Surveillance data and preprocessing

We obtained the weekly time series of influenza data from week 27, 2010 to week 26, 2019 in Xi'an city, which is a temperate climate city in Northwest China (Fig. 1a) with about 10.2 million residents. The influenza sentinel hospitals and network laboratory, established by the Chinese National Influenza Center, aim to conduct influenza surveillance by reporting weekly influenza-like illness (ILI) cases and detecting the sampled specimens. Xi'an has five influenza sentinel monitoring hospitals (Xi'an Children's Hospital, Xi'an Central Hospital, Xi'an No.1 Hospital, Xi'an No.4 Hospital, Xi'an No.12 Hospital) and one influenza network laboratory (laboratory biosafety II or above and qualified to identify influenza subtypes and virus isolation), dispersed through four densely populated districts. Each week, the sentinel hospitals would use an Internet-based platform to report number of ILI cases (defined as patients whose body temperature ⩾38 °C with either cough or sore throat), divided into five age groups (i.e. 0–4, 5–14, 15–24, 25–59 and ⩾60 years). Moreover, every sentinel hospital was required to provide influenza network laboratory with 20 specimens (5–15 specimens before 2014 [33, 34]) of nasopharyngeal or throat swab per month (from April to September), and same specimens per week (from October to next March). The specimens of ILIs were collected from patients within 3 days of the appearance of symptoms and have not received any antiviral treatment, and then were transported in viral transport media at 4 °C to influenza surveillance network laboratory to isolate and identify influenza virus and subtypes [34]. In our model, we define week 27 to next week 26 as the whole influenza season. For example, week 27, 2010 to week 26, 2011 was denoted as the 2010/11 influenza season.

Fig. 1. Time series of influenza surveillance data in Xi'an, China between 2010/11 and 2018/19 influenza seasons. (a) Map of China, marking Xi'an city (red rectangle in the left sub-figure), showing five influenza sentinel hospitals (red dot in the middle and right sub-figure). (b) The weekly time series of influenza infections. (c) The seasonal proportion of influenza subtypes in the positive influenza infections of ILIs.

We obtained the weekly influenza infections by processing ILI cases and their influenza-positive rate in different age groups (Fig. 1b). The number of reported ILI cases and influenza-positive rates s differ in five age groups in the nine influenza seasons (see Appendix Fig. S1). Denote c _i,t as the weekly reported ILI cases in i-th (i = 1, 2, …5) age group at time t, and θ _i as the influenza-positive rate in i-th (i = 1, 2, …5) age group during each influenza season (July to next June, about 1700–3700 specimens were tested in each influenza season) (see Appendix Fig. S1C). The total weekly influenza infections C _t were calculated as follows.

(1)$$C_t = \mathop \sum \limits_{i = 1}^5 c_{i, t} \times \theta _i.$$

Model construction

We proposed a compartmental model to describe the transmission of seasonal influenza with different influenza subtypes (Fig. 2a). The population was divided into six compartments: susceptible individuals (S _t), vaccinated individuals (V _t), latent infections but not yet infectious (E _t), infectious individuals with (I _s,t) and without (I _a,t) symptoms, and recovered individuals (R _t). The model is given as follows:

