key: cord-0779457-3yat87kq
authors: Mochan, Ericka; Sego, T. J.; Gaona, Lauren; Rial, Emmaline; Ermentrout, G. Bard
title: Compartmental Model Suggests Importance of Innate Immune Response to COVID-19 Infection in Rhesus Macaques
date: 2021-05-26
journal: Bull Math Biol
DOI: 10.1007/s11538-021-00909-0
sha: 484efbd1dc811985c8aabe04d409f0056126a024
doc_id: 779457
cord_uid: 3yat87kq

The pandemic outbreak of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has quickly spread worldwide, creating a serious health crisis. The virus is primarily associated with flu-like symptoms but can also lead to severe pathologies and death. We here present an ordinary differential equation model of the intrahost immune response to SARS-CoV-2 infection, fitted to experimental data gleaned from rhesus macaques. The model is calibrated to data from a nonlethal infection, but the model can replicate behavior from various lethal scenarios as well. We evaluate the sensitivity of the model to biologically relevant parameters governing the strength and efficacy of the immune response. We also simulate the effect of both anti-inflammatory and antiviral drugs on the host immune response and demonstrate the ability of the model to lessen the severity of a formerly lethal infection with the addition of the appropriately calibrated drug. Our model emphasizes the importance of tight control of the innate immune response for host survival and viral clearance. SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at 10.1007/s11538-021-00909-0.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has continued to be a worldwide health crisis since late 2019. Symptoms of coronavirus disease-2019 (COVID-19) include flu-like symptoms (cough, sneeze, fatigue, sore throat, etc.) as well as fever, shortness of breath, and even digestive symptoms. Most COVID-19 cases are mild and resolve on their own, but around 20% of cases lead to severe pathologies like pneumonia and acute respiratory distress syndrome (ARDS) .

ARDS and other excess symptoms and mortality are caused by two primary mechanisms: direct lung injury by viral infection, and overstimulation of inflammation causing a cytokine storm. Cytokine storm in COVID-19 has been shown to cause a high expression of many pro-inflammatory cytokines by effector cells Yang et al. 2020; Zhang et al. 2020) .

The balance and regulation of pro-and anti-inflammatory mediators in the wake of COVID-19 infection is paramount to the resolution of the infection. Proinflammatory regulators are key to driving the initial immune response to the virus, but COVID-19 can easily overstimulate this response, leading to a cytokine storm and excess damage to the epithelium Hu et al. 2020) .

Anti-inflammatory drugs can be used to help mitigate this excess inflammatory response. Proper calibration of the initiation, duration, and dosage of this drug therapy is key to controlling the inflammation. Providing this drug for too long can dampen the host's ability to heal from the infection, but too little drug can be insufficient to inhibit the inflammation enough to prevent cytokine storm (Nile et al. 2020; Ye et al. 2020; Zhang et al. 2020) .

Here we present an ordinary differential equation (ODE) model of the inflammatory response to COVID-19 infection. Mathematical modeling has long been used to describe the immune response to lung infection (Baccam et al. 2006; Hernandez-Vargas et al. 2014; Miao et al. 2011; Mochan-Keef et al. 2015; Mochan et al. 2018 Mochan et al. , 2014 Price et al. 2015; Saenz et al. 2010; Smith et al. 2011) . Our model focuses on the importance of the inflammatory response to viral infection, using feedback from the pro-inflammatory mediators to suppress viral growth and induce incidental damage to the epithelium. The model also features a distinction between the lower and upper respiratory tracts, allowing us to capture spatial differences in viral distribution and growth. Loss of smell or taste is an early potential marker of COVID-19 infection, and the virus has been shown to cause significant changes to mucociliary action in the upper respiratory system (Koparal et al. 2020) . Lower respiratory symptoms such as coughing and respiratory distress indicate a severe infection and substantial damage of the host's epithelium . The structure of our model allows us to separate these types of symptoms to investigate these dynamics in a variety of simulated infections.

This model can simulate multiple phenotypes in respiratory viral infection: asymptomatic infection, mild infection with fast recovery, severe but survivable infection, lethal infection due to excess viral titers, and lethal infection due to excess inflammation. We demonstrate the parameter bifurcations in which these modes may occur, and we fit this model to SARS-CoV-2 data from rhesus macaque infection (Munster et al. 2020) . Viral shedding and spread in macaques have been shown to be similar to that in human subjects (Josset et al. 2012; Munster et al. 2020) , and macaques were identified as a model organism for MERS coronavirus infections (Falzarano et al. 2013) . Previous work has also shown macaques accurately replicated mild COVID-19 infection in humans ).

This novel model of the intrahost immune response to respiratory viral infection follows the time-evolution of 9 variables representing virus (V U and V L ) and healthy (H U and H L ) and infected (I U and I L ) epithelial cells in the upper and lower respiratory tracts, as well as global levels of pro-inflammatory (F) and anti-inflammatory (G) mediators, and damage (D). The model consists of 9 differential equations with 24 parameters representing the interactions between these components. Figure 1 summarizes the interactions represented by the model, and Eqs. (1)-(9) detail the model.

