Seasonal influenza: Modelling approaches to capture immunity propagation

Seasonal influenza poses serious problems for global public health, being a significant contributor to morbidity and mortality. In England, there has been a long-standing national vaccination programme, with vaccination of at-risk groups and children offering partial protection against infection. Transmission models have been a fundamental component of analysis, informing the efficient use of limited resources. However, these models generally treat each season and each strain circulating within that season in isolation. Here, we amalgamate multiple data sources to calibrate a susceptible-latent-infected-recovered type transmission model for seasonal influenza, incorporating the four main strains and mechanisms linking prior season epidemiological outcomes to immunity at the beginning of the following season. Data pertaining to nine influenza seasons, starting with the 2009/10 season, informed our estimates for epidemiological processes, virological sample positivity, vaccine uptake and efficacy attributes, and general practitioner influenza-like-illness consultations as reported by the Royal College of General Practitioners (RCGP) Research and Surveillance Centre (RSC). We performed parameter inference via approximate Bayesian computation to assess strain transmissibility, dependence of present season influenza immunity on prior protection, and variability in the influenza case ascertainment across seasons. This produced reasonable agreement between model and data on the annual strain composition. Parameter fits indicated that the propagation of immunity from one season to the next is weaker if vaccine derived, compared to natural immunity from infection. Projecting the dynamics forward in time suggests that while historic immunity plays an important role in determining annual strain composition, the variability in vaccine efficacy hampers our ability to make long-term predictions.

The data compared the ILI rates of people with chronic diseases (Chronic Disease = 1) against people without chronic diseases (Chronic Disease = 0).

Contributors:
1. Simon de Lusignan Director, guarantor for these data, assisted with clinical knowledge 2. Rachel Byford Put together structures for the data extraction.

Ana Correa
Carried out the quality assessment of the output. 4. Chris McGee Carried out data extraction of ILI rates (covering week 1 2000 to week 52 2016). 5. Julian Sherlock Carried out data extraction of ILI rates and data quality checks (covering week 1 2017 to week 35 2018). 6. Sameera Pathirannehelage Carried out data extraction of ILI rates (covering week 1 2017 to week 52 2018). 7. Ivelina Yonova liaison with practices, and review of the original data request.

Acknowledgements:
Practices who have agreed to be part of the RCGP RSC and allow us to extract and used health data for surveillance and research. Filipa Ferreira (programme manager), and other members of the Clinical Informatics and Health Outcomes Research Group at University of Surrey. Apollo Medical Systems for data extraction. Collaboration with EMIS, TPP, In-Practice and Micro-test CMR supplier for facilitating data extraction. Colleagues at Public Health England.

Respiratory Virus RCGP Surveillance
We sourced data pertaining to the weekly percentage of sentinel virology samples that were influenza positive from a subset of general practices in the RCGP Weekly Returns Service, which submitted respiratory samples for virological testing from patients presenting in primary care with an ILI.
These data were displayed in figures within PHE annual influenza reports (for the 2009/10 to 2017/18 influenza seasons inclusive) [5,6], and PHE Weekly National Influenza reports (for the 2013/14 influenza season onwards) [7]. . As an exception, the data curve for the entire 2009/10 influenza season was available within the 2010/11 PHE annual influenza report [5].
In all seasons, we assumed weeks 36-39 and weeks 21-39 had an influenza positivity of 0%.
To inform positivity values for weeks preceding week 20 for the 2013/14 influenza season onwards that were not illustrated within the PHE annual influenza report, we used PHE Weekly National Influenza reports [7], which provided United Kingdom GP sentinel swabbing scheme sample positivity on a weekly basis. Note that we gave precedence to figures explicitly stated within annual reports over positivity quantities within the weekly PHE influenza summaries.

Circulating strain composition
Here, we expand on our use of FluNet [8] to inform the circulating influenza strain distribution per influenza season.
FluNet is a global web-based tool for influenza virological surveillance first launched in 1997. The virological data entered into FluNet, e.g. number of influenza viruses detected by subtype, are critical for tracking the movement of viruses globally and interpreting the epidemiological data.
FluNet reports weekly influenza surveillance data that are provided from over 140 National Influenza Centres of the Global Influenza Surveillance and Response System, national influenza reference laboratories, and WHO regional databases [9]. Per the WHO Global Epidemiological Surveillance Standards for Influenza guidance, influenza testing is conducted on specimens collected from persons presenting for medical care at participating surveillance sites who meet a clinical definition for influenza-like illness (defined as an acute respiratory infection with measured fever of ≥ 38 • C, and cough, with onset within the last 10 days), or severe acute respiratory infection (defined as an acute respiratory infection with history of fever or measured fever of ≥ 38 • C, and cough, with onset within the last 10 days, and requiring hospitalization) [9]. Influenza is confirmed by accepted laboratory diagnostic methods and actively reported by the reference laboratory to FluNet [10].
We computed the empirical circulating strain distribution in each influenza season (from the 2009/10 influenza season onward) using data from FluNet for the United Kingdom. A salient feature of the influenza type data is the presence of samples whose exact subtype/lineage were not determined (Fig. A). We assumed the fraction of undetermined samples ascribed to each subtype/lineage matched the proportions observed for the set of samples where strain-specific information were available (Fig. B) (a)

Vaccine uptake
For the 2009/10 influenza season, we acquired vaccine uptake profiles from survey results on H1N1 vaccine uptake amongst patient groups in primary care [11].
For the 2010/11 influenza season onward, we sourced seasonal influenza vaccine uptake figures from Public Health England official statistics [6,7]. Taken from PHE Weekly National Influenza Reports, these data comprised uptake profiles per target vaccination age group (e.g. at-risk under 65 years, 65 year and over, 2 year-olds, 3 year-olds) for each influenza season. Populationlevel measures were acquired through weighting each uptake curve by the respective population proportion. Age-partitioned population distribution estimates for England were obtained from the Office of National Statistics (ONS) [12,13].
All age, overall population influenza vaccine coverage has been on the rise, increasing from 21.0% to 26.0% between the 2009/10 and 2017/18 influenza seasons (Fig. C).
We constructed daily uptake rates from either weekly or monthly uptake values via linear interpolation.
The temporal resolution of the empirical vaccine uptake data was dependent upon the age category, stratified by: (i) Non-school aged children (weekly coverage figures); (ii) school-aged children (monthly coverage figures). Vaccination coverage (%) constructed based upon GP practices reporting weekly to Immform (the influenza vaccine uptake monitoring programme from Public Health England).
For the groupings of at-risk under 65 years and 65 years and above, vaccine uptake profiles were available from the 2012/13 influenza season onwards. For the 2010/11 and 2011/12 influenza seasons we invoked the uptake profile from the 2012/13 influenza season.
Collection of cumulative weekly vaccine uptake for 2 year-olds and 3 year-olds started in 2013/14.
(ii) School-aged children For school-aged children, monthly vaccination uptake estimates were available. These uptake estimates specified the proportion of children in England who received the influenza vaccine via school, pharmacy or GP practice.

Vaccine efficacy
We informed population level vaccine efficacy estimates for each historical influenza season from publications detailing end-of-season age adjusted seasonal influenza vaccine effectiveness for adults and children in preventing laboratory-confirmed influenza in primary care in the United Kingdom [14][15][16][17][18][19][20]. In influenza seasons where an equivalent publication were not available, we used mid-season or provisional end-of-season age adjusted seasonal vaccine efficacy estimates from Public Health England reports [21,22].
The 2009/10 influenza season was an exception. Due to the adjusted seasonal influenza vaccine efficacy being -30% (-89%, 11%) [14], we instead used the pandemic vaccine to inform efficacy against A(H1N1)pdm09, with the effectiveness against all other types set to zero.
The all age, strain-specific adjusted vaccine efficacy estimates gathered from the literature are displayed in Table A. In the United Kingdom, since the incremental introduction of the universal childhood influenza vaccine programme began in the 2013/14 influenza season, an intra-nasally administered live attenuated influenza vaccine (LAIV) has been administered to children. Adult age classes have generally been offered an inactivated influenza vaccine (IIV). Therefore, due to using a population-averaged efficacy measure, in influenza seasons where IIV and LAIV formulations were administered to different age classes and/or risk groups the overall efficacy estimate is a combination of IIV and LAIV effectiveness.
With the empirical data not providing individual estimates for each influenza A subtype and influenza B lineage, we invoked a series of assumptions to produce the strain-specific vaccine efficacy quantities used within our study (Table 1).
We outline herein the set of enacted assumptions.
Approach to address absent influenza A efficacy data When there were estimates present for overall influenza A efficacy, but estimates for only one of the A(H1N1)pdm09 or A(H3N2) subtypes, for the subtype with no efficacy data we assumed its efficacy matched the overall influenza A efficacy estimate.
Case two: Data absent for one subtype and overall influenza A efficacy Applicable to the following influenza seasons: 2010/11, 2011/12.
In those influenza seasons where there was an efficacy estimate for one of the two influenza A subtypes, but efficacy estimates for the other influenza A subtype and against influenza A overall were absent, we assumed the efficacy for the subtype with no relevant empirical data available was equal to that of the other influenza A subtype (for which there was an available estimate). In those influenza seasons where there was an overall influenza B efficacy estimate, but efficacy estimates for both influenza B lineages were absent, we assumed the vaccine efficacy for both lineages matched the overall influenza B efficacy estimate. : The adjusted seasonal influenza vaccine efficacy in the 2009/10 influenza season was -30% (-89%, 11%) [14]. We therefore used the pandemic vaccine to inform efficacy against A(H1N1)pdm09, with the effectiveness against all other types set to zero. †: Mid-season estimate of seasonal influenza vaccine effectiveness from [21]. Low incidence throughout the 2013/14 influenza season meant reliable end-of-season estimates for the vaccine efficacy could not be attained (Personal communication, Public Health England). ‡: Provisional end-of-season influenza vaccine effectiveness results. ' 2 Complementary details of the modelling approach 2.1 Epidemiological model specifics  Compartments situated within the shaded region signify vaccinated groups, with vaccination occurring at rate µ. Transmission events lead to movement from a susceptible state (S) to latent (E,infected but not yet infectiousness). Latency is lost at rate γ 1,m . Infected (I) transition to recovered (R) at rate γ 2 . Equivalent formulation for each strain m.

Approach to address absent influenza B efficacy data
The epidemiological model is a deterministic, non-age, multi-strain structured compartmental based model capturing demography, influenza infection status (with susceptible-latent-infectedrecovered, SEIR, dynamics) and vaccine uptake. In addition, we assumed disease transmission to be frequency-dependent.
We let E X m , I X m and R X m denote the proportion of the population that have vaccination status X ∈ {N, V } (indexed by N for non-vaccinated, V for vaccinated), and that are latent, infectious and recovered (as a result of natural infection) with respect to strain m.
Exposure to influenza virus in the previous influenza season, through natural infection or vaccination, modulated current influenza season susceptibility. Tracking immunity derived from natural infection and vaccination separately required ten distinct exposure history groups, with the susceptibility to strain m for a given exposure history group h encoded into the susceptibility array f (h, m) (Fig. 3).
To incorporate exposure history from the previous influenza season, we let S X,h denote the proportion of the population that have (current influenza season) vaccination status X ∈ {N, V } and exposure history h that are susceptible to all strains. With ten exposure history groups h in place, a total of 20 susceptibility states were used: ten S N,h states, accounting for susceptibles not vaccinated in the current influenza season stratified by exposure history grouping; ten S V,h states, tracking susceptibles who have been administered the vaccine in the current influenza season whilst retaining exposure history group information.
The epidemiological model has the following formulation, with time dependencies dropped: where γ 1,m is the strain-dependent rate of loss of latency, γ 2 the rate of loss of infectiousness, α m the strain-specific vaccine efficacy, and µ = with ν the rate of vaccination in the population.
To maintain a constant population size, birth and death rates were equal (B = D). We set the mortality rate based upon a gender-averaged life expectancy for newborns in England, computed from male and female specific statistics from ONS data [23]. We computed the gender-averaged life expectancy for newborns in England to be approximately 81 years, with the corresponding daily mortality rate being D = 1 81×365 . The strain specific force of infection, λ m , satisfies λ m = β m X I X m , where β m is the transmission rate for strain m, implicitly comprising the contact rate and the transmissibility of the virus (the probability that a contact between an infectious person and a susceptible person leads to transmission).

Between-season exposure history group mappings
At the beginning of each influenza season (1st September), we enacted the following exposure history group assignments: where X represents vaccination status, with N corresponding to non-vaccinated and V vaccinated (bars on vaccination state symbols relate to vaccination status in the previous influenza season).

Parameter inference
The aim of our parameter estimation was to select parameter sets generating influenza-attributed ILI GP consultation model predictions that best resemble the empirical data.
In full, we sought to estimate: transmissibility of each influenza virus strain (q m ), modified susceptibility to strain m given infection by a strain m type virus the previous influenza season (a), carry over cross-reactivity protection between influenza B lineages (b), proportion of prior influenza season vaccine efficacy propagated to the current influenza season (ξ), and an ascertainment probability per influenza season ( t ).
We selected an Approximate Bayesian Computation (ABC) approach, where results are first simulated with a parameter set, after which it is rejected or accepted based on a divergence measure or statistic.
We sought parameter sets that: (i) adhered to our infection prevalence temporal profile check; the peak in influenza infection could not occur after February (i.e. during March through August inclusive) in any influenza season; (ii) minimised our deviance statistic, measuring the correspondence of the model predicted influenza-attributed ILI GP consultations versus the observed data (Equation 1).
To update the parameter sets we employed an adaptive population Monte Carlo ABC algorithm [24], an ABC scheme combined with particle Monte Carlo methods to minimise the number of simulations.
Outlining the inference scheme in brief, we initially drew 10,000 particles (i.e. 10,000 parameter sets) from the prior distributions, simulated and calculated their statistics. From these firstgeneration particles, each with equal weight, we generated 10,000 additional particles. Next, we sampled new particles according to the weight of the first generation particles and perturbed according to a multivariate normal distribution with twice the weighted variance-covariance matrix (computed using analytic weights) of the first generation particles. After simulation of the new particles, we computed their statistics and new weights. Of the now 20,000 particles, we retained the 10,000 particles with the lowest statistics.
The procedure was repeated iteratively, except all subsequent perturbations were carried out with a local perturbation kernel. Local perturbation kernels assign a different covariance matrix for each particle of the previous population. The class of perturbation kernel we selected was the multivariate normal kernel with optimal local covariance matrix (OLCM) [25]. We ran the process for at least 600 generations.
The final generation of 10,000 particles represents a sample from the posterior distribution. To compute posterior parameter distribution summary statistics, parameter values were sorted in ascending order, x 1 , . . . , x 10,000 , with associated normalised weightsw 1 , . . . ,w 10,000 (sum of all weights equal to one). The weighted median was defined as the element x k such that: The 95% credible intervals were calculated as We additionally utilised the retained particles for forward simulations.

Fitting to subsets of the empirical data
To delve into model robustness, we sought parameter fits using subsets of the historical influenza season data. Categorically, we fit to two datasets. The first covered four influenza seasons spanning 2012/13-2015/16. The second also included the 2016/17 influenza season, bringing the total number of influenza seasons to five.
Invoking the adaptive population Monte Carlo ABC scheme upon each dataset, we obtained 10,000 parameter sets following the completion of 1,100 and 1,000 generations respectively (Figs. L and M). To view detailed summary statistic per parameter, we refer the reader to Table B.
We witness agreement with the statistical fit to the full dataset on the following aspects : (i) influenza A transmissibility exceeding type B virus transmissibility; (ii) amongst the three mechanisms for immunity propagation, retained natural infection in the previous influenza season having the greatest influence; (iii) little support for immunity propagation originating from influenza B cross-reactivity (b attaining values near 1); (iv) quantitatively similar ascertainment probabilities.
Whilst outcomes fitting to a range of influenza seasons generally exhibited qualitative consistency, there were notable quantitative differences. Transmissibility magnitudes were raised, with the transmissibility of the A(H3N2) subtype exceeding that of the A(H1N1)pdm09 subtype (whereas in the complete fit we inferred similar transmissibility values for the type A viruses). The elevated transmissibility levels were counteracted by the enlarged effect of prior infection and vaccination reducing susceptibility. Explicitly, while the parameter distributions inferred fitting to the period 2012/13-2017/18 indicated little residual vaccine immunity and an estimated 20% drop in strain-specific susceptibility as a result of natural infection in the previous influenza season, when fitting to the shorter time window, susceptibility at the outset of an influenza season against a given strain type saw a near 55% reduction if naturally infected in the prior influenza season by a virus of that class, and roughly 40% of prior influenza season vaccine efficacy was retained (Table B).
In a similar fashion to the dataset encompassing influenza seasons ranging from 2012/13 up to and including 2017/18, we found strong, qualitative correspondence between simulations of the paramterised model and the data (Figs. N and O).

Parameter identifiability
We generated and used a synthetic dataset to verify the parameter identifiability capabilities of our ABC inference scheme. To mimic the format of our empirical data, the synthetic data corresponded to strain-specific influenza positive ILI GP consultations rates per 100,000 population from the 2012/13 to 2017/18 influenza season (inclusive). The parameter values we used for producing the synthetic are listed in Table C.
When fitting the model to the synthetic data, we employed our inference scheme in an equivalent manner as for fitting to the empirical data. We ably recovered the parameter values from which the synthetic data had been generated (Fig. P and Table C). Furthermore, simulated outcomes from each parameter set had strong agreement with the synthetic data (Fig. Q).   Stratified by influenza season, we present back-to-back stacked bars for 1,000 simulation replicates, each using a distinct parameter set representing a sample from the posterior distribution when fitting our mathematical model to the synthetic data. Each influenza season is topped by a thicker stacked horizontal bar plot, corresponding to the synthetic data values. The left side depicts the cumulative total of ILI GP consultations attributable to type A influenza per 100,000 population (red shading denoting the A(H1N1)pdm09 subtype, orange shading the A(H3N2) subtype). In an equivalent manner, the right side stacked horizontal bars present similar data for type B influenza (cyan shading denoting the B/Victoria lineage, dark blue shading the B/Yamagata lineage). Overall, simulated outcomes from each parameter set had strong agreement with the synthetic data.  Fig R. Quantification of the projected dominating influenza A subtype and influenza B lineage by influenza season, up to 2029/30. The left half depicts the proportion of simulations in which the major subtype contributor (>50% of cases) to the overall type A influenza incidence was A(H1N1)pdm09 (red bars) or A(H3N2) (orange bars). In an equivalent manner, the right half depicts the proportion of simulations in which the major lineage contributor (>50% of cases) to the overall type B influenza incidence was B/Victoria (cyan bars) or B/Yamagata (dark blue bars). Constructed from simulation runs where vaccine efficacy against each strain were randomly sampled from the empirical distribution (totalling 1,000 replicates). For all forward simulated influenza seasons, vaccine uptake matched that of the 2017/18 influenza season.

Model extension: Immunity propagation across multiple influenza seasons
In our extended model, we altered the immunity propagation model component. The vaccination, epidemiological and observation model components were all unchanged.
We display a modified directed acyclic graph, which includes the revised immunity propagation component structure and incorporation of the data streams (Fig. S). Dotted arrows indicate relationships between prior influenza season epidemiological outcomes and immunity propagation factors. Solid arrows indicate within season processes. Circled capitalised letters indicate the relationships connecting the variables or data involved. These relationships are: process A, propagation of immunity as a result of exposure to influenza virus in the previous influenza season (through natural infection or vaccination); process B, modulation of current influenza season virus susceptibility; process C, estimation of influenza case load via the SEIR model of transmission; process D, ascertainment of cases through ILI recording at GP.

Details of the expanded immunity propagation mechanism
As before, we tracked immunity derived from natural infection and vaccination separately. We kept our three, original susceptibility modifying factors: (i) modified susceptibility to strain m given infection by a strain m type virus the previous influenza season, denoted a; (ii) carry over cross-reactivity protection between influenza B lineages, denoted b (as infection with one influenza B virus lineage may be beneficial in protecting against subsequent infection with either influenza B virus lineage [27]); (iii) residual strain-specific protection carried over from the prior influenza season influenza vaccine, denoted c m . We again mandated that 0 < a,b,c m < 1.
To enhance the flexibility of the immunity propagation component, we introduced and fit an additional parameter, δ, representing the proportion of those who both began the influenza season in an exposure history group linked to natural infection (Fig. 3: rows 2-5, 7-10) and were unexposed to influenza virus during the current season (who at the end of the current influenza season remained susceptible), that retained their pre-existing natural infection acquired immunity. Those keeping immunity arising from natural infection were mapped to the relevant prior infection exposure history group dependent upon current influenza season vaccination status. The remaining proportion (1 − δ) transitioned in the same manner as in the original model, reverting to either the naive exposure history group or vaccinated only exposure history group (Fig. T).
For clarity, there was no mechanism in the model framework to confer vaccine-induced immunity beyond one influenza season (i.e. vaccine-induced immunity could be retained for, at most, a single additional influenza season). As before, the collection of exposure history groupings and associated strain-specific susceptibilities were consolidated into a single susceptibility array (Fig. S, process B; Fig. 3).

Between-season exposure history group assignments
At the beginning of each influenza season (1st September), we enacted the following exposure history group assignments: At the beginning of a generic influenza season y, those susceptible may be divided into two categories based on exposure history: (i) susceptibles in either the naive exposure history group , with unmodified susceptibility to infection, or in the vaccinated only exposure history group (labelled S); (ii) all those that are susceptible and in an exposure history group involving natural infection (labelledŜ). During the influenza season, those that were initially susceptible may be infected (state I). Of those remaining in stateŜ at the end of influenza season y, a proportion δ are retained in stateŜ at the beginning of influenza season y + 1. The remainder (1 -δ) transition to state S. Note, δ = 0 corresponds to no residual immunity being retained due to natural infection beyond one subsequent influenza season, while δ = 1 coincides with residual immunity being kept in all future influenza seasons until the occurrence of a new infection event.

Fitting the extended model to the data
We performed parameter inference using an equivalent ABC scheme to that used to fit the original model. For the extended model, we fit the following collection of parameters (Table D): transmissibility of each influenza virus strain (β m ), modified susceptibility to strain m given infection by a strain m type virus the previous influenza season (a), carry over cross-reactivity protection between influenza B lineages (b), residual protection carried over from the prior influenza season vaccine (ξ), the fraction unexposed in the prior influenza season that retain immunity (δ), and an ascertainment probability per influenza season ( y ).
We amassed 10,000 parameter sets (representing a sample from the posterior distribution) to determine credible values of the model parameters (Table E). Samples were obtained after completion of 750 generations of the adaptive-population Monte Carlo ABC scheme, at which point the generation-by-generation tolerance level updates were minor ( Fig. U(a)). The error values attained for the parameter sets inferred for the extended model closely matched the error values for the parameter sets inferred for the original (less complex) model. The inferred model parameters for the extended model were well defined with Gaussian-shaped histograms ( Fig. U(b)).
Comparing the parameter fits to those for the original model, we maintained the main outcome of propagation of immunity being weaker if vaccine derived, compared to natural immunity from infection. As before, we inferred little carry over of prior influenza season vaccine efficacy, with the majority of the posterior distribution for ξ massed near 0.
We found several other similarities amongst the transmission parameters and ascertainment probabilities when comparing the two model fits. Virus transmission rates were larger for the two influenza A subtypes than the corresponding estimates for the two type B lineages, and B/Yamagata transmissibility estimates exceeded those for B/Victoria. Furthermore, for each Discrepancies did arise amongst parameters tied to immunity processes. First, a decrease in value for parameter a, corresponding to the modified susceptibility to a given influenza type that results from natural infection by a prior influenza virus of that type (median values, extended model vs original model: 0.7168 vs 0.7883). These predictions correspond to an almost 30% reduction in susceptibility in the extended model compared to little more than 20% in the original model. Second, posterior distributions for b imply a similar amount of modulation to susceptibility via influenza B cross-reactivity. This is in stark contrast to the less complex model, where there was little support for carry over of type B influenza cross-reactivity from one influenza season to the next (median values: 0.7214 vs 0.9703).
Inspecting the posterior distribution for δ, the majority of samples attained values greater than 0.9. The inferred magnitude implies long term retention of immunity acquired from natural infection (assuming no further vaccination or infection event to alter the immunity landscape), in agreement with prior studies supporting infection acquired immunity waning over a time scale exceeding a year [30][31][32]).
Performing 1,000 independent simulations of the extended model with parameter sets drawn from the ABC inference procedure, we found a similar fit to the data with the extended model (Fig. V) compared to the less complex, original model (Figs. 5-6). Thus, despite the added complexity, the extended model did not give noticeable improvements in correspondence of model outputs with the data (relative to fits with a model using the simpler immunity propagation setup).

Parameter identifiability
We again generated and used a synthetic dataset to verify the parameter identifiability capabilities of our ABC inference scheme under the revised model framework. The synthetic data corresponded to strain-specific influenza positive ILI GP consultations rates per 100,000 population from the 2012/13 to 2017/18 influenza season (inclusive), to mimic the format of our empirical data. The parameter values we used for producing the synthetic are listed in Table F.
Once more, we ably recovered the parameter values from which the synthetic data had been generated ( Fig. W and Table F). Furthermore, simulated outcomes from each parameter set had strong agreement with the synthetic data (Fig. X).   Posterior predictive distributions for influenza positive GP consultations per 100,000 population, using the parameter sets generated fitting the extended model to the synthetic data. Stratified by influenza season, we present back-to-back stacked bars for 1,000 simulation replicates, each using a distinct parameter set representing a sample from the posterior distribution when fitting our mathematical model to the synthetic data. Each influenza season is topped by a thicker stacked horizontal bar plot, corresponding to the synthetic data values. The left side depicts the cumulative total of ILI GP consultations attributable to type A influenza per 100,000 population (red shading denoting the A(H1N1)pdm09 subtype, orange shading the A(H3N2) subtype).
In an equivalent manner, the right side stacked horizontal bars present similar data for type B influenza (cyan shading denoting the B/Victoria lineage, dark blue shading the B/Yamagata lineage). Overall, simulated outcomes from each parameter set had strong agreement with the synthetic data.