In silico dynamics of COVID-19 phenotypes for optimizing clinical management

Undefirstanding the underlying mechanisms of COVID-19 progression and the impact of various pharmaceutical interventions is crucial for the clinical management of the disease. We developed a comprehensive mathematical framework based on the known mechanisms of the SARS-CoV-2 virus infection, incorporating the renin-angiotensin system and ACE2, which the virus exploits for cellular entry, key elements of the innate and adaptive immune responses, the role of inflammatory cytokines and the coagulation cascade for thrombus formation. The model predicts the evolution of viral load, immune cells, cytokines, thrombosis, and oxygen saturation based on patient baseline condition and the presence of co-morbidities. Model predictions were validated with clinical data from healthy people and COVID-19 patients, and the results were used to gain insight into identified risk factors of disease progression including older age, co-morbidities such as obesity, diabetes, and hypertension, and dysregulated immune response1,2. We then simulated treatment with various drug classes to identify optimal therapeutic protocols. We found that the outcome of any treatment depends on the sustained response rate of activated CD8+ T cells and sufficient control of the innate immune response. Furthermore, the best treatment –or combination of treatments – depends on the pre-infection health status of the patient. Our mathematical framework provides important insight into SARS-CoV-2 pathogenesis and could be used as the basis for personalized, optimal management of COVID-19.

Other potential treatments include anti-in ammatory drugs and anti-thrombotic agents. Because of the systemic in ammation seen in many patients, anti-in ammatory drugs have been tested, including anti-IL6/IL6R therapy (e.g., tocilizumab, siltuximab) and anti-JAK1/2 drugs (e.g. barcitinib). It is not clear if these drugs will be effective as stand-alone treatments, particularly after the recent failure of tocilizumab in a phase III trial 3,[24][25][26] . In addition, given that a common complication of COVID-19 is the development of coagulopathies with microvascular thrombi potentially leading to the dysfunction of multiple organ systems 4,5 , anti-thrombotic drugs (e.g., low molecular weight heparin) are being tested. Recognizing the interactions of COVID-19 with the immune system 27 , the corticosteroid dexamethasone have been tested, showing some promising results. Given the large range of patient comorbidities, disease severities, and variety of complications such as thrombosis, it is likely that patients will have heterogeneous responses to any given therapy, and such heterogeneity will continue to be a challenge for clinical trials of unselected COVID-19 patients 28 .
Here, we developed a systems biology-based mathematical model to address this urgent need. Our model incorporates the known mechanisms of SARS CoV-2 pathogenesis and the potential mechanisms of action of various therapeutic interventions that have been tested in COVID-19 patients. In previous work, we have exploited angiotensin receptor blockers (ARBs) and angiotensin converting enzyme inhibitors (ACEis) for the improvement of cancer therapies and developed mathematical models of the renin-angiotensin system in the context of cancer desmoplasia [29][30][31][32] . Using a similar approach, we developed a detailed model that includes lung infection by the SARS-CoV-2 virus and a pharmacokinetic/pharmacodynamic (PK/PD) model of infection and thrombosis to simulate events that take place throughout the body during COVID-19 progression (Figure 1 and Supplementary   Figure 1). The model is rst validated against clinical data of healthy people and COVID-19 patients and then used to simulate disease progression in patients with speci c co-morbidities. Subsequently, we present model predictions for various therapies currently employed for treatment of COVID-19 alone or in combination, and we identify protocols for optimal clinical management for each of the clinically-observed COVID-19 phenotypes.

