Comparative dynamics, seasonality in transmission, and predictability of childhood infections in Mexico

Summary

The seasonality and periodicity of infections, and the mechanisms underlying observed dynamics, can have implications for control efforts. This is particularly true for acute childhood infections. Among these, the dynamics of measles is the best understood and has been extensively studied, most notably in the UK prior to the start of vaccination. Less is known about the dynamics of other childhood diseases, particularly outside Europe and the United States. In this paper, we leverage a unique dataset to examine the epidemiology of six childhood infections – measles, mumps, rubella, varicella, scarlet fever and pertussis – across 32 states in Mexico from 1985 to 2007. This dataset provides us with a spatio-temporal probe into the dynamics of six common childhood infections, and allows us to compare them in the same setting over the same time period. We examine three key epidemiological characteristics of these infections – the age profile of infections, spatio-temporal dynamics, and seasonality in transmission – and compare them with predictions from existing theory and past findings. Our analysis reveals interesting epidemiological differences between the six pathogens, and variations across space. We find signatures of term-time forcing (reduced transmission during the summer) for measles, mumps, rubella, varicella, and scarlet fever; for pertussis, a lack of term-time forcing could not be rejected.

INTRODUCTION

The dynamics of directly-transmitted, acute infections are well-described by the Susceptible-Infected-Recovered (SIR) model [1]. Among acute childhood infections, the epidemic dynamics of measles are perhaps the best understood since it has received extensive study, most notably in England and Wales prior to the start of vaccination. This body of research has shown that the basic epidemic clockwork of measles is driven by the aggregation of children during school terms [2]. Measles dynamics are also tuned by birth rates as modulated by vaccination rates [3]. If both of these factors are encompassed into simple mechanistic models that reflect the SIR, long-term prediction of incidence is possible [4]. Less is known about the dynamics of other childhood diseases, particularly outside Europe and the United States. Here, we leverage a unique existing dataset to examine the epidemiology and incidence dynamics of six childhood infections – measles, mumps, rubella, varicella, scarlet fever and pertussis – across 32 states in Mexico from 1985 to 2007. The availability of data for each state provides 32 similar populations for comparing varying social and demographic contexts across the six pathogens. We test the extent to which the dynamics of these common childhood infections align with predictions from theory and past findings.

Measles, mumps, rubella, and varicella are viral infections, while scarlet fever and pertussis have a bacterial aetiology. All six are directly transmitted airborne diseases that have highest incidence in children. The pathogens, however, have several important life-history differences. There is considerable variation in the duration of generation times (latent period plus infectious period), transmission rates, the degrees of immunity conferred, and the existence of asymptomatic carriers across the six infections; the main differences and similarities are summarized in Table 1. Our focus is on examining the role of three key epidemiological characteristics of these six common childhood infections in the context of existing theory

Table 1. Comparison of the six infections [5]

Age profile of infection

The mean age of infection and the age-specific force of infection (FOI) are core epidemiological parameters. The variation of FOI with age has been documented for several childhood diseases [5] and has implications for control efforts [6]. Elimination following vaccination is easier if the infection is intrinsically less transmissible in older age groups than it would be if transmission were constant across all ages [5].

In the basic SIR model, at equilibrium, the average age of infection is inversely proportional to the FOI (and the underlying rate of disease transmission) [5]. Thus, we would expect denser states, and those with higher birth rates, to have a lower average age of infection due to higher rates of disease transmission. The mean age of infection could also be associated with socioeconomic factors. For instance, if school attendance rates vary across states, we would expect to see a corresponding variation in the mean age of infection; if fewer children are concentrated in schools, then the transmission rates may be lower, resulting in higher mean ages of infection. Accordingly, we present analysis of how the estimates of the mean age of infection vary with space, demography and socioeconomic factors.

Spatio-temporal patterns in incidence

Given life-history differences for the six pathogens, we may expect to see differences in the infection dynamics. Previous work has suggested that diseases with large values of the basic reproduction number, R ₀, tend to exhibit large resonant power (cycles with annual or multiennial periods), while diseases with small amplitude of seasonality or small R ₀ tend to exhibit significant non-resonant power (which can occur at any frequency) [7]. For instance, measles (large R ₀) epidemics are typically strongly annual or biennial, whereas pertussis (small amplitude of seasonality) often exhibits significant non-resonant power. Pertussis also tends to exhibit longer period cycles of 3–4 years due to a longer generation time [8]. Multi-annual cycles with varying periodicity have also been observed for mumps [9], rubella [5], and scarlet fever [10]; annual cycles have been observed for varicella incidence [11]. The dynamics of these pathogens have rarely been compared in the same setting over the same time period, which is possible with this unique dataset.

We also examine the epidemic synchrony between states and the timing of epidemics across states. Previous work has suggested that the strong synchronizing effect of seasonally forced transmission should lead to epidemics that are completely synchronized across large coupled areas [12]. However, travelling epidemic ‘waves’ have been observed for measles in England and Wales, arising as a result of repeated invasions from endemic core areas to the periphery [13]. This is thought to be a result of the coupling of large and small population centres. Travelling waves have also been observed for pertussis in the United States in the 1950s, although no similar spatial pattern was detected during the more recent re-emergence of the pathogen [14]. While our scale of observation (state-level) is too large to detect core-periphery dynamics, we can nonetheless examine the synchrony across states for each disease. Comparing across diseases allows us to comment on the extent to which the strength of seasonal forcing is correlated with synchrony, as predicted by previous theory [12].

