4Flu - an individual based simulation tool to study the effects of quadrivalent vaccination on seasonal influenza in Germany

Before 1985, circulating Influenza B viruses were closely related to the precursor of the subsequent B/Yamagata (B/Yam) lineage [1]. Since 2001, viruses of the distinct B/Victoria (B/Vic) lineage have been detected in clinical samples from various parts of the world. The two Influenza B lineages are antigenically distinct [1–3]; genome sequence analysis revealed that they substantially differ as to the gene encoding the surface hemagglutinin [1, 2] and that they continue to diverge [4, 5]. Viruses of both B lineages now co-circulate together with the Influenza A viruses A(H3N2) and A(H1N1) [6]. The incidence of Influenza B infection varies substantially from year to year [6]. Surveillance data of the seasons from 2001/02 to 2010/11 show that Influenza B is responsible for 1 to 60% (mean 23%) of influenza cases in Europe and for <1 to 44% (mean 24%) in the United States, respectively [6]. Influenza B infections cause severe disease in all age groups [7] and their clinical symptoms and outcomes are similar to those of Influenza A infections [8].

Influenza vaccine design is intended to anticipate which variants will circulate in the subsequent transmission season. The antigenic composition of the vaccine is revised twice a year by the World Health Organization (WHO) and adjusted to the antigenic characteristics of circulating influenza viruses in the Northern and Southern hemisphere. Current trivalent influenza vaccines (TIV) contain antigens of A(H3N2) and A(H1N1), but of only one Influenza B lineage [9]. The efficacy of these influenza vaccines varies with the age of the vaccinees and with the degree at which the vaccines match the circulating influenza strains. In case of a good match, their efficacy can range from 70% to 90% in young adults, but it is lower in very young children and in elderly [9]. A recent meta-analysis of randomized controlled vaccine efficacy trials showed that TIV confers a 52% protection (95% confidence interval (CI): 19 to 72%) against the B-lineage not contained in the vaccine, compared to 77% (95% CI: 18 to 94%) for the lineage which matches the vaccine [10]. Co-circulation, fluctuating prevalence, and differing geographical distributions of the two B lineages have challenged global recommendations for vaccine composition [3]. In the US, the Influenza B lineage in TIV failed to match the predominant circulating B lineage in five out of ten seasons between 2001 and 2011 [11]. In 2012, WHO recommended for the first time vaccine strains from both B lineages, paving the way for the development of quadrivalent influenza vaccines (QIV) which contain strains of two B lineages and two A subtypes [12].

A variety of deterministic and stochastic models have been proposed to address various aspects of seasonal [13–16] or pandemic influenza [17–21] or to predict the size and timing of seasonal waves [22–24]. We have decided to use an individual-based stochastic model which allows for maximum flexibility in modeling immunity and vaccination strategies [25] and, thus, circumvents common oversimplifications of deterministic models [26]. We have developed the simulation tool 4Flu to compare the effect of trivalent and quadrivalent vaccination on the incidence of influenza infection in Germany. 4Flu is based on an individual-based stochastic model, which allows for the independent and simultaneous transmission of the four influenza variants A(H1N1), A(H3N2), B/Vic and B/Yam. We use the official German age-distribution for simulation years which lie in the past, and official demographic predictions for the future [27]. Individuals are connected by a dynamically changing contact network which is synchronized with the German POLYMOD contact matrix [28]. The simulation model further allows for the introduction of new drift variants and for mismatched vaccine design against new drift variants. Special care has been taken to provide a realistic immunity model which allows for maternal protection of newborns, for booster effects and waning of existing immunity and for cross-immunizing effects of infections and vaccinations. Vaccination either uses TIV or QIV. Each simulation runs for several decades and reports the annual incidence of influenza infections. Simulations are typically repeated one thousand times to average out random effects.

We present an individual-based stochastic simulation model for the simultaneous and independent transmission of four influenza strains (A(H1N1), A(H3N2), B/Vic, B/Yam). Each simulation runs for 50 years, with each simulation year starting on July 1st: the first 10 years are used to initialize the contact network, the next 20 years for initializing age-dependent infection and immunity patterns in the population, and the final 20 years for comparing vaccines.

Naturally acquired immunity

After infection, individuals acquire a temporary immunity which prevents further infections; this immunity is lost continuously at a constant rate. After losing immunity, individuals become fully susceptible again. If an immune individual gets in contact with an infectious one, or if an immune individual is vaccinated, the duration of his or her immunity can be extended by a booster vaccination or infection (see details in Appendix 4) without rendering the individual infectious. As simulation models have shown that cross-immunizing events caused by co-circulating B lineages can play an important role in the transmission dynamics of seasonal influenza [14, 26], we have also incorporated this feature in 4Flu: infection with any one of the two B lineages has a 60% probability to lead to a cross-immunizing event against the other B lineage [10]. We do not consider specific cross-protection between the two Influenza A strains or among Influenza A and B, but examine the possible effect of a short-term unspecific protection in a sensitivity analysis. Newborn individuals are protected by maternal antibodies against the influenza strains against which their mother is immune (irrespective of whether her immunity was acquired by infection or vaccination). This maternal protection lasts for 2 to 4 months [33, 34] and is assumed to prevent infection. The dynamics of acquiring and losing immunity is schematically depicted in Figure 1.

Vaccination

