Competition between RSV and influenza: Limits of modelling inference from surveillance data

Respiratory Syncytial Virus (RSV) and Influenza cause a large burden of disease. Evidence of their interaction via temporary cross-protection implies that prevention of one could inadvertently lead to an increase in the burden of the other. However, evidence for the public health impact of such interaction is sparse and largely derives from ecological analyses of peak shifts in surveillance data. To test the robustness of estimates of interaction parameters between RSV and Influenza from surveillance data we conducted a simulation and back-inference study. We developed a two-pathogen interaction model, parameterised to simulate RSV and Influenza epidemiology in the UK. Using the infection model in combination with a surveillance-like stochastic observation process we generated a range of possible RSV and Influenza trajectories and then used Markov Chain Monte Carlo (MCMC) methods to back-infer parameters including those describing competition. We find that in most scenarios both the strength and duration of RSV and Influenza interaction could be estimated from the simulated surveillance data reasonably well. However, the robustness of inference declined towards the extremes of the plausible parameter ranges, with misleading results. It was for instance not possible to tell the difference between low/moderate interaction and no interaction. In conclusion, our results illustrate that in a plausible parameter range, the strength of RSV and Influenza interaction can be estimated from a single season of high-quality surveillance data but also highlights the importance to test parameter identifiability a priori in such situations.


Introduction
Respiratory Syncytial Viruses (RSV) and seasonal influenza viruses cause large burdens of respiratory disease, including in young children (Lafond et al., 2016;Shi et al., 2017).RSV was recently identified as the primary cause of hospitalisation for severe paediatric pneumonia (Pneumonia Etiology Research for Child Health (PERCH) Study Group, K. L. et al., 2019), particularly in the neonatal period.In the northern hemisphere both viruses cause pronounced annual winter epidemics peaking between October and March (Bloom-Feshbach et al., 2013).
Evidence from epidemiological and biological studies implies there is competitive interaction between influenza and RSV (Opatowski et al., 2017;Velasco-Hernández et al., 2015;Mak et al., 2012;Pascalis et al., 2012;Yang et al., 2012;Walzl et al., 2000).The biological mechanism for competition is activation of the innate "antiviral response" by infection that can inhibit further or subsequent infection (Walzl et al., 2000;Ascough et al., 2018;Lee et al., 2018), resulting in a period of cross-protection during and after infection.Mouse studies have shown this effect, where following influenza infection or live attenuated influenza vaccination (LAIV), RSV replication/severity was decreased (Walzl et al., 2000;Lee et al., 2018).Within-host animal studies, both in vitro and modelling, have shown that the growth rates of the viruses can be affected by other viruses present (Pinky and Dobrovolny, 2016;Shinjoh et al., 2000).The duration of this cross-reactive response is debated, varying from "short-term" (Ferguson et al., 2005), less than two weeks (Laurie et al., 2015) or up to 3 months (Kelly et al., 2010).Influenza epidemics caused by different strains are thought to exhibit competitive exclusion (Opatowski et al., 2017;Ferguson et al., 2003;Kucharski et al., 2016), and for RSV and influenza syndromic surveillance has shown shifts in the seasonal incidence peaks of RSV following abnormal (pandemic or early) influenza seasons (Mak et al., 2012;Hirsh et al., 2014;Meningher et al., 2014;Gröndahl et al., 2014;Casalegno et al., 2010;van Asten et al., 2016), which suggest this mechanism may not only occur but can substantially alter the epidemiology of influenza and RSV.There is, however, little evidence that links the strength of competition between RSV and influenza within a host to observed population dynamics.Understanding the dynamics is critical for predicting the effects of alteration of their ecological balance, for example through vaccination programs, and is the motivation for this study.
Influenza vaccination rates, especially in key transmission groups could disrupt transmission, potentially leading to effects on interacting viruses.In the UK, the RSV epidemic usually precedes the influenza epidemic by one or two months, so reduced influenza transmission as a result of childhood influenza vaccination may not affect RSV transmission dynamics.However, the competitive pressure exhibited by RSV on influenza may become highly relevant soon.The only RSV vaccine candidate yet that completed Phase 3 trials, the maternal vaccine, Novavax, demonstrated only partial efficacy that the Advisory Committee for Immunization Practices in the US deemed insufficient to warrant licensure (Novavax Announces, 2021).However, the RSV vaccine pipeline contains a number of Phase 1 and 2 candidates that aim to protect children in part by limiting RSV circulation.As such, these future RSV vaccines have the potential to decrease the competitive pressure on influenza and thereby increase influenza as an unintended consequence, both in children and other age groups, as children are a key driver of transmission.These impacts will need to be considered as part of their cost benefit proposition preceding routine use.
Mathematical modelling is an important tool for testing mechanisms and hypotheses of epidemiologically significant RSV and influenza competition, such as the hypothesis that they competitively interact.Models offer an opportunity to mechanistically combine observations from surveillance data and extrapolate beyond the observed.However, in the case of RSV and influenza competition the identifiability of model parameters from viral surveillance data is uncertain.Hence, we conducted a simulation study to test whether parameters can be backinferred from a range of realistic model-generated scenarios that include only partial observation of the infection dynamics from surveillance-like data.