Seasonality in transmission

Many childhood infectious diseases exhibit a marked seasonal pattern in incidence. In the SIR model, for the oscillatory behaviour to persist over a long period, as is observed in incidence data for many pathogens, seasonal ‘forcing’, or variation, in host births [15] and/or in transmission is necessary. Seasonal fluctuations in transmission for directly transmitted pathogens can occur for several reasons. These include changes in host behaviour [16], annual variation in environmental conditions [17, 18], and variations in host susceptibility [19]. Environmental conditions, such as temperature and humidity, may affect pathogen survival outside the host and the dispersal of infectious particles [20]. However, the link between environment factors and disease dynamics has not been shown conclusively for any of the pathogens considered here.

For childhood infections, there is evidence to suggest that changes in host behaviour, rather than environmental conditions, may be the key driver of seasonality in transmission [21, 22]. Seminal work on measles in England and Wales, demonstrated that the aggregation of children during school term time was a strong driver of the observed seasonality in incidence [2, 6, 16]. For measles in Niger, the strength of the seasonal forcing of transmission was shown to be uncorrelated with rainfall patterns but correlated with population size [22]. This is consistent with the hypothesis of changes in host behaviour (in this case, due to agricultural cycles) driving seasonality in transmission. The evidence for other pathogens, however, is less conclusive. A reduction in the transmission parameter during the summer months has been estimated for other childhood infectious diseases including mumps, varicella, and pertussis [16, 23]. Metcalf et al. [24] found evidence for school term-time forcing for viral infections (measles, mumps, varicella), but not for bacterial infections (diphtheria, scarlet fever, pertussis) in Copenhagen before the start of vaccination.

To examine seasonality in transmission, we fit a mechanistic model to the incidence data to estimate the unobserved transmission rates. Comparing across different pathogens, with varying transmission rates and mean ages of infection, over the same time period, allows us to better understand the drivers of seasonality for childhood infections. If term-time forcing were driving the seasonality in transmission, then we would expect to see consistently reduced transmission during the long summer holidays.

In addition to quantifying these three key epidemiological characteristics for each infection, we examine whether they are related to the states’ demographic and socioeconomic conditions. Socioeconomic factors may affect disease dynamics for a variety of reasons. For instance, they may affect host susceptibility (children in poorer states are more likely to suffer from malnutrition) or host aggregation patterns (states differ in schooling rates, spatial aggregation patterns, and agricultural practices). Only a few studies have closely examined these hypotheses. Metcalf et al. [25] found no association between a states’ socioeconomic characteristics and rubella dynamics in Mexico. On the other hand, scarlet fever epidemics in England and Wales in the nineteenth century were strongly correlated with wheat prices, possibly as a result of malnutrition [10]. Here, we examine the associations between socioeconomic factors and disease dynamics as a first step towards better understanding the role of these drivers in shaping the observed dynamics.

We first estimate the mean age of infection and the age-specific FOI from the age-structured incidence data. Next, we examine spatio-temporal patterns in the incidence data by examining the periodicity of the epidemic cycles, the epidemic synchrony between states, and the timing of epidemics across states. Last, we fit a mechanistic model to the incidence data for each disease to estimate the unobserved transmission rates, and examine the feasibility of making predictions using this model. We discuss how our results compare with predictions from previous theory and findings. Our analysis reveals interesting epidemiological differences between the six pathogens. With the exception of pertussis, we detect signatures of term-time forcing for all the pathogens.

METHODS

Data

Incidence data for the six infections was obtained from the Yearly Morbidity Report (2008) of the Directión General de Empidemiología, Mexico [26]. The dataset contains monthly counts for the number of infected individuals between 1985 and 2007 for each of the 32 states of Mexico for each infection (Supplementary Figs S1–S6). Figure 1 shows the incidence time-series for all states combined. For measles, data is missing for 1989; for scarlet fever, data is missing for 1985. We also have state-level yearly incidence data between 1990 and 2007 broken down by age groups of 5–10 years (<1, 1–4, 5–14, 15–24, 25–44, 44–64, ⩾65 years) for all the infections, except for measles for which no age-structured data is available. Widespread vaccination for measles and pertussis began around 1991, and reached relatively high levels of coverage [27, 28]. The MMR vaccine was introduced in 1998 and supplanted the monovalent measles vaccine [27]. Routine vaccination for varicella is not recommended in Mexico, and there is no vaccine for scarlet fever.

Fig. 1. Reported incidence for all states in Mexico from 1985 to 2007 for (a) measles, (b) mumps, (c) rubella, (d) varicella, (e) scarlet fever, and (f) pertussis. Approximate vaccination start date is indicated with the red line for measles, mumps, rubella, and pertussis. The monovalent measles vaccine and the pertussis vaccine were introduced in 1991; the MMR vaccine (which supplanted the monovalent measles vaccine) was introduced in 1998.

