Abstract
Molnupiravir is an antiviral medicine that induces lethal copying errors into SARS-CoV-2 RNA. Molnupiravir reduced hospitalization in one pivotal trial by 50% and had variable effects on reducing viral RNA levels in three separate trials. We used mathematical models to simulate these trials and closely recapitulated their virologic outcomes. Model simulations suggest lower antiviral potency against pre-omicron SARS-CoV-2 variants than against omicron. We estimate that in vitro assays underestimate in vivo potency 7-8 fold against omicron variants. Our model suggests that because polymerase chain reaction detects molnupiravir mutated variants, the true reduction in genetically intact viral RNA is underestimated by ∼0.5 log in the two trials, which included participants with omicron. Our results reinforce past work suggesting that in vitro assays are unreliable for estimating in vivo antiviral drug potency and suggest that virologic endpoints for respiratory virus clinical trials should be catered to the drug mechanism of action.
Introduction
Molnupiravir is an antiviral drug that induces errors in the SARS-CoV-2 genome, which typically renders the virus unable to replicate further.1 In the randomized double-blinded MOVe-OUT trial, which enrolled unvaccinated individuals when delta, mu, and gamma variants of concern (VOC) were circulating, molnupiravir reduced hospitalization by 50% and viral load after treatment by 0.3 log relative to placebo.2 In the platform adaptive PANORAMIC trial, which enrolled vaccinated individuals when omicron VOCs were circulating, molnupiravir did not lower hospitalization rates, which were only 1% in both arms, but lowered viral load after treatment by 0.94 log relative to usual care.3 In the platform adaptive PLATCOV trial, which enrolled low-risk individuals when omicron VOCs were circulating, molnupiravir lowered viral load after treatment by 1.09 log relative to ususal care.4 Taken together, these trials demonstrate that molnupiravir has both clinical and virologic efficacy which varied across trials and viral variants.
Overall, use of molnupiravir has been lower than that of nirmatrelvir / ritonavir based on lower reduction in hospitalization in MOVe-OUT relative to the EPIC-HR trial for nirmatrelvir / ritonavir.5 A concern has also been raised that molnupiravir’s mechanism of action could generate novel mutants that persist after cessation of treatment,6 and then spread in the population.7 Nevertheless, PANORAMIC and PLATCOV results suggest high potency, and molnupiravir is still considered in individuals in whom nirmatrelvir / ritonavir is contraindicated, and in combination with other drugs in immunocompromised hosts.8, 9 There is currently no explanation for the disparate antiviral effects in MOVe-OUT versus PANORAMIC and PLATCOV. Moreover, the fact that polymerase chain reaction (PCR) detects drug-altered viral RNA molecules6 has not always been considered in the analysis of trial outcomes.
We previously used clinical trial simulation to reproduce results from nirmatrelvir / ritonavir trials for SARS-CoV-2.10, 11 We first validated a viral immune dynamic model (VID) against a very large prospective cohort of infections that included multiple VOCs.11 We used diverse virologic output from this model to create simulated cohorts for the control arms of trials. We then layered pharmacokinetic (PK) and pharmacodynamic (PD) models for nirmatrelvir / ritonavir on the VID models to simulate treatment arms.10 This approach recapitulated mean viral load reduction in the EPIC-HR and PLATCOV trials, as well as individual viral load trajectories in PLATCOV. The validated model was then used to explain the high frequency of virologic and concurrent symptomatic rebound with use in the community,12 despite very low levels of virologic and symptomatic rebound in the EPIC-HR trial.13, 14 Model output suggests that extending therapy from 5 to 10 days would nearly eliminate rebound,10 a result confirmed with modeling of separate data by another group.15
Here we expand this approach to develop a new joint VID-PK-PD model to account for the unique mechanism of action of molnupiravir. We simulate by fitting the model to results from the MOVe-OUT, PLATCOV and PANORAMIC trials. Our results suggest that quantitative viral PCR likely underestimates the reduction in non-mutated viral RNA and therefore the true potency of molnupiravir during omicron infections.
Results
Viral immune dynamic, pharmacokinetic and pharmacodynamic clinical trial simulation models
We previously described our viral immune dynamic (VID) model that was fit to diverse viral loads from 1510 SARS-CoV-2 infected individuals in the National Basketball Association cohort.11, 16 The model assumes a finite number of susceptible cells. An eclipse phase delays viral production by infected cells. In keeping with an early innate immune response, susceptible cells become refractory to infection in the presence of infected cells but also revert to a susceptible state at a constant rate. Infected cells are cleared by cytolysis and delayed acquired immunity, which is activated in a time-dependent fashion (Fig 1a). We used a mixed-effect population approach implemented in Monolix to estimate model parameters.17
To reproduce levels of molnupiravir, we used a three-compartment pharmacokinetic (PK) model (Fig. 1b). Using Monolix and the mixed-effect population approach, we estimated parameter values by fitting the model to the average plasma concentration of healthy subjects.18 The model closely recapitulated observed drug levels following multiple doses of 50, 100, 200, 300, 400, 600, and 800 mg given twice daily for five days (Fig S1, Table S1). The estimated value for the transition rate from plasma to peripheral compartment (was dose-dependent in the form of. All other PK parameters were dose independent. For the pharmacodynamic (PD) model, we assumed drug efficacy follows a Hill equation with respect to concentration. We parameterized the model using in vitro efficacy data collected at different concentrations (details in Materials and Methods, Fig S2, Table S2).19
We combined the VID and PKPD models by using treatment efficacy to convert non-mutated virus to mutated virus, both of which are assumed to be detected with polymerase chain reaction (PCR) assays, given the low probability of drug-induced mutations in the PCR primer region. This assumption is based on the observed drug-induced mutation of approximately 1 in 2000 sites.20 Given the average length of most PCR primers of ∼25 base pairs, the chance of the primer remaining unmutated after treatment is (1999/2000)25, or 98.76%. A limitation of viral load data from the included clinical trials is that it lacks early pre-symptomatic endpoints to estimate the viral expansion slope. To further train the model, we included 1023 Omicron-infected participants from the NBA cohort to inform rates of viral upslope in the trials.11 We first fit the combined model to individual viral load data from 149 low-risk, symptomatic vaccinated participants infected with omicron VOCs in the PLATCOV trial (65 treated and 84 controls) and from 80 high-risk, symptomatic vaccinated participants infected with omicron VOCs in the PANORAMIC trial (38 treated and 42 controls) (Fig 2a, 2b, Fig S3, S4, S5, S6). We next fit the combined model to trial endpoint data (mean viral load drop from baseline) reported in three randomized, controlled trials: PLATCOV (Fig 3a, 4a-d),4 PANORAMIC (Fig 3b, 4e-h),3 and the MOVe-OUT trial with 1093 high-risk unvaccinated symptomatic individuals infected with pre-omicron VOCs (549 treated + 544 placebo, Fig 3c, 4i-k).2 All model fitting was using Monolix with non-linear mixed effects approaches described in the Materials and Methods.
Model fitting to individual viral load trajectories in PLATCOV and PANORAMIC
For each participant, we defined the in vivo EC50 as the plasma drug concentration required to inhibit viral replication by 50% and the potency adjustment factor (paf) as the ratio between the in vivo and in vitro EC50.21, 22 To estimate the paf, we fit the combined VID-PKPD model to individual viral load data of both arms of PLATCOV and PANORAMIC trials and Omicron infected individuals in the NBA cohort using the population mixed effect approach in Monolix. We achieved good model fit to individual viral load trajectories in the control and treatment arms of PLATCOV (Fig 2a, Fig S3, S4) and PANORAMIC (Fig 2b, S5, S6). The model projected higher levels of total detected SARS-CoV-2 RNA in most participants relative to non-mutated viral RNA (Fig 2a,b). We estimated a range of individual paf values with similar mean and median values estimated for both trials (Fig 2c). These values suggest that in vivo potency of molunpiravair is on average 7-8 fold higher than in in vitro estimates. Each participant had an estimated paf less than one indicating that enhanced potency is necessary to include to optimize fit to the data.
Model fit to trial virologic endpoint data from PLATCOV, PANORAMIC, and MOVe-OUT
As a second approach, we assessed whether a virtual cohort strategy where control participants are modeled from pre-existing cohorts can predict virologic trial endpoints. This approach is necessary for situations where individual viral load data are not available as with MOVe-OUT and is important as it demonstrates that the model can reproduce the primary virologic endpoint of the study. We simulated virtual cohorts using the combined VID-PKPD model and fit results to viral load decay from baseline in the three trials. For each trial arm, we randomly selected 400 individuals from the NBA cohort with the closest matching viral variant, symptom, and vaccine status and used their estimated individual viral dynamic parameters in simulations. To address variability in timing of baseline viral load measurement relative to infection, we randomly assigned all individuals an incubation period selected from a variant-specific gamma distribution found in the literature.23, 24 Treatment start day was randomly selected from a distribution based on observed enrollment windows in the three trials. Due to the lack of individual PK data, the same estimated population PK parameters were used for all simulated treated individuals (Table S1). PD parameters were randomly selected from a log-normal distribution with estimated mean and standard error (Table S2).
Our model closely reproduced kinetics of viral decay in PLATCOV in control (Fig 3a, 4a) and treatment arms (Fig 3a, 4b) and estimated a paf=0.13 (Fig 4c) similar to our mean estimate using individual fits (Fig 2c). The model also predicted individual-level variability in virologic responses observed in PLATCOV, including instances of increased viral load following therapy (Fig 4d). We compared simulated and actual distributions of viral load change among trial participants in the control and treatment arms. On most post-treatment days, simulated and actual distributions were not statistically dissimilar. Wider distributions of observed versus simulated viral load change were noted on post-randomization days 1 and 2 for control and treatment (Fig 4d), likely due to noise in viral load data from oral swabs: differences of 1–2 logs were often noted between replicates collected from PLATCOV participants at equivalent timepoints, particularly on days 1 and 2.10
Similarly, our model closely reproduced kinetics of viral decay in PANORAMIC in control (Fig 3b, 4e) and treatment arms (Fig 3b, 4f) and estimated a paf=0.19 (Fig 4g) similar to our mean estimate using individual fits (Fig 2c). The model also predicted individual-level variability in virologic responses observed in PANORAMIC, including instances of increased viral load following therapy (Fig 4h). We compared simulated and actual distributions of viral load change among trial participants in the control and treatment arms. On all post-treatment days other than day 2 control, simulated and actual distributions were not statistically dissimilar (Fig 4h). This likely indicates less noise in viral load data from nasopharyngeal swabs collected in PANORAMIC relative to oral swabs in PLATCOV.
Finally, the model reproduced kinetics of viral decay in MOVe-OUT in control (Fig 4i) and treatment arms (Fig 4j) but estimated a higher paf=2.64 (Fig 4k). The higher paf maps to the far less substantial viral load reduction in MOVe-OUT relative to the other two trials which in turn may be explained by less potency against pre-omicron variants which has been observed experimentally.25
As a further validation step, we performed counterfactual simulations which assess viral kinetics of control study participants assuming treatment and treatment participants assuming placebo/usual care. Counterfactual control simulations slightly overestimated late viral loads for PLATCOV (Fig S7a) and PANORAMIC (Fig S8a). This may be because therapy suppresses acquired immune responses, which is not captured in our model.15 Counterfactual treatment simulations fit the data well for PLATCOV (Fig S7b) and PANORAMIC (Fig S8b). Simulations occasionally predicted viral rebound following treatment (Fig S7c,d and S8c,d).
Estimates of reduction in fully mutated viruses versus non-mutated SARS-CoV-2 RNA
We used our optimized model with solved paf to project the trajectory of non-mutated viral RNA during treatment relative to values measured with PCR, which detects viral RNA with drug-induced mutations.6 In PLATCOV (Fig 5a,d,f) and PANORAMIC (Fig 5b,d,f), owing to higher drug potency, total SARS-CoV-2 viral RNA on treatment exceeded non-mutated viral RNA by ∼0.48 and 0.59 log respectively, suggesting that measured endpoints underestimate the drug’s true antiviral effect. However, these differences did not achieve statistical significance, perhaps because estimated total and non-mutated viral RNA levels converge at drug trough. In MOVe-OUT (Fig 5c,d,f), there was no significant difference between total SARS-CoV-2 viral RNA on treatment and non-mutated viral RNA.
In all 3 trials, the models suggest that non-mutated viral loads during treatment may be lowest at drug peak, reflecting the short plasma half-life of molnupiravir. Therefore, the cumulative effective of drug is best estimated with viral area under the curve which accounts for highly variable drug activity over time due to short drug half-life. By this estimate, the reduction in non-mutated viral RNA far exceeds that of total measured viral RNA (Fig 5e). In the case of MOVe-OUT, this may explain why significant clinical benefit is associated with only a marginal reduction in observed decline in total viral RNA. This also suggests that there might be utility to measure viral load after drug peak and trough for agents with short half-life as viral loads may differ by a full order of magnitude according to drug level.
Differences in trial participants and model parameters
We also compared features of each trial as they related to model predictions by assessing the viral dynamic range in each trial (Fig S9a). Control participants in PLATCOV had lower mean viral loads throughout the course of infection relative to PANORAMIC and MOVe-OUT (Fig S9b). Given PLATCOV and PANORAMIC enrolled participants with the omicron variant, we surmise that these differences relate to demographic differences in study participant demographics, slightly shorter estimated time to treatment in PANORAMIC versus PLATCOV estimated by our model (t0 in Fig S10), and/or characteristics of the PCR assay used in the studies. The trials also employed different limits of detection which impacted observed reductions in viral load (Fig S9b). Model parameters were largely equivalent between studies and between treatment and control arms, reflecting the flexibility of the model (Fig S10). The parameter governing transition of susceptible cells to a refractory state was higher in PLATCOV relative to PANORAMIC which likely was necessary to achieve lower viral loads overall and may reflect the younger study participants (Fig S10).
Antiviral potency, viral load assessment and trial design all impact observed antiviral reduction
We next combined PK and PD models to assess the average efficacy of the drug during days 0-5 in all three trials (Fig 6a,b) and noted an efficacy of 53% in MOVe-OUT (Fig 6c). The efficacy of molnupiravir in PLATCOV (94%) and PANORAMIC (95%) was similar to that of nirmatrelvir in PLATCOV (94%) and EPIC-HR (82%) owing to a much lower paf for molnupiravir relative to nirmatrelvir / ritonavir (Fig 6c).
These potencies mapped to different reductions in viral load relative to placebo/usual care. Total viral RNA reduction in PANORMIC and PLATCOV exceeded that in MOVe-OUT owing to lower paf (Fig 6d), but also due to a larger viral dynamic range (defined as the distance from baseline viral load to the lower threshold PCR (LOD) (Fig S9b)), which allows for a greater observed reduction in viral load. PANORAMIC used LOD = 100 (imputed as 50), and PLATCOV had LOD≌18 (copies/ml), versus LOD = 500 copies/ml for MOVe-OUT. PANORAMIC also had much higher average starting viral loads (7.4 log10 (copies/ml)) versus PLATCOV (5.8 log10 copies/ml) and MOVe-OUT (6.8 log10 copies/ml) (Fig S9b). Molnupiravir approached maximal possible total viral load viral RNA reduction in PANORAMIC and PLATCOV, whereas protease inhibitors could still achieve greater viral load reduction at lower paf (Fig 6d) as recently observed with ensitrelvir.26
The greater possible reduction in total viral load for protease inhibitors relative to molnupiravir owes to different mechanisms of action. Model projected reductions in non-mutated viral RNA reduction in PANORMIC and PLATCOV approximated viral load reductions observed in EPIC-HR and PLATCOV on nirmatrelvir / ritonavir (Fig 6e), suggesting that PCR detection of mutated viruses underestimates true molnupiravir potency and that a more potent mutagenic agent could accrue further virologic benefit.
Optimization of molnupiravir therapy to avoid viral rebound
Instances of viral rebound were observed in PLATCOV and PANORAMIC and have been observed following molnupiravir treatment.4, 27 We analyzed higher doses and prolonged therapy and noted that as with nirmatrelvir, prolonging therapy is a better method to prevent rebound than increasing dose (Fig S11).
Discussion
We recapitulated the virologic results of three clinical trials of molnupiravir with our combined clinical trial simulation VID/PK/PD models. Model output highlights key differences in viral load reduction between the trials and identifies mechanistic differences to explain this. The MOVe-OUT trial was associated with significantly less reduction in viral load between treatment and control arms than the other two trials. Accordingly, the drug efficacy (0.53) was lower in this trial, and our estimate for paf was 2.64, signifying marginally lower potency in vivo than in vitro. This result is compatible with 0.4-fold lower median EC50 values for molnupiravir in vitro against omicron relative to prior variants. 25 The higher paf in MOVe-OUT permitted drug troughs below the EC50, limiting potency throughout the dosing interval. In PLATCOV and PANORAMIC, our model suggests drug levels remain above the in vivo EC50 throughout the dosing interval though potency does fluctuate according to drug level.
The estimated drug efficacy in PLATCOV and PANORAMIC was considerably higher (0.94 and 0.95, respectively), and the paf was estimated to be 0.14 and 0.13, indicating greater potency in vivo than in vitro. We have applied our clinical trial simulation technique to multiple drugs for SARS-CoV-2,10 HSV-2,21, 22 and HIV,28 and this is the first time we have identified this trend. Allowing molnupiravir to be more potent in vivo in our modeling of PLATCOV and PANORAMIC was necessary to capture the much greater reduction in total viral RNA relative to off-treatment in these studies (1.09 log and 0.94 log respectively versus 0.3 log in MOVe-OUT).
A key outcome of our analysis is that SARS-CoV-2 PCR likely underestimates molnupiravir potency because it still detects drug-mutated viral RNA.6 This appears to be most significant when antiviral potency is higher, as in PANORAMIC and PLATCOV, leading to a 0.48-0.59 log underestimation of reduction in non-mutated virus. Our results suggest that use of standard PCR for assessing SARS-CoV-2 levels may lead to underestimation of drug potency. Multi probe assays as have been used for the HIV reservoir may improve specificity for viruses that remain intact and replication competent.29 Viral culture is potentially a useful metric but lacks sufficient sensitivity and precision and is too labor intensive to serve as a viable trial endpoint.12
Our results suggest that viral loads may vary according to drug level given molnupiravir’s short plasma half-life. A sub-study within future trials comparing viral loads between drug trough and peak would be useful for the field. This could validate our model’s prediction that even small reductions in viral RNA may be associated with substantial reductions in total viral area under the curve (a surrogate for infection surface area) particularly with a short half-life drug. Even minor reductions in viral load may be associated with substantial clinical benefit in this case. These fluctuations may be less evident if the intra-cellular half-life of the drug is longer or if PK measures in our model underestimate true drug levels in PANORMIC due to older age or impaired renal clearance. Model accuracy would be improved if viral load and PK data were available from the same patients in similar trials.
Another key practical outcome is that as for nirmatrelvir, extension of molnupiravir therapy to 10 days is likely to prevent rebound though our simulations do not suggest any benefit from increasing dose.10 This suggests that prolonging therapy or using longer half-life agents is ideal for treating SARS-CoV-2.30
Each trial represented a unique set of issues for model fitting. In MOVe-OUT, because the treatment arm mean viral kinetics curve differed only slightly from the control arm, the model without drug provided reasonable fit to the treatment arm. Nevertheless, the paf was identifiable for this trial indicating that the model was able to detect and specify the very limited potency of the drug. The fact that the drug’s potency and clinical efficacy appears to have increased with introduction of the omicron variant demonstrates a massive challenge for the therapeutics field: as with vaccines, trials performed when prior variants were circulating may prove less relevant as new variants continually emerge. A priority should be retesting existing agents against newly emerging variants in small nimble trials such as PLATCOV, with viral load endpoints.
For PLATCOV, the model for the treatment arm matched the trial data very precisely and identified the paf. The drug achieved nearly maximal observed viral load reduction in this trial. We identified a similar trend for PANORAMIC. It is notable that the model had the flexibility to account for different viral loads between these trials by predicting more rapid innate immune responses in PLATCOV.
We arrived at similar estimate for paf in PLATCOV and PANORAMIC, while using two separate methods: fitting to individual viral loads and fitting to mean viral load reduction trial endpoints. This suggests that our approach using in silico cohorts and fitting to population level outcomes produces reliable results.10 In many cases, it remains challenging for academics to obtain individualized data from industry sponsored trials. Therefore, it is important that the endpoint fitting approach should be considered when this is the only data available.
Finally, our results highlight challenges in trial design associated with the selection of PCR assays and limits of detection. Different limits of detection were selected for each trial which in turn impacts the degree of viral load decrease that can be observed. Initial viral loads were notably higher in PANORAMIC and MOVe-OUT than PLATCOV. The equivalent viral loads between PANORAMIC and PLATCOV may reflect a more sensitive PCR in the PANORAMIC study as past immunity has consistently predicted lower viral loads and more rapid viral elimination for omicron variants.11, 31 On the other hand, PANORMIC viral loads could have been higher than in PLATCOV despite both enrolling omicron infections due to an older and less healthy population, Ideally, equivalent internationally standardized PCR quantitation would be used across all trials.
Our study has a couple of limitations. The estimated paf is based on the in vitro assay data against delta variant in Calu-3 cells. In vitro EC50 is sensitive to assay conditions, including cell type, the variant of concern, the multiplicity of infection (MOI), and specific lab. In general, the inability of in vitro assays to match in vivo conditions makes it an unreliable proxy for the in vivo potency of the drug and explains the necessity of incorporating the paf parameter when simulating an antiviral clinical trial.
In our model, we assumed all mutated viruses are noninfectious. While most drug-induced random mutations reduce the fitness/viability of the virus, the molnupiravir mutation signature has been detected in circulating variants, especially in regions where the drug has been used most commonly.7 Furthermore, in the PANORAMIC study, no statistical difference was observed between the rate of culturable samples in the treatment and control arms which implies some mutated viruses might be transmissible.32
PK parameters were estimated using the data from the plasma concentration of healthy individuals with the age range of 19-60 (mean 39.6) years old. The clearance rate of renally cleared drugs often increase with age.33 This implies that the paf may be larger in an older population, such as in PANORAMIC participants. Further, we used the plasma concentration in the PD model to calculate the drug efficacy. However, using the drug’s intracellular concentration with a longer half-life, represented by the central compartment of the PK model, would also likely lead to a larger estimated paf.
In summary, we further demonstrate the utility of clinical trial simulation using models that capture drug PK and PD, as well as infection dynamics. In the case of molnupiravir, our results suggest that final viral endpoints should be adjusted based on the drug’s mechanism of action.
Materials and Methods
Study Design
We developed a viral dynamics model recapitulating viral load data collected from symptomatic individuals in the NBA (National Basketball Association) cohort.11 We used a two-compartment model to reproduce the PK data of molnupiravir.18 For clinical trial simulation, we constructed a virtual cohort by randomly selecting 400 individuals from the NBA cohort, trying to match the trial populations regarding vaccine status and history of infection. Due to lack of individual PK data, we used estimated population PK parameters for all individuals in the virtual cohort. The PD parameters for each individual were randomly selected from a log-normal distribution with their estimated means and standard errors. We fit the combined viral dynamic and PKPD model to the average change in viral load from the baseline of the control and treatment arms of three previously published molnupiravir clinical trials.2-4, 34 Comparing our model to the control arms validated our viral dynamic model and demonstrated how well our virtual cohorts represent the trial control arms. We used the average data from the treatment arms to estimate the potency adjustment factor (paf) by maximizing the R2 of the fit.10 We also fit to individual viral load trajectories in PLATCOV and PANORAMIC using the mixed-effect population approach implemented in Monolix35-37 and obtained both individual paf values and a population distribution.
Viral load data
The NBA cohort dataset published by Hay et al consists of 2875 documented SARS-CoV-2 infections in 2678 people detected through frequent PCR testing regardless of symptoms.16 We used the viral load data from 1510 infections in 1440 individuals that had at least 4 positive quantitative samples to estimate viral dynamic parameters. We used parameter sets estimated for the symptomatic subpopulation of these individuals to construct virtual cohorts.11
Clinical trial data
We used viral load data from three molnupiravir clinical trials. MOVe-OUT by Jayk Bernal et al. included 544 and 549 symptomatic high-risk individuals in the control and treatment arms, respectively.2 We obtained the average change in viral load data of the control and treatment arms as shared in the supplementary material of their published manuscript in Table S6. Nasal viral load was measured using PCR assay on days 0, 3, 5, 10, and 14 after the treatment start day and adjusted by the baseline viral load. The lower limit of detection (LOD) of 500 copies/ml were used in this trial. The treatment was started within five days of symptom onset.
PLATCOV by Schilling et al. is an open-label, randomized, controlled adaptive trial with 85 and 58 symptomatic, young, healthy individuals in the control and molnupiravir treatment arms, respectively.4 The oropharyngeal samples from each participant were collected daily on days 0 through 7 and on day 14 after the treatment start day, and viral load was measured using PCR assay. We used the individual viral load data published by the authors. From PLATCOV, we averaged over the two oral samples collected from each individual and calculated viral load drop from baseline or used the individual-level viral load data. We used the maximum LOD reported in the published data (∼1.26 log).
PANORAMIC is a platform adaptive randomized, controlled trial with 42 and 38 symptomatic, vaccinated individuals with at least one risk factor in the control and treatment arms, respectively.3 Nasal viral load was measured using PCR and samples were collected on days 0 through 6 and on day 13. We used individual viral load data shared by the authors and adjusted by baseline viral load to obtain mean viral load drop from the baseline. Mean days since symptom onset at baseline (sd) were 2.4 (0.78) for the treatment arm and 2.5(1.12) for the control arm. The lower limit of detection of 100 copies/ml (imputed as 50 copies/ml) was used.
In all three trials, the study participants were treated with 800mg molnupiravir twice per day, for five days.
PK data
Mean plasma concentration data of molnupiravir were obtained by digitizing Figure 3 of the phase I trial by Painter et al. using WebPlotDigitizer.18 The data used in this paper belong to the multiple-ascending part of the phase I trial where six participants were given 50, 100, 200, 300, 400, 600, and 800 mg of molnupiravir twice daily for 5.5 days, and their plasma concentrations were measured after the first and last doses.
PD data
The data on drug efficacy was obtained from experiments performed at the University of Washington. The efficacy of molnupiravir was measured against the delta variant of SARS-CoV2 in Calu-3 cells (human lung epithelial). Briefly, Calu 3 cells were treated with varying concentrations of molnupiravir prior to infection with SARS-CoV-2 (delta isolate) at a multiplicity of infection of 0.01. Antiviral efficacy and cell viability (of non-infected cells treated with drugs) were assessed as described.9, 10 There were five replicates per condition, pooled from 2 independent technical experimental repeats (one experiment with triplicate conditions, one experiment in duplicate conditions).
Viral dynamics model
We used our model of SARS-CoV-2 dynamics to model the viral load of symptomatic individuals with SARS-CoV-2 infection.11 Our model assumes that susceptible cells (S) are infected at rate βVS by SARS-CoV-2 virions. The infected cells go through a non-productive eclipse phase (IE) before producing viruses and transition to becoming productively infected (IP) at rate κIE. When encountering productively infected cells, the susceptible cells become refractory to infection (IR) at the rate ϕIPS. Refractory cells revert to a susceptible state at rate ρR. The productively infected cells produce virus at the rate πIP and are cleared at rate δI representing cytolysis and the innate immune response that lacks memory and is proportional to the amount of ongoing infection. If the infection persists longer than time τ, then cytotoxic acquired immunity is activated, which is represented in our model by the rate mIP. Finally, free virions are cleared at the rate γ. Of note, this model, previously proposed by Ke et al. was selected against other models based on superior fit to data and parsimony.38 The model written as a set of differential equations has the form, To estimate parameter values, we fit the model to viral load data from the NBA cohort using a mixed-effect population approach implemented in Monolix. Details on the model selection and fitting process can be found in Owens et al.11
We start the simulations with 107 susceptible cells. The initial value of the refractory cells is assumed to be zero since the interferon signaling is not active prior to infection. We further assume there are no infected cells (eclipse or productive) at the beginning of the infection. We fix the level of inoculum (V0) at 97 copies/ml for each individual.
To resolve identifiability issues, we fixed two parameter values, setting the inverse of the eclipse phase duration to κ = 4, and the rate of clearance of virions to γ = 15.11
PK model
We used a two-compartmental PK model which includes the amount of drug in the GI tract (AGI), the plasma compartment (Ap), and the respiratory tract (AL). The drug is administered orally, passes through the GI tract and gets absorbed into the blood at the rate κa. The drug then transfers from the blood into the peripheral compartment (or the lung) at the rate κPL. The metabolized drug transfers back into the plasma at the rate κLP from where it clears from the body at the rate κCL. The model in the form of ordinary differential equations is written as, We used Monolix and a mixed-effect population approach to estimate the parameters and their standard deviations. With the initial condition of (AGI = Dose, Ap = 0, AL = 0); we fit to the plasma concentration data where Vol is the estimated plasma volume. Details on parameter values and the error model provided in Table S1.
PD model
For the pharmacodynamics model we used Hill equation, , where is the drug concentration in plasma, Emax is the maximum efficacy, n is the Hill coefficient, and EC50 is the drug concentration in plasma required for 50% efficacy. We used least-squared fitting to obtain the three parameters (Emax, n, and EC50) and their standard deviations. The average drug efficacy is measured using, Where tstart and tend are the treatment start day and end day, respectively.
Combined PKPD and VL models
The plasma concentration of molnupiravir obtained from the PK model is used in the PD model to obtain time-dependent efficacy. Since molnupiravir imposes lethal mutations during the virus replication process, in our model, a portion of produced viruses, measured by ϵ(t), are mutated (Vm) and therefore assumed to be non-infectious, with the addition that most detected viral RNA pre-treatment is also non-infectious. The production rate of non-mutated viruses is decreased by a factor of (1-ε (t)). Equation 1e is written as, Total viral load (v+ vm) is used to fit the trial data.
Fitting the combined model to individual viral load data in the PLATCOV trial
We used the population mixed-effect approach and Monolix to estimate each individual’s viral dynamics parameters and the potency adjustment factor (paf). Due to the lack of data from the initial phase of infection in the PLATCOV and PANORAMIC trials, we include the data from Omicron-infected individuals in the NBA cohort to inform the model about the initial phase of infection. We fixed the PK parameters to the estimated population values (Table S1), and the PD parameters other than EC50 to the in vitro estimated population values (Table S2). We used the study category (NBA vs PLATCOV and PANORAMIC) as a covariate for tO (timing of infection) and τ (timing of the adaptive immune response) since the first recorded positive test is likely much later for the clinical trials. In the NBA study, samples were collected almost daily regardless of symptoms often leading to pre-symptomatic detection, while in the PLATCOV study, the baseline measurement occurred after symptom onset, trial enrollment and consent.
Construction of a virtual cohort
To generate a cohort for our simulated clinical trials, we randomly selected 400 individuals (for each arm of the simulated trials) from the unvaccinated symptomatic subpopulation of the NBA cohort for MOVe-OUT and the vaccinated, Omicron infected subpopulation for PLATCOV and PANORAMIC and used their individual viral load parameters estimated by fitting our viral dynamics model to the data. A sample size of n=400 (out of 822 vaccinated individuals with Omicron infection) was used to mimic a large-scale clinical trial and maintain relatively low overlap between virtual cohorts used in each arm of the simulations and between different simulations. Since the time of symptom onset is not available for all individuals in the NBA data, we randomly drew an incubation period for each individual from gamma distributions with variant-specific parameters estimated by Gamiche et al.39 The start of treatment relative to symptom onset was randomly selected according to a uniform distribution for MOVe-OUT and
PLATCOV, and a logit normal distribution for PANORAMIC with limits of [0,5] days and mean and standard deviation reported in the PANORAMIC trial for control and treatment arms. The same population PK parameters were assigned to each individual. The relevant dose in each scenario was added to the AGI compartment (the absorption equation) of the PK model (eq 2a) at each dosing timepoint (t=0, 0.5, 1, 1.5, …., 4.5 days). For each dose, the appropriate PK parameter were used (Table S1). PD parameters were also randomly drawn from a log-normal distribution with the estimated mean and standard deviation. The standard deviation of the PD parameters represents the accuracy of the assays and not individual variability.
Potency adjustment factor (paf)
The paf is defined as, We estimated the paf by maximizing R2 when fitting the change in viral load of the treatment arm of our simulation to the data from day 0 to day 7 of the treatment arm of the clinical trial.
Measuring rebound probability
A viral load rebound in the treatment arm was defined when the viral load at any time after treatment ended exceeded the viral load at the end of the treatment by 1 log. In the control group, viral rebound was defined in patients who had at least two peaks with minimum height of 1000 copies/ml in their viral load trajectories and the second peak was 1 log higher than its preceding local minimum.
Software
All the model fittings in this study were performed using Monolix version 2023R1. The data analysis and simulations were performed using Python 3.9.12.
Data Availability
The data analyzed in this work was previously published by Hay et al. and Schilling et al. and is available on GitHub at https://github.com/gradlab/SC2-kinetics-immune-history and https://github.com/jwatowatson/PLATCOV-Molnupiravir/tree/V1.0 Pharmacodynamics data is available on GitHub at https://github.com/sEsmaeili/MolnupiravirModeling
https://github.com/gradlab/SC2-kinetics-immune-history
https://github.com/jwatowatson/PLATCOV-Molnupiravir/tree/V1.0
Data availability
The data analyzed in this work was previously published by Hay et al. and Schilling et al. and is available on GitHub at https://github.com/gradlab/SC2-kinetics-immune-history and https://github.com/jwatowatson/PLATCOV-Molnupiravir/tree/V1.0
Pharmacodynamics data is available on GitHub at https://github.com/sEsmaeili/MolnupiravirModeling
Code availability
All codes and materials used in the analysis are available on GitHub https://github.com/sEsmaeili/MolnupiravirModeling
Funding
This work is supported by NIAID R01AI77512-01 (JTS, SP).