The seasonality of three childhood infections in a pre-industrial society without schools

Background: The burden of many infectious diseases varies seasonally and a better understanding of the drivers of infectious disease seasonality would help to improve public health interventions. For directly transmitted highly-immunizing childhood infections, the leading hypothesis is that seasonality is strongly driven by social gatherings imposed by schools, with maxima and minima during school terms and holidays respectively. However, we currently have a poor understanding of the seasonality of childhood infections in societies without schools and whether these are driven by human social gatherings. Here, we used unique nationwide data consisting of >40 epidemics over 100 years in 18th and 19th century Finland, an agricultural pre-health care society without schools, to (i) quantify the seasonality of three easily identifiable childhood infections, smallpox, pertussis and measles and (ii) test the extent to which seasonality of these diseases is driven by seasonal social gatherings.

Methods: We quantified the seasonality of transmission using time series Suscpetibel-Infected-Recovery models, wavelet analyses and general additive mixed models.

Results: We found that all three infections were seasonal and the seasonality patterns differed from those in industrialized societies with schools. Smallpox and measles showed high transmission in the first half of the year, but we could not associate this with seasonal human gatherings events. For pertussis, however, transmission was higher during social gathering events such as New Year and Easter.

Conclusions: Our results show that the seasonality of childhood infections is more variable than previously described in other populations and indicate a pathogen-specific role of human social aggregation in driving the infectious disease dynamics.

Funding: Academy of Finland (278751, 292368), Nordforsk (104910), the Ehrnrooth Foundation, the Finnish Cultural Foundation, the University of Turku Foundation and the Doctoral Programme in Biology, Geography and Geology, University of Turku.

