Mathematical Modeling of RNA Virus Sensing Pathways Reveals Paracrine Signaling as the Primary Factor Regulating Excessive Cytokine Production

: RNA viruses, such as inﬂuenza and Severe Acute Respiratory Syndrome (SARS), invoke excessive immune responses; however, the kinetics that regulate inﬂammatory responses within infected cells remain unresolved. Here, we develop a mathematical model of the RNA virus sensing pathways, to determine the intracellular events that primarily regulate interferon, an important protein for the activation and management of inﬂammation. Within the ordinary di ﬀ erential equation (ODE) model, we incorporate viral replication, cell death, interferon stimulated genes’ antagonistic e ﬀ ects on viral replication, and virus sensor protein (TLR and RIG-I) kinetics. The model is parameterized to inﬂuenza infection data using Markov chain Monte Carlo and then validated against infection data from an NS1 knockout strain of inﬂuenza, demonstrating that RIG-I antagonism signiﬁcantly alters cytokine signaling trajectory. Global sensitivity analysis suggests that paracrine signaling is responsible for the majority of cytokine production, suggesting that rapid cytokine production may be best managed by inﬂuencing extracellular cytokine levels. As most of the model kinetics are host cell speciﬁc and not virus speciﬁc, the model presented provides an important step to modeling the intracellular immune dynamics of many RNA viruses, including the viruses responsible for SARS, Middle East Respiratory Syndrome (MERS), and Coronavirus Disease (COVID-19).


