A Prototype QSP Model of the Immune Response to SARS-CoV-2 for Community Development.

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) pandemic requires the rapid development of efficacious treatments for patients with life-threatening coronavirus disease 2019 (COVID-19). Quantitative systems pharmacology (QSP) models are mathematical representations of pathophysiology for simulating and predicting the effects of existing or putative therapies. The application of model-based approaches, including QSP, have accelerated the development of some novel therapeutics. Nevertheless, the development of disease-scale mechanistic models can be a slow process, often taking years to be validated and considered mature. Furthermore, emerging data may make any QSP model quickly obsolete. We present a prototype QSP model to facilitate further development by the scientific community. The model accounts for the interactions between viral dynamics, the major host immune response mediators and tissue damage and regeneration. The immune response is determined by viral activation of innate and adaptive immune processes that regulate viral clearance and cell damage. The prototype model captures two physiologically relevant outcomes following infection: a "healthy" immune response that appropriately defends against the virus, and an uncontrolled alveolar inflammatory response that is characteristic of acute respiratory distress syndrome. We aim to significantly shorten the typical QSP model development and validation timeline by encouraging community use, testing, and refinement of this prototype model. It is our expectation that the model will be further advanced in an open science approach (i.e., by multiple contributions toward a validated quantitative platform in an open forum), with the ultimate goal of informing and accelerating the development of safe and effective treatment options for patients.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

The complex interactions between the in‐host viral dynamics, immune response and tissue damage underlying disease heterogeneity in coronavirus disease 2019 (COVID‐19) are not yet fully understood.

WHAT QUESTION DID THIS STUDY ADDRESS?

Development of a robust quantitative systems pharmacology (QSP) model of COVID‐19 will bring novel medicines to patients faster. Here, we propose a prototype model for further development by the scientific community.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

The model qualitatively captures a “healthy” immune response and a pathophysiological uncontrolled inflammatory response, characteristic of severe COVID‐19 cases. Exploration of model parameters highlights the importance of an appropriate immune response in mediating a balance between tissue damage and viral clearance.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

A prototype model for community development will accelerate the development of a validated QSP model of COVID‐19 and hence support the development of new therapeutics for patients.

Coronavirus disease 2019 (COVID‐19) is caused by the severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2) virus, a novel coronavirus that emerged in 2019. Similar to other respiratory coronaviruses (such as SARS‐CoV and Middle East respiratory syndrome coronavirus), human‐to‐human transmission is primarily via aerosolized respiratory droplets. The functional receptor for both SARS‐CoV and SARS‐CoV‐2 is believed to be the angiotensin converting enzyme 2 receptor. Although the primary infectible cell of SARS‐CoV‐2 is alveolar type 2 cells (AT2), their expression of angiotensin converting enzyme 2 is relatively low, suggestive of multiple factors driving infectibility of AT2 cells and the potential role of nonrespiratory infections contributing to the observed pathogenesis.

Detailed descriptions of the pathogenesis of COVID‐19 have appeared in recent reviews. , Briefly, the host antiviral response can be divided into two categories. First, induction of cellular programs that act to suppress the virus is regulated by type I and III interferons and their downstream genes. Second, a leukocyte response, including phagocytosis of the apoptotic debris by antigen‐presenting dendritic cells and macrophages, leads to T‐cell induction via antigen presentation. Infected cells can then be killed by antigen presenting cell (APC)‐recruited CD8+ cytotoxic T‐cells (CTLs).

It is suggested that the heterogeneity in response to infection, ranging from asymptomatic to death, could be due to a varying and sometimes inappropriate inflammatory response driven by limited induction of the type I and type III interferons and “exuberant” inflammatory cytokine production, including IL‐6 and TNF‐⍺. As a result, alveolar epithelial cells undergo inflammatory cell death, , which, in turn, compromises lung microvasculature and epithelial barrier function leading to pulmonary edema and ultimately hypoxia.

One of the common and potentially fatal comorbidities of COVID‐19 is acute respiratory distress syndrome (ARDS), which can develop secondary to the viral‐mediated injury. ARDS is characterized by an acute lung abnormality alveolar injury produces diffuse alveolar damage, resulting in impaired gas exchange, decreased lung compliance, and increased pulmonary arterial pressure. The etiology of ARDS may include sepsis (most commonly), pneumonia, and over 60 other causes. ,