Virus begins in both the upper and lower respiratory tract, in accordance with the experimental conditions. Macaques were infected orally, intranasally, and intratracheally to ensure the virus reached both compartments immediately upon infection. Virus invades healthy epithelial cells in the respiratory tract, creating a population of infected epithelial cells, which produce and release more virions to continue the Fig. 1 Graphical summary of the interactions represented in the model. Infection begins when virus is introduced in the upper and lower respiratory tracts, which are treated as two compartments in the model. Epithelial cells in each compartment do not move, but virus can progress from the lower compartment to the upper along the mucociliary elevator. Virus infects epithelial cells, which in turn produce and release free virus. Infected cells also initiate the stimulation of pro-inflammatory mediators. Inflammation causes excess damage of the epithelium and initiates the anti-inflammatory response to contain the inflammation infection. Free virus is removed by one of three mechanisms: (1) internalization by epithelial cells (parameter γ VH ), (2) clearance by nonspecific immunity (parameters a VU and a VL ), and (3) clearance by antibodies (parameters α VU and α VL) . In addition, virus moves from the lower respiratory tract to the upper as mucus and cilia in the lungs mechanically push the pathogen upwards to remove it from the host via the upper respiratory tract (Koparal et al. 2020) .

The presence of infected cells and cellular damage activates the pro-inflammatory arm of the immune response. Virus production is inhibited by the presence of pro-inflammatory mediators, which limit the ability of the infected cells to actively release virions. Inflammation causes a subsequent increase in anti-inflammatory mediators, which limits the release of more inflammation to protect the host against cytokine storm.

The pro-inflammatory mediators also cause excess damage to the host, killing epithelial cells in addition to the virus. Unlike target cell-limited models, epithelial cells can regenerate in this model to allow the host to heal the damage caused by the infection. In our model, survival is characterized by a return to the initial healthy state, meaning all healthy cells have regenerated, virus is cleared in both compartments, and inflammation and damage heal. Death is defined as a depletion of healthy cells in either compartment.

We validate the interactions and parameters presented in the model by fitting to experimental data from Munster et al. (Munster et al. 2020) . These authors infected four rhesus macaques with an inoculum of 4.6 log 10 TCID 50 /ml SARS-CoV-2 in both the upper and lower respiratory compartments and measured the progression of their disease for 21 days post-infection (dpi). Four other rhesus macaques were also infected but killed at day 3 post-infection, and we did not include these early measurements in our results, focusing only on those who survived a full 3 weeks after infection.

Measurements for viral RNA levels in both the nasal cavity and lungs were used for upper and lower compartment virus measurements, respectively. Inflammation is fit to data for IL-15, a pro-inflammatory cytokine which activates natural killer (NK) cells and T cells after virus infection . Damage, as a measure of host health, is fit to clinical symptoms scores, which were measured daily on the macaques. To represent the disparity in size of the two compartments, the total number of cells in the upper respiratory tract is taken to be about a third of the number of cells in the lower respiratory tract, so H U (0) = 0.33 and H L (0) = 1. Figure 2 shows the outcome of fitting the model to this data. Viral titers peak at 1 dpi, followed by about 2 weeks of slow decay in the nasal passages until viral clearance about 18-21 dpi. Because this peak is only about one order of magnitude higher than the inoculum, there is very little damage and infected cell creation in the upper compartment.

In the lower compartment, there is much more damage and infected cell generation, as the virus grows by about three orders of magnitude within the first 2 days. Within about 1 week post-infection, the virus has almost completely depleted from the lower compartment, and the healthy cells have returned to their baseline level, signifying that the host has survived the infection.

Inflammation (F) also peaks early, in accordance with the experimental data; only the day 1 data were statistically significant for the pro-inflammatory cytokines and chemokines in the macaques (Munster et al. 2020) . Damage (D) and the anti-inflammatory signals (G) peak after the pro-inflammatory and slowly return to their baseline values.

Clinically, COVID-19 has been shown to have particularly adverse effects on patients of advanced age or with preexisting conditions Wang et al. 2020; Yang et al. 2020; Zhou et al. 2020) . We evaluated the ability of the Page 6 of 26 model to predict host-specific dynamics, such as strength and speed of the immune response, by changing select model parameters.

Elderly patients have experienced more severe symptoms and outcomes than younger hosts with COVID-19 infection Verity et al. 2020 ). An aged host will typically replicate epithelial cells at a slower rate than a young, healthy host (Tomasetti et al. 2019) , controlled in our model by parameter b HD . Immunosenescence, dysregulated immune function related to the aging process, is also a major contributor to morbidity and mortality in the elderly population (Franceschi et al. 2000a; Ginaldi et al. 2001) . Adaptive immunity and innate immunity are impacted by aging. Aged hosts have been shown to experience a delayed onset and prolonged activation of inflammation, leading to excess symptoms and damage to the epithelium in a phenomenon termed "inflamm-aging" (Franceschi et al. 2000b; Ginaldi et al. 1999 ).

