Mechanistic Models of COVID-19: Insights into Disease Progression, Vaccines, and Therapeutics.

The COVID-19 pandemic has severely impacted health systems and economies worldwide. Significant global efforts are therefore ongoing to improve vaccine efficacies, optimize vaccine deployment, and develop new antiviral therapies to combat the pandemic. Mechanistic viral dynamics and quantitative systems pharmacology models of SARS-CoV-2 infection, vaccines, immunomodulatory agents, and antiviral therapeutics have played a key role in advancing our understanding of SARS-CoV-2 pathogenesis and transmission, the interplay between innate and adaptive immunity to influence the outcomes of infection, effectiveness of treatments, mechanisms and performance of COVID-19 vaccines, and the impact of emerging SARS-CoV-2 variants. Here, we review some of the critical insights provided by these models and discuss the challenges ahead.

Introduction

COVID-19, caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) virus, has infected more than half a billion individuals worldwide, resulting in the death of over 6.2 million people as of 26th April 2022 (WHO). SARS-CoV-2 is a small, single-stranded RNA virus that primarily infects the respiratory system. The clinical outcome of infection is heterogeneous, ranging from cure without symptoms to severe disease culminating in death [1]. Nearly two dozen COVID-19 vaccines have been approved for clinical use, with several hundreds of candidates in clinical and preclinical development [2]. The approved vaccines are powerful in conferring protection, particularly against severe disease [3]. Simultaneously, new antivirals such as Paxlovid and nirmatrelvir have been shown to significantly reduce COVID-19 deaths in clinical trials [4,5]. As SARS-CoV-2 continues to evolve, enormous efforts are therefore underway to understand how SARS-CoV-2 variants may impact vaccine and antiviral treatment effectiveness [6,7].

Over the past two years, SARS-CoV-2 viral dynamics models that quantitatively analyze viral load evolution in patients (Figure 1 ), and quantitative systems pharmacology (QSP) models that integrate systems biology and human physiology with clinical pharmacology (Figure 2 ) and can thus predict the system's response to therapeutics, have been developed. These models are playing a key role in understanding SARS-CoV-2 infection as well as in identifying optimal vaccine dosing regimen and antiviral treatment strategies. In this review, we first discuss insights gained into COVID-19 disease progression and the influence of host immune responses on suppressing infection. Next, we focus on studies that have modelled COVID-19 therapeutics, followed by descriptions of vaccine QSP frameworks for vaccine dose optimization and predicting efficacy. Finally, we conclude by outlining modelling efforts that would enhance our ability to combat current and future pandemics.

Figure 1

Statistical and mechanistic models of SARS-CoV-2 viral dynamics. (a) Schematic of typical predictions from a statistical model (line) calibrated to patient data (circles) capturing the proliferation phase, τ, and the clearance phase, τ, of SARS-CoV-2 infection. (b) Schematic of a mechanistic model of within-host SARS-CoV-2 viral dynamics, describing the interactions between target cells, infected cells, virions, immune response, and immune-mediated cells refractory to infection, and the mode of action of different drugs in clinical development and application. Therapeutics include antivirals such as entry and protease inhibitors that interfere with the virus life cycle, monoclonal antibodies that neutralize the virus and prevent infection, and immune modulators such as steroids and Janus kinase inhibitors that may prevent severe disease. (c) Representative fits with a mechanistic viral dynamics model to longitudinal viral load data from two patients are shown; adapted from Ref. [26].

Figure 2

Mechanistic quantitative systems pharmacology (QSP) models of SARS-CoV-2 vaccines. (a) Defining the scope of QSP modelling frameworks; adapted with permission from Ref. [59]. (b) Schematic of the prototypical Certara vaccine simulator, a QSP model incorporating antigen physiologically-based pharmacokinetics (PBPK) upon either mRNA or adenovirus vaccine administration, immune dynamics and resulting antibody as well as cellular responses, and human population attributes including age, sex, race, ethnicity, HLA genetics, and physiological characteristics of special populations such as immunocompromised, pregnant, and lactating women. (More details in [43]). (c) Predictions of the optimal prime-boost interval based on IgG responses upon Moderna vaccine administration; adapted with permission from Ref. [43].

