key: cord-0881577-7f4idbzi
authors: Sanche, S.; Cassidy, T.; Chu, P.; Perelson, A. S.; Ribeiro, R. M.; Ke, R.
title: A simple model of COVID-19 explains disease severity and the effect of treatments
date: 2021-12-02
journal: nan
DOI: 10.1101/2021.11.29.21267028
sha: f30ab2d51ab27e90a17fd597473b12f4f43331ca
doc_id: 881577
cord_uid: 7f4idbzi

Considerable effort was 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. The model 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). The former signals the presence of pathogens while the latter signals danger such as hypoxia or the lack of nutrients. Analyses show that SARS-CoV-2 infections can lead to a self-perpetuating feedback loop between DAMP expressing cells and inflammation. It identifies the inability to quickly clear PAMPs and DAMPs as the main contributor to hyperinflammation. The model explains clinical findings and the conditional impact of treatments 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.

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. 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 severity of symptoms, the need for ventilation and deaths. (6) (7) (8) This suggests that a sustained or exaggerated inflammatory response (hyperinflammation) may play an important role in determining disease severity. 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 provided 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) (12) (13) (14) (15) (16) (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 damage cells by the innate immune response is of the utmost importance to reach a state of resolved of inflammation. Despite its simplicity, 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.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint

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.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint 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 0 (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 ( 0 ) become activated ( 1 ) 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 ( 1 ) 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 1 . 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 1 1 . 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 ( 0 ) at rate λ( 00 − 0 ), where D00 is their homeostatic level.

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 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint and activated phenotypes, respectively. In our model, the rate at which 0 cells become 1 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 0 have a homeostatic level 00 , that they maintain through the recruitment of resting immune cells over the course of infection at rate . We assumed 0 cells are long-lived compared to the duration of acute infection, and that activated cells decay at rate . The role of 1 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 1 1 . 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 downregulated 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 1 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 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint

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: 2 , 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 modelled using the indicator function > (t). 

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 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint after the peak; iv) peak 1 (used as a proxy measure for peak inflammation); v) the time of peak 1 ; vi) 1 at 60 days post-infection.

We further computed the hyperinflammation index, which takes a value of 1 if the value of log10( 1 ) 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- 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 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021.

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.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; Table 1 A study demonstrated that viral production from infected cells was maintained at high levels for up to 6 days in vitro. ( An initial range of values was chosen based on bifurcation analyses. The range was refined to ensure that the around 5-20% of simulations lead to hyperinflammation.

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. Crossvalidation 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) Resting innate immune cells 00 (see Table 1) 1 Activate innate immune cells 0 cells . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint

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

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 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint steady states appear/disappear following a saddle-node bifurcation as shown in the bifurcation plots (Fig S1 and S2) .

The bifurcation analysis identified important parameters that dictate the existence of a hyperinflammatory steady state: i) 1 , which represents the effect of the innate response ( 1 ) 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 1 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 ( Fig S1B and Fig 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.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint

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.

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 give rise to the Disease Score distribution are illustrated in Fig. 2 . They are categorized in groups by inflammation trajectory: i) the inflammation marker 1 decreased after peak inflammation or ii) the inflammation marker 1 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) .

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint

Resolved inflammation (orange), Hyperinflammation (pink).

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

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint 

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 1 , , 0 , , and

. Lower values of 0 (type-I IFN secretion and/or response) and 1 (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), . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint 20 (activation rate), 1/ (survivability of activated innate immune cells) and (damage rate due to inflammation) were associated with hyperinflammation and a high Disease Score.

