Calibration and Validation of a Mechanistic COVID-19 Model for Translational Quantitative Systems Pharmacology - A Proof-of-Concept Model Development for Remdesivir.

With the ongoing global pandemic of coronavirus disease 2019 (COVID-19), there is an urgent need to accelerate the traditional drug development process. Many studies identified potential COVID-19 therapies based on promising nonclinical data. However, the poor translatability from nonclinical to clinical settings has led to failures of many of these drug candidates in the clinical phase. In this study, we propose a mechanism-based, quantitative framework to translate nonclinical findings to clinical outcome. Adopting a modularized approach, this framework includes an in silico disease model for COVID-19 (virus infection and human immune responses) and a pharmacological component for COVID-19 therapies. The disease model was able to reproduce important longitudinal clinical data for patients with mild and severe COVID-19, including viral titer, key immunological cytokines, antibody responses, and time courses of lymphopenia. Using remdesivir as a proof-of-concept example of model development for the pharmacological component, we developed a pharmacological model that describes the conversion of intravenously administered remdesivir as a prodrug to its active metabolite nucleoside triphosphate through intracellular metabolism and connected it to the COVID-19 disease model. After being calibrated with the placebo arm data, our model was independently and quantitatively able to predict the primary endpoint (time to recovery) of the remdesivir clinical study, Adaptive Covid-19 Clinical Trial (ACTT). Our work demonstrates the possibility of quantitatively predicting clinical outcome based on nonclinical data and mechanistic understanding of the disease and provides a modularized framework to aid in candidate drug selection and clinical trial design for COVID-19 therapeutics.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

Currently, there is an increasing number of coronavirus disease 2019 (COVID‐19) drug candidates demonstrating promising results in the nonclinical stage, but there is no mechanistic COVID‐19 model that can translate these nonclinical findings into clinical outcomes in a quantitative (predicting the values of clinical end points) and independent (no clinical efficacy data available for model adjustment) manner.

WHAT QUESTION DID THIS STUDY ADDRESS?

Is it possible to use a mechanistic model for predicting the clinical outcome of a COVID‐19 therapeutic based on nonclinical (in vitro) pharmacology data and our current knowledge of the underlying biology of the disease?

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

This study shows the possibility of using in silico models for translating nonclinical pharmacology data to clinical outcome by using remdesivir as a proof‐of‐concept example. It also provides a well‐calibrated and validated, modularized model to perform similar nonclinical‐clinical translations for other potential COVID‐19 therapies.

HOW MIGHT THIS CHANGE CLINICAL PHARMACOLOGY AND THERAPEUTICS?

Given the time‐consuming nature of clinical studies and the urgent need for treatments of COVID‐19, the full mechanistic model can play a key role in quickly identifying drug candidates that not only showed promising nonclinical results, but also have the best chances of demonstrating clinical efficacy. It can also aid in the subsequent clinical trial design and dose selection by simulating the drug actions in a quantitative and mechanism‐informed manner.

In the first dozen years of the 21st century, there have been 2 outbreaks of highly pathogenic coronaviruses, severe acute respiratory syndrome coronavirus (SARS‐CoV) in 2002 and middle east respiratory syndrome coronavirus (MERS‐CoV) in 2012, wherein spread of animal viruses to humans resulted in mild to fatal respiratory illness in humans. The year 2019 saw the initial outbreak of a new coronavirus, SARS‐CoV‐2, which subsequently spread rapidly across 7 continents and became a serious threat to all human beings around the world.