Introduction
Respiratory infections are a constant threat to public health with deadly infections often characterized by cytokine storms, i.e., overly aggressive innate immune responses that result in severe, unnecessary lung tissue inflammation. Influenza causes an annual 3-5 million infections and 290,000-650,000 deaths [1]. The disease presents a significant pediatric burden, causing 374,000 annual hospitalizations of children under 1 year old, globally [2]. Currently, the new Severe Acute Respiratory Syndrome of the novel Coronavirus (SARS-COV2), responsible for the disease COVID-19, has emerged, resulting in over 113,000 deaths in the US to date. Both influenza and SARS-COV2 are RNA viruses, and the data to date suggests that cytokine storms are a common feature of both viruses during severe infections [3]. Immune responses can help or hinder an organism's ability to overcome an infection, and excessively, inflammatory responses, like cytokine storms, can cause greater tissue damage, higher mortality, and slow recovery [4,5]. Vaccination is effective for protecting public health against seasonal influenza; however, when new strains unexpectedly emerge, such as the 2009 novel pandemic H1N1 virus or the 2019 SARS-COV2 virus, new treatment strategies that can be implemented rapidly, and preferably independently of the specific virus, are needed. Immunomodulatory treatments that aim to reduce inflammation while still managing virus growth are a promising approach to protecting against emergent disease, but several fundamental questions on how unnecessarily aggressive immune responses emerge remain unknown. Mathematical modeling can help quantify the kinetics of the interactions that define the immune system, revealing the interactions that are most likely to be responsible for unnecessarily aggressive responses and potential targets, to interfere with immunity to ensure healthy virus clearance.
The first step to occur during an immune response is the detection of the pathogens, leading to the early, localized, innate immune response [6]. This response to viral infection leads to the production of type I interferons (IFNs). Interferons serve to establish an antiviral state by activating and inducing Mx proteins, RNA-activated protein kinase, and the 2-5A system [7]; they also regulate other immune responses, by acting on natural killer cells, T cells, B cells, dendritic cells, and phagocytic cells [8]. The presence of influenza virus is primarily sensed by cytoplasmic retinoic acid-inducible gene 1 (RIG-I) and endosomal Toll-like receptors 7 and 9 (TLR) [9,10]. RIG-I senses viral RNA in the cytoplasm [11], but is antagonized by many influenza A viruses' nonstructural protein I (NS1) to varying, strain-specific magnitudes [12,13]. TLR7 is free of this antagonism [14,15], and is activated after the influenza envelope has been degraded by endosomal proteases. SARS-CoV's N protein has been implicated in the inhibition of Type-I interferon production via antagonism of RIG-I [16,17]. This suggests RIG-I as a common viral sensor protein and a common target of antagonism for RNA viruses. The activation of either sensor leads to the phosphorylation of interferon regulatory factor 7 (IRF7, IRF7P) and the production of interferons (IFN), to act as a signaling cytokine. IFN induces secondary messenger molecules, ultimately leading to the induction of immune modulation and antiviral genes [18]. IFN is secreted from the infected cell and sensed through the Janus kinase/signal transducer and activator of transcription pathway (JAK/STAT), in both an autocrine and paracrine manner. The JAK/STAT pathway leads to the induction of a broad family of IFN-stimulated antiviral genes [19], as well as the supplementary (autocrine) or novel (paracrine) production of IFN. These IFN-stimulated genes cause cell death through apoptosis, necroptosis and pyroptosis [20], slow viral replication within the cell [21], regulate the infiltration and activity of key innate immune cells to clear the infection, and help initiate the adaptive immune response [22].
The current models of innate immune response to RNA virus infection lack major intracellular components or lack important biological interactions, limiting their applicability to understanding how severe inflammation emerges. Extensive molecular pathway maps exist [6,23], but they currently lack mathematical description to support simulating the immune response. Some models of the intracellular innate immune response have incorporated RIG-I and TLR activity, but consider their effect to be constant, independent of the viral load, and the models are inherently unstable, complicating their use [24,25]. These models used ordinary differential equations (ODEs), a common approach in systems biology, after their demonstrable success in analyzing the robustness of biological signaling [26][27][28], the highly dynamic behaviors of NF-kB [29], and ultrasensitive cell fate binary responses [30]. ODEs allow for interpolation of the dynamics between a finite number of time points at which data has been measured, based on hypotheses of the mechanisms regulating the system's components.
In this study, we construct a novel ODE model to simulate the intracellular innate immune response of human bronchial epithelial cells (HBECs) to influenza A infection and use the model to determine the interactions that most affect cytokine production. This model was constructed with computational expense for parameterization and eventual agent based model implementation in mind, with a minimum number of ODE's and parameters that capture the dynamics of interest. The model is numerically stable under realistic conditions and non-stiff, enhancing its reproducibility and reducing computational cost. The model incorporates a viral growth model [31] and the proportionality of sensor protein activity to vRNA levels in the cytoplasm, the first such integration of cell dynamics and viral replication. The feedback of interferon production on viral replication through the interferon stimulated gene (ISG) family [32] is included. A literature search for data to perform parameterization produced viral titers [33,34] and the time-series of RNA data in HBECs [35]. RNA data originated from Shapira et al.'s 2009 work elucidating a network of viral-host interactions via genome wide expression profiling. Viral titers came from Ramos et al.'s 2013 work on the polyadenylation stimulating factor 30 (CPSF30) binding function of the NS1 protein. These consist of subsets, in which competing, parallel pathways were inhibited, allowing for improved identifiability of the model parameters; first, a wild-type A/Puerto Rico/8/1934 Influenza A (PR8) infection, in which RIG-I is assumed fully antagonized and TLR is fully active; and second, an NS1 knockout PR8 strain which has both TLR and RIG-I activity [35]. The antagonism of RIG-I in wild-type PR8 infection is shown to drastically alter infection outcomes. Additionally, some parameters were sourced from or bounded by their respective values in previous models. This work establishes the first cell-level model of interferon signaling induced by influenza infection that can be used to compare host responses between infections with different influenza viruses, antagonism motifs and different RNA viruses. Paracrine signaling is demonstrated to produce the majority of HBEC's cytokine response to influenza infection, while the initial sensor protein pathways are shown to serve as an ignition for said paracrine signaling.

