How Interactions during Viral–Viral Coinfection Can Shape Infection Kinetics

Respiratory viral infections are a leading global cause of disease with multiple viruses detected in 20–30% of cases, and several viruses simultaneously circulating. Some infections with unique viral copathogens result in reduced pathogenicity, while other viral pairings can worsen disease. The mechanisms driving these dichotomous outcomes are likely variable and have only begun to be examined in the laboratory and clinic. To better understand viral–viral coinfections and predict potential mechanisms that result in distinct disease outcomes, we first systematically fit mathematical models to viral load data from ferrets infected with respiratory syncytial virus (RSV), followed by influenza A virus (IAV) after 3 days. The results suggest that IAV reduced the rate of RSV production, while RSV reduced the rate of IAV infected cell clearance. We then explored the realm of possible dynamics for scenarios that had not been examined experimentally, including a different infection order, coinfection timing, interaction mechanisms, and viral pairings. IAV coinfection with rhinovirus (RV) or SARS-CoV-2 (CoV2) was examined by using human viral load data from single infections together with murine weight-loss data from IAV-RV, RV-IAV, and IAV-CoV2 coinfections to guide the interpretation of the model results. Similar to the results with RSV-IAV coinfection, this analysis shows that the increased disease severity observed during murine IAV-RV or IAV-CoV2 coinfection was likely due to the slower clearance of IAV-infected cells by the other viruses. The improved outcome when IAV followed RV, on the other hand, could be replicated when the rate of RV infected cell clearance was reduced by IAV. Simulating viral–viral coinfections in this way provides new insights about how viral–viral interactions can regulate disease severity during coinfection and yields testable hypotheses ripe for experimental evaluation.