Damage to the epithelium can be sensed by nearby cells, which send signals to upregulate the innate immune response (Chen et al. 2018 ). In our model, this behavior is governed by parameter FD . Figure 3 shows the effect of changing both b HD and FD over several orders of magnitude, reflecting changes that may occur in the host with age. The figure demonstrates the end behavior of each variable in the model, showing the value of these variables 200 dpi. When there is very low stimulation of the inflammatory response, virus can grow for the first 3-4 days post-infection, rather than reaching its peak within a day, as in the calibrated parameters of Fig. 2 . Without the tight control of the time to the virus peak, the infection will become lethal.

As FD increases, the inflammatory response gets stronger, and the initial growth of the virus becomes so low that the host can resolve the infection quickly. However, when FD is on the order of 10 1 , the infection becomes chronic. Virus does not deplete in the lower compartment, which leads to sustained inflammation and higher symptoms and sickness scores. There is also a narrow area of parameter space in which the upper virus will not clear, leading to chronically depressed healthy cells in the upper respiratory tract.

Aging also likely impacts parameter a HF , which controls the depletion of healthy cells by inflammation. Increases in this parameter would simulate the excess damage caused by overactive inflammation. Figure 4 demonstrates the effect of gradually increasing this parameter, which essentially divides the parameter space into three sections. At low values of a HF , the virus in the lower compartment cannot be cleared, as the inflammatory response is not stimulated quickly enough to control the viral growth. This leads to a permanent loss epithelial cells in the lower respiratory tract, and a sustained inflammatory response over 2 months post-infection.

At intermediate values of a HF , the stimulation of inflammation is high enough to contain the virus growth but also low enough to resolve within 60 dpi. As a HF rises to high levels, inflammation will grow very quickly to a high magnitude, causing a substantial amount of damage and continued positive feedback to the inflammatory response, likely leading to hospitalization and serious outcomes. Depending on the value of b HD , this high damage may or may not be accompanied by depletion of epithelial cells in the upper and lower respiratory tracts.

There is also a narrow area in which the infected cells in the upper compartment reach a nonzero steady state. The virus in the upper compartment does not decrease quickly enough and rebounds after about 15 dpi, resulting in damped oscillations and a nonzero steady state for H U , I U , and V U .

The value of parameter a HF is clearly a significant factor in the resolution of each simulation. In Fig. 5 , we demonstrate how changing only this parameter can result in three distinct regimes of behavior in the F − H l plane. When a HF = 0.015 (red line), the infection clears and the host resolves the infection fully, as in our calibrated parameter set shown in Fig. 2 . Raising a HF to 0.025 creates a chronic infection, with a high level of inflammation and permanently depressed number of epithelial cells in the lower compartment (blue line). Finally, raising a HF again to 0.035, the simulation predicts inflammation-induced death of the host. The inflammation reaches a high peak value and depletes all the cells in the lower compartment (black line).

Sensitivity to age-related changes to the immune response. Parameters b HD and FD control two important effects of immunosenescence: age-related slowing of cellular replication, and excess inflammation and damage known as inflamm-aging. At low levels of b HD , the healthy cells will deplete and lead to death of the host. As b HD Increases, the healthy cells can regenerate quickly enough to offset the viral load and transform the infection from lethal to survivable. High levels of FD can initiate a strong inflammatory response to eradicate the virus from the lower compartment even with low levels of b HD . In the middle ranges, the inflammation reaches a chronic state of consistently elevated F, D, and G, and a nonzero steady state for H and I

Sensitivity to cell replication rates and death rates. At high levels of b HD , the epithelial cells can regenerate quickly, but more damage accumulates along the way. As a HF increases, we see a narrow region in which the system can recover regardless of the magnitude of b HD . If a HF is either too high or too low, the host will sustain too much damage to recover

Sensitivity of the lower respiratory tract to inflammation-induced damage. Changing only the parameter a HF , which controls the excess damage caused to the epithelium by the inflammatory response, can change the outcome of the infection. Low levels of a HF will lead to resolution of the infection and survival of the host (red line). As a HF increases, the severity of the infection also increases, becoming chronic (blue line) or lethal (black line)

We next explored the paired sensitivity of some parameters that differ between the upper and lower compartments. First, we studied the impact of changing the a VU and a VL parameters, which control the strength of the adaptive immune response against the virus in the upper and lower compartments, respectively.

As Fig. 6 shows, high levels of VU are required to clear the virus in the upper compartment. Limiting the adaptive response in that compartment will lead to a sustained viral titer and chronic infection. The lower compartment is not sensitive to the adaptive response in these scans. As shown in Fig. 2 , the calibrated parameters lead to a fast depletion of virus in the lower compartment, with all virus cleared from the lungs by day 9 post-infection. That clearance is more dependent on the innate and inflammatory response to virus than the adaptive, as the innate immune responses are felt much earlier in the infection progression.