Model Construction
Upon initial viral entry to the cell, there is an eclipse phase, during which the virus enters the cell, traffics to the nucleus, and starts the replication process [36]. Analogously, a whole-organism infection with influenza A observes an eclipse phase before significant viral titers and immune response are seen [32,37]. As the virus replicates within the cell, the concentrations of single-stranded viral RNA (vRNA) and NS1 in the cytoplasm rise until cell death. Intuitively, the magnitude of action of the sensor protein TLR should increase proportionally to the vRNA levels. RIG-I's activity level will also be proportional to the vRNA level, but inversely proportional to the NS1 concentration.
A model of IFN production, virus sensing, and JAK/STAT feedback was built to reflect these mechanisms. The initial model was unstable; such models fail to capture important aspects of the dynamics of shutdown and steady-state transitions. The complete model had degradation rates for environmental type-I interferons (IFNe), interferon regulatory factor 7 and its phosphorylated form (IRF7, IRF7P) added to assist in stability. For each modeled interaction, mass action was assumed if no other information was available on the reaction kinetics, unless an existing literature model had specific evidence for the representation of the step in some other form, in order to maintain a minimum number of parameters. The IFN Hill type kinetic was selected based on in vivo data from mice, which showed that IFN has an ultrasensitive response to viral load. Equations (7) and (8) were added to incorporate the classical model of virus kinetics [31,32], to model virus concentration within the cell, infected cell counts, and type-I interferons' effects on viral replication. The use of viral counts, [V], in Equation (2), makes the sensor proteins' actions proportional to virus level. RIG-I's action is modeled as a mass-action kinetic, while TLR is modeled separately, as a Hill-type kinetic. All Type-I Interferons are considered indistinguishable in the complete model, as are the different STAT species. The transport of interferons from the cytoplasm to the extracellular space is now represented with a first-order mass transfer kinetic. Qiao et al. [24] implemented Michaelis-Menten (MM) for this transport step; this was rejected due to a lack of evidence for any specific MM-type mechanism, and the redundancy of limiting maximum IFN concentrations with the introduction of degradation rates. Finally, the binding of interferons and the activation of the JAK/STAT pathway [24,25] was reduced to constants in a form reminiscent of a Hill kinetic. The production of cellular species were made proportional to the living cell population, [P].

Data Sources
Three primary literature sources were used for data to estimate model parameters. First, micro array gene expression data [35] of two influenza strain (PR8 and an NS1-knockout PR8) time-course experiments in human bronchial epithelial cells (HBECs) were used to fit IFN, STATP, and IRF7 gene expression. Second, viral titers of wild-type PR8 influenza in human lung adenocarcinoma epithelial (A549) cells [34] and NS1-knockout PR8 [33] were used to fit viral titers. Viral titers would ideally be obtained with the same cell type, time points, and infection methodology as the micro-array data; however, the immune response similarity of A549 cells and HBECs [38] justifies this approach.

ODE Simulation and Sensitivity Analysis
Julia v1.3 was used to simulate the ODE model with the DifferentialEquations v6.11.0, ParameterizedFunctions v4.2.1, and DiffEqParamEstim v1.12.0 packages. The ODE system is solved with the non-stiff solver VERN7. A Sobol method global sensitivity analysis was conducted to determine the degree of control that each parameter exerted on the system. Julia's DiffEqSensitivity v6.7.0 was used for this analysis.

Model Parameterization
Since the ODE system relies on 15 unknown parameters, simple regression methods are insufficient to successfully parameterize the model. Instead, a parallel temping Markov chain Monte Carlo method (PT MCMC) was implemented in Julia v1.3. To initialize the parameters, a literature search and manual fitting methods were conducted. The literature search provided estimated decay rates for STATP [39], IRF7 [40], and IRF7P [41]. A manual fitting gave estimates and stability-based bounds for cell death (k 61 ), viral replication (k 71 ), and nonspecific viral clearance (k 73 ) rates. Estimates from the Qiao model [24] were used to initialize the remaining parameters. Since parallel tempering results in faster convergence than single-chain methodologies [42], PT MCMC was run with 1 million iterations with three parallel chains, for a total of three million samples per fitting attempt. The sum squared error minimized by the MCMC was: The left-hand portion of the objective function determines the error of the system dynamics, where O ij is the experimentally observed log-fold change of Species i at Time j [35], relative to a control RNA level of the same species and time point. Two points from biological replicates are available for each i,j pair. Intracellular interferon and IRF7 were directly tracked by their RNA levels. Because phosphorylated STAT cannot be measured using microarray, the RNA levels of Interferon-induced GTP-binding Protein Mx1 (MX1), which is induced primarily by the action of STATP [6], were used as a proxy for STATP levels. Next, n is a weighting factor. E k is the normalized literature viral titer estimate at Time k [33], and V k is the normalized calculated viral titer at Time k.