Vaccinations are performed in October and November of each simulation year [37], whereby a percentage of individuals is vaccinated, using either TIV or QIV. As shown in Table 1, the vaccinated percentage depends on the age of the individuals, on their risk status and on their vaccination status [38]: individuals who were vaccinated in the previous simulation year have twice the chance of being vaccinated again, compared to previously unvaccinated individuals (see Appendix 5 for details). An age-dependent percentage of previously susceptible vaccinees become immune after vaccination (vaccine efficacy) [39–41]. In case of TIV vaccination, we also consider the possibility of cross-immunizing effects against the B lineage which is not included in the vaccine: the age-dependent vaccine efficacy is multiplied with 0.6 to determine what percentage of vaccinees is protected against the B lineage which is missing in the vaccine. Successful vaccination can also extend the duration of pre-existing immunity (“booster-vaccination”; see Appendix 5 for details). When a new drift variant has just been introduced, the vaccine may not match this variant (vaccine mismatch). In this case, the age-dependent vaccine efficacy is multiplied by 0.6 (calculated from [39]). The decision of whether the vaccine is matched or mismatched to a new drift variant [42] is made identically for the TIV and the QIV scenario. Vaccine efficacy depends on the age of the vaccinee, but is assumed to be independent on the virus subtype or lineage. We use the same vaccine efficacy for TIV and for QIV (Table 1) [39–41], but distinguish between age classes and individual risk status [43]. For convenience, we use the same age-dependent values of the vaccination coverage in every simulation year. In the first simulation year, a total of 26.8% of the population is vaccinated. Due to demographic changes (Figure 2) and due to a high vaccination coverage in the elderly (Table 1), the overall vaccination coverage reaches 29.4% at the beginning of the evaluation period and continues increasing until the end of the simulation when it reaches 33.4%. Vaccination-derived immunity is lost much more quickly than naturally acquired immunity [44].

Parameter	Baseline	Reference
Basic reproduction number R 0	1.575	Calibration
Maximum seasonal transmission factor	1.43	[31]
Day of maximum seasonal transmission	Dec. 21st	[31]
Duration of the latent period	2 days	[32]
External infection probability	0.0003/year	Assumed
Duration of the infectious period		[32]
- children (age 0–17 years)	4 days
- adults (age 18 years and above)	2 days
Duration of maternal protection	2 - 4 months	[33, 34]
Immunity loss rate after infection (irrespective of influenza subtype or lineage)	1/(9.13 years)	Calculated from [31]
Average circulation time per drift variant		Calculated from [35] and [3]
- A(H1N1), B/Vic, B/Yam	7.0 years
- A(H3N2)	3.5 years
Vaccination coverage		[38]
- no risk group, 0.5-2 years of age	19.2%
- no risk group, 3–6 years of age	22.4%
- no risk group, 7–10 years of age	23.6%
- no risk group, 11–15 years of age	11.0%
- no risk group, 16–59 years of age	16.9%
- no risk group, 60 years of age or older	48.8%
- risk group, 0.5-59 years of age	33.0%
- risk group, 60 years of age or older	64.9%
Revaccination preference factor	2	Assumed
Probability of mismatched vaccine design when a new drift variant occurs	40%	Calculated from [3, 35, 42]
Vaccine efficacy (well-matched vaccine; irrespective of influenza subtype or lineage)
- 0–1 year of age	45%	[40]
- 2–5 years of age	39%	[40]
- 6–15 years of age	69%	[40]
- 16–64 years of age	73%	[39]
- 65 years of age or older	58%	[41]
Cross protection after infection		Calculated from [39]
- percentage of individuals who are immunized against a B lineage when they are infected (or booster-infected) with the other B lineage (lineage cross protection)	60%
- percentage of individuals who were immune against the previous drift variant, who are still protected against the new one (drift cross protection)	60%
Cross protection after vaccination		Calculated from [39]
- vaccine efficacy multiplication factor for B lineage not contained in TIV (lineage cross protection)	0.6
- vaccine efficacy multiplication factor for vaccinations with drift mismatch (drift cross protection)	0.6
Immunity loss rate after vaccination	1/(1.8 years)	Calculated from [44]
Percentage of the population with elevated risk		[43]
- newborn individuals	3.0%
- age 0 to 15 years	6.0%
- age 16 to 59 years	14.2%
- age 60 years and above	47.1%

Figure 2 shows how the contact matrix changes over time: Figure 2a shows the state of the contact network before initialization: the 5-year age classes of the German POLYMOD matrix are slightly disturbed by the German demography, but are still clearly visible. After about five years of initialization, the contact structure has changed to mostly allocate contacts to individuals of the same age, although (by construction of the algorithm) the summary statistics of the number of “outgoing” contacts per 5-year age group remain identical to the POLYMOD matrix. After ten years of network initialization, the contact network has converged to a steady state (Figure 2b). After that time, the immunologic and epidemiologic initializations are started and the demography of the simulated population is allowed to change according to the official German demography. The resulting age distributions of the German population are shown on the right hand sides of the contact matrix visualizations in Figure 2b-f.

As the number of children keeps declining and the number of elderly grows in the German population, the frequency of contacts becomes less focused on children and young adults and tends to be more spread out over the ages. In 2012, contacts are most frequent in juveniles and then decline over age (Figure 3b). Taking into consideration the longer duration of infectiousness of children and juveniles and considering the susceptibility status of their contacts, we can calculate how many individuals are infected on average by a single infectious individual of a given age as shown in Figure 3c. The median time until an infected individual passes on the infection (generation period) is 3.3 days in our simulations.