Distribution overlap may be discriminated from multivariate models. Inflammation trajectory groups R: Resolved inflammation, H: Hyperinflammation.

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 .

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint 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. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint and DAMPs ( 0 and 1 ) and the immune cell activation and recruitment rate ( and ). In particular, combined low values of 0 and 1 and high values of and led to a dramatically increased chance of hyperinflammation (73.2% risk of hyperinflammation).

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

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 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint 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 is 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) .

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint

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.

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 revealed 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)

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org /10.1101 /10. /2021 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 3C).The model offers an explanation to the slower viral load decay after peak infection in cases of severe COVID-19 through parameter 1 ; a lower 1 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 1 value was the most important predictor of hyperinflammation 

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 response (lower 0 ) leads to a higher likelihood of hyperinflammation (see Fig 4C, 5A and 5B ).

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 0 values) had a significant impact on inflammation. IFN by limiting the rate and shear number of cells that carry PAMPs also constrains inflammation.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org /10.1101 /10. /2021 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 1 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 ( 0 + 1 ) 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.

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 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint improves outcomes for those with more severe symptoms. (51, (53) (54) (55) (56) 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.(53) 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.(57) More recently, a reduction of around 50% in the risk of hospitalization was observed after administering the antiviral molnupiravir. (58) The studied cohort consisted of individuals that had mild/moderate symptoms and were not expected to be hospitalized within 48 hours of randomization.(58) Our results suggest this latter criteria may be crucial to ensure the beneficial effect of antivirals on disease severity.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint

The model we propose is unique, but many of the effects it includes are present in other withinhost 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) (13) (14) (15) However, models differed in terms of the components of the inflammatory, innate or adaptive immune response they include. (11) (12) (13) (14) (15) (16) (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) (12) (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. 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 antiinflammatory cytokines. (11, 12) 

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 1 , 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 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint , 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. 

The authors have declared that no competing interests exist.

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 2, 2021. ; https://doi.org/10.1101/2021.11.29.21267028 doi: medRxiv preprint

Mild or moderate Covid-19

Thrombotic complications of COVID-19

Risk factors for severe and critically ill COVID-19 patients: a review

Vander Heide RS. Pulmonary and cardiac pathology in African American patients with COVID-19: an autopsy series from New Orleans. The Lancet Respiratory Medicine

Pathophysiology of SARS-CoV-2: the Mount Sinai COVID-19 autopsy experience

An inflammatory cytokine signature predicts COVID-19 severity and survival

Longitudinal analyses reveal immunological misfiring in severe COVID-19

Progress on role of cytokine storm in exacerbation of coronavirus disease 2019 (COVID-19)

Immunopathology of Hyperinflammation in COVID-19

Imbalanced host response to SARS-CoV-2 drives development of COVID-19

COVID-19 virtual patient cohort suggests immune mechanisms driving disease outcomes

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

Modeling Within-host SARS-CoV-2 Infection Dynamics and Potential Treatments

A quantitative model used to compare within-host SARS-CoV-2, MERS-CoV, and SARS-CoV dynamics provides insights into the pathogenesis and treatment of SARS-CoV-2

Immune-Viral Dynamics Modeling for SARS-CoV-2 Drug Development

Modeling SARS-CoV-2 viral kinetics and association with mortality in hospitalized patients from the French COVID cohort

Potency and timing of antiviral therapy as determinants of duration of SARS-CoV-2 shedding and intensity of inflammatory response

Detection of viral infections by innate immunity

DAMPs, PAMPs, and LAMPs in immunity and sterile inflammation

Clinical implications of cellular stress responses

Processes of sterile inflammation

Release mechanisms of major DAMPs

Immunothrombosis in acute respiratory distress syndrome: cross talks between inflammation and coagulation

Hypoxia and inflammation

Neutrophil extracellular traps contribute to immunothrombosis in COVID-19 acute respiratory distress syndrome

Neutrophil extracellular traps in COVID-19

Neutrophil extracellular traps in pulmonary diseases: too much of a good thing

TLRmediated activation of NK cells and their role in bacterial/viral immune responses in mammals

Natural killer cell cytotoxicity: how do they pull the trigger?

NKG2D ligands are expressed on stressed human airway epithelial cells

Type 1 interferons as a potential treatment against COVID-19

In vivo kinetics of SARS-CoV-2 infection and its relationship with a person's infectiousness

Modeling withinhost dynamics of influenza virus infection including immune responses

New features of the software MatCont for bifurcation analysis of dynamical systems

Numerical bifurcation analysis of a class of nonlinear renewal equations

Transit and lifespan in neutrophil production: implications for drug intervention

Timing of antiviral treatment initiation is critical to reduce SARS-CoV-2 viral load

SARS-CoV-2, SARS-CoV, and MERS-CoV viral load dynamics, duration of viral shedding, and infectiousness: a systematic review and meta-analysis. The Lancet Microbe

Improvements to and limitations of Latin hypercube sampling

Kinetics of SARS-CoV-2 infection in the human upper and lower respiratory tracts and their relationship with infectiousness

Morphogenesis and cytopathic effect of SARS-CoV-2 infection in human airway epithelial cells

Mathematical model of a three-stage innate immune response to a pneumococcal lung infection

Within-host models of high and low pathogenic influenza virus infections: The role of macrophages

NK cells in the human lungs

Classification and regression trees. Routledge

Random forests

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

Kinetics of viral load and antibody response in relation to COVID-19 severity

Impaired type I interferon activity and inflammatory responses in severe COVID-19 patients

Functional exhaustion of antiviral lymphocytes in COVID-19 patients

Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19

Remdesivir for the treatment of Covid-19-preliminary report

Efficacy evaluation of early, low-dose, shortterm corticosteroids in adults hospitalized with non-severe COVID-19 pneumonia: a retrospective cohort study

Effect of hydroxychloroquine in hospitalized patients with COVID-19: preliminary results from a multi-centre, randomized, controlled trial

Effect of systemic glucocorticoids on mortality or mechanical ventilation in patients with COVID-19

Corticosteroid treatment of patients with coronavirus disease 2019 (COVID-19)

A review of potential suggested drugs for coronavirus disease (COVID-19) treatment

Merck and Ridgeback's Investigational Oral Antiviral Molnupiravir Reduced the Risk of Hospitalization or Death by Approximately 50 Percent Compared to Placebo for Patients with Mild or Moderate COVID-19 in Positive Interim Analysis of Phase 3 Study