Structural Identifiability
The structural identifiability of an ODE system is a prerequisite to parameter estimation and the subsequent application of the model. A common problem in ODE representations of biological systems is a lack of identifiable parameters, leading to non-unique parameter sets [43]. Structural identifiability analyses were carried out with structural identifiability taken as extended-generalized observability with lie derivatives and decomposition [44] (STRIKE-GOLLD) in MATLAB R2019a. These analyses were done under two sets of conditions-perfect identifiability and practical identifiability. Under perfect identifiability conditions, all seven species are assumed to be perfectly observed (i.e., measured directly by experiment). Under practical identifiability, only IFN, STATP, IRF7, Cells, and Virus were considered to be observable, which reflects the availability of data under which the model was trained. Structural identifiability results are available in Appendix A.

Interparameter Correlation
Parameter correlation in MCMC training results was tested using Pearson's correlation coefficient from SciPy v1.4.1 in Python 3.6.8. Significant correlations were considered as those with a correlation coefficient > ±0.5 [45].

ODE Model
The model consists of seven ordinary differential equations (states) with three fixed parameters and 15 unknown parameters. The complete model is given in Equations (2)-(8) and illustrated graphically in Figure 1.

ODE Model
The model consists of seven ordinary differential equations (states) with three fixed parameters and 15 unknown parameters. The complete model is given in Equations (2)-(8) and illustrated graphically in Figure 1.
Equations (2) [46], it is assumed that the entire cell population is producing new viral particles at the start of the trial.
[P] will thus vary from 1 to 0 and is unitless. State 7, [V], represents a virus concentration normalized to the maximum Equations (2) [35] results in 99.3% of target cells becoming infected [46], it is assumed that the entire cell population is producing new viral particles at the start of the trial.
[P] will thus vary from 1 to 0 and is unitless. State 7, [V], represents a virus concentration normalized to the maximum amount observed; [V] = (virus concentration)/(max virus concentration). These values can take on molarity, PFU, or similar matching units; [V] is unitless and can vary from 0 to 1. The viral titers [33,34] are in units of PFU mL −1 , and the maximum measure was the 24-h observation of wild-type PR8. A [V] value > 1 is possible for viral strains with higher peak viral loads than wild-type PR8 and at model timespans greater than 24 h.

MCMC Parameterization
During parameterization, the system was simulated out to 48 h post infection (HPI); only the first 36 h are shown here for clarity. All accepted parameter fits require the system to be stable, returning to zero after the infection has run its course. In this system, the multiplicity of infection is 5 [35], leading to 99.3% of all cells being initially infected [46]. Shapira et al.'s [35] experiments were carried out in vitro, without any immune cell presence. Thus, the model does not incorporate immune cell presence or phagocytosis, and it was assumed that all cells will die solely because of viral effects. The steady state for the model is complete cell culture death and the eventual decay of all species. The parameter set with the lowest SSE is shown in Figure 2.
Processes 2020, 8, x FOR PEER REVIEW 6 of 16 [33,34] are in units of PFU mL −1 , and the maximum measure was the 24-h observation of wild-type PR8. A [V] value > 1 is possible for viral strains with higher peak viral loads than wild-type PR8 and at model timespans greater than 24 h.

