A Novel Viral Infection Model to Guide Optimal Mpox Treatment

In 2022, worldwide mpox outbreaks have called attention to mpox virus infection and treatment opportunities using the drugs cidofovir and tecovirimat, which target different stages of in-host viral proliferation, respectively production and shedding. We propose a new model of in-host viral infection dynamics that distinguishes between the two stages, so as to explore the distinct effects of the two drugs, and we analyse the model properties and behaviour. Reducing the model order via timescale separation is shown to lead to the classical target-cell limited model, with a lumped viral proliferation rate depending on both production and shedding. We explicitly introduce the effect of the two drugs and we exemplify how to formulate and solve an optimal control problem that leverages the model dynamics to schedule optimal combined treatments.

by inhibiting envelope protein p37; although non-enveloped mature internal viruses can also leave the host cell at cell death, enveloped virus is thought to be the major contributor to cell-to-cell transmission [6], [7], [8].
Mathematical models have been successfully proposed to describe the in-host evolution of infectious diseases and plan the optimal treatment of viral infections, including HIV and COVID-19 [9], [10], [11], [12], [13], [14], [15]. The classic target-cell limited model [16], [17] of viral infection describes the interplay among three populations (uninfected target cells; infected productive cells; virus) and lumps the production of new virus within infected cells and its subsequent shedding into a single process. However, this prevents the exploration of the distinct effects of cidofovir and tecovirimat on mpox infection.
The main contribution of this letter is to propose and analyse a new model structure that explicitly shows the effect of reducing the viral production rate p and shedding rate η.
• Our new model (Section II) considers five populations: uninfected target cells, C U ; latent cells, C L , that are infected but not yet producing virus; infected productive cells, C P ; intracellular virus V I ; enveloped extracellular virus, V E , able to infect cells. We show that reducing the model via timescale separation techniques yields a target-cell limited model with a lumped constant for extracellular virus production that depends on both p and η; yet, we consider the full model to capture the specific effect of each drug and exemplify how therapies can be planned accordingly.
• We thoroughly analyse the model behaviour, characterise the disease-free equilibrium set and assess its stability, also in relation to the reproduction number (Section III).
• Considering mpox dynamics when cidofovir and tecovirimat can be infused, and possible parameter values from the literature, we formulate and solve an optimal control problem (Section IV), to demonstrate how our model could be used to plan antiviral treatments that optimally combine both drugs. Virus reproduction cycle corresponding to system (1). Uninfected cells C U are infected by external virus V E with rate constant β. The resulting latent cells C L are transformed with rate constant γ into virus producing cells C P , which produce internal virus with rate constant p and die with rate constant δ. Internal virus V I is shed with rate constant η, thus producing external virus, which is cleared with rate constant c. Internal virus is also lost, due to the death of productive cells, at rate f δ .

II. MPOX VIRUS INFECTION MODEL
where all the parameters are positive: β (t N) −1 represents the infection coefficient; γ t −1 is the rate constant of progression of an infected cell from latent to productive; δ t −1 is the death rate constant of productive cells; p t −1 is the viral production rate constant; η t −1 is the viral shedding rate constant; c t −1 is the clearance rate constant of extracellular virus. Vital dynamics of uninfected cells, death of latent cells and clearance dynamics of intracellular virus are assumed to be negligible over the considered time scale. The dynamics of V I involve production, shedding, and viral loss occurring due to productive cell death events at a rate f δ (C P ) ≥ 0.
Virus lost in such events is not enveloped, thus it does not add to the external virus population [8]. The reaction of the immune system is implicitly accounted for by the parameters: the immune response reduces β and increases δ and c. The two drugs cidofovir and tecovirimat would reduce p and η, respectively. To find an expression for the reduction of intracellular virus caused by productive cell death, we consider the virus population inside a single productive cell, V I,C (t). Experiments on isolated host cells have shown that, after the latent period, internal virus particles are produced at a constant rate p [18] and are then shed from the cell through the process of extracellular enveloped virus formation, which requires membrane protein p37 [8], [19]. The shedding process is modelled as proportional to the amount of internal virus, with rate constant η. Therefore, the internal virus population evolves as: The average lifespan of a productive cell is δ −1 . Given the solution V I,C (t) = p η 1 − e −ηt to (2), the average number of intracellular virus copies lost due to one cell death is V I,C δ −1 = p η 1 − e − η δ , corresponding to the intracellular virus population in the average productive cell at the end of its life. Multiplying this term by the death rate δC P for all productive cells yields f δ (C P ) = p η 1 − e − η δ δC P . Hence, By defining we can rewrite (3) so that the full model, with state The state variables, associated with population densities, cannot become negative, or arbitrarily large.
Proposition 1 (Positivity and Boundedness): Proof: All the system parameters are positive, including p r . In fact, since η δ > 0, also ψ η δ is positive: it is monotonically increasing, with lim η δ →0 ψ η δ = 0 and lim η δ →∞ ψ η δ = 1. Then, the state variables cannot become negative, because for all i, when Variables C U , C L and C P are individually upper bounded by κ, because they cannot be negative. Then,

