A simple model of COVID-19 explains disease severity and the effect of treatments.

Considerable effort has been made to better understand why some people suffer from severe COVID-19 while others remain asymptomatic. This has led to important clinical findings; people with severe COVID-19 generally experience persistently high levels of inflammation, slower viral load decay, display a dysregulated type-I interferon response, have less active natural killer cells and increased levels of neutrophil extracellular traps. How these findings are connected to the pathogenesis of COVID-19 remains unclear. We propose a mathematical model that sheds light on this issue by focusing on cells that trigger inflammation through molecular patterns: infected cells carrying pathogen-associated molecular patterns (PAMPs) and damaged cells producing damage-associated molecular patterns (DAMPs). The former signals the presence of pathogens while the latter signals danger such as hypoxia or lack of nutrients. Analyses show that SARS-CoV-2 infections can lead to a self-perpetuating feedback loop between DAMP expressing cells and inflammation, identifying the inability to quickly clear PAMPs and DAMPs as the main contributor to hyperinflammation. The model explains clinical findings and reveal conditions that can increase the likelihood of desired clinical outcome from treatment administration. In particular, the analysis suggest that antivirals need to be administered early during infection to have an impact on disease severity. The simplicity of the model and its high level of consistency with clinical findings motivate its use for the formulation of new treatment strategies.

Introduction

COVID-19 symptoms severity differs wildly between infected individuals. Some individuals are infected without experiencing many of the characteristic symptoms, such as fever, coughs, body aches and the loss of taste or smell[1]. At the other end of the spectrum, a substantial minority will experience more extreme symptoms, such as acute respiratory distress syndrome and thrombotic complications that can lead to organ failure and death[1,2]. What distinguishes individuals experiencing more severe symptoms has been extensively studied[3]. These studies have identified a set of risk factors associated with severe COVID-19 such as older age, obesity, diabetes and past or present cancer[3]. COVID-19 severity likely depends on both the trajectory of the viral infection and the trajectory of the inflammatory response. An association was found between endothelial cell expression of angiotensin-converting enzyme 2 (ACE2), the receptor for severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) entry to host cells, and the presence of microthrombi in major organs such as the lung, heart, brain, and liver[4,5]. This suggests that the spread of the infection may be responsible for damage to vital tissues and organ dysfunction. It also underlines the necessity of an appropriate innate and adaptive immune response to limit the spread of the infection. To this effect, inflammation plays a crucial role by coordinating the immune response. However, elevated inflammatory markers in COVID-19 patients (IL-1β, IL-2R, IL-6, IL-8, IL-10, TNF-α, to name a few) have been associated with increased symptom severity, the need for ventilation and deaths[6-8]. This suggests that a sustained or exaggerated inflammatory response (hyperinflammation) may play an important role in determining the severity of disease outcomes. In particular, an inappropriate innate immune response was pointed out as a significant contributor to the hyperinflammatory state in COVID-19[9,10]. Despite the large effort by the scientific community, there are still many unknowns with regards to the mechanistic drivers of severe COVID-19. In turn, this knowledge gap hampers our ability to find new treatment strategies aiming to improve clinical outcomes.

Our main objective was to provide a simple quantitative framework to understand the pathogenesis of severe COVID-19 and to determine the importance of potential mechanisms. We aimed to provide a model that is simple and yet adequately captures the main clinical findings. Many within-host models of SARS-CoV-2 infection have been published[11-17] with a wide range of model complexities (from less than 10 to more than 80 parameters). These models mostly differ in terms of the complex interactions between the viral infection and the immune response that they include in their formulation (see Discussion for details). The large difference between model formulations suggests that identifying the key elements having an impact on clinical outcomes is a difficult task.

The model we formulated focuses on cells that trigger inflammation through molecular patterns: infected cells carrying pathogen-associated molecular patterns (PAMPs) and damaged cells producing damage-associated molecular patterns (DAMPs). We show that the clearance rate of infected and damaged cells by the innate immune response is of the utmost importance to reach a state of resolved inflammation. Our model can explain the following findings: (i) severe COVID-19 tends to be accompanied by hyperinflammation, (ii) those with severe COVID-19 generally experience a similar viral trajectory as mild cases, albeit with a slower viral load decay after the peak, (iii) the complex and conditional effect of antivirals and corticosteroids on disease severity, (iv) an inefficient type-I IFN response is associated with severe COVID-19, and (v) generation of bystander cell damage and infective removal of these cells are a critical component of severity. Note that this last point is reminiscent of clinical observations that, for example, less cytotoxic NK cells and higher levels of neutrophil extracellular traps (NETs) are associated with severe COVID-19. Overall, the simplicity of the model we propose along with its high level of consistency with clinical observations suggest it is an adequate framework for the study of COVID-19 pathogenesis and the effect of therapy.

Methods

The model

Our goal was to formulate a model that is simple enough to guide intuition, yet complex enough to allow relating with clinical outcomes. A schematic representation is provided in Fig. 1. The model is described below.

Figure 1