MCMC Parameterization
During parameterization, the system was simulated out to 48 h post infection (HPI); only the first 36 h are shown here for clarity. All accepted parameter fits require the system to be stable, returning to zero after the infection has run its course. In this system, the multiplicity of infection is 5 [35], leading to 99.3% of all cells being initially infected [46]. Shapira et al.'s [35] experiments were carried out in vitro, without any immune cell presence. Thus, the model does not incorporate immune cell presence or phagocytosis, and it was assumed that all cells will die solely because of viral effects. The steady state for the model is complete cell culture death and the eventual decay of all species. The parameter set with the lowest SSE is shown in Figure 2. The parameterization of the ODE model was accomplished with a parallel tempering Markov chain Monte Carlo (MCMC). The first 10 3 iterations comprised burn-in, a period where the MCMC algorithm was searching parameter space for potential local minima. This can be seen in the acceptance ratio ( Figure 3A) and sum squared error ( Figure 3B). The parameterization of the ODE model was accomplished with a parallel tempering Markov chain Monte Carlo (MCMC). The first 10 3 iterations comprised burn-in, a period where the MCMC algorithm was searching parameter space for potential local minima. This can be seen in the acceptance ratio ( Figure 3A) and sum squared error ( Figure 3B). The MCMC algorithm would ideally accept new parameters 23% of the time [47]; for this parameterization, the acceptance ratio of the primary chain was 19%. This was considered acceptable. Further hyperparameter tuning to obtain exactly 23% acceptance was possible, but not pursued. The overall fit of the model, as quantified by the sum squared error (SSE) (Section 2.4), is shown in Figure  3B. After burn-in, no trend towards fit improvement is observed, despite high-temperature chain exploration for other local minima, thus, the MCMC fitting algorithm is considered to be sufficiently converged. Two additional optimizations of the same length were run with randomized starting parameter values, which converged to parameterizations with the same SSE value and ranges as Run 1. The parameterization with the lowest SSE was used. The values for this parameterization are given in Table 1, along with their origin.  The MCMC algorithm would ideally accept new parameters 23% of the time [47]; for this parameterization, the acceptance ratio of the primary chain was 19%. This was considered acceptable. Further hyperparameter tuning to obtain exactly 23% acceptance was possible, but not pursued. The overall fit of the model, as quantified by the sum squared error (SSE) (Section 2.4), is shown in Figure 3B. After burn-in, no trend towards fit improvement is observed, despite high-temperature chain exploration for other local minima, thus, the MCMC fitting algorithm is considered to be sufficiently converged. Two additional optimizations of the same length were run with randomized starting parameter values, which converged to parameterizations with the same SSE value and ranges as Run 1. The parameterization with the lowest SSE was used. The values for this parameterization are given in Table 1, along with their origin. The parameter distributions are shown in Figure 4. The lack of normalcy was a possible indication of interparameter correlation. This correlation was tested using the Pearson correlation coefficient. Nine significant correlation pairs were identified, and are shown in Figure 5. The magnitude of correlations for all parameter combinations is summarized in Appendix B, Table A1. The parameter distributions are shown in Figure 4. The lack of normalcy was a possible indication of interparameter correlation. This correlation was tested using the Pearson correlation coefficient. Nine significant correlation pairs were identified, and are shown in Figure 5. The magnitude of correlations for all parameter combinations is summarized in Appendix B, Table A1.    The parameter distributions are shown in Figure 4. The lack of normalcy was a possible indication of interparameter correlation. This correlation was tested using the Pearson correlation coefficient. Nine significant correlation pairs were identified, and are shown in Figure 5. The magnitude of correlations for all parameter combinations is summarized in Appendix B, Table A1.   High correlation and non-normal distributions were not unexpected, and are a frequent challenge in fitting nonlinear, biological signaling systems. Several model features, namely multiple parameters affecting the same species and the presence of feedback loops by nature, necessitate a correlated random walk to maintain or improve the model's fit. This manifests as a lower than expected acceptance ratio, which can be overcome by hyperparameter tuning if undesirably slow parameter space exploration results.

Model Validation by Predicting Response to Infection Using a NSI Knockout Influenza Virus
Once the model was trained, a validation study on a nonstructural protein 1 (NS1) knockout strain of PR8 influenza (dNS1PR8) was conducted. Since NS1 was assumed to be fully antagonizing the action of the sensor protein RIG-I during the initial training, the kinetic parameter (k 11 ) associated with RIG-I was unfit and assumed to be zero (unitless). It was manually estimated at 10 5 (unitless) for the dNS1PR8 simulation, based on the same SSE function as the MCMC parameterization. All other parameters maintained their values from the wild-type PR8 training. An ensemble of the top 1000 model parameterizations from the wild-type training were simulated with the nonzero k 11 term and plotted against the NS1 knockout data [35]. The validation model simulation results are shown in Figure 6. High correlation and non-normal distributions were not unexpected, and are a frequent challenge in fitting nonlinear, biological signaling systems. Several model features, namely multiple parameters affecting the same species and the presence of feedback loops by nature, necessitate a correlated random walk to maintain or improve the model's fit. This manifests as a lower than expected acceptance ratio, which can be overcome by hyperparameter tuning if undesirably slow parameter space exploration results.