Our model is calibrated to observed data from Germany. As our simulation results vary strongly (coefficient of variation CV = 91.5% for R ₀  = 1.575), we use the median of the simulated incidence when evaluating simulation results: for each value of the basic reproduction number R ₀ between 1.5 and 1.6 (step size 0.025), we ran 1,000 to 2,000 simulations and calculated the median of the infection incidence in the 2006/07 transmission season for young adults. A basic reproduction number of R ₀  = 1.575 yields a median infection incidence of nearly exactly 10.6% as was observed among German health care workers and controls in that season [45]. Figure 4a shows the large fluctuation of annual influenza waves in a single simulation in which TIV is used throughout, whereby all infections with the two A strains and the two B lineages are added up. Figure 4b selectively shows the incidence with B/Yam for the same simulation. Years in which B/Yam is not contained in TIV are shown as dark grey bars. New B/Yam drift variants are introduced in 2015, 2022, 2028 and 2031, shown as vertical dashed lines. In 2028, the B/Yam lineage in the vaccine does not match the newly introduced drift variant, indicated by a light grey background. During two seasons, large B/Yam outbreaks occur, dwarfing the smaller waves in the other seasons. Starting with the year 2013, simulations branch into a TIV scenario and a QIV scenario (Figure 4c). The cumulative numbers of infections are calculated in each of these scenarios and are used to calculate the average annual number of infections prevented by QIV as compared to TIV.

Averaged annual simulation results for the baseline set of parameter values are given in Table 2. As the simulated population is exactly one thousandth of the German population, these results have to be multiplied with 1,000 to obtain estimates for the whole of Germany. The annual number of influenza infections in Germany can be reduced from 8,943,000 to 8,548,000, i.e. by 395,000 infections per year without increasing the vaccination coverage if QIV is used instead of TIV. The prevented infections constitute 4.3% (95% confidence interval (CI): 4.1-4.5%) of all influenza infections which still occur under TIV vaccination. The simulated fraction of Influenza B infections among all infections varies widely from one season to the next; even if averaged over the 20 years of each simulation, the 95% reference interval for the percentage of Influenza B infections among all influenza infections ranges from 24.6 to 51.3% (median 37.3%) in the TIV scenario, and from 21.3 to 49.1% (median 34.3%) in the QIV scenario. Switching from TIV to QIV prevents on average 11.2% (CI: 10.7-11.8%) of all Influenza B infections which still occur under TIV vaccination. In 53% of 10,000 simulation years, the B lineage not contained in TIV caused more infections than the B lineage which was contained in TIV.Additionally to lowering the expected number of infections, QIV furthermore reduces the number of seasons with rather extreme Influenza B outbreaks (Figure 5).

Age group	Annual mean number of infections using TIV	Annual mean number of infections using QIV	Annual mean number of prevented infections	Percentage prevented
0-15	2152	2073	79	3.6
	(2140–2164)	(2061–2085)	(75–83)	(3.4-3.8)
16-60	5470	5247	223	4.0
	(5433–5506)	(5211–5282)	(212–235)	(3.8-4.2)
61+	1321	1228	93	6.9
	(1311–1331)	(1219–1237)	(90–97)	(6.7-7.2)
all	8943	8548	395	4.3
	(8885–9001)	(8491–8604)	(376–414)	(4.1-4.5)

The risk that in Germany more than 10 million Influenza B infections occur in a single transmission season is about 13.4% for the TIV scenario and 11.2% for the QIV scenario (Figure 5a), meaning that the expected duration between such extreme events can be increased from 7.5 years (TIV) to 8.9 years (QIV). For 15 million infections per year, the difference is even more pronounced (21.3 years for TIV, 29.2 years for QIV). The risk of strong Influenza B seasons can be reduced by up to 27% if QIV is used (Figure 5b).

We address parameter uncertainty in various sensitivity analyses: the basic reproduction number, R ₀ , is varied from 1.2 to 2.0 (baseline value: 1.575) while the values of the other parameters are kept at their baseline values. As R ₀ summarizes the contagiousness of the disease, the average annual number of infections in Germany grows with increasing R ₀ from 3,354,000 to 12,973,000 per year in the TIV scenario, and from 3,029,000 to 12,587,000 in the QIV scenario (Figure 6a-c and Appendix 8); the average number of infections annually prevented by QIV rises from 325,000 infections for R ₀  = 1.2, to about 386,000 infections for R ₀  = 1.5, and, thereafter, remains about constant on a level of 395,000 infections. Due to the rising denominator, the percentage of prevented infections falls with increasing R ₀ from 9.5 to 3.0% (all infections), or from 21.6 to 7.8% (Influenza B infections).

In a further sensitivity analysis, the rate at which immunity is lost is varied such that the duration of naturally acquired immunity lasts on average 4 to 12 years against A(H1N1) and the B lineages, and for 3.2 to 7.1 years against A(H3N2) (see Appendix 4 for details). The average annual number of infections drops with increasing immunity duration from 11,267,000 to 6,735,000 infections per year in the TIV scenario, and from 10,751,000 to 6,453,000 in the QIV scenario (Figure 7a-c and Appendix 8). The average number of infections annually prevented by QIV correspondingly drops from 517,000 (4.6% of all infections, or 11.7% of influenza B infections) to 282,000 infections (4.1% of all infections, or 10.7% of influenza B infections).

Following the parameter choice of Vynnycky et al.[31], we perform a sensitivity analysis in which we only increase the duration of naturally acquired Influenza B immunity from the baseline value of 6 years to 12 years while keeping the immunity duration of Influenza A at the baseline value. This reduces the annual number of infections from baseline 8,943,000 to 8,098,000 in the TIV scenario, and from baseline 8,548,000 to 7,814,000 in the QIV scenario. The average number of infections annually prevented by QIV correspondingly drops from baseline 395,000 (4.3% of all infections, or 11.2% of all Influenza B infections) to 284,000 infections (3.5% of all infections, or 10.8% of all Influenza B infections). Increasing the duration of Influenza B immunity shifts the 95% reference intervals of the percentage of Influenza B infections among all influenza infections towards lower values, reaching 16.3-46.6% (median 30.5%) for the TIV scenario, and 13.7-43.3% (median 27.5%) for the QIV scenario.

