A Computational Framework for Exploring SARS-CoV-2 Pharmacodynamic Dose and Timing Regimes

: Emerging diseases—and none as recently or devastatingly impactful toward humans as COVID-19—pose an immense challenge to researchers concerned with infectious disease. This study is tasked with expanding the computational probe of treatment regimes in a differential equations-based model of the SARS-CoV-2 host–virus interaction. Parameters within the model are tweaked to simulate dose speciﬁcations. Further, parametric variations are introduced in a timed manner to infer the importance of dose timing. Arming in silico testing, and eventually, clinical testing, with abundant information on simulated therapeutic regimes is the overall contribution of this pharmacodynamic model; thus, a wide range of dose and timing combinations are examined. Therapeutic interventions that block viral replication inhibit viral entry into host cells, and vaccination-induced antibodies are all studied alone and in combination. Especially during early detection, exhaustive parameter sweeps of well-suited within-host models are often the ﬁrst step in the clinical response to a novel disease. immune response; pharmacodynamic model MSC: 92-10; 92B05


Introduction
Model-based analyses play a very important role in informing the human struggle with infectious diseases. This has become increasingly important as the emergence of the SARS-CoV-2 virus has renewed urgency in epidemiological research. The terrible burden of COVID-19 has been mitigated expediently with the adaptation of antiviral therapies and the inclusion of novel immunization techniques. In this context, in-host modeling is one way to discuss the efficacy of immune response to therapy at the cellular level. An in-host model allows researchers to probe a patient's response to viral exposition and therapies based on timing, dosage, individual mitigating factors or patient-specific co-morbidities.
Mathematical in-host models allow researchers to model as many in-vivo metrics as the modeler prescribes. Exerting control over metrics describing rates of viral replication, antiviral agents and other key host-virus interactions allows conclusions concerning the utility of various treatment regimes to be analyzed by the pharmacodynamic model. Increasing the predictive performance of such differential equation-based models is important as clinicians face uncertainty in treating novel diseases. Within-host dynamical models are arguably the first and best weapon against novel pathogens-serving as an intermediary between therapeutic discovery and/or adaptation and clinical intervention.
Differential equation-based models have been used to estimate the time-to-infection, to differentiate secondary infections from imported strains, and to inform potential therapies [1,2]. Other models intensively represent the interactions between virus and components of the host's immune response, with a focus on T-cell interactions [3,4]. Bifurcation analyses have been employed to conclude how extremely rapid SARS-CoV-2 replication overcomes T-cell response [5]. Sanche et al. provide a framework that focuses on viral interactions leading to disease-progressing inflammation [6]. The interaction of angiotensinconverting enzyme inhibitors and angiotensin II receptor blockers with COVID-19 is also modeled thoroughly by within-host models [7]. In this paper, the effects of three COVID-19 therapies-which emerged as top prospects in combating the disease during our early understanding of it-are studied: Remdesivir [8], Host-cell entry inhibitors [9] and antibody/vaccination treatments. To simulate Remdesivir treatment, which inhibits viral transcription rate, we throttle the rate in our model, which describes the replication of infected cells. To simulate host-cell entry inhibitor antiviral therapy types, we throttle the rates that describe how the virus interacts with healthy cells to produce latently infected (termed as latent cells from hereon) and then infected cells. Finally, to model vaccination, we introduce an alternate antibody agent in the model to capture the presence of vaccineinduced antibodies. The inclusion of vaccination efficacy in the framework is crucial for interpreting patient-specific simulations. Chemaitelly and Abu-Raddad [10] study waning immunity provided by vaccination as a function of time since vaccination and many other personal patient characteristics such as: age, sex, place of birth and existing co-morbidities. Similar host-virus models have probed the effectiveness of other potential treatments, such as hydroxychloroquine [11] or favipiravir [12]. Additionally, the importance of patientspecific characteristics is highlighted by the simulation of adverse neurological reactions to certain COVID-19 treatment regimens in stroke patients [13] or as a function of a patient's geographic location [14]. Jenner et al. explore patient-specific outcomes using within-host models to describe individual patients' predisposition to produce infected cell-derived interferon [15]. Ordonez-Jimenez et al. describe patient-specific reactions as a function of initial infection severity [16]. Finally, it must be noted that pharmacodynamic models are subject to limitations imposed by modeling choices. Venisse et al. [17] discuss the need for existing host-virus models for SARS-CoV-2 to reduce the possibility of unreliable conclusions and call urgently for the development of further pharmacodynamic models. The results presented herein widen explored therapeutic parameter regimes while being still limited by computational cost. Further collaboration with clinical studies being needed to bridge the gap between purely computational probes and clinically informative numerical simulation.
This work provides a framework to probe large treatment doses and timing regimes in the context of treating emerging diseases. In the case of the novel coronavirus, key in-host viral dynamics have been studied [18]. This study sweeps various parameter combinations to provide insights into how dose, therapy combination, and vaccination status can affect an individual's response to treatment. The novelty of large parameter sweeps, including timed administration simulation, is a useful tool for computational, clinical integrative treatment design. Further, the inclusion of vaccination-induced antibodies and the vaccination efficacy parameter, η, allow for explicitly patient-specific modeling as a function of time since vaccination. Additionally, a stability analysis to conclude the neutral stability of the system is performed.