Model Validation by Predicting Response to Infection Using a NSI Knockout Influenza Virus
Once the model was trained, a validation study on a nonstructural protein 1 (NS1) knockout strain of PR8 influenza (dNS1PR8) was conducted. Since NS1 was assumed to be fully antagonizing the action of the sensor protein RIG-I during the initial training, the kinetic parameter (k11) associated with RIG-I was unfit and assumed to be zero (unitless). It was manually estimated at 10 5 (unitless) for the dNS1PR8 simulation, based on the same SSE function as the MCMC parameterization. All other parameters maintained their values from the wild-type PR8 training. An ensemble of the top 1000 model parameterizations from the wild-type training were simulated with the nonzero k11 term and plotted against the NS1 knockout data [35]. The validation model simulation results are shown in Figure 6. The validation case lends credence to the underlying model structure; a different combination of sensor proteins and a modified strain of influenza can be modeled only by introducing a new term for the previously antagonized RIG-I. The validation simulation shows several key differences. First, the interferon and extracellular interferon peaks are larger in magnitude and occur faster than in the wild type, in good agreeance with microarray data. Second, the viral load is 96% lower than the wild type, which agrees with the viral titer data [33]. Finally, the simulation predicts that only 20% of The validation case lends credence to the underlying model structure; a different combination of sensor proteins and a modified strain of influenza can be modeled only by introducing a new term for the previously antagonized RIG-I. The validation simulation shows several key differences. First, the interferon and extracellular interferon peaks are larger in magnitude and occur faster than in the wild type, in good agreeance with microarray data. Second, the viral load is 96% lower than the wild type, which agrees with the viral titer data [33]. Finally, the simulation predicts that only 20% of infected epithelial cells die, regardless of the simulation's time frame. While this qualitatively agrees with the much lower lethality of the dNS1PR8 strain, it does not reflect the biological expectation of certain cell death after the infection of all epithelial cells, regardless of viral strain. Overall, the validation study suggests that the model structure is sound and is capable of novel predictions.

Sensitivity Analysis Reveals IRF7 Phosphorylation as Critical Step
A Sobol global sensitivity analysis was conducted, allowing each parameter to vary over the same parameter space that the MCMC algorithm explored in the initial training. The resulting system sensitivities are shown in Figure 7.
infected epithelial cells die, regardless of the simulation's time frame. While this qualitatively agrees with the much lower lethality of the dNS1PR8 strain, it does not reflect the biological expectation of certain cell death after the infection of all epithelial cells, regardless of viral strain. Overall, the validation study suggests that the model structure is sound and is capable of novel predictions.

Sensitivity Analysis Reveals IRF7 Phosphorylation as Critical Step
A Sobol global sensitivity analysis was conducted, allowing each parameter to vary over the same parameter space that the MCMC algorithm explored in the initial training. The resulting system sensitivities are shown in Figure 7. The most sensitive parameters are those which exert the most control over the innate immune response. Notably, the initial sensing of the virus' presence via TLR (k12) contributes a relatively small proportion of the overall system response; most of the immune response originates from the paracrine signaling pathway. Parameters k42 and k51 correlate to the phosphorylation of IRF7 to IRF7P, and the induction of additional IRF7 by the action of IRF7P, respectively, which dominate the paracrine signaling pathway. Thus, the sensitivity analysis suggests that the model's outcomes are strongly controlled by these two interactions.

Simulating Varying Levels of RIG-I Antagonism Reveals Robust Sensor Protein Action
Next, an in silico knockdown study of RIG-I was performed. This study is comparable to varying the production and effectiveness of the NS1 protein across several influenza virus strains [12,13]. This was done by varying the k11 parameter in 25% increments from dNS1-PR8's lack of antagonism (0% RIG-I knockdown, k11 = 10 5 ) to wild-type PR8's total antagonism (100% RIG-I knockdown, k11 = 0). As shown in Figure 8, any RIG-I activity above zero (0% to 75% knockdown) showed significant reductions in viral load and target cell lethality. 0% knockdown yielded the largest magnitude immune response, as quantified by extracellular IFN, however, 25% through 75% knockdowns showed a robust immune response of near equal magnitude despite the reduction in RIG-I activity. This suggests that RIG-I is robust against viral antagonism and plays a vital role in initializing the host's cytokine response to viral infections. The most sensitive parameters are those which exert the most control over the innate immune response. Notably, the initial sensing of the virus' presence via TLR (k 12 ) contributes a relatively small proportion of the overall system response; most of the immune response originates from the paracrine signaling pathway. Parameters k 42 and k 51 correlate to the phosphorylation of IRF7 to IRF7P, and the induction of additional IRF7 by the action of IRF7P, respectively, which dominate the paracrine signaling pathway. Thus, the sensitivity analysis suggests that the model's outcomes are strongly controlled by these two interactions.