Similar to other coronaviruses, SARS‐CoV‐2 initiates the infection process by binding to the host cells via its viral spike protein and then fusing with the cell membrane and releasing viral RNA into the cell. After genome replication and transcription, newly translated structural proteins and replicated RNA genomes are assembled in the endoplasmic reticulum‐Golgi compartment, where new virus particles are formed and exit the host cell through exocytosis to complete the growth cycle. The variability in initial viral load and the strength of host immune responses can lead to differences in inflammatory and immune responses and, consequently, different courses of the disease. Typically, more severe pneumonia in patients with coronavirus disease 2019 (COVID‐19) is linked with higher initial viral load, more virus replication, massive inflammatory cell infiltration, and elevated pro‐inflammatory cytokine responses, leading to a known phenomenon called “cytokine release syndrome.” Despite greater knowledge and awareness of this disease since the outbreak, and widespread prevention and control policies set by most countries on the planet, the COVID‐19 pandemic continues throughout all areas of the world and has accounted for more than 200 million confirmed cases and the deaths of 4.0 million people as of August, 2021. Therefore, there is an urgent need to speed up the traditional drug discovery process to develop effective therapeutic agents against COVID‐19.

One way to accelerate the drug development process is to utilize nonclinical (in vitro and/or in vivo) data for better clinical trial design and dosage selection. This translational science strategy is important for drug repurposing, where the effort is to identify COVID‐19 therapies that have already been proven safe and approved to treat other diseases in humans. However, a major challenge is the translatability from nonclinical to clinical settings due to various factors, such as the lack of immune responses for in vitro assays and potential species differences for in vivo assays. Some drugs with promising nonclinical data have failed to show efficacy against COVID‐19 in clinical studies, such as hydroxychloroquine and lopinavir. Up until now, despite a large number of studies having been done for potential therapeutics against SARS‐CoV‐2, the antiviral drug remdesivir is the only approved treatment in the United States. In addition, the US Food and Drug Administration (FDA) has authorized emergency use of other drugs and biological products to treat or prevent COVID‐19.

Remdesivir can be intracellularly metabolized into its active form, nucleoside triphosphate analog (TP), that inhibits the RNA‐dependent RNA polymerase (RdRp) of SARS‐CoV‐2 in vitro, and has been shown to improve time to recovery for patients with COVID‐19 in the Adaptive COVID‐19 Treatment Trial (ACTT). However, even with remdesivir, the reported in vitro SARS‐CoV‐2 inhibition potency (half‐maximal effective concentrations (EC50s)) is highly inconsistent and it is unclear if these in vitro findings can be quantitatively translated to clinical efficacy. Indeed, to date, we are unaware of studies showing a quantitative prediction of COVID‐19 clinical trial outcomes based on in vitro or other nonclinical data.

In this work, we propose a framework to bridge the gap between nonclinical and clinical end points by combining mechanistic understanding of pharmacokinetics (PKs)/pharmacodynamics and nonclinical data into a human disease model of COVID‐19 with physiologically detailed immune responses. Through a stringent calibration/validation strategy, we show that it is possible to independently predict the clinical outcome (time to recovery) of the remdesivir clinical trial (ACTT) through a model‐informed drug development approach based, primarily, on in vitro data.

METHODS

Detailed methods can be found in the supplementary document and only essential procedures are described here.

COVID‐19 model construction

The modular structure of the COVID‐19 model is based on the human immune response model for influenza infection published previously with a few additional components, such as IL‐6 response, T cell apoptosis within lymph nodes, and antibody response, etc. The full structure of the model can be found in Figure . In order to illustrate the model clearly, we limited Figure to showing only the essential processes and avoided showing the reaction details. Reaction formulations, reaction rates, parameters, and differential equations are listed in full detail in the supplementary document.

Figure 1

Schematic diagram of the mechanistic in silico COVID‐19 model. The in silico model includes two main submodels: (a) COVID‐19 disease model and (b) pharmacological model. The disease model consists of three submodules: virus life cycle, immune response module in lung, and immune response module in the lymphatic compartment. The pharmacological model describes the pharmacokinetics processes that translate drug dose to clinical exposure and then links it to the disease model. In this figure, remdesivir is shown as an example, whose pharmacological model is linked to the disease model at the point of virus replication due to remdesivir's known effect of suppressing SARS‐CoV‐2 replication. Drugs with other mechanisms of action may have other connection points. For example, pharmacological models of inflammation suppressors (e.g., IL‐6 antibody) can be linked to the disease model at the point of IL‐6 production. Multiple pharmacological modules (submodels) could be added to simulate drug combinations. COVID‐19, coronavirus disease 2019; EC50, half‐maximal effective concentration; Epi, epithelial cell; IFNα, interferon alpha; IL‐6, interleukin 6; IgG, immunoglobulin G; IgM, immunoglobulin M; SARS‐CoV‐2, severe acute respiratory syndrome coronavirus 2. This figure is a high‐level summary of essential processes connected by sequential arrows, with solid lines ending in filled circles representing cell death induced by immune responses (for example cytotoxic T cells killing infected epithelial cells). More details about reaction interconnections, rates, parameters, and equations are provided in the Supplementary Document.