Model structure
We developed an age-stratified deterministic compartmental transmission model for RSV and Influenza with interactions (Fig. 1 and Supplement Section 2).The population could be Susceptible (S), Infectious (I), Protected (P) or Recovered (R) for each of RSV and influenza viruses.We simulated one season so we did not consider potential loss of immunity, and current estimates for RSV immunity lasts less than a year (Weber et al., 2001), and we take influenza immunity into account by fitting the percentage susceptible at the start of the season (see below).There were separate transmission and recovery rates for each virus (subscripts RSV and INF), and i and j denote age groups.Susceptible individuals were infected at rates λ INF, i and λ RSV, i and enter the I compartment.They recovered at rates γ INF and γ RSV , and entered the P compartment where they were no longer infectious.In P, individuals were fully protected against homologous re-infection and also had some cross protection against the second virus.Loss of cross-protection occurred at rate ρ.Infection with a second virus was less likely in the I and P classes and occurred at a rate reduced by (σ).The key parameters determining interaction are therefore the strength of competition (σ) and the rate of loss of cross-protection (ρ).Compartments I RSV , i P INF,i and I RSV,i R INF,I , as well as P RSV,i I INF,I and R RSV,i I INF,I were combined as they were effectively identical when modelling only one season.
The model was stratified into 5 age categories: infants: 0− 1 years, pre-school-aged children: 2− 4 years, school-aged children: 5− 15 years, adults: 16− 64 years, and older adults: aged 65 + .Age-dependent contact patterns relevant to the transmission of infections are highly age heterogenous (Mossong et al., 2008), and we used social contact patterns (including both physical and verbal contacts) in England from POLYMOD (Mossong et al., 2008), a European wide contact study in 2005/6, and in the socialmixr R package (Funk, 2018).We calculated forces of Infection, λ RSV,i and λ INF,i , from the baseline transmission rates β INF and β RSV and the mixing parameters as: where α ij is the contact rate between groups i and j and I INF,j and I RSV,j are the proportion infected with Influenza and RSV in age group j.See Supplement Section 1 for model equations.
We modelled one year from the start of the respiratory virus season, and initiated the model with a proportion of each age group susceptible to influenza set from serological data (Baguelin et al., 2013) (Table 2) and the rest in S RSV R INF .RSV immunity to re-infection may last less than a year (Weber et al., 2001), therefore we considered the population to be N.R. Waterlow et al. susceptible to RSV at the start of the season.However, RSV susceptibility differs with age (Henderson et al., 1979), and therefore we reduced the susceptibility to the same range as in other models (Moore et al., 2014) by decreasing the infection rate by a susceptibility parameter, τ i (Table 2).
RSV was seeded at time η RSV when one individual from the fully Susceptible class (S RSV S INF ) becomes infected (I RSV S INF ).Influenza infections are introduced at a rate of 0.1 cases per day from S RSV S INF to S RSV I INF , starting on day η INF .Influenza introduction assumptions differ from those of RSV as with a single introduction the influenza epidemic was supressed for the whole season at certain parameter values, which is not seen in UK surveillance.See Supplementary section 6 for further details.
An observation process layer converted infections to detected cases using a binomial distribution.The number of detected cases is assumed to follow a binomial distribution as follows: where X is the number of detected cases in n infections and x virus,i,t is the cases detected for each virus, age group, and timestep.The proportion detected was different for each age group and virus (Table 2).We implemented the model in R (R Core Team, 2018) and C++ using the Rcpp (Eddelbuettel and Francois, 2011) and deSolve (Soetaert et al., 2010) packages.