We obtained demographic information and socioeconomic indicators for each state from the Instituto Nacional de Estadística y Geografía. In our analysis, we used population size, population density, and crude birth rate (averaged over the time period) as our demographic variables for each state. We used the average of the estimated poverty rates (measured as food poverty and asset poverty) in 1990, 2000 and 2010 as our main socioeconomic indicator.

Age profile of infection

We estimated the age-specific FOI for mumps, rubella, varicella and scarlet fever from the age-structured data (for all states combined) using the ‘catalytic framework’ described in detail in [29]. There is no age-specific data available for measles, and the age-specific data for pertussis is too sparse for accurate estimation of the mean age of infection. We used a piecewise constant model to estimate the age-dependent FOI, assuming a constant FOI within pre-determined age intervals. We assumed a binomial distribution for the proportion infected at each age, and numerically maximized the likelihood to estimate age-dependent values of FOI that best fit the age-structured prevalence data. Using the estimated age-specific FOI, we can also estimate the mean age of infection for each disease for all states combined.

We also used the age-specific incidence data to estimate the mean age of infection for mumps, rubella, varicella and scarlet fever for each state separately. While the approach described above works well for all states combined, the maximum-likelihood estimates are difficult to obtain when the data is sparse. To estimate the state specific mean ages of infection, we estimated the cumulative distribution function by fitting a cubic smoothing spline to the cumulative proportion of cases at each age (calculated by summing the data over age and dividing by the total sum). We estimated the derivative of the fitted curve to obtain an estimate of the probability density by age (scaled to sum to one). The mean age of infection can then be calculated as the numerical mean of the probability density curve.

We used ordinary least squares to estimate associations between the mean age of infection and the district's population size, population density, crude birth rate, and poverty rates. We examined the associations of mean age of infection with each of these factors separately, but also included them as multiple predictors in the same regression model. We also examined any spatial gradient in the estimated mean age of infection.

The estimates of the average age of infection also allows us to calculate a rough estimate of the basic reproduction number, R ₀, for each disease in every district. Assuming a constant FOI and Type I mortality (constant survivorship over age until life expectancy is reached), R ₀ = μ ⁻¹ /A, where μ is the per capita birth rate and A is the mean age of infection [5]. Although this estimate of R ₀ may be biased due to many factors [30], it provides an independent estimate that can be compared with the time-series susceptible-infected-recovered (TSIR) estimate of R ₀, to ensure that the model fit is producing reasonable estimates.

Spatio-temporal analysis

We used wavelet analysis, described in detail in [31], to examine the period of oscillations in the incidence data over time. Unlike Fourier analysis, wavelet analysis can capture the non-stationarity in the period of recurrent outbreaks over time. This allows us to describe how the period of incidence cycles vary over time and across diseases. We log(x + 1)-transformed the incidence time-series to dampen the amplitudes. We used a Morlet wavelet and applied the wavelet transform to the transformed incidence data to compute the power as a function of period and time. Statistical significance was tested against red noise background (signal noise produced by random-walk motion).

We also examined the synchrony in recurrent epidemics between states for each disease. We defined epidemic synchrony between two states as the Pearson correlation coefficient between the two incidence time-series. The incidence time-series was square-root-transformed to stabilize the variance. We examined how synchrony was related to distance between states (as measured by the latitude and longitude of the centre of the state) by estimating spatial correlation functions [13, 32]. We used a nonparametric spline function to estimate the correlation functions, and computed confidence limits using bootstrapping (with 500 resamples) [33]. For diseases where vaccination was implemented during the observed timeframe we computed the spatial correlation functions separately for pre-vaccination and vaccination-era data. For each disease, we also computed the cross-correlation coefficient between all states and the Federal District (Mexico City), to see if there were any lags between the incidence time-series.

To examine any spatial gradient in the timing of epidemics, we computed the ‘centre of mass’ of the incidence time-series. The centre of mass represents the average month in the year in which the most cases occurred over the entire observation time period. Bootstrapped confidence intervals were computed with 1000 resamples, by computing the maximum-likelihood estimate for the mean of a von Mises distribution [34] for each resampled dataset. We examined the correlation between the centre of mass for each state and the latitude and longitude of the centre of the state, and repeated the analysis for each disease.

The TSIR framework

To estimate seasonal transmission rates we used the TSIR model developed by [4]. The TSIR model is a stochastic, discrete time analogue of the classic SIR model. The time-series of number of infected and the number susceptible are described by a set of difference equations. At each time step, the expected number of infected cases is given by

(1) $$E[I_{t + 1}] = \displaystyle{{\beta _{s\;} I_t^\alpha \; S_t} \over {N_t}},$$

where E[I _t+1] is the expected number of infected individuals one infection generation time in the future; I _t and S _t are the numbers infected and susceptible at time, t, respectively; N _t is the total population size at time, t; β _s is the seasonal transmission factor; and α captures heterogeneities in mixing, not captured by the seasonality, and the effects of discretization.

Assuming all individuals eventually become infected, the number of infected individuals tracks the births, B, and the number of susceptible individuals is defined by:

(2) $$S_{t + 1} = S_t + B_t - I_t + u_t.$$

