Compartmental Model Suggests Importance of Innate Immune Response to COVID-19 Infection in Rhesus Macaques.

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.

Introduction

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) (Yang et al. 2020).

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 (Li et al. 2020; 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. Pro-inflammatory 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 (Guan et al. 2020; 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, 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 (Yang et al. 2020). 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 (Shan et al. 2020).

Results

Overview of Modeling Structure

This novel model of the intrahost immune response to respiratory viral infection follows the time-evolution of 9 variables representing virus (VU and VL) and healthy (HU and HL) and infected (IU and IL) 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.

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

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 infection. Free virus is removed by one of three mechanisms: (1) internalization by epithelial cells (parameter γVH), (2) clearance by nonspecific immunity (parameters aVU and aVL), 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.

Fits to Experimental Data: Macaque Model of SARS-CoV-2 Infection

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 log10 TCID50/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 (Verbist et al. 2012; Verbist and Klonowski, 2012). 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 and .

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.

Fig. 2

Fits to experimental data from macaque SARS-CoV-2 infection. Line of best fit shown in purple generated with MATLAB data fitting function fminsearchbnd. Data points with standard deviations are represented by the black triangles with error bars. Virus in the nasal cavity and IL-15 (used to fit the Inflammation variable) were measured at 1,3, 5, 7, 10, 12, 14, 17, and 21 days post-infection. Virus in the lungs was measured only at 1, 3, and 5 days post infection, and clinical symptoms scores (used to fit the Damage variable) were measured each day from 0 to 21 days post-infection

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.

Explorations of Parameter Sensitivity

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

Predicted Effects of Aging on Host Immune Response Dynamics

Elderly patients have experienced more severe symptoms and outcomes than younger hosts with COVID-19 infection (Liu et al. 2020; 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 . 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 . Figure 3 shows the effect of changing both and 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.

Fig. 3

Sensitivity to age-related changes to the immune response. Parameters and 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 , the healthy cells will deplete and lead to death of the host. As Increases, the healthy cells can regenerate quickly enough to offset the viral load and transform the infection from lethal to survivable. High levels of can initiate a strong inflammatory response to eradicate the virus from the lower compartment even with low levels of . 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

As 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 is on the order of , 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 , 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 , 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.

Fig. 4

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

At intermediate values of , the stimulation of inflammation is high enough to contain the virus growth but also low enough to resolve within 60 dpi. As 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 , 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 , and .

The value of parameter 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 plane. When (red line), the infection clears and the host resolves the infection fully, as in our calibrated parameter set shown in Fig. 2. Raising 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 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).

Fig. 5

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

Spatial Disparities in Immune Response Dynamics

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 and 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 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.

Fig. 6

Sensitivity to adaptive immune responses. At high levels of , 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

Spatial differences also exist in the inflammatory response between the two compartments. Parameters and control the suppression of the viral growth by inflammation in the upper and lower respiratory tracts, respectively. High values of 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 and I_U subplots). Results of this parameter scan are demonstrated in Fig. 7.

Fig. 7

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

Drug Dosage Response

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.

The drug dose is simulated with a step function:

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

Anti-Inflammatory Drug Administration

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.

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

We performed parameter sweeps with , , and to determine the optimal calibration of the drug dosage. First, time of onset was varied between 1 and 5 days post-infection, with and 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 virus-driven 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.

Fig. 9

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 , as the duration increases, the infection will either remain chronic or heal fully. When is above , 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

We also explored the impact of changing the amount and duration of drug delivery while keeping time of onset constant. With dpi, we varied smax from 100 to 103 and 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 and 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 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.

Fig. 10

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

Antiviral Drug Administration

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:

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.

Fig. 11

Simulation of death due to excess viral load. Parameters in the model were adjusted to create a simulation in which viral titers remain elevated for a long period of time. This creates a lethal infection, as shown by the depletion of healthy cells in the upper compartment and high number of infected cells in that compartment. Parameter values used to simulate these trajectories are provided in Table S3

We then performed parameter sweeps with the drug dosage parameters to determine the appropriate calibration for the antiviral drug (Fig. 12). As in the anti-inflammatory simulations, we keep and and vary only 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 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.

Fig. 12

Effect of antiviral drug on simulated virus-induced death of the upper respiratory compartment. 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, falls below its threshold level and cannot regenerate, leading to death of the host. b Parameters and 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 or is too low, the system cannot recover

Discussion

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 pro-inflammatory 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 anti-viral 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 anti-inflammatory 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.

Methods

Experimental Data

Experimental data used for model calibration were downloaded from Munster et al. (Munster et al. 2020). In these experiments, 8 macaques were challenged with TCID50/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).

Model Structure and Equations

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 (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 (). 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).

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 aD. Anti-inflammatory mediators are created directly in response to pro-inflammatory signals, and they deactivate over time as the infection resolves.

Inflammation in the Absence of Virus: The H-F Subsystem

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: where and

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.

Fig. 13

Nullclines and fixed points for H-F subsystem. The H–F subsystem was obtained by making steady-state assumptions for all other variables. The H nullclines (green) and F nullclines (orange) denote the behavior of the H–F subsystem. The subsystem has three stable fixed points, circled on the graph: (1) at H = 0 and F = 0, (2) at H = 1 and F = 0, and (3) at nonzero H and F, representing chronic inflammation

Parameter Fitting

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 generated from the experimental conditions. Virus began at 4 × 105 TCID50/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.

Table 1

Biological interpretation of model parameters

Parameter	Meaning	Value
γ_V	Replication rate of the virus	3.459E9
κ_VFU	Inhibition of viral production in upper compartment by F	0.9044
γ_VH	Uptake of free virus by epithelial cells	3.1792
a_VU	Clearance of upper respiratory virus by nonspecific immunity	6.0624
a_V2	Saturation of nonspecific immunity	6.4517E-4
α_Vu	Clearance of upper respiratory virus by adaptive immunity	0.3077
β	Transport rate of virus from lower to upper compartment	5.4065E-6
w	Saturation of viral transport	8.5439E4
ε_V	Threshold for viral regeneration	1E-5
b_HD	Regeneration of healthy epithelial cells	4.3910E4
ε	Allee threshold for healthy cell regeneration	0.05
γ_HV	Rate of creation of infected cells	7.6059E-8
a_HF	Damage caused to healthy cells by inflammation	0.0099
a_IF	Damage caused to infected cells by inflammation	0.0126
a_I	Death of infected cells	4.4981
κ_VFL	Inhibition of viral production in upper compartment by F	0.7658
a_VL	Clearance of lower respiratory virus by nonspecific immunity	431.59
α_VL	Clearance of lower respiratory virus by adaptive immunity	0.0014
γ_FD	Stimulation of inflammation by presence of cellular damage	160.0381
b_FI	Stimulation of inflammation by presence of infected cells	0.0078
k_FG	Inhibition of inflammation by anti-inflammatory signal	0.3467
a_F	Decay rate of inflammation	10.2506
a	Inhibition of damage	3.193E-05
a_D	Decay/removal rate of damage	0.1503
γ_FG	Stimulation of anti-inflammatory signal by inflammation	0.1179
a_G	Decay rate of anti-inflammatory mediator	0.0835

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.

Table 2

Biological interpretation of drug dosage parameters

Parameter	Meaning	Bounds
s_max	Maximum concentration of drug dosage	100–103
t_on	Day of onset of drug delivery	0–14
t_dur	Number of days of drug delivery	0–30

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 13 kb)