Simulated data
We generated simulations to resemble data collected through surveillance systems in the UK (Reeves et al., 2017), such as the Respiratory Datamart System in England and Wales (Zhao et al., 2014) and other Public Health England surveillance systems (Weekly national flu reports, 2021; Respiratory infections, 2021).This provides laboratory test results from routinely tested clinical respiratory samples from a range of respiratory viruses.The proportion detected varies by age-group for RSV (Table 1), as younger infants are more likely to present with severe symptoms (Ohuma et al., 2012).Output from the model is weekly number of positive tests in the under-five population for RSV and influenza.
We generated simulations with parameter values from the literature and if unavailable we calibrated the values to realistic ranges (Table 1).
Across simulations we varied σ (strength of interaction), for which we used 11 different values, and ρ (the rate of loss of cross-protection), for which we used 5 different values.This resulted in 55 combinations of σ and ρ and we simulated 5 replicates of each.

Parameter estimation
We assumed that the observed cases followed a Poisson distribution with likelihood: where θ is the modelled detected cases of RSV and influenza in the two youngest age groups, x j is the observation and n is the total number of observations.We fitted only to the lowest 2 age groups to represent where the majority of samples for RSV are taken from and detected.We fitted the model to simulated data using Metropolis Hastings Markov Chain Monte Carlo (MCMC) sampling.Estimated parameters were transmission rates (β interaction parameters (ρ, σ) and season start times (η RSV , η INF ).For each scenario, we ran two chains with 450 000 iterations as burn in followed by a further 250 000 iterations.For chains that did not converge, we extended the chains for a further 250 000 iterations iteratively until convergence was reached or a total of 1 200 000 iterations were run.We used weak priors, and the priors for β RSV and β INF were calculated from R 0 values, assuming no interaction (see Supplement Section 3).We adapted the shape of the proposal distribution during burn in, from 5000 accepted proposals to a further 300000 proposals, to take correlation between parameters into account by allowing the covariance matrix for proposal distributions to change.Parameter limits are defined in Supplement Section 7 and β RSV , β INF , Δ RSV1 , Δ RSV2 , Δ INF and ρ were sampled on a log scale to improve mixing where the parameter values were very low.
We assessed MCMC convergence via the Gelman-Rubin statistic, which compares the within-chain variance to the between-chain variance for each parameter.Scenarios with a statistic greater than 1.1 we deemed as practically unidentifiable from simulated data.We also calculated the Pearson correlation coefficient between each two estimated parameters and assessed how these changed with the values of the interaction parameters (σ and ρ), in order to further understand difficulties with parameter estimation.
We compared the inferred parameter estimates to the simulated parameter values to determine inaccuracy and imprecision of the fit, where inaccuracy is defined as the difference between the median value of the posterior distribution and the true value, and imprecision is the range between the 95 % credible intervals (95 % CI).We present results from one replicate set of simulations in Results and others are given in the supplement.

Epidemic profiles
Altering the strength or duration of cross-protection did not notably affect the timing or shape of the RSV epidemic (Fig. 2), due to the higher transmission rate and earlier start of RSV in our scenarios.However, increasing the strength or duration of interaction delayed the influenza peak.The total number of influenza infections in the youngest two age groups did not change (percentage difference <1% between σ = 1 and σ = 0) with the strength of interaction (Supplement Section 8).Increasing the duration of cross-protection resulted in an 11 % lower total number of infections from the shortest (2 days) to longest (40 days) duration of cross-protection (Supplement Section 8).Plots showing the epidemic curves for each infectious compartment are shown in in Supplement Section 9.

Correlation analysis
The most strongly correlated parameters were consistently the transmission rate for RSV (β RSV ) with the RSV season start time (η RSV ) and the transmission rate for influenza (β INF ) with the detection rate for influenza (Δ INF ) and the start of the influenza season (η INF ) (Fig. 3A).The correlation between parameters changed dependent on the values of the interaction parameters, an example of which is shown in Fig. 3B, where the correlation coefficient between the strength of interaction (σ) and the proportion of influenza cases detected (Δ INF ) varies depending on the values of σ and ρ.As the strength of interaction decreases (as σ→0), the correlation between the strength of interaction and the proportion of influenza cases detected becomes more positive.The correlation changes across the interaction parameters for other parameter combinations are shown in in the Supplement Section 10, and matrices for individual simulations are in supplement section 12.

Inferring the strength of cross-protection (σ)
Across simulations, the imprecision and inaccuracy of the estimated strength of cross-protection (σ) varied (Fig. 4), with the imprecision ranging from 0.15 to 0.66 (where 1 is poor precision) and average imprecision decreasing as the duration of protection (ρ) increased.We did not observe a trend in the inaccuracy of the parameter estimates and they ranged from 0 to 0.24.However, the lowest value tested (σ = 0.01)  was overestimated in each simulation and the highest value tested (σ = 0.99) was underestimated, showing that the extreme values are less well estimated.This may be due to the true value being very close to the limit.

Inferring the duration of cross-protection (1/ρ)
The imprecision of the estimated duration of cross protection ranged from 7 to 87 days (Fig. 5).Estimates were generally less precise when the period of cross-protection is longer (Fig. 5).In 70 % of our simulations   where σ interaction reduced the transmission rate by no more than 10 % (i.e.σ = 0.1 or 0.01) the duration of protection estimates exceeded an imprecision of 50 days.For scenarios assuming stronger competition estimates were much more precise.Indeed, one would expect that once the strength of competition is negligibly small the duration of such protection would be largely irrelevant.ρ estimates increased for simulations generated with longer duration of protection (smaller ρ).

Variation between replicates
Each parameter set was used to generate a further four replicate simulations of the observation process.Of these 275 simulations (5 replicates of 55 parameter combinations), 259 reached convergence at the cut-off point (see Supplement Section 11).The true parameter values were included in the 95 % CI in the majority of replicates (Fig. 6).The true value of ρ was not included in the 95 % CI in 6 simulations (2%), whereas the true value of σ was not included in the 95 % CI in 31 simulations (11 %).These simulations were more concentrated in areas with extreme interaction strengths (0.99 and 0.01) and very short duration of protection.We conclude from this that the stochastic variation in the simulation of the observation can occasionally result in difficulty estimating the true value of the parameter.

Discussion
We tested whether a transmission model including competitive interaction between RSV and influenza is identifiable from a single season of simulated high-quality surveillance data.We determined that it is possible to re-estimate strength and duration of interaction in most tested scenarios, although often imprecise due to large credible intervals, but that there are some areas of parameter space where posterior estimates are potentially misleading, particularly when the strength of interaction is assumed to be low or the duration of interaction short.However, we only estimated the parameters from information from one season at a time, and without strong priors.
While we are the first to test robustness of RSV and influenza competition inference, other identifiability studies, e.g. on Rift Valley Fever, have previously highlighted the importance of robustness testing to avoid misleading conclusions stemming largely from insufficient power of the data to inform the model parameters of interest (Kao and Eisenberg, 2018;Tuncer et al., 2016).Structural and practical identifiability analysis have also been used to select appropriate models, given the data available (Tuncer et al., 2018;Roosa and Chowell, 2019), for example a study that evaluated six different Zika models, and the identifiability of parameters within each (Tuncer et al., 2018).
Our analysis shows that there are potentially misleading results at extreme competition values, and it is almost impossible to get a "null estimate" for the strength of competition from this study.Evidence from mouse models suggest that the duration of RSV cross protection following influenza infection may last more than two weeks (Walzl et al., 2000;Hamilton, 2021) which, under the assumption that the duration of cross-protection is non-differential to the initiating virus, may suggest that the imprecision of our estimates at short durations of cross protection is unlikely to be a key risk for inference.However it may not be possible to distinguish such competition from no competition in our model, if the strength of the competition is low.We deliberately used uninformative priors for this parameter in order to be able to fully explore its identifiability, however, subsequent work may further improve precision of estimates by including prior estimates based on published evidence.This may also reduce the correlation of estimated parameters, which has challenged convergence in our simulations.Similarly, mouse models have suggested strong modulation of the RSV immune response if preceded by an influenza infection (Walzl et al., 2000), which may suggest that difficulties in our inference in scenarios that assume very small amounts of competition may not be the most relevant.
For parameter combinations where the simulated parameter value could not be re-estimated, we found that despite the relatively high assumed sample size stochastic noise from the observation model can occasionally result in incorrect estimates.This implies that inference based on a single season may be misleading purely because of the observational process associated with surveillance, however, including multiple seasons of observation should limit problems stemming from the observation process alone and further increase accuracy of estimates.
We assumed that the RSV-influenza interaction was bidirectional; particularly we assume that the strength and duration of interaction that influenza exhibits on RSV is the same as vice versa.Given that the proposed mechanisms for interaction are not virus specific this seems reasonable, and is supported by studies looking at the shift in RSV epidemics following the early 2009 influenza pandemic (Mak et al., 2012;Hirsh et al., 2014;Gröndahl et al., 2014;Casalegno et al., 2010).However, the RSV epidemic in the UK typically precedes influenza and similarly we only investigate such scenario.Therefore, in this work we can only estimate the competition of RSV on influenza dynamics and do not have power to estimate the other direction.Hence our results are applicable for considerations around RSV vaccine introduction but should be treated cautiously for any studies interested in the impact of Influenza on the transmission dynamics of RSV.This model did not include multiple strains of either RSV or influenza, which could have an impact on the interaction dynamics, as the interaction may differ between strains.Including strains would significantly increase the complexity of the model (see review on strain interaction models (Kucharski and Gog, 2012)), which we think would have rendered it unidentifiable.In addition, the aim was to assess the practical identifiability of the model parameters that govern viral interaction from routine surveillance data, and in many scenarios the surveillance data does not record strain type.The biological mechanism underpinning the period of cross-immunity is that of viral infection-induced protection, which is potentially induced by many viruses so may not be specific to RSV and influenza, and may not differ between influenza subtypes/strains.
We ran the model for one season at a time in order to reduce the complexity, as in other influenza models (Baguelin et al., 2013).We captured influenza immunity from previous seasons in the proportion of individuals susceptible for influenza at the start of the season, and RSV immunity is considered to last less than a year (Weber et al., 2001), so we simplified to a single season but included the major multi-season effects.A further sensitivity analysis could be to vary the susceptibility of individuals to influenza at the start of the year, in order to simulate different dominant influenza strains.We have however not included this, as our aim here was to look at the identifiability of parameters, and these differences would be taken into account when fitting to surveillance data from different seasons.Ideally, we would fit to multiple seasons of surveillance data, in order to account for variations by year.In the model we assumed a constant, age-dependent observation rate, as in other influenza models (Magal and Webb, 2018).Time varying reporting rates would substantially hinger inference, in fact a previous study comparing model fit of age-dependent vs time and age-dependent reporting rates concluded it was not possible to prefer one model terms of fit alone (Dorigatti et al., 2012), so we assume age-dependent only reporting rates for simplicity.We did not include additional seasonal effects in the model.While no or small effects have been reported for RSV and seasonal factors (Tian et al., 2017;Hogan et al., 2016), there is stronger evidence for the impact of climatic factors on influenza transmission, particularly ambient temperature and absolute humidity (Shaman et al., 2010;Shaman and Kohn, 2009;Lowen and Steel, 2014).While this does not affect our results on identifying parameters from simulated data, it should be noted as a potential confounder when estimating these parameters using surveillance data.
Further data may help to identify parameters in the model where it currently has difficulties.Data on the frequency of co-infections would allow us to use stronger priors for the strength of interaction, as well as providing an informative data source to fit the model to.Further information on the circulation of RSV, as opposed to only clinical cases, is also important due to the current uncertainty in infection numbers.In addition, surveillance systems would ideally provide daily data on RSV and influenza cases, giving us more granularity and potentially allowing us to identify all areas of parameter space.
Behavioural changes may also impact respiratory viral circulation, after infection with a virus (staying inside while recovering), or largescale behavioural change due to restrictions (social distancing measures in response to the SARS-CoV-2 pandemic).Whilst we do not investigate these mechanisms in this paper, such changes can have drastic impacts, such as the largely absent 2020 influenza season in Australia (Sullivan et al., 2020).
Overall, this study shows that in principle interaction parameters can be estimated from high quality surveillance-like data using mathematical models, although the precision and accuracy of the estimates varies depending on the scenario and stochasticity in the surveillance data.More power to reliably infer parameters may be available if fitting multiple seasons.It also highlights the importance of validating complex models, especially in light of the rapid development of models in emergency situations, which can have large impacts on public policy (Panovska-Griffiths, 2020).

Fig. 1 .
Fig. 1.Model diagram for RSV and Influenza (INF).Individuals could be either Susceptible (S), Infected, (I), Protected (P) or Recovered (R) to either virus.Following infection, (which occurred at rate λ RSV,i and λ INF,i ), recovery occurred at a constant rate (γ RSV and γ INF ), and the population entered the P state.Here they are immune to the virus they were infected by and protected to a varying extent (σ) against infection from the second virus.This protection waned at rate ρ, and the population entered the R compartment.In the R compartment the population was immune to the virus it was infected by, but not the other virus.We ran the model for one season and compartments I RSV, i P INF,i and I RSV,i R INF,I, were combined, and P RSV, i I INF,I and R RSV,i I INF,I were combined, because they are effectively identical.Parameters were: age susceptibility to RSV infection (τ i ), For clarity, age structure is given only by the subscript (i), for further details see supplement section 2.

Fig. 2 .
Fig. 2. Mean weekly incidence of observed cases in under 5 s (sum of age groups 0-1 and 2-4) from simulations with A) varying σ values and a fixed protection duration of 10 days (ρ = 0.1), and B) varying ρ values, and a fixed σ of 0.5.Simulations were run and sampled 1000 times for each parameter set and the shaded windows are the 95 % quantiles for each week.In both A and B the top panel shows the observed cases for RSV, and the lower panel the cases for Influenza.

Fig. 3 .
Fig. 3. A) Mean Pearson correlation coefficient between parameters.B) Correlation coefficient between σ ( strength of cross-protection) and Δ INF (start day of influenza).This is shown for 1 simulation, but the patterns were similar for all (Supplementary Section 12).

Fig. 4 .
Fig. 4. A) Estimated σ values for simulations with different σ and ρ values.Median value and 95 % CI are shown.The black line is the simulated (true) value of σ in each case.B) Imprecision of σ estimates calculated as the 95 % quantile range.C) Inaccuracy of the σ estimates, calculated as the difference between the posterior median and the true value.
N.R.Waterlow et al.

Fig. 5 .Fig. 6 .
Fig. 5. A) Estimated 1/ ρ values for simulations with different σ and ρ.Lines represent 95 % quantiles of the posterior sample and the circle represents the median value.The black line shows the true 1/ ρ value in each case.B) Imprecision of ρ estimates calculated as the 95 % quantile range.C) Inaccuracy of the ρ estimates, calculated as the difference between the posterior median and the true value.

Table 1
Parameter values used for generating simulations.

Table 2
Demography and susceptibility input used for model simulations.