Here, u _t is additive noise, with E [u _t ] = 0. By reducing the susceptible population, vaccination has the same effect as lowering births. Vaccine coverage levels were assumed to be 80% for measles, mumps and rubella, and 90% for pertussis. Accordingly, the number of births in the TSIR model was discounted by 0·8 for measles, mumps, and rubella, and 0·9 for pertussis following the introduction of vaccination. Our main model results also hold over a range of assumed vaccination rates for these infections.

If the number of susceptible individuals fluctuates around a mean, then $S_t = \; \bar {S} N_t + \; Z_t$ , where $\bar {S} $ is the average proportion of susceptible individuals in the total population, N _t , and Z _t is the unknown deviation around the mean number of susceptible individuals. We can rewrite the susceptible difference equation in terms of the deviations, Z _t , and iterate successively with an initial starting condition, Z ₀. This yields:

(3) $$\mathop \sum \nolimits_{k = 0}^{t - 1} B_k = - Z_0 + \mathop \sum \nolimits_{k = 0}^{t - 1} I_k + Z_t + u_t.$$

Assuming u _t is small, Z _t can be estimated as the residuals from the locally varying regression of the cumulative number of births on the cumulative number of cases. We can then rewrite equation (1) as:

(4) $$E[I_{t + 1}] = \displaystyle{{\beta _{s\;} I_t^\alpha \; \left( {\,\bar {S\;} N_t + Z_t} \right)} \over {N_t}}.$$

The mean proportion of susceptible individuals, $\bar{S\;}$ , can be estimated using marginal profile likelihoods from estimating equation (4) for a range of values of $\bar{S\;}. $ We used a quasi-Poisson generalized linear model to estimate equation (4). We constrained the estimated $\bar{S\;}$ to have an upper limit of 0·8 (i.e. 80% of the population). Conditional on the estimated $\bar {S\;} $ , the seasonal transmission rates, β _s , were estimated using a quasi-Poisson generalized linear model with monthly or biweekly-specific factors (and no intercept). Since the TSIR model is not age-structured, here we are estimating the ‘effective’ seasonal transmission rates that encompass age structure [2].

We first corrected the incidence time-series for under-reporting using methods described in [4]. Assuming that all children eventually become infected in the absence of vaccination, the reporting proportion can be estimated as the slope of the regression line relating cumulative number of births to cumulative number of cases (see [4] for details). To allow for local variation in the regression of cumulative number of births on cumulative number of cases, we used a cubic smoothing spline. This allows for temporal variation in the reporting proportion. Having corrected the incidence data for under-reporting, we then fitted the TSIR model to the mumps and pertussis time-series using monthly time steps, which is close to the estimated generation times for these pathogens [5]. As a robustness check, we also used biweekly time steps for mumps (generation time is estimated to be between 16 and 26 days [5]), and the main results remained unchanged (Supplementary Fig. S10). For measles, rubella, varicella and scarlet fever, we fitted the TSIR model to the incidence time-series using biweekly time steps to match their shorter generation times. To convert the observed monthly incidence into biweekly time steps we used linear interpolation to produce observations for 24 time points in the year. This method maintains the peaks in the incidence data, although our results are robust to using cubic interpolation. Assuming that the time step used in the TSIR model is the true average time from infection to recovery, the estimates of the seasonal transmission rates, β _s , correspond to the seasonally varying basic reproductive number, R ₀. We fit the TSIR model to the incidence data for all states combined, as well as for each state separately.

Although our main focus was to analyse seasonal patterns in transmission, we also explored the ability to predict long-term incidence dynamics using the TSIR parameter estimates. For each disease, we fitted the TSIR model to the reported cases after correcting for under-reporting. We used the corresponding estimates of the seasonal transmission rates to forecast incidence for the observed time period. Predictions were made by sampling the incidence, I _t+1, from a negative binomial distribution:

(5) $$I_{t + 1} \sim NB\left( {\displaystyle{{\beta _{s\;} I_t^\alpha \; S_t} \over {N_t}},\; I_t} \right),$$

with mean $E[I_{t + 1}] = \displaystyle{{\beta _{s\;} I_t^\alpha \; S_t} / {N_t}}$$ and shape parameter I _t . We used the observed I ₀, corrected for under-reporting, and estimated S ₀ as the initial conditions, and allowed the simulation to continue until the end of the observed time-series. To assess the accuracy of the predicted incidence time-series, we calculated the normalized root mean squared error (NRMSE) between the predicted incidence and the actual incidence (normalized by dividing the root mean squared error by the average incidence over the observed time period).

RESULTS

Age profile of infection

The mean age of infection estimated from age-structured case data for all states combined was 6·47 for mumps, 6·52 for rubella, 6·78 for varicella, and 5·46 for scarlet fever (Fig. 2). The age-specific FOI was highest in the 4–14 years age group for mumps and scarlet fever. This is in the same range as what has been found for other childhood diseases in the past [5, 29]. The estimated mean age of infection for varicella was within the 5–14 years age range that has been estimated in the past [5, 11]. The FOI, however, was highest in the 1–4 years age group for varicella. Somewhat surprisingly, for rubella, the highest FOI was in the 0–1 years age group, although the confidence intervals for the estimates are largest for rubella. The 4–14 years age group has the second highest FOI for rubella.