In this study we had two aims. First, we identified whether the three childhood infections show 134 seasonality (Aim 1). Second, we tested whether weeks with seasonal social gatherings have 135 increased infectious disease transmission (Aim 2). To accomplish these two aims, we used three 136 methods: (i) we estimated the seasonality of transmission rates from time-series Susceptible-137 Infected-Recovery (tSIR) models, (ii) to provide an independent confirmation of the tSIR results 138 in (i) we repeated the analysis using wavelet analyses, and (iii) we tested whether weeks with 139 seasonal social gatherings have increased infectious disease transmission using general additive 140 models (GAMs). 141 2.2.2 Method 1: Identifying infectious disease seasonality using tSIR-models 142 To estimate variation in disease transmission per era, we used the time-series susceptible-143 infected-recovered (tSIR) model as developed in Bjørnstad et al. (2002). To fit the tSIR model to 144 the data, we followed previously established techniques which are implemented in R version 145 4.0.2 (R Core Team 2020) with the function 'estpars' of the package 'tsiR' (Becker and Grenfell 146 2017). In brief, in the tSIR model the population is divided into Susceptibles S, Infected I and 147 Recovered R. The parameter of interest is the transmission factor β, which captures the rate at 148 which individuals transfer from the S to I and hence transmit the infection. The time series 149 approach divides time in discrete steps which reflect the diseases' generation time, which is two 150 weeks for smallpox and measles and four weeks for pertussis (Anderson and May 1991). It is 151 possible that there is a time lag between onset of disease and death, which we did not include in 152 the tSIR models. While having an estimate of such time lag from historical pre-industrial 153 societies is difficult, we note that tSIR models with generation times that were two weeks longer 154 than those presented here, a time possible lag in infants as indicated by some studies on 155 smallpox and pertussis in 20 th century industrialized societies (Downie et al. 1969;Centers for 156 Disease Control and Prevention 2002), gave results consistent with those presented here. 157 To estimate β, we first reconstructed the susceptible dynamics. At each time step, the number of 158 susceptibles St fluctuates around a mean Sm such that St = Sm + Zt. Based on pre-vaccination data 159 from the United Kingdom, Sm was fixed at 0.035 (Bjørnstad et al. 2002;Becker and Grenfell 160 2017). Note however that analyses with other values of Sm (0.015-0.15) gave conclusions 161 consistent with those presented here (results not shown). The deviation of susceptibles Zt+1 at 162 each time step is then given by: 163 -Zt+1=Zt+Bt-ρtIt 164 Bt and It are the number of births and infected respectively as given by the data. ρt is a time 165 varying correction factor for underreporting. However, because our data reports only deaths, ρt 166 here is a combined measure of underreporting and fatality rate. It can be estimated using the 167 regression between cumulative births and cumulative observed cases. Given these parameters, 168 (1) . CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. ; https://doi. org/10.1101org/10. /2021 the dynamics of the deviation susceptibles Zt+1 can be estimated by the residuals of the 169 regression model in (1), which for small populations as in pre-industrial Finland, is best 170 captured by a Gaussian process (Caudron et al. 2015). 171 Following the tSIR model, the infection dynamics at each time step can be described by: 172 log(E[It+1])=log(β)+ α log(It) +log (Sm+Zt) 173 where the expected number of infected individuals E[It+1] is given by the data. α is a correction 174 factor that accounts for the discretization of a continuous time process and for inhomogeneous 175 population mixing. We here fixed α=0.97, a common value that impacts transmission dynamics 176 independently of epidemic size (Bjørnstad et al. 2002;Caudron et al. 2015;Becker and Grenfell 177 2017 To provide an independent statistical test of the seasonality of infectious diseases, we 185 performed wavelet analyses. In brief, wavelet analysis decomposes time series data using 186 functions (wavelets) simultaneously as a function of both period and time, thereby making it 187 possible to detect changes in periodicity, such as annual seasonality, over time (Cazelles et al. 188 2007). We pooled the data into monthly intervals and tested periods between six and 18 189 months with 12 months showing annual seasonality. To determine statistical significance, we 190 compared the periodicity of the data with that of 1000 'white noise' datasets with significance at 191 the 5% level. We performed wavelet analyses following (Cazelles et al. 2007;Roesch and 192 Schmidbauer 2018), in which we log-transformed (+1) the monthly data and fitted a Morlet 193 wavelet with the function 'analyze.wavelet' using the package 'WaveletComp' (Roesch and 194 Schmidbauer 2018). To statistically test whether transmission rates were higher during weeks of seasonal social 198 gatherings events we started off with a trivariate (smallpox, pertussis, measles) general additive 199 model (GAM) using the function 'gamm' from the package 'mgcv ' (Woods et al. 2016;Woods 200 2017). We determined the 'statistical significance' of terms and models based on (i) the model 201 selection approach using the second order Akaike's Information Criterion (AICc) and (ii) classic 202 (2) . CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 p-values. Model fitting should be considered as a continuum for which alternative models within 203 4 ΔAICc are plausible and become increasingly equivocal up to 14 ΔAICc, after which they 204 become implausible (Burnham and Anderson 2002;Burnham et al. 2011). We also checked the 205 terms p-values at p=0.05. 206 In our models, the dependent variable was the infection's transmission rates estimated from the 207 tSIR-model. To capture the interannual variance in the transmission of infection, we estimated 208 the transmission coefficients from tSIR analyses in four-year intervals. We chose four-year 209 intervals because it was the shortest time interval between epidemics (Fig. 1, Briga et al. 210 unpublished), but note that analyses with any intervals up to 16 years (25% of the sample size) 211 gave conclusions consistent with those presented here (results not shown). We also ran all 212 analyses with four-year running means, which also gave conclusions consistent with those 213 presented here (results not shown). To avoid that some years influenced the results more than 214 others, we standardized all transmission rates to an annual mean of zero and, to avoid that the 215 results be driven by a few influential some seasonal maxima, we log-transformed (+1) all 216 transmission rates. 217 We had two types of predictor variables in our GAM models. First, we included 'social', a six-218 level factor that reflected whether the week belonged to either one of the five aforementioned 219 seasonal social gathering events or was 'asocial'. Second, because human seasonal social 220 gatherings correlate with seasonal climatic cycles (e.g., agricultural activities are inevitably 221 climate bound and social activity may be greater during summer) we tested whether the effect 222 of social gathering events on transmission remained statistically significant when including 223 gradual cyclic changes by adding the predictor variable 'daylength' as a smoothed term using 224 cubic regression splines (note that any spline gave conclusions consistent with those shown 225 here). We obtained the daylength data from 1750 until 1850 at the location of Helsinki from the 226 US National Oceanic and Atmospheric Administration (NOAA) daylength calculator at 227 https://gml.noaa.gov/grad/solcalc/index.html. 228 We started off our GAM models with a trivariate analysis to capture a possible covariation 229 between dependent variables. However, because multivariate GAMs do not allow inclusion of 230 random terms or temporal autocorrelation of the time series, we performed follow-up 231 univariate analyses including temporal autocorrelation as an auto-regressive factor of the order 232 one (corAR1; -17.4<ΔAICc<-79.7) and year as a random intercept even though it provided a 233 somewhat worse model fit (ΔAICc=+2 for all infections). The co-variation between the 234 transmission rates of all three infections remained low (<0.1 in the best fitting model in Table  235 S1) and conclusions from multivariate and univariate analyses were always consistent, though 236 . CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. ; https://doi.org/10.1101/2021.10.08.21264734 doi: medRxiv preprint multivariate models were often more conservative than univariate models (Tables S1 & S2), and 237 we formulated our conclusions based on the most conservative approach. We tested a variety of 238 smoothing functions, which gave consistent results and here we present all models using cubic 239 splines 'cs'. Model residuals fulfilled all assumptions as controlled with the functions 'resid' and 240 'gam.check' from the package 'mgcv ' (Woods et al. 2016;Woods 2017).  We found evidence of seasonality for all three infections in both eras (Fig. 2, Fig. S1, Fig. S2). For 251 smallpox and measles, transmission (i.e. β in the tSIR models) was generally higher in the first 252 half of the year than the second, in both eras ( Fig pre-vaccine era, transmission in the first half of the year was 6% higher than the era mean, 254 while in the vaccine era the change was larger at 19% higher than the era mean ( Fig 2B, Fig. S1B). Some of the seasonal maxima were larger than the 95% CI around the eras´ mean 264 transmission ( Fig. S1B), indicating statistically significant seasonality and this was confirmed 265 through wavelet analyses (Fig. S2B). Hence, all three infections were seasonal, but the 266 seasonality of pertussis was different from that of smallpox or measles. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. ; https://doi.org/10.1101/2021.10.08.21264734 doi: medRxiv preprint Next, we tested the hypothesis that the transmission of childhood infections was higher during 269 weeks of seasonal social gathering. In our study population, we identified five periods with 270 seasonal social gathering events (grey weeks in Fig. 3), for a total of 32 out of 52 weeks (62%) 271 having recurring social gathering events. These data require statistical testing for at least three 272 reasons. First, we separated the dataset in two eras, and while this may be justified from a 273 public health perspective, it remains unclear whether the introduction of the smallpox vaccine 274 and the changing demography of the growing population affected the seasonality of childhood 275 infections. Second, differently from school-term studies, where schooling occurs in relatively 276 standardized conditions, the social gathering events in this study differed for example in their 277 importance, group size, whether they were celebrated indoor vs. outdoor and some might 278 increase transmission more than others. Third, seasonal social gatherings can correlate with 279 seasonal changes in climate and we should control to what extent social events are confounded 280 with seasonal changes in climate. To address these challenges, we estimated tSIR-based 281 transmission rates in four-year intervals and tested statistically whether transmission rates 282 were higher during social gathering events using GAMs. 283 The results did not support increased transmission of smallpox and measles during weeks with 284 social gathering events (smallpox: ΔAICc=+9.6, p=0.85, model 37 in Table S1, ΔAICc=+9.4, 285 p=0.82, Table S2; measles: ΔAICc=+4.7, p=0.09, model 57 in Table S1; ΔAICc=+1.8, p=0.13, Table  286 S2). Furthermore, support for a gradual seasonal cycle was inconclusive at best, i.e. models with 287 the daylight term were within 4 AICc of models without and the p-values for the daylight term 288 in all these models were >0.05 (smallpox: ΔAICc=+2.2, p=0.96, model 21 in Table S1, 289 ΔAICc=+1.8, p=0.12 in Table S2; measles: ΔAICc=-1.9, p=0.12, model 53 in Table S1, ΔAICc=+2.0, 290 p=0.20 in Table S2). For neither infection did we find statistical support for era-specific 291 seasonality (ΔAICc=+2.0, p=0.99, Table S2A & C), era-specific social gathering effects 292 (ΔAICc>+10.1, p>0.54, Table S2A & C) or era-specific daylight effects (ΔAICc=+4.0, p>0.34, Table  293 S2A & C). 294 In contrast, pertussis showed increased transmission at the end of December and early January, ΔAICc=-18.7, p=0.0046, model 53 in Table S1; univariate: ΔAICc=-7.4, p=0.01, Table S2; Fig. 3B, 299 Fig. S3A), but the latter two transmission periods were better described by a gradual cycle (i.e. a 300 smooth term) rather than an abrupt event (multivariate: ΔAICc=-114.8, p<2e -16 , model 53 in 301   Table S1; univariate: ΔAICc=-23.3, p<e -8 , Table S2 is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. ; https://doi.org/10.1101/2021.10.08.21264734 doi: medRxiv preprint summer-winter gradient we would expect from a uniquely climate driven seasonal 304 phenomenon. The winter minimum in the smooth term did not cause the statistically significant 305 increase in transmission during New Year as this increase remained statistically significant in a 306 model without seasonal cycle (multivariate: ΔAICc=-11.8, p=0.0017, models 63 vs. 64 in Table  307 S1; univariate: ΔAICc=-8.8, p=0.004 in Table S2B). We found no statistical support for era-308 specific seasonality (ΔAICc=+34.0, p=0.97, Table S2B), era-specific social gathering effects 309 (ΔAICc=+35.1, p=0.94, Table S2B) or era-specific daylight effects (ΔAICc=+37.1, p=0.01, Table  310 S2B). Thus, pertussis showed three seasonal maxima during New Year, in late March/April and 311 in August/September, and the New Year maximum was abrupt, while the March/April and 312 August/September maxima were gradual. 313