Spatial differences also exist in the inflammatory response between the two compartments. Parameters k VFU and k VFL control the suppression of the viral growth by Fig. 6 Sensitivity to adaptive immune responses. At high levels of a VU , the host can clear the virus from both compartments and return to a baseline healthy state. Without a strong adaptive response in the upper compartment, the virus will remain elevated for months after infection, indicating a chronic infection state inflammation in the upper and lower respiratory tracts, respectively. High values of k VFL are beneficial to the lower respiratory tract, as the inflammation depletes the virus and allows the healthy cells to heal, regardless of the strength of the inflammatory response in the upper compartment. However, this immediate depletion is not necessarily a benefit to the upper respiratory tract. Without the lung pathogen present to initiate a strong immune response, the virus in the upper compartment will not be able to clear, and the host will experience a chronic upper respiratory infection (red and yellow region in H U and I_U subplots). Results of this parameter scan are demonstrated in Fig. 7 .

We also tested the model's ability to simulate drug intervention, both for infections that are lethal due to excess virus, and for infections lethal due to excess inflammation. We adapted our equations to include a drug dosage curve which can be customized based on the day the drug intervention begins, the length of time the drug is administered, and the maximum strength of the drug.

Sensitivity to inflammation-driven suppression of virus growth. Elevated suppression of viral growth in either compartment will not necessarily lead to a better outcome for the host. Excess effects from inflammation can lead to a chronically inflamed state in the upper respiratory compartment, as damage increases and feeds back to the pro-inflammatory response The drug dose is simulated with a step function:

The height of the step function, s max , represents the maximum efficacy of the drug given. The drug is administered at time t on and continues for t dur days.

We first tested whether our model could replicate the effects of administering an anti-inflammatory drug to a patient experiencing cytokine storm. We adapted the inflammation equation to include the drug dosage as an additional anti-inflammatory element:

We then adjusted the parameters of the model from those used in the fitting procedure to create a simulation in which the model predicts the host to die from excess inflammation (Fig. 8) . Inflammation peaks around 800, much higher than the peak of about 70 we see in the F trajectory for the calibrated parameter set (Fig. 2) . Damage and anti-inflammatory signals also remain elevated longer and peak at much higher values. This prolonged inflammatory response causes excess damage to the epithelium in both compartments. The healthy cells in the upper compartment are completely eradicated, and in the lower compartment, while the cells will eventually return to their baseline level, the damage is sustained for weeks.

We performed parameter sweeps with s max , t on , and t dur to determine the optimal calibration of the drug dosage. First, time of onset was varied between 1 and 5 days post-infection, with s max = 50 and t dur = 21 days. When the drug is given between 0 and 2.3 dpi, the inflammation is initially depressed too much, which results in extended growth of virus in the lower compartment, which will then lead to virusdriven death of the host.

When administered between days 2.3 and 4.6, the infection becomes survivable for the host. The inflammation peaks much lower than in the other simulations, on a level like that in the calibrated parameter set. The damage immediately begins to heal, and the healthy cells return to their baseline level in both compartments. These experiments demonstrate the timing of the drug dosage is critical to control of the disease, especially in patients with cytokine storm and ARDS.

Administering the drug too late (beyond day 4.6) means the inflammation will not be controlled soon enough, and the cells in the upper compartment will die out regardless of the strength of the drug administered. An example trajectory from each of these three regions of behavior is shown in Fig. 9a .

We also explored the impact of changing the amount and duration of drug delivery while keeping time of onset constant. With t on = 3 dpi, we varied s max from 10 0 to 10 3 and t dur from 0 to 30 days. The results of this experiment are plotted in Fig. 9b . The parameter set used in this experiment will always lead to

clearance of the virus in the upper compartment. However, depletion of cells in the upper compartment may occur anyway if the duration of the drug delivery is too short. As the duration increases, healthy cells may either reach a steady-state value below their baseline or return fully to their healthy baseline state. When the epithelium reaches its chronically damaged state (red areas on H U and H L plots), this corresponds to a nonzero steady state for F, D, and G, indicating a permanently elevated inflammatory response.

Inflammation-induced death can also impact the lower compartment. We again adapted the calibrated parameters to simulate an infection in which the virus is cleared, but there is chronic inflammation (parameter values provided in Table S2 ), this time focusing on the health of the lower respiratory tract. Here we introduce a long-term anti-inflammatory drug. As in the previous simulations, the time of drug administration is the most significant factor in the host's likelihood of survival. As Fig. 10 demonstrates, delaying the onset of the drug delivery to either day 2 (green line) or day 3 (blue line) will lead to a full recovery of the host. When delivered only 1 dpi (red line), the virus will still kill the host, as Fig. 8 Simulation of death due to excess inflammation. Parameters in the model were adjusted to create a simulation in which the virus is cleared from the host quickly, but inflammation and damage remain elevated long after the virus is eradicated from the body. This excess inflammation leads to anosmia and chronic upper respiratory symptoms, as shown by the depletion of healthy cells in the upper compartment, prolonged depression of healthy cells in the lower compartment, and high levels of damage. Parameters used to generate this simulation are given in Table S1