Fig. 2. Top row: Estimated force of infection (FOI) by age for (a) mumps, (b) rubella, (c) varicella, and (d) scarlet fever for all states combined. Bottom row: Cumulative proportion of cases by age for (e) mumps, (f) rubella, (g) varicella, and (h) scarlet fever for all states combined. For the same age, the different points correspond to different years. The piecewise fit is shown in red. Mean age of infection, A, for each disease was calculated from the fitted model.

The mean age of infection estimated separately for each state, by fitting a smoothing spline to the cumulative proportion of cases by age, were lower than the mean age of infection estimated for all states combined using the catalytic framework. This is likely because there is less year-to-year variation in the cumulative proportion of cases by age within states (as opposed to when we look at all states together), and the spline function is overfitting to the initial increase in the cumulative proportion of cases from the 0–1 to 1–4 years age group. The mean age of infection, estimated separately by state, was not significantly associated with population size, density, or the crude birth rate. However, we found significant geographical variation in the mean age of infection (Figs 3 and 4). The mean age of infection estimated separately for each state was significantly associated with the longitude of the centre of the state for mumps [A = 12·39 (1·81) + 0·06 (0·02) longitude; P = 0·001, r ² = 0·27, s.e. in parentheses], rubella [A = 11·32 (2·78) + 0·08 (0·03) longitude, P = 0·008, r ² = 0·19, s.e. in parentheses], varicella [A = 14·31 (1·94) + 0·10 (0·02) longitude, P = 0·00001, r ² = 0·46, s.e. in parentheses], and scarlet fever [A = 8·35 (2·02) + 0·05 (0·02) longitude, P = 0·023, r ² = 0·13, s.e. in parentheses]. There were also significant associations of mean age of infection with latitude for all four infections, although the strength of the association was weaker. The association between latitude/longitude and the mean age of infection remained significant, and of similar magnitude, even after controlling for the population size and density. The mean age of infection was significantly positively associated with (logged) poverty for mumps (estimated association = 0·67, P = 0·005), rubella (estimated association = 0·95, P = 0·008), and varicella (estimated association = 0·76, P = 0·014).

Fig. 3. Estimated mean age of infection plotted against longitude for each state. Ordinary least squares regression fit and 95% confidence intervals are shown in red.

Fig. 4. Estimated mean age of infection by state. Darker green colours indicate a higher mean age of infection. The key indicates the upper age bound for each colour category. Red circle indicates the location of Mexico City.

The estimated R ₀ (averaged over all states) was 6·89 (variance = 0·97) for mumps, 13·26 (variance = 18·4) for rubella, 9·92 (variance = 3·87) for varicella and, 11·83 (variance = 3·68) for scarlet fever. We observe the same spatial gradient in R ₀ (controlling for the birth rate and population density) as we do for the mean age of infection, with states in the east having lower R ₀ compared to states in the west.

Spatio-temporal dynamics

Timing of outbreaks

Figure 5 shows the median number of cases by month, for consecutive pairs of years (pre-vaccination) for each disease for all states combined. Measles and pertussis had the greatest variance in monthly incidence over the observed time period. Measles and varicella, and to some degree, rubella, show clear annual patterns, while the other pathogens exhibit a slightly more complicated pattern. For all of the infections, incidence fell during the summer school holidays. With the exception of mumps, all infections had peak incidence in the spring around April. For mumps, we observed biannual peaks. There was a second peak just at the start of summer holidays, with the second peak being higher than the April peak every other year. For scarlet fever, the annual peak incidence switched from April to June every other year.

Fig. 5. Median number of cases by month (and interquartile range) for consecutive pairs of years in the pre-vaccination time-series for each disease (for all states combined). Approximate school holidays are shown in grey. n, number of years of data used for each infection.

All diseases had similar epidemic centres of mass, around April. Measles incidence had the earliest centre of mass (mid-March), and rubella incidence had the latest centre of mass (late April). Overall, we found no evidence of a relationship between the epidemic centre of mass and demographic factors. However, centre of mass (CoM) was significantly negatively correlated with longitude (Fig. 6) for measles [CoM = −0·85 (0·58) – 0·02 (0·01) longitude, P = 0·001, r ² = 0·29, s.e. in parentheses], varicella [CoM = 0·92 (0·18) – 0·004 (0·001) longitude, P = 0·016, r ² = 0·15, s.e. in parentheses] and, scarlet fever [CoM = −0·22 (0·58) – 0·02 (0·01) longitude, P = 0·007, r ² = 0·19, s.e. in parentheses]. The centre of mass was significantly positively associated with poverty (even after controlling for population size and longitude) for rubella (estimated association = 0·08, P = 0·029), varicella (estimated association = 0·05, P = 0·044), and scarlet fever (estimated association = 0·18, P = 0·032). Combined with our earlier findings, this result indicates that states that had higher rates of poverty, on average, had higher mean ages of infection (i.e. lower rates of disease transmission) and experienced epidemics later in the year (Supplementary Fig. S7).

