Interactions among common non‐SARS‐CoV‐2 respiratory viruses and influence of the COVID‐19 pandemic on their circulation in New York City

Abstract Background Non‐pharmaceutical interventions (NPIs) and voluntary behavioral changes during the COVID‐19 pandemic have influenced the circulation of non‐SARS‐CoV‐2 respiratory infections. We aimed to examine interactions among common non‐SARS‐CoV‐2 respiratory virus and further estimate the impact of the COVID‐19 pandemic on these viruses. Methods We analyzed incidence data for seven groups of respiratory viruses in New York City (NYC) during October 2015 to May 2021 (i.e., before and during the COVID‐19 pandemic). We first used elastic net regression to identify potential virus interactions and further examined the robustness of the found interactions by comparing the performance of Seasonal Auto Regressive Integrated Moving Average (SARIMA) models with and without the interactions. We then used the models to compute counterfactual estimates of cumulative incidence and estimate the reduction during the COVID‐19 pandemic period from March 2020 to May 2021, for each virus. Results We identified potential interactions for three endemic human coronaviruses (CoV‐NL63, CoV‐HKU, and CoV‐OC43), parainfluenza (PIV)‐1, rhinovirus, and respiratory syncytial virus (RSV). We found significant reductions (by ~70–90%) in cumulative incidence of CoV‐OC43, CoV‐229E, human metapneumovirus, PIV‐2, PIV‐4, RSV, and influenza virus during the COVID‐19 pandemic. In contrast, the circulation of adenovirus and rhinovirus was less affected. Conclusions Circulation of several respiratory viruses has been low during the COVID‐19 pandemic, which may lead to increased population susceptibility. It is thus important to enhance monitoring of these viruses and promptly enact measures to mitigate their health impacts (e.g., influenza vaccination campaign and hospital infection prevention) as societies resume normal activities.


| INTRODUCTION
Viral respiratory infections are one of the leading causes of disease in humans. Each year, numerous respiratory viruses co-circulate in the population, causing substantial public health burden 1 and economic loss. 2,3 Previous studies have suggested that respiratory viruses may interfere with and change the risk, timing, or natural history of infection of one another. 4 For instance, in 2009, seasonal epidemic of respiratory syncytial virus (RSV) in Israel was temporarily delayed due to the A(H1N1) pandemic. 5 Potential mechanisms including competitions within hosts (e.g., infecting cells) and population-level interactions have been proposed to explain such virus interactions. 4 However, the specific interactions among different respiratory viruses and the impact on their collective transmission dynamics have not been well characterized.
Before the COVID-19 pandemic, influenza was the foremost public health concern among all respiratory infections. As such, much research effort has been devoted to understand the transmission dynamics of influenza viruses and their interactions with other respiratory viruses. 4 In contrast, many other infections such as human endemic coronaviruses have received far less attention. In addition, previous studies tend to ignore the different subtypes of respiratory viruses and only examine interactions at the level of genus. However, subtypes from the same virus group may have different seasonality (e.g., the four subtypes of parainfluenza viruses) and competitions within genus tend to be more intense (e.g., influenza A[H1N1] and A [H3N2]). 6 As such, combining all subtypes of a virus regardless of their circulation patterns may mask the true interactions.
Following its emergence in late 2019, the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has spread to 214 countries and territories, causing the coronavirus disease 2019 (COVID- 19) pandemic. 7 The widespread prevalence of SARS-CoV-2 may affect the circulation of other respiratory viruses, via virus interactions. In addition, during the COVID-19 pandemic, non-pharmaceutical interventions (NPIs) such as social distancing, school closure, travel ban, and mask-wearing that aimed to mitigate SARS-CoV-2 transmission had also limited the transmission of other respiratory viruses. Seasonal respiratory viruses such as influenza were found to be at low circulation during the 2020 respiratory virus season amid the COVID-19 pandemic. 8,9 However, to what extent the COVID-19 pandemic has impacted the circulation of common respiratory viruses and potential differences by virus remains unclear.
In this work, we analyzed incidence data for seven groups of respiratory viruses in New York City (NYC) before and during the COVID-19 pandemic. We first used elastic net regression to identify potential interactions among the viruses at the subtype level (13 in total). We further hypothesized that strong interactions should allow models incorporating the relationship to more accurately predict the incidence of related respiratory viruses. To examine the robustness of the found interactions, we thus built and compared the performance of Seasonal Auto Regressive Integrated Moving Average (SARIMA) models with and without the interactions. Lastly, we used the best-performing models to estimate the impact of the COVID-19 pandemic on circulation of each of the 13 respiratory viruses.