Boundedness of the trajectories is ensured.
We can reduce the model by exploiting two time-scale separation arguments. If cells start producing virus very soon after infection, the dynamics of the latent cell population C L can be neglected and γ C L = βC U V E . Also, if the intracellular virus population is assumed to be in a quasi-steady state, then V I = p η ψ η δ C P can be plugged into the equation ofV E . The resulting reduced-order model with state has the same form as the classical target-cell limited model, but the production rate p r of extracellular virus depends on both the shedding rate η and the production rate p, as in (4).
To be able to capture the specificity of each drug action at best, henceforth we focus on the complete model (5).

III. ANALYSIS OF THE IN-HOST INFECTION DYNAMICS
We provide a qualitative analysis of the dynamics of the complete system (5). Similar results are known to hold for the target-cell limited model, as proven, e.g., in [9].

A. Equilibria, Reproduction Number, Critical Cell Number
As shown in Proposition 1, system (5) is positive, since R 5 ≥0 is an invariant set for the system, and its trajectories are ultimately bounded in the compact and convex set X ∞ = {x ∈ R 5 : 0 ≤ x ≤ x max }. This ensures the existence of a steady state in X ∞ [20], [21]. Only two types of steady states are possible: the trivial equilibrium, where all state variables are zero; and the disease-free equilibrium, where only C U has a strictly positive value, while all other variables are zero. The equilibrium set is therefore We consider the set of meaningful initial conditions If the infection occurs at time t 0 , it leads to a discontinuity in V E (t): when t < t 0 , all variables are zero apart from the strictly positive C U (healthy state); for t = t 0 , all state variables are zero apart from the strictly positive C U and V E . The infection progression can be related to the in-host basic reproduction number R 0 (defined as the number of host cells infected by a single infected cell when a small initial viral load targets a population of C U,0 > 0 uninfected cells), which can be computed as the spectral radius of the next-generation matrix [22], [23]. For system (5), we have The critical number C * U of uninfected cells, defined as the value of C U,0 that solves (9) for R 0 = 1, is To perform a local stability analysis, we evaluate the system Jacobian matrix at the equilibriumx = [C U 0 0 0 0]: The block-triangular matrix Jx is singular and its spectrum is completed by the eigenvalues of A, whose characteristic polynomial has all positive coefficients apart from the constant term a 0 = γ η(δc − p r βC U ). The Routh-Hurwitz table reveals that the disease-free equilibrium is unstable, due to the presence of a positive real eigenvalue, if a 0 < 0, i.e., R 0 > 1 or, equivalently,C U > C * U . WhenR 0 < 1, or equiva-lentlyC U < C * U , matrix A is Hurwitz, but the zero eigenvalue prevents us from drawing conclusions on the stability of the equilibria based on the system linearisation.