Model Description
The model includes SARS-CoV-2 infection, the renin angiotensin system (RAS), in ammatory and anti-in ammatory cytokines, innate and adaptive immune cells, and factors involved in the coagulation cascade ( Fig. 1). SARS-CoV-2 enters the cell by docking to ACE2, a key component of the RAS. ACE2 can be membrane-bound or soluble, and it regulates in ammation by converting Ang II to Ang 1-7 and Ang I to Ang 1-9; as opposed to Ang I and Ang II, which lead to in ammation. Ang 1-7 and Ang 1-9 have anti-in ammatory effects. Intracellular virus initiates in ammatory pathways through toll-like receptors and NFkB, which produces interferons and other in ammatory cytokines. The viral antigens, along with in ammatory cytokines, cause activation of naïve T cells, creating virus-speci c T effector cells. T cell activation is controlled by viral antigen strength and the presence of PD-L1 / PD-1 inhibition. We combine in ammatory cytokines into a single variable, but explicitly account for IL6 production via the trans pathway in epithelial and endothelial cells and the canonical pathway in immune cells. In the presence of in ammatory cytokines and virus, neutrophils can produce neutrophil extracellular traps (NETs).
Because the virus can infect endothelial cells, we also consider viral dissemination via the blood stream, and the possibility of systemic infection and thrombosis. We include the major organs in a PK/PD model, with physiological blood ow patterns explicitly modeled. Infection of endothelial cells, combined with high levels of in ammatory cytokines in the plasma, can result in thrombosis. Damage to virally-infected endothelial cells and the production of NETs can exacerbate the thrombosis, and microthrombi can enter the blood stream to accumulate in other organs, including the brain, heart and lung. We use a simpli ed model of the coagulation pathways, assuming that formation of microthrombi is proportional to the number of infected endothelial cells, the presence of neutrophil NETs, and the level of in ammatory cytokines. Transport of oxygen from the alveolar space to the blood vessels in the lung is calculated using a modi ed diffusion model, which accounts for damage-induce thickening of the alveolar membrane 33 .

Results
Model validation with data of healthy people and various COVID-19 phenotypes.
Model predictions were validated with clinical data from healthy humans and severe COVID-19 patients [34][35][36] . Figure  2A Subsequently, we focused on the clinically observed phenotypes of COVID-19. We used the model to predict the time evolution of critical disease variables, namely the virus load, levels of IL6 and other pro-in ammatory cytokines, the formation of microthrombi in the lung, the numbers of neutrophils, macrophages and activated cytotoxic CD8 + T cells, and the blood oxygen saturation (SpO2). The COVID-19 patient phenotypes considered include young patients (age<35), older patients (age>65) who more frequently require hospitalization, female patients, as well as patients with co-morbidities such as hypertension, obesity and diabetes, and patients with a dysregulated immune response, whose condition is characterized by high levels of pro-in ammatory cytokines ( Figure 3). Supplementary Table 2 presents the model parameters that were modi ed from the baseline values to account for these various patient phenotypes. In all phenotypes, the viral load increases during early lung infection with divergent trajectories seen after day 5, depending on the levels of activated T cells. In younger patients (age<35) with a healthy immune system, the sustained recruitment of T cells results in a reduction in viral load and in ammation, as well as a decrease in neutrophils and macrophages. All these effects cause a signi cant reduction in the formation of thrombi and restoration of oxygenation.
In general, the simulations resulting in poor outcome were due to increased baseline in ammation or more active innate immune response combined with less effective adaptive immunity. This was the case -in varying degreesfor older patients and those with diabetes, obesity, hypertension, and dysregulated immunity ("hyper-in amed"; see supplement Table 2 for the parameters varied in each population). In the latter case, the adaptive immune response is intact and clears the virus, but the innate immune response is sustained, perpetuating in ammation, thrombosis and hypoxemia.
An outstanding question is why males tend to have more severe COVID-19 disease compared with females 38 .
Proposed mechanisms include the higher CD4/CD8 T cell ratio in females 39 and androgen-induced differences in susceptibility to viral entry into the cells due to higher TMPRSS levels in males 40 . To simulate COVID-19 in an older female, we decreased the production of naïve T cells and the rate of virus entry into the cells 41 . The simulations predict lower viral load compared with the older males, with correspondingly less in ammation and hypoxemia.
These results suggest that reduced viral entry into cells due to lower TMPRSS2 level can explain the improved outcomes in females-even if the adaptive immune response is not as vigorous 38 .