Elevations in several coagulation biomarkers also have been associated with disease severity and prognosis in patients with COVID‐19. , , Whereas the tight coordination between hemostatic and inflammatory responses is a part of innate immunity, dysregulation in this feedback system leads to thrombo‐inflammation, as observed in cases of sepsis‐associated disseminated intravascular coagulation (DIC). , Although there are some differences, as in DIC, coagulopathy in COVID‐19 is likely triggered due to injury to the endothelial cells by endotoxins or inflammatory cytokines. Endothelial stimulation triggers the activation of the coagulation pathway, which is sustained and amplified by thrombin generation. Platelet activation and endothelial cell stimulation also further lead to production of cytokines and growth factors, including IL‐6.

To quantitatively understand the complex interactions among in‐host viral dynamics, tissue damage, and the immune response, a mathematical model is required. Quantitative systems pharmacology (QSP) models are mechanistic models of pathophysiology and treatment and can support decision making throughout the drug discovery and development process. Hence, a QSP model of the immune response of COVID‐19 has clear applications for understanding the pathophysiology and supporting clinical trial design and analysis for novel treatments.

Specifically, the initial goal of the model presented here is to study factors determining the manner in which the immune response influences severity and duration of SARS‐CoV‐2 infection and potentiates the development of a severe pathophysiological hyperinflammatory state (e.g., in ARDS) observed in certain infected individuals. Accordingly, our mathematical model accounts for the interactions between processes determining viral infection and shedding, the activation of the major host immune response mediators and tissue damage and regeneration (Figure ). In this prototype version of the model, we do not consider the response of the coagulation pathway.

Figure 1

Model diagram of viral dynamics coupled with immune response. Model diagram of severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2) dynamics coupled with immune response. The model schematic shows lung/alveolar space where virus particles infect healthy alveolar cells (AT2), which results in a positive feedback where infected cells produce additional virus particles. The infected cells produce cytokines to activate the innate immunity, which in turn leads to tissue damage of healthy alveolar cells, inducing additional innate immunity response. Although the immune response aids the clearance of infected cells and virus particles, the neutrophil activation can also damage AT2 and AT1 (not shown). The positive feedback, due to infected and damaged cells producing proinflammatory cytokines that enhance immune activation, is closely interlinked with the adaptive immunity. In addition, the anti‐inflammatory signaling via regulatory T (Treg) cells serves to balance the pro‐inflammatory response. A number of physiological and clinical biomarkers including surfactant protein D (SP‐D), ferritin, C‐reactive protein, neutrophils, M1 macrophages, cytotoxic lymphocytes, TH17, TH1, Treg and dendritic cells, IFN, TNF and interleukins IL‐1, IL‐2, IL‐6, IL‐10, IL‐12, and IL‐17 are found in the systemic circulation. ROS, reactive oxygen species.

The disease pathophysiology of viral replication and immune response is complex, making the development of therapeutic intervention challenging. One way to overcome this is the utilization of QSP models to leverage and incorporate existing mechanistic knowledge and data to make forward predictions. The typical development time for models of this scope is at least a year and probably significantly longer to fully report. To shorten this timeline, we have decided to break our traditional cycle for development and publication of QSP models. We are publishing a prototype version of this model that passes a series of unit tests, outlined below, that are designed to demonstrate that the model structure is sufficient for qualitatively capturing a range of plausible responses (Figure ). It is our hope and expectation that researchers across academia and industry find this model a useful basis for further development and application to the discovery, development, and utilization of safe and effective treatments for COVID‐19. We strongly encourage expedient publication of subsequent work based on this. We note that other mechanistic models, with different applications due to their smaller or distinct focus, emerged as we reached this stage of development. , These alternative approaches present a future opportunity for collaboration and sharing of data and findings.

Figure 2

Past and future development plan. We developed a set of ordinary differential equations (ODEs), specified preliminary parameters sets and devised qualitative unit tests to assess whether the model’s structure captures the required behavior (green). The results from these efforts are described in this paper. We are building a database for quantitative testing and seeking external feedback for critical additions to the model (orange). Next, we will incorporate new emerging data, build virtual populations, and ensure that updated model passes the quantitative tests before applying this quantitative systems pharmacology (QSP) model to questions of interest (red).

