In Silico Evaluation of Paxlovid’s Pharmacometrics for SARS-CoV-2: A Multiscale Approach

Paxlovid is a promising, orally bioavailable novel drug for SARS-CoV-2 with excellent safety profiles. Our main goal here is to explore the pharmacometric features of this new antiviral. To provide a detailed assessment of Paxlovid, we propose a hybrid multiscale mathematical approach. We demonstrate that the results of the present in silico evaluation match the clinical expectations remarkably well: on the one hand, our computations successfully replicate the outcome of an actual in vitro experiment; on the other hand, we verify both the sufficiency and the necessity of Paxlovid’s two main components (nirmatrelvir and ritonavir) for a simplified in vivo case. Moreover, in the simulated context of our computational framework, we visualize the importance of early interventions and identify the time window where a unit-length delay causes the highest level of tissue damage. Finally, the results’ sensitivity to the diffusion coefficient of the virus is explored in detail.


Introduction
Even with extensive vaccination, COVID-19 will most probably not be eradicated from the human populations; thus, new options for therapy need to be explored. This paper concentrates on the mathematical evaluation, assessment, and computation-based simulation of a promising antiviral drug, Paxlovid [1], which is essentially nirmatrelvir co-packaged with ritonavir. Nirmatrelvir is a protease inhibitor that is active against M pro : inhibition of the SARS-CoV-2 main protease renders it incapable of processing polyprotein precursors, preventing virus production. Ritonavir is given as a pharmacokinetic enhancer: it slows down nirmatrelvir's metabolism, allowing a twice daily administration regimen. Paxlovid (also known as PF-07321332) is orally bioavailable, and it has excellent in vivo safety profiles [2][3][4].
Modeling has enhanced our understanding of the dynamics of viral spread, and it played an instrumental role in developing successful therapies for chronic viral infections such as HIV and HCV [5]. Mathematical models have proved to be indispensable tools in overcoming the challenges posed by the SARS-CoV-2 pandemic, too. In terms of investigating cellular-level antiviral dynamics, one of the most modern, state-of-the-art approaches consists of considering each component's physical dimension and capturing them within the framework of a spatial multiscale model accordingly. These systems incorporate size in a particular manner: rather than operating with a simple numerical value, what changes in these models-depending on the size of a given biological entity-are the mathematical tools themselves that are applied to grasp these variables on different scales.
More formally, building on our previous work [6], we define a hybrid mathematical model by merging (i) a partial differential equation representing local virus concentration, (ii) an agent-based model describing target cells in the lung and their three possible states (uninfected, infected, and dead), and (iii) a partial differential equation representing nirmatrelvir concentration. Naturally, the respective parts are closely and meaningfully intertwined: each considerable interaction and feedback process connecting these separate biological participants is given formal definition and appears in the implementation. A detailed motivation, definition and construction of this model type and its various advantages are discussed in [6]-here, we limit ourselves to highlighting the role of including crucial spatial mechanisms: unlike other classical models such as the ODE (ordinary differential equations) approach, the present system is defined in both space and time, and consequently, it is able to capture highly significant-and in nature inherently spatial-physical phenomena such as virus diffusion.
Our goal here is to provide an assessment of Paxlovid based on mathematical simulations. We explore questions such as what happens if nirmatrelvir is taken without ritonavir, or how do expectations for treatment outcome change if tablets are taken with some delay. We emphasise that performing the analogous in vivo experiments on an actual person or animal would be either simply impossible or unethical; at the same time, the corresponding in silico experiment can be conducted in just a few minutes at an ideally low cost.
While various other software packages implement hybrid mathematical concepts [7,8], we highlight that our implementation of the proposed multiscale system modeling Paxlovid is based on an adaptation of the free and open source library, HAL (Hybrid Automata Library) [9].

The Hybrid PDE-ABM Model
As described in Section 1, the main multiscale framework is defined via forming meaningful bridges between two important and fundamentally different modeling techniques: continuous partial differential equations and discrete agent-based models.
In the following, let Ω denote the area we are considering. For an in vivo experiment, this would mean a small part of the lung tissue, while in the case of an in vitro experiment it would be the area of a single well in a laboratory plate. Section 2.1.2 introduces the discrete part of our hybrid model accounting for epithelial cells, whereas Sections 2.1.2 and 2.1.3 discuss the continuous component modeling the virus (V) and drug (N) concentration. The cornerstone observation motivating this separation is that viruses and drug molecules are several magnitudes smaller than epithelial cells [10]. Then, we present the parametrization of the model in Section 2.2. Finally, Section 2.3 contains details of the implementation.