Virus dynamics models of COVID-19 disease progression and immune response

Like many acute viral infections, multiple clinical [8], [9], [10], [11], [12] and non-clinical studies [13] (some with very high temporal resolution) have shown that COVID-19 is mainly characterized by two phases (Figure 1a): a proliferation phase, where viral loads at sites of infection increase exponentially due to abundant target cells until it attains a peak, and a clearance phase, where a lack of target cells and/or mounted antiviral host immune responses lower and clear the virus from the system [9]. The relative timescales of the two phases, τ and τ (Figure 1a), are determined by viral infectiousness and strength of immune responses, and are signatures of transmissivity and severity of infection, therefore influencing the choice of treatment. For example, an individual in the proliferation phase may need longer isolation and closer medical monitoring due to a higher risk of complications than a person in the clearance phase. Therefore, mathematical modelling of longitudinal virus load data has been employed to not only quantitatively determine in vivo disease dynamics at the individual and population levels, but also the interplay between virus-host immune responses and the resulting pathology.

Statistical models are a class of models employed for quantifying SARS-CoV-2 virus dynamics, where virus load time series measured using RT-qPCR are fit to piecewise linear regression or other mathematical models to estimate parameters underlying individual and population virus trajectories (Figure 1a) [9]. While such phenomenological models were useful for deployment in the initial stages of the pandemic with limited available data, a glaring limitation is a lack of mechanistic insight into disease progression, co-evolving immune dynamics, heterogeneity of responses, and estimates of efficacies of various therapeutics and vaccines. More mechanistic virus dynamics models [8,14] (Figure 1b, c) largely address the above limitations and have provided valuable insights into the progression of many acute and chronic infections caused by HIV [8,15], hepatitis C virus (HCV) [16], influenza [17,18], and Zika virus [19] among others. Example model fits to SARS-CoV-2 virus trajectories from two patients are shown in Figure 1c.

The basic model of acute viral infection, also known as the ‘target cell limited model’, predicts COVID-19 disease progression to be a result of interactions between SARS-CoV-2 virions, infected cells, and the availability of uninfected target cells. The effects of innate and adaptive immune responses are lumped into the viral dynamic parameters such as the rate constants of infection and death of infected cells. In these models, after the proliferation phase, viral load declines until complete clearance due to the lack of target cells for new infection. Such models have quantitatively described and compared the within-host dynamics of MERS, SARS-CoV-1, and SARS-CoV-2, and predicted that SARS-CoV-2 had a shorter time from the symptom onset to the acute infection viral load peak compared to the other two coronaviruses [20]. Subsequent models [21] have accounted for the known delays between cellular infection and viral production by including an eclipse phase for infected cells, allowing for more accurate estimates of SARS-CoV-2 viral dynamic parameters and the within-host reproductive ratio that determines the ability of the virus to establish infection in an individual. Modelling efforts to quantify viral dynamic parameters pertaining to the SARS-CoV-2 variants-of-concern [22,23] and to link viral kinetics to an individual's disease transmission probability are also underway [24], [25], [26].