Pharmacological model

The pharmacological model includes two sequential parts, remdesivir compartmental PK model and remdesivir intracellular metabolism model. The compartmental PK model includes the central and peripheral compartments and is used to convert remdesivir intravenous dose to plasma concentration. Each compartment is considered to be well‐mixed, with a uniform concentration throughout. Distribution, redistribution, and elimination taking place between the compartments follow first‐order kinetics and are characterized by their rate constants. By using a protein bindingvalue12.1%, total remdesivir plasma concentration is converted to free plasma concentration that is assumed to be equal to the extracellular remdesivir concentration of remdesivir in the lungs. The free fraction of extracellular remdesivir enters the cell through diffusion and can either be degraded or converted to the active ingredient nucleoside TP. It is of note, that the metabolism is a multistep process involving multiple intermediate metabolites. We lumped all these reaction steps into a single step to produce the final metabolite TP, as this is the active metabolite whose concentration was measured in vitro.

Computational methods

All simulations are run in R version 4.0.2 (www.R‐project.org) through the deSolve package (http://desolve.r‐forge.r‐project.org) so that model equations can be implemented in C to accelerate the simulation speed. The FDA High Performance Computing cluster was used for large scale parallel computation for simulating virtual populations. Model calibration was conducted through covariance matrix‐adapted evolution strategy through the R package CMAES (https://CRAN.R‐project.org/package=cmaes). The full simulation code is also available in the GitHub repository https://github.com/FDA/Mechanistic‐COVID‐19‐Model).

RESULTS

Design of the overall COVID‐19 model structure

Expanded based on the human immune response model for influenza virus that we previously developed, our in silico COVID‐19 model (Figure ) includes the disease model that represents the virus infection process as well as human immune responses, and a pharmacological model that represents the PKs of potential therapies. The disease model starts with the virus life cycle that includes three main steps: virus infection, virus replication in different subpopulations of lung epithelial cells, and virus release. The two immune response modules recapitulate the innate and adaptive immune responses in the lungs and lymphatic compartments, respectively. In response to virus infection, the innate immune system in the lung compartment produces interferon alpha (IFNα) in order to convert free (uninfected) epithelial cells to refractory epithelial cells that are resistant to virus infection and also promote the death of virally infected cells. In addition, cytokines are released that may lead to inflammation and damage of both healthy and infected tissues. Whereas many pro‐inflammatory cytokines can be released, IL‐6 was found to be the key cytokine that consistently shows statistically significant correlation with disease severity of COVID‐19 (Table ) and, hence, was used in our disease model to represent cytokine response to the virus infection. The immune response module (in the lymphatic compartment) captures various essential processes such as the activation and migration of antigen presenting cells from the lungs to draining lymph nodes, the conversion of naive T cells to mature helper (CD4+) and cytotoxic (CD8+) T cells, the migration of cytotoxic T cells from the lymph nodes to the lungs to attack infected cells, and the maturation of B cells for antibody (immunoglobulin G (IgG) and immunoglobulin M (IgM)) production and release (Figure ).

Calibration of the disease model

Clinical time course of virus titer in nasopharyngeal swab and plasma profiles of IFNα, lymphocytes, IL‐6 and IgM and IgG antibodies were used simultaneously for parameterizing the disease model (model calibration; Figure ). Due to the different clinical outcome of patients with severe vs. mild symptoms, , and inherent heterogeneity of the disease, we estimated two (sub)populations of parameters, one for patients with severe and another for mild disease.

Figure 2

Disease model calibration. Left column is mild case. Right column is severe case. Points: clinical data; error bar: clinical variability (1 times SD); black line: best fitting curve. Gray band: predicted population behavior (90% confidence interval). Y axes have the same scale for each row and X axes have same scale and are shared in each column.

Viral titer

The mean values of the viral load time course for patients with severe and mild disease were well‐matched by the disease model (Figure ). As the initial viral loads (amount of virus entering the body at day 0) cannot be measured directly in clinical studies, they are usually treated as an unknown quantity and estimated along with other parameters. , Interestingly, our disease model estimated that, on average, patients with severe disease had a higher initial viral load compared with the patients with mild disease (mean value: 99890 (copy num/mL) vs. 8,926 (copy num/mL)). This is consistent with the notion that the initial viral load is one factor that contributes to subsequent disease severity. It is important to note that, on the individual level, some simulated patients had relatively low initial virus load but subsequently developed into a more severe disease status (Figure ), suggesting that other factors (e.g., degree of immune response) also play a role in determining the disease severity.

Interferon alpha

Clinical data reveal that patients with COVID‐19 with low IFNα responses are associated with mechanical ventilation and poorer health outcomes. Consistent with this, our disease model, being calibrated with clinical data, illustrates that the simulated patients in the mild disease population have higher IFNα levels than those in the severe disease population (Figure ). Moreover, it has been reported that SARS‐CoV‐2 can delay the response of the innate immune system, resulting in right‐shifted (delayed) IFNα peak compared with the peak of virus, a crucial mechanism for the virus to evade early immune responses. Our simulations confirmed this pattern (Figure compare panel c to a). In contrast, the influenza infection model we previously developed showed the IFNα peak occurs before the virus titer peak after influenza virus infection, providing a plausible explanation to why SARS‐CoV‐2 can cause much more severe diseases than influenza (Figure ).

Interleukin‐6

In patients with COVID‐19, a hyperinflammatory state known as cytokine storm is associated with the disease severity and IL‐6, a common component of cytokine storm, is a well‐documented predictor of clinical severity. Figure demonstrate that, consistent with the reported IL‐6 clinical time course, both patients with mild and severe simulated disease have an elevated IL‐6 level, with patients with severe disease peak the IL‐6 level reaching about 3‐fold higher levels than patients with mild disease.

Lymphopenia

Based on longitudinal analysis of peripheral lymphocytes after SARS‐CoV‐2 infection, the development of lymphopenia in patients with COVID‐19 was mainly related to the significantly decreased absolute counts of T cells in peripheral blood. Although the exact mechanisms behind SARS‐CoV‐2 induced lymphopenia are elusive, possible mechanisms include (direct or indirect) induction of apoptosis of T cells and enhanced lymphocyte redistribution (migration of lymphocytes from peripheral blood to the lungs or lymphoid organs). Accordingly, our model implements the mechanism of virus‐dependent death of T cells in the lymph node and migration of cytotoxic T cells to the lungs. We used the T cell percentage decrease in the lymph node to represent the overall degree of lymphopenia and were able to replicate clinical data showing a more significant degree of lymphopenia in patients with severe vs. mild disease (Figure ).

IgG and IgM

IgG and IgM clinical time courses were used to represent the antibody response in the human body after SARS‐CoV‐2 infection. IgM and IgG appear earlier, and their titers are significantly higher in patients with severe disease compared with patients with mild disease. Our disease model was able to capture this trend (Figure ).

Calibration of remdesivir PK model

Next, we complemented the COVID‐19 disease model with a pharmacological model to evaluate clinical efficacy of COVID‐19 therapies. As a proof‐of‐concept example, we selected remdesivir, one of the first therapies approved by the FDA for COVID‐19. Due to the complex PKs of remdesivir, the pharmacological model includes two sequential components (Figure ): remdesivir compartmental PK model and remdesivir intracellular metabolism model.

Figure 3

Full pharmacological model. Pharmacokinetic (PK) model includes the central and peripheral compartments used to convert remdesivir (RMD) intravenous infusion dose to total plasma concentration. A scaling factor of 0.12 is used to convert total plasma RMD to free plasma RMD. An intracellular metabolism model converts extracellular free remdesivir to intracellular active metabolite (triphosphate).

A two‐compartment PK model was used to convert remdesivir intravenous (i.v.) dose to plasma concentration. By assuming a drug free plasma concentration of remdesivir in equilibrium with the tissue (extracellular) concentration, the compartmental PK model is linked to the intracellular metabolism model, which describes the diffusion of remdesivir across the cell membrane and conversion of the parent drug remdesivir to its active metabolite nucleoside TP inside the cell. The compartmental PK model was parameterized by fitting to published plasma profiles of remdesivir under different i.v. infusion schemes (Figure ). The intracellular metabolism model was developed by fitting to the time course of intracellular TP concentrations after incubating 1 μM remdesivir with human airway epithelial (HAE) cultures in vitro (Figure ). Because the human lung TP time course data are not available (and unlikely to be obtainable), we used the full remdesivir pharmacological model (including the compartment PK model and intracellular metabolism model) to predict TP concentrations in a monkey's lungs following a single i.v. (slow bolus) dose. Our model predicted a similar trend as monkey experimental data, where after initial i.v. infusion and plasma distribution, remdesivir diffused from the extracellular compartment into the intracellular compartment and was converted to the active metabolite TP, leading to a fast increase of intracellular TP concentration, followed by a gradual decrease of TP due to intracellular degradation (Figure ).

Figure 4

(a) Calibrating of the pharmacokinetic (PK) model. Fitting remdesivir plasma concentration after 2 hours of intravenous administration of 6 different single loading doses (3, 10, 30, 75, 150, and 225 mg). Black line: best fitted curve. Gray band: predicted population, points with error bars: clinical data. (b) Intracellular nucleoside triphosphate (TP) concentration following in vitro incubation with the parent drug remdesivir. Black points: the intracellular TP concentration in primary human airway epithelial (HAE) cultures after incubating with 1 μM remdesivir. Black line: fitted curve based on intracellular metabolism model, gray band: predicted population.

Calculation of in vitro SARS‐CoV‐2 inhibition potency of remdesivir's active metabolite

Because remdesivir needs to be converted to its active metabolite TP in vivo to inhibit SARS‐CoV‐2 replication, , it is conceivable that the in vitro virus inhibition potency (EC50) measured for the parent drug remdesivir, as generally reported in literature, , , may not reflect the drug's in vivo potency. Indeed, it has been speculated that the vastly different SARS‐CoV‐2 EC50 for remdesivir as measured in different in vitro cell cultures (Vero E6, Calu3 2B4, and HAE Cultures) may be due to various tissues' different capabilities of converting remdesivir to TP inside the cell. We hypothesized that, if the potency was calculated based on intracellular concentrations of TP, the EC50s may better reflect the characteristic potency of the drug and be more consistent among different cell cultures. To test this hypothesis, we constructed intracellular metabolism models for three cell cultures (Vero E6, Calu3 2B4, and HAE), and then used the predicted intracellular TP concentrations, instead of the nominal extracellular remdesivir concentrations, to recalculate the virus inhibition potency (Supplemental Methods, Figure ). The resulting EC50 values (Table ) are indeed much closer among different cell cultures, suggesting the tissue‐dependent intracellular metabolism capability may explain up to 90% of the apparent variability in SARS‐CoV‐2 inhibition potency among different cells/tissues in vitro. Because HAE cultures are widely considered the in vitro cell culture that is closest to in vivo human lung environment, we used the EC50 calculated by HAE intracellular TP concentration (0.01 μM, with a Hill coefficient 0.64, see Supplementary Methods for calculation) as the in vivo inhibition potency for the model (Figure ).

Table 1

In vitro EC50 variability across cell cultures for remdesivir vs. its active metabolite nucleoside TP

Remdesivir				Nucleoside TP
	Vero E6	Calu3	HAE		Vero E6	Calu3	HAE
Plaque assay EC50 (μM)	1.65	0.28	0.01	Plaque assay EC50 (pmol/million cells)	1.46	0.46	0.1
Ratio	165x	28x	1x	Ratio	14.6x	4.6x	1x
Genome copy EC50 (μM)	1.49	0.6	NA	Genome copy EC50 (pmol/million cells)	1.17	1.24	NA
Ratio	2.48x	1x		Ratio	0.94x	1x

Using plaque assays to measure efficacy, the largest cross‐culture differences (EC50 fold change) is 165 times if using remdesivir concentrations to calculate EC50s, compared to 14.6 times if using triphosphate metabolites concentrations. Using genome copy assays, the cross‐culture differences (EC50 in Vero E6 cells/EC50s in Calu3) is 2.48 times if using remdesivir concentrations to calculate EC50s, compared with 0.94 times if using triphosphate metabolites concentrations.

Calu3, 2B4 human lung adenocarcinoma cells; EC50, half‐maximal effective concentration; HAE, primary human airway epithelial cell culture; NA, not available; TP, triphosphate; Vero E6, African green monkey kidney cells.

Clinical endpoint calibration

Next, we calibrated the model to reproduce primary clinical end points used by clinical trials testing COVID‐19 therapeutics. We selected time‐dependent recovery percentage, which was used as the primary end point by a double‐blind, randomized, placebo‐controlled trial of i.v. remdesivir in adults hospitalized with COVID‐19 (ACTT). Importantly, only the placebo arm's data were used to calibrate the model to reproduce the disease progression profiles in patients with severe and mild disease. The remdesivir arm's data were left aside for independent validation of the full model (see the next section about model validation).

Because the remdesivir trial (ACTT) defined recovery as hospital discharge or hospitalization without requiring supplemental oxygen or ongoing medical care, we investigated what pathophysiological conditions could differentiate hospitalized vs. non‐hospitalized patients with COVID‐19. Interestingly, quantitative computed tomography scan and 3D reconstruction of the lungs in patients with COVID‐19 revealed that the mean lesion percent relative to the whole lung volume is ~ 4–7% in non‐severe or early‐stage disease patient groups. , Accordingly, we operationally defined those patients whose lung epithelial cell damage percent (as predicted by the disease model) was reduced below 5% (relative to the total number of lung epithelial cells in the model) as “recovered” patients. We calibrated the simulated patient populations so that initially (at day 0 in Figure ) both severe and mild disease patient populations have no patient in the “recovered” status (all simulated patients' lung epithelial cell damage percentage > 5% initially; Figure ). As time goes by, there is a steady increase in percentage of patients in the “recovered” status, with the time‐dependent recovery curve quantitively matching the clinically observed curve (Figure ). Importantly, such an increase of recovered patient percentage reflects the natural progression of the disease, (virus being eliminated by immune responses and tissue recovered from inflammation damage), but not the effects of remdesivir, as only the placebo arm's data are being used for model calibration.

Figure 5

Calibrating (a, b) and validating (c, d) the model for the primary endpoint (time to recovery) used in the remdesivir trial. For model calibration, placebo arm time‐dependent percentage of recovered patients in the mild (a) and severe (b) disease groups were used to adjust the model parameters. For model validation, the calibrated model was used to independently predict remdesivir arm (RMD) data for mild (c) and severe (b) patients without any parameter adjustments. Of note in the Adaptive Covid‐19 Clinical Trial (ACTT) study mild/moderate disease was defined by a SpO2 > 94% and respiratory rate < 24 breaths per minute without supplemental oxygen requirement, whereas severe disease was defined as meeting one of the following criteria: requiring invasive or noninvasive mechanical ventilation, requiring supplemental oxygen, an SpO2 ≤ 94% on room air, or tachypnea (respiratory rate ≥ 24 breaths per minute). These definitions are generally aligned with other studies and also aligned with our virtual populations of mild and severe disease. Black points: Kaplan–Meier estimation of clinical data. Error bar: 95% confidence interval of clinical data. Black line: simulated curve based on subjects in model population. Gray band: predicted uncertainty quantification. X axis: time (in days) since the start of the trial. Y axis: percentage of recovered patients. Day 0, is the day that subjects were admitted into the clinical trial. Table shows the median recovery day and 95% confidence interval.

As expected, the increase of recovered patient percentage is faster in the mild disease patient population (median time to recovery 5 days in clinical data and 6 days in model simulation) vs. severe disease patient population (median time to recovery 18 days in both clinical data and model simulation). Of note, the ACTT study stratified patients into a moderate‐to‐mild stratum and a severe stratum, using definitions generally aligned with other studies, and hence also aligned with our virtual populations of patients with mild and severe disease.

In addition to the moderate‐to‐mild vs. severe stratum, the ACTT study also did a more granular stratification according to patients' baseline ordinal clinical score. We calibrated separate virtual populations to match the placebo arm recovery time from each stratum (Figure ). When the baseline ordinal score is 4, 5, 6, and 7, the simulated median time to recovery is 6 (6 in clinical data), 10 (9 in clinical data), 21 (20 in clinical data), and 26 (28 in clinical data) days, respectively.

Model validation by predicting remdesivir efficacy end point

Finally, the full model (with the pharmacological model for remdesivir and the disease model for COVID‐19) was combined to predict the primary efficacy end point (time‐to‐recovery) of the remdesivir trial ACTT. Even though the placebo arm of this trial was used to calibrate the model (see the previous section about clinical end point calibration), the remdesivir arm data were never used during model calibration. We ran the full model by using the same dosing scheme as the remdesivir trial, where remdesivir was administered intravenously as a 200‐mg loading dose on day 1, followed by a 100‐mg maintenance dose administered daily on days 2 through 10. Our predicted curve of time‐dependent increase of recovered patients quantitatively agrees with the clinical data in the remdesivir arm for both severe (median time to recovery 11 days in both clinical data and model prediction) and mild populations (median time to recovery 5 days in clinical data and 6 days in model prediction; Figure ). For the more granular stratification, the model was able to predict the time to recovery in the remdesivir group quantitatively as well (Figure ). When the baseline ordinal score is 4, 5, 6, and 7, the predicted median time to recovery is 6 (5 in clinical data), 8 (7 in clinical data), 17 (15 in clinical data), and 28 (29 in clinical data) days, respectively. It is worthwhile to reiterate that we did not adjust any parameters or model structures based on remdesivir arm data, and hence the model validation is completely independent.

DISCUSSION

We have developed a mechanistic in silico COVID‐19 model aimed at predicting clinical outcomes based on nonclinical data and shown through a proof‐of‐concept example that this model can independently and quantitatively translate nonclinical pharmacology data into the clinical endpoint of time to recovery of the remdesivir clinical trial ACTT. Using a modularity concept, the whole model can be divided into smaller modules or submodels (virus life cycle, lung immune responses, lymph node immune responses, pharmacology, etc.) to facilitate model development and calibration. The individual modules can be reconnected to form model variants or reused in new models, such as to link other classes of COVID‐19 therapeutics to the disease (sub)model for clinical outcome prediction, making this an attractive strategy for quantitative systems pharmacology modeling.

In the disease model part, viral life cycle and spread of the virus among cell populations are supplemented by incorporating the innate and adaptive immune response to the viral kinetics of SARS‐CoV‐2 in the lung and lymphatic compartments. The disease model parameters were calibrated to reproduce key characteristics of longitudinal clinical data from patients with mild and severe COVID‐19. Although there is no universal definition for patients with mild vs. severe COVID‐19, most clinical studies , , separated these 2 subpopulations due to their different clinical presentations, with the need for supplemental oxygen therapy as one of the most prominent distinguishing factors. , Despite the simplifications made in the model, such as using IL‐6 as a surrogate for all pro‐inflammatory cytokines and lumping or ignoring complex immune responses (e.g., antigen cross‐presentation and immune memory), the disease model captures typical clinical responses to mild and severe SARS‐CoV‐2 infections. The inclusion of human immune responses and simultaneous fitting to a series of interacting clinical variables (viral titer, cytokines, antibodies, lymphocytes reduction, etc.) are key differences between our approach and some other work modeling viral kinetics. , These models use viral titers as the only clinical output to track disease progression, and usually assume viral titers will not decline until all target cells are consumed. This contrasts with the recognized effects of immune systems on both virus clearance , and (with overactive immune response) disease deterioration. As an illustrative example, our model was able to show that, for some patients, disease severity is not positively correlated with the initial viral load (Figure ), which has also been reported clinically.

Among the many drugs being researched, we have focused on remdesivir, the only approved COVID‐19 therapy by the FDA, because of the availability of comprehensive information in different phases of studies (in vitro, in vivo, and human clinical trial data). The substantial difference of in vitro EC50s across cell cultures between remdesivir as a prodrug and its active metabolite nucleoside TP led us to design a pharmacological model to capture the transition pathway from the prodrug to active metabolite TP inside cells. We found that the dosing scheme used in ACTT achieved a predicted lung intracellular TP concentration fluctuating between ~ 0.07 and ~ 0.35 μM for an average patient during maintenance dosing (data not shown), resulting in an in vivo inhibition percentage between 77% and 90%. This is consistent with the observed remdesivir efficacy in ACTT, and also raises an interesting question whether alternative dosing schemes would decrease the time‐to‐recovery even further. One limitation of our current study is that PK drug–drug‐interactions is not considered, as some participants in ACTT received hydroxychloroquine, which was shown in in vitro studies to have the potential of affecting PKs of remdesivir. However, sensitivity analyses performed by ACTT suggested that this may not affect the efficacy of remdesivir in a clinical setting.

One hallmark of our modeling framework is to use (sub)populations, instead of single “average” patients, to represent the disease status of patients with mild and severe COVID‐19, respectively. This helps to capture the heterogeneity within and between the 2 subpopulations. The mild and severe subpopulations are calibrated to reproduce the clinical time courses of various biomarkers (virus titer, plasma IL‐6, etc.) from different studies, , , , making them useful to simulate the baseline disease progression for COVID‐19 (e.g., forming the basis to mimic the placebo group in a clinical study for a specific COVID‐19 treatment), although parameters may have to be re‐evaluated for different variants. In addition, the model's capability of predicting efficacy end points from the remdesivir trial suggests it might have credibility to predict other clinical scenarios for remdesivir, such as other dosing routes, alternative dosing schemes, etc. On the other hand, whether such credibility can be generalized to predicting clinical outcomes of other COVID‐19 treatments awaits further investigation.

In summary, the mechanistic prediction framework presented in this work is the first proof‐of‐concept model development example to illustrate that it is possible to quantitatively predict clinical outcomes based on nonclinical (in vitro) data for COVID‐19 therapies. This paves the way for quickly selecting candidate drugs that have shown promising results in the early nonclinical stages of drug discovery and can facilitate the design of clinical studies accordingly. The modular design of the model also makes it possible to simulate drug combination effects based on monotherapy data, without the need to systematically recalibrate the disease model. Such a mechanism‐based, model‐informed drug development tool may potentially play an important role in our urgent pursuit of safe and effective COVID‐19 therapies.

FUNDING

This project was supported by the Research Participation Program at CDER, administered by the Oak Ridge Institute for Science and Education (ORISE) through an interagency agreement between the US Department of Energy and the FDA.

CONFLICT OF INTEREST

The authors declared no competing interests for this work. This report is not an official US Food and Drug Administration guidance or policy statement. No official support or endorsement by the US Food and Drug Administration is intended or should be inferred. As an Associate Editor for Clinical Pharmacology & Therapeutics, David G. Strauss was not involved in the review or decision process for this paper.

AUTHOR CONTRIBUTIONS

Z.L. and M.S. wrote the manuscript. Z.L., M.S., J.L.W., K.E.H., R.R., J.F., and D.G.S. designed the research. Z.L. and M.S. performed the research. Z.L., M.S., A.C., J.M., J.Z., X.H., and G.A. analyzed the data.

Appendix S1

Click here for additional data file.