Another highly influential parameter is the degree of B lineage cross protection: in a sensitivity analysis, we increase the percentage of individuals who are immunized against an Influenza B lineage when they are successfully vaccinated or infected with the other B lineage, using values from 0% up to the baseline value of 60% (Appendix 8). The higher the probability of B lineage cross protection is, the fewer Influenza B infections occur (Figure 8a-c). Whereas for the baseline value, the median percentage of Influenza B infections among all infections is 37.3% (TIV scenario), it increases to 54.1% if no B lineage cross protection is considered. Accordingly, the number of infections which can annually be prevented by QIV increases from 395,000 (baseline value) up to 2,180,000 infections if lower values are considered. If the degree of B lineage cross protection decreases over time from initially 60% to a lower value which is reached at the end of the simulation, the additional benefit of using QIV increases from the baseline result of 395,000 infections per year (obtained for a constant level of 60% B lineage cross protection) to 559,000 (decreasing the cross protection from 60% to 50%), to 706,000 (40%) and to 853,000 (30%), respectively (average values of 1,000 simulations for each parameter set). Our model does not explicitly consider cross immunity between the A strains or among Influenza A and B, but we have run a sensitivity analysis where we consider a short-term unspecific immunity which prevents infection with any influenza virus within the first weeks following an infection. As shown in Appendix 8, such an indirect protection, which may last for one to four weeks, has only a marginal effect on the results. If anything, it may slightly increase the difference between the QIV and the TIV scenario.

In order to keep the age distribution of “outgoing contacts” synchronized to the POLYMOD contact matrix despite demographic changes, contacts have to be removed or added daily. Although we keep the number of changes to a minimum, only 0.36% of the childhood contacts of 10 year old children at the beginning of the simulation are still among their contacts at the end of the simulation when they are 50 years old. We have examined the influence of contact turnover on our results in a sensitivity analysis in which we have either increased the contact turnover by replacing at random 10% of all contacts each year, or decreased contact turnover by updating contacts only once a year instead of daily (resulting in larger and longer lasting deviations of the simulated age-dependent contact distribution from the POLYMOD matrix). The effect of contact turnover on the number of infections additionally prevented by QIV is marginal (Appendix 8).

Our results indicate that annually an average number of 395,000 influenza infections (representing 11.2% of all Influenza B infections) can additionally be prevented without increasing vaccination coverage by switching from TIV to QIV. Although the composition of TIV during the evaluation period of our simulations is determined at random, our finding that in 53% of seasons, the B lineage not contained in TIV causes more infections than the B lineage which actually is contained in the vaccine is practically identical to findings from the US (the Influenza B lineage in TIV failed to match the predominant circulating B lineage in five out of ten seasons between 2001 and 2011) [11]. It had to be expected that this percentage should even be higher if the vaccination coverage was higher, yet this finding frequently is regarded as an indicator of having chosen the wrong B lineage in TIV.

As QIV is designed to additionally prevent Influenza B infections, our results strongly depend on the assumptions made on the transmissibility and on the duration of acquired immunity against Influenza B. Cross immunization after infection with an Influenza B lineage or after successful vaccination with TIV increases Influenza B immunity in the population. As we show in sensitivity analyses (Figure 8a-c and Appendix 8), our baseline value of 60% lineage cross protection leads to a highly conservative estimate of the effects of QIV. Using lower cross protection values, the additional benefit of using QIV instead of TIV can be up to 5.5 times as high as the baseline result of preventing 395,000 infections per year (Figure 8a-c). As the antigenic overlap of the two B lineages is clearly disadvantageous to their transmission, mutations which increase the genetic difference between the two lineages and, thus, reduce the degree of cross protection should be selected for. To check the effect of this hypothetical antigenic drift, we have continually lowered the degree of B lineage cross protection over time in a sensitivity analysis. If the degree of B lineage cross protection drops from initially 60% to finally 30%, the additional benefit of using QIV increases up to an average of 853,000 infections per year. It must be noted that we restrict cross-immunizing vaccine effects to TIV and omit them when individuals receive QIV, because it has been shown that increasing antigen amount does not enhance the efficacy of influenza vaccines [57, 58] and because there is no convincing evidence that increasing the amount of antigen enhances the efficacy of influenza vaccine [59].