Epithelial Cells
One of the most important modeling decisions in [6] was defining epithelial cells as discrete agents in an agent-based model (ABM). The diameter of epithelial cells is relatively significant [11,12] and consequently, in terms of mathematical conceptualization, it is natural to approach cells as separate entities and follow their respective states on an individual level.
The discrete state space of target cells is defined precisely as in [6]. For completeness, we recall some of the fundamental technical details; namely, we construct a twodimensional ABM state space by introducing a lattice of k 1 × k 2 agents representing epithelial cells (k 1 , k 2 ∈ N). Cells are identified by means of the corresponding agent's place in the grid, or formally, by the Finally, by setting the Ω i,j notation for the open set occupied by the (i, j)-th cell, we havē Regarding cell states in the context of the ABM space, the main concept is rather straightforward: each agent has three potential states. The latter is formally captured by the state function s i,j (t) representing that in this modeling framework, an epithelial lung cell is considered to be either uninfected, infected, or dead: T, if the (i, j)-th cell is alive and uninfected at time t I, if the (i, j)-th cell is infected at time t D, if the (i, j)-th cell is dead at time t.
We note that the uninfected, susceptible, and target (cell) expressions are used interchangeably in the context of viral dynamics: they all refer to living cells that are susceptible to SARS-CoV-2 infection but are (for the time being) free from it.
Concerning state dynamics, the transition rules are set to mimic the biological phenomenon in question, the complete list is as follows: • All living uninfected cells are susceptible target cells to virus infection; • Since the time frame of infection is relatively short, cell birth and cell division are ignored; • Infection is not reversible: an infected cell can not become a healthily functioning uninfected cell again; • Viral infection itself is the only reason for cell death, i.e., death related to any other natural cause is not accounted for (considering the 17-month half-life of lung epithelial cells [13], natural apoptosis may be ignored over the course of a 5-day Paxlovid treatment); • The uninfected→infected state change: a target cell may become infected depending on the local virus concentration at the given cell. Infection itself is randomized and it occurs with a probability of P I (for more details, see [6]); • The infected→dead state change: Analogously to infection, death is governed by a stochastic model with probability P D .

Virus Concentration
As discussed, virus concentration V(t, x, y) is described as a variable that is continuously changing in both space and time, and as such, it is formally described by means of a partial differential equation (PDE): where D V stands for the virus diffusion coefficient, µ V represents the viral clearance rate, N is the local concentration of nirmatrelvir (i.e., the active antiviral component of Paxlovid), η N is the efficacy function of nirmatrelvir, while g i,j denotes the viral source term for the (i, j)-th infected cell. Equation (1) accounts for the following modeling assumptions.
(i) Virus particles spread across the domain primarily via diffusion. (ii) A non-specific, non-adaptive, simplified immune system is assumed which clears virions with a constant rate. (iii) A local nirmatrelvir concentration of N(t, x, y) reduces virus production from infected cells by a ratio of η N (N(t, x, y)); for more details, see Section 2.1.3. (iv) Infected cells generate new virus particles in a process that is formally described by the g i,j source functions: In general, any reasonable f i,j (t, x, y) function may be allowed in the above formula (for more details, see [6,14]). We adopted the standard simplification commonly used in the field of viral dynamics: analogously to [15,16], a constant virus budding rate is assumed. In particular, we used the estimate f i,j = 3.72 × 10 −3 copies/(mL · minute · cell); see [17].