Materials and Methods
The within-host SARS-CoV-2 viral model is probed with potential therapeutic interventions. The effect of widely used Remdesivir, a monophosphoroamidate prodrug that is utilized to inhibit the viral transcription rate, is modeled by controlling the parameter associated with viral replication. The potential use of host-cell inhibitors is modeled by controlling the parameters that inhibit the virus's ability to catalyze the transformation of healthy cells to latent cells. A summary of the key viral interactions in the within-host model, first described by Sadria and Layton [18], is given in Figure 1. The entire in-host model is governed by Equations (1)- (12). The system is solved employing MATLAB function ode15s: a variable-step, variable-order multistep solver based on numerical differential formulae of order 1-5 by default. The model's components are summarized in Table 1. The parameters and their values are displayed in Table 2. Original parameter values are borrowed from precursor work modeling influenza infection in-host model also derived from clonal selection theory, law of mass-action, interaction characteristics and natural cell creation-death balances [19]. Model parameters are then tuned to fit observational COVID-19 infection data. In brief, the model describes interactions between viral infection, respiratory epithelial cells, and the host's immune system. The model allows the control of key parameters to simulate therapeutic intervention. Model parameters are scaled such that viral load clearance is realized in about three weeks, which is quite typical in COVID-19 patients.
Ultimately, a modeling paradigm to include vaccination is also introduced. The model in Figure 1 already includes a general antibody agent, A. Sadria and Layton present an Here, β A is the rate by which antibodies are introduced through plasma donor treatment. A proportional change in Equation (1) would be observed to reflect additional viral elimination due to A . The case of β A = 0 and initial value A = A 0 can be abstracted as exposure to some initial concentration of antibodies through vaccination. While in the case A 0 = 0, β A = 0 is the rate of incoming antibodies through infusion.
An additional modeling modification is introduced to examine the effect that newly developed vaccines have on SARS-CoV-2 within-host interactions. To abstract vaccine efficacy, parameter η is introduced. η is a scaling parameter that affects the rate by which vaccine-induced antibodies are effective: allowing the study of the patient outcome as a function of waning immunity. η = 0 implies no vaccination or complete loss of immunity, while η = 1 implies that vaccine-induced antibodies, X, are equally as effective as regular antibodies, A. Thus, to model patient outcome, when vaccines are introduced, Equation (1) becomes:  [20,21]. A brief description of the parameters is given.

Parameter
Description Virus reproduction rate in infected cells γ VA = 309. 6 Virus elimination rate by antibodies Rate by which virus enters healthy cells α V = 0.85 Natural virus degradation rate Interferon production rate by antigen-presenting cell δ F = 1000 Interferon production rate by infected cell β FH = 8.5 Rate by which interferon binds to healthy cells α F = 4 Interferon natural decay rate β EM = 4. 15 Rate by which effector cells are produced by antigen-presenting cells β EI = 1.36 Rate of effector cell death by infected cell Effector cell natural death rate β PM = 5.75 Plasma cell production rate α P = 0.2 Plasma cell natural death rate β A = 0.0225 Antibody production rate by plasma cells γ AV = 73. 1 Rate by which antibody binds to virus α A = 0.0215 Antibody natural death rate r = 0.000015 Rate of antibody specificity change

