Divergent COVID-19 Disease Trajectories Predicted by a DAMP-Centered Immune Network Model

COVID-19 presentations range from mild to moderate through severe disease but also manifest with persistent illness or viral recrudescence. We hypothesized that the spectrum of COVID-19 disease manifestations was a consequence of SARS-CoV-2-mediated delay in the pathogen-associated molecular pattern (PAMP) response, including dampened type I interferon signaling, thereby shifting the balance of the immune response to be dominated by damage-associated molecular pattern (DAMP) signaling. To test the hypothesis, we constructed a parsimonious mechanistic mathematical model. After calibration of the model for initial viral load and then by varying a few key parameters, we show that the core model generates four distinct viral load, immune response and associated disease trajectories termed “patient archetypes”, whose temporal dynamics are reflected in clinical data from hospitalized COVID-19 patients. The model also accounts for responses to corticosteroid therapy and predicts that vaccine-induced neutralizing antibodies and cellular memory will be protective, including from severe COVID-19 disease. This generalizable modeling framework could be used to analyze protective and pathogenic immune responses to diverse viral infections.

of the main text provides a schematic of the interactions and processes considered for our mathematical model, which represents an abstraction and simplification of the myriad of immune events occurring during a viral infection. It considers the interaction of four compartments: 1) Virus-infected cells, CV, 2) Tissue/organ/organism dysfunction that drives inflammation and serves as an indicator of host health, D, 3) Innate inflammatory components, I, and 4) Adaptive Immune components including antigen presenting cells, CD8 T cells, and virusspecific antibodies, A. The mean field dynamics are described by differential equations (1)-(4). Each compartment represents a variety of mediators sharing a common role in the response.
The virus-infected cell population, CV, is scaled to have units of 10 6 cells, such that a value of = 10 −6 corresponds to one infected cell. The variables I and A have arbitrary units given that they each represent various cell types, cytokines, chemokines, and other mediators considered to be a part of the respective response. The units of D are also arbitrary as it is not linked to one specific biomarker. Parameters are described in Table S1, some of which are based on our previous work (37,91).

Virus-infected Cells, CV
Viral-cell infection dynamics are often modeled in terms of compartments for the free/extracellular viral particles, target or susceptible cells, viral infected cells, and sometimes 'unproductive' viral infected cells (33)(34)(35)(36) where the CV, P and S are fractions and the decay and interaction rates can depend on other compartments. We reduced these dynamics to a single compartment by assuming 1) viral particle clearance from the vicinity of cells is fast compared to their transmission rate between cells and thus the viral dynamics can be approximated to be in quasi-steady state, P = (kV/kP)CV; and 2) that the intrinsic death rate of susceptible cells is slow compared to the time course of the infection, giving S + CV = 1 (i.e. S + CV is conserved) and thus we can set S = 1 -CV . This results in a single compartment reduced dynamics for the virus-infected cells CV given as: We use this as the basis for our CV compartment for the virus-infected cell population that increases as a result of the process of the virus hijacking the cellular machinery of its target cell to produce virions intracellularly that are released into the extracellular environment to infect new target cells. This construct along with our scaling choice implies that time zero in our model corresponds to the time when there are C × 10 6 infected cells already present. We model this increase of the virus-infected cell population, CV, with the density dependent logistic growth term: Logistic growth is bounded above, which for this population can be best attributed to a limited number of target cells for the virus to infect (33). This type of growth starts off exponential and then slows as the population nears a carrying capacity, 1 .
The adaptive response, A, has dual action on the virus-infected cell population. A inhibits the growth of CV via virusspecific antibodies by preventing additional target cells from becoming infected. This effect is modeled by dividing the growth term of CV by 1 + ( ) 2 : (1 − ) 1 + ( ) 2 However, this mechanism of inhibition is only used when considering vaccination simulations (including simulations of existing SARS-CoV-2-specific immunity). If ≪ 1 then inhibition is weak until A has grown considerably; and when A is zero, as it would be at the beginning of a naive infection, there is no inhibition of virusinfected cell growth. In all the core scenarios presented, we assume that A directly kills CV, representing the cytotoxic effects of CD8 T cells. This is modeled as a simple mass action interaction between A and CV occurring at The two interactions for the induction of I by D and the self-recruitment/self-upregulation of I, were modeled with terms that followed the convention used in previous modeling work by (37) and (73)  . The function R given by equation (5) consists of terms for the induction of I by D and the self-recruitment/self-upregulation of I: + . Inhibition of these interactions by the adaptive response, A, is also incorporated by dividing by 1 + ( ) 2 (i.e. multiplying by the function 1 (1, ) given after equations (1)-(4) above) . Using the above, the following term in equation (3) represents these events: In addition, we utilized the same functional form used in (37)