METHODS

Our approach was enabled by significant adaptation of a prior model of inflammatory bowel disease. The model we are proposing is outlined below and in accompanying Supplementary Materials. Briefly, the model consists of 49 state variables and 206 parameters. An exhaustive list of the states and parameters of the model is provided in the Supplementary Material s. Code is available in full at https://github.com/openP​fizer/​covid​19-immun​e-model.

Viral infection dynamics

We developed a set of ordinary differential equations to describe the dynamics of SARS‐CoV‐2 viral load, based on previously published models of influenza infection. , , In our model, uninfected but susceptible AT2 cells are infected by SARS‐CoV‐2 to form productively infected cells, which shed viable virus particles that might be measured in a sputum sample or nasal wash. In contrast to prior models of viral dynamics, which assume a static initial pool of susceptible target (AT2) cells that cannot be replenished, we account for the regeneration of AT2 cells, intended to be a phenomenological representation of the activation of wound healing mechanisms upon alveolar cell depletion.

Immune response dynamics

Equations from Rogers et al. were adapted to model immune cells that are specifically activated by the presence of virus particles, infected cells, and by alveolar tissue damage due to excessive inflammation. Although the immune response to respiratory viral infections involves multiple compartments, for simplicity we only consider alveolar and plasma compartments.

Both free virus particles and infected cells activate the innate and subsequently adaptive components of host immune response. , Our model specifically describes the virus and infected cell‐induced activation/maturation of alveolar‐resident macrophages, neutrophils, and dendritic cells (DCs). Moreover, we consider the DCs to be the primary APCs and, as such, responsible for the activation of the adaptive immune cell populations. The adaptive immune cell populations considered in the model include the CD8+ CTLs and the CD4+ Th1, Th17, and regulatory T cells. Although the CD8+ CTLs are critically involved in the clearance of infected cells, the CD4+ Th1 and Th17 cells are assumed to potentiate and sustain their activation by maintaining a permissive inflammatory milieu through the secretion of proinflammatory cytokines together with DCs, macrophages, and neutrophils. The CD4+ regulatory T cells are assumed to be the primary negative regulators of the immune response through the secretion of anti‐inflammatory cytokines, such as IL‐10 and TGF‐β. As an example, Eq. 1 describes the activation of CD8+ CTLs as considered in the model.

Where,

Type I IFN induction of CTL activation

IL‐12 induction of CTL activation potentiated by IL‐2

IFNγ induction of CTL activation

IL‐10 and TGFβ‐mediated inhibition of CTL activation

The production of CD8+ CTLs is activated by viral epitope‐responsive mature DCs, and further induced by IL‐12, IL‐2, IFN‐γ, type I IFN, and inhibited by IL‐10 and TGF‐β concentrations. Additionally, the ability of IFN‐γ and IL‐6 to negatively regulate each other’s activity is also incorporated. The clearance rate of CTL is determined by a nonspecific death/deactivation rate (β), and intercompartmental transport rate . Given the uncertainty of viral infection‐induced differential recruitment, we assumed equal immune cell transit rates between alveolar and plasma compartments.

Viral clearance and antiviral dynamics

The simulated viral dynamics are influenced by type I IFN, released by virus‐infected cells and mature DCs, which inhibits the generation of infected cells, thus, implicitly accounting for effects of IFN‐stimulated gene expression products that are known to confer viral resistance. Additionally, virus‐activated CTLs contribute to viral clearance through their cytotoxic effects on infected cells, thereby reducing the rate of infectious viral shedding. Finally, we assume that DCs, macrophages, and neutrophils, upon viral recognition, contribute to overall viral clearance through the phagocytosis of viral particles and apoptotic infected cells.

Virus and immune‐induced damage

As mentioned above, in contrast with most other models of the immune response to respiratory viral infections, we explicitly account for damage and inflammatory death of AT2 cells that might occur due to both viral infection as well as the effects of proinflammatory cytokines. Such immune‐induced damage is particularly important in facilitating and sustaining pathophysiological outcomes associated with cytokine storms. ,

