Data sources
The original data collection was approved by the Swedish Ethical Review Authority. The current study does not involve research with human subjects, and written informed consent was waived by the Yale Institutional Review Board.
Data on the pediatric population in Northern Stockholm were obtained from the online base Statistics Sweden (https://www.statistikdatabasen.scb.se). The catchment area of the hospital includes the northern districts (Swedish: kommuner) of the Stockholm region and approximately 51% (51.8% in 1998 and 50.7% in 2018) of the population residing in central Stockholm. The number of children (0-17 years) living in the catchment area increased from 201,334 in 1998 and to 265,268 in 2018. Of these, 5.6% were under 1 year old, 23.1% were aged 1–4 years, 28.9% were aged 5–9 years, and 42.4% were aged 10–17 years.
Astrid Lindgrens Children’s Hospital at Karolinska Hospital Solna, is a tertiary referral university hospital with a pediatric intensive care unit (with 6–10 beds during this time frame) and an extracorporeal membrane oxygenation center (with 2–4 beds). The neonatal intensive care unit is a regional referral unit for babies from gestational age 22+ weeks. The emergency department had approximately 48,000-50,000 visits/year in the later time periods. The number of hospital beds varied between 150-190 over the study period.
The number of live births in the population residing in the hospital’s catchment area was approximately 10,750 in 1998 and 14,832 in 2018. The birth clinic at Karolinska University Hospital in Solna has approximately 2500 live births per year and prioritizes high-risk births. There are 3 more birth clinics at other hospitals in Stockholm, and another, smaller birth clinic opened and closed during the study period.
Stockholm experienced the first wave of the 2009 Influenza A (H1N1) pandemic at the end of summer 2009, with a total of 15 pediatric admissions at Astrid Lindgren Children’s Hospital from mid-July until mid-September. After a three-week hiatus, pediatric admissions resumed in early October, peaking in mid-November before declining sharply, with the last admission occurring in early December. In total, 77 children were hospitalized during the second wave of the 2009 Influenza A (H1N1) outbreak. The impact of the pandemic on hospitalized children in Stockholm has also been described in relation to co-circulating viral infections and within the context of a 16-year span of influenza epidemics52,53.
Clinical data
Positive viral findings of RSV and influenza virus among pediatric patients (under 18 years old) were provided by the Karolinska University Hospital Virology unit and matched to medical records. Virus detection was performed on samples from nasopharyngeal aspirates. In the first period of the study (1998–2007), viral detection of RSV and influenza virus was performed using immunofluorescence. In the second period of the study (2008–2018), viruses were detected with a rapid-antigen test (RSV only) or through PCR (RSV and influenza virus). Improved sensitivity of viral testing in the second study period leads to an increased number of cases among older children, as they have lower viral loads compared to young children. From the matched medical records, the following data were extracted for each virologically confirmed case: the duration of hospitalization, age and sex of the patient, underlying diseases (i.e., healthy, asthma, neuromuscular disease, immunosuppression, prematurity, chronic respiratory diseases, metabolic and gastrointestinal diseases), and disease severity. Based on medical records, only cases with positive findings and ongoing clinical signs of infection were included in the study. If the same child was admitted more than once during a 4-week span and still tested positive for RSV or influenza virus, only the first episode was included. Only cases where the primary address of the patient was within the Astrid Lindgren Children’s Hospital catchment area were included.
Demographic data
We used the smooth.spline function (with 10 degrees of freedom) implemented in R (version 4.3.2) to interpolate weekly birth rates. Within our transmission model, we divided the <1 year age class into 12-month age groups to more accurately capture aging among this age class. The remaining population was divided into 9 age classes: [1,2) years, [2,3) years, [3,4) years, [4,5) years, 5–9 years, 10–19 years, 20–39 years, 40–59 years and 60 years old and above. We estimated the net rate of immigration/emigration for each age group (detailed in S1 Text) to produce a rate of population growth and age structure similar to that of northern Stockholm. Data on age-specific contact rates were obtained from22 specifically, we used the POLYMOD contact matrix from the Netherlands, which had a similar contact pattern as Sweden.
Climatic data
The climatic variables used in this study were weekly temperature (in Celsius) and relative humidity (as a percentage) from July 1998 to June 2018. Weekly averages were calculated from the daily data. To incorporate the climate data into the RSV transmission model, we normalized the data to between −1 and 1.
The center of gravity and the intensity of RSV activity
The center of gravity of RSV activity for each season (\({G}_{s}\)) was measured as the mean epidemic week, with each week weighted by the weekly number of admissions (\({Y}_{s,w}\)), such that
$${G}_{s}={\sum }_{{{\rm{w}}}\in [1:52]}w\times {Y}_{s,w}/{\sum }_{w\in [1:52]}{Y}_{s,w}$$
(1)
where \(w\) is an index for the week of each epidemic season, \(s\). The center of gravity is a measure of the mean epidemic week, indicating the average timing of RSV activity. This measure has been used repeatedly in the literature for RSV and other pathogens9,22,54,55,56. The intensity of RSV activity for each season was determined by the maximum weekly number of hospitalizations. The rationale for using the maximum values is to capture the contrast between large and small epidemics observed in the data.
Wavelet analysis
We obtained the timing of RSV epidemics in each season based on phase decomposition obtained from wavelet analysis57,58. In the wavelet analysis, we used the 0.8–1.5 year and 1.5-2.5 year periodicity bands from the wavelet spectrum to extract weekly phase angles for the annual period and the biennial period, respectively.
Dynamic model description
Here, we used an age-stratified SIS (Susceptible-Infectious-Susceptible) model for RSV transmission dynamics. The model was proposed by Pitzer et al. 22 to study the environmental drivers of the spatiotemporal dynamics of RSV in the US. The model assumed individuals are born with protective maternal immunity, which wanes exponentially, leaving the infants susceptible to infection. We assumed a progressive build-up of immunity following up to four previous infections, based on data from birth cohort studies. Following infection with RSV, individuals develop partial immunity, reducing the rate of subsequent infection and relative infectiousness of the following infections. We also assumed that subsequent infections have a shorter recovery time compared to the primary infection. The model was described by a system of ordinary differential equations; see S1 Text for details.
As a sensitivity analysis, we also used a two-strain mechanistic model (Fig. S9) to assess its compatibility with our original single-strain model. We used an SIRS (Susceptible-Infectious-Recovered-Susceptible) model for influenza transmission dynamics and integrated the model into the RSV transmission dynamic model described above.
In the model, influenza dynamics are coupled to the RSV transmission model (Fig. 3) through two alternative interference pathways. First, we allow prior infection with one virus to change susceptibility to the other by a factor \(\theta\). Second, we allow coinfected individuals to have altered infectiousness at a rate \(\xi\) compared to infection with either virus, capturing changes in viral shedding when both pathogens are present.
We calibrated the two-strain model to both weekly RSV and influenza admissions from July 2009 to June 2016. We assumed the RSV parameters were fixed at their pre-pandemic estimates (see Parameter estimation below). We estimated six influenza transmission parameters: a seasonal amplitude parameter, a seasonal offset parameter, a baseline transmission rate parameter, a reporting fraction parameter to scale infections with influenza to hospitalized cases, duration of influenza immunity, duration of influenza infection, and four viral interference parameters. We estimated these parameters using a maximum likelihood approach, assuming the number of hospitalizations during each week was Poisson distributed with a mean equal to the model-predicted number times the estimated reporting fraction. Other parameter values are provided in S1 Table. The full two-strain model structure, its detailed description, and model results are included in the Supplementary Materials.
Parameter estimation
For the single-strain model without viral interference, we first fitted the transmission dynamic model for RSV to weekly RSV admissions from July 1998 to June 2008 (i.e., 11 seasons before the influenza pandemic). We estimated four parameters: a seasonal amplitude parameter (\(\alpha\)), a seasonal offset parameter \((\phi )\), a baseline transmission rate parameter (\({\beta }_{0}\)), a reporting fraction parameter (\(f\)) to scale RSV low respiratory disease cases to hospitalized cases. For the model testing the climate hypothesis, we estimated two additional parameters, including the seasonal amplitude parameters for temperatures (\({\alpha }_{{Temp}}\)) and relative humidity (\({\alpha }_{{RelHum}}\)). The force of infection is given by:
$$\lambda={\beta }_{0}(1+\alpha \cos (2\pi {vt}-\phi )+{\alpha }_{{Temp}}\times {Temp}+{\alpha }_{{RelHum}}\times {RelHum}){I}^{*},$$
(2)
where \({Temp}\) and \({RelHum}\) are normalized data of temperatures and relative humidity, and \({I}^{*}\) denotes all infection states. Note that we assumed \({\alpha }_{{Temp}}=\,{\alpha }_{{RelHum}}=0\) when testing the birth rate hypothesis. We estimated these parameters using a maximum likelihood approach, assuming the number of hospitalizations during each week was Poisson-distributed with a mean equal to the model-predicted number times the estimated reporting fraction. Other parameter values for the model were adopted from ref. 22, and they are provided in S1 Table. We seeded the model with one RSV-infected individual in each age group except the <1 year group, then used a burn-in period of 47 years to ensure the model reached a quasi-equilibrium steady state.
For the models with viral interference, we fitted the transmission dynamic model for RSV to weekly RSV admissions from July 1998 to June 2016 (i.e., all seasons across the study period for which we had influenza admissions data). We first estimated the parameters for RSV transmission (i.e., \({\beta }_{0},\alpha,\phi,f\)) and viral interference parameters using a maximum likelihood approach. We also estimated climatic parameters (i.e., \({\alpha }_{{Temp}},\,{\alpha }_{{RelHum}}\)) when all three factors (the observed birth rate change, climatic factors, and viral interference) were included. We determined the best model using Akaike Information Criterion (AIC) scores. To further quantify the uncertainty of viral interference parameters in the best-fitting model, we then applied a Sampling-Importance-Resampling method to estimate credible intervals (CIs) of those parameters. We used Latin Hypercube Sampling (LHS) to generate representative samples from a wide range of values for the parameter space \(\varPhi=({\xi }_{s},{\xi }_{p},{\xi }_{{sp}},\eta )\). We drew 100,000 samples from a uniform distribution \(U(-5,-1)\) for \({\log }_{10}({\xi }_{s},{\xi }_{p},{\xi }_{{sp}})\), and from a uniform distribution \(U(1,\,20)\) days for \(\eta\). Note that we drew samples using a log10-scale for the \(\xi\) parameters because we had no prior knowledge of the magnitude of these parameters. To enhance computational efficiency, we used informative priors for the RSV transmission parameters (i.e., \({\beta }_{0},\alpha,\phi,f\)), setting their mean values to the estimates from the maximum likelihood approach and assuming a small variance, such that \({\beta }_{0} \sim N(8.51,\,0.5)\), \(\alpha \sim N(0.27,\,0.1)\), \({\beta }_{0} \sim N(3.24,\,0.5)\), \(f \sim U({\mathrm{0,1}})\). Then, we generated forward simulations using the sampled parameter sets and fitted them to weekly RSV admissions from 1998 to 2016. We calculated the log-likelihoods of the model under each parameter set, assuming the number of hospitalizations during each week was Poisson-distributed with a mean equal to the model-predicted number times the estimated reporting fraction. We then normalized the log-likelihoods to assign weights to each parameter set and resampled (with replacement) 10,000 parameters from the joint distribution based on the weights.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