The simulated absolute and relative numbers of Influenza B infections fluctuate strongly from year to year. Even if they are averaged over 20 simulation years, the percentage of Influenza B infections varies from 24.6 to 51.3% (95% reference interval in the TIV scenario). This range is contained in the range of reported values on clinical Influenza B cases, yet the median percentage of Influenza B infections (37.3%) is higher than the reported average of 23% of clinical cases who are caused by Influenza B [6]. This difference between the simulated percentage of Influenza B infections and the reported percentage of Influenza B cases may partly be due to a different fraction of Influenza A and B cases being reported: it has been shown that Influenza A infections more frequently lead to high fever which may lead to more doctoral visits and finally to more case reports [60]. During a mixed A(H3N2) and B season in the Netherlands, the proportion of laboratory confirmed Influenza A was higher in hospital samples than in GP sentinels [61]. This indicates an underestimation of the incidence with types which tend to cause mild illness when laboratory data include hospital samples. With increasing severity, the percentage of strains which cause milder infections should increasingly be underestimated (e.g. from diagnosis of acute respiratory infections, ARI, to diagnosis of influenza like illness, ILI, or hospitalization). Influenza infections are diagnosed by using Flu-rapid tests. Most of these tests can identify Influenza A and B, but the proportion of selective tests which only test for Influenza A is unknown. This might cause a slight underestimation of the true Influenza B incidence. Vynnycky et al.[31] use in their simulation study different immunity durations for Influenza A (6 years) and Influenza B (12 years). We have addressed this in a sensitivity analysis in which we double the duration of Influenza B immunity (unlike Vynnycky et al. we use the resulting duration of 12 years only for naturally acquired immunity and not for immunity conferred by vaccination). As a consequence, the median percentage of Influenza B infections among all influenza infections drops to 30.5% in the TIV scenario, and the average number of infections prevented by QIV drops to 284,000 per year (preventing 3.5% of all influenza infections, or 10.8% of all Influenza B infections, respectively). If, on the other hand, we assume that the B lineage cross protection is less than 60%, the number of Influenza B infections increases, and, consequently, QIV can prevent more infections (Figure 8a-c). Without any B lineage cross protection, we expect to see more Influenza B than A infections (54.1 vs. 45.9%) in the TIV scenario, because TIV only immunizes against one B lineage, but against both A subtypes. In that extreme case, QIV would additionally prevent over two million infections per year. Although the annual incidence of infections strongly depends on R ₀ , the number of prevented infections is nearly constant for R ₀  > 1.5.

The seasonal influenza waves depicted in Figure 4a show a much higher variability than in the country-wide statistics published by the Robert Koch Institute [42]: the coefficient of variation of the simulated annual number of influenza infections is in the order of 100% (TIV scenario). This high variability mainly results from using a population size of only one thousandth of the German population: geographic differences in immunity patterns and transmission reduce the variability of the reported influenza incidence. When adding up the annual incidences of ten randomly picked simulations (to mimic the combined result of ten locations with independent transmission), the coefficient of variation reduces to about 30%. In reality, much more locations ought to be combined, but their epidemiology must not be assumed to be independent.

Various contact networks have been proposed to study the spread of influenza [62–65]. In our simulations, we use a novel contact network which is based on contact information gathered for the German population [28], yet which allows for dynamic changes during the ageing process of the population. As was shown in a sensitivity analysis, the network turnover had only negligible effect on the difference between the TIV and QIV scenario. It has been suggested to use bipartite networks which explicitly consider the place of contact (school, work place, household) [66]. Although the POLYMOD data include some information on the place of the contact [28], many additional assumptions would be required to explicitly consider contact place in our simulations. Some of the contacts which regularly occur at school or work place [67] are covered by the age-dependency of contacts in the underlying POLYMOD matrix [28], but our contact network may be too rigid to allow for random encounters (e.g. during public transportation).

In our simulation studies, we make several simplifying assumptions: (1) the basic reproduction number is assumed to be equally high for all four influenza strains; differences between Influenza A and B are only due to the higher drift frequency of A(H3N2), to the missing B lineage in TIV, and to B lineage cross protection. (2) Infected children are as contagious as adults, yet they have a longer infectious period (cf. Table 1) and have more contacts than adults - specifically with other children, resulting in a larger number of secondary infections (Figure 3c). (3) We apply the same reduction factor 0.6 to the vaccine efficacy which cross-protects against the B lineage which is not contained in TIV, and to a vaccine design mismatch; likewise, we assume that only 60% of previously immune individuals are immune against a newly introduced drift variant. Tricco et. al.[10] report for vaccinees from 18 to 49 years of age that the vaccine efficacy against Influenza B is 77% if the lineage is contained in TIV and 52% if it is not contained in TIV, respectively. The ratio 52/77 = 0.675 is similar to our ratio of 0.6, calculated from a Cochrane review for any type of mismatch. Langley et al. show in a vaccination study with children that the sero-conversion rate with respect to the B lineage which is not contained in TIV is half as high as the sero-conversion rate with respect to the lineage which is contained in the vaccine [59] (cf. Appendix 6). Using a B lineage cross protection of only 50%, leads to a larger than baseline average annual number of infections prevented by QIV (Figure 8c). (4) School holidays have been shown to have a strong impact on the transmission of influenza [62], but we did not want to further complicate our model by changing contact patterns during such time periods, as different school holiday periods are used in different German states and as they change from year to year. We also do not distinguish between week days and weekends in our simulation. Using averaged social contacts throughout the week has been shown not to affect the major outcome of influenza simulation studies [63]. (5) Although vaccination coverage in Germany varies strongly between different locations and may also slightly vary over time [68], we use a constant age-dependent vaccination coverage throughout our simulations [38] which is assumed to represent a country-wide average. (6) We do not explicitly consider human travel in our simulations [20, 21, 47], but we apply a very small “external” infection rate which is present throughout the year. (7) Individuals who were vaccinated against A(H1N1) pdm09 in 2009 were more likely to be revaccinated in 2010 with an odds ratio (OR) of 6.02 in The Netherlands [69], yet in our simulations, the OR of revaccination ranges from 2.2 to 6.5 (for details, see Appendix 5). If the reported OR is also representative for non-pandemic years in Germany, too few individuals are revaccinated in our simulations and, thus, the effects of vaccinations may be overestimated. (8) We use the same average immunity duration after booster infections as after infections of previously susceptible individuals as it is not known whether this immunity duration should be larger or smaller than after infection. (10) Although great care has been given to use realistic assumptions on the immunity acquired by infection, the parameters associated with immunity are quite uncertain. We have addressed this uncertainty in various sensitivity analyses, which explore the effect of the duration of immunity and the degree of heterologous immunity between the B lineages. It has been difficult to obtain quantifiable evidence for cross-immunizing effects between the A strains or between Influenza A and B, yet an indirect protection may be caused by an increased production of interferons, tumor necrosis factor alpha and beta and other cytokines as has been reported after infection with influenza and various other viral infections [70–73]. As we did not want to incorporate non-influenza infections, we have restricted the study on the effects of unspecific protection to a sensitivity analysis which showed only a minor effect of indirect protection (Appendix 8).

