Epidemic model

The epidemic model (Supplementary Figure S1) followed the age-structured susceptible-exposed-infectious-removed (SEIR) construct [Reference Yang, Li and Shaman8, Reference Yang12], per Equation (1)):

(1) $$ \left\{\begin{array}{l}\frac{d{S}_i}{dt}=-{\lambda}_i(t)\frac{S_i^{m_1^i(t)}}{N_i}+b(t)\left(1-c(t)\right)N{\mathbf{A}}_{\mathbf{i}=\mathbf{1}}+\frac{M}{p(t)}{\mathbf{A}}_{\mathbf{i}=\mathbf{1}}-{v}^r(t){\delta}_{i-1}{\mathbf{A}}_{\mathbf{i}=\mathbf{2}}{S}_{i-1}\\ {}\hskip4.12em -{v}_i^c(t)+{\delta}_{i-1}{S}_{i-1}-{\delta}_i{S}_i-{\mu}_i(t){S}_i+{r}_i(t){f}_i^S(t)\\ {}\frac{d{E}_i}{dt}={\lambda}_i(t)\frac{S_i^{m_1^i(t)}}{N_i}-\alpha {E}_i+{\delta}_{i-1}{E}_{i-1}-{\delta}_i{E}_i-{\mu}_i(t){E}_i+{r}_i(t){f}_i^E(t)+{h}_i(t)\\ {}\frac{d{I}_i}{dt}=\alpha {E}_i-\gamma {I}_i+{\delta}_{i-1}{I}_{i-1}-{\delta}_i{I}_i-{\mu}_i(t){I}_i+{r}_i(t){f}_i^I(t)\\ {}\frac{d{R}_i}{dt}=\gamma {I}_i+{v}^r(t){\delta}_{i-1}{S}_{i-1}{\mathbf{A}}_{\mathbf{i}=\mathbf{2}}+{\delta}_{i-1}{R}_{i-1}-{\delta}_i{R}_i-{\mu}_i(t){R}_i+{v}_i^c(t)+{r}_i(t){f}_i^R(t)\\ {}\frac{dM}{dt}=b(t)c(t)N-\frac{M}{p(t)}-{\mu}_1(t)M\end{array}\right. $$

where $ t $ is time in days. $ i $ ( $ i=1,2,3,4 $ ) represents 4 age groups: <1, 1–14, 15–50, and > 50 years. For Group $ i $ , $ {N}_i $ is the population size, which is divided into those who are susceptible, $ {S}_i $ ; exposed (latently infected but not yet infectious), $ {E}_i $ ; infectious, $ {I}_i $ ; or recovered or immunised, $ {R}_i $ . M represents infants with maternal immunity. $ \alpha $ is the rate of progression from latent infection to infectiousness and $ \gamma $ is the recovery rate. $ {\delta}_i $ and $ {\mu}_i $ are the aging and natural death rates ( $ {\delta}_0 $ and $ {\delta}_4 $ were set to 0).

The force of infection, $ {\lambda}_i(t) $ , is given by

(2) $$ {\displaystyle \begin{array}{c}{\lambda}_i(t)=\sum \limits_{j=1}^4{\beta}_{ij}(t){I}_j^{m_2^j(t)}.\end{array}} $$

where $ {m}_2^i(t) $ is the mixing ratio of infectious individuals in Group i; $ {m}_1^i(t) $ is the corresponding mixing ratio of susceptible individuals. These exponents represent the extent of inhomogeneous mixing, with the value 1 representing homogeneous mixing. The exponents are further scaled by $ {k}_i(t) $ to test hypotheses related to differential mixing among migrants and infants (see Sections 4.4 and 4.6 of the Supplementary Materials), per

(3) $$ {\displaystyle \begin{array}{c}{m}_1^i(t)={m}_1{k}_i(t);\\ {}\hskip-0.52em {m}_2^i(t)={m}_2{k}_i(t)\end{array}} $$

The transmission rate, $ {\beta}_{ij}(t) $ , accounts for seasonality and school term-time mixing per

(4) $$ {\beta}_{ij}(t)=\left\{\begin{array}{c}\begin{array}{ll}{\beta}_{ij}\left[1+{b}_{season}\mathit{\cos}\left(\frac{2\pi }{365}\left(t-\phi \right)\right)\right]\frac{1+{b}_1 Term(t)}{b_{term}}& \mathrm{if}\;\left(i,j\right)=\left(2,2\right)\\ {}{\beta}_{ij}\left[1+{b}_{season}\mathit{\cos}\left(\frac{2\pi }{365}\left(t-\phi \right)\right)\right]& \mathrm{otherwise}\end{array}\end{array}\right. $$

where $ {\beta}_{ij} $ is the mean transmission rate from Group $ j $ to Group $ i $ ; $ {b}_{season} $ is the amplitude of sinusoidal forcing, and $ \phi $ represents the day of a year when the sinusoidal forcing reaches the maximum. As in previous work [Reference Keeling and Rohani13], for school-age children (here Group 2), we included a school term-time forcing $ \frac{1+{b}_1 Term(t)}{b_{term}} $ to model the impact of congregation during school terms. $ {b}_1 $ is the amplitude of the school term forcing, and $ Term(t) $ is given by

(5) $$ {\displaystyle \begin{array}{c} Term(t)=\left\{\begin{array}{ll}1& \mathrm{if}\;t\in {T}_{school\ terms}\\ {}-1& \mathrm{if}\;t\in {T}_{school\ break}\end{array}\right..\end{array}} $$

Based on the school schedule in China, $ {T}_{school\ break} $ includes a summer break that lasts for 7 weeks and ends on August 31st, and a winter break that lasts for 5 weeks and ends on the 14th day after the Lunar New Year (LNY). $ \frac{1+{b}_1 Term(t)}{b_{term}} $ averages 1 in a calendar year.

The epidemic model includes 4 age groups and thus needs 16 (=4 × 4) parameters for all $ {\beta}_{ij} $ ’s. To reduce the number of parameters, the $ \boldsymbol{\unicode{x03B2}} $ matrix is formulated using 6 parameters:

(6) $$ {\displaystyle \begin{array}{c}\boldsymbol{\beta} =\left[\begin{array}{cccc}{\beta}_1& {\beta}_1& {\beta}_5& {\beta}_4\\ {}{\beta}_1& {\beta}_2& {\beta}_6& {\beta}_4\\ {}{\beta}_5& {\beta}_6& {\beta}_3& {\beta}_4\\ {}{\beta}_4& {\beta}_4& {\beta}_4& {\beta}_4\end{array}\right].\end{array}} $$

Here, we assumed the transmission rate among infants (aged <1) equals that between infants and older children (aged 1–14), that is, $ {\beta}_1 $ . We also assumed all transmission rates associated with Group 4 are the same (i.e., $ {\beta}_4 $ ) because the oldest age group contributes less to measles dynamics and also has less social contact.

The basic reproductive number, $ {R}_0 $ , in this age-structured model is related to the $ \boldsymbol{\unicode{x03B2}} $ matrix per the next generation matrix method [Reference Keeling and Rohani13]:

(7) $$ {\displaystyle \begin{array}{c}{R}_0= eige{n}_{max}\left(\boldsymbol{\beta} \frac{1}{\gamma}\right)\end{array}} $$

where $ eige{n}_{max}\left(\right) $ is the function that gives the largest eigenvalue of the matrix. We reparametrised the model based on Equation (7) and included $ {R}_0 $ for estimation by setting $ {\beta}_1=1 $ and estimating the relative magnitude of $ {\beta}_k $ ( $ k=2,\dots, 6 $ ).

Relating to births and maternal immunity, $ b(t) $ is the birth rate per data, and $ c(t) $ is the proportion of infants with maternal immunity (set to $ \frac{R_3(t)}{N_3} $ , i.e., the immunity level of the age group of new mothers). $ p(t) $ is the mean duration of maternal immunity, and varied by hypothesis (Section 4.5 of the Supplementary Materials). $ {\mathbf{A}}_{i=1} $ and $ {\mathbf{A}}_{i=2} $ are indicator functions that were set to 1 for $ i=1 $ and $ i=2 $ , respectively; otherwise, these functions were set to 0. In Line 1 of Equation (1)), for infants (Group 1), $ b(t)\left(1-c(t)\right)N{\mathbf{A}}_{\mathbf{i}=\mathbf{1}} $ represents newborns without maternal immunity; $ \frac{M}{p(t)}{\mathbf{A}}_{\mathbf{i}=\mathbf{1}} $ represents infants who are initially protected by maternal immunity gradually becoming susceptible as this protection wanes.

For vaccination, we included both routine and catch-up vaccination. For simplicity, the routine vaccination in the model only includes 1 dose at age 1; this roughly captures the impact of the 2-dose measles vaccination programme, in which the first dose is administered at 8 months and the second at 18–24 months of age [Reference Ma, An, Hao, Cairns, Zhang, Ma, Cao, Wen, Xu, Liang, Yang and Luo6]. To account for vaccine effectiveness, we computed the immunization rate, $ {v}^r(t) $ , combining vaccination rates and effectiveness for both doses. In Lines 1 and 4 of Equation (1)), for 1- to 14-year-olds (Group 2), $ {v}^r(t){\delta}_{i-1}{\mathbf{A}}_{\mathbf{i}=\mathbf{2}}{S}_{i-1} $ represents the number of children turning age 1 (i.e., $ {\delta}_1{S}_1\Big) $ and getting immunised. $ {v}_i^c(t) $ represents the number of individuals acquiring immunity from catch-up vaccination (see Supplementary Materials for details).

We also accounted for migration. Specifically, $ {r}_i(t) $ represents the net number of migrants (i.e., the change in population size excluding births and deaths, due to a lack of detailed data) in Group $ i $ ; $ {f}_i^{\left[.\right]}(t) $ represents the fraction of migrants to Group i with a given disease status as specified by the superscript. The setting for these migration-related terms varied by hypothesis (Sections 4.1 and 4.3 of the Supplementary Materials). Lastly, $ {h}_i(t) $ models migration-related case importation (see specific settings in Section 4.2 of the Supplementary Materials).

Model inference (parameter estimation)

To estimate the model state variables ( $ {S}_i $ , $ {E}_i $ , $ {I}_i $ , $ {R}_i $ , $ M $ ) and parameters ( $ \alpha $ , $ \gamma $ , $ {m}_1 $ , $ {m}_2 $ , $ {b}_1 $ , $ {b}_{season} $ , $ \phi $ , $ \eta $ , $ {R}_0 $ , and $ {\beta}_2 $ to $ {\beta}_6 $ ), we ran the epidemic model in conjunction with a modified particle filter [Reference Yang, Karspeck and Shaman14]. Briefly, we initialised the model-filter system at the start of 1995, that is, before the large influx of migrants. The initial model state and parameter values (i.e., prior distributions) for the particles (i.e., model replica) came from our previous work [Reference Yang, Li and Shaman8] and the literature [Reference Mistry, Litvinova, Pastore, Chinazzi, Fumanelli, Gomes, Haque, Liu, Mu, Xiong, Halloran, Longini, Merler, Ajelli and Vespignani15] (see Supplementary Materials, Section 3). The system then estimated the model state sequentially from 1995 to 2014 via prediction-update cycles. In the prediction stage, the system integrates the particles forward stochastically with a daily time-step according to the epidemic model (Equation (1)); this generates the prediction, or the prior, for the next time step. In the update stage, the system computes the likelihood for each particle, and combines it with the prior to compute the posterior per Bayes’ rule.

To compute the likelihood, we modelled the reported incidence in Group $ i $ during reporting period $ t-\varDelta t $ to $ t $ ( $ \varDelta t $ is 1 year during 1995–2004 and 1 week during 2005–2014 based on data availability), Y_i,t, using a Gaussian distribution:

(8) $$ {\displaystyle \begin{array}{c}{Y}_{i,t}\sim N\left(\eta {C}_{i,t},{\sigma}_{i,t}^2\right).\end{array}} $$

where $ {C}_{i,t} $ is the corresponding model-estimated measles cases (including unreported cases) and $ \eta $ is the reporting rate (estimated by the filter). Following our previous work for Beijing [Reference Yang, Li and Shaman8], we assumed the reporting rate increased linearly before 2005 to reflect the improvement in disease surveillance. We set the annual increase to 0.0037, corresponding to a 20% increase from 1951 to 2004. $ {\sigma}_{i,t}^2 $ is the observational error variance. During 1995–2004, only yearly incidence combining all ages (Y_t) was available. Accordingly, we aggregated model estimated cases across all ages and, given the large year-to-year variation, simply scaled the observational error variance $ {\sigma}_t^2 $ to the reported incidence for the same year in addition to an arbitrary baseline of 100, per

(9) $$ {\displaystyle \begin{array}{c}{\sigma}_t^2=100+{Y}_t\times 3.\end{array}} $$

For 2005–2014, weekly incidence data were available for the 4 age groups. We used these data for inference and computed the observational error variance as

(10) $$ {\displaystyle \begin{array}{c}{\sigma}_{i,t}^2=100+\frac{\sum \limits_{k=0}^4{Y}_{i,t-k}}{5}\times 3.\end{array}} $$

As a sensitivity analysis, we tested smaller observational errors, that is, replacing the baseline of 100 with 50 and scaling factor of 3 with 1.

During filtering, we applied space reprobing, a technique to explore the state space to a larger extent to prevent the particles from being trapped in a sub-state space (a challenge known as particle impoverishment) [Reference Yang and Shaman16]. We reprobed the parameter space during 1995–2004 when only yearly data were available, and during 2005 when the filter switched to using weekly data. In addition, the net migration estimates used in the model likely underestimated the influx of susceptible migrants, particularly during 2011–2014 when the estimates were low (Supplementary Figure S2); as such, we also reprobed $ {S}_3 $ (i.e., susceptibles aged 15–50, the age group where most migrants were in) during the first four weeks after the LNY (i.e., the likely timing of worker migration) in 2011–2014.

Rationale and description

(5) Shorter duration of maternal immunity

Recent serological studies have showed that maternal measles antibodies wane earlier in infants born to mothers who acquired immunity via vaccination than those via natural infections [Reference Leuridan, Hens, Hutse, Ieven, Aerts and Van Damme9]. In addition, a substantial proportion of measles cases occurred among infants (<1 year; on average 22% during 2005–2016; Table 2). We thus hypothesised that, after decades of mass vaccination, the duration of maternal immunity may have been shortened; this in turn may have rendered more infants susceptible prior to receiving their first vaccine dose at 8 months of age and contributed to measles persistence. To test this, we set $ p(t) $ in Equation (1), the duration of maternal immunity, to 180, 150, or 120 days (Imm180, Imm150, and Imm120).

Table 2. Age distribution of measles cases by year and cases combing all ages in Beijing from 2005 to 2019

^a During 2017–2019, measles cases by age groups were not available.

Comparison procedure

Given the large number of combinations across the six sets of hypotheses, we compared the aforementioned hypotheses in two steps that focussed on migrant-related and infant-related hypotheses, respectively. In Step 1, we tested all combinations of the four sets of migrant-related hypotheses (n = 81; Figure 2), while setting the two sets of infant-related hypotheses to the baseline scenarios, that is, Imm180 and InfantMixBase. In Step 2, setting the migrant-related hypotheses at the most plausible scenario identified from Step 1, we further tested all combinations of the two sets of infant-related hypotheses (n = 9; Figure 3).

Figure 2. Comparing migrant-related hypotheses based on medians of log-likelihood. We formulated four sets of migrant-related hypotheses (increased migrant mixing [MigMixBase/MigMixMedium/MigMixHigh], seeding intensity [BaseSeed/MediumSeed/StrongSeed], higher migrant susceptibility [MigSusBase/MigSusMedium/MigSusHigh], and timing of migrant influx [MigConst/MigLNY12/MigLNY4]), and tested all combinations of the four sets (n = 81, i.e., the number of cells in the Fig). For each model (or hypothesis combination), we conducted 100 model inference runs using a model-filter system, and show the median of the 100 log-likelihood estimates in red. For infant-related hypotheses, all models tested here assumed Imm180 and InfantMixBase as baseline scenarios. See Table 1 for a summary of all hypotheses.

Figure 3. Comparing infant-related hypotheses based on medians of log-likelihood. We formulated two sets of infant-related hypotheses (increased infant mixing [InfantMixBase/InfantMixMedium/InfantMixHigh] and duration of maternal immunity [Imm180/Imm150/Imm120]), and tested all combinations of the two sets (n = 9, i.e., the number of cells in the Fig). For each model (or hypothesis combination), we conducted 100 model inference runs using a model-filter system, and show the median of the 100 log-likelihood estimates in red. For migrant-related hypotheses, all models tested here assumed MigMixBase, BaseSeed, MigLNY4, and MigSusHigh, because the previous comparison step identified this combination as one of the most plausible scenarios. See Table 1 for a summary of all hypotheses.