Skip to main content
Download PDF

Data-driven interdisciplinary mathematical modelling quantitatively unveils competition dynamics of co-circulating influenza strains

Journal of Translational Medicine volume 15, Article number: 163 (2017) Cite this article

Abstract

Background

Co-circulation of influenza strains is common to seasonal epidemics and pandemic emergence. Competition was considered involved in the vicissitudes of co-circulating influenza strains but never quantitatively studied at the human population level. The main purpose of the study was to explore the competition dynamics of co-circulating influenza strains in a quantitative way.

Methods

We constructed a heterogeneous dynamic transmission model and ran the model to fit the weekly A/H1N1 influenza virus isolation rate through an influenza season. The construction process started on the 2007–2008 single-clade influenza season and, with the contribution from the clade-based A/H1N1 epidemiological curves, advanced to the 2008–2009 two-clade influenza season. Pearson method was used to estimate the correlation coefficient between the simulated epidemic curve and the observed weekly A/H1N1 influenza virus isolation rate curve.

Results

The model found the potentially best-fit simulation with correlation coefficient up to 96% and all the successful simulations converging to the best-fit. The annual effective reproductive number of each co-circulating influenza strain was estimated. We found that, during the 2008–2009 influenza season, the annual effective reproductive number of the succeeding A/H1N1 clade 2B-2, carrying H275Y mutation in the neuraminidase, was estimated around 1.65. As to the preceding A/H1N1 clade 2C-2, the annual effective reproductive number would originally be equivalent to 1.65 but finally took on around 0.75 after the emergence of clade 2B-2. The model reported that clade 2B-2 outcompeted for the 2008–2009 influenza season mainly because clade 2C-2 suffered from a reduction of transmission fitness of around 71% on encountering the former.

Conclusions

We conclude that interdisciplinary data-driven mathematical modelling could bring to light the transmission dynamics of the A/H1N1 H275Y strains during the 2007–2009 influenza seasons worldwide and may inspire us to tackle the continually emerging drug-resistant A/H1N1pdm09 strains. Furthermore, we provide a prospective approach through mathematical modelling to solving a seemingly unintelligible problem at the human population level and look forward to its application at molecular level through bridging the resolution capacities of related disciplines.

Background

Over the past 100 years, the three major influenza pandemics, namely the 1918 pandemic, 1957 pandemic, and 1968 pandemic, claimed millions of human lives [1, 2]. A number of authors used simulation models to explore the optimal strategies of containing a potential influenza threat [3,4,5,6,7,8]. The results not only provided insights into the control strategies of influenza pandemics but also shed light on the resolution of some enigmas of the seasonal influenza epidemics. One common interest pertinent to seasonal epidemics and pandemic emergence concerns the competition dynamics of co-circulating strains.

Given the ever worldwide dissemination of drug-resistant seasonal A/H1N1 H275Y variants during the 2007–2009 influenza seasons, the detection of drug-resistant A/H1N1pdm09 variants was of growing concern [9,10,11,12,13,14,15,16,17]. On one hand, surveillance systems tried to track and analyze the evolution and epidemiology of influenza viruses for better vaccine strain selection and timely antiviral susceptibility assessment [18]. On the other hand, nevertheless, the transmission dynamics of drug-resistant seasonal A/H1N1 H275Y variants was still mysterious at the human population level.

According to Yang et al. [19], A/H1N1 H275Y variants were seldom isolated in Taiwan until September 2008 when they were found emerging in 14.3% of A/H1N1 isolates. As A/H1N1 subtype progressively prevailed, the proportion of A/H1N1 H275Y strains ultimately reached 100% by the end of year 2008. The vast majority of A/H1N1 H275Y strains belonged to clade 2B-2, which dominated the 2008–2009 influenza season in Taiwan. Worthy of note, however, were the preceding strains belonging to A/H1N1 clade 2C-2, which emerged in June 2008 but disappeared before the peak of the 2008–2009 influenza season.

We therefore constructed a transmission model to simulate the competition dynamics of co-circulating clades 2C-2 and 2B-2 during the 2008–2009 influenza season in Taiwan. Virological factors, epidemiological observations, social factors, and climatic factors were incorporated into the model. Annual influenza vaccination program was set accordingly to simulate the real situation. As to antivirals, oseltamivir was the major antiviral stockpile then but was documented as rarely prescribed before May 2009 in Taiwan [19, 20]. Taking the advantage of the antiviral-free scenario, we studied the competition dynamics between the preceding clade 2C-2 and the succeeding clade 2B-2 in a quantitative way.

Methods

Virological data

We intended to develop a heterogeneous dynamic transmission model and ran the model to fit the weekly A/H1N1 influenza virus isolation rate through an influenza season. We inspected the laboratory data of the national influenza surveillance system coordinated by Taiwan Centers for Disease Control from week 1 (starting on January 1) of 2007 through week 26 (ending on June 28) of 2009. The surveillance system covered approximately 75% of the 352 basic administrative units throughout the northern, central, southern, and eastern regions of Taiwan [19, 21, 22]. Considering the inherent variability of the weekly data, we adopted the aggregate weekly A/H1N1 influenza virus isolation rate (IVIR) of each week, computed from all positive isolations divided by all collected specimens of the three consecutive weeks, for the model to simulate (see Additional file 1: Table S1). To accommodate to the time frame of the model, i.e., 30 days a month and 12 months a year, we further reframed the time frame of the aggregate weekly data so that 52 weeks a year were transformed into 360 days a year. April 1 was arbitrarily defined as day 1 of the model.

Epidemiological evidence