Adaptive Immune Response, A
During a viral infection, the initial innate immune response will induce an adaptive immune response for a more specific and effective response against the virus-infected cells. As mentioned earlier, the variable A representing the adaptive response serves to directly kill virus-infected cells (i.e. cytotoxic capability of CD8 T cells); and inhibit the infection of other cells (via antibodies). The adaptive response is also responsible for signaling that inhibits the inflammation of the innate immune response as the more specific response is developed. However, if the innate inflammatory response becomes/is quite strong, it can serve to inhibit the initiation/recruitment of the adaptive response. The adaptive equation also encompasses the role of antigen presenting cells as the component of the response that provides the initiation of the adaptive response by the innate.
The priming of the adaptive response by innate inflammation assumes that a naive A population, 0 , is primed by I initially but, as I becomes larger, it inhibits this induction of A: The 0 population is modeled as 0 = 0 − 0 0 − 0 , which includes a source term, a natural death term, and a priming term by I. This equation is then set to quasi-steady state so that the formula for 0 becomes 0 = 0 0 + and is inserted in place of 0 , giving the term: The adaptive response includes CD8 T cells that after being primed will undergo clonal expansion. We include a term for this in equation (4) of A which assumes that expansion is dependent on the CV population to indicate immune signaling from virus-infected cells and sensing by the T cells. Also, we assume that expansion of A will have a self-saturating effect to account for the self-regulation of the response. The functional form for saturation uses a hill function with hill coefficient 2 to emphasize a steeper threshold-like response to self-regulation: ⋅ ⋅ 2 2 + 2 . As with the other variables, A, has a natural decay term, − . Also, a source term, , is included to explore an existing source of virus-specific memory cells if considering prior exposure to re-infection or a level of immunity built up over time with a vaccine. Thus, A is modeled as:

Other Modeling Notes and Considerations
We are not considering a direct negative effect from innate immune factors on the virus-infected cell population.
Instead the presence of virus-infected cells will induce signals via D that will initiate innate immune mediators, I, which in turn initiate specific adaptive immune mediators that can directly kill the virus-infected cells.
= Table S1: Model parameter descriptions, values, and notes. As specified above, the values given are those used to simulate the seven core scenario using the ODE system (1)-(4).
With a parameter set specified, including initial conditions for the variables, the equations are solved with established numerical algorithms for estimating solutions of ordinary differential equations. These algorithms are standard in computational mathematics software such as MATLAB and XPPAUT, both of which were utilized in the analysis and simulation of the model. Parameter values are admittedly difficult to estimate for a model such as this, given the lumped nature of the variables and parameters and the lack of data available for SARS-CoV-2. Many of the model parameters govern the rates of various processes represented in the model (e.g., growth rates, decay rates, activation/priming rates) while others govern the strength of inhibition of one variable on another (e.g., xvainv, xiainv). Still, other parameters determine the levels of variables at which their effects on processes are at halfsaturation (e.g., xdv, xdv, xia, xag). In any of these roles, the value of any one parameter is essentially tied to an individual's genetic variation, co-morbidities or even lifestyle history.
A baseline set of parameter choices were based on previous work and known inflammatory dynamics (37) as well as estimates. The differences between archetypes were due to changes in only four parameters (kdi; xdi; xiainv; and kai), while the others were fixed. However, the advantage of modeling is the ability to explore the parameter space to determine the myriad of model behaviors both transiently and in the long-term. Thus, many parameter sets can be explored around the baseline set of parameters. Time courses of the model variables can be generated for each parameter set and together form an interpretable picture of the pathogenesis of the immuno-inflammatory response to the infection, showing both the transient evolution of the variables as well as an overall outcome.
Multiple parameter sets can lead to similar manifestations of a response that fall under one of the four major archetypes presented in this manuscript. We chose seven parameter sets that provided a representation of the possibilities and which all arose from changes made to one parameter set such that the overall differences between the number of parameter values differing from one scenario to the next is on the order of three. For instance, the pathways in the model that lead from the Mild recovery scenario to the Moderate recovery scenario to the Severe recovery scenario are all with respect to a change in two parameters (namely the parameters responsible for inducing D via I: kdi and xdi); though, we acknowledge this is not the only pathway to achieve this progression. Below we describe these parameter changes among the different scenarios relative to the baseline parameter set used for the Mild Scenario, where these particular modifications do not represent unique ways to achieve these various scenarios/archetypes. Recurrence Scenario -Recrudescence Archetype (Figure 6). The parameter xiainv controls the strength of A to inhibit the inflammatory response, I. Thus, an increase in this parameter will increase the inhibitory effect on I by A.