Fig. 6. Centre of mass against longitude for (a) measles, (b) mumps, (c) rubella, (d) varicella, (e) scarlet fever, and (f) pertussis. Ordinary least squares regression fit and 95% confidence intervals are shown in red.

Periodicity

Figure 7 shows the wavelet power spectra for the six infections. Measles dynamics are erratic, exhibiting annual, biennial and multiennial cycles. Rubella (pre-vaccination), varicella and scarlet fever had a dominant annual periodicity. The annual periodicity of rubella disappeared after the introduction of vaccination. Mumps incidence exhibited a mixture of annual and multiennial cycles, particularly in the earlier years. The 4-year period was dominant for pertussis. This observed difference between pertussis and other similar childhood diseases such as measles has been described in detail in [8], and the longer cycles for pertussis are likely to be a result of differing response to dynamical noise as a result of a longer infectious period. Pertussis and rubella incidence exhibited significant non-resonant power, while varicella exhibited very large resonant power. This is consistent with previous work that has suggested that diseases with large R ₀ (such as measles and varicella) typically exhibit large resonant power (annual or multi-annual dynamics), while diseases such as pertussis (small amplitude of seasonality) and rubella (small R ₀) also exhibit non-resonant power due to transient dynamics [7].

Fig. 7. Wavelet power spectra of the incidence time-series for all states combined for (a) measles, (b) mumps, (c) rubella, (d) varicella, (e) scarlet fever, and (f) pertussis. Lighter shades of grey indicate greater power.

Synchrony

To characterize spatial patterns we examined the synchrony between states for each disease, which reflects correlations between states in both amplitude and timing of epidemics (Fig. 8). Epidemic synchrony declines with distance for measles, mumps, rubella but the decline is not statistically significant at the 95% level. For scarlet fever, the epidemic synchrony drops significantly with distance. Measles epidemics became less synchronized with the introduction of vaccination, a pattern that has also been documented for cities in the UK [35].

Fig. 8. Estimated spatial correlation functions showing epidemic synchrony (Pearson correlation coefficient) between states as a function of distance (as measured by the latitude and longitude of the centre of the state) for (a) measles, (b) mumps, (c) rubella, (d) varicella, (e) scarlet fever, and (f) pertussis, along with 95% confidence intervals. Red indicates pre-vaccination and blue indicates post-vaccination era epidemics (1991 onwards for measles and pertussis; 1998 onwards for mumps and rubella). The dashed horizontal line represents the average correlation across all states.

Rohani et al. [35] found a reverse pattern, i.e. greater synchrony following the introduction of vaccination, for pertussis epidemics in the UK. Interestingly, in Mexico, mumps and rubella epidemics became more synchronized in the vaccination era, while the synchrony for pertussis was not significantly different pre- and post-vaccination. Synchrony in varicella epidemics remains relatively high and constant with distance. This could be the result of the strong seasonality in varicella transmission (discussed below) synchronizing epidemics across the country.

The cross-correlation coefficients between all states and the Federal districts confirmed that epidemics were synchronized across Mexico for all diseases, except pertussis. The zero time-step lag was highest for all districts for measles and varicella, for 30/31 districts for mumps, and for 25/31 states for rubella and scarlet fever. For pertussis, the zero time step lag was the most common (for 6/31 states), followed by a lag of two months (for 5/31 states; epidemic occurring two months earlier in the Federal District).

Seasonal variation in transmission

Table 2 shows the results of fitting the TSIR model to the incidence data for each of the six infections, for all states combined. There was a great deal of variation in the estimated $\bar {S} $ for each infection, with the lowest estimated for pertussis ( $\bar {S} $  = 0·05) and the highest estimated for scarlet fever and varicella ( $\bar {S} $  = 0·8). The estimated reporting proportion was low for all six infections, particularly for scarlet fever and pertussis. A few studies, looking at measles reporting, have estimated the reporting proportion in Mexico to be as low as 3% during the 1980s [27]. Our estimated reporting proportion for some infections, such as pertussis, is quite a bit lower. Errors in our estimates of the reporting proportion and the mean proportion susceptible will affect the estimated magnitude of the seasonal transmission rates. This makes it difficult to compare the magnitude of the seasonal transmission rates across diseases and states, as the uncertainty in the estimated $\bar {S} $ and reporting proportion will vary. For scarlet fever and varicella, $\bar {S} $ was estimated to be 0·8 which is at the upper bound of our constraint on $\bar {S}. $ These estimates are most likely biased upwards, and has the effect of downward biasing the magnitude of the estimated seasonal transmission rates.

Table 2. Infection characteristics estimated by the TSIR model. The approximate generation time of each infection determines the discrete time step used in the TSIR model

TSIR, Time-series susceptible-infected-recovered; CI, confidence interval.

$\bar {S} $ is the estimated mean proportion of susceptibles in the population (95% CIs in parentheses, derived from the profile likelihood using the χ ² distribution with 1 degree of freedom). Mean reporting proportion is the mean of the time-varying reported rates estimated using locally varying regression of cumulative births on cumulative cases (standard deviation of the time-varying rates in parentheses). Mean β _s is the mean value of the seasonal transmission rates estimated via the TSIR model (maximum standard error in parentheses). The R ² values (on a log scale) for the TSIR fit are measured as 1 – [residual deviance/null deviance].