E cacy of current treatments on COVID-19 progression
To identify optimal treatment protocols, we simulated the effects of the various treatments currently being investigated for COVID-19. We rst simulated a patient of age >65 years (assumed to have baseline in ammation and impaired T effector cell function) and performed simulations of various treatments given at early (day 3, SpO2 92%) and late (day 7, SpO2<90%) stages of the disease. For the purpose of the simulations, day 0 is the time of arrival of the virus in the lung, so signi cant symptoms are expected starting around day 3. Simulating treatment on day 3 represents early treatment (e.g., upon hospital admission); later treatment (day 7) could be due to delayed diagnosis/hospital admission or in response to disease progression. Treatments that were considered include: heparin (anti-coagulant), anti-viral therapy, dexamethasone, ARBs, ACEi, hrACE2, anti-IL6 and anti-IL6R treatment.
Unsurprisingly, we found that heparin decreases thrombus formation and improves oxygenation when administered during early stages of disease, but has no effect on viral load. The anti-viral remdesivir is effective in reducing viral load in early stages, but less effective in reducing thrombus formation and improving oxygenation. Interestingly, dexamethasone is helpful when started later, but early administration can prevent the production of activated T cells, making it di cult to reduce the viral load. This mirrors the data from the RECOVERY trial, where dexamethasone was demonstrated to improve outcomes in patients requiring ventilatory support, but not in those with milder symptoms. These simulations also identify IL-6 directed therapies as potential therapeutics that may bene t patients if started early in the disease course. As far as inhibitors of the RAS are concerned, both ARBs and ACEi show modest bene ts only when administered early, whereas the model predicts a signi cant bene t from the use of hrACE2 on day 7, but not day 3.
Optimizing clinical management of COVID-19 by combining therapies Because of the diverse mechanisms involved in COVID-19 pathogenesis, clinical management currently involves combination of multiple therapies in an effort to optimize therapeutic outcomes. To investigate the effects of combination treatments, we performed simulations combining various therapeutic approaches in pairs, using the best treatment time for each single treatment (i.e., day 3 or day 7, Table 1). We also accounted for older patients and those with co-morbidities and dysregulated immune responses. Table 2 presents a summary of the combined treatments for older patients (>65 years); the results for the other COVID-19 phenotypes are presented in Supplementary Tables 3-6. According to the model predictions, effective clinical management of older patients (with some baseline in ammation and impaired adaptive immunity) involves combination of heparin (day 3) with dexamethasone (day 7), which can improve oxygenation and decrease micro-thrombosis signi cantly; however, these treatments have little effect on viral load. For patients with obesity or hypertension, the combination of heparin with dexamethasone is again bene cial compared to other treatments. In addition, the use of anti-IL6/anti-IL6R therapy could also be considered in both these patient populations. For diabetic patients, the combination of heparin with dexamethasone can improve oxygenation and thrombus formation similar to the other scenarios. Other treatments predicted to be bene cial for this phenotype is the combination of dexamethasone with hrACE2 or anti-IL6/anti-IL6R. These combination therapies in diabetic patients can treat hypoxemia and coagulation but also reduce the viral load. For patients with dysregulated immune systems, early blockade of IL6 combined with dexamethasone on day 7 give optimal results; anti-IL6/anti-IL6R combined with anti-viral therapy and/or heparin could also be considered.