B. Asymptotic Behaviour
Here, we show that the equilibrium subset is the smallest set that is attractive in X (the distance from X * eq of all trajectories originated in X converges to zero as t → ∞) and is Lyapunov stable (for all > 0, there exists ζ > 0 such that, if the distance of x(t 0 ) from X * eq is smaller than ζ , then the distance of x(t) from X * eq is smaller than for all t ≥ t 0 ). Hence, X * eq is the smallest asymptotically stable equilibrium set, with domain of attraction X .
All trajectories starting in X converge to a state in X eq . Proposition 2 (Asymptotic Behaviour): For system (5), if x(t 0 ) ∈ X at some time t 0 , then lim t→∞ x(t) = x eq ∈ X eq , with lim t→∞ C U (t) = C ∞ ∈ [0, C U (t 0 )).
Proof: BeingĊ U (t) < 0 ∀t ≥ t 0 when x(t 0 ) ∈ X , C U is monotonically decreasing, hence lim t→∞ C U (t) = C ∞ is a constant value in the interval [0, C U (t 0 )). As a consequence, lim t→∞ĊU (t) = 0 and, from (5a), We can show that X * eq is an attractive set, since C ∞ ∈ [0, C * U ], and is actually the smallest. Theorem 1 (Attractivity of X * eq ): For system (5), X * eq in (12) is the smallest attractive equilibrium set in X . Proof: To prove attractivity, we show that C ∞ ∈ [0, C * U ]. Summing the first three equations of system (5) yieldsĊ U + C L +Ċ P = −δC P . Rearranging (5d) gives Dividing both sides of (5a) by C U yieldṡ Integrating and taking the limit for t → ∞ of both sides, since where R 0 , defined as in (9), and Q 0 = p r β δc (C L,0 + C P,0 ) + Moreover, X * eq is the smallest attractive set in X : from two different initial conditions, the system converges to two different equilibria in X * eq , hence neither single states in X * eq nor subsets of X * eq are attractive in X ; for more details, see the proof of [9, Th. 3.1] for a model of the form (6).
We can also prove that X * eq is Lyapunov stable. Theorem 2 (Lyapunov Stability of X * eq ): For system (5), the equilibrium set X * eq in (12) is Lyapunov stable. Proof: Given the generic equilibriumx = [C U 0 0 0 0] with C U ∈ [0, C * U ), we consider the candidate Lyapunov function which is continuously differentiable, zero at the equilibrium and positive for all other x ∈ R 5 ≥0 . The Lyapunov derivative wheneverC U ≤ C * U . Hence, V is a Lyapunov function for system (5), proving stability of all equilibria in X * eq . Fig. 2. Asymptotic value C ∞ of C U as a function of the initial condition C U,0 , when β = e and all other parameters are set to 1, resulting in C * U = 1. All the statements in Proposition 3 hold for V E,0 = 10 −5 . As the initial amount of external virus V E,0 is increased, C ∞ (C * U )< C * U , but all statements except for (ii) still hold.
In view of Theorems 1 and 2, X * eq is asymptotically stable. Theorem 3 (Asymptotic Stability of X * eq ): For system (5), X * eq in (12) is the smallest asymptotically stable equilibrium set, with domain of attraction X .
The domain of attraction cannot be the whole nonnegative orthant: all states with C U = 0 or V E = 0 are equilibria when C L = C P = V I = 0.
The asymptotic number of uninfected cells enjoys properties akin to those proven in [9] for a system of the form (6). We state and prove the properties, visualised in Fig. 2.
Proposition 3 (Asymptotic Value of C U ): For system (5) starting from the initial condition [C U,0 0 0 0 V E,0 ] ∈ X , with V E,0 > 0 small enough, denoting by C ∞ (C U,0 ) the asymptotic value of C U with initial condition C U,0 , we have: Proof: From the proof of Theorem 1, we know that C ∞ = −C * U W 0 (−R 0 e −(R 0 +Q 0 ) ). Given an initial condition with C L,0 = C P,0 = V I,0 = 0 and V E,0 ≈ 0, it is Q 0 ≈ 0 and thus C ∞ ≈ −C * U W 0 (−R 0 e −R 0 ). Then: −1 , 0), the composed function is increasing: Remark 1 (Transient Behaviour): The evolution of (5) is different from that of (6), analysed, e.g., in [9], for which, if R 0 ≤ 1,V E (t) ≤ 0 for all t > t 0 when only C U (t 0 ) and V E (t 0 ) are non-zero. For some parameter choices, model (5) can yield an increase in V E even when R 0 < 1. Two typical behaviours of system (5), shown in Fig. 3, correspond to monotonic healing (V E (t) ≤ 0) when R 0 < R < 1 is small enough, and acute infection spreading otherwise. The different transient evolution suggests the importance of considering the full model to plan treatment, whenever the dynamics of C L and V I cannot be safely neglected.  of R(t). Left: R 0 <R < 1; V E decreases monotonically. Right: R 0 > 1; V E first decreases, reaches a minimum, then grows up to a maximum (and keeps growing whileR < R(t)< 1) before finally decreasing to zero.