Data on viral-viral coinfections are somewhat limited, but some studies have begun evaluating the outcomes of coinfection with commonly observed viral pairs (e.g., IAV and RSV [26][27][28][29][30][31], IAV and PIV [32], IAV and RV [28,[33][34][35], RSV and RV [28], RSV and HMPV [36], IAV and severe acute respiratory syndrome (SARS)-CoV-2 [37][38][39][40][41][42][43], RSV and SARS-CoV-2 [43], and RV and SARS-CoV-2 [44]). Collectively, these studies observed diverse outcomes of viral-viral coinfections with some interactions resulting in enhanced spread of one or both viruses within the respiratory tract, while other mechanisms seem to work in an inhibitory manner. For example, PIV-2 enhanced cell-to-cell fusion through the expression of its surface glycoproteins, which boosted viral spread between cells, and increased IAV titers but not PIV titers [32]. On the other hand, infection with IAV limited a concurrent RSV infection by promoting intracellular competition for proteins or amino acids needed for the successful replication of both viruses within cell cultures [11,26]. In vivo infections in animal models support the exclusion of RSV by IAV and suggest that RSV prior to IAV decreases disease severity [27,29,30]. A similar competitive exclusion was observed in 3D tissue cultures during coinfection with RSV and HMPV, where HMPV was inhibited without any effect on RSV [36]. The reduction was accompanied by higher Type I and III interferon (IFN) responses [36]. IFN-mediated effects were also implicated in RV-IAV coinfection where RV-induced IFN protected against subsequent IAV infection within differentiated airway cell cultures [34].
The relative timing between viruses and the order in which the viruses infect the host seem to contribute to differing disease outcomes [27,29,33,37,39,40,42]. Interestingly, the exclusion effect during RSV-HMPV coinfection was more robust during a concurrent infection compared to that in an infection where HMPV followed RSV after 2 days. This is in contrast to other viral-viral coinfections where infections separated by 2 to 5 days had more robust effects [33,37,42]. For example, RV attenuated IAV-mediated disease severity and reduced IAV titers when RV infection had occurred 2 days before IAV, but the effect was reduced during simultaneous infection [33]. Conversely, animals coinfected with IAV 2 days before RV experienced greater disease severity [33]. Similar outcomes occurred in animals coinfected with IAV 3 days before SARS-CoV-2 [37].
These empirical studies illuminate the breadth of interactions that lead to diverse outcomes of respiratory viral-viral coinfections and the need for mathematical methods that could dissect complex, time-dependent, and potentially nonlinear mechanisms, as they do for viral-bacterial coinfection [45]. One study on viral-viral coinfections suggested that a faster-replicating virus would outcompete other viruses in a scenario where two viruses were competing for epithelial cells [46]. However, it is possible that viruses infect different cells or infect different areas of the respiratory tract. In addition, as noted above, they may inhibit or enhance other processes (e.g., replication rates and/or immune responses) and ultimately modulate disease. Thus, expanded modeling infrastructures are needed and are a focus of this study. Unfortunately, most studies on viral-viral coinfections lack quantitative information on viral loads and/or host immune responses that are needed to effectively use mathematical approaches, which typically include fitting a mechanistic model to data. However, weight loss, which is a measure of disease severity, is tracked in most murine studies, and our recent work showed that mathematical models can accurately connect animal weight loss to infection kinetics during monoinfections and viral-bacterial coinfections [47,48]. These links allow for us to better interpret weight-loss data, and afford the ability to assess mechanisms with limited data to modeling studies such as this one.
To address gaps in understanding viral-viral coinfections, we assessed RSV, RV, and SARS-CoV-2 coinfections with IAV using two mathematical models. We predicted potential underlying mechanisms of viral interference and cooperation and assessed how different viral orders, timings, and pairings affected the infection dynamics and disease severity. The results provide important insights into divergent outcomes, in addition to generating novel hypotheses regarding why certain viral orders enhance or reduce disease severity.

Data for RSV-IAV Coinfection in Ferrets
Viral-load data were digitized from a study where ferrets were monoinfected or coinfected with a long strain of RSV and/or influenza A/ Tasmania/2004/2009 (A[H1N1]pdm09; IAV) [27]. Briefly, groups of 4 ferrets were intranasally infected with 3.5 log 10 50% tissue culture infectious dose (TCID 50 ) of IAV, 5.0 log 10 plaque-forming units (PFU) of RSV, or IAV followed 3 days later by RSV. Viral RNA copy number per 100 µL of nasal wash was measured daily for 14 days post infection (pi).

Data for IAV, RV, and SARS-CoV-2 Infections in Humans
Viral-load data were digitized from studies where humans were experimentally or naturally monoinfected with IAV [49], RV [50], or SARS-CoV-2 [51]. For each, we chose one patient for ease of investigating potential coinfection dynamics. For IAV, human volunteers were experimentally infected intranasally with 4.2 log 10 TCID 50 of influenza A/Hong Kong/123/77. Nasal washes were collected daily for 7 days, and infectious viral titers were determined via TCID 50 . The used data here were from Patient 4 due to this individual having clear viral growth, peak, and decay phases. For RV [50], 14 human volunteers were intranasally infected with 2.4 log 10 TCID 50 of RV, which was reported as the geometric mean. Nasal washes were collected daily for 5 days. For SARS-CoV-2 [51], the data were from naturally infected patients. Viral loads were sampled from throat swabs and measured in RNA copies/mL. All samples were taken approximately 2-4 days after the symptoms. We used Patient 8 due to their clear viral-load dynamics, and assumed that the infection had initiated 5 days before the onset of symptoms.

Data for IAV Coinfection with RV or SARS-CoV-2 in Mice
Weight-loss data were digitized from a study where BALB/c mice had been intranasally infected with 7.6 × 10 6 TCID 50 of RV1B and/or 100 TCID 50 of influenza A/Puerto Rico/8/1934 (PR8) [33]. Coinfections were initiated simultaneously or sequentially at a 2 day intervals (IAV-RV or RV-IAV) [33]. Weight loss was measured daily for 14 days.

Mathematical Model of Viral Monoinfection
To describe the dynamic interactions between epithelial cells and the virus during monoinfection, we used the viral kinetic model in Equations (1)-(4) (reviewed in [45,52]). Briefly, in the model, target (epithelial) cells (T) are infected by the virus (V) at a rate βV per day. Once virus is internalized, the cell undergoes an eclipse phase (E) during which infected cells do not yet produce the virus. The cells then transition to the infectious phase (I) at a k rate per day. Productively infected cells are cleared at a δ rate per day. The virus is produced at a p rate per cell per day and cleared at a c rate per day.

Target Cell Competition
To model viral-viral coinfections, we expanded the model in Equations (1)-(4) using two hypotheses ( Figure 1). The first hypothesis assumed that two viruses (V 1 and V 2 ) competed for target cells (T) ('target cell competition model'; Equations (5)-(8)) [46]. In this model, a single equation for target cells (T) was used where each virus could infect these cells at rates β 1 V 1 per day and β 2 V 2 per day. All other equations were equivalent to those in the monoinfection model. Subscripts i = 1, 2 denote each virus.  (5)- (8) and Equations (9)- (12). In the 'target cell competition' model, two viruses (V 1,2 ) interact indirectly by competing for target cells (T) [46]. In the 'target cell partitioning' model, two viruses do not interact, and each has their own pool of target cells (T 1,2 ). In each model, target cells become infected by the virus at rates β i V i , where subscript i = 1, 2 denotes the specific rates to V 1 and V 2 , respectively. Infected cells enter an eclipse phase (E i ) and transition to producing the virus at rate k i . Productively infected cells (I i ) produce the virus at rate p i and are cleared at rate δ i . The virus is cleared at rate c i . Direct interactions were implemented by increasing and/or decreasing one or more of the rates using functions α(V i ) and/or ζ(V i ), respectively, due to the other virus.

Target Cell Partitioning
The second hypothesis assumed that each virus had its own pool of epithelial cells to infect (T 1 and T 2 ) because viruses may preferentially infect certain cell types [53][54][55][56][57] or be present in a different areas of the respiratory tract ('target cell partitioning model'; Equations (9)- (12)). All other equations remained the same but were distinct for each virus, resulting a total of 8 equations. Subscripts i = 1, 2 denote each virus.

Modeling Viral-Viral Interactions
To assess the effect of one virus on another, we used functions that enhanced (α(V i ) (Equation (13)) or inhibited (ζ(V i ); Equation (14)) a particular infection process (i.e., rates of viral infection (β i ), viral production (p i ), infected cell clearance (δ i ), or viral clearance (c i )).
Parameter κ (per RNA/mL or TCID 50 /mL) is the strength of the interaction.

Quantifying the Relative Change in Total Virus
To quantify changes in the viral loads and the total viral burden as a consequence of an interaction during coinfection, we calculated the relative change in total virus (i.e., the area under the curve (AUC)) of viral load using Equation (15), where V coinf AUC is the AUC for the coinfection and V single AUC is the AUC for the monoinfection. The AUC was calculated using the Python function scipy.integrate.trapz.

Quantifying Disease Severity
To quantify the percentage of the lung infected by the virus that related to animal weight loss [47], we calculated the cumulative area under the curve (CAUC) of the infected cell dynamics [47] using the Python function scipy.integrate.cumtrapz.

Parameter Estimation
For model fits to the ferret data, parameters were estimated using a nonlinear mixedeffect modeling (NLME) and stochastic approximation expectation minimization (SAEM) algorithm implemented in Monolix 2019R1 [58]. In the NLME approach, each individual parameter is drawn from a log-normal distribution and written as θ i = θe η i , η i = N (0, ω 2 i ), where θ denotes the median value of the parameter in the population, and η i denotes the random effect that accounts for the interindividual variability of the parameter within the population. Interindividual variability was allowed for all the estimated parameters with the assumption of no correlation and applying an additive residual error model for log 10 viral loads. For model fits to the human data, parameters were estimated using scipy.optimize.minimize in Python. Here, only a single patient monoinfected with a virus was used.
The initial number of target cells (T 0 ) was set to 5 × 10 7 cells for ferrets and 2 × 10 8 cells for humans. Similar to our previous studies [47,59,60], we fixed the initial number of infected cells (E 0 ) to 3.1 × 10 3 cells for IAV infection and 1.0 × 10 5 cells for RSV infection in ferrets [27] and 1 × 10 2 cells for all infections in humans. We considered other values of E 0 and found no significant differences in estimated parameters, which is consistent with our prior studies [47,59,60]. The initial number of productively infected cells (I 0 ) and the initial free virus (V 0 ) were set to 0.
The duration of the eclipse phase (1/k) for each virus was kept within a biologically feasible value, and set to 4.8 h for IAV and to 8.0 h for RSV, RV, and SARS-CoV-2. For the monoinfection model (Equations (1)-(4)), estimated parameters included the rates of viral infection (β), viral production (p), viral clearance (c), and infected cell clearance (δ). The rate of viral infection (β) was allowed to vary between 1 × 10 −9 and 1.0 RNA −1 d −1 or TCID −1 50 d −1 , and the rate of viral production (p) was allowed to vary between 1 × 10 −3 and 1 × 10 3 RNA/cell/d or TCID 50 /cell/d. The rate of infected cell clearance (δ) was given a lower limit of 1 × 10 −2 d −1 and an upper limit of 1 × 10 3 d −1 . For the coinfection models, we first simulated the monoinfection until the day of coinfection and then employed the coinfection models (Equations (5)- (8) or Equations (9)-(12)) while incorporating the enhancement or inhibition functions (Equations (13) and/or (14)). The strength of interaction (κ) was estimated for each scenario.
Fit quality was assessed using the Akaike information criterion with small sample size correction (AIC c ). The model with the lowest AIC c was considered the best, and ∆AIC c ≤ 2 was considered statistically equivalent [61].

Model-Predicted Mechanisms of RSV-IAV Coinfection
We used data from ferrets infected with a long strain of RSV and/or influenza A/Tasmania/2004/2009 (A[H1N1]pdm09; IAV) that had their viral loads measured until 14 days post infection (pi) [27]. When ferrets were inoculated with RSV followed by IAV 3 days later, morbidity was reduced, and IAV titers were slightly lower. To begin examining the interactions between RSV and IAV during coinfection that resulted in these dynamics, and establish the baseline parameter values for use in our mathematical models, we first fit the monoinfection model (Equations (1)-(4)) to the data from IAV-or RSVinfected ferrets (Table 1, Figure 2A). This showed a robust fit to each data set and yielded a faster rate of infection (β RSV = 1.0 × 10 −5 (RNA/100 µL) −1 d −1 vs. β IAV = 3.3 × 10 −6 (RNA/100 µL) −1 d −1 [IAV]) and slower rate of viral production for RSV (p RSV = 3.0 × 10 −2 RNA/100 µL/cell/d vs. p IAV = 2.5 × 10 1 RNA/100 µL/cell/d). Table 1. Best-fit parameters for the monoinfection model. Best-fit parameters obtained from fitting the monoinfection model (Equations (1)-(4)) to viral titers from ferrets intranasally infected with IAV at 3.5 log 10 TCID 50 or with RSV at 5.0 log 10 PFU [27]. Parameters are reported as the population median and the standard deviation of the associated random effect (ω). The initial numbers of target cells (T 0 ) and infected cells (E 0 ) were fixed to the indicated values, and the initial number of productively infected cells (I(0)) and the initial viral amount (V(0)) were set to 0.
Thus, to examine whether increases or decreases in the rates of viral infection, production, clearance, or infected cell clearance could better explain the data, we refit the models together with Equations (13) and/or (14). In total, we evaluated 16 scenarios for a single interaction and, on the basis of those results, up to 18 scenarios for dual interactions where each virus affected the other (Tables A1 and A2). When assuming that the two viruses competed for target cells (Equations (5)-(8)), single interactions that resulted in improved fits (i.e., lower AIC c ) included an RSV-induced reduction in the rate of IAV infected cell clearance (ffi − IAV ) or in the rate of IAV clearance (c − IAV ) ( Table 2; Figure A1A,B). Comparatively, some interactions within the target cell partitioning model (Equations (9)-(12)) led to a rebound of RSV ( Figure A2), which were excluded from consideration. The best suggested mechanisms under this hypothesis were either a decrease in the IAV infection rate (β − IAV ) by RSV together with a decrease in the RSV production rate (p − RSV ) by IAV (Table 2, Figure 1C) or an increase in the rate of RSV infected cell clearance by IAV (δ + RSV ; Table A2, Figure A1C). When allowing for dual interactions, the target cell competition model suggested that there were two sets of mechanisms that provided fits with similar AIC c values as the case when a single interaction was considered ( Table 2). Similar to the single interaction results, both sets of mechanisms included a RSV-induced reduction in the rate of IAV infected cell clearance (ffi − IAV ). This was paired with either an IAV-induced reduction in the rate of RSV production (p − RSV ; AIC c value of 175.7 (lowest); Figure 2B) or an IAV-induced increase in the rate of RSV clearance (c + RSV ; AIC c value of 175.8; Figure A1B). Allowing for dual interactions in the target cell partitioning model suggested an RSV-mediated reduction in the rate of IAV infectivity (fi − IAV ) coupled with an IAV-mediated reduction in rate of RSV production (p − RSV ; AIC c value of 139.3; Table 2, Figure 2C). The remaining single and double interactions for both models provided fits that were not statistically justifiable (Tables A1 and A2). Table 2. Best-fit parameters of the predicted mechanisms of RSV-IAV coinfection. Best-fit parameters from simulating the coinfection models with no interactions ('no interaction') or fitting the target cell competition model ('competition'; Equations (5)-(8)) or the target cell partitioning model ('partitioning'; Equations (9)-(12)) with the Equations (13) and/or (14) to viral loads from animals infected with RSV followed by IAV after 3 days. An NLME modeling approach was used and only the strength of interaction (κ) was estimated. Parameters are reported as the population median (κ) with standard deviation of the associated random effect (ω κ ). Fit quality is reported as log-likelihood (-2LL), AIC c , and standard deviation of the residual error (σ). The resulting relative change in total virus (∆V AUC ) is provided for each scenario.  (13) and (14)) to viral titers from ferrets infected with RSV followed by IAV after 3 days (IAV, white squares; RSV, white circles). Dynamics of the (B) target cell competition model (δ − IAV and p − RSV ) or (C) target cell partitioning model (β − IAV and p − RSV ) are shown along with the model schematic. Heatmaps are the relative change in total viral burden (i.e., ∆V AUC ; Equation (15)) evaluated for a range of interaction strengths (κ = 1 × 10 −9 to 1 × 10 2 (RNA/100 µL) −1 ) and infection intervals (0 to 11 days). The best-fit κ for a coinfection at 3 days is denoted by a white star.

Effect of Infection Timing and Interaction Strength in RSV-IAV Coinfection
Because the order, infection interval, and strength of interaction can influence coinfection dynamics and result in diverse disease phenotypes, we sought to better understand how these metrics alter RSV-IAV coinfection. To achieve this, we evaluated the relative change in total viral burden for each virus (i.e., ∆V AUC ; Equation (15)) across a wide range of infection intervals (0 to 11 days) and interaction strengths (κ; 1 × 10 −9 to 1 × 10 2 (RNA/100 µL) −1 ) in each model. Here, we focused on the best-fit models with the lowest log-likelihood (Table 2; i.e., δ − IAV and p − RSV in the target cell competition model and β − IAV and p − RSV in the target cell partitioning model). For the predicted interactions from the earlier analyses, differing the interaction strength (κ) had the strongest effect when the two infections were separated by shorter intervals (i.e., <3 days; Figure 2). When the interval was <3 days in the target cell competition model, a reduction in the RSV burden (up to ∆V AUC = −1) and prolonged IAV infection (up to ∆V AUC = 15) were observed for a higher interaction strength than the best fit δ − IAV value (white stars in Figure 2B). In the target cell partitioning model, the RSV burden was again reduced but the range of interaction strengths was narrower ( Figure 2C). Consistent with the data when the IAV infection was initiated 5 or 7 days after RSV [27], the simulations showed that intervals >3 to 4 days after RSV infection resulted in minimal changes in both models ( Figure 2B,C). However, only the target cell partitioning led to uninterrupted IAV infection or unchanged IAV viral burden for longer intervals between IAV and RSV infection ( Figure 2C).

IAV Coinfection with RV
Animals infected with IAV and RV at the same time, or RV two days before IAV, yielded lower weight loss and milder disease severity compared to animals infected with IAV alone (Figure 3). On the other hand, animals infected with IAV 2 days before RV underwent significantly higher weight loss that led to death of all animals by 7 days pi ( Figure 3; [33]). To assess the potential mechanisms during IAV coinfection with RV that could lead these empirical observations and provide potential translation to human infection, we first fit the monoinfection model (Equations (1)-(4)) to viral loads from human volunteers infected with IAV [49] or RV [50] (Table 3, Figure 3A). The model fit resulted in different infection kinetic rates for each virus where rate of IAV infection was lower (β IAV = 2.5 × 10 −6 (TCID 50 ) and the rate of IAV production was higher (p IAV = 3.0 × 10 −1 /TCID 50 /cell/day vs. p RV = 2.9 × 10 −3 /TCID 50 /cell/day (RV)).
We next used these parameters in the coinfection models (Equations (5)-(8) or Equations (9)-(12)) with or without interaction (Equation (13) or/and Equation (14)) to predict which mechanisms could lead to the distinct disease outcomes observed in the experimental study. Because viral loads were not measured, but weight loss in the infected animals was measured for IAV-RV and RV-IAV coinfections, we compared the estimated cumulative area under the curve (CAUC) of the infected cell dynamics to the weight loss [47], qualitatively matching the magnitude and timing of change. The CAUC of the combined infected cells dynamics of each virus from the coinfection models without any interactions could not recapitulate the weight-loss dynamics for any interval or order, confirming that interactions were occurring.

Simultaneous or Sequential RV-IAV Coinfection
When testing different interaction mechanisms that enhanced or inhibited one virus within the target cell competition model, the model predicted that the mechanism that could lead to the reduced disease severity observed in RV-IAV coinfection (simultaneous or separated by 2 d) was an IAV-mediated decrease in the rate of RV infected cell clearance (δ − RV ; Figure 3B,C). An intermediate signal strength was required for the simultaneous infection (κ δ − RV = 7 × 10 −5 /(TCID 50 /mL); Figure 3B) and a larger signal strength was required for a coinfection separated by 2 days (κ δ − RV = 1/(TCID 50 /mL); Figure 3C). In the two scenarios, this could produce similar reductions in the CAUC of the infected cells as the weight-loss patterns in coinfected animals ( Figure 3B,C). In addition, when IAV and RV were initiated simultaneously, the model-predicted kinetics showed significant reductions in IAV titers compared to IAV monoinfection ( Figure 3B). This was accompanied by a small increase in RV titers around the peak that was due to the IAV-mediated decrease in the rate of RV infected cell clearance. However, in contrast to the simultaneous infection, the model indicated that RV significantly reduced IAV titers when RV was initiated two days before IAV ( Figure 3C).
The reduced rate of RV-infected cell clearance (δ − RV ) was also identified by the target cell partitioning model ( Figure 4A,B). However, for a simultaneous infection, this needed to be coupled with an RV-mediated increase in the rate of IAV infected cell clearance (δ + IAV ; κ δ + IAV = 1 × 10 2 /(TCID 50 /mL); Figure 4A), an increase in the rate of IAV clearance (c − IAV ; κ c + IAV = 1 × 10 −3 /(TCID 50 /mL);), a decrease in the rate of IAV production (p − IAV ; κ p − IAV = 6 × 10 1 /(TCID 50 /mL); Figure A3B), or a decrease in the rate of IAV infectivity (β − IAV ; κ β − IAV = 6 × 10 1 /(TCID 50 /mL); Figure A3C) to achieve a consistent CAUC of the infected cells with the reduced weight loss. In all cases, the combined CAUC of the infected cells was reduced to a level below that of an IAV single infection and accompanied with a complete suppression of IAV titers without affecting RV titers.
For RV-IAV coinfection, no single interaction could reproduce the reduced disease severity. Thus, we did not consider dual interactions. However, because viral infections could initiate and/or modify host responses (e.g., Type I interferon, macrophages, and neutrophils) that were not included in our model, and this could translate into a reduced number of susceptible cells that is not automatically created by the target cell partitioning hypothesis, we examined the effect indirectly by reducing the initial number of target cells that were available for the second virus, as in prior studies [47,60]. For RV-IAV coinfection, reducing the initial number of target cells by 1 log 10 (i.e., T 0 = 2 × 10 7 cells for IAV compared to T 0 = 2 × 10 8 cells for RV; Table 3) was sufficient to reduce the combined CAUC of the infected cells compared to IAV monoinfection ( Figure 4B). However, the estimated CAUC of the infected cells deviated from the experimental results at later time points, where the model suggested similar but delayed IAV titers ( Figure 4B).

IAV-RV Coinfection
The mechanism that could lead to the increased disease severity observed in IAV-RV coinfection (separated by 2 days) within the target cell competition model was an RV-mediated decrease in the rate of IAV-infected cell clearance (δ − IAV ; Figure 3D). The required signal strength was large (κ δ − IAV = 3/(TCID 50 /mL)). Despite the significant model-predicted reduction in RV titers, the small increase in IAV titers was sufficient to create an increase in the combined CAUC of the infected cells that aligned with the increased weight loss observed in coinfected animals ( Figure 3D).
In contrast, the target cell partitioning hypothesis alone (i.e., no interactions) led to a higher combined CAUC of the infected cells and alignment with the observed increase in disease severity ( Figure 4C). This resulted in similar viral loads as those with the monoinfection for both viruses. However, including an IAV-mediated increase in the rate of RV production (κ p + RV = 1 × 10 −6 /(TCID 50 /mL); Figure 4C) or the rate of RV infectivity (κ β + RV = 1 × 10 −6 /(units/mL); Figure A3D) resulted in an earlier increase in the CAUC of infected cells, which matched the timing of the deviation in weight loss slightly better. In the two scenarios, the predicted IAV titer dynamics were similar and unchanged from a monoinfection, and the predicted RV titer dynamics had a similar shape, but were much higher when the production rate was increased (5.51 log 10 TCID 50 /mL vs. 4.63 log 10 TCID 50 /mL). Other possible mechanisms identified by this model included a reduction in the initial number of target cells (i.e., T 0 = 2 × 10 7 cells for RV compared to 2 × 10 8 cells for IAV; Table 3) coupled with either a reduction in the rate of IAV infected cell clearance by RV (c − IAV ; κ = 10/(TCID 50 /mL); Figure A3E) or in the rate of RV infected cell clearance by IAV (δ − RV ; κ = 10/(TCID 50 /mL); Figure A3F). In the first scenario, the model suggested that IAV titers would remain high for an extended period of time, which created an extended, flat viral-load peak. In the latter case, the model indicated that there would be no changes to IAV titers, and that RV had a slower increase with a lower peak.

IAV Coinfection with SARS-CoV-2
Animals infected with IAV followed 3 days later by SARS-CoV-2 resulted in increased weight loss and more severe disease severity compared to animals infected with IAV or SARS-CoV-2 alone (Figure 3) [37]. To examine the potential interactions between these two viruses, we employed the same approach as above. The results suggested similar mechanisms for enhanced disease severity for the two coinfection models, although the quantitative dynamics was distinct between the models. That is, the target cell competition model predicted a slightly higher combined CAUC of the infected cells when SARS-CoV-2 reduced the rate of IAV infected cell clearance (δ − IAV ) at the signal strength κ δ − IAV = 2 × 10 −2 /(RNA/mL) ( Figure 3E) whereas the target cell partitioning without any interaction led to a significantly higher combined CAUC of the infected cells ( Figure 4E). Reducing the initial number of target cells (i.e., a log 10 reduction [T 0 = 2 × 10 7 cells]) available to SARS-CoV-2 coupled with a reduction in the rate of SARS-CoV-2-infected cell clearance (δ − CoV2 ) by IAV (κ δ − CoV2 = 10/(RNA/mL)) or a reduction in the rate of IAV infected cell clearance (δ − IAV ) by SARS-CoV-2 (κ δ − IAV = 1 × 10 −2 /(RNA/mL)) led to a higher combined CAUC of infected cells. The predicted viral-load dynamics was distinct between the two models. In the target cell competition model, there was a significant reduction in SARS-CoV-2 titers and a delay in the resolution of IAV titers. In contrast, the target cell partitioning model suggested that the SARS-CoV-2 infection was simply delayed.

Discussion
Respiratory coinfections with multiple viruses are becoming more recognized clinically, particularly in light of the SARS-CoV-2 pandemic. Experimental studies have begun illuminating the outcome heterogeneity, which seems to rely on numerous factors such as the viral pairing, order, and timing of each infection, and specific immune factors. Although viral and immune dynamics during viral-viral coinfections is only beginning to be defined, mathematical models are useful to predict and narrow the spectrum of potential mechanisms, guide new experiments, and help in interpreting clinical, experimental, and epidemiological observations. Our analysis of different viral coinfection scenarios suggests that only a small subset of mechanisms could lead to the alterations in viral loads and/or disease severity observed in animal models (summarized in Table 4).
When a new virus invades within a few days before or after influenza, innate immune responses may be disrupted. The early regulation of Type I IFNs, macrophages, and/or neutrophils, among other innate immune responses, by a virus could impact the dynamics of subsequent infections. Although we did not assess these immune responses directly, we modeled them indirectly in various ways. The target cell competition model automatically assumes that fewer cells are available for infection by the second virus, which could emulate a protective mechanism that reduces the possible infection size for the coinfecting virus. In the target cell partitioning model, we decreased the number of target cells available for the second virus, which was used to mimic lower doses [47,60]. Both approaches assume that some cells are protected or otherwise unavailable for infection, which can be interpreted as the IFN-mediated protection of susceptible cells or an immediate clearance of the virus upon infection due to the activation of macrophages and/or neutrophils by the first virus. We observed the latter phenomenon in an experiment where CD8 T cells were depleted before infection [47]. In that case, innate responses were activated by the CD8 T cell death, which resulted in the partial clearance of the inoculum that emulated a reduced dose, and led to fewer infected cells and thus less virus. This automatically reduced inflammation and weight loss [47], which was similar to the reduced weight loss observed during RV-IAV coinfection [33]. Allowing for fewer susceptible cells within the target cell partitioning model was sufficient to explain the reduced weight loss in RV-IAV coinfected animals, potentially indicating a similar mechanism. While neither our model nor the data were sophisticated enough to specify the exact mechanism, higher IFN-β was detected at 2 days post coinfection, and IAV titers trended slightly lower compared to in the monoinfection [33]. A follow-up study suggested that the protection of IAV-mediated disease by RV was dependent on IFN and that this contributed to the control of neutrophilic inflammation [35]. These data align nicely with our predictions, which may help in connecting the underlying reasoning (i.e., fewer cells becoming infected) with downstream consequences (i.e., reduced immune activation and inflammation). Table 4. Summary of model-predicted mechanisms resulting in increased or decreased disease severity during viral-viral coinfection. Summary of model-predicted mechanisms that resulted in altered viral loads (RSV-IAV [27]) or disease severity as quantified by the CAUC of the infected cells [47,48] and measured by animal weight loss (IAV-RV and RV-IAV [33], and IAV-CoV2 [37]). IAV-infected cell clearance reduced by RSV and RSV production reduced by IAV Figure 2B IAV infected cell clearance reduced by RSV Figure A1A IAV clearance reduced by RSV Figure A1A IAV infected cell clearance reduced by RSV and RSV clearance increased by IAV Figure A1B IAV RV 0 Decreased [33] RV infected cell clearance reduced by IAV Figure 3B RV IAV 2 Decreased [33] RV infected cell clearance reduced by IAV Figure 3C IAV RV 2 Increased [33] IAV infected cell clearance reduced by RV Figure 3D IAV SARS-CoV-2 3 Increased [37] IAV infected cell clearance reduced by SARS-CoV-2 Figure 3E Partitioning RSV IAV 3

Model First Virus
Decreased [27] IAV infectivity reduced by RSV and RSV production reduced by IAV Figure 2C RSV-infected cell clearance increased by IAV Figure A1C IAV RV 0 Decreased [33] RV infected cell clearance reduced by IAV and IAV infected cell clearance increased by RV Figure 4A RV infected cell clearance reduced by IAV and IAV clearance reduced by RV Figure A3A RV infected cell clearance reduced by IAV and IAV production reduced by RV Figure A3B RV infected cell clearance reduced by IAV and IAV infectivity reduced by RV Figure A3C RV IAV 2 Decreased [33] Reduced number of target cells for IAV Figure 4B IAV RV 2 Increased [33] RV production increased by IAV Figure 4C No interaction Figure 4C RV infectivity increased by IAV Figure A3D Reduced number of target cells for RV and IAV infected cell clearance reduced by RV Figure A3E Reduced number of target cells for RV and RV infected cell clearance reduced by IAV Figure A3F IAV SARS-CoV-2 3 Increased [37] SARS-CoV-2 production increased by IAV Figure 4D No interaction Figure 4D SARS-CoV-2 infectivity increased by IAV Figure A4A Reduced number of target cells for SARS-CoV-2 and IAV infected cell clearance reduced by SARS-CoV-2 Figure A4B Reduced number of target cells for SARS-CoV-2 and SARS-CoV-2-infected cell clearance reduced by IAV Figure A4C In addition to potential effects on early host responses, one of the most common mechanisms defined by our analysis was altered rates of infected cell clearance, which may indicate an effect on virus-specific CD8 T-cell responses [47]. In some cases, there was a negative effect, while in others, there was a positive effect (Table 4). Further, there were some scenarios where each virus affected the other virus' infected cell clearance rate. This may indicate variation in the number, composition, and/or function of epitope-specific T cells, which was observed for other viral pairs (e.g., lymphocytic choriomeningitis (LCMV) and Pichinde (PICV) viruses) [62]. For scenarios where the rates are reduced, it could also indicate that having fewer available target cells from preinfection would have automatically reduced the number of T cells needed to clear the infection. Data defining host responses are important to in helping to distinguish between these possibilities, particularly because CD8 T cells have to be significantly reduced to have a robust impact on viral loads, and their efficacy is dependent on the number of infected cells [47].
Effects on other infection processes such as the rates of viral infectivity or production were detected for some coinfection scenarios, although this was rarely the sole mechanism (Table 4). Only in IAV-RV and IAV-CoV2 coinfections within the target cell partitioning model were these potential single interactions, and they each resulted in similar viral load dynamics ( Figures 4C,D, A3D and A4A). This is because the type of model used here cannot typically distinguish between the effects of these two processes [63]. There is some evidence that the infectivity of SARS-CoV-2 is increased by IAV but not RSV within cell cultures [41,43]. The underlying mechanism driving this remains unknown, but other studies with IAV and PIV have shown enhanced infection rates where PIV increased cell-to-cell fusion and, thus, spread of IAV [32]. However, another study found that IAV titers were increased while SARS-CoV-2 titers were decreased during IAV-CoV2 coinfection in mice [42]. The discrepancy between these studies, in which one found there was an increase in the number of infected cells, while the other found reduced viral loads, may have been due to the results being obtained in vitro vs. in vivo or to another interaction (e.g., IFN suppression of SARS-CoV-2 or clearance by macrophages). Our analyses suggest that SARS-CoV-2 titers would be reduced within the target cell competition model ( Figure A3E,F) or when the number of susceptible cells was reduced within the target cell partitioning model ( Figure A4B,C). This may indicate a role for the IAV-activated innate immune response and/or a lower effective dose of SARS-CoV-2.
A limitation of this work is that we used minimal data that were of different types (i.e., viral loads or weight loss) to predict potential mechanisms during viral coinfection. The quantity and type of data used are important because altered viral loads (e.g., as in RSV-IAV coinfection), immune cells, or cytokines do not always directly equate to differences in disease severity. We previously showed this phenomenon where we found that the number of infected cells and inflammation were nonlinearly correlated to disease severity [47]. This was particularly evident in one experiment where CD8 T cells were depleted, which resulted in only small reductions in viral loads, yet large reductions in weight loss. However, as mentioned above, the early immune activation led to a predicted lower effective dose (i.e., fewer cells becoming infected) and, thus, reduced disease severity. Similarly, in RSV-IAV coinfection, IAV titers were slightly increased, yet less weight loss occurred [27]. Thus, the higher viral loads later in infection may be insignificant with respect to severity. This highlights that, while some mechanisms may occur and alter viral loads, they could be distinct from those that yield disease outcomes. Our results from matching the qualitative differences in weight-loss data, which is a measure of disease severity, for IAV coinfection with RV or SARS-CoV-2 may better represent potential mechanisms with measurable differences in outcome. However, many of these also led to predicted differences in viral loads. Some information about mechanism may be able to be deduced from the timing of when weight loss begins to deviate from the monoinfection. In several scenarios, this occurred directly after the initiation of the secondary infection, which suggests that the environment created by the first virus had immediate effects.
The mechanisms suggested by our analysis occasionally differed depending on the underlying model hypothesis (i.e., whether viruses competed for epithelial cells) and, in some cases, resulted in different predicted viral load kinetics. Because most respiratory viruses can infect various types of airway epithelial cells, and replicate in the upper and lower respiratory tracts, it is conceivable that each virus would have ample cells to infect. However, by chance or due to the airway structure, they may enter the same region and interact on a local level. This may lead to cells being coinfected with both viruses, which we did not model explicitly. We indirectly modeled the potential effects of coinfected cells by assuming that the rates of infection and/or production could be different. Interestingly, cellular coinfection was detected during simultaneous infection with RSV and HMPV, where coinfected cells were possible but in small numbers and less likely in the presence of IFN [36]. The same may be true during other coinfections with viruses that are sensitive to IFN antiviral responses. To model the impact of coinfected cells and potential heterogeneity in their prevalence, agent-based models may be better suited than those used here. However, it is important to determine whether and for how long these local effects have a measurable impact during an in vivo infection.
Examining data from viral-viral coinfections using mathematical models allowed for us to reduce the number of possible underlying mechanisms that could result in altered viral load kinetics and/or disease severity. Although the models were relatively simple and lacked investigation into specific host immune responses, the analysis provides the infrastructure to integrate immunological models of higher complexity once data becomes available. Models for some immune responses during respiratory virus infections are already being developed and validated with experimental data (reviewed in [45,64]). Some of the insight from those studies was integrated here and aided our ability to interpret the small amount of experimental data currently available. However, establishing better methods that could predict disease severity (e.g., as in [47]) is critical. Our ability to assess the contribution and timescales of different mechanisms to infection kinetics and outcome should increase as more temporal viral load, immunologic, and pathologic data become available. In addition, the generated hypotheses should aid in experimental design, ultimately leading to a more complete understanding of respiratory virus coinfection. Funding: This research was funded by NIH grants AI125324 and AI139088 (AMS), and P20GM104420 (CHR).

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results. The complete set of best-fit parameters obtained from fitting the coinfection models ('competition', Equations (5)-(8); 'partitioning', Equations (9)-(12)) with the Equations (13) and/or (14) to viral loads from animals infected with RSV followed by IAV after 3 days is in Tables A1 and A2. The top models are highlighted in gray and the dynamics are in Figures 2 and A1. Models with a predicted RSV rebound ( Figure A2) were excluded. Table A1. Full parameter results from fitting the target cell competition model to RSV-IAV coinfection in ferrets. Parameters from simulating ('no interaction') or fitting the target cell competition model ('competition'; Equations (5)-(8)) with the Equations (13) and/or (14) to viral loads from animals infected with RSV followed by IAV after 3 days. An NLME modeling approach was used and only the strength of interaction (κ) was estimated. Parameters are reported as the population median (κ) with standard deviation of the associated random effect (ω κ ). Fit quality is reported as log-likelihood (-2LL), AIC c , and standard deviation of the residual error (σ). The best models and predicted mechanisms are highlighted in gray and summarized in Table 1.    (8)) with the Equations (13) and/or (14) to viral loads from animals infected with RSV followed by IAV after 3 days. An NLME modeling approach was used and only the strength of interaction (κ) was estimated. Parameters are reported as the population median (κ) with standard deviation of the associated random effect (ω κ ). Fit quality is reported as log-likelihood (-2LL), AIC c , and standard deviation of the residual error (σ). The best models and predicted mechanisms are highlighted in gray and summarized in Table 1. Although some models resulted in the lower AIC c values, they were excluded because they produced a viral load rebound.   4)) and fit of the coinfection models (solid lines; Equations (5)- (8) or Equations (9)-(12)) with the interaction functions (Equations (13) and (14)) to viral titers from ferrets infected with RSV followed by IAV after 3 days (RSV, white circles; IAV, white squares). Heatmaps are the relative change in total viral burden (i.e., ∆V AUC ; Equation (15)) evaluated for a range of interaction strengths (κ = 1 × 10 −9 to 1 × 10 2 RNA/100 µL) and infection intervals (0 to 11 d). The best-fit κ for a coinfection at 3 days is denoted by a white star.   4)) or the target cell partitioning model with the indicated mechanism (solid lines; Equations (9)- (12)) to viral titers from ferrets infected with RSV (circles) followed by IAV (squares) after 3 days. Each mechanism included a reduction in the rate of RSV infectivity (β − RSV ), which was paired with either decreased rates of RSV infected cell clearance (δ − RSV ), IAV infectivity (β − IAV ), or IAV clearance (c − IAV ). Because the RSV dynamics rebounded, these were excluded as possible mechanisms.  [37]. Dynamics of IAV-CoV2 coinfection with (A) SARS-CoV-2 infectivity increased by IAV (β + CoV2 ), (B) reduced number of target cells for SARS-CoV-2 (T − 0 ) and IAV infected cell clearance reduced by SARS-CoV-2 (δ − IAV ), or (C) reduced number of target cells for SARS-CoV-2 (T − 0 ) and SARS-CoV-2-infected cell clearance reduced by IAV (δ − CoV2 ).