Effect of anti-inflammatory drug on simulated cytokine storm. a Early drug dosing tempers the inflammatory response too quickly, leading to a large peak in the virus growth and depletion of epithelial cells in both compartments. When the drug is given too late, although the virus can be cleared from the system by the immune system, the excess damage overtakes the system and the healthy cells will deplete from at least one compartment. Given at day 7, the drug has the same effect on the virus as the day 9 simulation, but the damage remains low enough that the host will survive the infection, with epithelial cells returning to baseline in both compartments.

b Sensitivity to magnitude and duration of drug delivery. With too short a duration of drug administration, the infection remains lethal, regardless of the strength of the drug. At low to medium values of s max , as the duration increases, the infection will either remain chronic or heal fully. When s max is above 10 2 , the virus in the lower compartment will not be regulated enough by the remaining inflammation, and it will reach a high steady-state value, leading to death of the host the immune response will be too suppressed to destroy the virus. This results in the virus destroying all the cells in the lower compartment. If the drug is instead delayed to day 4 (black line), it is already too late to overcome the inflammation and the host will experience a chronic infection that never fully resolves. This permanent deterioration of the epithelial cells in the lungs could simulate the chronic respiratory problems experienced by many of those who contract COVID-19.

We next tested whether the model could simulate the action of administering an antiviral drug to a patient overwhelmed by high viral titers. While high viral titers are not a common cause of death in COVID-19 infection, excess virus can commonly be lethal in other respiratory viral infections, such as influenza (De Jong et al. 2006; Price et al. 2015; Toapanta and Ross, 2009 ). Using the same drug delivery step function, we adapted the virus equations in both compartments:

Effect of anti-inflammatory drug on simulated virus-induced death of the lower respiratory compartment. If the anti-inflammatory is administered to the host either 2 or 3 dpi, the virus can be cleared fully by the inflammatory response which can then be quelled by the drug; the epithelial cells will fully recover. If the drug is provided too early, the host will die due to an insufficient immune response, and if the drug is given too late, the infection becomes chronic, leading to permanently depressed levels of healthy cells in the lungs. Parameter values used to simulate these trajectories are provided in Table S2 We again began with the parameters calibrated to the macaque data, but changed a small subset of parameters to simulate virus-induced death. Parameters used to generate this lethal scenario are listed in Table S3 . Figure 11 shows what this initial lethal simulation looks like. Healthy cells in the upper compartment are depleted after about 1 week due to the high sustained viral titer in the upper respiratory tract.

We then performed parameter sweeps with the drug dosage parameters to determine the appropriate calibration for the antiviral drug (Fig. 12) . As in the antiinflammatory simulations, we keep s max = 50 and t dur = 21 and vary only t on from 1 to 9 dpi. When the drug is provided by 5 dpi, the healthy cells in the upper compartment may begin to regenerate and the host can eventually survive the infection, though the reduction of the virus occurs very slowly. However, after 5 dpi, the healthy cells have already fallen below their threshold value of 0.05, and it is impossible for them to recover. Figure 11a shows example trajectories in each of these two regions of behavior, and Fig. 12b demonstrates the sensitivity of this system to changes in the magnitude and duration of the drug delivery. Only H U demonstrates any sensitivity to the antiviral drug; all other variables will always return to 0 using these parameters. The upper respiratory tract will only heal given a sufficiently strong drug over at least 10 days; if the drug is stopped too early, the host will clear the virus from the compartment but will not regain the epithelial cells, resulting in chronic anosmia. a If the drug is administered within the first week post-infection, the cell damage can be reversed and the host can survive, though the eradication of the virus will be slow and upper respiratory symptoms will persist for several weeks.

After 5 dpi, H U falls below its threshold level and cannot regenerate, leading to death of the host.

b Parameters s max and t dur were scanned while keeping time of onset at a constant 5 days. Other parameters are identical to those used to create the trajectories in Fig. 11 , where host death is defined by loss of epithelial cells in only the upper compartment. Sensitivity to drug delivery is also contained to the upper compartment. Healthy cells in the upper respiratory tract will recover only with a prolonged, strong drug delivery. When either s max or t dur is too low, the system cannot recover