Simulating Varying Levels of RIG-I Antagonism Reveals Robust Sensor Protein Action
Next, an in silico knockdown study of RIG-I was performed. This study is comparable to varying the production and effectiveness of the NS1 protein across several influenza virus strains [12,13]. This was done by varying the k 11 parameter in 25% increments from dNS1-PR8's lack of antagonism (0% RIG-I knockdown, k 11 = 10 5 ) to wild-type PR8's total antagonism (100% RIG-I knockdown, k 11 = 0). As shown in Figure 8, any RIG-I activity above zero (0% to 75% knockdown) showed significant reductions in viral load and target cell lethality. 0% knockdown yielded the largest magnitude immune response, as quantified by extracellular IFN, however, 25% through 75% knockdowns showed a robust immune response of near equal magnitude despite the reduction in RIG-I activity. This suggests that RIG-I is robust against viral antagonism and plays a vital role in initializing the host's cytokine response to viral infections. Processes 2020, 8, x FOR PEER REVIEW 11 of 16 Figure 8. Simulations at varying levels of RIG-I knockdown (% knockdown or % kd of the k11 parameter). Here, 0% knockdown means zero NS1 antagonism, matching the dNS1PR8 strain results. Moreover, 100% kd is equivalent to total antagonism via NS1, matching the wild-type PR8 results.

Sensor Protein and JAK/STAT Originated Interferon Production
This model incorporates the production of IFN through both sensor protein action and the paracrine JAK/STAT pathway. These contributions to IFN production were isolated in silico, to determine the relative contribution of each pathway under different conditions. Figure 9 demonstrates these contributions for simulated wild-type PR8 and dNS1PR8 influenza infections. PR8 simulation showed almost complete antagonism of the sensor protein signaling, while dNS1PR8 simulation revealed a dynamic interplay between sensor proteins and paracrine IFN production, based on the infection stage. This analysis suggests that paracrine signaling is the major contributor to IFN production, especially in the presence of RIG-I antagonism, and that TLR activity alone is insufficient to trigger a strong immune response to infection.

Discussion
In this paper, an ODE model of the intracellular innate immune system's response to influenza infection was developed and used to evaluate system properties associated with viral lethality and Here, 0% knockdown means zero NS1 antagonism, matching the dNS1PR8 strain results. Moreover, 100% kd is equivalent to total antagonism via NS1, matching the wild-type PR8 results.

Sensor Protein and JAK/STAT Originated Interferon Production
This model incorporates the production of IFN through both sensor protein action and the paracrine JAK/STAT pathway. These contributions to IFN production were isolated in silico, to determine the relative contribution of each pathway under different conditions. Figure 9 demonstrates these contributions for simulated wild-type PR8 and dNS1PR8 influenza infections. PR8 simulation showed almost complete antagonism of the sensor protein signaling, while dNS1PR8 simulation revealed a dynamic interplay between sensor proteins and paracrine IFN production, based on the infection stage. This analysis suggests that paracrine signaling is the major contributor to IFN production, especially in the presence of RIG-I antagonism, and that TLR activity alone is insufficient to trigger a strong immune response to infection.

Figure 8.
Simulations at varying levels of RIG-I knockdown (% knockdown or % kd of the k11 parameter). Here, 0% knockdown means zero NS1 antagonism, matching the dNS1PR8 strain results. Moreover, 100% kd is equivalent to total antagonism via NS1, matching the wild-type PR8 results.

Sensor Protein and JAK/STAT Originated Interferon Production
This model incorporates the production of IFN through both sensor protein action and the paracrine JAK/STAT pathway. These contributions to IFN production were isolated in silico, to determine the relative contribution of each pathway under different conditions. Figure 9 demonstrates these contributions for simulated wild-type PR8 and dNS1PR8 influenza infections. PR8 simulation showed almost complete antagonism of the sensor protein signaling, while dNS1PR8 simulation revealed a dynamic interplay between sensor proteins and paracrine IFN production, based on the infection stage. This analysis suggests that paracrine signaling is the major contributor to IFN production, especially in the presence of RIG-I antagonism, and that TLR activity alone is insufficient to trigger a strong immune response to infection.

Discussion
In this paper, an ODE model of the intracellular innate immune system's response to influenza infection was developed and used to evaluate system properties associated with viral lethality and