Linear Stability Analysis
In order to discuss the stability of the solutions under the conditions prescribed in Section 2, a linear stability analysis [22][23][24] is performed. Let the system of differential Equations (1)-(11) be represented by F(x). The model is evaluated at equilibrium, where F(x eq ) = 0. Small perturbations are introduced around such equilibrium and variables are replaced where x = x eq + ∆x. Thus, d∆x dt is approximated by F(x eq ) + J∆x. This implies where J is the Jacobian matrix of F where x = x eq is given by Figure 2 Thus, the system is determined to be at least neutrally stable, with each eigenvalue containing real components taken to be equal to or less than zero. In other words, for small changes in initial parameters, one can still reasonably expect the robustness of equilibrium solutions.

Results
Viral load within the simulated host throughout infection is modeled as a function of therapy type, treatment efficacy and treatment timing. Treatment efficacy is described on a [1, 0] scale-vaccine efficacy describes a boost in antibody effectiveness, while simulated therapies represent a decrease in viral replication rate, Remdesivir or host-cell entry inhibition with a decrease in viral infectiousness rates, γ HV and γ V H . For instance, in the modeling of Remdesivir and entry-inhibition therapy, efficacy = 0.9 means the associated rate constant has been reduced to 90% of its normal value; thus, γ HV is 90% of its original value, resulting in a lower cell infection rate by the virus. An efficacy of 1 implies no administered treatment, while lower values of efficacy correspond to the treatment regime with either lower cell infection rates or lower rates of entering the healthy cells by the virus. In other words, here, efficacy denotes the reduction in the associated rate constants and does not directly equate to the overall efficiency of the treatment. The simulations are 25 days in duration. The viral load on day 1 is set to 30,000 parts per mL, simulating a relatively early detection of the virus. Efficacy is most readily abstracted as treatment dosage, and the two terms are used interchangeably henceforth. However, efficacy quantification implies a more general approach to therapeutic modeling, which could include an individual patient's predisposition to certain treatment types or to model other factors that may affect in-host interactions.

Singular Therapies
First, the dependence of maximum viral load on treatment efficacy is examined. Remdesivir treatment efficacy of 0.9-decreasing the viral transcription rate by 90%corresponds to the efficacy of about 5 µM Remdesivir [8]. Efficacy is varied around this value from 0.75 to 0.99 in steps of 0.01.
Host-cell entry inhibitor therapy is modeled for a range of inhibition efficacy. This therapy is modeled by toggling rate constants, which govern the effect the virus's presence has on producing latent cells. Practically, many potential therapies that inhibit viral entry into host cells could be modeled by the wide range of efficacy values chosen as 0.75 or 0.90 and values between 0.50 and 0.80 in steps of 0.01 to produce contours.
Finally, the effects of vaccination effectiveness are also discussed through the modified system. Chosen values are in the range between 0 and 1 in steps of 0.01. Figure 3 clearly shows the patient's simulated outcome is dependent upon the paradigm of modeled treatment efficacy. In both the case of Remdesivir and host-cell entry inhibition, increased efficacy results in decreased maximum viral load. Additionally, the time until clearance metric, the time passed to clear two-thirds of the initial viral load, is also dependent on Remdesivir dosage.   Most interestingly, the effect of vaccine antibodies seems to have a cut-off, by which the introduced viral load is cleared immediately, whereas less capable antibodies result in nearly full-term infection duration (as shown in Figure 5). Further study into this phase transition will help uncover the understanding of how waning immunity affects in-host interactions. The rate at which vaccine effectiveness is lost and the efficacy of vaccination against SAR-CoV-2 variants are of paramount concern as we continue developing novel treatments. Comparing the exploratory study of parameter choice to emerging models of waning vaccine immunity [25] can help inform therapy regimens as dependent upon vaccination status or time since the last booster.