Recovery
In the Recurrence scenario (of the Recrudescence archetype), this parameter was increased from the value of xiainv=0.1 in the Mild Scenario to a value of xiainv=0.5 and, in addition, the parameter xdi=1.0 in Mild was decreased to 0.5 (See above discussion on this parameter). completely, leading to a situation in which virus-infected cells are not eliminated yet remain at lower, perhaps undetectable levels that continue to drive a lengthy immune response having signs of recovery only to be followed by reprise. This gave rise to the scenario name has it reflected descriptions by some health professionals of patients who struggle with the disease progression/symptoms for months at a time.

Long-lasting Disease Scenario -
Uncontrolled Infection Archetype (Figure 8). As in the Long-lasting Disease scenario, the Uncontrolled Infection scenario decreased the parameter kai from a value of 3.0 in the Mild Scenario to kai =0.1. However, all other parameters were the same as that for the Mild case. This gave rise to a more severe instance of an inadequate immune response (relative to that in the Long-lasting Disease scenario) which causes uncontrolled growth of the virus-infected cell population as well as sustained inflammation.

Parameters involved in treatment Simulations. Simulations of corticosteroid therapy modified parameters kii and
kid to have values of zero during the therapy window as described in the text. These two parameters are rates of induction of inflammation, I, respectively by the innate inflammatory response itself through self-enhancement and by damage-induced mechanisms. Thus, setting these two parameters to a value of zero during a window of time essentially shuts off any production of I completely for that window. After the treatment window ends, these two parameters are restored to their pre-treatment values.
Vaccination Protocol Implementations: short-term vaccine effects and existing immunity.
We implemented two different vaccination-type scenarios in the core scenarios. One, termed 'short-term vaccination,' looks at the short-term effects of a vaccine that is adjuvant-based and where a COVID infection happens withing 2 weeks of the last dose. The other, termed "existing immunity," looks at considering already established, existing level of virus-specific immunity at the time of infection.
For the short-term vaccination, a function ( ) = − = −3 was defined to represent an exponentially decaying stimulus of the I response and was incorporated within the function R in the I equation to emulate induction of I by an adjuvant which decays over time. The induction of I leads to the priming of A, the specific adaptive response to the virus. If CV is present, the A response will be magnified in a CV-dependent manner as described in model development. The decay period of the adjuvant stimulus function was chosen with initial condition at Day 0 normalized to a value of 1.0 and exhibiting exponential decay at a rate such that the stimulus is negligible by Day 3.
The parameter, xvainv, governs the strength of A to inhibit the growth/expansion of the CV population. This parameter is set to 0 when not implementing a vaccination protocol, thus only allowing for direct (e.g., cell-mediated) killing of CV in the model without vaccination. However, for vaccination we assume xvainv >0 (specifically, xvainv =1.0) to allow the vaccine to be responsible for the induction of antibodies needed to prevent virus from infecting more cells, and therefore inhibiting the expansion of the virus-infected cell population, CV, in addition to the direct killing of CV by A. There is no build-up mechanism/equation of long-term memory yet modeled and thus the vaccination effects are assumed to be viable for the short-term, which drove the decision to initiate an infection 14 days after the second vaccine dose was simulated initiated (with 8x105 virus-infected cells). The two vaccine doses administered 28 days apart was chosen to emulate current vaccine protocols (94).
To simulate an established build-up of immunity, we utilized the parameter Ae in the A equation to represent an existing source of A and set this parameter to a nonzero amount (Ae =0.05) and used xvainv =1.0 both chosen so that the effects would not be unrealistically too weak or too strong (data not shown) in inhibiting CV. Thus, at time zero when an infection is initiated (with 8x105 virus-infected cells), A levels would not be zero (as in the other scenario) but instead (0) = and therefore able to respond immediately to the viral threat.  Table S2: Exploration of the impact of increasing initial viral load on the time to peak infection. The days to peak CV level as a function the initial CV were calculated for each of the seven core modeling scenarios and for eight different initial doses of CV ranging from the equivalent of one virus-infected cell to one million virus-infected cells. Time to virus positivity was calculated from incidental COVID-19 patients as 11.3±5.3 (n=8; range 4-22) days as discussed in the main text. From these data, we calibrated to a CV peak time of 11.1±4.5 days (range 8.5-19.1 days) in subsequent simulations based on an initial CV dose equivalent to 8x10 5 virus-infected cells (bolded row).