Discussion
COVID-19 is a complex disease that can affect multiple organs and is characterized by extremely heterogenous presentation. As with other critical care syndromes such as septic shock, and non-COVID-19 ARDS disease severity results from a complex interplay of viral replication, adaptive and maladaptive immune response and patient comorbidities. Such heterogeneity is a challenge to the development of effective therapies as candidate approaches may elicit different responses depending on the phase of illness or patient comorbidities. Because of the complexity of the underlying mechanisms and interactions among the various components involved in the disease, the response of COVID-19 patients to any treatment is not intuitive and for that reason we developed a highly sophisticated mathematical model to provide insights into the underlying mechanisms and predict optimal treatment strategies.
A main conclusion of our study is that disease progression and outcome of any treatment largely depends on the response rate of activated CD8 + T cells and subsequent control of the innate immune system response. A sustained activation of CD8 + T cells along with the control of the populations of neutrophils and macrophages can lead to a decrease in viral load and in ammation. This, in turn, will improve arterial oxygen saturation levels due to the decreased endothelial damage and micro thrombosis, and limited formation of neutrophil extracellular traps (NETs) 42 . In line with these conclusions about the underlying mechanisms of the disease, our model predicts that anti-viral and anti-in ammatory drugs that were rst employed to treat COVID-19 might have limited e cacy, depending on the stage of the disease progression. Furthermore, our simulations suggest that an optimal approach would be to enhance the adaptive immune response in the early stages while limiting harmful in ammation in the later stages of the disease. We also found that addition of heparin to the treatment regimen of COVID-19 can improve therapeutic outcomes. In case of patients with co-morbidities (obesity, diabetes, hypertension) or dysregulated immunity, the treatment regimen could further include anti-in ammatory drugs (e.g., anti-IL6/anti-IL6R) and RAS inhibitors and hrACE2 for the diabetic patients in particular.
Our mathematical modeling framework is subject to certain limitations. The model is able to account for the currently-known components and interactions involved in disease pathogenesis but results in a large number of model parameters. Most of the values of these parameters were taken from the literature. However, because COVID-19 is a new disease, many of the parameter values are not known yet. For these parameters, it was necessary to make reasonable assumptions about their values, as summarized in Supplementary Table 1. We further validated the model predictions quantitatively with available clinical data for healthy and infected humans and qualitatively based on clinical observations for the different COVID-19 phenotypes. Our mathematical framework can be further re ned as new mechanisms and more data become available and can be adapted for other diseases similar to COVID-19.
Another limitation of the model is that it does not account for adverse effects of the various treatments considered in this study, which could affect the therapeutic outcome. For instance, heparin use might be associated with hemorrhage and RAS inhibitors can reduce blood pressure. Even though incorporation of adverse effects could be important, they would not affect the basic conclusions of the study and were omitted to avoid adding further complexity to the model. In conclusion, this study presents a mathematical representation of the known mechanisms of COVID-19 and could be utilized by the scienti c community as a useful tool for further understanding the disease and investigating the bene ts of treatments on a patient-speci c basis. To study how the virus affects systemic events, such as in ammation and thrombosis, we coupled the lung model with a PK/PD model of viral infection. The PK/PD model represents the major organs of the body as compartments connected in an anatomical manner by the blood and lymphatic circulations (Supplementary Figure 1). The organ compartments are then further subdivided into vascular and extravascular sub-compartments, and each organ has a draining lymph node compartment. Brie y, the main transport processes for the biochemical species include: i) convective and diffusive transport across capillary walls, ii) reversible, non-saturable, nonspeci c binding in the extravascular compartments, iii) reversible, and saturable, speci c binding of virus to endothelial cells, and probabilistic infection of bound cells. Any infected cell generates microthrombi, which enter the circulation and can accumulate in target organs including brain, heart and lung. Each species can be produced, bound or degraded within each compartment. Blood ow is distributed according to known ows through each organ, and the overall mass balance allows calculation of the concentration of each species in each compartment over time. The PK/PD model includes the following key processes involved in the tra cking of viruses: i) transport from the lung via the systemic circulation, ii) entry into endothelial cells via binding to ACE2, iii) entry into the cell and replication, iv) exit from the cell and entry into the blood circulation. The important outputs of the model are the level of thrombosis in each organ and the viral load. When coupled with the lung model, we are able to analyze the dynamics of these readouts in light of local lung pathologies and predicted cytokine levels. Data availability. All data supporting the ndings of this study are available in the paper and the Supplementary Information.

Methods
Code availability. The COMSOL code will become available in Zenodo upon acceptance of the manuscript.
Declarations Tables   Table 1. Summary of model predictions for the therapeutic outcome of currently employed treatments initiated either in the beginning (day 3) of the disease or later (day 7). The table presents the results at the end of the simulation (day 20). Values have been normalized to the corresponding initial values except for SpO2.    Dynamics of disease progression. Time evolution of key model variables: viral load, IL6 concentration, proin ammatory cytokines concentration, activated T cells, microthrombosis in the lung, saturated oxygen levels, neutrophils and macrophages for young (<35 years old) and old (>65 years old) patients, females and patients with co-morbidities: hypertension, obesity and diabetes, as well as for hyper-in amed patients and for anti-androgen treatment. Values have been normalized to the corresponding initial values except for S(Po2).