Combination Therapy
Here, multiple treatments are applied contemporaneously for the entire treatment window (Figures 6 and 7), or in time-dependent dosing schedules (Figure 8). The results of combination therapy tests are of paramount importance to this study-the specific case study of COVID-19 allows modelers to inform clinical studies by describing the interaction of various dose and treatment timings to provide individual treatment suggestions. This is especially important when existing drugs are being repurposed for emerging diseases and when new treatments are being developed; thus, the important interactions between novel treatment techniques beg large-scale investigation. The proposed generic treatment efficacy model suggests a course of action for clinicians.  Figure 3a, we can see how the introduction of host-cell therapy greatly affects a patient's viral load dependence on Remdesivir dosage. The shift from a linearly shaped curve in Figure 3 to the S-shaped curve in Figure 6 suggests interesting phase transitions when studying concurrent treatment regimens. This is an important conclusion in that it suggests to clinicians that treatment plans do not respond linearly when varying dosage and timing-again begging further study before administering treatment.  Figure 7 displays the varying dependence of viral load on both treatments. One important conclusion suggests that the effect Remdesivir has on patient maximum viral load is dependent on the efficacy of concurrent host-cell entry inhibitor treatments. The nonlinearity further suggests the importance of broad-ranging parametric sweeps.
Additionally, timing and treatment efficacy are varied between dosing regimens in Figure 8. The relationship between timing a 5µM Remdesivir treatment with a hostcell entry treatment is displayed. Both time to viral peak and maximum viral load are dependent heavily on dose timing. Where the dose is constant, timing is still shown to produce important results in viral load and clearance rate. Next, the combination of vaccination and both proposed therapeutics is examined. This highlights again the effect that individual patient characteristics have on treatment outcome. In this case, the patient-specific information is vaccination status. However, broadly speaking, this model could be extrapolated to include other patient-specific predispositions to treatment.
Many important conclusions can be drawn from Figures 9 and 10. First, the introduction of treatment greatly affects the range of efficacy by which vaccines produce near-instantaneous clearance. The apparently step-wise phase transition dependency on treatment dosage offers insight into the level of treatment needed for a patient depending on time since vaccination.  Finally, the results of various timed treatment combinations in key vaccination efficacy regimes are explored to highlight the differing outcomes. In other words, depending on patient predisposition, the "best" treatment-in terms of maximum viral load-is not always the same treatment timing and combination choice. This is most apparent in Figure 11, where identical treatment regimes are tested for various vaccination efficacies, and the mixed results for producing maximum viral load are reported.  Figure 11 further examines how optimal treatment timing depends upon vaccine efficacy. The effects on viral load in the first five days of treatment are observed in Figure 12. The host's viral load is shown to depend closely upon dose timing. Figure 11, in conjunction with the maximum viral loads reported in Figure 11 and time until clearance reported in Figures 9 and 10, illuminates the diversity of outcomes attainable by small differentiations in dose and timing. The diversity of outcomes further stresses the utility of pharmacodynamic models in integrative clinical treatment regimens.
(a) (b) Figure 12. Viral load as a function of time for various treatment timing configurations of host-cell entry inhibitor followed by Remdesivir with patient vaccine efficacies of (a) 0.1 and (b) 0.15. Colored lines describethe time before the entry-inhibitor to Remdesivir therapy swap. Blue, black, orange, green and red correspond to one, two, three, four and five days, respectively.

Conclusions
The in-host model of SARS-CoV-2 infection is shown to be sensitive to parameter conditions describing the dose, timing and efficacy of various COVID-19 therapies. This correlation is drastically important to understanding how COVID-19 patients should be treated. COVID-19, a novel disease with rapidly evolving treatments, requires exploratory studies of treatment regimens.
The results presented include modeling the effect of Remdesivir or host-cell entry inhibition on a patient's viral load. The clear dependence on small increments in treatment efficacy underscores the importance of in-host models in informing clinical explorations.
Further investigation of how different vaccinations can interact with different therapies needs only minor tweaks to the modeling paradigm and available data on vaccines and outcomes. Furthermore, the parameters in Table 2 remain widely unchanged as they represent the average patient. However, systemic parametric changes could be made to describe underlying differences in patient populations. Pre-existing conditions and co-morbidities could be added to the model through informed parameters. Importantly, in-host COVID-19 therapy models must continue to adapt to challenges presented by treatment design and patient condition modeling.
Funding: This work is partially supported by NSF grant CBET-1802588. Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.