Unless individuals have pre-existing cross-reactive T-cell immunity to SARS-CoV-2 due to prior exposure to other infections, it takes about a week or two post-infection to mount effective T cell responses against SARS-CoV-2. Many models have added an explicit T cell compartment to understand how T cell dynamics affects SARS-CoV-2 infection and disease progression. Here, antigen presentation on infected cells activate cognate CD8+ T cells, which upon differentiation into effector cytotoxic phenotypes kill the infected cells. Inclusion of this mechanism in virus dynamics models has been helpful in describing the rapid viral load declines observed in some individuals during the clearance phase. Recently, a minimal model consisting of the essential interactions between infected cells, CD8+ T cells, and innate immune response was developed to understand the diverse outcomes of SARS-CoV-2 infection, from clearance without symptoms to severe illness followed by death [27]. By analyzing data from patients with different degrees of infection severities, the model hypothesized that variations in the timing and strengths of the innate and the CD8 T-cell responses may be a factor that potentially explains the observed spectrum of SARS-CoV-2 infection outcomes. The role of the humoral response, the other component of adaptive immunity, in resolving SARS-CoV-2 infection in non-vaccinated individuals is still unclear. By modelling and analyzing longitudinal data of viral loads and SARS-CoV-2-specific IgG and IgM antibodies in patients, a recent model suggested that the surge in IgG levels occurs 5-10 days post symptom onset and had a minor impact on viral clearance [28].

Several classes of immunomodulatory agents, such as interferons, steroids such as dexamethasone, kinase inhibitors such as baricitinib, and interleukin inhibitors such as tocilizumab, are either approved or in clinical development. The choice, timing, combination, and duration of these treatments are likely to depend on the stage and severity of the infection. For instance, in hamster models challenged with SARS-CoV-2, early prophylactic interferon treatment conferred protection whereas later treatment did not [29]. Future modelling efforts, especially larger QSP models [30], [31], [32] incorporating viral dynamics and associated innate and adaptive immune responses at various anatomical sites while accounting for host factors such as comorbidities and age, may accurately predict how different arms of the immune system would interact to determine diverse infection outcomes and guide rational optimization of treatments with immune modulators.

An obstacle with developing larger disease models incorporating the immune response is that, often, detailed kinetics are only available for viral load measurements while the spatio-temporal dynamics of innate and adaptive immune arms, which are mechanistically intertwined with pathogen and disease dynamics, are harder or nearly impossible to obtain, especially from deeper tissues. Therefore, this is a common issue with many disease models and perhaps a grand challenge in the field of QSP and virus dynamics modelling. Three approaches are usually taken to ameliorate this issue. Firstly, modelling insights are often borrowed from similar diseases. For example, influenza virus dynamics models were very useful for developing SARS-CoV-2 virus dynamics models since both are respiratory diseases with qualitatively similar viral kinetics in the respiratory tract [8]. The second approach is to employ measurements of serum immune biomarkers which could potentially reflect the disease state of the individual even if the events occur in solid tissues, and then constrain mechanistic models to reproduce these biomarkers. With COVID-19, the dynamics of many serum immune biomarkers such as antibodies (IgM/IgG/IgA), cytokines (interferons, interleukins, chemokines…), and immune cells (B cells, T cells, monocytes…) have been measured with unprecedented resolution, including in deeper, immunologically-active sites such as draining lymph nodes, thus yielding opportunities for parameterization of large disease models. Careful calibration-qualification cycles with independent training and testing datasets can be performed to prevent model overparameterization, and robust QSP workflows [33], [34], [35] are available that may guide modelling efforts in this regard. Further, non-linear mixed-effects approaches for parameter estimation can be employed to compare and rank multiple complex QSP models, each based on mutually exclusive scientific hypotheses, against multiple streams of pathogen-immune dynamics data. The best ranked QSP model, based on its ability to capture the data as well as parsimony in model structure and number of parameters, can then be chosen (we have previously employed this approach for selecting HIV virus dynamics models [15]). Lastly, experiments with animal models may yield immune response dynamics even in solid tissues, but this must be carefully interpreted given that such models typically do not capture all aspects of human infection.

Models of COVID-19 therapeutics