Discussion
In this paper, an ODE model of the intracellular innate immune system's response to influenza infection was developed and used to evaluate system properties associated with viral lethality and IFN production dynamics. This model considers signaling components not available in previous work, uses unique data for parametrization and validation that perturbs distinct aspects of the pathways, and is numerically stable and well suited for continued adaption in future multiscale modeling efforts. The data used to inform the model were collected from human bronchial epithelial cell infections, with two influenza strains in vitro. Both strains induce an innate immune response, but the sensor protein RIG-I's activity is antagonized by only one strain, leading to vastly different interferon production and peak viral load. The model presented here provides estimates of the rates that regulate intracellular responses to RNA virus infection. The model can be used to assess how distinct RNA viruses impact IFN production, a key early step in activating the immune system, and is a valuable platform for determining how intracellular immune signaling may be distinctly regulated between influenza and coronaviruses.
A sensitivity analysis of the model revealed that IRF7 and IRF7P have significant control over the innate immune response. This strong single-protein control of the system suggests an area of further experimental investigation; viruses may try to impede these reactions as well as initial sensor protein action, to limit the innate immune response mounted against the invader. Moreover, pre-stimulation of the TLR4 pathway has been shown to lead to an earlier induction of IRF7P production and increased protection from deadly influenza infection [48]. The systems-level analysis here suggests that IRF7 is a potent target for the immune-targeted treatment of severe respiratory infection, both as a means of increasing host immune response and as a target for interference for the mitigation of cytokine storms.
By isolating paracrine and sensor protein originating production of IFN in silico, paracrine signaling is revealed to be responsible for most cytokine production and, thus, immune response. Early and strong sensor protein action serves to ignite this feedback loop in a dNS1PR8 strain. In both PR8 and dNS1PR8 strains, TLR's activity is too slow and of insufficient magnitude to significantly alter infection trajectory. When active, RIG-I has a profound effect on the peak viral titer. Effects would be more profoundly distinguished in low multiplicity of infection or in vivo cases, where uninfected neighboring cells would produce IFNs solely because of paracrine signaling. RIG-I is robust against antagonism, with 25-100% activity all reducing viral titer by about 97%. Full antagonism via NS1 is necessary for the virus to reach maximal peak. This level of antagonism is not achieved immediately, leaving a window for RIG-I to act before enough NS1 is produced inside the cell and its antagonism sets in. An area of future work is the relaxation of this simplifying assumption and the incorporation of NS1 production, based on viral load and the exploration of different RNA viruses' antagonism mechanisms.
Developing a minimal ODE model that minimized the computational expense of performing MCMC optimization in high dimensional space necessitated several simplifying assumptions. All type-I interferon species were grouped together in a single IFN equation. Several species' decay rates were estimated from literature [39][40][41], rather than being included in the MCMC optimization. RIG-I was assumed to be completely antagonized by the wild-type PR8 strain, which is supported by investigations of the NS1 protein [12,13]. The wild type simulation is stable at all times, since all simulated cells die. A validation study was performed by predicting a dNS1PR8 strain of influenza, using the same parameter set from training in the wild-type infection. The predictions of the validation case lend support to the model's capability to capture the interactions of interest without overfitting, and suggest that the model can predict responses outside of the training data. This validation case predicts only 20% cell lethality, despite total initial infection ( Figure 6, Cell Population), which leads to a mathematically unstable chronic inflammation state. This places limits on model interpretation beyond 24 h for NS1 knockout, or other reduced severity strains, since an infection in vitro is expected to be fully lethal in the absence of immune response and cellular regeneration. However, simulated cell death proportion acts as an indicator for virus severity.
This model does not stand alone in literature, and indeed must serve as a basis of further advancements. The action of drugs can be simulated by modifying only the relevant term(s) which the drug is thought to affect. Similarly, other strains or species of viral infection can be modeled by changing only the viral replication and clearance parameters. Both possibilities would be significantly computationally cheaper than the reparameterization of the entire model for each new system. If other parameters must change to reasonably fit a new data source, this could indicate previously unknown effects of the drug or infection on the species modeled herein, like NS1's well known antagonistic effects. The degree of this antagonism on RIG-I massively controls cytokine production trajectory and infection outcomes, as demonstrated here in silico by varying levels of RIG-I knockdown. This insight applies to all RNA viruses, with varying mechanisms of antagonism. Such predictions generated by the ODE model will provide insight into infection trajectory, disease outcome, and their manipulation by intervention.