Data
Serology and viral activity in King County, WA
Our serology study was initiated in November 2021 as an add-on to the larger Seattle Flu Alliance surveillance study (SFA), which began in 2018 and carried out hospital and community-based surveillance to monitor the circulation of 26 respiratory pathogens45. The purpose of the serology arm of SFA was to document the immunity debt that may have accumulated in the first two years of the COVID-19 pandemic and explore changes in antibody levels coinciding with putative pathogen rebound in the 2021–22 winter. Cross-sectional sera were collected from convenience residual blood samples of 999 immunocompetent children seeking medical care, evenly distributed between the <5 yo and 6–10 yo age groups, through Seattle Children’s Hospital. We also sourced convenience residual sera from a blood donor bank in Seattle, Bloodworks Northwest, which collected specimens from 509 adults prospectively, 274 of which provided repeat samples46. We collected the same number of samples in two adult age groups, 20–64 and over 65 yo. We defined three sampling periods to titrate immunity at key points of the pandemic, after balancing the availability of samples and seasonality of pathogen circulation (Table 1). The first period, occurring in the summer of 2020, followed a mostly regular respiratory season in 2019–2020, with only samples from children <11 yo being collected. In subsequent periods, sera was collected from both adults and children. The second period in fall 2021 followed a year and a half of repressed pathogen circulation due to COVID-19 interventions and was presumed to mark a trough in immunity. The third period in spring-summer 2022 followed a respiratory virus season in which multiple pathogens rebounded (Fig. 1A).
Sera was analyzed using the multiplex meso scale discovery (MSD) electrochemiluminescence immunoassay and laboratory analysis was performed at the Fred Hutchinson Cancer Center21. The V-PLEX Respiratory Panel 3 IgG Kit has been used increasingly since the pandemic, primarily to characterize antibody responses to SARS-CoV-2 vaccination and infection47,48. We used the V-PLEX COVID-19 Respiratory Panel 3 IgG Kit to analyze antibody concentration in arbitrary units [AU] per milliliter against Flu A/Hong Kong/2014 H3, Flu A/Michigan/2015 H1, Flu A/Shanghai/2013 H7, Flu B/Brisbane/2008 HA, Flu B/Phuket/2013 HA, HCoV 229E Spike, HCoV-HKU1 Spike, HCoV-NL63 Spike, HCoV-OC43 Spike, MERS-CoV Spike, RSV Pre-Fusion F, SARS-CoV-1 Spike, SARS-CoV-2 N, SARS-CoV-2 S1 RBD, and SARS-CoV-2 Spike antigens (SARS-CoV-2 based on ancestral strain). We transformed the antibody concentration to log10 units. The King County sample collection and this study were approved by the Institutional Review Board of Seattle Children’s Hospital. This study was granted a waiver of consent since it used residual clinical samples and existing clinical data. The authors did not have access to identifiable information.
To provide context for serology, we sourced detailed respiratory surveillance data from Seattle Children’s Hospital (SCH) from January 2018 to January 2023. SCH performed 105,207 PCR tests among children <21 yo meeting criteria (symptom onset within the past 7–10 days and at least one respiratory symptom, e.g., cough, sore throat, shortness of breath) in the inpatient and outpatient setting during this time period. Children were tested for influenza and RSV using the Cepheid Xpert XpressR (Sunnyvale, CA) panel and for influenza A/H3N2, A/H1N1, HCoV 229E, NL63, OC43, and HKU1 using the Biofire FilmArrayR multi-pathogen respiratory viral panel (Biomerieux, Salt Lake City, UT). We presented monthly tests percent positive for these pathogens to visualize endemic respiratory pathogen circulation throughout the study period and additionally derived a COVID-19 time series from the King County COVID-19 Dashboard49. Because we were interested in how serology could improve understanding of the influenza rebound in the post-COVID-19 pandemic period, we obtained more resolved age-specific data on influenza healthcare encounters from Public Health—Seattle and King County (PHSKC). The dataset comprises age-specific (0–4, 5–19, 20–64, over 65 years) emergency department visits and hospitalization data reported by PHSKC to the Washington State Department of Health Rapid Health Information Network (RHINO)50. The number of King County hospitals that report to RHINO increased from 6 to 21 from 2017-2022. Given limited information on how the underlying surveillance catchment area changed, and to obtain a stable time series of influenza healthcare encounters, we assumed that reporting increased linearly with the number of reporting hospitals. This means that the time series was upscaled by the fraction of reporting hospitals relative to the final set of 21.
Serology in South Africa
Because we used a novel multiplex serological assay to assess immune dynamics in the King County, WA, population, we looked for confirmation of our findings with an independent dataset and assay. Accordingly, we sourced serologic data from a 3-year household cohort of influenza in South Africa22,51. We analyzed paired sera collected from 1684 individuals through the Prospective Household cohort study of Influenza, Respiratory syncytial virus, and other respiratory pathogens community burden and Transmission dynamics (PHIRST) in a rural and urban setting located in Mpumalanga and North West Province, South Africa, 2016–2018. Detailed study protocols have been published elsewhere22,51. Longitudinal samples were collected at study enrollment and before and after the influenza season annually in 2016, 2017, and 2018 for enrolled households, including 375 children and 653 adults (Table 1). Antibody levels were measured using the well-established hemagglutinin inhibition assays against influenza strains circulating during the study period, namely A/South Africa/2517/2016 EPI_ISL_230453, A/Singapore/INFIMH-16-0019/2016 EPI_ISL_225834, B/South Africa/R3037/2016 EPI_ISL_231726, and B/South Africa/R5631/17 EPI_ISL_17008503. Antibody concentration was transformed into log2(x/5) units. The South Africa PHIRST protocol was approved by the University of the Witwatersrand Human Research Ethics Committee (Reference 150808), and the US CDC’s Institutional Review Board relied on the local review (#6840). The protocol was registered on http://clinicaltrials.gov in 2015 (Reference NCT02519803) and participants provided individual written consent or assent prior to enrollment.
Modeling age-specific immunodynamics of respiratory pathogens
We first used simple descriptive statistics to assess changes in antibody concentration levels for individuals in King County by age group and pathogen throughout the study period. We summarized the serological data by 7 age groups: <1 yo, 1–2 yo, 3–4 yo, 5–10 yo, 18–49 yo, 50–64 yo, and 65+ yo, and calculated the percent change in average antibody levels between the 2020, 2021, and 2022 time points. Two-sample Kolmogorov–Smirnov tests were performed to assess whether the distribution of antibody measurements by age band and pathogen changed across annual time points (Supplementary Tables 2–3).
Next, we used a Bayesian hierarchical model developed to make epidemiological inferences from individual-level serological data using the R package serosolver20. The model estimates the joint posterior distribution of (1) underlying antibody kinetics following infection, including antibody boosting and waning, (2) infection histories for all individuals throughout the study period, and (3) the population infection probability over time, which is estimated independently for each age group. In short, the model estimates the set and timing of infections that are most consistent with an individual’s antibody profile, accounting for dynamic antibody responses following infection and measurement variability. We separately fit the model to the King County serological data using cross-sectional data from children <5 years and paired sera from adults for RSV, SARS-CoV-2 N, the seasonal coronaviruses HCoV 229E, NL63, OC43, and HKU1, and influenza subtypes A/H3N2, A/H1N1, B/Victoria, and B/Yamagata (Table 1). We did not model the 5–10 yo age group cross-sectional data due to the challenge of distinguishing between high initial antibodies and antibody boosts following exposure28. We also excluded infants <6 months old to ensure that our results were not biased by transplacental maternal antibody dynamics in young infants52.
We additionally fit the same model to serological samples paired across time for individuals from South Africa for three age groups, aligning with the age groups represented in the King County serological data, <5 yo, 5–10 yo, and 18+ yo, and for each influenza subtype: A/H3, A/H1, B/Victoria, and B/Yamagata (Table 1). Model outputs for the 10–18-year-old population in South Africa are additionally provided in Supplementary Table 4. The model was fit to two samples per individual in the South African dataset.
We adopted a simplified version of the serosolver modeling framework, where the observed antibody level Y for individual i at a given time t for age group a is the sum of antibody boosting \({\mu }_{s}\) that has accumulated since birth and wanes at rate r following exposure. We assume that estimated parameters are pathogen- and subtype-specific (p). An individual’s expected antibody level can then be modeled as:
$${Y}_{i,a,p,t}={\varSigma }_{j}{Z}_{i,j}({\mu }_{s,a,p}(1-r(t-j)))$$
(1)
where the distribution of observed antibody levels for a given age group and pathogen is normally distributed (truncated given the upper and lower limit of detection) with mean \({Y}_{i,a,p,t}\) and the set of antibody kinetics parameters to be estimated is defined by \(\varTheta=\{{\mu }_{s},r,\sigma \}\), where \(\sigma\) is the estimated measurement error. Individual infection histories are represented as vectors of binary latent states \({Z}_{i,j}\), indicating whether an individual was infected (\({Z}_{i,j}\) = 1) or not infected (\({Z}_{i,j}\) = 0) in a given time period j. We estimated infection histories back to birth, with the exception of SARS-CoV-2 antigens, which are estimated from 2020 to 2022. The model is run at an annual timescale, meaning an individual can be infected a maximum of once per season. Each infection event \({Z}_{i,j}\) is modeled by a Bernoulli trial that is conditional on the time-varying population probability of infection \(\phi\).
The serosolver framework allows for the implementation of different priors on the probability of infection \(\phi\) depending on the specific disease context. We implemented a beta prior on the probability of infection in each time period, which is appropriate when an individual’s time-varying probability of infection is governed by the population-level force of infection; as is the case for endemic respiratory pathogens20. In this model framework, the per-capita attack rate is beta-distributed, and the total number of infections experienced by each individual follows a binomial distribution. To set the priors, we inferred the distribution of influenza population attack rates from PCR-confirmed data collected through the PHIRST study. Annual attack rates were right-skewed with a mean of 10% standard deviation of 13%; we used this to define our infection probability prior for all endemic respiratory pathogens across both study locations. Thus the population probability of infection P(ϕ) is defined by a beta distribution in which we have set alpha = 0.37 and beta = 3.5 (mean: 0.10, standard deviation: 0.13) (Supplementary Fig. 4). Because the King County dataset reported influenza titers at the subtype but not strain level (given use of the MSD multiplex assay) while the South Africa dataset reported strain-specific HI titers over a 3-year period with limited antigenic changes, we did not incorporate an antigenic map into the analysis. Analyses considering longer time periods and antibody titers specific to different strains can include this information in the serosolver framework.
The joint estimation of all parameters given observed individual-level antibody data is defined in Eq. 2:
$$P(Z,\phi,\theta {|Y})\alpha {\prod }_{i=1}^{n}({\prod }_{t=1}^{{tmax}}f({Y}_{i,t}|{Z}_{i,},\theta )){\prod }_{t=1}^{{tmax}}P({Z}_{i}|\phi )P(\phi \left)\right.)P(\theta )$$
(2)
The model uses an adaptive Markov Chain Monte Carlo (MCMC) Metropolis-Hasting algorithm to infer \(P(Z,\phi,\theta {|Y})\). We performed each model run with a minimum of three chains and 500,000 iterations per chain and assessed post-run model convergence by ensuring \(\hat{R}\) was <1.1 for all parameters in each model run. Model diagnostics are provided in Supplementary Fig. 5.
Population-level influenza transmission model
Next, we integrated our serologic findings into a mechanistic transmission model to assess how age-specific immunodynamics may impact transmission, particularly in the post-COVID-19 period when endemic pathogens rebounded. We focused this analysis on influenza for which the age dynamics of the rebound in 2021–23 was particularly poorly explained19 and for which we had more serological data. We considered two variants of a susceptible-exposed-infected-recovered-susceptible (SEIRS) model, a null model in which we assumed uniform duration of immunity across the population, and a tiered immunity model, in which we assumed that children <5 yo experienced a different waning rate than older individuals (to be estimated), inspired by our serological analysis. We calibrated these models to monthly age-specific influenza healthcare encounter data for King County, as reported to the RHINO’s syndromic surveillance system, which collates visits based on discharge diagnosis codes. The full set of model equations is given by Eqs. 3–6:
$$\frac{d{S}_{i}}{{dt}}=\mbox{-}{{{{\rm{\lambda }}}}}_{{{{\rm{i}}}}}\,{S}_{{{{\rm{i}}}}}+(v/{r}_{i}){R}_{i}$$
(3)
$$\frac{d{E}_{i}}{{dt}}={{{{\rm{\lambda }}}}}_{{{{\rm{i}}}}}\,{S}_{{{{\rm{i}}}}}\mbox{-}{{{\rm{\delta }}}}{E}_{i}$$
(4)
$$\frac{d{I}_{i}}{{dt}}={{{\rm{\delta }}}}{E}_{i}\mbox{-}{{{\rm{\gamma }}}}{I}_{i}$$
(5)
$$\frac{d{R}_{i}}{{dt}}={{{\rm{\gamma }}}}{I}_{i}\mbox{-}(v/{r}_{i}){R}_{i}$$
(6)
where S is the susceptible population; E is exposed; I is infectious; R is recovered, and i represents a given age cohort (monthly age cohorts for children ≤ 5 yo and 5-year age cohorts for persons aged 5–75 yo). 1/γ is the infectious period, which is set to 2.27 days, and 1/δ is the latent period, which is set to 2 days based on literature53,54,55,56,57.
To integrate age-specific immunodynamics in this model, we explored different age-specific immunity durations. Let \(1/v\) be the baseline duration of immunity and ri the proportion change in duration of immunity for age i compared to baseline. We fix \(1/v\) = 4 years based on prior literature17,58,59. In the tiered immunity model, we allow the proportion change in immunity duration ri to vary in children under 5 yrs, while retaining the longer immunity duration in older individuals (i.e., ri = 1 for all ages above 5 yrs, so that duration of immunity is 4 years for individuals >5 yo). We then consider a range of ri values that capture putative durations of immunity in children <5 yo that are both longer and shorter than in older individuals and estimate the value of ri that is most consistent with the incidence data at the calibration stage. In the uniform immunity model, ri is set at 1 no matter the age group (i.e., immunity lasts on average 4 years for all ages).
\({\lambda }_{i}(t)\) i the influenza force of infection on a susceptible individual in age class i and j is given by Eq. 7:
$${\lambda }_{i}\left(t\right)={c\left(t\right)\beta }_{0}(1+{\beta }_{1}\cos \left(\frac{2\pi t}{12}+\phi \right))\frac{1}{{N}_{i}}\sum\limits_{j=1}^{75}{M}_{i,j}{\omega }_{j}{I}_{j}(t)$$
(7)
Here t is time in months, λi(t) is the monthly transmission rate, β0 is the baseline transmission coefficient, and β1 and ϕ represent the amplitude and phase shift respectively that capture the seasonality in transmission in King County, to be estimated from the data. Mi,j describes the contact matrix between individuals in age groups j and i, where Mi,j is measured in contacts per day per person, and children experience a higher absolute number of contacts than adults. We used an expanded version of the contact matrix originally described by Mossong et al., which describes population mixing patterns for several European countries and aggregated daily contacts to monthly to align with the model’s temporal scale60. ωj represents different infectiousness by age, with children under 5 years assumed to be 2x more infectious than older age groups61. c(t) is the time-varying strength of the pandemic-related control periods, representing the percent reduction in transmission due to NPIs, where 0 <c(t) < 1 and c(t) is to be estimated at the calibration stage.
Our model generates the number of new infections, which is given by δ\({E}_{i}\), while our observations are healthcare encounters. To connect the two, we write that Hi~ Poisson(hiδ\({E}_{i}\)), where Hi is the number of reported healthcare encounters at time t in age group i, δ\({E}_{i}\), is the number of new infections at time t, and hi is the proportion of infections in this age group that result in reported hospitalizations, where hi is to be estimated from the data. The full parameter table can be found in Supplementary Tables 5, 6.
We calibrated the uniform and tiered immunity models in two steps. In the first step, we calibrated the model to age-specific influenza infections for the pre-pandemic period (October 2017–March 2020). In the absence of observed population-level influenza infection time series, we created synthetic time series based on reported healthcare encounters. We estimated that in a given influenza season, 25–35% of children and 15% of adults are infected29,62,63,64. Using the population size of King County in 2021, we scaled up the observed age-specific monthly healthcare encounter time series to align with realistic incident infection curves, creating a synthetic infection time series that reflects both the true seasonality trends captured by the healthcare encounter time series and realistic population attack rates. To calibrate the baseline uniform immunity model, we fit the transmission parameters (βo, β1, and ϕ) and age-specific hospitalization rates hi. In the age-specific tiered immunity model, we fit one additional parameter, the relative duration of immunity in children vs adults ri, to test whether there was support in the data for a different (shorter or longer) duration of immunity in children.
In the second calibration step, we fit two control periods beginning in March 2020 to estimate the effect of NPIs on influenza transmission in the first two years of the pandemic. We were particularly interested in how these NPI affected the influenza rebound in the 2022–2023 season, where influenza hospitalization burden was higher than expected, particularly in children. We allowed for two different NPI periods to account for differences in overall transmission impacts of the first year of the pandemic (March 2020–March 2021), when more stringent NPIs were implemented in King County, and the second year (April 2021–April 2022), ending when Omicron-related precautions were subsiding. All parameters were estimated using maximum likelihood. We also calculated the root mean square error (RMSE) on age-specific hospitalization time series to assess the fit of the uniform and tiered immunity models (RMSE = 3792 vs 3920, respectively). Lastly, we performed a Chi-square test to analyze changes in the age distribution of observed influenza healthcare encounters between the pre- and post-pandemic periods. All analyses were conducted in R Studio Version 2024.12.1+563 (2024.12.1+563)65.
Ethics and inclusion statement
We have complied with all relevant ethical regulations and use only de-identified data. The research process included local researchers from both South Africa and King County, WA, with all parties participating in the study design, data ownership, and authorship.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