4Flu is a novel simulation tool which describes the spread of influenza viruses using a dynamically evolving contact network, reflecting German demography, with different vaccinations strategies. The results indicate that the current vaccination with TIV in Germany, although covering less than 30% of the total population [38], substantially reduces the number of influenza infections compared to no vaccination. Replacing TIV with QIV would further reduce Influenza B incidence by 395,000 to 2,180,000 influenza infections per year, depending on the degree of B lineage cross protection.

The observed and predicted demographic distributions of Germany are given in 1-year cohorts for every year until 2060, together with annual birth rates [27]. Our simulations start with the demographic distribution of 1993. Each individual is given a random birth day; in each simulation year, individuals are added according to the predicted number of births. If the number of individuals of a cohort declines from one year to the next, individuals are chosen to die during the simulation year; if it increases, individuals are introduced as “immigrants”. Birth, death and immigration take place at random times during the year. On each birthday, an individual’s age grows by one year. Individuals are also given a risk status which influences their probability of being vaccinated and which may influence the clinical outcome of an infection. Some individuals have an increased risk when they are born, others acquire a risk status over time. In order to avoid sudden transitions from one age group to the next one, we assume that the percentage with risk changes linearly within each age group, whereby we keep the total percentage in the age group identical to the literature value. Due to births, ageing and death, the risk fractions in the simulations keep drifting away from the specified values given in Table 1. At each birthday, it is checked whether there are currently too many or too few individuals with risk status in the cohort of the individual who is now one year older; if needed, the individual’s risk status is changed in order to reduce this deviation.

In order to transmit the infection from one individual to another one, individuals are connected in a dynamic contact network. The average number of contacts depends on the age of the individual and on the ages of the contact persons, and is given by the POLYMOD matrix [28]. Before initializing the contact network, each individual “initiates” on average 5 contacts with others, whereby the age distribution of these contacts is proportional to the numbers given by the POLYMOD matrix. These contacts are classified as “outgoing contacts” for the source individual and as “incoming contacts” for the contacted individual (i.e. on average, every individual has 10 contacts). For each newborn individual, a random person of the age-group from 20 to 40 years is chosen (denoted as the “mother”) who initiates a contact to the newborn individual. If the mother is connected to other individuals (which usually is the case), the newborn initiates a contact to the individual whose age is closest to that of the mother; thus newborns are initially connected in triangles. Immigrants are incorporated into the existing contact network in the same way as newborns. Due to births, immigrations, ageing and deaths, the contact distribution of “outgoing” contacts keeps drifting away from the values specified by the POLYMOD matrix. In order to correct this, the distribution of “outgoing” contacts is compared at the end of each simulation day with the POLYMOD matrix assuming a population average of 5 “outgoing” contacts per individual. If too many contacts between two age groups exist, superfluous contacts are chosen at random and removed; if too few contacts exist, random individuals are chosen from the required age groups and connected. In a sensitivity analysis this “repair mechanism” is only performed annually instead of daily (resulting in lower network turnover); or an additional random percentage of contacts is removed and replaced by new random contacts (resulting in higher network turnover). In order to let the contact network adjust to the demographic structure of Germany, we run the simulation for 10 years (“contact network initialization”) without changing the demographic structure of Germany before introducing immunity and infections in the simulations.

An infected individual does not necessarily infect all of his or her contacts; the expected number of secondary infections in a completely susceptible population is given by the basic reproduction number R ₀ . In order to calculate the infection probability per contact from the basic reproduction number, the “next generation matrix” is calculated which combines the contact information in the population and the duration of contagiousness [32]. As we start the simulations with 80,976 individuals, the contact network can be described by an 80,976 × 80,976 matrix which is mostly filled with zeros. The next generation matrix is defined by the product of the network matrix and the average duration D _a of contagiousness (which depends on the age a of the infectious individual). The largest eigenvalue E of this matrix is calculated by power iteration (Von Mises iteration). The probability that a given contact of an infectious individual of age a is infected by that individual is then given by p _a  = R ₀ D _a /E[30]. During the first 5 years of the initialization of the contact structure, the largest eigenvalue of the next generation changes markedly, but after that it remains stable, indicating that the contact network has converged to a steady-state. As influenza transmission is assumed to vary seasonally, we multiply the transmission probability with a seasonal factor: p _a (1 + 0.43cos(2π(t-173)/365)). The seasonality is chosen such that the transmission reaches its maximum around Christmas where it is 43% higher than the average value [31]. When an individual becomes contagious, a random decision is made for each contact whether an infection event occurs. Such an event can either lead to the infection of a susceptible individual or booster existing immunity. The time point of infection is sampled at random from the infectious individual's period of infectiousness, assuming a uniform distribution. Additional to the transmission of infection between individuals, all individuals (irrespective of their age or immunity status) are constantly exposed to a very low infection rate from the “outside” of the population. Transmission of the four influenza variants occurs independently, thus, on very rare occasions, individuals may be infected with more than one variant at the same time.