Mathematical models combining virus dynamics with pharmacokinetic/pharmacodynamic (PK/PD) models of therapeutics have been extensively used to estimate the efficacy and identify the mechanism(s) of action of therapeutics (Figure 1b). Models achieve this by introducing an efficacy factor varying between 0 and 1 that lowers viral production or de novo infection rate constants depending on the mode of drug action [14]. For instance, Pfizer's Paxlovid and nirmatrelvir, which work by inhibiting the activity of the SARS-CoV-2 protease [4,5], would be mechanistically modelled as lowering the viral production rates in the model (Figure 1b). Similarly, Merck's molnupiravir, which may inhibit viral replication by increasing viral mutation rates [4], would render a fraction of progeny virions non-infectious [8,36]. Anti-SARS-CoV-2 neutralizing monoclonal antibody drugs such as REGEN-COV (casirivimab and imdevimab), sotrovimab, bamlanivimab, Evushield (tixagevimab and cilgavimab), and etesevimab, block infection of target cells and simultaneously induce phagocytosis and clearance of virus. We expect that analyzing viral kinetics in patients treated with these antivirals in the future will: (i) quantify the treatment efficacies, (ii) clarify the in vivo dominant mechanism(s) of action of these antivirals, (iii) optimize dosage, as shown recently for the combination of bamlanivimab and etesevimab using population PK/PD and viral dynamics modeling [37,38], and, (iv) predict optimal antiviral combinations with high genetic barriers [8] that may lower the probability of emerging drug-resistant variants, akin to combination therapies for other infectious diseases caused by viruses such as HIV and hepatitis B/C.

Several models have predicted the need for early administration before the viral peak of therapeutics that blocked new infections or lowered viral production, as most cells are likely to be infected by the time viral load peaks, and late administration is likely to inhibit only a small proportion of cellular infection [20,21,32,39]. The current estimates of the critical therapeutic efficacy, which is the minimum drug efficacy above which infection cannot be established, is ∼80-95%. Until effective SARS-CoV-2 therapeutics are developed, one strategy would be to repurpose available drugs for treatment. However, analysis of SARS-CoV-2 virus dynamics in patients treated with different repurposed drugs revealed the efficacies were far below the critical efficacy [21]. Therefore, combination treatments that display synergistic effects are desirable [40]. We recently developed a model of SARS-CoV-2 entry into target cells and predicted that targeting two host proteases, which mediated SARS-CoV-2 entry into target cells via independent pathways, could be synergistic [41]. Such frameworks can be advantageously combined with mathematical models of drug resistance to elucidate optimal combination therapies that would simultaneously maximize synergy and minimize the probability of escape variants.

Application of QSP models for vaccine development

With a raging pandemic, no definitive sterilizing cure (yet), and substantial vulnerable populations worldwide, regulatory-approved COVID-19 vaccines are critical for safely suppressing infections and mortality. In the current scenario with multiple variables such as limited vaccine supplies, availability of multiple vaccine types with varying levels of protection across populations, evolving SARS-CoV-2 variants, and uncertainty regarding the durability of elicited immune responses [42], vaccine biosimulations, including with QSP [43,44], mathematical [12,[45], [60]], and epidemiological [46] models, can deliver potent intelligence to support decision-making for maximal vaccine deployment. For example, early clinical trials with the Pfizer/BioNTech (NCT04368728), Moderna (NCT04470427), and Oxford-AstraZeneca (NCT04400838) vaccines were conducted with either 3-week, 4-week, or 4-6-week intervals, respectively, between the first (prime) and second (boost) doses. However, in December 2020, the UK approved a 12-week prime-boost interval to maximize the first-dose vaccine coverage across the population [47], against recommended intervals based on clinical trials and thus sparking debate among experts that a longer prime-boost interval (which was untested in clinical trials at that point) could potentially lead to decreased immune responses. Around that time, our group at Certara deployed a novel, proprietary, vaccine QSP model [43] (Figure 2b), which was extensively calibrated with pre-clinical and clinical mRNA vaccine data, to predict that a longer 7-8-week prime-boost interval would elicit maximal antibody titers, and reassuringly, a 12-week interval would still yield higher titers than a 3- or 4-week interval (Figure 2c). Strikingly, ∼6 months after this prediction by the QSP model, the PITCH human study [48] was performed by employing an extended dosing interval of 6-14 weeks with the Pfizer/BioNTech vaccine, and this study quantitatively confirmed the model prediction [49] highlighting the power of QSP frameworks to robustly forecast complex what-if scenarios in the face of uncertainty. Mechanistic modelling also predicted that along with higher titers, the neutralization potential of the elicited antibodies may also be improved with delayed boost dosing (and/or lowered prime dose amount) due to altered B cell selection stringency in germinal centers [45]. Next, we discuss other important questions about which QSP predictions can help ameliorate uncertainty and thus significantly contribute to decision-making.