Drug Concentration
Paxlovid (also known as PF-07321332) is an an orally administered SARS-CoV-2 main protease inhibitor [2][3][4]. It is essentially a combination of two different drugs: nirmatrelvircapable of effectively blocking virus production in infected cells-acts as its main antiviral component, while ritonavir serves to slow down the metabolism of nirmatrelvir to maintain significantly higher concentrations of the participant responsible for M pro -inhibition. We emphasize that nirmatrelvir and ritonavir are not only separate entities as acting components: the corresponding drugs themselves are packaged in individual, separate tabletsthis means that the theoretical possibility to take, for example, nirmatrelvir only (without the beneficial effect of ritonavir) is readily available. The remaining part of the section is dedicated to formulate the above statements in the context of the mathematical framework.
Nirmatrelvir concentration is modeled as a continuous variable and is denoted by N(t, x, y). We highlight that-both for simplicity and because of the apparent lack of clinical data-we do not explicitly introduce the analogous R(t, x, y) function for ritonavir concentration. Instead, we focus only on two specific scenarios: ritonavir is either taken as instructed (i.e., 100 mg of ritonavir every 12 h), or it is not taken at all. Formally, we introduce the boolean r to mathematically grasp the above concept: r = true, if ritonavir is taken following official regimen, false, if ritonavir is not administered at all.
Hence, this boolean is responsible for controlling N through the metabolismrelated descriptors.
Two further anatomical details need to be taken into account before we can formulate the governing equations for drug concentration.
First, capillary density is very high in the lung. For example, in the case of rats, there are about 11 epithelial cells per a single alveolus [18]. Considering that there are approximately 40 capillary loops per alveolus [19], this gives circa 4 capillary loops per epithelial cell. Consequently, because of the abundant presence of neighboring capillaries for a single cell, it is natural to assume a completely homogeneous drug distribution in the alveolar epithelium. Therefore, we work with N(t) instead of N(t, x, y), and the equation describing nirmatrelvir concentration becomes an ODE instead of a PDE.
Second, in order to reproduce the characteristic local concentration curves observed in clinical data (we refer to Figure 2A in [20]), we apply a standard pharmacokinetic (PK) dual compartment approach-akin to that utilized in [21]-assessing antiviral therapy targeting SARS-CoV-2. The latter model consists of a central compartment (e.g., stomach) responsible for first-level drug metabolism and a peripheral one (for the purpose of this manuscript the lung) containing the target site of nirmatrelvir. Technically, we introduce an additional c(t) function representing drug concentration at the central compartment-this is the the amount of nirmatrelvir that is already present in the patient's system but is not yet locally available at the level of the lung's epithelial cells. In the context of the two-compartment model, N(t) corresponds to the peripheral compartment's nirmatrelvir concentration. Similarly to the case of N(t), an ODE is used to describe c(t).
The complete system for nirmatrelvir concentration is then formally described by the following set of equations: where S represents the nirmatrelvir source function in the body, corresponding to a twicedaily administration regimen (the time unit being τ = 1 min): The choice of K(r) is discussed in the following section. We note that the above set of equations holds primarily for in vivo scenarios. In case of in vitro experiments, one might, for example, consider a simpler, constant presence of nirmatrelvir.

Parametrization
The configuration of the stochastic ABM state space and the PDE layer describing SARS-CoV-2 infection had been given in [6]. Here, we set the parameter values that are related to the (new) calibrated layer representing Paxlovid-based antiviral therapy.
• Drug removal rates: µ c (r), µ N (r). As we do not have direct information on the µ c , µ N coefficients, we deduce them indirectly by using frequently measured nirmatrelvir blood concentration values communicated in [20]. Of course, this argument raises the question of whether it is reasonable to use blood concentration values to estimate local drug concentrations in the lung-the validity of this approach is reassured by the results of [22]. The µ c (r) and µ N (r) coefficients were set using Wolfram Mathematica [23]. We note that in the process, we also benefited from a priori information on ritonavir from [20,24]: we used that in non-ritonavir-boosted cases, the active component is apparently metabolized 3-4 times faster. The Wolfram Mathematica notebook is available in our public Github repository [25]. • Efficacy: η N . In our model, η N is defined by means of a Hill function-the parameters of the latter are set precisely to obtain an efficacy of 50% when the drug concentration takes the value of EC 50 for nirmatrelvir w.r.t. SARS-CoV-2 (the latter parameter is approximately 62 nM according to [3]). Formally, η N is defined as For simplicity, we use the notation N(t) for drug concentration whether it is understood in nanomolars or in nanogramms per milliliter. Our implementation internally takes care of conversions when necessary due to data arriving from different sources.
The most important parameter values are summarized in Table 1. Table 1. Parameter configuration is primarily based upon best fit to actual data communicated in [20]. Previously existing parameters are defined in [6].

Implementation
The present work is a direct continuation of [6], and hence, its technical foundations and principles remain unchanged. For the sake of compactness, we avoid repetitive detailshere, we limit ourselves to summarizing the extended structure of our updated software with the help of the flowchart in Figure 1. The chart is divided into two main columns. On the left, we have the program flow itself from start to finish, while the three additional boxes on the right contain further details the respective functions. After setup and initialization, the execution is controlled by a time loop: nirmatrelvir concentration, virus concentration, and cell states are updated for each time step (in our case τ = 1 min). 1.
Nirmatrelvir concentration is calculated first (highlighted with a green frame): the most fundamental details of the process are given on the right-hand side, highlighted with an identical color. There is a clear distinction between in vitro and in vivo cases: while the first scenario is implemented assuming a constant drug concentration, the latter utilizes a two-compartment PK system.