| Virus surveillance data
The virus surveillance data were collected from a subset of NYC participating laboratories that test for multiple respiratory viruses in addition to influenza and RSV. Tested viruses included adenovirus (Adv), coronavirus (CoV), human metapneumovirus (HMPV), rhinovirus (RV), parainfluenza (PIV), RSV, and influenza virus (IV) overall and by subtype ( week the next year) increased over time with the expansion of testing ( Figure S1). Previous studies have used the percent positivity (i.e., the ratio of positive samples to total samples) to account for such changes in testing. However, during the summer when fewer respiratory viruses are in circulation and fewer people seek testing, the much smaller sample sizes ( Figure S1) tend to result in very high percent positivity for some viruses (e.g., RV in Figure S2), which may not reflect the true circulation levels of these viruses in the population. Thus, to account for the testing time trend and better represent the viral circulation levels, here, we adjusted the weekly incidence by multiplying the ratio of the number of samples tested during the season to that number in season 2015-2016 (the first season of data collection). For the COVID-19 pandemic period (Week 10 of 2020 to Week 20 of 2021; note the first COVID-19 case was reported in NYC during Week 10 of 2020 10 ), the number of samples tested each week fluctuated substantially from week to week (in particular, there were initial increases followed by a large drop during early weeks of the COVID-19 pandemic period; see Figure S1). To account for this short-term fluctuation, we adjusted the incidence during COVID-19 pandemic period week by week relative to the corresponding week during season 2015-2016. A comparison of the adjusted weekly incidence using this method and the percent positivity for each virus is shown in Figure S2.

| Selection of key viral interactions during the pre-COVID period
We used elastic net regression models 11 to preliminarily identify, for each virus, the set of other viruses consistently included as interacting viruses, during Week 40 of 2015 to Week 9 of 2020, that is, before the first COVID-19 case was reported in NYC. 10 To avoid spurious correlation due to a common winter-time seasonality shared by some viruses, 12 we first used a linear regression model to identify and remove the seasonal trend for each virus. The model took the following form: where Y is the adjusted weekly incidence (see Section 2.1) and Week_of_Year is an indicator variable for each week of the calendar year (1:52 for annual seasonal cycle and 1:104 for biennial cycles). For each virus, we fitted both annual and biennial cycles and used the adjusted R 2 to determine the most likely seasonal cycle for each virus for removal of seasonal trend.
We then regressed on the detrended time series (i.e., the residuals after removing the seasonal trend) with an elastic net penalty: where Y i, detrended is the detrended weekly incidence for a given virus of interest, i, and X j, detrended is the detrended weekly incidence for other viruses (i.e., for any j ≠ i); β j is the corresponding regression coefficient. Similar to lasso (i.e., least absolute shrinkage and selection operator) and ridge regressions, elastic net shrinks regression coefficients by imposing a penalty on their size. However, instead of penalizing by the sum of absolute coefficients (L1-lasso penalty) or the sum-of-squares coefficients (L2-ridge penalty), elastic net regression penalizes with both L1 and L2 and the penalty function is formulated as where λ controls the amount of shrinkage, α controls the distribution between L1 and L2, and β j represents the regression coefficients that minimize the penalty (i.e., formula 3). Since elastic net shares traits of both ridge and lasso regression, while it selects covariates like lasso, it also allows coefficients of correlated covariates to shrink together and provide a more stable selection result.
We fitted 500 elastic net regressions with tenfold cross-validation and pooled all interactions selected at least in half of the 500 runs for further testing (see the next section). Thirteen models were developed, one model each for Adv, HMPV, RV, RSV, and IV and for each of four subtypes of CoV and PIV. For influenza, due to the more erratic circulation pattern of different types and subtypes and short study period with available data (i.e., 5 years), we combined all types and subtypes.

| Testing the identified interactions
To further test the identified interactions from the elastic net regression models, we examined if inclusion of any of the found interactions in an SARIMA model with exogenous variables (SARIMAX) model 13 where y 0 is a differenced time series with a seasonal differencing degree of 1 and non-seasonal differencing degree d; p is the order of autoregressive (AR) model; q is the order of moving average (MA) model; ϕ i (i = 1, …, p) are the coefficients for the AR terms; Þare the coefficients for the MA terms; and ε is the error term. Parameters d, p, and q were selected based on model fit, and the values from the best-fit model for each virus are shown in Table S1.
Similarly, after the seasonal differencing, the SARIMAX model was formulated as where

| Model validation
We tested the SARIMA and SARIMAX models using data prior to the (either SARIMA or SARIMAX) were first fit to the training subset; the trained models were then used to estimate the weekly incidence for each virus during the testing period for out-of-fit model validation.
We compared model performance based on relative Root Mean Square Error (rRMSE) during the testing period. In addition, as the weekly incidence tended to be low, we also evaluated the models based on the cumulative incidence during the testing period and the 95% prediction interval; if the observed cumulative incidence fell within the 95% prediction interval, the model was deemed accurate.