Do the elderly mount less effective immune responses upon vaccination? Extensive simulations with QSP models calibrated on age-stratified clinical trial data suggested that the elderly (65–85 years) elicited reduced antibody titers [43]. The reduction was exacerbated at lower vaccine doses potentially due to lower antigenic stimulation, and with time due to waning immune responses [43]. Thus, multiple boosters with the current vaccine doses may be necessary for the elderly, especially with the rise of SARS-CoV-2 variants that partially escape pre-existing immune responses. On the other hand, increasing the prime and boost vaccine (antigen) dose amounts in the elderly were predicted to elicit comparable antibody titers to that within younger individuals, suggesting that age-stratified dose optimization – higher doses for the old and lower doses for the young – may prove to be beneficial in resource-constrained settings. Further, QSP biosimulations suggest that antibody titers upon vaccination dip below protective levels within a months across both young and old populations [43], thus perhaps necessitating periodic booster shots for all, similar to annual influenza vaccines.

Similar predictions with QSP models can in principle also be made for special populations where detailed clinical trials may be impractical due to recruitment challenges, or unethical, such as patients with impaired hepatic or renal systems, but for whom swift go/no-go vaccine policy decisions must be made. Such populations include neonates, pediatrics populations (children and juveniles), pregnant and lactating women, immunocompromised adults, patients on immune-suppressive therapies, individuals with cancer and/or chronic diseases like HIV and hepatitis B/C, patients who have recently undergone various surgeries, individuals on specific comedications, diabetics, and others. Since QSP biosimulations deal with computer-generated virtual patients whose characteristics and responses are constrained by mechanistic models and real-world data, they can rapidly enable virtual clinical trials (VTs) on these populations. VTs are generally much faster and cost-effective compared to clinical trials, and are often used to streamline, optimize, and sometimes potentially replace these trials. The speed factor in performing VTs and exploring myriad combinations of principal variables is pivotal for COVID-19 in the context of optimizing the vaccine supply chain and ensuring global vaccine equity. An example of such a scenario with urgent consequences for global vaccine coverage is evaluating heterologous vaccine regimens, where different vaccines are administered for the prime and subsequent booster shots. Limited clinical trials for a few vaccine combinations are available, and they show promising but non-intuitive results such as an impact of the order of vaccine administration – the Oxford-AstraZeneca (prime) + Pfizer/BioNTech (booster) combination may be better than Pfizer/BioNTech + Oxford-AstraZeneca, and both mixed vaccine regimen may be better than administering two homologous doses of the Oxford-AstraZeneca vaccine [50]. However, all possible vaccine combinations with varying doses and prime-boost intervals cannot be realistically tested across unvaccinated populations with large-scale clinical trials. Here, mechanistic vaccine QSP models [43] which incorporate quantitatively calibrated antigen physiological-based pharmacokinetic (PBPK) modeling, and thus account for differences across vaccine types (e.g., mRNA vaccines vs. adenovirus vector vaccines) as well as individual patient characteristics (Figure 2b), can enable rigorously performed VTs to evaluate all heterologous vaccine combinations in silico and thus greatly streamline decision-making. Finally, mechanistic mathematical models can also be employed to integrate available experimental data and generate testable predictions on updating the antigen in SARS-CoV-2 vaccines to match that of arising variants, such as Omicron [60].

Predicting the efficacy of COVID-19 vaccines

