Estimating the impact of influenza on the epidemiological dynamics of SARS-CoV-2

Stringency index data

Country-level time series of the stringency index were available from the Oxford COVID-19 Government Response Tracker, developed at the University of Oxford and described elsewhere (Hale et al., 2020). Briefly, the stringency index provides an aggregate measure of the number and of the strictness of non-pharmaceutical control measures implemented by governments in response to the COVID-19 epidemic. The stringency index is defined as the average of nine normalized ordinal variables, which quantify the strength (e.g., recommended or required) and the scope (e.g., targeted or general) of closure and containment measures (8 variables) and of health measures (1 variable). The resulting index allows quantifying the strength of control measures in a systematic way, on a scale ranging from 0 (no interventions) to 100 (maximum number and maximal intensity of control measures). Of note, however, the stringency index does not quantify the impact of control measures, which likely varied across countries (Flaxman et al., 2020). In formulating our model, we therefore modeled the relationship between the stringency index and the relative reduction in SARS-CoV-2 transmission using a non-decreasing function, whose parameters represented the impact of control measures and were estimated from the data.

Influenza incidence data

Virological data on the weekly numbers of samples tested and of samples positive to any influenza virus were available from the FluNet database, compiled by the WHO (Fig. S1A). Parallel syndromic data on the weekly incidence rate of influenza-like illnesses (ILI) were available from the FluID database, also compiled by the WHO (Fig. S1B). These data were deemed high-quality and used in a previous study on influenza forecasting in the countries considered here (Kramer & Shaman, 2019). The weekly incidence rate of influenza was then calculated as the product of ILI incidence and of the fraction of samples positive to any influenza virus (Fig. 1B). Because the magnitude of influenza incidence thus calculated varied markedly across countries (e.g., as a result of different surveillance systems and case definitions), we rescaled each time series by its average during the period of co-circulation of influenza and SARS-CoV-2 (Fig. 1B). The resulting time series was therefore dimensionless and equalled 1 when influenza incidence equalled its average value during that period.

COVID-19 mortality data