| Estimating the impact of COVID-19 pandemic on circulation of non-SARS-CoV-2 viruses
We used the validated models (SARIMA or SARMIAX) to generate counterfactual estimates for each virus during the pandemic periodthat is, the expected cumulative incidence should there be no pandemic. To enhance model performance, we refitted the validated models using data during the entire pre-COVID period (i.e., through Week 9 of 2020) and used them to predict the incidence for each virus during the COVID-19 period (here Week 10 of 2020 to Week 20 of 2021). We then compared the model counterfactual estimates of cumulative incidence during the COVID-19 period (C counterfactual ) to the observations (C observed ) to estimate the impact of the COVID-19 pandemic. We computed the percent reduction in cumulative incidence due to the COVID-19 pandemic for each virus as 3 | RESULTS

| Respiratory virus circulation before and after the introduction of SARS-CoV-2
During the pre-COVID period, influenza viruses were the most commonly detected in our dataset (up to 400 adjusted case counts per year (Figures 1 and 2). In contrast, RV cases were detected year-round and tended to have two comparable epidemics each year-one in the winter and one in the summer (Figure 1). PIV cases were also detected throughout the year, but the outbreak patterns were less obvious ( Figure 3). For the same virus group, different subtypes tended to alternate in circulation and recurred biennially with irregular peaks (e.g., the four coronaviruses in Figure 2; PIV-1 and PIV-2 in Figure 3).
In addition, among the four coronaviruses, weekly case counts were highly correlated for virus pairs belonging to different genera (r = 0.82 between CoV-OC43 and CoV-229E and 0.76 between CoV-NL63 and CoV-HKU). In contrast, the circulation of different influenza types and subtypes and PIV-4 appeared less regular, with few cases detected in some years.

| Potential virus interactions
The elastic net regression combining with the SARIMAX model selection identified several potential associations among the respiratory viruses (Tables 1 and S2). However, the patterns of found interactions were not readily clear (Table S2) (Table S2 and Figure S2).
Further examination using the time series models indicated that, for three coronaviruses (i.e., CoV-NL63, CoV-HKU, and CoV-OC43), RV, PIV-1, and RSV, the inclusion of the identified interacting viruses in the SARIMAX model generated more accurate out-of-fit estimates than the SARIMA model (Table S3). However, this improvement was not substantial (2.31-15.85% reduction in rRMSE, see Table S3), likely because epidemics of these respiratory  Table 1) were similarly accurate for both time series models. As such, below we present results from both models.

| Impact of the COVID-19 pandemic on non-SARS-CoV-2 viruses
Most respiratory viruses included here appeared to have had lower circulation during the COVID-19 pandemic than would be expected (comparing the observed incidence and model counterfactual estimates in Figures 1-3). This is likely due to the continued NPIs period-the observed cumulative incidence during this period fell outside the model predicted 95% intervals and the estimated mean reduction was around 70-90% for these viruses (Table 2). In contrast, Adv and RV appeared to be less affected ( Figure 1A,D). Incidence of Our analysis found multiple potential interactions among the 13 viruses and that most associations were positive, before the COVID-19 pandemic (Table S2)  For RV, even though viral activities were also lower during the COVID-19 period, case increases were observed during the summer of 2020 when NPIs were relaxed ( Figure 1D). Similarly, continued transmission of RV was reported during and after the 2009 influenza A(H1N1) pandemic. 22 Rhinovirus infections occur in most age groups, and infection of one serotype confers little immune protection against others. [23][24][25] This wider infection demographics and the large breadth of RV serotypes (around 160 discovered by 2018 24 ) likely facilitated its transmission locally in the population.
Our study has several limitations. First, the data analyzed here are a subset of all tests done in NYC (i.e., only those from laboratories using the expanded respiratory panel tests) and thus may not be fully representative of the entire population. Second, even though the selection criteria have not changed during the study period, underlying patient characteristics may differ among specimens tested before and during the COVID-19 period, due to changing medical seeking behaviors in response to COVID-19 (e.g., people with mild respiratory symptoms may be more likely to seek testing at the early stage of the pandemic due to concern of COVID-19); this in turn may temporally change the composition of underlying sample population. Third, fewer specimens were tested each week during the early phase of the pandemic due to limited testing supplies and human resources; this reduced sampling likely increased model uncertainty. Fourth, the iden-

PEER REVIEW
The peer review history for this article is available at https://publons. com/publon/10.1111/irv.12976.

DATA AVAILABILITY STATEMENT
The virus surveillance data were used with permission under a Data Use and Nondisclosure Agreement between the New York City Department of Health and Mental Hygiene and Columbia University.