What is the relationship between immune responses evoked by vaccination and protection against infection? To answer this, immune correlates of protection must be identified for vaccine-induced immune responses [51]. While cellular responses such as primed T cells undoubtedly confer some protection, neutralizing antibody (NAb) titers, usually the first line of defense against SARS-CoV-2 viruses and infection in vaccinated individuals, appears to be a robust immune correlate of vaccine efficacy from statistical analyses of clinical and epidemiological data [51], [52], [53]. Higher the NAb titers elicited by a vaccine, higher the protection accorded to the population immunized by that vaccine [52]. A testable mechanistic link between NAb titers and protection, which could potentially unravel the mechanistic underpinnings of COVID-19 vaccine efficacies and enable robust efficacy predictions, was recently proposed by us [12,54]. First, we hypothesized that inter-individual variability in NAb responses upon vaccination can be described by assuming that the elicited NAb repertoire within an individual was a randomly sampled subset of a shape space constituting all possible NAbs. Next, we constructed such a shape space from analyses of in vitro dose-response neutralization curves of ∼80 anti-SARS-CoV-2 NAbs and then recapitulated the reported NAb repertoires of convalescent patients by sampling from that subset, thus validating the shape space. Finally, we developed a within-host model of SARS-CoV-2 infection which quantitatively captured human viral kinetics in vivo, coupled individual NAb responses after vaccination (sampled from the shape space) to virus dynamics, performed VTs of vaccinated populations with a spectrum of immune and virus characteristics, and predicted the overall vaccine efficacies. Our model predictions quantitatively captured the efficacies of major vaccines from multiple clinical trials and thus offer a framework for assessing and including other immune correlates of protection, thus opening avenues for optimizing vaccine deployment.

Outlook: Combining models of within-host disease progression, therapeutics, and vaccines with population epidemiological models

We anticipate that the optimal development and deployment of therapeutics and vaccines, especially in times of urgent need such as global pandemics, would be greatly accelerated by integrated, semi-mechanistic/statistical, population-based models. Such multiscale frameworks would seamlessly combine models of basic vaccine immunobiology [45,55], individual-scale infection models including within-host models of pathogen-immune dynamics and disease progression [30,31], action of anti-pathogen and immunomodulatory therapeutics, vaccine QSP models accounting for antigen PBPK and inter-individual variations [12,43,44,56], and population epidemiological models [46,57,58] that would enable real-world efficacy predictions in the context of spreading and evolving disease as a function of virus/patient characteristics, social policies, resource constraints, and human behavior. The key challenges for parameterizing such large QSP models are mainly the requirements for detailed mechanistic knowledge, seamlessly integrating multiscale phenomena from cellular to population levels, extensive literature surveys or experiments to obtain parameter estimates, and robust workflows [33], [34], [35] that curate datasets for efficient calibration-validation cycles. The wealth of fine-resolution and diverse data generated by massive consortium efforts during the COVID-19 pandemic is expected to offer a unique opportunity to calibrate and validate such mega-models. Subsequently, these frameworks may enable static or dynamic predictions of waning anamnestic immune responses and protection against future disease while accounting for inter-vaccine differences, history of infection and vaccination in convalescent individuals, heterologous vaccinations, special populations, individual genetics, pharmacogenomics, combination therapies in vaccinated and unvaccinated individuals, and epidemiological aspects such as optimal supply/deployment of vaccines in the face of resource constraints and rapidly transmitting disease, pathogen evolution, impact of interventional measures such as lockdowns and vaccination strategies, and changes in human behavior. Such simulations may enable both ‘nowcasting’ and forecasting in the face of uncertainty, and better guide not only drug/vaccine discovery and dose optimization, but also holistic policies towards the appropriate administration of vaccines and therapeutics for effectively combating future diseases and pandemics.

Funding

No funding.

Competing interests

R.D. was a former employee, and A.M.K., P.H.v.d.G. are current employees of Certara QSP who were involved in the development of Certara's vaccine QSP model. All the authors declare no conflict of interest.