Cross-immunity between strains explains the dynamical pattern of paramyxoviruses

We perform wavelet decomposition of all time series from the years 2002 to 2014 to analyze the data in time–frequency domain. The continuous wavelet transformations of the RSV and HPIV time series are shown in SI Appendix, Fig. S2. The significant peaks of RSV periodicity are around the 52-wk band, whereas peaks of HPIV-1 and HPIV-2 activity are around the 104-wk band. However, there is high power for HPIV-2 around the 52-wk band in the period from 2010 to 2013. In contrast, HPIV-3 does not exhibit any consistent period, although there is significant high power around 52 and 104 wk in the periods of 2008–2010 and 2010–2012, respectively.

The wavelet coherence (WC) (SI Appendix, Fig. S3) and cross-wavelet transformation (XWT) (SI Appendix, Fig. S4) between pairs of strains reveal how strongly they are correlated in the time and frequency domain. Both XWT and WC of RSV with HPIV serotypes 1 and 2 (SI Appendix, Fig. S3 A and B) indicate that HPIV serotypes 1 and 2 are strongly correlated with RSV in the 52-wk band and lead RSV throughout the entire period. In contrast, the WC between RSV and HPIV-3 (SI Appendix, Fig. S3C) shows several significant spots of correlation along the 52-wk band with RSV leading HPIV-3. WC between HPIV-1 and HPIV-2 (SI Appendix, Fig. S3D) shows consistent significant power in the 104-wk band in a perfect antiphase relation. The wavelet analysis of hMPV data also shows high power in the 52-wk band, and cross-wavelet analysis with RSV exhibits in-phase relation in the 52-wk band. Thus, XWT and WC between RSV and HPIV serotypes are strong preliminary indicators of interactions among these pathogens.

We extend the susceptible–infected–recovered (SIR) compartmental model (45) to include two different disease interaction mechanisms: cross-immunity and convalescence (detailed descriptions are given in Materials and Methods and SI Appendix). Cross-immunity acts by reducing the rate at which a host that has recently recovered from one pathogen may become infected with another (38). Convalescence acts by temporarily removing individuals after infection because of behavioral modification and thus, transiently reducing the susceptible pools of the infections (43). Our models include seasonally forced transmission dynamics and temporary immunity of acquired infection. Studies suggest that reinfection by RSV and HPIV is likely in infants and the elderly, although it remains largely undetected in the latter (19, 20). A detailed description and schematic (SI Appendix, Fig. S1) of the models are given in SI Appendix.

Both models can exhibit a variety of dynamics depending on the strength of interactions (46). For example, when the two viruses are in a similar parameter regime with no interaction, they often have annual and in-phase oscillations. Increasing the strength of cross-protection can induce oscillations with higher-order cycles and out of phase dynamics (SI Appendix, Fig. S5). Generally, either cross-immunity or convalescence has little influence on the interepidemic period (SI Appendix, Fig. S6). These theoretical results indicate that both mechanisms may generate broad dynamical features of our paramyxovirus case studies. In the next section, along with the no interaction model (as “null model”), we fit these two models to the time series of case reports to select the most parsimonious model to explain RSV, HPIV, and hMPV dynamics.

We begin parameter estimation by fitting the single-disease seasonally forced SIR model (i.e., no interaction model) for each strain separately [a detailed discussion of the observation model, baseline parameter values (SI Appendix, Table S1), parameter estimation, and model selection methodology is given in Materials and Methods and SI Appendix]. The single-disease SIR model gives a good fit to timing and size of RSV peaks, a poor fit to the biennial dynamics of HPIV-1 and HPIV-2 (the goodness of fit being negative), a moderately accurate fit for HPIV-3 (SI Appendix, Fig. S7 and Table S2), and good predictive power for hMPV (results in SI Appendix).

The two-disease cross-immunity model fits well to the paired datasets with RSV and each serotype of HPIV (Fig. 2), describing the high annual peak of RSV and the relatively small biennial peaks of HPIV (Fig. 2 A and B). The estimated strength of cross-protection (SI Appendix, Table S2) suggests that strong cross-protective immunity created by RSV on HPIV serotypes 1 and 2 shapes the biennial abundance pattern in these serotypes. The model also captures the behavior of the pair of RSV and HPIV-3, describing relatively lower annual peaks of HPIV-3. The cross-immunity model also captured the dynamics of the paired dataset of RSV and hMPV (SI Appendix, Fig. S9).

The two-disease convalescence model, in contrast, broadly fails to simultaneously explain the dynamics of RSV and serotypes of HPIV (SI Appendix, Fig. S8), although it exhibits a biennial oscillation at low convalescence rate and the annual peaks of HPIV-3 at a higher convalescence rate (SI Appendix, Table S2). This model does provide a relatively better fit for the hMPV data. However, the error profile of the convalescent parameter θ for all three estimations in SI Appendix, Fig. S13 indicates that the convalescence mechanism has less influence in explaining the paired RSV–HPIV datasets.