We developed a compartmental differential equation model of intrahost immune dynamics of COVID-19 infection. The 9-equation model is fitted to experimental data measured on rhesus macaques infected with SARS-CoV-2 in both the upper and lower respiratory tract. Our model captures the dynamics of the immune response to the viral infection, as well as the movement of virus from the lower to upper compartment along the mucociliary elevator. All four macaques survived the infection, clearing the virus from both compartments and returning the proinflammatory mediators to their baseline levels. The model can qualitatively and quantitatively replicate these dynamics. Systemic symptoms like fever, aches, chills, and malaise are primarily attributed to pro-inflammatory cytokines like IL-6 and type I interferon. These global symptoms have been shown to peak early during other respiratory viral infections like influenza (Doyle et al. 1998; Fritz et al. 1999; Wadowsky et al. 1995) . Systemic symptoms are closely tied to inflammation, both in the structure of the model and in experimental results. Pathology and imaging performed on the macaques showed lesions and pulmonary edema during infection, indicating extensive damage to the lungs as a result of the virus. Clinical scores, used as a proxy for lung damage in our data fits, peak around 5 dpi. These scores amalgamate a variety of clinical parameters, such as weight loss, fever, respiration rate, and pulmonary edema. In the experiments, all macaques were fully recovered between 9 and 17 dpi. Our model shows a steady decline in lung damage from day 3 onward, with health recovered nearly to baseline by day 21 post-infection. Importantly, the model can also replicate lethal infection dynamics, including infection that become lethal due to excess viral titers in either compartment, or infections that are lethal due to cytokine storm. Given that about 20% of COVID-19 infections have been severe, with about 5% of infections becoming lethal, accurate replication of these dynamics is vital to the translatability of the model.

All parameter explorations demonstrate the importance of a tight calibration of the immune response for survival of the host. If the inflammatory response cannot react quickly and with a large magnitude, the initial rise of the virus will exceed the capacity of the immune response and quickly exhaust the epithelial cells. With too strong of an immune response, however, the excess inflammation will cause high levels of damage, corresponding to more severe symptoms and possibly inflammation-induced death of the host.

We have also shown that the model demonstrates the ability to transform formerly lethal infection scenarios into survivable scenarios with the addition of an antiviral or anti-inflammatory drug. Our drug dosage curves verify that the time of onset of the drug delivery is the most important factor to improving the outcome of the disease in severe or lethal infection conditions (Sego et al. 2020) . Chances of recovery from multi-organ dysfunction after COVID-19 have been low, so prompt and adequate application of drug dosage is critical to survival of the patients (Siddiqi and Mehra, 2020).

Corticosteroids have been suggested as a possible treatment for severe viral infection which induces cytokine storm (Brun-Buisson et al. 2011; Kudo et al. 2012; Liu et al. 2016; Zhang et al. 2020) . However, there are relatively few studies in which steroid treatments have been tested on human patients infected with influenza or other respiratory viral infections. Ultimately, effective treatment must be optimized to promote viral clearance while limiting inflammation and tissue degradation. Previous studies have shown that, to be effective in immunomodulation after SARS infection, steroids should be administered early and continued for at least 10 days post-infection (Carter 2007) .

Studies of corticosteroid treatment of viral infection have not shown significant impacts on mortality rates (Bernard 2002) . Our anti-inflammatory drug simulations align with these results; the window of opportunity to change an infection from lethal to nonlethal is narrow, and in many cases, the anti-inflammatory drug actually made the infection worse as the virus could not be controlled without an early, strong inflammatory response. Recent studies suggest lower doses of steroids are more effective than the traditional short-term, high-dosage steroid treatment of the past (Annane et al. 2004; Bernard 2002) . Some have suggested that steroids are most effective as a pre-treatment, which is impractical in human subjects (Carter 2007) .

Some authors have suggested that steroid therapy should be combined with oseltamivir (anti-viral) treatment (Carter 2007; Chotpitayasunondh et al. 2005) , as steroids can be ineffective in clearing the virus from the host. In our model, the antiviral drug again must be administered early in the infection progression in order to save the host, but the window of opportunity is certainly wider than in the antiinflammatory simulations.

Future extensions of this work may include expanding upon the immune components included explicitly in this model to incorporate more inflammatory markers, such as macrophages or chemokines, as the current model is naturally a simplification of the complex interactions and dynamics that comprise the inflammatory response. Adaptive immune elements such as T cells or antibodies may also be included explicitly in future iterations of the model. The biological accuracy of these analyses is limited by the sparsity of the data available to which to fit the model. Ideally, future iterations of the model will incorporate these immune dynamics explicitly, fitting the model to detailed longitudinal data for viral infections.

While macaques have been used as models of the human immune response (Chiou et al. 2020; Falzarano et al. 2013; Shan et al. 2020) , future work may also include fitting the model to human patient-level data for COVID-19 or other severe respiratory viral infections. The model could be used to demonstrate important differences in host immune responses, clarifying reasons why some patients experience far more serious outcomes than others after COVID-19 infection. In addition, our drug treatment simulations can eventually be matched to experimental data for antiviral or anti-inflammatory drug treatment of SARS-CoV-2 infection to validate our predicted outcomes. The model does replicate many expected qualitative behaviors of limiting inflammation and diminishing viral titer growth, but these conclusions would certainly be improved with experimental data to which to compare the model's behavior.