IV. OPTIMAL TREATMENT WITH ANTIVIRAL DRUGS
We demonstrate how our novel model could be used to support optimal mpox treatment. Interest in mpox has only very recently arisen and therefore the available data are still scarce. Our model parameters are thus not identified from time series, due to the lack of specific clinical data, but inferred from the literature, also on the closely related Vaccinia virus [24]. Here, we exemplify the general methodology and the insight that our novel model could provide within an optimal control framework for therapy design, once the most appropriate parameter values for the specific clinical case have been chosen. The resulting optimal therapy depends on the parameter values, which are subject to large uncertainties and variability across patients. Even though all the model parameters are specific for each patient, given the typical infection and treatment duration of about 14 days [5], [30], the time required for data collection on the single patient would lead to late therapeutic interventions; a possible practical approach towards personalised treatment is to plan specific therapies for different patient categories, taking into account, e.g., age, sex, weight, specific risk factors.
We assume that the antiviral drugs cidofovir (inhibiting intracellular virus production, thus reducing p) and tecovirimat (inhibiting intracellular virus excretion, thus reducing η) can be infused daily (for 12 hours) through inputs u 1 and u 2 respectively. Each drug administration cannot exceed a maximum tolerated dose. The drug concentrations D i (t), i = 1, 2, in the body exhibit first-order decay, as is common in pharmacokinetic modelling: Each drug inhibits its target process through a Michaelis-Menten function of D i , with half-inhibitory concentration h i [31]. Hence, the drug-dependent rate constants are: Our goal is to minimise the viral load V E over a finite horizon, both throughout the dynamic evolution (with an integral cost) and at the final time (with a terminal cost), and possibly also V I throughout the horizon. This maximises the amount of surviving healthy cells C U . We also wish to minimise the side effects caused by the two drugs [32], [33]. With weights w 1 , w 2 , w 3 , α 1 , α 2 ≥ 0 and cost functional Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.  = 1200 [mg day −1 ], assuming constant infusion for 12 hours followed by 12 hours without infusion. The initial conditions are set as C U,0 = 5 · the optimal control problem can be formulated as: where the input is piecewise constant: u i (t) = u i,j for t ∈ [t j , t j+1 ], with j = 1, . . . , N − 1, t 1 = 0, t N = T = 14 days and t j+1 − t j = 12 hours ∀j; and it is u i,j = 0 for j even. We numerically solve the optimal control problem with a multiple shooting method using CasADi [34] through its MATLAB interface. We assume a treatment period of 14 days. The chosen parameter values are discussed in the caption of Fig. 4, which compares the effects of the optimal treatment schedules obtained with different choices of the cost weights.
The four scenarios illustrate the type of conceptual insight that we can expect from embedding our model within an optimal control approach. In the clinical practice, tecovirimat is preferred over cidofovir [5], [31], because it has a stronger efficacy for the same dose (in fact, we set h 1 h 2 ) and milder side effects (in fact, we choose the cost weights as α 2 = 0.001 < α 1 = 0.01). Even though all the advantages of tecovirimat are reflected in our chosen parameters, our numerical results show that the optimal choice is not to administer tecovirimat only (as is typically done in practice), but to use (small) doses of cidofovir in synergy with it.
To allow a complete comparison, not only between the considered treatment scenarios, but also between the treated and untreated evolution, all the bottom panels of Fig. 4 also report, with dashed lines, the dynamic evolution of the system with the same parameters, but without any treatment (u 1 = u 2 = 0). Without treatment, in 15 days C U would be reduced from its initial value of 5·10 5 to 5·10 4 (while the viral load values would have a peak around 10 7 ), leading to tissue necrosis and to the patient's death. In Fig. 4A (w 1 = 100, w 2 = 1, w 3 = 0), both controls are saturated, although cidofovir (u 1 ) is only given for the first two days, while tecovirimat (u 2 ) for the whole therapy horizon; however, even though V E is brought below one, V I decreases very slowly. When a fast reduction of the amount of intracellular virus is an important objective, as in Fig. 4B (w 1 = 100, w 2 = 1, w 3 = 10), an even stronger use of cidofovir becomes necessary. Conversely, in Fig. 4C (w 1 = 10, w 2 = 1, w 3 = 0), cidofovir (u 1 ) is only given for the first day at the maximum dose and for the second day at a small dose, while tecovirimat (u 2 ) is given for the whole therapy horizon, at the maximum dose for the first two days and at smaller doses afterwards. In Fig. 4D (w 1 = 10, w 2 = 1, w 3 = 0.1), a slightly faster reduction of V I is required (barely noticeable in logarithmic scale) and is obtained with higher doses of cidofovir and smaller doses of tecovirimat.
It is worth stressing that, to successfully reach the treatment goal, it is not necessary to drive V E and V I to zero: once the viral load goes below a safe threshold, the patient's immune system will take care of it and suppress the infection.

V. CONCLUDING DISCUSSION
We have proposed a new model of in-host virus infection dynamics that explicitly captures internal virus production and shedding, to account for the respective effects of the antiviral drugs cidofovir and tecovirimat, candidates in the treatment of mpox. The classic target-cell limited model is recovered from the proposed model as a quasi-steady-state approximation; we have shown that our model shares many features with the target-cell limited model, such as equilibria and asymptotic behaviours, but can have different transient behaviours. With model parameters taken from the medical literature, we have showcased the general methodology to formulate and solve an optimal treatment control problem, illustrating how our model can be used to support the optimal scheduling of antiviral drug therapies. Future work includes quantitative model validation and parameterisation based on infection time-series data from in vitro experiments and/or patients suffering from mpox infection, so as to inform real therapies through optimal control approaches with tailored parameter values. Since previous studies have effectively used the target-cell limited model with constant coefficients to fit real viral infection data (e.g., for HIV, West Nile virus and Influenza [35], [36], [37]), we expect that our proposed model, which is more detailed and thus more flexible, will be able to match real data even with constant coefficients.