Schematic representation of the model described by Eqs. (1–7). Target cells () transition to a productively infected state () after successful infection by virions () at rate β. Virions are cleared at per capita rate c, while new virions are produced by infected cells at rate p. Target cells become refractory to infection () at rate (we assume target cells are exposed to a concentration of type-I IFN that is proportional to and that puts the cells into an antiviral state). Resting innate immune cells () become activated () at rate σ (I + J), where I and J are the number of infected and damaged bystander cells, respectively (we assume the extent of PAMP and DAMP signaling is proportional to +. Also, activated immune cells () die at per capita rate . Damaged bystander cells () are generated from an extensive proinflammatory response at a rate that is a Hill function of the number of activated immune cells . Infected cells die due to viral cytopathic effects at rate and damaged cells die from their injury at rate . The clearance of these cells also occurs by the action of activated innate immune cells at rate . The effect of the adaptive immune response is modeled by adding a constant term κ2 to the clearance of infected cells at time τ after infection. Finally, homeostatic processes allow replenishment of the population of resting cells at rate , where D00 is their homeostatic level.

Modeling the course of infection in the absence of immune responses

We describe the course of the infection in the absence of immune response through cells susceptible to viral infection () that can become infected () by contact with virions (). We assumed this process follows mass action with rate . Infected cells die at constant rate due to viral cytopathic effects. Infected cells produce virions at rate , which are cleared at per capita rate .

Modeling inflammation and the innate immune response

The model focuses on molecular patterns that initiate an inflammatory response. Pathogen associated molecular patterns (PAMPs) trigger an inflammatory response by signaling the presence of pathogens[18]. An inflammatory response can also be initiated in the absence of pathogenic infection through recognition of damage associated molecular patterns (DAMPs)[19]. Various conditions promote DAMP expression including hypoxia, low levels of glucose and amino acids, exposure to heat, physical stress or exposure to toxic molecular products[19-22]. Many of the same conditions that promote DAMP expression can be induced by inflammation itself, for example, an excessive presence of neutrophil extracellular traps (NETs) released during inflammation is linked to immunothrombosis[23]. This can contribute to hypoxia and nutrient deprivation that can lead to DAMP expression[24]. Overall, a positive feedback loop between inflammation and DAMP expression may be triggered[22,24]. Interestingly, patients with severe COVID-19 have increased levels of NETs[25,26]. NETs and platelet dysregulation may be specific to COVID-19 and both are associated with lung microthrombi, conditions that can ultimately favor an inflammatory response in the lungs[27].

The inflammatory response signals a need to eliminate infected or damaged cells and protect the rest of the organism from the perceived danger. Infected cells carrying PAMPs and uninfected cells producing DAMPs are represented by and in the model, respectively. In the organism, these molecular patterns promote the downstream recruitment and activation of numerous innate immune cells such as neutrophils, monocytes, macrophages and natural killer (NK) cells[28]. For simplicity, we lumped all cells that take part directly or indirectly (through cytokine signaling) in cytolytic or phagocytic activities into two model compartments, and representing resting and activated phenotypes, respectively. In our model, the rate at which cells become is proportional to the amount of PAMPs and DAMPs, themselves assumed proportional to the number of cells carrying these molecules, i.e. . We further assumed resting cells have a homeostatic level , that they maintain through the recruitment of resting immune cells over the course of infection at rate . We assumed cells are long-lived compared to the duration of acute infection, and that activated cells decay at rate . The role of in our model is to promote cytolytic and phagocytic activities to rid the system of PAMPs and DAMPs. Due to the innate nature of the modeled response, we assumed the effect to be similar on the decay of both I and J at per capita rate . This was partly motivated by the behavior of NK cells, which target injured epithelial cells through stimulation of receptor NKG2D, as well as infected cells having down-regulated major histocompatibility complex class 1 (MHC-I) molecules or upregulated MHC class I polypeptide–related sequence A (MICA) or sequence B (MICB)[29,30].

We assumed that high levels of inflammation can induce significant cellular stress and promote the expression of DAMPs in bystander cells. This is represented by the Hill term in Eq. (5) for the rate of generation of uninfected cells expressing DAMPs with being the maximum damage rate, being the concentration leading to half of the maximum damage rate, and is the Hill coefficient. We assumed damaged cells decay at rate < in the absence of an innate immune response.

Finally, we assumed an additional innate immune response, e.g., a type-I interferon response. Type-I interferons signal danger to neighboring cells, making the latter refractory to infection ()[31]. For simplicity, we approximated the amount of interferon in the microenvironment to be proportional to the number of infected cells[32]. Accordingly, we modeled the rate at which susceptible cells become refractory (R) by . Refractory cells expected to naturally revert to a susceptible state after 1/.

Modeling the adaptive immune response

The model mainly focuses on the innate response. To avoid unnecessary complexity, the adaptive immune response is represented by a single term that includes two parameters: , which is the maximum decay rate of infected cells due to the adaptive immune response, and τ, which represents the time post-infection that the adaptive response takes effect. This approach is similar to that in Pawelek et al.[33]. This response is thus modeled using the indicator function (t).

Structural analysis of the model

A substantial minority of individuals with COVID-19 experience a state of sustained high inflammation or hyperinflammation[9]. We evaluated whether the model allowed the existence of such a hyperinflammatory state. We used as a marker of inflammation. Through a bifurcation analysis, we searched for the existence of a stable steady state of resolved inflammation ( = 0) and the existence of a second stable steady state of hyperinflammation where is maintained well above 0. This type of analysis also allows identifying basins of attraction, i.e. regions in the space of variables that lead to specific inflammation trajectories over time. We used analytical techniques combined with the numerical bifurcation software Matcont[34], a Matlab software package designed for the analysis of equations such as Eqs. (1–7)[35,36]. Matcont uses continuation techniques to follow equilibria and performs normal form computations to classify bifurcation points.

Set of virtual markers

To analyze model behavior, a set of virtual markers were computed. These markers were used to compare model predictions to clinical observations and to investigate determinants of severe COVID-19. These markers are: (i) the peak viral load (peak VL, maximum value of V); (ii) the time of peak viral load; (iii) the difference between the viral load at its peak and the viral load 5 days after the peak; (iv) peak (used as a proxy measure for peak inflammation); (v) the time of peak ; (vi) at 60 days post-infection.

We further computed the hyperinflammation index, which takes a value of 1 if the value of log10 at the end of the simulated period was more than 99% of its peak value, and 0 otherwise. This resulted in the categorization of each simulation into inflammation trajectory groups: (i) resolved inflammation (R) or (ii) hyperinflammation (H). Finally, we computed the Disease Score, defined as the total number of cells (I + J) that died over the 60-day period post-infection. This score is meant to describe disease severity, with higher scores representing more severe COVID-19.

Identifying a space of realistic parameter values

The next step consisted in identifying a space of parameter values for which model predictions are consistent with a minimal set of clinical observations to allow in silico investigation of the model. In particular, acceptable parameter values should result in: (i) peak viral loads (VL) values between 4 and 10 log10, (ii) peak viral loads achieved 2 to 14 days following infection, and (iii) an innate immune response during the course of infection (the activation of at least 1% of all resting innate immune cells). These are referred to as conditions (i)–(iii). Condition (i) was chosen to represent peak viral loads observed using nasopharyngeal swabs in a population of individuals infected by the virus[7,37]. Condition (ii) was chosen based on reports of peak viral loads occurring around the time of symptom onset and symptom onset primarily occurring within 14 days of infection (median 4–5 days)[14,16,32,38]. Finally, we chose condition (iii) to ensure reaches high enough values for an observable effect of on either or .

The determination of the space of parameter values that satisfy the above conditions was done iteratively. First, we established bounds of values for each parameter using literature estimates when available, results from the bifurcation analysis to ensure that hyperinflammation was achievable, and preliminary simulations when no estimate was available (see Table 1). We then selected n = 100,000 vectors of parameters from the parameter space using a Latin hypercube approach to minimize the chances of unexplored multidimensional subspaces[39]. Simulations were performed for each selected vector of parameters using initial conditions listed in Table 2 to predict infection and inflammation trajectories from day 0 (day of infection) to day 60. We identified regions in the parameter subspace leading to unacceptable results based on the conditions (i)-(iii) listed above. We subsequently refined the parameter space and repeated the procedure until we were satisfied that a randomly selected vector of parameters would lead to a high likelihood of satisfying acceptance criteria (i) to (iii) (> 70% acceptance probability). Table 1 describes all parameters, the resulting space of parameter values along with references from the literature when applicable.

Table 1

Model parameters, the range of explored values used for the in silico investigation along with justifications and references.

Par	Description (unit)	Range of explored values	Justification and References
β	Infectivity rate (virion−1 day−1)	[− 12;− 6] (log10 scale)	β values of 1.9E-6 and 6.6E-7 were used in Ke et al. for the upper and lower respiratory tract, respectively[40]. A value of 5.2 × 10–6 was used in Kim et al.[14]. A larger range of values was used to account for variability between individuals and uncertainty in parameter values due to a possible correlation with parameter p
p	Virion production rate (virions cell–1 mL−1 day−1)	[–1;4] (log10 scale)	p values of 51.4 per swab per day and 0.35 per mL per day were used in Ke et al. for the upper and lower respiratory tract, respectively[40]. A larger range of values was used to account for variability between individuals and uncertainty in parameter values due to a possible correlation with parameter β. Accordingly, it was inefficient to sample β and p independently. It was determined through linear regression that a value of log_10p determined by  6.8–1.0*log_10β + U, where U is a uniform random variable with range [ 0.4;0.5], greatly enhanced the likelihood of peak viral loads being between 4 and 10 log10 and occurring between day 2 and 14 post-infection .
δ_I	Death rate from viral cytopathic effects (day–1)	[0.05;0.1]	Jenner et al. used an infected cell death rate of of 0.014 day–1 in the preprint version of the paper [11]. The lower bound was reviewed to allow a greater probability of satisfying the simulation acceptance conditions (see Methods). A study demonstrated that viral production from infected cells was maintained at high levels for up to 6 days in vitro[41]
c	Virion clearance rate (day–1)	[10;30]	Values explored in Goncalves et al. as well as Ke et al. were between 5 and 20 day–1[32,37]
κ_0	Rate of transition to IFN induced refractory state (cell–1 day–1)	[− 8;− 5] (log10 scale)	A value of 1.3 × 10–6 was estimated in Ke et al.[32]
δ_R	Refractory state reversion rate (day–1)	[− 4;− 2] (log10 scale)	A value of 0.0044 day–1 was estimated in Ke et al.[40]
τ	Time delay of adaptive immune response post-infection (days)	[7;40]	The range was chosen to match the variability between individuals in time of viral clearance post-infection, with a median of around 25 days[7,40]
κ_2	Effect of adaptive immune response on the clearance of infected cells (day–1)	[2;6]	This range was chosen to ensure viral clearance is achieved shortly after τ days
m	Hill coefficient for bystander cell damage	3	A value of 3 ensures a steep progression of the damage rate as a function of D_1. In other words, we assumed that the damage due to inflammation is not substantial unless inflammation reaches high levels
ν	Maximum bystander damage rate (cells day–1)	[7.5;8.5] (log10 scale)	An initial range of values was chosen based on bifurcation analyses. The range was refined to ensure that around 5–20% of simulations led to hyperinflammation
K	D_1 concentration leading to half the bystander damage rate (cell)	[1.5 × 106;2 × 106]	An initial range of values was chosen based on bifurcation analyses. The range was refined to ensure that around 5–20% of simulations lead to hyperinflammation
δ_J	Death rate from cell damage	[0.01;0.05]	Chosen such that the death rate for damaged cells is slower than for infected cells
D_00	Resting immune cells homeostatic constant (cells mL–1)	106	A value of 106 was used for alveolar macrophages in Smith et al.[42]. Values ranging from 4 × 105 and 3 × 107 were used for various populations of innate immune cells in Jenner et al.[11]. A single value was used since we were not interested in investigating the impact of variability in D_00 on the severity of COVID-19
λ	Resting immune cells replenishing rate or recruitment rate (day−1)	[− 2;1] (log10 scale)	A comparative value of 0.22 was used in Jenner et al. for monocytes[11]
σ	Innate immune cell activation rate (cell–1 day–1)	[− 8.5; − 7.5] (log10 scale)	A value of the order of 1 × 10–6 was used in Jenner et al.[11]. The range of values was adjusted to ensure that at least 1% of all resting innate immune cells activated over the course of infection and that high inflammation levels (D_1 > 1 × 106) did not occur prior to peak VL
δ_D	Activated innate immune cell average death rate (day–1)	[0.1;0.3]	A value of 0.3 was used for activated macrophages in Jenner et al.[11]. A value of 0.04 was used for activated macrophages in Pawelek et al.[43].Natural killer cells are shown to have a turnover of around 2 weeks[44]. In Sadria et al., a value of 0.2 was used to represent the natural death rate of effector cell[13]
κ_1	Effect of innate immune response on the clearance of PAMP and DAMP expressing cells (cell–1 day1)	[− 8;− 4] (log10 scale)	An initial range of values was chosen based on bifurcation analyses. The range was refined to ensure that around 5–20% of simulations lead to hyperinflammation

Par Parameters.

Table 2

Initial conditions for all simulations.

Variable	Description	Initial value
T	Target cells	4.8 × 108 cells[40]
R	Refractory cells	0 cells
I	Infected cells	10 cells
V	Virions	0 virions
J	Damaged cells	0 cells
D_0	Resting innate immune cells	D_00 (see Table 1)
D_1	Activate innate immune cells	0 cells

Simulating clinical observations

We sampled a large number of sets of parameter values (n = 1,000,000) from the space of parameter values described in Table 1 using a Latin hypercube approach[39]. Accordingly, the marginal distribution of values being sampled for each parameter was uniform across the considered range (Table 1). For each selected vector, simulations were performed using initial conditions listed in Table 2 to predict infection and inflammation trajectories from day 0 (day of infection) to day 60. Analyses were only performed on the simulations satisfying the acceptance criteria (conditions (i)–(iii)). We first investigated the distribution of parameter values and virtual markers using histograms. We used violin plots to study bivariate associations between parameters and inflammation trajectory groups (resolved inflammation or hyperinflammation). Multivariate analyses were performed using decision trees[45]. We used the hyperinflammation index as the variable being predicted by model parameters. We first obtained a single tree. Cross-validation and deviance plots were used to guide the choice of the optimal tree[45]. To account for uncertainty around the formulation of a single tree, we also performed a Random Forest analysis, generating 100 trees, and reported the mean decrease in GINI index (a measure of decreased node impurity from choosing a parameter for tree splits, i.e. its ability to discriminate inflammation trajectory groups)[46].

Simulating the effect of treatments on COVID-19 severity

Finally, we evaluated if the model could replicate clinical findings regarding the treatment of COVID-19. For this purpose, we selected 10,000 simulations leading to lower Disease Severity scores and the same number of simulations leading to higher scores from the sample of accepted simulations. To allow comparison with clinical data, we assumed the former group represents mild/moderate disease, while the latter represents severe COVID-19. To simulate the use of corticosteroids, we assumed a reduction of by 50% for a period of 10 days. To simulate the use of potent antivirals, we decreased parameters or to 1% of their original value at peak viral load for the remainder of the simulation. For each simulation, we modified parameter values at the time of peak viral load as it is estimated that the peak is reached within the few days following symptom onset[14,32,38]. We also simulated treatment administration a day prior to the expected viral load peak to study the effect of early treatment. The difference in log10 Disease Score (with vs without treatment) was computed and reported.

Results

The model allows two stable steady states (Hyperinflammation and resolved inflammation) under realistic parameter values

To understand the general dynamics of the model given by Eqs. (1)–(7) and the types of infection outcomes predicted by the model, we first performed a bifurcation analysis using baseline parameter values (see Supplementary Material for detail). Interestingly, the analysis shows that under certain biologically plausible parameter values, there exists bistability in the system with both a stable hyperinflammatory state and a stable resolved inflammatory state. In more technical terms, there are three equilibria in the model: a stable high-inflammation state corresponding to hyperinflammation, an unstable equilibrium with non-zero, but low, inflammation and another stable equilibrium corresponding to resolved infection and resolved inflammation. The first two steady states appear/disappear following a saddle-node bifurcation as shown in the bifurcation plots (Figs. S1 and S2).

The bifurcation analysis identified important parameters that dictate the existence of a hyperinflammatory steady state: (i) , which represents the effect of the innate response ( on the clearance of cells carrying PAMPs and DAMPs, and (ii) the parameter that dictates the amount of bystander cell damage resulting from inflammation (Fig. S1). Our analysis shows there is a threshold value of beyond which the hyperinflammatory state ceases to exist (Fig. S1A). Further, hyperinflammation can only occur if there is enough immune driven inflammation beyond a threshold value of (Figs. S1B and S2B).

When bistability exists in the system, inflammation trajectories could either converge to one stable state or the other over time. What determines the long-term inflammation trajectory is the amount of bystander cell damage caused by inflammation over the course of the infection. Mathematically, this is represented by a saddle-node bifurcation. The unstable lower branch of equilibria acts as a separatrix between the hyperinflammatory and resolved-inflammation states and so defines the basin of attraction of the hyperinflammatory state (see Fig. S1). Thus, hyperinflammation can only occur if the infection induced inflammation is severe enough to force the system across this separatrix and into this basin.

Overall, these results suggest that the clearance of cells carrying PAMPs and DAMPs and the amount of bystander cell damage from inflammation are important determinants of disease outcomes. When the clearance rate of infected and damaged cells is sufficiently high or bystander cell damage is low, infection leads to non-severe outcomes.

Model simulations of a virtual cohort of infected individuals

We next simulated the model by sampling n = 1,000,000 sets of parameter values across biologically plausible ranges (Methods). We use these simulation results to analyze the different viral load and inflammatory response trajectories in the population, such that key determinants of disease outcomes can be further identified. In our analysis, we only included simulations satisfying the acceptance criteria that are consistent with broad patterns seen in clinical studies (described in the Method section as conditions i-iii) (accepted simulations, n = 739,465 of 1,000,000 or 73.9%). The distribution of parameter values that led to accepted simulations are presented in Supplementary Material (Fig. S3). Most parameters preserved their sampling distributions, i.e., the distribution remained uniform across the considered ranges after discarding simulations that do not satisfy the acceptance criteria. Further, there was little correlation between almost all parameters. However, there was a strong correlation between and ( 0.99). Consequently, parameters and could not be independently sampled in order to output acceptable simulations (Fig. S4). Table 1 describes the sampling strategy used to account for this dependency.

A strong association between hyperinflammation and the disease score

To represent disease severity, we define the Disease score as the total number of infected and bystander cells and that died over a 60-day period. Simulations that give rise to the Disease Score distribution are illustrated in Fig. 2. They are categorized in groups by inflammation trajectory: (i) the inflammation marker decreased after peak inflammation or (ii) the inflammation marker kept increasing or was maintained at high level following peak inflammation (hyperinflammation). There was a direct link between the presence of hyperinflammation and the Disease Score (see Fig. 2a). Overall, 13.5% of accepted simulations exhibited hyperinflammation. Hyperinflammation and a high Disease Score were both associated with a higher number of cells that died following injury from inflammation (Fig. 2b).

Figure 2

(a) Distribution of Disease Scores and (b) distribution of the total number of bystander cells that died by inflammation trajectory groups: Resolved inflammation (orange), Hyperinflammation (pink).

Similar viral load but divergent inflammation trajectories

Viral load and inflammation trajectories over time by inflammation groups are presented in Fig. 3. Simulations resulting in both resolved inflammation or hyperinflammation exhibited similar viral load dynamics. The distribution of peak viral load and the time to reach this peak after infection largely overlapped between groups. However, those with resolved inflammation tended to have a faster VL decay after peak VL (Fig. 3c). The most remarkable difference between groups pertained to the dynamics of the inflammation marker . Peak levels of inflammation were higher for simulations resulting in hyperinflammation. For those with resolved inflammation, peak inflammation was generally observed around the time of peak viral load while for those with hyperinflammation, peak inflammation was observed much later (Fig. 3h).

Figure 3

Viral load and inflammation trajectory characteristics by inflammation trajectory groups. (a) Viral loads over the course of infection. The shaded area corresponds to the 10th and 90th percentiles of the viral loads, while the curve represents the median. (b) The distribution of peak viral loads, (c) the VL decay from peak infection to 5 days after peak infection and (d) the time of occurrence of peak VL after infection. (e) Inflammation trajectories by inflammation trajectory groups. The shaded area corresponds to the 10th and 90th percentiles of , while the curve represents the median. (f) the distribution of peak , (g) the distribution of at 60 days post-infection and (h) the time of occurrence of peak  after infection. In orange and represented by the symbol R, Resolved inflammation. In pink and represented by the symbol H, Hyperinflammation.

Association between hyperinflammation and characteristics of the innate immune response

Next, we compared the distribution of parameter values between inflammation trajectory groups. Figure 4 shows violin plots for the 15 parameters that were allowed to vary between simulations. The most striking differences between groups are observed for parameters and . Lower values of (type-I IFN secretion and/or response) and (cytolytic and phagocytic activities of innate immune cells) were associated with a greater risk of hyperinflammation and a high Disease Score. Similarly, higher values of (innate immune cell recruitment rate), (activation rate), (survivability of activated innate immune cells) and (damage rate due to inflammation) were associated with hyperinflammation and a high Disease Score.

Figure 4

Distribution of the model parameters by inflammation trajectory groups. Distribution overlap may be discriminated from multivariate models. Inflammation trajectory groups R: Resolved inflammation, H: Hyperinflammation.

Prediction of the risk of hyperinflammation from characteristics of the innate immune response

We used regression tree analysis to reveal the discriminatory importance of parameters from a multivariate perspective. The regression trees attempted to discriminate simulations leading to hyperinflammation from those leading to a resolved inflammatory state. Results are shown in Fig. 5.

Figure 5

Regression tree analysis results. (a) Single optimal tree for the prediction of hyperinflammation from model parameters. The tree reads from left to right. At each labeled node, simulations either go up if the value for the associated parameter is higher than a threshold determined by the procedure (threshold not shown, see Supplementary Figure S7), or down otherwise. Branch length represents the amount of classification error explained by the node. At each terminal node, the percentage of simulations as well as the risk of hyperinflammation within members of the node are reported. (b) Parameter importance based on the GINI index. A greater mean GINI decrease indicates a parameter that is more discriminatory.

The tree has a root (left side), branches, nodes (where branches separate) and terminal leaves (right side). Simulations enter at the root and separate at nodes. At each node, one parameter and one threshold value were chosen by the regression algorithm based on their ability to separate simulations into more homogeneous groups in terms of inflammation trajectory. Simulations are directed toward the lower branch if the value for the chosen parameter for the simulation is smaller than the threshold, or toward the upper branch otherwise. This resulted in 5 terminal leaves that partitioned all simulations into 5 groups. Longer branches following a node signify that the node allowed a better separation of the two inflammation trajectory groups.

Figure 5 illustrates the multivariate conditions favoring hyperinflammation. The most important parameters pertain to type-I IFN response, the ability of the system to clear cells carrying PAMPs and DAMPs ( and ) and the immune cell activation and recruitment rate ( and ). In particular, combined low values of and and high values of and led to a dramatically increased chance of hyperinflammation (73.2% risk of hyperinflammation).

Differences in the effect of treatments in terms of treatment type, administration time and predicted disease severity

To validate and demonstrate the utility of our model, we simulated the use of corticosteroids and antivirals in infected individuals and compared the model results with clinical findings. First, we modeled corticosteroid treatment by assuming the treatment leads to a reduction of by 50% for a period of 10 days following peak infection. The in silico administration of corticosteroids had a remarkably different effect on the Disease Score depending on the inflammation trajectory groups (Fig. 6). Among those simulations where resolved inflammation is predicted in the absence of treatment (orange plot), corticosteroids were often detrimental (32% chance of an increase in Disease Score). Such a detrimental effect was not observed among those for which hyperinflammation was predicted in absence of treatment (pink plot). Greater improvements were observed in the latter group: 23% had a greater than 0.5 log10 decrease in Disease Score, compared to only 0.7% in the former group. There were only small differences in the effect of corticosteroids across the investigated timing of drug administration (at peak viral load in Fig. 6a and a day prior in Fig. 6b). Corticosteroid treatment also had an effect on the slope of viral load decay, which was substantial for those simulations of mild/moderate disease and lesser among simulations of severe COVID-19 (Fig. S8).

Figure 6

Violin plots of the effect of virtual treatment on the Disease Score. (a) represents administration of corticosteroids while (b) represents antiviral drug administration. Negative values represent improvements while positive values represent the worsening of symptoms. Note there were no clear difference between a reduction in or in the simulation of antivirals so the figure applies to both cases. Orange denotes the effect of treatment among individuals who would have resolved inflammation in the absence of treatment, whereas pink denotes the effect of treatment among individuals who would have had hyperinflammation in the absence of treatment.

Simulations of the effect of antivirals were performed by decreasing parameters or to 1% of their original value at peak viral load. Modifying or had a very similar effect on Disease Scores None of the simulations revealed a substantial increase in Disease Score. However, the impact of antivirals was very different between inflammation trajectory groups and between times of drug administration. Earlier drug administration led to large improvements in Disease Score for both groups (Fig. 6b). However, administration at peak viral load led to much reduced improvements in Disease Score (Fig. 6a), in particular for those simulations leading to hyperinflammation in absence of treatment. In this case, infection had already driven inflammation to high levels by the time peak viral load was reached (although peak inflammation was not reached until much later), which translated into a number of damaged cells that corresponded to being within the basin of attraction of hyperinflammation (Fig. S1).

Discussion

In this work, a mathematical model was formulated to represent the within-host dynamics of COVID-19 infection and inflammation. The objective was to provide a quantitative explanation for the range of COVID-19 symptom severity among individuals and to reveal the discriminatory importance of modeled mechanisms. The hypothesis we explored was that high levels of inflammation in COVID-19 may produce a significant amount of damage to uninfected cells. These cells would then produce DAMPs that further stimulate the inflammatory response. This model produced predictions that are consistent with clinical observations.

Hyperinflammation and disease severity

The model predicts that those with higher Disease Scores had substantially higher levels of inflammation (Fig. 2a). Further, peak inflammation was not reached until much later in those simulations leading to hyperinflammation (Fig. 3h). These modeling results are consistent with clinical findings that non-survivors tend to experience hyperinflammation and increasing levels of inflammatory biomarkers up to the time of death[7,9,47]. Our analysis provides an explanation of these findings. In particular, the numerical bifurcation analysis suggests that viral infection can push the immune system into a self-sustaining high inflammatory state that can persist well beyond the resolution of infection. This hyperinflammatory state causes additional damage to uninfected cells, consequently leading to much higher Disease Scores. Due to the strong association between hyperinflammation and the Disease Score, we used those simulations that led to hyperinflammation as an in silico description of severe COVID-19. The percentage of simulations exhibiting hyperinflammation (13.5%) closely matched the reported proportion of the infected population with severe COVID-19 (14%), further motivating this decision[1].

Viral load dynamics and disease severity

In terms of viral load, the main difference that was reported between severe COVID-19 and those with milder disease pertained to the slope of VL after disease onset[7]. In particular, a faster VL decay was observed for those with milder disease after symptom onset, as measured using nasopharyngeal swabs[7]. This finding was also reported for VL sampled using other means and from various physiological compartments[48]. Interestingly, disease severity does not correlate strongly with peak viral loads[7,48]. Our simulation results are consistent with these findings (Fig. 3b and c).The model offers an explanation to the slower viral load decay after peak infection in cases of severe COVID-19 through parameter a lower value both results in slower clearance of productively infected cells (see Eq. 3) and increases the risk of hyperinflammation (see Fig. 5). In fact, having a low value was the most important predictor of hyperinflammation (Fig. 5b). In the model, low values lead to the persistence of and cells, thereby prolonging PAMP and DAMP signaling and its downstream impact on the inflammatory response. The heightened inflammatory response promotes the generation of damaged cells , strengthening DAMP signaling. Our analysis suggests that the ensuing feedback loop is the hallmark of hyperinflammation and severe COVID-19.

Characteristics of the innate immune response and disease severity

Many clinically observed associations were found between components of the innate immune response and disease severity. Among them, a dysregulated IFN response has been repeatedly associated with severe COVID-19[10,49]. Comparatively, our model suggests that a weaker IFN antiviral response (lower leads to a higher likelihood of hyperinflammation (see Figs. 4c, 5a,b). One of the roles of type-I IFN is to limit the number of target cells that can be infected[31]. Although it did not have a big impact on peak VL or the time to reach this peak in our simulations, an efficient type-I IFN response (higher values) had a significant impact on inflammation. IFN by limiting the rate and the shear number of cells that carry PAMPs also constrains inflammation.

Poor NK-cell cytotoxic ability was also linked with severe COVID-19[50]. NK-cells are important actors in the innate immune response that can target both infected and damaged cells and release molecules that precipitate their apoptosis[29,30]. In the model, poorer clearance of cells having PAMPs and DAMPs by the innate immune response is represented by lower values, the most important predictor of hyperinflammation and severe COVID-19 (Fig. 5b). A poorer response from NK cells could hence lead to severe COVID-19 by enabling prolonged PAMP and DAMP signaling. Patients with severe COVID-19 also have increased levels of neutrophil extracellular traps (NETs)[25,26]. These can lead to immunothrombosis and the generation of damaged cells[23,24]. In the model, higher generation of damaged cells is represented by larger values of the parameter (Eq. 5). The bifurcation analysis suggests that larger values of lowers the threshold for the number of damaged cells that are required to reach a hyperinflammatory state (Fig. S1), facilitating severe COVID-19.

Finally, clinically defined severe COVID-19 has been associated with higher abundance and activation of proinflammatory macrophages[51]. In the model, those with higher disease scores had a higher number of innate immune cells () and a higher proportion of these cells were activated. The parameters that dictate innate immune cell recruitment () and activation () were both strongly associated with hyperinflammation and having a high Disease Score.

Effect of corticosteroids and antivirals on disease severity

We evaluated if the model could replicate clinical findings regarding the treatment of COVID-19. We used the model to simulate the treatment of corticosteroid, such as dexamethasone, or an effective antiviral. The model predictions are consistent with many clinical findings, suggesting that the model we developed here will be a useful tool to understand SARS-CoV-2 pathogenesis and predicting the impact of treatment. Indeed, observational studies revealed the use of corticosteroids can lead to more severe symptoms in those with milder disease, but generally improves outcomes for those with more severe symptoms[51-55]. Further, it is reported that those with milder COVID-19 generally experience slower viral load decay under corticosteroid treatment, an effect that was not found to be statistically significant among those with severe disease[52]. This latter result was also observed in our simulations (Fig. S8). One of the roles of inflammation is to stimulate cytolytic and phagocytic activities. By lowering this ability among those who experienced milder disease, the use of corticosteroids may lead to slower viral clearance. Hence, corticosteroids could have both a negative effect (slower clearance of PAMP carrying cells) and a beneficial effect (slower rate of bystander cell damage) on disease pathogenesis in those with mild disease. This could result in corticosteroids sometimes improving, sometimes worsening disease severity. Comparatively, the net beneficial effect of corticosteroids on those with severe COVID-19 may be the result of a smaller downstream impact of the drug on already less efficient NK cells[50].

For antivirals, simulations suggest that early administration is crucial for antivirals to impact disease severity, particularly for those that would have experienced severe COVID-19 in absence of treatment (Fig. 6). Our results also suggest the existence of an inflammation threshold beyond which antivirals may be unable to prevent hyperinflammation. It also suggests that early administration could reduce disease severity by preventing the infection from driving inflammation across the threshold. Comparatively, many clinical trials failed to show a substantial effect of antivirals[56]. More recently, a reduction of around 50% in the risk of hospitalization was observed after administering the antiviral molnupiravir[57]. The studied cohort consisted of individuals that had mild/moderate symptoms and were not expected to be hospitalized within 48 h of randomization[57]. Our results suggest this latter criteria may be crucial to ensure the beneficial effect of antivirals on disease severity.

Fitting clinical markers of disease severity

We used the disease score as a proxy measure of disease severity in our model. Quantitatively fitting our model by comparing the disease score with usual clinical disease severity scores may be inappropriate as it would require an established relationship between these quantities. In addition, clinical score data may be too sparse to reliably estimate the parameters in our model. Instead, biologically plausible parameter ranges were used to investigate gross changes in disease severity. Accordingly, changes that give rise to sufficiently different distributions of disease severity, such as those effected by some variants of concern, may inform parameters characterizing the variants. For example, the delta variant may cause overall higher viral loads[58], which would lead to stronger immune activation on average, explaining the higher severity of this variant.

Other within-host models of SARS-CoV-2 infections

The model we propose is unique, but many of the effects it includes are present in other within-host models of SARS-CoV-2 infection. Despite differences in model formulations, there is a general agreement across models about the necessity for early antiviral administration[12-15]. However, models differed in terms of the components of the inflammatory, innate or adaptive immune response they include[11-17]. Some of the models included an effect of type-I interferon, either as a promoter of the pro-inflammatory response or helping susceptible cells resist infection[11-13]. One of these models concluded, as per our analysis, that an inefficient type-I IFN response can lead to accentuated tissue damage[11]. In comparison, our model simultaneously explains the more important clinical associations between severe disease and biomarkers of the innate immune response. It distinguishes itself by the inclusion of both a positive feedback loop between damaged cells and inflammation and an effect of innate immune cells on both infected and damaged cells. This latter effect highlights the importance of cells capable of clearing both types of cells in the pathogenesis of severe COVID-19.

Clinically, it has been shown that many individual characteristics, such as age and immune markers, e.g. lymphocyte counts and neutrophil counts, are associated with COVID-19 disease severity[59].This complexity contrasts with the simplicity of the model we propose. While it is likely that many factors are involved in severe disease development, our model represents a high-level description of some of the physiological pathways, in particular, the innate immune response pathways, that may be involved in disease progression. This is consistent with the data from Li et al. and Lasso et al. that many immune markers predict disease severity better than other features, such as demographic characteristics[59,60]. Although we investigated more complex models, we decided against the modeling of individual cytokines or cells, as they often exhibit overlapping functions, and because some of these functions have been poorly studied leading to uncertainties in parameter values. However, some of the more complex models reported in the literature did give rise to interesting hypotheses that may warrant further investigation, such as the role of monocyte-to-macrophage differentiation or the role of anti-inflammatory cytokines[11,12]. Further, more complex models have provided complementary insight on the conditional impact of immunomodulatory treatment on disease severity that are consistent with the results reported herein[61,62].

Conclusion

Our analysis revealed key aspects of the innate immune response that dictate inflammation trajectories and disease severity. The most important parameter suggested by bifurcation and decision tree analyses was , representing the ability of the system to rid itself of cells carrying PAMPs and DAMPs. The analysis suggested that when this parameter is high enough, hyperinflammation can be avoided. The bifurcation analysis also suggested that small values of , representing the amount of bystander cell damage due to inflammation has a similar effect. In other words, the ability of the innate immune response to target PAMPs and DAMPs carrying cells appears to be key to the determination of COVID-19 severity. This suggests that therapies that specifically target aspects of the innate immune response may prove beneficial in comparison to broadly acting anti-inflammatory agents. The model also underlined the role of DAMPs in maintaining high levels of inflammation later in the course of infection. It suggests DAMPs may be an interesting therapeutic target for COVID-19. When exploring such novel treatment strategies, the model presented here could provide a means of exploring timing of treatment and dose effects in silico. Hopefully, a better understanding of the pathology of SARS-CoV-2 will lead to decreased mortality in this and similar diseases.

Supplementary Information.