However, we are most interested in the shape of the seasonal transmission rate which is unaffected by the estimated $\bar {S} $ and reporting proportion. The transmission rates, estimated for all states combined, exhibit a marked decline during the summer for mumps, rubella, varicella, and scarlet fever (Supplementary Figs S8 and S9). We estimate a small decline in transmission for measles; however, the estimates for measles are highly uncertain as we only had a few years of incidence data prior to the start of vaccination. We did not detect any clear patterns in the transmission rate for pertussis.

We also estimated monthly transmission rates for all diseases separately by state (Fig. 9). For pertussis, for 16 out of the 32 states, we were unable to estimate at least one seasonal transmission parameter because there were not enough reported cases for that month. Furthermore, some estimated transmission rates for pertussis were unreasonably large, and we excluded all estimated transmission rates above 50 for our analysis (nine states). The yearly average transmission rate (averaged across all states) was 6·83 for measles, 10·58 for mumps, 4·20 for rubella, 2·19 for varicella, 3·14 for scarlet fever, and 14·54 for pertussis (Supplementary Fig. S11). The surprisingly low estimate of the average transmission rate for varicella may be due to biases in our estimate of $\bar {S}. $

All of the diseases, with the exception of pertussis, exhibited a decline in transmission during the summer months, which is consistent with school term-time forcing. School holidays in Mexico consists of 2 weeks at Christmas, 2 weeks at Easter, and summer holidays from July to mid-August [25]. Transmission (averaged across all states) peaked in December for measles, varicella and pertussis; February for rubella and scarlet fever; and in March for mumps. Transmission was lowest in August for measles, mumps, and rubella; July for varicella and scarlet fever; and June for pertussis. There was no significant spatial gradient in the magnitude of transmission or the coefficient of variation in the transmission parameter. For rubella and varicella, we found a negative significant association between the coefficient of variation and density (rubella: estimated association = −0·01, P = 0·011; varicella: estimated association = −0·01, P = 0·028). For rubella, we also found a significant correlation with poverty (estimated association = −0·18, P = 0·004), but not for any of the other infections.

Fig. 9. Mean centred seasonal transmission coefficients from the time-series susceptible-infected-recovered fit for (a) measles, (b) mumps, (c) rubella, (d) varicella, (e) scarlet fever, and (f) pertussis for each state (black lines). The red line indicates the seasonal transmission coefficient averaged over all states. Approximate school holidays are shown in grey.

Forecasting epidemics

Finally, although our main aim with the TSIR analysis was to examine seasonal patterns in transmission, we also explored the ability to predict long-term dynamics using the TSIR framework (Fig. 10). Our results, for all states combined, suggest that it may be possible to make reasonable predictions for varicella (NRMSE = 0·33) and scarlet fever (NRMSE = 0·37). The more complex dynamics of mumps (NRMSE = 0·44), rubella (NRMSE = 0·67), and pertussis (NRMSE = 1·38) are not captured as well by the TSIR framework; for measles the model fit is poor since we fit the TSIR model to only a few years of pre-vaccination incidence data (NRMSE = 3·40).

Fig. 10. Actual cases corrected for underreporting (in black) and 500 stochastic forward predictions (in red) from the time-series susceptible-infected-recovered fit for (a) measles, (b) mumps, (c) rubella, (d) varicella, (e) scarlet fever, and (f) pertussis for all states combined.

We also examined the ability to predict long-term dynamics of mumps, rubella, varicella, and scarlet fever at the state-level. The results varied widely across states. The NRMSE averaged across states was lowest for varicella (mean NRMSE = 0·75), followed by mumps (mean NRMSE = 0·92), scarlet fever (mean NRMSE = 1·01), rubella (mean NRMSE = 1·60), and measles (mean NRMSE = 4·47). The NRMSE was not significantly correlated with the state's demographic and socioeconomic indicators (P > 0·05 for all associations).

DISCUSSION

Our analysis reveals interesting epidemiological variations across the six infections, but our results are consistent with predictions from a seasonally forced SIR model. All infections, except pertussis, exhibited a seasonal pattern in transmission with a drop during the summer months. This is consistent with the well-defined term-time forcing of directly transmitted childhood infections observed elsewhere, most notably for pre-vaccination measles in the UK [2]. Unlike Metcalf et al. [24] who found increased transmission of bacterial infections (including pertussis and scarlet fever) during the summer, we found a reduction in transmission for scarlet fever; for the other bacterial infection (pertussis) we found no change in transmission during the summer.

The four pathogens – mumps, rubella, varicella, and scarlet fever – for which we had age-structured data, all had similar mean ages of infection (around 6 years). Basic compulsory education in Mexico starts at the age of 6, although almost 100% of children are enrolled in formal care and pre-school by age 5 [36]. Net enrolment rate in primary school was 94·3% for Mexico in 1985, and has remained fairly steady (lowest rate was 93·5% in 2006) [37]. On the other hand, fewer than 10% of children aged <3 years are enrolled in any form of formal care. Thus, the estimated mean age of infection, which is between the ages of 5 and 7 years, is consistent with the hypothesis of high transmission in schoolchildren during term time driving the epidemic dynamics.