(2)$$\left\{{\matrix{ {\displaystyle{{dS_t} \over {dt}} = {-}\kappa_tS_t-\lambda_tS_t + \mu_tN_t + \varepsilon R_t + \varepsilon V_t-\nu_tS_t} \hfill \cr {\displaystyle{{dV_t} \over {dt}} = \kappa_tS_t-( {1-\rho } ) \lambda_tV_t-\varepsilon V_t-\nu_tV_t} \hfill \cr {\displaystyle{{dE_t} \over {dt}} = \lambda_tS_t + ( {1-\rho } ) \lambda_tV_t-\delta E_t-\nu_tE_t} \hfill \cr {\displaystyle{{dI_{s, t}} \over {dt}} = ( {1-p} ) \delta E_t-\gamma I_{s, t}-\nu_tI_{s, t}} \hfill \cr {\displaystyle{{dI_{a, t}} \over {dt}} = p\delta E_t-\gamma I_{a, t}-\nu_tI_{a, t}} \hfill \cr {\displaystyle{{dR_t} \over {dt}} = \gamma I_{a, t} + \gamma I_{s, t}-\varepsilon R_t-\nu_tR_t} \hfill \cr } } \right.$$

where λ _t is the force of influenza infection at time t (i.e. the probability of susceptible population being infected by asymptomatic or asymptomatic infected person) and it is given by

(3)$$\lambda _t = \varphi _{s, t}\displaystyle{{I_{s, t}} \over {N_t}} + \varphi _{a, t}\displaystyle{{I_{a, t}} \over {N_t}}$$

Fig. 2. Flow chart of influenza transmission model and best model fit. (a) The total population is divided into six compartments (S _t, V _t, E _t, I _s,t, I _a,t and R _t denote the number of susceptible, vaccinated, exposed but not yet infectious, infectious with and without symptoms, and recovered individuals at time t, respectively). The force of infection for susceptible is denoted as λ _t, which involves influenza subtypes of different seasons. We assume that λ _t varies periodically as a sine function annually. More details are described in Materials and methods section. (b) Model calibration and data fitting based on weekly time series influenza infections between 2010/11 and 2018/19 influenza seasons. (c) Model validation based on the Pearson correlation coefficient of simulated and observed annual influenza infections.

Here φ _s,t, φ _a,t are the transmission rates when contacting with symptomatic and asymptomatic infections, respectively, which are assumed to change periodically in different influenza seasons based on the periodic characteristics of the influenza epidemic. We used sinusoidal function to express φ _s,t, φ _a,t as follows:

(4)$$\left\{{\matrix{ {\varphi_{s, t} = \phi_{{\rm season}} \times \left({1 + A \times \sin \left({\displaystyle{{2\pi t} \over {52}} + \omega } \right)} \right)} \hfill \cr {\varphi_{a, t} = q \times \varphi_{s, t}} \hfill \cr } } \right.$$

where ϕ _season represents the baseline transmission rate decided by the transmission rate of different influenza subtypes and their proportion in each influenza season, given by

(5)$$\phi _{{\rm season}} = \beta _{{\rm H}1{\rm N}1}\theta _{{\rm H}1{\rm N}1, } + \beta _{{\rm H}3{\rm N}2}\theta _{{\rm H}3{\rm N}2} + \beta _{{\rm type}\,{\rm B}}\theta _{{\rm type}\,{\rm B}}$$

Here β _H1N1, β _H3N2, β _type B are the transmission rates of H1N1, H3N2 and B type, respectively, while θ _H1N1, θ _H3N2 and θ _type B are the proportion of these subtypes in each influenza season (Fig. 1c). q is the relative transmissibility for asymptomatic infections compared with symptomatic infections. A (⩽1) is the seasonal amplitude, and ω is the phase shift in the sinusoidal function. p is the fraction of asymptomatic infection. μ _t denotes the natural birth rate and ν _t denotes the natural death rate. Denote 1/δ, 1/γ and 1/ɛ as the average period of latency, recovery and immunity protection, respectively. κ _t represents the vaccination rate of UIV for susceptible individuals at time t and ρ denotes the effectiveness of UIV. Our model did not consider the effects of seasonal influenza vaccine due to very low seasonal influenza vaccine coverage (~2% in China, <0.5% in Xi'an).

Parameter estimation and model calibration

We obtained the following four parameters from the published literatures. The mean incubation time for influenza is 2 (1–4) days (1/δ = 2) [35]. The fraction of asymptomatic infection was chosen as P = 40% (30–50%) [Reference Furuya-Kanamori36]. The natural birth rate (μ _t) and death rate (ν _t) were obtained from the Xi'an Bureau of Statistics [37] (see Appendix Table S2). The other parameters were estimated by fitting the model (equations (2)–(5)) to the weekly time series of influenza cases from week 27, 2010 to week 26, 2019 (Fig. 2b) using the Nonlinear Least Square (NLS) method, which minimised the residual error of model-simulated and reported influenza infections ((1 − p)δE _t − C _t)². These initial estimates are used as prior information for carrying out the Bayesian Markov Chain Monte Carlo (MCMC) procedure with an adaptive Metropolis-Hastings (M-H) algorithm and post estimates can be obtained [Reference Haario38]. The algorithm was run for 10 000 iterations with a burn-in (some iterations at the beginning of an MCMC run are thrown away) of 3000 iterations, and we used the rest 7000 iterations to derive the mean value and 95% confidence intervals (CI) of parameters as shown in Table 1. We conducted Latin Hypercube Sampling (LHS) 1000 times [Reference Helton and Davis39] in estimated parameters and its 95% CI to run the dynamic model procedure and obtain the 95% CI of influenza incidence cases.

Table 1. Prior information of estimated parameters based on references or assumptions and its' post estimate values using MCMC methods

^a The average population size served by each hospital is about 110 000 (total 100 hospitals in Xi'an city with population size 11 million), so we assumed the population size covered by 5 sentinel hospitals was about 550 000.

We calculated the Pearson correlation coefficient r and coefficient of determination R² [Reference Taylor40, Reference Du41] to evaluate the goodness of fit between model-simulated and reported influenza infections in each influenza season (Fig. 2c). All the procedures and analyses were implemented by MATLAB R 2019b.

Construction of scenarios

We considered two vaccination patterns about the vaccination time based on WHO's recommendations [32]. First, WHO recommends that the ideal time to implement the influenza vaccination is in fall, before the influenza season begins (November in the northern hemisphere), so we choose the time of vaccination as between September and October, i.e. 2-month pattern. Second, WHO also recommends that people could be vaccinated at any time during the influenza season to prevent more infections, so we choose another vaccination pattern that people can be vaccinated between September and the following February, i.e. 6-month pattern, which is closer to the current situation of real-world influenza vaccination.

Based on NIAID's influenza research programme [18] and the average effectiveness of the seasonal influenza vaccines across different subtypes [Reference Belongia7], we defined the UIV with a 75% effectiveness as high effectiveness and a 50% effectiveness as low effectiveness. The current vaccine coverage in China is as low as 1.5–2.2% [Reference Yang11] and we consider it would probably increase to 10% after UIV rollout. Thus, we defined a 10% coverage as low coverage and a 30% coverage as high coverage.

We projected the number of influenza infections in each influenza season under the following five scenarios (Fig. 3): (1) the no UIV scenario (baseline scenario); (2) the 2-month vaccination pattern with low UIV coverage rate (10%) and effectiveness (50%) scenario; (3) the 6-month vaccination pattern with low UIV coverage rate (10%) and effectiveness (50%) scenario; (4) the 2-month vaccination pattern with high UIV coverage rate (30%) and effectiveness (75%) scenario; (5) the 6-month vaccination pattern with high UIV coverage rate (30%) and effectiveness (75%) scenario. Previous study has reported the time window of antibody response is about 2 weeks after vaccination [42], so we assumed there was 2 weeks delay of immunity protection of vaccinated individuals.

Fig. 3. The simulated influenza infections for five constructed scenarios between 2010/11 and 2018/19 influenza season. The solid black line means influenza infections in the no vaccine scenario. The dotted red line means influenza infections in the 6-month vaccination pattern with low UIV coverage rate and effectiveness scenario. The solid red line means influenza infections in the 2-month vaccination pattern with low UIV coverage rate and effectiveness scenario. The dotted blue line means influenza infections in the 6-month vaccination pattern with high UIV rate and effectiveness scenario. The solid blue line means influenza infections in the 2-month vaccination with high UIV coverage rate and effectiveness scenario.

Sensitivity analysis

We expended the range of UIV coverage rate (0–50%) and effectiveness (0–75%) to perform sensitivity analysis (Fig. 4). Under the 2-month vaccination pattern, we calculated the percentage of averted influenza infections with varied UIV effectiveness and coverage in each influenza season, and plotted them as a function of UIV effectiveness and coverage. We also gave a special example that defined 50% of averted influenza infections as a threshold to evaluate the effects of varied UIV effectiveness and coverage. A similar plot for the 6-month vaccination pattern was shown in Appendix Fig. S2.

Fig. 4. Contour plots about the percentage of averted infections as a function of universal influenza vaccine (UIV) coverage and effectiveness with 2-month vaccination pattern from 2010/11 to 2018/19 influenza seasons. The solid black isoclines indicate the threshold that the percentage of averted infections is 50% which is just for show. The dashed black lines correspond to the minimal vaccine effectiveness and vaccine coverage rate when the percentage of averted infections is 50%.