Additionally, we include a population of type I alveolar (AT1) cells, which are assumed to be relatively resistant to viral infection but might perpetuate a hyperinflammatory state by undergoing inflammatory cell death (e.g., pyroptosis and necroptosis/necrosis) as a result of the viral induction of the immune response.

Clinical biomarkers

We modeled multiple disease‐relevant clinical biomarkers. , We incorporate three biomarkers commonly monitored in hospitalized cases of COVID‐19: CRP, which is a general marker of inflammation; and ferritin and SP‐D, which are leakage products of alveolar cell damage. , In our model, we assume IL‐6 induces the hepatic production of CRP, whereas ferritin and SP‐D are produced by the death of damaged alveolar cells. All three biomarkers are assumed to be released into systemic circulation from their respective sites of production.

Model development plan and unit tests

The size and scope of the model is such that for efficient development we decided to report its content at a relatively early point in its development (Figure ). At this stage of our workflow, we have refined the model (green boxes, Figure ) such that it passes key qualitative unit tests.

Unit Test 1: Nominal case of healthy response

Unit test 1 (U1) Purpose: Demonstrate engagement of the innate immune system and pathogenesis of COVID‐19 that might typically be seen in patients with moderate symptom severity.

U1 Criteria: Simulated viral peak at approximately 4 days, viral clearance (viral load < 100 copies/mL) within 30 days, and an engagement of innate immunity. Post‐viral clearance, alveolar cells, activated immune cells, and cytokines return to baseline.

U1 Simulation Method: The model was calibrated to fit the criteria of U1, a physiological representation of a “moderate” immune response to COVID‐19 infection, which resolved upon viral clearance. U1 was used as the foundation for the other unit tests.

Unit test 2: No innate immune recognition of virus

Unit test 2 (U2) Purpose: Demonstrate the contribution of innate and adaptive immune mediators on viral clearance by uncoupling the contribution of the immune response from the other clearance mechanisms accounted for in the model. This nonphysiological example is an important control for model development and ensures the functional role of the immune response in this model.

U2 Criteria: Without innate immunity viral load peak is higher and clearance is slower than U1.

U2 Simulation Method: U1 plus silencing of innate immune activation (no activation of dendritic cells, macrophages, or neutrophils and effect of type I IFNs).

Unit test 3: Sterile inflammatory response

Unit test 3 (U3) Purpose: Demonstrate that there is a damaged cell‐induced inflammatory response independent of the virus‐induced immune response.

U3 Criteria: Damaged cells induce a sterile inflammatory response. Post‐clearance of damaged cells, the alveolar cells, activated immune cells, and cytokines return to baseline.

U3 Simulation Method: U1 with an initial condition of damaged AT1 and AT2 cells at the start of the simulation, but no viral inoculum.

Unit test 4: Sustained inflammatory response

Unit test 4 (U4) Purpose: Demonstrate a sustained inflammatory response characterized by prolonged immune activation even after the virus is cleared.

U4 Criteria: Sustained immune response after virus clearance that leads to alveolar cell damage.

U4 Simulation Method: U1 with increased cytokine‐induced damage to healthy alveolar cells at baseline and after virus exposure.

Sensitivity and strength of response to infection

To determine the sensitivity and strength of response to infection, the half‐saturation constants and the maximal innate immune activation rates of virus, infected cells, and damaged cells were explored. The explored sensitivities and strength of response to infection were 2 and 1 order of magnitude above and below that of the nominal case (U1).

Virtual population analysis

One thousand plausible patients were generated using a Latin hypercube sampling approach. The immune response parameters bounds were one magnitude above and below that of the nominal case (U1) to explore a wide range of parameter space. The viral dynamic parameter bounds were 50% above and below that of the nominal case (U1) to avoid nonphysiological viral dynamics.

RESULTS

Viral and cell damage dynamics

The critical role of viral and immune‐induced tissue damage in facilitating and sustaining pathophysiological response under normal (healthy) immune response, no innate immune recognition of virus, sterile and sustained inflammatory response conditions (as specified by U1–U4 tests) is shown in Figure . As required by U1, the viral load peaks around 4 days, with a peak of around 10 log10 viral RNA molecules/mL and is significantly cleared within ~ 30 days. The simulated time to peak in viral load is in qualitative agreement with clinical observations, which suggest that SARS‐CoV‐2 achieves peak viral loads relatively quickly, roughly coinciding with symptom onset, and with an estimated incubation period of about 4–7 days. Simulated viral loads in the nominal unit test are qualitatively comparable to those typically observed in nasopharyngeal swab samples in patients with COVID‐19, as reported in a recent meta‐analysis of over 2,200 samples in 315 patients by Gastine et al. Moreover, peak viral loads are in agreement with data that have been shown to range between 8 log10 and 11 log10 RNA molecules/mL , , in severe hospitalized cases.