Discussion 314
In this study we quantified the seasonality of three childhood infections in a pre-industrial 315 society without schools, 18 th and 19 th century Finland. For all three infections, we detected 316 statistically significant seasonality. For smallpox and measles, transmission was higher in the 317 first half of the year than in the second half, but we found no evidence that this seasonality was 318 associated with social gathering events. In contrast, pertussis transmission increased at (i) the 319 end of December and early January, (ii) in March and April and (iii) in August and September. 320 The first two maxima were associated with the social gathering events New Year and Easter. 321 Here we discuss the implications of our results for the seasonality and dynamics of infectious 322 diseases in general and of childhood infections in particular. 323 A large body of work on the seasonality of measles in pre-vaccine Europe and in contemporary 324 industrialized societies found higher transmission in winter, spring and autumn and minima 325 during spring and summer breaks, indicating that seasonality is driven by school term-forcing 326 (Bjørnstad et al. 2002;Metcalf et al. 2009;Mahmud et al. 2017;Klinkenberg et al. 2018). Our 327 results show that two hundred years earlier, the seasonality of measles was very different from 328 that in contemporary industrialized European societies with school terms. Measles showed 329 increased transmission in the first half of the year without association with seasonal social 330 gathering events. We believe that the lack of association with social gathering events is because 331 Finland in the 18 th and 19 th century can best be described as a meta-population, with different 332 regions of Finland regularly being disconnected from one another and from neighboring 333 countries. In this meta-population, infections would often experience local or even nation-wide 334 extinctions (Fig. 1), only to be re-introduced from other regions in Finland or from abroad. In 335 such a scenario, the seasonality of measles could be driven by larger-scale human movement, 336 with the increased early-year transmission (Fig. 2)  is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. ; https://doi.org/10.1101/2021.10.08.21264734 doi: medRxiv preprint seasonal human movements driving the seasonality of measles have been described, for 339 example, in Niger, a contemporary largely agricultural based society (Ferrari et al. 2010;Bharti 340 et al. 2011). We suggest that the seasonal dynamics of measles in our population could reflect a 341 similar phenomenon. 342 We know almost nothing on the seasonality of smallpox epidemics (Krylova and Earn 2020). A 343 recent study on the seasonality of smallpox mortality in London indicated a weak seasonal 344 trend from summer in the 1700's that gradually moved towards winter in the 1800's (Krylova 345 and Earn 2020). Our results quite differ from those in London in the same period, with smallpox 346 transmission being predominantly in the first half of the year, a seasonality consistent with that 347 of measles (Fig. 2). Interestingly, in the 18 th and early 19 th smallpox in London was endemic 348 (Duncan et al. 1994;Krylova and Earn 2020). This is different from the situation in the Finnish 349 meta-population and we believe that the seasonality of smallpox is driven by the same dynamics 350 as that of measles, namely early-year introductions from neighboring regions through large-351 scale human movement. 352 In contrast to the predictable seasonality of measles, the seasonal dynamics of pertussis in 353 contemporary industrialized societies are more stochastic, with a seasonality that does not 354 match periods of school-term aggregation well (Farizo et al. 1992;Rohani et al. 2002;355 Skowronski et al. 2002;Metcalf et al. 2009). In North America there appears to be a dominant 356 seasonality during August, September or in winter (Skowronski et al. 2002;Fisman et al. 2011). these earlier studies in that we found increased transmission around New years and in late 360 March/April, which matched the New Year and Easter social gatherings respectively. However, 361 we also found a maximum from August until October (Fig. 3) and this is consistent with many 362 previous studies, even in contemporary societies with a relatively high vaccination coverage 363 (Farizo et al. 1992;Skowronski et al. 2002;Metcalf et al. 2009). In our study, the timing of that 364 maximum did not match any seasonal social gathering events, which occurred either before or 365 after the maximum in transmission (Fig. 3)  It remains puzzling as to why local social gathering events did not drive the seasonality of 370 measles and smallpox, but did affect the seasonality of pertussis. The difference between these 371 childhood infections may be caused by several factors and here we bring forward two potential 372 mechanisms. First, the degree of immunization differs between infections, with immunity to 373 . CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. There are at least three limitations to our study. First, the causes of death in historical 389 records were registered by parish priests, who may or not have the know-how to identify 390 infections. To this end, we focused our analyses on three easily identifiable causes of death. 391 Second, individuals might be more likely to die from infections in some seasons rather than 392 others, irrespective of infection incidence. In the Northern hemisphere, the general pattern is 393 that mortality gradually increases in winter (Rau 2007). The seasonal infection maxima in our 394 study did not match that pattern and hence are unlikely to be driven by general seasonality in 395 mortality. To what extent there is a possible discrepancy between the seasonality of incidence 396 and mortality in our historical data remains yet unclear, but we are gathering historical 397 incidence data to investigate the question. Secondly, there can be seasonal variation in access to 398 health care or registration of information on infectious diseases. We do not believe this to be the 399 case here, as health care was virtually non-existent and because deaths were recorded by the 400 parish priests, who were locally-based. However, even though we believe that seasonal 401 variation in registration of information was limited, it remains possible that for example parish 402 priests were postponing recording information when they were busy with other activities, such 403 as New Year or Easter. Interestingly though, for a subset of the data the time lag between death 404 date and burial date during these two gathering events did not differ from other periods of the 405 year (results not shown). 406 To conclude, here we provide a quantification of the seasonality of three childhood infections 407 smallpox, pertussis and measles in a society for which data are rarely available, an 18 th and 19 th 408 . CC-BY 4.0 International license It is made available under a perpetuity.
is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. ; https://doi.org/10.1101/2021.10.08.21264734 doi: medRxiv preprint century pre-industrial European nation. Our results show that the seasonality of infections was 409 very different from that in 20 th century industrial societies and emphasize a pathogen-specific 410 role of local social movements as a driver of infectious disease dynamics. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. ; https://doi.org/10.1101/2021.10.08.21264734 doi: medRxiv preprint Table 1. Overview of the infections in this study, with their sample sizes before and after the 568 start of the smallpox vaccination program, respectively the pre-vaccine era (1750-1801) and 569 vaccine era  is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021  is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021     is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. ;https://doi.org/10.1101https://doi.org/10. /2021 Supplementary Information S1: Infectious disease seasonality per vaccine era based on 618 tSIR models 619 620 is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The copyright holder for this this version posted October 10, 2021. ; https://doi.org/10.1101/2021.10.08.21264734 doi: medRxiv preprint Fig. 3. 633 Table S1. Model selection of tri-variate general additive model of tSIR-based transmission rates for the three infections. Predictor variable 'social' 634 captures whether or not week belonged to one of five seasonal social gathering events (six-level factor for five social gathering events and asocial 635 weeks). To account for a possible climatic seasonal dynamic in the transmission of infection we included daylight as a cubic smoothed term shown in 636 the whether or not week belonged to one of five seasonal social gathering events (six-level factor 642 for five social gathering events and asocial weeks). To account for a possible climatic seasonal 643 dynamic in the transmission of infection we included daylight as a cubic smoothed term shown 644 in the table as s(daylight). For each infection, the best model is shown in italic. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint