Forecasting the 2017/2018 seasonal influenza epidemic in England using multiple dynamic transmission models: a case study

Background Since the 2009 A/H1N1 pandemic, Public Health England have developed a suite of real-time statistical models utilising enhanced pandemic surveillance data to nowcast and forecast a future pandemic. Their ability to track seasonal influenza and predict heightened winter healthcare burden in the light of high activity in Australia in 2017 was untested. Methods Four transmission models were used in forecasting the 2017/2018 seasonal influenza epidemic in England: a stratified primary care model using daily, region-specific, counts and virological swab positivity of influenza-like illness consultations in general practice (GP); a strain-specific (SS) model using weekly, national GP ILI and virological data; an intensive care model (ICU) using reports of ICU influenza admissions; and a synthesis model that included all data sources. For the first 12 weeks of 2018, each model was applied to the latest data to provide estimates of epidemic parameters and short-term influenza forecasts. The added value of pre-season population susceptibility data was explored. Results The combined results provided valuable nowcasts of the state of the epidemic. Short-term predictions of burden on primary and secondary health services were initially highly variable before reaching consensus beyond the observed peaks in activity between weeks 3–4 of 2018. Estimates for R0 were consistent over time for three of the four models until week 12 of 2018, and there was consistency in the estimation of R0 across the SPC and SS models, and in the ICU attack rates estimated by the ICU and the synthesis model. Estimation and predictions varied according to the assumed levels of pre-season immunity. Conclusions This exercise successfully applied a range of pandemic models to seasonal influenza. Forecasting early in the season remains challenging but represents a crucially important activity to inform planning. Improved knowledge of pre-existing levels of immunity would be valuable.


B Appendix: Model Details
In what follows, all references to weeks will be to the ISO standard definition for week numbers, and a shorthand will be employed, where 2018w4 refers to ISO week 4 of 2018.

B.1 Stratified Primary Care (SPC) model
The stratified primary care model has been presented elsewhere as the parallel-region model of Birrell et al (2016). In summary, a deterministic model is used to approximate influenza transmission resulting in a time series of daily numbers of new infections. These infections then link into a disease reporting model which accounts for the proportion of the infections that are observed through consultation with a general practitioner (GP) for influenza-like illness (ILI) as well as the time period from infection to the reporting of the GP consultation. Data on GP ILI consultations are not disease specific and therefore are 'contaminated' with consultations for other infections. These noninfluenza related consultations are referred to as background ILI cases. Virological (swab) positivity data is therefore required to understand the fraction of the daily ILI consultations that are attributable to influenza.

B.1.1 Transmission Dynamics
The dynamics are as detailed in Birrell et al (2016) with the exception that we only work here with a single age-group (due to insufficient information on age-group prior immunity). The model is stratified into five regions: Greater London, North (combining NHS regions North West, North East and Yorkshire and Humberside), Midlands and the East (combining the NHS regions Midlands and East of England), South East and South West. For a map of the regions see https://www.hee.nhs.uk/about/how-we-work/your-area.
Here U V W is a school-holiday effect equal to 1 when # $ is a time outside of a school holiday period and U during school holidays, X`, 2 is the basic reproductive number, and Z 2 is the population size, both in the region 0. The initial conditions for the system of equations (1) are detailed in Birrell et al (2016), but it bears pointing out that the initial pool of susceptible individuals is given by 1 2 (# Y ) = aZ 2 , where parameter a gives the initial fraction of the population that are susceptible to infection. Throughout, a value of a = 0.9 will be used to represent a belief that at most people will be susceptible to at least one of the circulating influenza viruses.
In each time interval [# $%& , # $ ), a number of new infections is calculated as forming the input to the disease reporting model.

B.1.2 Disease Reporting Model
The disease reporting component of the SPC model accounts for • the fraction of infections that lead to a GP consultation; • the delay time from infection to the reporting of the GP consultation.
• the level of GP consultations for ILI not directly attributable to influenza infection.
It is assumed that a fraction of the infections, e, will develop ILI symptoms. Of those who develop symptoms, it is assumed that a further fraction, f gh , will consult with their GP. Where a consultation occurs, the time between infection and the reporting of the consultation is governed by a discretised gamma distribution, with parameters i and j, such that k(l; i, j) gives the probability that the time from infection to the reporting of the healthcare event is l.# (assuming that infections and disease reporting occur at the same point in each time interval). The number of GP consultations attributable to influenza infection is therefore given by In 2018w1, when we only have a single region and do not directly incorporate the virological data into the analysis, the ILI+ data are assumed to be imperfect observations of the v 2 RPR^( # $ ), and therefore are assumed to be distributed with (3) as their mean.
However, from 2018w2 onwards, the full GP ILI consultation data are used with virological swabbing data. As the data arrive daily, there is also a strong within-week pattern to the consultations: there is little to zero reporting of consultations at weekends, whilst many will delay their consultation until the following Monday, the busiest day of the week.

B.1.2.1 Background Consultations
The background consultations give a level of GP consultation for ILI that would exist in the absence of influenza activity. These consultations are assumed to be highly seasonal, with a peak in midwinter. To account for this pattern, we assume a model for Ö 2 (Ü), the weekly rate of non-influenza ILI consultations in region 0 during week Ü, that can be expressed as: the type of seasonal model used for the endemic component in the HHH model of Held, Höhle and Hofmann, 2005. All three parameters, â 2 , ä 2 and å 2 , will be estimated, but will be done so with informative priors placed on them. These informative priors were found by taking the posteriors obtained when fitting the epidemic/endemic HHH model to GP In-Hours ILI consultation data over the period 12/1/2015 to 1/10/2017. If ê 2,ë is the GP ILI consultation data in week Ü, then: where Z 2,ë is the catchment (observable) population in region 0 in week Ü, an average of the daily catchment populations, and ï 2 is the region-specific parameter of the epidemic component of the HHH model.

B.1.3 Likelihood
Denote • û 2,V W to be the random variable representing the observed number of ILI consultations in region 0 at time # $ , with observed value ü 2,V W . • † 2,V W to be the random variable representing the observed number of swabs testing positive for the presence of influenza out of a sample of size °2 ,V W in region 0 at time # $ , with observed value ¢ 2,V W .
It is assumed that the ILI consultation data are distributed according to the negative binomial distribution: where the Z 2,V W is the size of the catchment population on day # $ , and ¶ 2 is a dispersion parameter such that varßû 2,V W ® = (1 + ¶ 2 )íû 2,V W .
Similarly, the swab positivity data are distributed † 2,V W~ Bin ö°2 , These distributional assumptions help us to compile a likelihood. If we denote Φ to be a vector of all the model parameters, then the likelihood is For ease of presentation, equations (4)-(6) assume that the data arrive at intervals of .#. In practice, data are daily, whereas .# = 1 2 ⁄ days. In which case, in (4) and (5) the v 2 MLM^( # $ ), Ö 2 (# $ ) and v 2 (# $ ) are summed to give daily totals, and the product in (6) is over days, not intervals.

B.1.4 Model Evolution
Throughout the influenza season, there were changes in the available information (see Appendix A: Data). Therefore, as with the other models, some increase in the model complexity and detail was possible over time. The table below indicates how this changed for the SPC model.

B.2 Strain Specific SIR model
The model used here has an SIR structure, fitted to both the weekly ILI and virological data from RCGP. It is assumed that there are three circulating strains each acting independently.
The SIR model is specified in continuous time by the system: where 1 ≠ (#) represents the number of individuals susceptible to strain j at time #. Similarly, : ≠ (#) and X ≠ (#) represent the numbers of infected and recovered individuals respectively at time # and for strain j. Z is the total population size (Z = 1 ≠ + : ≠ + X ≠ ), AE ≠ is the strain specific rate of transmission and ; is the daily rate of recovery. In practice, this system of equations is solved using a first-order Euler approximation with time-steps of .# = 1 day.
We denote the number of ILI attributable to influenza strain j in week Ü as :™:^, ≠ (Ü), and the corresponding level of non-influenza ILI as :™: % (Ü). The total number of ILI cases in week Ü is then the sum of the two components i.e. :™:(Ü) = :™: % (Ü) + ∑ :™:^, ≠ (Ü) ≠ , where these components result from with e ≠ the daily rate of becoming symptomatic due to influenza and e % the corresponding rate of developing ILI symptoms when not infected with influenza.
As in the SPC model, the background influenza-negative ILI rates increase during the winter (Fleming & Elliot, 2008). To account for this, e % was assumed to change over time as follows: where e ≥ is the maximum value of the (log) influenza negative ILI rate, ¥ is the amplitude of oscillation, τ gives the timing of the peak and ª governs the period of fluctuation. Note that the rate of influenza positive ILI is assumed to be proportional to the rate at which individuals recover from infection (see first line of Equation (8)). This introduces a delay between becoming infected and developing symptoms, visiting the GP and being diagnosed. For simplicity, we assumed that this delay is on average the same as the delay between becoming infected and recovery.

B.2.1 Likelihood
Again, assume that the observed number of ILI consultations, ü ë MLM , is the realisation of a random variable û ë MLM , with the following distribution: where the weekly catchment population is Z ë , and º is the probability with which someone with ILI is diagnosed, i.e. this is a combination of the probability that the infected cases consult the GP and the GP correctly diagnoses ILI.

B.2.2 Model evolution
Two different models have been fitted over time as new data became available. The table below describes the scenarios considered and the dates from which the analyses were carried out. Scenario 2 involves a time-dependent influenza negative ILI rate, while in Scenario 1 the influenza negative ILI rate was kept constant.

2018-01-28
Figure S2 Estimated number of infected over time for the different strains by analysis date and type of model (see Table  above).

B.3 ICU model
The ICU model was initially formulated to analyse data collected by USISS during the 2012/13, 2013/14 and 2014/15 influenza seasons. For a detailed description of the model we refer to Corbella et al. (2018) and we report below only the key elements and assumptions.

B.3.1 Transmission model
We assume that the transmission dynamics can be approximated by a continuous-time deterministic SEEIIR model: where 6 and ; relate to the mean latent and infectious periods respectively as in Equation (2).
The transmission rate AE(#) is assumed to change over time as follows: The factor U ∈ (0,2) allows for the possibility that rates of infectious contact might change as a result of school holidays. Beyond this, it is assumed that contact and transmission are homogeneous over age, geographical region and influenza sub-type. At 2018w8 a number of alternative models were considered, models assuming: • a constant transmission rate: AE(#) = AE Y • only a Christmas holiday effect: otherwise. • transmission dependent on climate variables: AE(#) assumed to be a linear function of either the daily absolute humidity or of the daily temperature, denoted below by ò(#) The competing models were used to analyse all the data up to 2018w8 inclusive. The Deviance Information Criterion (DIC) (Spiegelhalter et al, 2002) was then used to compare the possible models. The best-fitting model according to the DIC was the school-holiday model (Equation (11)). The other models were not considered further.

B.3.2 Reporting model
Denote by: • ∆(Ü) the number of new infections during week Ü; • f Mƒ≈ the probability of Intensive Care Unit (ICU) admission given infection; • k Mƒ≈|M«» (…) the probability of … weeks elapsing between infection and ICU admission; • K(Ü) the probability of detection given ICU admission at week Ü.

B.3.3 Inference
We assume the latent and infectious periods have fixed means K L = 2 and K M = 3, respectively. Moreover, we assume that the weekly probability of detection K(Ü) is known and equal to the known fraction of the total population within the catchment of reporting hospitals that reported data over the total monitored population, i.e., where Z ë is the total population within the catchment of reporting hospitals in week w and N is the full population of England.
Like the other modelling approaches, inference is carried out within the Bayesian paradigm: prior distributions for model parameters are assumed to be uniform within some plausible limits, as described in Corbella et al. (2018), with the exception of the initial susceptibility, a, and f Mƒ≈ which are distributed: a~Beta(5, 45) f Mƒ≈~L ogNormal(log(µ) = log(0.00012), 6 = 0.51).
Approximate posterior distributions are found using the Metropolis Hastings algorithm described in the appendix of Corbella et al. (2018).

B.4 Synthesis model
The synthesis model was proposed as a parsimonious attempt to utilise ILI, virological (combined as ILI consultations × swab positivity, to give a dataset labelled ILI + ), hospitalisation and ICU data to jointly infer and to predict the transmission dynamics of 2017/18 seasonal influenza outbreak. The schematic representation of the model is shown in the following figure: As with the SPC model, the reproductive number is linked to an initial rate of exponential growth, with rate -, Here, the mean infectious period ß 2 ; -® and the mean latent period ß 2 6 -® are both held fixed at 3.47 and 2.0 days respectively (Birrell et al 2011).

B.4.2 Disease and reporting processes
The disease and reporting component of the model is the component coloured blue in Figure S5. The following three parameters describe the proportion of infections that appear in the respective datasets: • f ghh : the proportion of infections that will lead to an influenza positive GP consultation for ILI, the product of the proportion of symptomatic infections (e) and the proportion that consult with a GP consultation given ILI (f gh ) (cf. Equation 3).
• f Ω : the proportion of infections that will lead to a hospital admission for ILI.
• f M|Ω : the proportion of hospitalisations for ILI that lead to an ICU admission.
The proportion infections leading to an ICU admission, f M is given by so that there is an implicit assumption that all ICU admissions are initially recorded as a hospitalisation.
In week Ü, the value of the ILI+ data ( In a slight variation from the ICU and SPC models, ∆ M (…) represents the weekly number of newly infectious (not infected) cases in week …. Correspondingly, k MLM gives the appropriate quantiles of a gamma delay distribution, with k MLM (Ü − …;•,•) describing the probability that the reporting of an ILI consultation occurs Ü − … weeks after an individual becomes infectious. The parameters governing this gamma distribution are those governing the distribution of the delay from symptom onset to GP consultations and are taken from Birrell et al (2011) to be i MLM = 3.06 and j MLM = 10.22 (when multiplied to a weekly, rather than daily, scale). In this study we simple assume that the incubation is equal to the latent period.
The weekly number of new hospital admissions at week w was similarly generated as Here k Ω is parameterised by i Ω = 0.708, j Ω = 1.813 from the assumed distribution of delay times from infectiousness onset to hospital admission.
The weekly number of new ICUs at week w was generated as Here k M|Ω is parameterised by i Mƒ≈ = 0.425, j Mƒ≈ = 2.163 from the assumed distribution of delay times from hospital admission to ICU admission. The parameters used for the above two delay distributions have been estimated on the basis of a simple analysis of USISS data from the 2012/13 and 2013/14 influenza seasons (not shown).

B.4.3 Inference
To reflect the over-dispersion in the weekly ILI+ numbers, the negative binomial likelihood is still assumed even though there is no guarantee that the ILI+ data, † ë MLM^, will be integer-valued. The observed number, ¢ ë MLM^, within a catchment population in week Ü of Z ë has probability density Here ¶ MLM is the dispersion parameter. Similarly, and more straightforwardly, the hospitalisation and ICU data, ¢ ë Ω and ¢ ë Mƒ≈ , also follow negative binomial distributions with dispersion parameters: ¶ Ω , ¶ Mƒ≈ respectively.
Assuming independence across the datasets and independent observation errors over time (conditional upon model parameters), the total likelihood given model parameter Φ is  Table S1). Again, posterior distributions are obtained via Markov Chain Monte Carlo simulations (MCMC). From these samples, we can obtain means and 95% credible intervals for model parameters.
In the followings we show the model predictions using the model calibrated at four different weeks: 2018w1, 2018w4, 2018w8 and 2018w12. These results were obtained from 5000 MCMC samples.    Table S1 gives a comparison of the assumed parameter values and, where the parameter is to be estimated, prior distributions for each of the four modelling approaches.    The main manuscript discusses the availability of pre-season (i.e. based on samples taken prior to week 40, 2017) serological data. These data are informative as they are indicative of the levels of pre-existing immunity in the population prior to the seasonal epidemic. In practice, this tells us the number of people in the susceptible states at time-0 (!(# $ ) in equations 1, 7, 10 and 12). The analyses of the pre-season serological data were complete by week 8 of 2018, towards the end of the modelling period under scrutiny in this manuscript. The modelling efforts included in the main manuscript were as completed in real-time, using the best available data at analysis time. The information that came from the analysis of the serological samples was therefore not quite timely enough to be properly understood within the required timeframe and therefore we study its impact here as a sensitivity analysis.

B.5 Prior Distributions and Model Assumptions
An initial assumption had been made that 0.9 was a suitable figure to describe the proportion of the population who were susceptible to infection with any of the circulating influenza strains. As can be seen from Table S1 this was the value for the model parameter π that has been used largely throughout. However, the analysis of the serological data revealed that 37.5% of samples returned titer values indicative of the presence of antibodies to all three circulating strains. This suggests that, for the single-strain models, a value of π=0.625 for the proportion of individuals susceptible to infection. This may still be an over-estimation of the proportion susceptible to infection, as the serology samples were taken in the summer of 2017, prior to the vaccination of ~14 million people in autumn 2018. If the vaccine is effective, this is likely to have a profound impact on population immunity if the vaccine is effective.
Different models handled initial susceptibility differently. Again, from Table S1 it can be seen that the SPC and Synthesis models use a fixed value, whereas the ICU models place a prior probability distribution on this value. This parameter cannot be directly identified from data, and therefore strong prior information is required, and this is reflected in the prior choice of the ICU model. For example, the ICU model uses a prior &~Beta(45,5), which has a mean of 0.9 and a small standard deviation of 0.04. In the sensitivity analysis, this is replaced by &~Beta(30,20) which has a mean of 0.6 and a standard deviation of 0.07, comfortably covering the 0.625 value used in the SPC and Synthesis model's sensitivity analysis.
The lowering of the initial susceptibility has the following impacts: 1. Leads to higher estimates of 2 $ (but not necessarily 2 3 ), and higher estimates of severity. In the models using GP ILI data, this manifests itself in a higher propensity to consult with a GP and in the models that use data on severe illness, it results in higher estimates of the ratio of ICU or hospital admissions per infection. 2. Induces lower estimates of attack rates, the cumulative number of infections caused by the epidemic (despite higher 2 $ ).
3. Affects the timing of the estimated peak in the number of underlying number of infections. Specifically, prior to any observed peak in the relevant data, lowering the initial susceptibility produces an estimated peak in incidence that occurs earlier. 4. Makes little difference to the estimated pattern of GP consultations/ICU admissions etc.
These are directly identified by the data and therefore are strongly anchored.
These findings are most clearly illustrated when looking at the Severity/ICU model. The SPC model, relying on direct data on GP consultations, has these sensitivities dampened by the presence of the background ILI consultations which can adapt to absorb some of the effects of changing the susceptibility parameter. Similarly, the multiplicity of data in the synthesis model and the competing strains in the SS model make these effects more difficult to clearly depict, though they are present. Figure S10 shows some parameter correlations obtained under the SS model using data up to the end of Week 1, 2018. Each plot shows the correlation of strain-specific parameters with the initial susceptibility. For susceptibilities >0.1, there is a clear negative correlation between 2 $ and π while there is no such obvious relationship for the effective reproductive number, 2 3 . Moving to an analysis that covers the full epidemic period, up to Week 12, based on the Severity/ICU model, Figure S11 shows a similar relationship. Here, this plot shows the posterior distributions of (top row) the susceptibility parameter π and the data dispersion parameter, 4 (middle row) the severity parameter, 5 ICU , and the school holiday effect on transmission, 9 (bottom row) reproductive numbers 2 $ and 2 3 , for each of the two priors placed on the susceptibility. The top left panel in each half shows that little is learned about the susceptibility parameter, whilst there are significant shifts in the parameters 5 ICU and 2 $ between the two halves.

Susceptibility Transmissibility
Figure S11 (On the left) Posterior distributions for parameters of interest obtained using a prior centred on π=0.9; (On the right) The same distributions obtained using a prior centered on π=0.6.
It stands to reason that, if the number of healthcare events (ICU admissions, GP consultations etc.) are anchored by data, yet there is a higher severity ratio, then it must mean there are consequently fewer infections. This also makes sense on the basis of similar values for the initial effective reproduction number 2 3 , but with a smaller susceptible population. The difference in the numbers infected are illustrated in the right-hand column of Figure S12 which shows the estimated cumulative numbers of infections over time. In the left-hand column of Figure S12 the plots show the estimated/forecast incidence per day. The error bars at the bottom of the plot indicate the estimated range for the timing of the peak in infection incidence. In 2018w1 (top-left panel) the estimated peak in incidence is slightly earlier when there is a lower susceptibility. The same phenomenon occurs for the analysis at 2018w2, the timing of the peak, after which the peak estimates largely coincide. In all models, prior to the peak in incidence, there is a tendency to forecast an earlier peak in incidence at the lower level of susceptibility. However, once the peak has occurred and is being retrospectively estimated, the estimated timing is independent of the initial susceptibility. Figure S12 Left-hand column: plots of the estimated number of daily infections from the SPC model at two levels of susceptibility. The grey-shaded area corresponds to the incidence of infections forecast beyond the time of analysis (2018w1, top; 2018w12, bottom). Right-hand column: plots for cumulative incidence (and thus forecast attack rates for the influenza season).
The SPC model can provide age-and region-stratified results. The sample sizes in the serological data were too small to be able to adequately estimate age-specific susceptibilities. Whilst this didn't prove a significant barrier to regional modelling, it prevented adequate age-specific inference. If the same susceptibility is assumed across age-groups, typical contact matrices that describe the pattern of infection across the age groups such as POLYMOD (Mossong et al, 2008) drive infection towards school-age children, a phenomenon not borne out in the data. To adequately account for this, the serological data could usefully be augmented by information on vaccine uptake amongst children to build a better picture of the population's immunity profile prior to any influenza season.
This study therefore identifies the importance of having pre-season serological data analysed in a timely fashion in readiness for any coming influenza season. It enables more accurate estimation of model parameters, can help to determine infection rates underlying the burden being placed on healthcare services and, if serological samples are representative of different geographical regions and across age groups, it can facilitate more detailed modelling at a finer granularity across the population.