Altogether, the data suggest that the dynamics are shaped, in part, by interactions among strains and the cross-immunity hypothesis (SI Appendix, Table S2). Confidence intervals, best-fitting parameter values, and AIC index for all models are given in SI Appendix, Figs. S15–S18 and Tables S2 and S5.

Although HPIV-1 and HPIV-2 are less common among hospitalized cases, they could have a strong interaction among the potentially many more less severe cases. The outbreaks of these two serotypes are separated by the outbreak of RSV in alternate years. Therefore, we further extended the SIR model to a three-strain model to check if including the full web of potential interactions better explains the dynamics of these three pathogens. However, the three-disease model does not explain the data better than the two-disease cross-immunity model (SI Appendix, Table S3), perhaps because the resultant model has many parameters relative to the total amount of data available.

To investigate further statistical limitations on inference, we compute the correlation matrix among parameters (obtained from inverting the Hessian matrix) (SI Appendix, Fig. S11). We also represent the two-parameter profile (${\text{ε}}_P-{\text{ε}}_R$) of normalized error from each paired dataset (SI Appendix, Fig. S12). As is to be expected from this type of mechanistic model, there is substantial multicollinearity among parameters.

The data presented in the paper (Fig. 1) are weekly case reports of respiratory virus results (RSV and serotypes of HPIV from 2000 to July of 2014 and hMPV from 2007 to July of 2014) from all IHC laboratories (22 hospitals and a medical group with more than 185 clinics) in the states of Utah and Idaho in the intermountain west with a population of ∼3 million in 2011 (56). During the study period, direct respiratory sampling (mainly saline-assisted nasopharyngeal aspiration) for respiratory viral testing was performed on children admitted in the clinic or hospitals with respiratory complaints. The test methods were enzyme immunoassay for RSV, direct fluorescent antibody strain, culture, and/or PCR. Test methods differ in sensitivity and turnaround time. Testing algorithms vary across IHC facilities, with some facilities using automatic reflex testing to more sensitive methods (e.g., direct fluorescent antibody to PCR), whereas other facilities require a clinician order for a more sensitive method to be performed. Results of testing are aggregated into a single patient encounter, even if multiple tests are performed. The RSV epidemic year was defined to be from July 1 of the year to June 30 of the following year. This time period was chosen to place the beginning date close to the middle of the interepidemic period, approximately 6 mo from the average historical peak of the seasonal epidemic. This study was reviewed by the Institutional Review Boards of IHC and the University of Utah and determined by both organizations to be exempt.

The full set of equations and the descriptions of all three models (cross-immunity, convalescence, and no interaction) are given in SI Appendix. Here, we outline how the interaction mechanisms are implemented in the model framework.

The precise immunological mechanisms of cross-protective immunity between strains of paramyxoviruses are not yet known. Our models incorporate a parameter to reflect the strength of cross-immunity between strains that is assumed to hold for the same duration as for the primary pathogen but with lesser effect. We also assume that partial cross-immunity acts as reducing susceptibility to other pathogens. Usually, the acquired immunity from infection protects the individual from original pathogen, but the same immunity may work for other antigenically related pathogens, and hence, the host will be rarely infected by such strains. Thus, two strains interact by conferring some protection to the host against other strains. Here, we assume that cross-protective immunity has the same duration for the related and original pathogens. In our main analysis, we follow clinical data that suggest that immunity lasts for around 1 y. However, in the analysis in SI Appendix, we allow this parameter to vary freely to be estimated from time series data.

This interaction between strains arises by competition for susceptible hosts and the temporary removal of individuals after infection because of quarantining or convalescence (43, 44, 53). In respiratory viral infections, such as RSV, child patients are often sent home or to the hospital until they are recovered or cured. We incorporate a parameter for the decreased likelihood of acquiring another infection during the convalescent period.

We estimate transmission rate, seasonality, interaction parameters, and reporting probability along with initial conditions for all different models using the time series data from the years 2002–2014. The reporting probability (i.e., observation model) is modeled in the following way: if $T_{i,t}$ is the total number of newly infected in week t for pathogen i, and f is the reporting probability, then weekly case notification is assumed to have been drawn from a Poisson distribution with mean $f{T}_{i,t}$. We estimate parameters of two-disease models using pairs of time series and parameters of three-disease models using data of RSV, HPIV-1, and HPIV-2. We minimize the negative of log likelihood to estimate on the best parameter values that describe the case data and use the Akaike Information Criterion (57) to evaluate the parsimony of competing models. Detailed descriptions of the likelihood function and model selection methodology are given in SI Appendix.

We find confidence intervals of estimated parameters using generalized cross-validation (SI Appendix, refs. 8 and 9). In this case, we generate 100 samples by randomly partitioning the entire dataset into two subsets (training set and test set) and validate parameter estimates on the test set. Three different partitions (70%, 80%, and 90%) were used in each case. Details are in SI Appendix.

We also use rms error to calculate the mean deviance for error profile. We calculate rms error from 10 sample model fits and average these to calculate the mean deviance. Details are given in SI Appendix.