In summary, our model can accurately replicate quantitative and qualitative behaviors of a host experiencing severe respiratory viral infection, such as that seen in the current COVID-19 pandemic. The model can simulate nonlethal infections, infections lethal due to excess viral load, or infections lethal due to excess inflammation, by changing only a few parameters. Drug therapy can also be added to the model equations to simulate either an antiviral or anti-inflammatory drug, depending on the cause of death to the host. These drugs can be used to transform a lethal simulation to a nonlethal one, provided the drug is calibrated correctly and administered early, in accordance with experimental studies.

Experimental data used for model calibration were downloaded from Munster et al. (Munster et al. 2020) . In these experiments, 8 macaques were challenged with 4 × 10 5 TCID 50 /ml SARS-CoV-2 via intranasal, intratracheal, and oral inoculation. Four of the macaques were predetermined to be killed at day 3 post-infection, and the remaining four were observed for 21 days post-infection. We used data from only the four macaques followed for 21 dpi to characterize the full time-course of the infection.

Viral loads in the nasal cavity and cytokine levels were measured on days 1, 3, 5, 7, 10, 12, 17, and 21 post-infection, while viral loads in the lungs were measured only on days 1, 3, and 5 post-infection. Clinical symptoms scores were measured daily. We used IL-15 as the marker for inflammation; all cytokines measured in the macaques were statistically significant only on day 1. Symptoms scores were used as a marker for damage, as increased damage to the epithelium is a sign of systemic symptoms (Schleimer et al. 2009 ).

Equations (1)-(9) describe the ODE model developed here. The upper respiratory tract dynamics are summarized in Eqs. (1)-(3). Virus is produced by infected cells, tempered by the presence of inflammatory mediators like type I interferon. Free virus is taken up by epithelial cells at rate γ VH to create new infected cells. Two immune responses are also featured in the model. We include a saturating term representing the nonspecific clearance of the virus at rate a VU (Reynolds et al. 2006 ) . This rate includes mechanisms such as mechanical clearance of virus from the lungs via mucociliary action, phagocytosis by resident macrophages, defensins, and nonspecific antibodies present in the lungs (Gomez and Prince, 2008) . These defenses react quickly to pathogen but are also easily overwhelmed by a large inoculum, so the term saturates at high levels of virus.

The model also includes a linear term representing specific clearance by adaptive immunity at rate α VU . Virus can also be added to the upper compartment as it travels along the mucociliary elevator from the lungs to the upper respiratory compartment, controlled by parameter β.

Healthy epithelial cells can become infected after interaction with virus. The inflammation can also lead to excess death of the healthy cells, through the release of reactive oxygen species after innate immune cells like neutrophils are activated. Healthy cells regenerate proportionally to the number of cells that have died ( H max − H U − I U ). Parameter ε represents a threshold value for healthy cells. If H falls below ε, this term becomes negative and the cells can no longer regenerate, thus stabilizing the dead state.

Equations (4)-(6) describe the dynamics of the lower respiratory tract. Most interactions are identical to those in the upper respiratory tract. Parameters controlling the strength of the immune response against the virus can take on different values in the lower compartment than the upper. The virus in the lower compartment can move to the upper compartment, but virus does not move in the opposite direction.

Dynamics for global variables (inflammation, damage, and anti-inflammation) are described in Eqs. (7)-(9). Epithelial cells sense damage nearby and release pro-inflammatory signals. Infected cells also release more type I interferon, further increasing the inflammatory signals in the respiratory tract. Inflammation is tempered by the presence of anti-inflammatory mediators (G).

(1)

Damage is created by the inflammation and the death of infected cells and is modeled with a Hill function. Damage can be healed, governed by parameter a D . Anti-inflammatory mediators are created directly in response to pro-inflammatory signals, and they deactivate over time as the infection resolves.

An important facet of this model is the ability to maintain a sustained inflammatory response, even in the absence of virus. This feature is important in that COVID-19 often leads to cytokine storm and excess inflammation long after the viral infection has been resolved. To ensure our model can replicate this behavior, we studied the H-F subsystem of the model. By making steady-state assumptions for the remaining variables in the system, we can construct the subsystem:

In Fig. 13 , we demonstrate the three stable fixed points of the H-F subsystem of the model. The point at F = 0, H = 1 shows the stability of the healthy state, and the point at F = 0, H = 0 shows the stability of the death state. There is also one fixed point at nonzero values for F and H, demonstrating a nonzero steady state for these variables. Note that there are two saddle equilibria whose stable manifolds separate the three stable equilibria into distinct basins. We do not depict them here as they will be quite different in the full system; the steady state reduction has been applied here to make the equilibria clearer.

We calculated a mean and standard deviation for each data point across the 4 variables with experimental data and fit to this pooled data. Initial conditions were = FD D − a F F 1 + k FG G generated from the experimental conditions. Virus began at 4 × 10 5 TCID 50 /ml in both compartments. All other variables were assuming to start at their natural initial conditions if no infection was present; infected cells, damage, inflammation, and anti-inflammation all begin at 0, and healthy cells begin at their maximum value. The best fit parameters were generated using MATLAB R2019a with function fminsearchbnd. Error was calculated as the sum of squared errors between the model output and the experimental data. Table 1 gives the list of parameters and their values when fit to the experimental data.