2.
Updating virus concentration values entails recalculating the current values considering both natural decay and inflow from infected cells taking virus diffusion into account.

3.
After the new drug and virus concentration values have been obtained at each cell, we are ready to update the cell states, i.e., consider potential cell infection and cell death-the former is highlighted in purple, the latter in gray. The schematic details are shown on the right with corresponding colors; we emphasize that the respective stochastic cores of these processes are very similar to each other. In both cases, the algorithm calculates the probability of either infection or death.
Our numerical simulations are based on a free and open source java software package, HAL (Hybrid Automata Library) [9]; our source code is publicly accessible in the Github repository [25].

Replication of In Vitro Pharmacometrics of Paxlovid
The initial step in identifying and testing clinically promising antiviral drugs consists in performing a great number of in vitro experiments evaluating their overall effects. We begin with this straightforward approach, too. In this first scenario, we simulate a series of experiments corresponding to in vitro cases with different nirmatrelvir concentrations-all these configurations are otherwise identical in every other aspect. We simulate the course of SARS-CoV-2 infection over the course of four days, and we compare our computergenerated predictions with real-life observations obtained by scientific experiments assessing nirmatrelvir. Specifically, we consider Figure 3D in [2]-here, the authors evaluate PF-07321332 inhibition for (among other viruses) SARS-CoV-2 in viral-induced CPE assays, and their results are given for a series of different drug concentration values. Figure 2 demonstrates a notable resemblance to Figure 3D in [2]. Key features of inhibition efficacy match in a reassuring way: the characteristic shape itself of the calculated curve looks identical to its clinical counterpart, and the numbers connected to the main concentration window (approximately between 10 and 300 nM) corresponding to tangible increase are essentially the same, too.   Figure 3D in [2].
While there is a clear match between clinical data and our calculated results, we highlight that there is a natural limit to accuracy due to a simple lack of data. Both simulated and real-life outcomes naturally depend on key features such as the number of days the experiment went on for or the complete resolution of the state space (i.e., the total number of cells). While the supplementary material of [2] suggests that the authors mostly considered time intervals corresponding to 3-5 days, several parameter values are either unknown by nature or have not been disclosed.

Exploring In Vivo Pharmacometrics of Paxlovid
In this section, we present and explain our most significant computational results representing simplified in vivo cases. We explore a series of scenarios with various configurations, one basic feature remains unchanged however in all of them: we assume that the most fundamental instructions given in Paxlovid's documentation [3] are followed at least for its nirmatrelvir component. Technically, this means that once patients start taking Paxlovid, they steadily take at least nirmatrelvir for 5 days straight, 1 dose every 12 h. Since Paxlovid is essentially nirmatrelvir co-packaged with ritonavir, technically, it is possible that a patient-either consciously because of an existing drug allergy or simply because of forgetfulness-takes only nirmatrelvir without the added benefits of ritonavir. This degree of freedom is allowed and investigated throughout the simulations. Some other combinations and scenarios were excluded due to lack of data; again, others were omitted simply because of the limited scope of the article.
Remark 1. The inherent, rather sharp distinction between the in vitro and in vivo clinical categories becomes notably smoother in the simulated context of our mathematical model. Figuratively speaking, we perform in vivo experiments "as if they were" in vitro in the sense that we have full control over (and full information on) which specific biological or anatomical processes are included and which ones are left out. Our framework is called hybrid because of the different mathematical theories it unites, but it proves to be hybrid in this point of view as well.