Our results suggest that there was substantial variation in the epidemiology and dynamics of the six childhood infections in Mexico, and we also found evidence for the role of demography and socioeconomic factors in determining dynamics, at least for some of the infections. All six pathogens had peaks in incidence in the spring and troughs in autumn. Schenzle [6] showed, using an age-structured model, that for countries where the school year starts after a long summer vacation, with relatively short interruptions throughout the year, peak measles incidence occurs around March, and then declines throughout the school year as the susceptible populations within schools are depleted. The incidence time-series for childhood infections in Mexico are broadly consistent with this result. This also aligns with our estimates of seasonal patterns of transmission, which are consistent with term-time forcing (see above).

Mumps and pertussis had the most complex dynamics with multiennial cycles. This could be the result of longer generation times interacting with the nonlinear dynamics of transmission, as has been suggested earlier [8]. Rubella and scarlet fever, which have shorter generation times, had annual epidemic cycles. Although the generation time for varicella is similar to that of mumps and pertussis, it has annual periodicity possibly because of higher transmission rates.

Our analysis also reveals interesting spatial patterns. We document lower mean ages of infection for all pathogens around Mexico City, and higher mean ages of infection in the states in the eastern coast. This is consistent with past findings for rubella epidemics in Mexico [25]. States with higher poverty rates also had a higher mean age of infection for mumps, rubella, and varicella. Poorer states also experienced epidemics later in the year compared to other states. The observed spatial gradient could be due to biased estimates of the age-specific FOI [38]. Another possible explanation could be regional variability in R ₀, for instance due to differences in schooling rates and practices across states.

While the longitudinal gradient in the mean age of infection can be partly explained by the spatial gradient in poverty, further work is needed to disentangle whether or not the spatial gradient is driven by socioeconomic factors or perhaps by environmental factors. Environmental conditions, such as temperature and absolute humidity, may affect pathogen survival outside the host and the dispersal of infectious particles [20]. Scarlet fever epidemics in England and Wales in the mid-nineteenth century were significantly correlated with dry conditions in the spring/summer; this was hypothesized to be the result of changes in host susceptibility and/or better pathogen survival in the environment due to low absolute humidity [10]. All six pathogens are fomite transmitted, but unlike influenza [17], the link between environmental factors and diseases dynamics has not been shown conclusively for any of the pathogens considered here. Future work will require a more complete analysis of various drivers using formal techniques such as those described in [39].

Somewhat surprisingly, we found that states that are further east experienced measles, varicella, and scarlet fever epidemic peaks earlier in the year compared to other states. If the higher mean age of infection observed in these states was being driven by a lower transmission rate, then we would expect the epidemics to occur later in the year compared to other states. This could be due to biases in our estimation in the mean age of infection and corresponding R ₀. However, further work is needed to disentangle biases in the data and possible socioeconomic and environmental drivers.

We found that epidemics were synchronized between states for all diseases, with the exception of pertussis. Our scale of observation was too large to detect potential lags in remote areas of the country. Unlike the other infections, pertussis epidemics were not synchronized across the country. The irregular dynamics of pertussis incidence, along with the lack of epidemic synchrony between states, suggests that pertussis epidemics in Mexico were strongly influenced by stochasticity, which is consistent with past findings [8].

Interestingly, despite differences in the dynamics and epidemiology of the infections, five of the six childhood infections had a similar seasonal pattern in transmission rates. The peak transmission varied from December to March for the six infections, but lowest transmission was consistently in July or August (for all infections, except pertussis). For two of the infections, rubella and varicella, we found that states with lower poverty and lower population density had higher variation in the transmission rate through the year. It is possible that term-time forcing is more dominant in those states because of greater proportions of children in school (for the states with less poverty) or because most of the transmission is occurring in schoolchildren (for less dense states). The results from the TSIR simulations suggest that the dynamics of varicella and to some extent, scarlet fever, can be predicted with reasonable accuracy using stochastic simulations. At the state-level the quality of predictions varied greatly, but was not significantly associated with the states’ demographic and socioeconomic indicators.

In conclusion, our results suggest that despite differences in dynamics and epidemiology, five of the six childhood infections studied here, behave in accordance with a seasonally forced SIR model. Term-time forcing is likely to be a major component of seasonal dynamics for measles, mumps, rubella, varicella and scarlet fever in Mexico. Pertussis shows more unpredictable dynamics, possibly due to life history differences such as waning immunity and a longer generation time. While there was a longitudinal gradient in the mean age of infection, the pattern and magnitude of seasonal transmission rates were uncorrelated with longitude. It is possible that our scale of observation (state-level) is too large to capture more erratic local dynamics. Nonetheless, these results suggest that epidemics were fairly synchronized across Mexico for five of the six infections. This has important implications for control efforts and the predictability of epidemics. For infections, such as pertussis, that show less clear patterns of term-time forcing, predicting dynamics becomes difficult. However, for most directly transmitted childhood infections, with varying degrees of transmissibility, generation times, and immunity conferred, term-time forcing appears to be a dominant driver of disease dynamics, leading to more predictable epidemic cycles.