The basic reproduction number of the model is calibrated such that the median percentage of young adults (16 to 60 years of age) who are infected during the simulation year 2006/07 closely matches the observed value of 10.6% of German health care workers and controls who had a serologically confirmed influenza infection or at least a fourfold rise in influenza antibody titer during the 2006/07 season [45]. The chosen baseline value of the basic reproduction number of 1.575 yields a median infection rate of 10.6% (calculated from 2,000 simulations).

Immunity is assumed to be lost at a constant rate, leading to an exponentially distributed duration of immunity. The average duration of immunity is shorter if the immunity is caused by vaccination than when it is caused by infection. Immunity can be boostered by another infection or by a successful vaccination. If a booster event occurs, a new exponentially distributed duration of immunity is calculated for the individual. If the new duration is longer than the existing immunity, the individual’s immunity is extended accordingly. Although we use the same average length of immunity when calculating the new “booster immunity” duration, using the maximum of the existing and the new immunity duration automatically leads to a larger average duration of immunity after a booster infection than after a non-booster infection. If an individual is infected or successfully vaccinated with a virus of one of the B lineages, a cross-immunizing event can occur which protects the individual also against the other B lineage (this is determined by calculating a random number). Such cross-immunizing events can immunize susceptible individuals or booster existing immunity. In most of our analyses, we keep the percentage of B lineage cross protection constant, but in a sensitivity analysis, we allow it to decrease linearly over time, starting with a value of 60% in 2001 when B/Vic is introduced and reaching the final value of e.g. 30% on the last simulation day.

A single infection with Influenza A(H1N1) is assumed to result in an immunity which lasts on average for 6 years [31]. As we assume that new A(H1N1) drift variants appear on average every 7 years and that 60% of the old immunity is active against new drift variants (see Appendix 6 for details), immunity loss due to drift alone would, therefore, result in an average immunity duration of 7 years/(1–0.6) = 17.5 years. In order to obtain an average immunity duration of 6 years, we assume that immunity is constantly lost at a rate λ while no drift variant is introduced. The average immunity duration D _I is then given by D _I  = ((1–0.6)/7 + λ) ⁻¹ . The value λ = 1/(9.13 years) which yields an average duration of 6 years is used for all four influenza strains (due to the higher drift frequency, the average immunity duration against A(H3N2) results in about 4.5 years).