In Silico Testing of Immediate Paxlovid-Based Intervention
We begin by simulating three basic scenarios and observing the respective outcomes. Figure 3 shows the course of SARS-CoV-2 infection assuming no antiviral intervention, Figure 4 follows a case where the patient takes nirmatrelvir only (i.e., the main acting component of Paxlovid, without the benefits of ritonavir), while Figure 5 represents the scenario where Paxlovid is taken exactly according to official instructions.  For the latter two cases, we plot total virus concentration and nirmatrelvir concentrations both at the first level of metabolism in the body and locally at the epithelial lung cells; the results are shown in Figure 6. In Figure 6, both the integrated virus and drug concentration values are noteworthy. Firstly, we highlight that nirmatrelvir concentration levels clearly correspond to Figure 2A in [20]-this means that our simulations (both in the ritonavir-boosted and in the nirmatrelvir-only case) are running with highly realistic nirmatrelvir concentration levels. Secondly, our computational results correspond to straightforward, basic expectations suggested by the packaging of Paxlovid. In more detail, on the one hand nirmatrelvir in itself seems to be insufficient to control the infection (which explains why Paxlovid does not simply consist of nirmatrelvir tablets), and on the other hand, ritonavir-boosted nirmatrelvir is apparently capable of stopping infection entirely (which is in accordance with the simple fact that Paxlovid is an authorized drug of great promise).

Evaluating the Effects of Treatment Delay
The previous section's premise was similar to a classical in vitro configuration-in this original default case, infection and treatment began simultaneously. In order to make our model more realistic, here, we introduce and explore a new degree of freedom: treatment delay.
We begin by exploring how simulated predictions seen in Figures 4-6 would change if Paxlovid tablets were given with a delay.
In particular, Figures 7 and 8 illustrate virus dynamical processes that are otherwise identical to the scenarios of Figures 4 and 5, respectively, except for a 36-h delay in initiating nirmatrelvir-based treatment (this also means that we follow these cases for an overall longer time period). The ritonavir-boosted scenario is particularly interesting. Although the first 36 h see uninhibited virus spread, Figure 8 confirms that Paxlovid can control infection relatively well even in this particular, less favorable scenario: after the first 2 days, there are almost no new cell infections at all, and the only detectable change between the last four subfigures is infected cells gradually turning dead.
Similarly to Figure 6, Figure 9 shows integrated virus concentration and nirmatrelvir concentration levels, but, naturally, considering a 36-h delay before Paxlovid is given.   Now, we are ready to move on to this section's main purpose, namely, investigating outcomes and eventual averted tissue damage rates for a series of delay values with respect to the default (i.e., no delay) case. Note that with no particular immune response, total tissue damage (i.e., the ratio of cells that are either infected or already dead) after 5 days reaches 100%-this means that the ratio of eventually remaining susceptible target cells after Paxlovid treatment corresponds precisely to the damage that is averted because of Paxlovid. Figure 10 illustrates the damaging effect of treatment delay from two different viewpoints. The first one, Figure 10a, considers averted damage for a series of scenarios where each scenario assumes a 12-h additional delay compared to the previous one. We highlight the sharp fall in effectiveness after a delay of 1.5 days: lack of timely Paxlovid-based antiviral intervention proves to be the most costly at this exact time window. For clarity, we note that the expression surviving cells refers to the fraction of initially susceptible cells that has not become infected by the end of the observation period.
The main idea of the second approach (Figure 10b) is to redefine the quantity on the horizontal axis: here, instead of linearly increasing delay times, the x axis follows initial damage rates (i.e., the level of damage that has been done until the moment treatment with Paxlovid is started). In other words, the latter approach depicts the relation between initial damage and averted damage. Surviving cells [%] (b) Figure 10. The damaging effect of treatment delay in two different approaches. Both subfigures illustrate the ratio of remaining uninfected target cells-the substantial difference between the two plots is the quantity measured on the horizontal axes. Subfigure (a) follows time, directly, on its x axes, while graph (b) depicts results w.r.t. initial damage rates. Results were calculated with the same fixed diffusion coefficient as used in [6], namely, D V = 0.2σ 2 /min. Due to the lack of precise clinical data-and consequently, the relative uncertaintyregarding the exact diffusion coefficient value of SARS-CoV-2, we explore the respective sensitivity of the results shown in Figure 10b. Specifically, Figure 11 illustrates the corresponding results in a heatmap for different diffusion values. Compared to the default scenario assuming D V = 0.2, the outcomes do not change substantially for even significantly higher D V values; however, there is a clear pattern suggesting that infection outcomes are expected to be more favorable if the diffusion coefficient is several magnitudes lower. Finally, we visualize the potency of Paxlovid in Figure 12: in this graph, we principally approach damage rates as areas. While this image is similar to Figure 10b in that the horizontal axis corresponds to initial damage, Figure 12 is ultimately structured differently. It distinguishes three types of damages and represents them as two-dimensional volumesnamely, we consider initial damage, damage after treatment initialization, and averted damage.
Naturally, the area between the x = y line-depicted in (dotted) red-and the horizontal axis corresponds to the level of tissue damage suffered until the moment of Paxlovidbased intervention, i.e., initial damage.
As our next step, we visualize the unavoidable damage that occurs after intervention begins: the (dotted) curve depicted with blue shows the further damage that takes place even after the patient starts taking Paxlovid. Evidently, the area between the red line and the blue curve is the visual representation of damage after treatment initialization. This rate of damage is especially high when soaring virus concentration values are combined with a significant fraction of susceptible target cells at the initialization time of Paxlovid treatment. The latter is explained simply by nirmatrelvir's mechanism of action: nirmatrelvir does effectively block virus production in infected cells, but it can not prevent target cells from getting infected, which is also apparent in the figure itself.
The third category, averted damage emerges in Figure 12 as the area between the blue curve and the horizontal line framing the graph from above (the latter naturally corresponds to the scenario where no medical intervention happens and full-scale damage takes place after 5 days). This shaded, light green area is the visual equivalent of the damage that is averted as a result of Paxlovid treatment, or in other words, the epithelial lung cells that are saved by this new M pro inhibitor. Similarly to numerous other antiviral drugs (targeting a large variety of viruses), the principle of 'the sooner the better' proves to hold in this case, too: if intervention happens right at the beginning, almost the entire cell population can be saved by Paxlovid in case of a SARS-CoV-2 infection.  The quantity on the horizontal axis (and the x = y line itself) represents the level of cell culture damage suffered until Paxlovid treatment begins, while data points depicted in blue show the unavoidable further damage that occurs after therapy commences. The shaded areas are a precise visual representation of initial damage (red), unavoidable post-intervention damage (blue), and averted damage (green). Evidently, the light green area represents those healthily functioning epithelial lung cells that were ultimately saved by Paxlovid.