Data on the daily number of deaths caused by SARS-CoV-2 (counted by date of death) were available from national public health public institutes, in Belgium (Belgian institute for health, sciensano)(2020), in Italy (Dipartimento della Protezione Civile, 2020), and in Spain (Instituto de Salud Carlos III, official data with historical corrections compiled by the media (DATADISTA, 2020). In Norway, the data were available from the worldwide database compiled by the European Center for Disease Control and Prevention (European Centre for Disease Prevention and Control, 2020). Following a previous study (Flaxman et al., 2020), and to avoid a possible bias caused by the dominance of deaths due to non-locally acquired infections early in the epidemic, we included observed deaths from the date after which the cumulative observed death count exceeded 10. Data points before that date were treated as missing and were assigned a conditional log-likelihood of 0, such that they did not contribute to the overall log-likelihood. The data were not further pre-processed, except in Italy, where a negative death count was reported on 24 June 2020 and was treated as missing and also assigned a log-likelihood of 0.

Model formulation

We formulated a variant of the standard Susceptible–Exposed–Infected–Recovered transmission model (Keeling & Rohani, 2008), using the method of stages to allow for realistic distributions of the latent, infectious, and onset-to-death periods (Lloyd, 2001; Wearing, Rohani & Keeling, 2005). Specifically, we assumed that the latent and infectious periods were Erlang-distributed with shape parameter 2 and mean 1/σ = 4 days and 1/γ = 5 days, respectively (Li et al., 2020). The resulting generation time T_g (i.e., the time from infection of a primary case to transmission to a secondary case) had a mean of 6.5 days and a coefficient of variation of 0.58 (see Fig. S2 for the full distribution and the details of the calculation), to Svensson (2007) and Camacho et al. (2011) consistent with empirical observations and with the values fixed in a previous modeling study (Flaxman et al., 2020; Bi et al., 2020). To model the impact of the gradual implementation of non-pharmaceutical control measures (e.g., border closure, school closure, lockdown), we mapped the stringency index (denoted by si(t)) to the time-varying relative reduction in transmission of SARS-CoV-2 (denoted by r_β(t)). Specifically, we used the following simple linear scaling function, with saturation: ${\displaystyle r_{\beta }\left(t\right)=min\left(1,b\times \frac{si\left(t\right)}{100}\right)}$

Here the parameter b quantifies how steeply the transmission rate of SARS-CoV-2 decreases as the stringency index increases. Hence, this parameter can be interpreted as a measure of the impact of non-pharmaceutical control measures on SARS-CoV-2 transmission. The deterministic variant of the model was represented by the following set of differential equations: ${\displaystyle \begin{array}{ccc}\hfill \dot{S}& \hfill =& \hfill -\lambda \left(t\right)S\\ \hfill {\dot{E}}_1& \hfill =& \hfill \lambda \left(t\right)S-2\sigma E_1\\ \hfill {\dot{E}}_2& \hfill =& \hfill 2\sigma \left(E_1-E_2\right)\\ \hfill {\dot{I}}_1& \hfill =& \hfill 2\sigma E_2-2\gamma I_1\\ \hfill {\dot{I}}_2& \hfill =& \hfill 2\gamma \left(I_1-I_2\right)\\ \hfill \dot{R}& \hfill =& \hfill 2\gamma I_2\end{array}}$

The force of infection (that is, the per capita rate at which susceptible individuals contract infection (Keeling & Rohani, 2008)), λ(t), was modeled as: ${\displaystyle \begin{array}{ccc}\hfill \lambda \left(t\right)& \hfill =& \hfill \beta \left(t\right)\frac{I_1+I_2}{N}\\ \hfill \beta \left(t\right)& \hfill =& \hfill {\beta }_0\left(1-r_{\beta }\left(t\right)\right){\beta }_F\left(t\right)\\ \hfill r_{\beta }\left(t\right)& \hfill =& \hfill min\left(1,b\times \frac{si\left(t\right)}{100}\right)\\ \hfill {\beta }_F\left(t\right)& \hfill =& \hfill max\left(0,1+{\beta }_{F}F\left(t\right)\right)\\ \hfill R_e\left(t\right)& \hfill =& \hfill \frac{\beta \left(t\right)}{\gamma }\times \frac{S\left(t\right)}{N}\end{array}}$

where R₀ represents the basic reproduction number of SARS-CoV-2, β₀ = R₀γ the basic transmission rate, R_e(t) the time-varying effective reproduction number, N the population size (assumed constant during the study period), and F(t) the renormalized time series of influenza incidence, incorporated as a covariate into the model (Fig. 1). With this formulation, the parameter β_F quantifies the impact of influenza on SARS-CoV-2 transmission: β_F > 0 if influenza increases transmission, β_F < 0 if influenza decreases transmission, and β_F = 0 if influenza has no impact on transmission (null hypothesis). More specifically, the average incidence of influenza during the period of co-circulation with SARS-CoV-2 corresponds to F(t) = 1, such that 1 + β_F represents the average relative variation of SARS-CoV-2 transmission associated with influenza. In writing the equations, we implicitly assume that the impact of influenza on SARS-CoV-2 transmission, if any, is short-lived and does not extend long after influenza infection.

Finally, we incorporated an observation model that related the dynamics of SARS-CoV-2 infection to that of COVID-19 mortality, taking into account the fact that only a fraction of infections results in death and that, among those, death occurs some time after symptom onset (Flaxman et al., 2020; Verity et al., 2020; Khalili et al., 2020). We assumed an average duration of pre-symptomatic of 2.5 days, resulting in an average incubation period of 6.5 days, in broad agreement with previous empirical studies (Khalili et al., 2020; Tindale et al., 2020). Hence, individuals in the first infected state (I₁) were considered pre-symptomatic, and the onset of symptoms was assumed to coincide with the transition from I₁ to I₂. The onset-to-death time was then assumed to be Erlang distributed with shape parameter 5 and mean 1/κ = 17.8 days (coefficient of variation of 0.45), the value estimated in a previous epidemiological study (Verity et al., 2020). In a sensitivity analysis, we also tested a mean onset-to-death time of 1/κ = 13 days, the lower bound estimated in a meta-analysis (Khalili et al., 2020). According to previous studies in European countries, the infection fatality ratio (IFR) typically ranged from 0.5% to 1% early in the epidemic (Salje et al., 2020; O’Driscoll et al., 2021; Pastor-Barriuso et al., 2020). We fixed the IFR to μ = 0.01 in the base model, but we considered an alternative value of 0.005 in a sensitivity analysis. Given those assumptions, the observation model was modeled by the following set of ordinary differential equations: ${\displaystyle \begin{array}{ccc}\hfill {\dot{Q}}_1& \hfill =& \hfill 2\gamma \mu I_1-5\kappa Q_1\\ \hfill {\dot{Q}}_{i=2,\dots ,5}& \hfill =& \hfill 5\kappa \left(Q_{i-1}-Q_i\right)\\ \hfill {\dot{D}}_M& \hfill =& \hfill 5\kappa Q_5\end{array}}$

Here D_M is the simulated number of daily deaths, modeled as an accumulator variable and reset to 0 at the end of each day. The observed number of daily deaths, D_O, was modeled using a negative binomial distribution with mean D_M and over-dispersion k_D (i.e.,  $\mathbb{V}\left(D_O|D_M\right)=D_M+k_D{D}_M^2$), a standard distribution used in previous modeling studies (Flaxman et al., 2020; King et al., 2015).

As in Flaxman et al. (2020), simulations were started 30 days before the date from which the cumulative observed death count first exceeded 10. At that date, we assumed that E₁(0) individuals had been exposed to SARS-CoV-2; other individuals were assumed susceptible to infection (i.e.,  S(0) = N − E₁(0)), and all other compartments were initialized to 0.

Stochastic variant and modeling of superspreading

The stochastic variant of the model was implemented as a continuous-time Markov process approximated via a multinomial modification of the τ-leap algorithm (He, Ionides & King, 2010), with a fixed time step of Δt = 10⁻¹ day. To model the effect of superspreading events—a key feature of SARS-CoV-2 transmission dynamics (Althouse et al., 2020)—, extra-demographic stochasticity was added to the transmission rate β(t). Specifically, as proposed by Kain et al. (2020), at every time step we drew a value of β₀ from a Gamma white noise distribution: ${\displaystyle {\beta }_0\sim {\Gamma }_{\mathrm{WN}}\left(\sigma =\sqrt{\frac{R_0}{\left(I_1+I_2\right)k\Delta t}},\mu =R_0\gamma \right)}$

with mean µand variance μσ². Here k represents the dispersion parameter of the Negative-binomial distribution for the individual reproduction number (with mean R₀ and variance $R_0+\frac{R_0^2}{k}$), as estimated in previous studies (Lloyd-Smith et al., 2005; Endo et al., 2020). As in Kain et al. (2020) and in keeping with empirical estimates from contact tracing studies of SARS-CoV-2 (Endo et al., 2020), we fixed k = 0.16.

Model estimation

Following the method presented in Flaxman et al. (2020), we estimated unknown model parameters using observed COVID-19 mortality data alone. Indeed, because of the initially limited testing capacity (typically reserved to severe cases or high-risk groups), mortality data were arguably more reliable than case data early in the epidemic in most countries (Flaxman et al., 2020). By incorporating known epidemiological parameters (key among those the onset-to-death time, the IFR, and the generation time), however, the method allows back-calculating infection rates from observed death rates. Hence, in addition to the dynamic of mortality, we also reconstructed the dynamic of infection and, as a validation, compared it to external epidemiological data—like cross-sectional seroprevalence estimates—when available.

The following five parameters were estimated from the data:

1. The basic reproduction number, R₀. According to a previous meta-analysis (Alimohamadi, Taghdir & Sepandi, 2020), this parameter was searched in the interval 1–10.

2. The impact of non-pharmaceutical control measures, b. The lower bound of the search interval of this parameter was fixed to 0.5, such that the maximal value of the stringency index (s = 100) corresponded to a minimal reduction of SARS-CoV-2 transmission of 50% (Flaxman et al., 2020).

3. The impact of influenza on SARS-CoV-2 transmission, β_F (search interval: ℝ).

4. The initial number of individuals exposed to SARS-CoV-2, E₁(0) (search interval: 0–10⁴).

5. The over-dispersion in death reporting, k_D (search interval: ℝ⁺).

A summary list of fixed and estimated model parameters is presented in Table 1.

All parameters were transformed to be estimated on the real line, using a log transformation for positive parameters and the extended logistic function $f\left(\theta \right)=log\frac{\theta -a}{b-\theta }$ for parameters constrained in the interval [a, b]. The maximum iterated filtering algorithm (MIF2, Ionides et al., 2015), implemented in the R (version 3.6.3) package pomp (King, Nguyen & Ionides, 2016; R Core Team, 2020) (version 2.7), was used to estimate model parameters. The R checkpoint package was used to freeze all the packages’ version at the date of April 3, 2020 (De Vries, 2020). The estimation was completed in several steps, starting with trajectory matching to identify good starting parameters for MIF2, followed by 100 independent runs of MIF2 to locate the maximum likelihood estimate (MLE). Each MIF run had 150 iterations with 5,000 particles, geometric cooling, and a random walk standard deviation of 0.1 for the initial condition E₁(0) and of 0.02 for the other parameters. The log-likelihood of every parameter set was calculated as the log of the mean likelihood of 5 replicate particle filters, each with 20,000 particles. The profile likelihood was calculated to verify the convergence of MIF2 and to derive an approximate 95% confidence interval for the parameter β_F (Raue et al., 2009)—the parameter of key interest in our study. For the other parameters, a parametric bootstrap was used to calculate approximate 95% confidence intervals, by re-estimating the parameters for each of 200 synthetic datasets simulated at the MLE (Domenech de Celles et al., 2018; Domenech de Celles et al., 2019). Compared with the profile likelihood, the parametric bootstrap requires less computation and was found to perform well in previous applications (Domenech deCelles et al., 2018; Domenech de Celles et al., 2019).

Sensitivity analyses

To verify the robustness of the parameter estimates, we conducted six sensitivity analyses (further detailed in the Supplementary Results). First, we estimated an extended model in which the reduction of SARS-CoV-2 transmission was allowed to scale non-linearly with the stringency index. Second, to test the possible presence of other variables confounded with influenza, we estimated a model that included an exponential trend in the transmission rate of SARS-CoV-2. Finally, we varied the fixed value of 3 parameters and re-estimated the parameters of the base model. Specifically, we tested two alternative values of the average generation time (𝔼(T_g) = 5 days and 𝔼(T_g) = 7.5 days), one alternative value of the infection fatality ratio (μ = 0.005), and one alternative value of the average onset-to-death period ($\frac{1}{\kappa }=13$ days).