In addition, there are three parameters used only in the drug dosage simulations, but not in the parameter fitting procedure. These parameters are summarized in Table 2 . Acknowledgements T.J.S. acknowledges funding from National Institutes of Health grants U24 EB028887 and R01 GM122424 and National Science Foundation grant NSF 1720625. BE acknowledges funding from the National Science Foundation 1951099. The authors would like to thank James Glazier for bringing this group together. 

Corticosteroids for severe sepsis and septic shock: a systematic review and meta-analysis

Kinetics of influenza A virus infection in humans

The International Sepsis Forum's controversies in sepsis: corticosteroids should not be routinely used to treat septic shock

Early Corticosteroids in Severe Influenza A/H1N1 Pneumonia and Acute Respiratory Distress Syndrome

A rationale for using steroids in the treatment of severe cases of H5N1 avian influenza

Inflammatory responses and inflammation-associated diseases in organs

Rhesus macaques as a tractable physiological model of human ageing

Human disease from influenza A (H5N1)

Fatal outcome of human influenza A (H5N1) is associated with high viral load and hypercytokinemia

Effect of rimantadine treatment on clinical manifestations and otologic complications in adults experimentally infected with influenza a (H1N1) virus

Treatment with interferon-α2b and ribavirin improves outcome in MERS-CoV-infected rhesus macaques

Human immunosenescence: the prevailing of innate immunity, the failing of clonotypic immunity, and the filling of immunological space

Inflamm-aging. An evolutionary perspective on immunosenescence

Nasal cytokine and chemokine responses in experimental influenza A virus infection: results of a placebocontrolled trial of intravenous zanamivir treatment

The immune system in the elderly

Immunosenescence and infectious diseases

Bacterial interactions with the airway epithelium. In: Proud D (ed) The pulmonary epithelium in health and disease

Clinical characteristics of coronavirus disease 2019 in China

Effects of aging on influenza virus infection dynamics

The cytokine storm and COVID-19

Messaoudi I (2012) Increased viral loads and exacerbated innate host responses in aged Macaques infected with the 2009 Pandemic H1N1 Influenza A virus

Assessment of mucociliary clearance as an indicator of nasal function in patients with COVID-19: a cross-sectional study

Systemic corticosteroids and early administration of antiviral agents for pneumonia with acute wheezing due to influenza A(H1N1)PDM09 in Japan

Molecular immune pathogenesis and diagnosis of COVID-19

Clinical features of COVID-19 in elderly patients: a comparison with young and middle-aged patients

The cytokine storm of severe influenza and development of immunomodulatory therapy

On identifiability of nonlinear ODE models and applications in viral dynamics

A three-tiered study of differences in murine intrahost immune response to multiple Pneumococcal Strains

A systems and treatment perspective of models of influenza virus-induced host responses

A mathematical model of intrahost pneumococcal pneumonia infection dynamics in murine strains

) Respiratory disease in rhesus macaques inoculated with SARS-CoV-2

COVID-19: pathogenesis, cytokine storm and therapeutic potential of interferons

The inflammatory response to influenza A virus (H1N1): an experimental and mathematical study

A reduced mathematical model of the acute inflammatory response: I. Derivation of model and analysis of anti-inflammation

Dynamics of influenza virus infection and pathology

Epithelium, inflammation, and immunity in the upper airways of humans: studies in chronic rhinosinusitis

A modular framework for multiscale multicellular spatial modeling of viral infection, immune response and drug therapy timing and efficacy in epithelial tissues: A multiscale model of viral infection in epithelial tissues

Infection with novel coronavirus (SARS-CoV-2) causes pneumonia in Rhesus macaques

Effect of 1918 PB1-F2 expression on influenza A virus infection kinetics

Impaired immune responses in the lungs of aged mice following influenza infection

Cell division rates decrease with age, providing a potential explanation for the age-dependent deceleration in cancer incidence

Functions of IL-15 in anti-viral immunity: multiplicity and variety

IL-15 participates in the respiratory innate immune response to influenza virus infection

Estimates of the severity of coronavirus disease 2019: a model-based analysis

Effect of experimental influenza A virus infection on isolation of Streptococcus pneumoniae and other aerobic bacteria from the oropharynges of allergic and nonallergic adult subjects

Clinical characteristics of 138 hospitalized patients with 2019 Novel Coronavirus-Infected Pneumonia in Wuhan, China

Clinical course and outcomes of critically ill patients with SARS-CoV-2 pneumonia in Wuhan, China: a single-centered, retrospective, observational study

The pathogenesis and treatment of the `'Cytokine Storm' in COVID-19

The use of anti-inflammatory drugs in the treatment of people with severe coronavirus disease 2019 (COVID-19)

Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study