Discussion
Even with worldwide vaccination programs, SARS-CoV-2 and its newly emerging variants represent an unprecedented global challenge. Consequently, new alternative treatment options are still very much needed. This paper yields a mathematical, computation-based evaluation of one of the most promising SARS-CoV-2 inhibitors to date, Paxlovid. We implemented and carefully calibrated a multiscale mathematical framework to serve as a small in silico laboratory where the basic features of Paxlovid can be replicated, explained, and further investigated. Our calculations correspond to clinical expectations remarkably well: we successully replicated the outcome of a real-life in vitro experiment in the simulated context of our model; moreover, both the sufficiency and the necessity of Paxlovid's two main components were verified by our computations for a simplified in vivo case. To further improve Paxlovid's assessment, we generated a heatmap investigating the results' sensitivity to the inherently vaguely specified virus diffusion coefficient.
The proposed hybrid model has its limitations. In its present form, the system operates on a simplistic two-dimensional grid, ignoring the complex 3D geometry of the lungs, which may introduce bias or delays in the predictions. The implementation of a biologically more realistic three-dimensional structure falls beyond the scope of this study and is subject of future research. The increase of dimension (and of lattice size) inevitably affects the computational load and, hence, requires additional, technical optimization of the code in order to achieve the desired performance. Similarly, the current assumption of a constant virus clearance rate is ignoring the intricacies of the immune system that is a major limitation when considering in vivo scenarios. While in real life, there is a significant virus release at burst, our model is averaging out the virus source over a time interval; hence, we work with constant virus production rates, similarly to [15,16]. Consequently, this limitation is responsible for a slight overestimation regarding advancement of virus release. We note that in the context of our model, it is straightforward to implement more sophisticated approaches as well (which would not be the case for example in an ODE-based system), and we also highlight that our source code includes a built-in option allowing to consider latency periods. However, that would require sufficiently detailed biological data for parametrization. Finally, we mention that our basic modeling approach to natural cellular life cycles, though being a simplification, is not expected to imply significant deviation from reality in the context of a 5-day long Paxlovid treatment: natural cell death is responsible only for circa 1% of cell population loss during this time frame according to the 17-month half-life of epithelial lung cells [13].
Despite the mathematical model's necessary simplifications and the short scope of this case study, we were able to visualize and verify the importance of early interventions; moreover, we identified the specific time window where delaying treatment initiation proves to be most costly. As for broader directions of future work, we highlight that such hybrid models and computational frameworks hold a great deal of promise with applications such as supporting clinical trials by means of in silico experiments. Computation-based evaluation and simulation of therapies not only can enhance optimization of treatments, but a further development of this technology could also serve to reduce the need for animal testing in the future.