As shown in Yang et al. [19], the 2007–2008 seasonal influenza epidemic, dominated by A/H1N1 clade 2B-1, peaked in January of 2008, and the 2008–2009 seasonal influenza epidemic, dominated by clade 2B-2, peaked in February of 2009. Besides, an earlier hillock of clade 2C-2 occurred during the period from September through October of 2008. To simulate the weekly A/H1N1 IVIR curve, the model hereby assumed the peaks of the 2007–2008 and 2008–2009 seasonal influenza epidemics were day 267 and day 304, respectively. For the 2008–2009 epidemic, as evidenced in Yang et al. [19], the model further assumed day 180 as the hillock peak, with a height between 7/192 and 10/192 the height of the epidemic peak, and assumed the crossover occurred on day 210.

Demographic data and age-specific contact rates

The model population adopted the mid-year population structure of the year, 2007 and 2008 respectively, and was classified into six age groups: 0–5, 6–12, 13–19, 20–39, 40–59, and ≥60 years (see Additional file 2: Table S2). The mid-year age-structured population sizes were obtained from Department of Household Registration, Ministry of the Interior (http://www.ris.gov.tw/zh_TW/346). Births and deaths were neglected because of the relatively short time scale that the simulation spanned. The model further adopted the normalized age-specific contact rates estimated in Wallinga et al. [23]. See Additional file 3: Appendix, for details.

Seasonality modulated by absolute humidity

According to Shaman et al. [24, 25], under the assumption of influenza transmission seasonality being driven by absolute humidity, the following equation was adopted to modulate influenza transmissibility in the model:

$$R_{\text{e}} = R_{{ 0 {\text{min}}}} + e^{{\left( { - 1 8 0\,\times\,{\text{SH}}\,+\,{\text{log }}\left( {R_{{ 0 {\text{max}}}}\,-\,R_{{ 0 {\text{min}}}} } \right)} \right)}}$$
(1)

where R e, effective reproductive number, defined as the actual average number of secondary cases per primary case during a time period, was adopted to estimate the transmission rate coefficient during the period. R 0max and R 0min represented the maximum and minimum daily basic reproductive numbers, respectively, and SH represented specific humidity and acted as proxy of absolute humidity. The meteorological variables used to derive specific humidity could be referred to Taiwan Central Weather Bureau (http://www.cwb.gov.tw/V7/climate/monthlyData/mD.htm). The relevant data from April 2007 through June 2009 were shown in Additional file 4: Table S3. See Additional file 3: Appendix, for details.

Transmission model construction

We proceeded to construct a transmission model with four dimensions taken into consideration, i.e., virological factors, epidemiological observations, social factors, and climatic factors [26]. We started on constructing an susceptible-exposed-infective-recovered (SEIR) compartmental model as the backbone [27]. The equations below described the deterministic behaviors of the model:

$$\frac{{{\text{d}}S\left[ i \right]\left( t \right)}}{{{\text{d}}t}} = - \lambda \left[ i \right]S\left[ i \right]\left( t \right)$$
(2)
$$\frac{{{\text{d}}E\left[ i \right]\left( t \right)}}{{{\text{d}}t}} = \lambda \left[ i \right]S\left[ i \right]\left( t \right) - \theta E\left[ i \right]\left( t \right)$$
(3)
$$\frac{{{\text{d}}I\left[ i \right]\left( t \right)}}{{{\text{d}}t}} = \theta E\left[ i \right]\left( t \right) - \alpha I\left[ i \right]\left( t \right)$$
(4)
$$\frac{{{\text{d}}R\left[ i \right]\left( t \right)}}{{{\text{d}}t}} = \alpha I\left[ i \right]\left( t \right)$$
(5)

where S[i], E[i], I[i], and R[i] represented the size of age group i in compartment S, E, I, and R, respectively (see Table 1 for details of the parameters). The relationship between R e and the λ[i] and the interaction between various age groups were illustrated in Appendix (see Additional file 3).

Table 1 Summary of model parameter values
Full size table

The criteria of a successful simulation were defined correspondingly for single-clade epidemic simulation of the 2007–2008 influenza season and for two-clade epidemic simulation of the 2008–2009 influenza season (see the “Results” section for details). The simulation spanned 15 months, covering the length of an annual influenza season and extending 3 months forward so that the observed emerging stage of the seasonal influenza epidemic could be included (see Additional file 3: Appendix, for details).

However, viewing the exploration that it would be impossible to reproduce the scenario of the 2008–2009 influenza season unless the preceding clade incurred reduction in transmissibility on encountering the succeeding clade, the model further tried probing into strain fitness in two ways. On one hand, since we were interested in whether the succeeding clade per se had fitness advantage over the preceding one, we defined relative fitness (RF) to quantify the ratio of transmissibility, in terms of R e, of the former over the latter as the expression below.

$$RF = R_{{{\text{e}}\;{\text{succeeding}}\;{\text{clade}}}} /R_{{{\text{e}}\;{\text{preceding}}\;{\text{clade}}}}$$
(6)

On the other hand, because we were also interested in whether the preceding clade suffered from transmissibility impairment when encountering the succeeding one, we defined competition cost (CC) to quantify the possible reduction in transmissibility of the former associated with the interaction between the two clades as the expression below.

$$CC = 1{ - }(R_{{{\text{e}}\;{\text{preceding}}\;{\text{clade/succeeding}}\;{\text{clade}}}} /R_{{{\text{e}}\;{\text{preceding}}\;{\text{clade}}}} )$$
(7)

Vaccination against circulating strains

In our model, clade X, Y, and Z represented A/H1N1 clade 2B-1, 2C-2, and 2B-2, respectively. All three clades were predominantly reported as low reactors to the contemporaneous annual vaccine strains [30]. We therefore assumed 50% of the vaccinees were protected from infection of the contemporaneous seasonal strains [28, 29] and acquired full immunity 15 days after vaccination [31]. For simplicity, full cross-immunity was further assumed for clades Y and Z, and people infected by either clade entered into compartment R without differentiation. We assumed that all those receiving vaccination were originally in compartment S rather than compartment E, I, or R. Hence, our model moved the vaccinated of each age group in the compartment S, adjusted by vaccine efficacy, directly into compartment R on a daily basis 15 days after vaccination [32, 33]. See Additional file 5: Table S4, for the progress of the annual influenza vaccination program.

Statistical analysis

The size of compartment I, under the assumption that the infectives could reflect the influenza virus isolation rate, was set to track the virological data as well as epidemiological evidence in pursuit of high correlation. Pearson method was used to estimate the correlation coefficient between the simulated epidemic curve and the observed weekly A/H1N1 influenza virus isolation rate curve.

Programming language and software environment

All mathematical modelling, simulations, calculations, graphical outputs and statistical analyses were performed with R version 3.0.1(https://www.r-project.org).

Results

Model

We constructed a heterogeneous dynamic transmission model [27]. The core of the model is schematically represented in Fig. 1.

Fig. 1
figure 1

The SEIR model. Individuals in the susceptible compartment (S) progress along the exposed compartment (E), infective compartment (I), and recovered compartment (R) at per capita rates λ, θ, and α, respectively, once they are infected with influenza. The subscript, W, W1, W2, and V, specifies the transition process involving each respective clade/strain. Individuals pass from compartment S directly into compartment R 15 days after receiving influenza vaccination at a proportion of ν according to the program agenda (see Additional file 5: Table S4). a Model for single-clade influenza season. b Model for two-clade influenza season. The model applied to each age group as well as the whole population

Full size image

Fitting single-clade influenza season and estimating effective reproductive numbers

We first ran the model under the context of the 2007–2008 influenza season and introduced clade X on April 1, 2007. A successful simulation was defined as follows: (a) the correlation coefficient of at least 95% to the observed weekly A/H1N1 IVIR curve and (b) the epidemic peak at most 3 days away from the one of the weekly A/H1N1 IVIR curve (day 267) (see “Methods” for details). Through fitting the simulation to the weekly A/H1N1 IVIR curve, the model found not only successful simulations but also the potentially best one. The ten best-fit simulations, ranked by the correlation coefficients, were shown in Table 2 (see Table 1 for details of the parameters). For the first parameter set, the model generated an epidemic curve of clade X that peaked on day 269 and was positively correlated to the observed weekly A/H1N1 IVIR curve with the correlation coefficient up to 96.83%. The monthly R e’s of the 15 months spanned by the simulation, corresponding to the months from April 2007 through June 2008, were estimated around 1.33, 1.32, 1.31, 1.31, 1.31, 1.31, 1.32, 1.33, 1.34, 1.35, 1.36, 1.34, 1.32, 1.32, and 1.31 (see Additional file 6: Table S5). The annual R e of clade X was hereby estimated around 1.32. The model therefore could successfully identify the best combination of R 0max and R 0min and evaluate the R e of the single-clade influenza season. In other words, in terms of the R e estimated, the first parameter set along with the actual data embedded in the model could establish an epidemic 96.83% correlated to the real epidemic of the 2007–2008 influenza season dominated by A/H1N1 clade 2B-1.

To illustrate the effect of combinations of parameter values on the simulation curve and the sensitivity of the model to the parameter values, we selectively chose R 0max, the maximum daily basic reproductive number [25], at 1.45, 1.55, and 1.65, and R 0min, the minimum daily basic reproductive number [25], at 1.20, 1.30, and 1.40, and ran the model generating the epidemic curves for the nine combinations (Fig. 2). The weekly A/H1N1 IVIR curve, serving as the reference, was scaled by the best-fit simulation (set 1 in Table 2) for visual presentation. Obviously, higher R 0max or R 0min shifted the epidemic curve left with higher peak and lower R 0max or R 0min otherwise. Besides, R 0min dominated the main course of the epidemic curve and R 0max further shaped the simulation. All the epidemic curves converged to the best-fit simulation, and the ten best-fit that the model reported were shown earlier in Table 2.

Fig. 2
figure 2

2007–2008 seasonal influenza A/H1N1 epidemic simulations. The model generated a simulated epidemic curve for each combination of R 0max and R 0min. The combinations 1 through 9 presented in the figure, in the form of (R 0max, R 0min), were (1.65, 1.40), (1.65, 1.20), (1.65, 1.30), (1.45, 1.40), (1.45, 1.20), (1.45, 1.30), (1.55, 1.40), (1.55, 1.20), and (1.55, 1.30). The combination 9 simulated the observed weekly A/H1N1 IVIR curve (A/H1N1, shown in black) best and was 96.83% positively correlated. Against the time horizon, the proportion of the total population being infective was plotted on the left vertical axis and the weekly A/H1N1 IVIR on the right

Full size image
Table 2 Parameter sets with corresponding characteristics of the ten best-fit simulations for the 2007–2008 influenza season
Full size table

Fitting two-clade influenza season and estimating effective reproductive numbers

We continued to run the model under the context of the 2008–2009 influenza season and introduced clade Y on June 1 and clade Z on September 1, 2007. A successful simulation was defined by the following five criteria: (1) the correlation coefficient of at least 95% to the observed weekly A/H1N1 IVIR curve, (2) the peak of clade Y at most 15 days away from day 180, (3) the peak of clade Z at most 3 days away from the peak of the weekly A/H1N1 IVIR curve (day 304), (4) the height ratio of the peaks of clade Z to Y in the range between 192:10 and 192:7, and (5) the curve of clade Y crossing over the one of clade Z on the day at most 3 days away from day 210 (see “Methods” for details). However, the model reported that no successful simulation returned until strain fitness was taken into consideration. The ten best-fit simulations, ranked by the correlation coefficients, were shown in Table 3 (see Table 1 for details of the parameters). Through fitting the simulation to the weekly A/H1N1 IVIR curve as well as the clade-based A/H1N1 epidemiological curves, the model found not only successful simulations but also the potentially best one.

Table 3 Parameter sets with corresponding characteristics of the ten best-fit simulations for the 2008–2009 influenza season
Full size table

For the first parameter set, the model generated an epidemic curve that represented the summation for clades Y and Z and was 96.96% positively correlated to the observed weekly A/H1N1 IVIR curve. The curve of clade Y peaked on day 167 and declined thereafter. The curve of clade Z peaked on day 307 with a height 21.93 times that of clade Y. The two curves crossed on day 210. For clade Z, the monthly R e’s of the 15 months spanned by the simulation, corresponding to the months from April 2008 through June 2009, were estimated around 1.65, 1.64, 1.63, 1.62, 1.62, 1.63, 1.63, 1.66, 1.69, 1.70, 1.66, 1.67, 1.66, 1.65, and 1.63 with the annual R e of around 1.65 (see Additional file 7: Table S6). As to clade Y, the annual R e would originally be equivalent to 1.65 but finally took on around 0.75 after the introduction of clade Z (see Table 3 and Additional file 7: Table S6). The model reported that the strain fitness of clade Z was equivalent to that of clade Y (RF Z  = 1.00) and the latter encountered a reduction of strain fitness of around 71% (CC Y  = 0.71) while competing with the former. The model therefore could successfully identify the combination of R 0max, R 0min, RF Z and CC Y that best simulated the 2008–2009 influenza season and evaluate the R e’s of the co-circulating clades. In other words, in terms of the R e estimated for each respective clade, the first parameter set along with the actual data embedded in the model could establish an epidemic 96.83% correlated to the real epidemic of 2008–2009 influenza season dominated by clade 2B-2 along with ever preceding clade 2C-2.

Relative fitness declaring both clades per se equivalent in transmissibility

To further illustrate the interaction between clades Y and Z and the sensitivity of the model to the RF Z and CC Y values, we would elaborate the simulation results in two ways. First, we selectively chose relative fitness of clade Z (RF Z ) from 0.80 to 1.20 at intervals of 0.05 while fixing R 0max at 2.00, R 0min at 1.61, and competition cost of clade Y (CC Y ) at 0.71, with the weekly A/H1N1 IVIR curve scaled by the best-fit simulation (set 1 in Table 3) as the reference. A spectrum of epidemic curves thereby generated by the model demonstrated the effect of RF Z on the simulation curve (Fig. 3). We noticed that higher RF Z shifted the simulated curve to the left with a higher peak and lower RF Z to the right with the peak lower. In fact, strains of clade Z with relative fitness of 1.00 (Strain Z 1.00) along with the preceding strains of clade Y suffering from competition cost of 0.71 (Strain Y 0.29) initiated a simulated epidemic that correlated to the observed one represented by the weekly A/H1N1 IVIR curve best. All the epidemic curves converged to the best-fit simulation, and the ten best-fit that the model reported were shown earlier in Table 3.

Fig. 3
figure 3

The effect of relative fitness of clade Z on 2008–2009 seasonal influenza A/H1N1 epidemic simulation. As clade Y assumed competition cost of 0.71 (Strain Y 0.29), with R 0max and R 0min fixed at 2.00 and 1.61 respectively, the model generated epidemic curves of clade Z of a series of relative fitness from 0.80 (Strain Z 0.80) through 1.20 (Strain Z 1.20). Among them, the Strain Z 1.00 simulated the observed weekly A/H1N1 IVIR curve (A/H1N1, shown in black) best and was 96.96% positively correlated. Against the time horizon, the proportion of the total population being infective was plotted on the left vertical axis and the weekly A/H1N1 IVIR on the right

Full size image

Competition cost explaining the dynamics between both clades

Secondly, we instead selectively chose CC Y between 0.00 and 0.90 while fixing R 0max at 2.00, R 0min at 1.61, and RF Z at 1.00, a rather intriguing pattern of the vicissitudes of two co-circulating clades was then unfolded (Fig. 4). Facing co-circulating clade Z equivalently fit to itself, if clade Y had not suffered much from the competition cost, i.e., if the CC Y had been less than 0.20, it could have been the dominant clade during the 2008–2009 influenza season. However, once clade Y encountered a competition cost of more than 0.30, clade Z became the dominant one during the 2008–2009 influenza season. Moreover, strains of clade Y with competition cost of 0.71 (Strain Y 0.29) coupled with the emerging strains of clade Z with relative fitness of 1.00 (Strain Z 1.00) presented a simulated epidemic that correlated to the observed weekly A/H1N1 IVIR curve best. Again, all the epidemic curves converged to the best-fit simulation.

Fig. 4
figure 4

The effect of competition cost of clade Y on 2008–2009 seasonal influenza A/H1N1 epidemic simulation. With R 0max and R 0min fixed at 2.00 and 1.61 respectively, while clade Z assuming relative fitness of 1.00 (Strain Z 1.00), the model generated the dynamic behavior of clade Y with competition cost between 0.00 (Strain Y 1.00) and 0.90 (Strain Y 0.10). The Strain Y 0.29 coupled with Strain Z 1.00 simulated the observed weekly A/H1N1 IVIR curve (A/H1N1, shown in black) best and was 96.96% positively correlated. Against the time horizon, the proportion of the total population being infective was plotted on the left vertical axis and the weekly A/H1N1 IVIR on the right

Full size image

From the implementations stated above, the model reported Strain Y 0.29 with Strain Z 1.00 as the best combination of clades Y and Z co-circulating during the 2008–2009 influenza season. Through step-by-step exploration (see “Methods” for details), with R 0max set between 2.13 and 1.99, R 0min between 1.63 and 1.57, RF Z between 1.03 and 0.97, and CC Y between 0.60 and 0.78 in the final step, there were 617 successful simulations reported among the final 13,965 simulations. In the 617 successful simulations, we observed that RF Z was in the range of 0.99 and 1.02 and CC Y was in the range of 0.65 and 0.76. The findings strongly converged onto the scenario that the competition cost encountered by clade Y was involved in the interaction between the two clades.

Discussion

On the basis of the results, we address three points pertinent to the 2008–2009 influenza season in Taiwan about strain competition in particular. First, in terms of R e, clade 2B-2 would be as fit as clade 2C-2, as evidenced by RF Z of 1.00 derived from the transmission model. Although there were compelling evidences that the clinical effects of A/H1N1 H275Y variants, which circulated globally during the 2007–2009 influenza seasons, were consistent with normal seasonal activity [34, 35], the transmission efficiency of the variants, nevertheless, was not well studied but in one animal study [36]. In our study, the model not only reported that clade 2B-2, carrying H275Y mutation in the neuraminidase, transmitted as efficiently as clade 2C-2 but also revealed that the transmission efficiency of clade 2B-2 could not sufficiently explain the vicissitudes of the two co-circulating clades. The model demonstrated that clade 2B-2 outcompeted clade 2C-2 for the 2008–2009 influenza season mainly as a result of the competition cost that the latter incurred in the interaction between them.

Secondly, according to Hurt and Bloom et al. [37, 38], H275Y mutation was typically associated with poor replication and transmission, and permissive mutation could enable oseltamivir-resistant A/H1N1 variants to restore viral fitness. However, given that oseltamivir was rarely prescribed before May 2009 in Taiwan, in the absence of selection pressure from the antiviral, the replacement of clade 2C-2 with clade 2B-2 could not be well explained by the combination effect of H275Y mutation with permissive mutation. In our study, had there been no competition cost (refer to the Strain Y 1.00 in Fig. 4), clade 2C-2 would have generated an earlier epidemic and, under the assumption of full cross-immunity, depleted compartment S to a level at which the real epidemic of clade 2B-2 would not have occurred. In terms of transmission dynamics, our model successfully adopted an ever innovative approach dissecting into strain interactions of influenza epidemics through a quantitative way. The embedded significance of the competition cost unveiled through our study probably plays an important role in deciphering the evolution trajectory of the influenza virus and deserves further research to discover its nature at the molecular level.

Thirdly, with higher-resolution information, integrating laboratory surveillance with modelling techniques could provide a good opportunity to simulate seasonal influenza epidemics specifically to the levels of type, clade, and even strain. Coherent with the findings of Shaman et al. [24, 25], absolute humidity could explain the theme of influenza seasonality at the population level in our study. However, to explore the competition of influenza strains, clade-based epidemiological data as well as weekly A/H1N1 IVIR were further exploited to excavate the underlying processes that might not be visible without skillful integration. The fact that our model could report for the 2008–2009 influenza season the best-fit simulations with correlation coefficient of up to 96% signifies the invaluable collaboration of, among others, laboratory, epidemiology, and mathematical modelling.

Through the transmission model we constructed, in terms of the R e’s estimated for clade 2B-2 and clade 2C-2, relative fitness (RF Z ) and competition cost (CC Y ) along with the basic reproductive number range (R 0max and R 0min) could establish an epidemic 96.83% correlated to the real epidemic of 2008–2009 influenza season. Mathematical modelling has the potential to probe the seemingly intractable complexity of infectious disease dynamics and offers valuable tools to dissect epidemiological surveillance. As illustrated in the “Results” section (see Figs. 3, 4), the vicissitudes of two co-circulating clades could be elucidated and, in a sense, be predicted if we could further study the nature of relative fitness and competition cost. Besides, in our study, witnessing the emergence of the succeeding clade 2B-2 in September 2008 with the recession of the preceding clade 2C-2 on day 167 and even the displacement of the latter by the former on day 210 (see Table 3), our model would argue that it could provide a promising opportunity predicting the epidemic trend for an influenza season and assist in decision making for public health preparedness and response.

Facing the progress in the emergence of drug-resistant A/H1N1pdm09 variants, Oh et al. [39] proposed an experimental ferret model to assess the viral fitness of influenza viruses in terms of within-host replication and between-host transmission. Petrie et al. [40] further found that the within-host replication fitness of the oseltamivir-resistant A/H1N1pdm09 variants was not compromised relative to that of related oseltamivir-sensitive strains. Hurt et al. [15] analyzed a widespread community cluster of oseltamivir-resistant A/H1N1pdm09 influenza in Australia and warned that widespread emergence of oseltamivir-resistant A/H1N1pdm09 strains might be more likely. At the human population level, we were therefore concerned with the between-host transmission fitness of the drug-resistant A/H1N1pdm09 variants for better preparedness strategy and response capacity.

Being a study of the aggregate of individual persons, our study reflected the summation of the transmission dynamics of co-circulating influenza strains in the human population. Thanks to the contribution of virological surveillance, epidemiological documentation, vaccination program recording, demographic registry, and meteorological observations, the competition dynamics of co-circulating influenza strains could be explored in a quantitative way to simulate the real epidemic. Worthy of note, because our model took the advantage of the antiviral-free background before May 2009 in Taiwan, the results of our study were implicitly assumed independent of selection pressure exerted by antivirals and de novo antiviral resistance generated during prophylaxis or treatment. However, as sagacious antiviral usage becomes indispensable part of the portfolio for influenza control, factors concerning selection pressure, de novo resistance, and cross-resistance of antivirals shall be further incorporated into the transmission model in the case that the corresponding data could be available in the real word to facilitate the communication among basic, clinical, and population sciences.

Although clade-based epidemiological data were adopted, the scope of our study was actually at the human population level. Our model made an implied assumption that the interaction of influenza strains, through transmission dynamics, would manifest itself in the trend of seasonal epidemics. Because of the limited resolution of the laboratory surveillance in our study, we did not intend to investigate viral fitness at the molecular level, which is though fundamentally important. Nevertheless, we identified the nature and the magnitude of the interaction, and found that A/H1N1 clade 2B-2 outcompeted clade 2C-2 as the latter encountered the former and suffered from 71% reduction in the transmission fitness. Furthermore, we found that clade 2B-2 would be 24.06% higher in transmission fitness than clade 2B-1. These findings are not only in a sense consistent with the results from the animal model of Bouvier et al. [36] but also provide quantitative assessment of the enhanced human transmissibility of oseltamivir-resistant strains during the 2007–2009 influenza seasons.

Conclusions

To our knowledge, this is the first time that the competition dynamics of co-circulating influenza strains was studied in a quantitative way at the human population level. Our study provides quantitative descriptions about the transmission efficiency of each circulating clade per se and brings to light the transmission dynamics of the A/H1N1 H275Y strains during the 2007–2009 influenza seasons worldwide. An interdisciplinary quantitative approach through mathematical modelling could probably facilitate the resolution process of the transmissibility of continually emerging drug-resistant A/H1N1pdm09 variants.

Abbreviations

A/H1N1pdm09:

pandemic (A/H1N1) 2009

CC:

competition cost

IVIR:

influenza virus isolation rate

R e :

effective reproductive number

RF:

relative fitness

References

  1. Nelson MI, Holmes EC. The evolution of epidemic influenza. Nat Rev Genet. 2007;8(3):196–205.

    Article  CAS  PubMed  Google Scholar 

  2. European Centre for Disease Prevention and Control. Pandemics of the 20th–21st centuries. http://ecdc.europa.eu/en/healthtopics/pandemic_preparedness/basic_facts/Pages/historical_pandemics.aspx. Accessed Nov 2013.

  3. Longini IM Jr, Halloran ME, Nizam A, Yang Y. Containing pandemic influenza with antiviral agents. Am J Epidemiol. 2004;159:623–33.

    Article  PubMed  Google Scholar 

  4. Patel R, Longini IM Jr, Halloran ME. Finding optimal vaccination strategies for pandemic influenza using genetic algorithms. J Theor Biol. 2005;234(2):201–12.

    Article  PubMed  Google Scholar 

  5. Ferguson NM, Cummings DA, Cauchemez S, Fraser C, Riley S, Meeyai A, Iamsirithaworn S, Burke DS. Strategies for containing an emerging influenza pandemic in Southeast Asia. Nature. 2005;437(7056):209–14.

    Article  CAS  PubMed  Google Scholar 

  6. Germann TC, Kadau K, Longini IM Jr, Macken CA. Mitigation strategies for pandemic influenza in the United States. Proc Natl Acad Sci. 2006;103(15):5935–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Lipsitch M, Cohen T, Murray M, Levin BR. Antiviral resistance and the control of pandemic influenza. PLoS Med. 2007;4(1):e15.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Halloran ME, Ferguson NM, Eubank S, Longini IM Jr, Cummings DA, Lewis B, Xu S, Fraser C, Vullikanti A, Germann TC, et al. Modeling targeted layered containment of an influenza pandemic in the United States. Proc Natl Acad Sci. 2008;105(12):4639–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Samson M, Pizzorno A, Abed Y, Boivin G. Influenza virus resistance to neuraminidase inhibitors. Antiviral Res. 2013;98(2):174–85.

    Article  CAS  PubMed  Google Scholar 

  10. Mai LQ, Wertheim HFL, Duong TN, van Doorn HR, Hien NT, Horby P. A community cluster of oseltamivir-resistant cases of 2009 H1N1 influenza. N Engl J Med. 2010;362:86–7.

    Article  CAS  Google Scholar 

  11. Storms AD, Gubareva LV, Su S, Wheeling JT, Okomo-Adhiambo M, Pan CY, Reisdorf E, St George K, Myers R, Wotton JT, et al. Oseltamivir-resistant pandemic (H1N1) 2009 virus infections, United States, 2010–11. Emerg Infect Dis. 2012;18(2):308–11.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Lackenby A, Moran Gilad J, Pebody R, Miah S, Calatayud L, Bolotin S, Vipond I, Muir P, Guiver M, McMenamin J, et al. Continued emergence and changing epidemiology of oseltamivir-resistant influenza A(H1N1)2009 virus, United Kingdom, winter 2010/11. Euro Surveill. 2011;16:1–6.

    Google Scholar 

  13. Hurt AC, Chotpitayasunondh T, Cox NJ, Daniels R, Fry AM, Gubareva LV, Hayden FG, Hui DS, Hungnes O, Lackenby A, et al. Antiviral resistance during the 2009 influenza A H1N1 pandemic: public health, laboratory, and clinical perspectives. Lancet Infect Dis. 2012;12:240–8.

    Article  PubMed  Google Scholar 

  14. Hurt AC, Hardie K, Wilson NJ, Deng YM, Osbourn M, Gehrig N, Kelso A. Community transmission of oseltamivir-resistant A(H1N1)pdm09 influenza. N Engl J Med. 2011;365:2541–2.

    Article  CAS  PubMed  Google Scholar 

  15. Hurt AC, Hardie K, Wilson NJ, Deng YM, Osbourn M, Leang SK, Lee RT, Iannello P, Gehrig N, Shaw R, et al. Characteristics of a widespread community cluster of H275Y oseltamivir-resistant A(H1N1)pdm09 influenza in Australia. J Infect Dis. 2012;206(2):148–57.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Meijer A, Jonges M, van Beek P, Swaan CM, Osterhaus AD, Daniels RS, Hurt AC, Koopmans MP. Oseltamivir-resistant influenza A(H1N1)pdm09 virus in Dutch travellers returning from Spain, August 2012. Euro Surveill. 2012;17:1–7.

    Google Scholar 

  17. Takashita E, Kiso M, Fujisaki S, Yokoyama M, Nakamura K, Shirakura M, Sato H, Odagiri T, Kawaoka Y, Tashiro M. Characterization of a large cluster of influenza A(H1N1)pdm09 viruses cross-resistant to oseltamivir and peramivir during the 2013–2014 influenza season in Japan. Antimicrob Agents Chemother. 2015;59(5):2607–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Ampofo WK, Al Busaidy S, Cox NJ, Giovanni M, Hay A, Huang S, Inglis S, Katz J, Mokhtari-Azad T, Peiris M, et al. Strengthening the influenza vaccine virus selection and development process: outcome of the 2nd WHO Informal Consultation for Improving Influenza Vaccine Virus Selection held at the Centre International de Conferences (CICG) Geneva, Switzerland, 7 to 9 December 2011. Vaccine. 2013;31(32):3209–21.

    Article  PubMed  Google Scholar 

  19. Yang JR, Lin YC, Huang YP, Su CH, Lo J, Ho YL, Yao CY, Hsu LC, Wu HS, Liu MT. Reassortment and mutations associated with emergence and spread of oseltamivir-resistant seasonal influenza A/H1N1 viruses in 2005–2009. PLoS ONE. 2011;6(3):e18177.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Liu CH, Wang JL, Su CP, Chuang JH, Chang CH, Lai MS. Oseltamivir use and outcomes during the 2009 influenza A H1N1 pandemic in Taiwan. BMC Public Health. 2013;13:646–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Shih SR, Chen GW, Yang CC, Yang WZ, Liu DP, Lin JH, Chiu SC, Chen HY, Tsao KC, Huang CG, et al. Laboratory-based surveillance and molecular epidemiology of influenza virus in Taiwan. J Clin Microbiol. 2005;43(4):1651–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Jian JW, Chen GW, Lai CT, Hsu LC, Chen PJ, Kuo SH, Wu HS, Shih SR. Genetic and epidemiological analysis of influenza virus epidemics in Taiwan during 2003 to 2006. J Clin Microbiol. 2008;46(4):1426–34.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Wallinga J, Teunis P, Kretzschmar M. Using data on social contacts to estimate age-specific transmission parameters for respiratory-spread infectious agents. Am J Epidemiol. 2006;164(10):936–44.

    Article  PubMed  Google Scholar 

  24. Shaman J, Kohn M. Absolute humidity modulates influenza survival, transmission, and seasonality. Proc Natl Acad Sci. 2009;106(9):3243–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Shaman J, Pitzer VE, Viboud C, Grenfell BT, Lipsitch M. Absolute humidity and the seasonal onset of influenza in the continental United States. PLoS Biol. 2010;8(2):e1000316.

    Article  PubMed  PubMed Central  Google Scholar 

  26. World Health Organization. Pandemic influenza risk management. WHO Interim Guidance. Geneva: World Health Organization; 2013.

    Google Scholar 

  27. Anderson RM, May RM. Infectious diseases of humans: dynamics and control. Oxford: Oxford University Press; 1991.

    Google Scholar 

  28. Hannoun C, Megas F, Piercy J. Immunogenicity and protective efficacy of influenza vaccination. Virus Res. 2004;103(1–2):133–8.

    Article  CAS  PubMed  Google Scholar 

  29. Clements ML, Betts RF, Tierney EL, Murphy BR. Serum and nasal wash antibodies associated with resistance to experimental challenge with influenza A wild-type virus. J Clin Microbiol. 1986;24:157–60.

    CAS  PubMed  PubMed Central  Google Scholar 

  30. Levy A, Sullivan SG, Tempone SS, Wong KL, Regan AK, Dowse GK, Effler PV, Smith DW. Influenza vaccine effectiveness estimates for Western Australia during a period of vaccine and virus strain stability, 2010 to 2012. Vaccine. 2014;32(47):6312–8.

    Article  PubMed  Google Scholar 

  31. Moldoveanu Z, Clements ML, Prince SJ, Murphy BR, Mestecky J. Human immune responses to influenza virus vaccines administered by systemic or mucosal routes. Vaccine. 1995;13:1006–12.

    Article  CAS  PubMed  Google Scholar 

  32. Diekmann O, Heesterbeek JAP. Mathematical epidemiology of infectious diseases: model building, analysis and interpretation. Chichester: Wiley; 2000.

    Google Scholar 

  33. Dorigatti I, Cauchemez S, Ferguson NM. Increased transmissibility explains the third wave of infection by the 2009 H1N1 pandemic virus in England. Proc Natl Acad Sci. 2013;110:13422–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Dharan NJ, Gubareva LV, Meyer JJ, Okomo-Adhiambo M, McClinton RC, Marshall SA, St. George K, Epperson S, Brammer L, Klimov AI, et al. Infections with oseltamivir-resistant influenza A(H1N1) virus in the United States. JAMA. 2009;301:1034–41.

    Article  CAS  PubMed  Google Scholar 

  35. Meijer A, Lackenby A, Hungnes O, Lina B, Van Der Werf S, Schweiger B, Opp M, Paget J, Van De Kassteele J, Hay A, et al. Oseltamivir-resistant influenza virus A (H1N1), Europe, 2007–08 season. Emerg Infect Dis. 2009;15(4):552–60.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Bouvier NM, Rahmat S, Pica N. Enhanced mammalian transmissibility of seasonal influenza A/H1N1 viruses encoding an oseltamivir-resistant neuraminidase. J Virol. 2012;86(13):7268–79.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Hurt AC. The epidemiology and spread of drug resistant human influenza viruses. Curr Opin Virol. 2014;8:22–9.

    Article  PubMed  Google Scholar 

  38. Bloom JD, Gong LI, Baltimore D. Permissive secondary mutations enable the evolution of influenza oseltamivir resistance. Science. 2010;328(5983):1272–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Oh DY, Hurt AC. A review of the antiviral susceptibility of human and avian influenza viruses over the last decade. Scientifica. 2014;2014:430629.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Petrie SM, Butler J, Barr IG, McVernon J, Hurt AC, McCaw JM. Quantifying relative within-host replication fitness in influenza virus competition experiments. J Theor Biol. 2015;382:259–71.

    Article  PubMed  Google Scholar 

Download references

Authors’ contributions

KC and BH conceived the study and developed the model. BH implemented the model and compiled the data. KC and BH analyzed and interpreted the data, and contributed to the development of the manuscript. KC was responsible for submission of the publication. Both authors read and approved the final manuscript.

Acknowledgements

We acknowledge the hard work of the national influenza surveillance network coordinated by Taiwan Centers for Disease Control.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

All data used, generated, and analyzed during the current study are included within the article and its Additional files 1, 2, 3, 4, 5, 6 and 7. Computer code is available upon request from the corresponding author on reasonable request.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Funding

KC was supported in part by MOST Grants 104-2221-E-002-046-MY3 and 103-2221-E-002-157-MY3, Taiwan.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information

Authors and Affiliations

  1. Department of Computer Science and Information Engineering, National Taiwan University, Taipei, Taiwan, ROC

    Bin-Shenq Ho & Kun-Mao Chao

  2. Public Health Bureau, Hsinchu, Taiwan, ROC

    Bin-Shenq Ho

  3. Taiwan Centers for Disease Control, Taipei, Taiwan, ROC

    Bin-Shenq Ho

  4. Graduate Institute of Biomedical Electronics and Bioinformatics, National Taiwan University, Taipei, Taiwan, ROC

    Kun-Mao Chao

Authors
  1. Bin-Shenq Ho
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Kun-Mao Chao
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Kun-Mao Chao.

Additional files

12967_2017_1269_MOESM1_ESM.docx

Additional file 1: Table S1. The weekly A/H1N1 influenza virus isolation rate (A/H1N1 IVIR) of the 2007-2008 and 2008-2009 influenza seasons. The surveillance system covered approximately 75% of the 352 basic administrative units throughout the northern, central, southern, and eastern regions of Taiwan. Considering that the original weekly data reflected the data variability of the national influenza surveillance network, we therefore adopted the aggregate weekly data, computed from all positive isolations divided by all collected specimens of the three consecutive weeks.

12967_2017_1269_MOESM2_ESM.docx

Additional file 2: Table S2. The mid-year population structure of the year 2007 and 2008. The mid-year age-structured population sizes were obtained from Department of Household Registration, Ministry of the Interior (http://www.ris.gov.tw/zh_TW/346). The model population was classified into six age groups: 0-5, 6-12, 13-19, 20-39, 40-59, and ≥ 60 years.

Additional file 3: Appendix. Transmission model construction description.

12967_2017_1269_MOESM4_ESM.docx

Additional file 4: Table S3. Meteorological data used to derive the monthly specific humidity during the 2007-2008 and 2008-2009 influenza seasons. Refer the meteorological variables used to derive specific humidity to Taiwan Central Weather Bureau (http://www.cwb.gov.tw/V7/climate/monthlyData/mD.htm).

12967_2017_1269_MOESM5_ESM.docx

Additional file 5: Table S4. Progress of annual influenza vaccination program from October 1 through December 30, year 2007 and 2008 respectively. The year 2007 influenza vaccination program started on October 1 and aimed at those aged 6 months through 2 years, the school children aged 6-7, and the elderly aged 65 years and above. From December 1 on, the program was open to the general population to maximize the vaccine utilization. As to the year 2008 influenza vaccination program, starting on October 1, it aimed at those aged 6 months through 3 years, the school children aged 6-9, and the elderly aged 65 years and above before December 1 and open to the general population thereafter

12967_2017_1269_MOESM6_ESM.docx

Additional file 6: Table S5. Parameter sets with the corresponding estimated monthly effective reproductive numbers and annual effective reproductive number of each simulation for the 2007-2008 influenza season. For a single-clade influenza season, the model found successful simulations and evaluated the effective reproductive numbers. The annual effective reproductive number was deduced from the geometric mean of the estimated monthly effective reproductive numbers. The range of the basic reproductive number was estimated accordingly.

12967_2017_1269_MOESM7_ESM.docx

Additional file 7: Table S6. Parameter sets with the corresponding estimated monthly effective reproductive numbers and annual effective reproductive number of each simulation for the 2008-2009 influenza season. For a two-clade influenza season, the model found successful simulations and evaluated the effective reproductive numbers. The annual effective reproductive number was deduced from the geometric mean of the estimated monthly effective reproductive numbers. Relative fitness of succeeding clade and competition cost of the preceding clade, as well as the range of the basic reproductive number, were estimated accordingly.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ho, BS., Chao, KM. Data-driven interdisciplinary mathematical modelling quantitatively unveils competition dynamics of co-circulating influenza strains. J Transl Med 15, 163 (2017). https://doi.org/10.1186/s12967-017-1269-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12967-017-1269-6

Keywords

Journal of Translational Medicine

ISSN: 1479-5876

Contact us