Figure 3

Viral and cell damage dynamics in unit tests. (a) Viral load with an initial inoculum equivalent to 1 million viral RNA molecules/mL in a nasal or sputum sample. (b) Infected cells formed from healthy type II alveolar (AT2) cells. (c) Angiotensin converting enzyme 2 (ACE2)‐negative healthy type I alveolar (AT1) cells assumed to be uninfectable. (d) Susceptible AT2 cells in response to viral infection. (e) Damaged AT1 cells scaled to the maximum of nominal case of 4.45 million damaged AT1 cells and (f) Damaged AT2 cells scaled to the maximum of nominal case of 5.16 million damaged AT2 cells. Median and maximal viral load measurements of nasopharyngeal swab samples from patients with coronavirus disease 2019 (COVID‐19) from Gastine et al. are depicted with green dashed lines. Data from Gastine et al. depicted in this figure was published by the authors (at https://github.com/ucl-pharmacometrics/SARS-CoV-2-viral-dynamic-meta-analysis).

In the case where there is no innate immune response, the model simulation results in a higher viral load and slower clearance of the virus, as required by U2, which is entirely mediated by the complete depletion of infectible AT2 cells in the terminal phase of the dynamics. In the sustained inflammatory case (U4), the viral dynamics are similar to those of the nominal case; however, in contrast to U1, there is a large increase in the damaged AT2 cell population that is sustained post‐viral clearance.

Immune cell and cytokine response

The degree of immune response across tests U1–U4 is defined by both plasma and lung cytokine dynamics (Figure ) and cell populations (Figure ). Although there is currently much uncertainty regarding the tissue level response, we anticipate being able to refine the model once these data emerge (either from clinical or representative preclinical experiments). The nominal case simulation (U1) gives rise to an immune response that resolves as the virus is cleared. Specifically, sufficient virus exposure induces the activation of innate immune cells, such as macrophages, DCs, and neutrophils (Figure a,b,e,f), which exhibit a peak in activity that aligns with the peak in viral load. The adaptive immune cells, activated by APC, peak shortly after the peak in innate immune cells (Figure c,d,g,h) with CD8+, cytotoxic T‐cells serving to clear infected cells. In addition to their clearance mechanisms of virus and infected cells, innate and adaptive immune cells also facilitate the production of cytokines (Figure a–h) that limit the infection and support the resolution of the immune response. For simplicity, we model virus‐activated and damage‐activated immune cells only, which are minimally present prior to viral exposure. In order to depict the qualitative correspondence of our simulations with clinically reported immune cell dynamics over the course of a viral infection we superpose plasma immune cell counts that are virus induced over corresponding immune cell counts that are typically observed in healthy volunteers (Figure ). , , The proinflammatory cytokines in the plasma and lungs follow a similar time course, with the increased cytokines levels in the lungs, reflected in systemic circulation (Figure a–h). We have superimposed median plasma cytokine concentrations measured in patients with severe COVID‐19 from studies by Lucas et al. and Vanderbeke et al. to show that simulated cytokine profiles in the nominal case (U1) are in qualitative agreement with a range of reported observations (Figure a–d). Furthermore, cytokine levels in the alveolar compartment are substantially higher than those in the plasma compartment, in agreement with reports comparing sputum, bronchoalveolar lavage, and plasma samples in patients with severe respiratory infections. , , Finally, the acute activation of the immune response leads to a sharp increase in CRP levels, whereas damage of alveolar cells by infection and immune mediators results in an increase in ferritin and SP‐D. Similar to the immune response dynamics, biomarker levels return to baseline upon resolution of the immune response.

Figure 4

Cytokine response in unit tests. (a) Increase in damage inducing IL‐6 levels in plasma in response to viral infection, commonly used as a marker of inflammation. (b) Type I IFN levels in plasma indicative of endogenous antiviral response. (c) IFNγ levels in plasma, indicative of CD8+ cytotoxic T‐cell (CTL) activation. (d) Immune suppressive IL‐10 levels in plasma indicative of regulatory T (Treg) activation. Corresponding levels of (e) IL‐6, (f) type I IFN, (g) IFNγ, and (h) IL‐10 at alveolar site of infection. Median cytokine levels from patients classified as having severe/critical COVID‐19 from Lucas et al. and Vanderbeke et al. are depicted with the dashed green lines. Raw data depicted in this figure from Lucas et al. was published by the authors, whereas data from Vanderbeke et al. Figure S2 was digitized from their publication using DigitzeIt (https://www.digit​izeit.de/). U1, unit test 1; U2, unit test 2; U3, unit test 3; U4, unit test 4.

Figure 5

Immune cell response in unit tests. (b) Increased levels of total dendritic cells (DCs) in response to detection of viral particles. (b) Time course of total neutrophil (N) levels in plasma indicative of innate immune activation. (c) Total CD4+ T cells in plasma, composed of Th1, Th17, and regulatory T (Treg) subsets. (d) Total CD8+ cytotoxic T‐cells (CTLs) in plasma. Nominal plasma levels of each cell type are superimposed in b–d to depict change in respective total cell count. Corresponding levels of (e) DC, (f) N, (g) CD4 + T cells, and (h) CTL at alveolar site of infection.

Conversely, if there is no innate immune recognition of the virus the simulations show an absence of activation of immune cells and release of cytokines (as required by U2). The response to sterile damage (U3) shows a mild induction of immunity; however, as expected, the response is a function of the extent of damage (not shown). The hyperinflammatory response shows a prolonged engagement of the immune response that extends beyond virus clearance (passing the criteria for U4). Specifically, the increase in proinflammatory cytokines, such as IL‐6 (Figure a,e), results in the increased accumulation of damaged alveolar cells (Figure e,f), which perpetuates a prolonged immune response (in contrast to the resolving immune response observed in U1). This prolonged immune engagement and increased immune‐induced alveolar apoptosis leads to substantial elevation of all simulated biomarkers, which is in semi‐quantitative agreement with clinical observations in patients with COVID‐19 with hyperinflammatory ARDS. ,

Sensitivity and strength of response to infection

To explore the behavior of the model further, we performed parameter sweeps across the sensitivity and strength of the innate immune response to the virus (Figure ). It has been observed that severe pathogenesis of COVID‐19 can be related to an escape of innate immunity and associated poor IFN response. , Our model qualitatively matches this observation, with increased viral‐load and increased infection time as the strength and sensitivity of the response is reduced (Figure ). Interestingly, and with potential relevance for identification of at‐risk patients for putative treatment, we note a bimodal response in IL‐6 and hence CRP (Figure ). Our finding suggests that a high level of inflammation can be induced both by a virus‐mediated innate immune cell activation that is too weak (bottom left corner, Figure ) or too strong (top right corner, Figure ). In the case of weak innate activation, the resulting low‐level of immune‐mediated viral clearance leads to high levels of intrinsic IL‐6 production by virus‐infected cells, caused by the elevated viral load and prolonged virus‐induced stimulation. Importantly, weak innate immune activation and delayed viral recognition by the immune system is hypothesized to lead to poor outcomes of COVID‐19. , Interestingly, the model predicts that the hyperinflammatory response to strong innate immune activation is triggered by an overly sensitive immune response, due to extensive alveolar cell damage despite enhanced viral clearance. Future work will investigate these dynamics and will explore how model and parameter uncertainty/variability, sufficient to explain observed clinical variability, affect these findings.

Figure 6

Model outcomes are sensitive to the strength of the innate immune response. Strength and sensitivity are defined as the maximum rate and the inverse of half‐saturation constant of immune activation relative to the nominal response (outlined in black), respectively. (a) Maximum viral load. (b) Time postinfection for viral road to reduce to 50% of peak. (c) Maximum circulating IL‐6. (d) Maximum circulating CRP. t1/2, terminal half‐life.

Pathophysiological heterogeneity

To evaluate the ability of our model to support the construction of a virtual population, we examined the response to resampling a subset of parameters that describe the relative strength of virus infectivity and the immune system (see Supplementary Material s).

The model is capable of representing cases of robust viral clearance, with viral loads decreasing to 100 RNA molecules/mL within 20 days postinfection, as well as cases where the viral infection is sustained indefinitely (Figure ). Similar to the viral dynamics, we also identified cases of AT2 recovery and the clearance of damaged cells, as well as sustained viral infection and cell damage (Figure ). Preliminary virtual population analysis show the Type I IFN and proinflammatory IL‐6 response correspond to the viral and cell damage dynamics (shown in Figure ). To compare the temporal dynamics of the simulated species depicted in Figure against clinically reported dynamics, we have superimposed published data on longitudinally obtained nasopharyngeal viral load measurements, plasma IL‐6, and IFN‐⍺ measurements from a cohort of patients moderate and severe COVID‐19. We note that whereas we depict the simulated trajectories with respect to time since viral infection, clinical data are almost always characterized with respect to the time since reported symptom onset. Despite substantial uncertainty in the relation between reported time since clinical symptom onset and time since viral infection (i.e., the incubation period), our preliminary virtual population simulations are in qualitative agreement with range of clinically reported observations. In future iterations, we will explicitly calibrate the model response to longitudinal studies characterizing both viral load and immune mediators, which will enable us to better align simulated trajectories with the observed temporal dynamics of these species. Moreover, we will more comprehensively tune the variability in the parameters used to generate virtual populations in order to better differentiate between mild and severe cases of COVID‐19 based on available data. Furthermore, our results are consistent with recent studies, which suggest that type I IFN activity in response to viral infection may be correlated to disease progression. , , The variability of type I IFN and IL‐6 response indicates that the current framework is capable of capturing an antiviral response ranging from mild to critical.

Figure 7

Preliminary virtual population captured variability of viral dynamics to evaluate mild and severe cases of coronavirus disease 2019 (COVID‐19). Shaded grey regions indicate 20th percentile—80th percentile interval, with nominal response (unit test 1 (U1)) superimposed (blue). Simulated time‐course of viral load with an initial inoculum ranging from 100 thousand to 10 million viral RNA molecules/mL in a nasal or sputum sample. (b) Time course of susceptible type II alveolar (AT2) cells in response to viral infection. (c) Time course of infected cells formed from AT2 cells. (d) Time course of infection‐resistant AT1 cells. (e) Time course of damaged AT1 cells due to immune activation and (f) damaged AT2 cells due to viral infection and immune activation. Red dots depict the 20th–80th percentile interval for viral load, type I IFN, and IL‐6 levels in patients classified as having moderate or severe COVID‐19 from Lucas et al. Raw data depicted in this figure from Lucas et al. was published by the authors as online Supplementary Material.

For completeness, the response range of all model species to variability in model input parameters is included in the Supplementary Material s.

DISCUSSION

We present a prototype QSP model of the pathogenesis of COVID‐19. QSP models often serve as evolving repositories of our current quantitative biological and pharmacological knowledge, and have been applied to improve the quality, cost, and speed of clinical development. , However, given the urgency in finding treatment options amidst continuously emerging data coupled with the typically slow development time for a QSP model, we deliberately modified our standard workflows and timelines in order to publish this model as a prototype for feedback and continued development by the scientific community in an open forum.

The presented prototype model is capable of simulating moderate (resolving) and severe (sustained inflammatory) responses to SARS‐CoV‐2 infection. We associate this sustained inflammatory response and damage with an ARDS‐like phenotype. By varying the strength and sensitivity of the response to infection, the model demonstrated that variability in response to infection could contribute significantly to the heterogeneous response in terms of viral and immune dynamics. Furthermore, by varying a small subset of model parameters, we confirmed the model will support development of a virtual population for close matching of clinical data. For this purpose, we include key clinical biomarkers, such as CRP and IL‐6, to link clinical observations to disease state at the tissue level, particular in response to treatment.

The model has been designed so that putative therapies, including anti‐inflammatory or antiviral therapies, are relatively straight forward additions to the model; however, we have refrained from inclusion of pharmacological interventions at this early stage to emphasize that this model requires further development and evaluation/validation prior to these applications. For instance, a preliminary incorporation of antiviral therapies based on an understanding of their mechanisms of action could be performed by accounting for their influence on several parameters influencing viral dynamics in the model (e.g., the rate of productive infection of AT2 cells by virus, the rate of viral shedding, or the endogenous viral clearance rate might be appropriately modified to capture the effects of antiviral therapies on viral entry, viral translation, and viral protein cleavage, respectively). However, if required, the model can be further expanded to incorporate additional mechanistical detail in order to more appropriately model antiviral or anti‐inflammatory therapies based on their specific mechanisms of action.

There are several areas where emerging data will lead to model refinements. In particular, quantitative tissue dynamics and concentrations are yet to be elucidated in literature. For instance, despite several clinical reports of lymphopenia in severe cases of COVID‐19, we do not explicitly account for the emergence of this phenomenon in our model, primarily due to a high degree of uncertainty around the causal mechanisms. Future work could include model evaluation of hypothetical causes of lymphopenia, including differential immune cell recruitment to the site of infection, and infection‐induced or excessive inflammation‐induced immune cell apoptosis. Moreover, although we account for adaptive immune CD4+ and CD8+ T cells, we have not explicitly included a virus‐specific antibody response, given that there is significant uncertainty around the time to onset and duration of such a response. We anticipate that accounting for an antibody response will enable further refinement of viral clearance to better match clinically observed dynamics. , However, we will note that the most severely ill patients at their time of hospitalization appear to have significant viral load—indicating they have not, as yet, mounted a successful antibody response. In this model, we have assumed well‐mixed compartments suitable for modeling with ordinary differential equations. However, spatial effects could lead to different conclusions about certain aspects of COVID‐19 pathophysiology. As such, a future direction for this work is to adapt this model into other frameworks that describe the pathophysiology of COVID‐19 at this level. Incorporating spatial effects could, for instance, provide additional insight into whether the clinical presentation of ARDS in COVID‐19 is atypical compared with other virus‐caused ARDS (as suggested by emerging and controversial data, there appears to be two types of ARDS, with decreased and nondecreased lung compliance).

The next stage of model development is quantitative unit testing to constrain the model using, ideally, clinical data. The existing data reported in the literature presents several challenges. First, it is biased toward hospitalized patients, which limits a comparison with milder pathophysiology. Second, for hospitalized cases, it is not always clear how to compare across groups because the time of infection onset is generally unknown (only the day of admittance). Hence, a future area of work is the development of virtual patients with variable admittance time relative to infection time.

Key model outputs capturing the viral and immune dynamics are clearly sensitive to parameters that control recognition of the virus and the strength of the response to it (Figure ). Furthermore, varying a small subset of key parameters (such as viral shedding, clearance, and endocytosis rates; the activation rate of innate immune cells by virus and infected cells, the rate of cell damage from cytokines, and the clearance rate of infected cells by CD8+ T cells) was able to drive a significant amount of the observed clinical heterogeneity (Figure ). Further data, and model development, will be required to validate the mechanistic sources of the heterogeneity in COVID‐19 pathogenesis.

There are a myriad of potential adaptations and extensions to this model. In particular, other distinct components of the pathogenesis could be included. In particular, the complex feedback between the coagulation and inflammation responses is beyond the scope of the modeling presented here, but quantitative studies of this coordination form an interesting future direction for better understanding of the pathophysiology of COVID‐19 and DIC.

Funding

This research is supported by Pfizer Inc.

Conflict of Interest

Pfizer Inc. supported the research by W.D., R.R., A.S., N.T., and C.J.M., and R.A., W.D., R.R., A.S., N.T., C.J.M., and R.A. were employees of Pfizer during the completion of this study.

Author Contributions

W.D., R.R., A.S., N.T., C.J.M., and R.A. wrote the manuscript. W.D., R.R., and R.A. designed the research. W.D. and R.R. performed the research. W.D., R.R., A.S., N.T., C.J.M., and R.A. analyzed the results.

Figure S1

Click here for additional data file.

Figure S2

Click here for additional data file.

Figure S3

Click here for additional data file.

Figure S4

Click here for additional data file.

Figure S5

Click here for additional data file.

Supplementary Material S1

Click here for additional data file.

Supplementary Material S2

Click here for additional data file.

Supplementary Material S3

Click here for additional data file.