The average duration of immunity after vaccination is calculated from a figure given by Ruf et al.[44] in which the sero-protection rate (SPR) (i.e. the fraction of vaccines with a heamagglutination inhibition titer of ≥ 40) is shown for various time periods after vaccination. Values of the percentage of vaccinees who were protected t months after vaccination were read from a figure using Aotsuka’s freeware package “Data Picker” M: (Shareware Version 1.2. http://hp.vector.co.jp/authors/VA019223/). We obtain the parameters of the function p ₀  + (100-p ₀ ) _* VE _*  e ^-αt/12 by least squares fit to the picked values. The SPR of vaccinees who were immune before vaccination (p ₀ ) is estimated to be 28.6% (95% CI from 23.5 to 33.7%), the fraction of individuals who were unprotected before vaccination, yet become protected by vaccination (VE) is estimated as 92.2% (95% CI from 84.9 to 99.6%), the immunity loss rate α is estimated to be 0.553 per year (95% CI from 0.394 to 0.729). The average duration of vaccination derived immunity (in the absence of drift) is given by 1/α = 1.81 years (95% CI from 1.37 to 2.54 years).

Each newborn individual inherits maternal antibodies against all influenza strains or lineages against which his or her mother is immune at the time of birth; maternal antibodies protect the newborn individual against infection for some months. The immunological and epidemiological statuses of immigrants also need to be initialized: this is done by choosing a random individual of similar age (preferably in the year cohort of the immigrant) and transferring his or her immunity status and infection status to the immigrant. In a sensitivity analysis, we additionally consider the possibility that infection causes a short-term unspecific protection which temporarily prevents infection with any one of the influenza viruses.

The vaccination coverage depends on the age and the risk status of individuals. Furthermore, people who were vaccinated in the previous year are more likely to be vaccinated again. The baseline revaccination factor f _rev of 2 means that previously vaccinated individuals are twice as likely to be chosen for vaccination; it can be interpreted as the “relative risk” (RR) of re-vaccination given a previous vaccination. At the beginning of each simulation year, each individual is assigned a “vaccination day”, irrespective of whether this individual will be vaccinated or not. The vaccination decision will be made on that day. If a percentage c _a,r of individuals of age a and risk category r is intended to be vaccinated in a simulation year, we first calculate the expected number of vaccinees c _a,r N _a,r , whereby N _a,r is the number of individuals in the population with age a and risk status r. If there are N ⁺ _a,r previously vaccinated and N ⁻ _a,r previously unvaccinated individuals, we sample (c _a,r N _a,r ) _* (f _rev N ⁺ _a,r )/(N ⁻ _a,r  + f _rev N ⁺ _a,r ) previously vaccinated and (c _a,r N _a,r ) _* N ⁻ _a,r /(N ⁻ _a,r  + f _rev N ⁺ _a,r ) previously unvaccinated individuals to be vaccinated. During the vaccination period, it is continuously memorized how many individuals of each combination of age group, risk status and previous vaccination status have already been vaccinated and how many of them have already made the “decision” not to be vaccinated. Due to demographic changes (births, ageing, deaths, and immigrations) and due to individual changes in the risk status, individuals may change their age group or risk status during the vaccination period. If this happens, the memorized vaccination statistics is changed accordingly, so that the exact numbers of yet vaccinated and of yet unvaccinated individuals are known at any time during the vaccination period.

whereby z = N ² _a,r /(N ⁻ _a,r  + f _rev* N ⁺ _a,r ). Using the standard parameter value of f _rev  = 2, we get OR = 2/(1-c _a,r ). For a vaccination coverage of 20%, the OR is 2.5, for 60%, it is 5.0.

Smith et al.[35] constructed an antigenic map of Influenza A(H3N3) data from 35 years of surveillance (1968 to 2003) and found 11 antigenically distant clusters; the average antigenic distance between the centers of consecutive clusters was 4.5 units (SD 1.3). As influenza vaccines are updated when there is an antigenic difference of at least 2.0, these clusters can be regarded as representing genetically relevant drift variants. As ten new drift variants occurred in the time period of 35 years, we get an average duration of D = 3.5 years for new A(H3N2) drift variants (95% CI from 2.2 to 6.8). No direct estimates for the duration between consecutive Influenza B drift variants are available, but as they change at a rate which is about two to three times lower than this [3], we assume the average sojourn time of Influenza B drift variants to be 7 years. Judging by the frequency of changes in vaccine composition from 1987 to 2000, the variability of A(H1N1) may be similar to that of Influenza B [3]; we, therefore, use the same average duration of 7 years for Influenza A(H1N1) drift variants. For each one of the two A subtypes and the two B lineages, an independent random number r is calculated at the beginning of each simulation year to determine whether the circulating variant is replaced by a new drift variant. If r < 1/D, a new drift variant is introduced.

The possibility of a vaccine design mismatch is only considered for years in which new drift variants appear. Of course, TIV vaccination can only be affected by an Influenza B mismatch, if the mismatched B lineage is actually contained in the vaccine. For the last 18 years, four mismatch events have been reported (A(H3N2) in 2003 and 2004, and A(H1N1) in 2006 and 2009 [42]). For this time period, we calculate that 18/3.5 = 5.14 new drift variants for A(H3N2) and 18/7 = 2.57 new drift variants for A(H1N1) appeared. For the B lineage which was included in TIV during the 18 years, we also calculate 18/7 = 2.57 new drift variants. As we expect altogether 10.28 drift events and observed 4 mismatch events, we calculate a mismatch probability of 40% per drift. We use this probability for all four influenza variants, including the two B lineages for which no vaccine design mismatch has been reported.

The immunity shared between a previously circulating variant and a new variant is not well known. We infer information on this shared immunity from vaccination studies summarized by Jefferson et al.[39]: well-matched TIV vaccination protects on average 73% of healthy adults, whereas mismatched TIV vaccination protects only 44%. We take the ratio of 44/73 = 60% also as indicator for the shared immunity between subsequent drift variants. 60% of the individuals who were immune against the previously circulating variant are also immune against the new drift variant, whereas the remaining 40% are susceptible to the new drift variant. For the sake of simplicity, we replace variants by new drift variants at the beginning of a simulation year (i.e. on July 1st). We do not allow co-circulation of the old and the new variant in the simulation, and we assume that a previously circulating variant will never be introduced again.

We furthermore use this estimate of 60% for shared immunity between the two B lineages: if an individual is infected with one B lineage, an additional immunization (or booster event) with respect to the other lineage takes place with a probability of 60%. Likewise, TIV vaccination can cause immunity (or booster events) against the B lineage which is not contained in the vaccine; the age-dependent vaccine efficacy is multiplied by 0.6 to obtain the efficacy against the missing lineage. We restrict this cross-immunizing vaccination effect to TIV because there is no convincing evidence that increasing the antigen amount enhances the efficacy of influenza vaccines [57–59]. Tricco et. al.[10] report for vaccinees from 18 to 49 years of age that the vaccine efficacy against Influenza B is 77% if the lineage is contained in TIV and 52% if it is not contained in TIV, respectively. The ratio 52/77 = 0.675 indicates a B lineage cross protection of 67.5%. Langley et al. examine in a vaccination study with children what percentages seroconvert with respect to the B lineage which is contained and which is not contained in TIV. As they use two different TIV vaccines which either contain B/Vic or B/Yam, they obtain four different sero-conversion rates with respect to the Influenza B lineages [59] (summarized in Table 3). We have developed a maximum likelihood model based on the assumption that either one of the two TIV vaccines causes a direct sero-conversion S with respect to the B lineage contained in the vaccine and an indirect effect S x with respect to the lineage not contained in the vaccine. The joint likelihood L of the observations is then given by $\mathit{L}=\left(\begin{array}{c}\hfill 870\hfill \\ \hfill 359\hfill \end{array}\right){\left(\mathit{Sx}\right)}^{359}{\left(1-\mathit{Sx}\right)}^{511}\left(\begin{array}{l}870\\ 622\end{array}\right){\mathit{S}}^{622}{\left(1-\mathit{S}\right)}^{248}\left(\begin{array}{l}877\\ 644\end{array}\right){\mathit{S}}^{644}{\left(1-\mathit{S}\right)}^{233}\left(\begin{array}{l}877\\ 262\end{array}\right){\left(\mathit{Sx}\right)}^{262}{\left(1-\mathit{Sx}\right)}^{615}$, omitting the constant factors, this simplifies to L ∝ (Sx)^{359 + 262}(1 − Sx)^{511 + 615} S ^{622 + 644}(1 − S)^{248 + 233}.