key: cord-0117579-sairo7d0 authors: Perez, Mara; Abuin, Pablo; Actis, Marcelo; Ferramosca, Antonio; Hernandez-Vargas, Esteban A.; Gonzalez, Alejandro H. title: Optimal control strategies to tailor antivirals for acute infectious diseases in the host date: 2021-06-17 journal: nan DOI: nan sha: 508d21edefb31ddf621d0cbb1a5afe9369f3c67e doc_id: 117579 cord_uid: sairo7d0 Several mathematical models in SARS-CoV-2 have shown how target-cell model can help to understand the spread of the virus in the host and how potential candidates of antiviral treatments can help to control the virus. Concepts as equilibrium and stability show to be crucial to qualitatively determine the best alternatives to schedule drugs, according to effectivity in inhibiting the virus infection and replication rates. Important biological events such as rebounds of the infections (when antivirals are incorrectly interrupted) can also be explained by means of a dynamic study of the target-cell model. In this work, a full characterization of the dynamical behavior of the target-cell models under control actions is made and, based on this characterization, the optimal fixed-dose antiviral schedule that produces the smallest amount of dead cells (without viral load rebounds) is computed. Several simulation results - performed by considering real patient data - show the potential benefits of both, the model characterization and the control strategy. Mathematical models of with-in infections can be used to characterize pathogen dynamics, optimize drug delivery, uncover biological parameters (including pathogen and infected cell half-lives), design clinical trials, among others. They have been employed to study chronic (i.e.: HIV [1, 2, 3] , hepatitis B [4, 5] , hepatitis C [6, 7] ) and acute (i.e.: influenza [8, 9, 10] , dengue [11, 12] , Ebola [13] ) infections. Currently, they are based on ordinary differential equations (ODE), which allows to analyze these systems employing mathematical and computational tools. This way, in-host basic reproduction numbers (R), stability analysis of equilibrium states, analytical/numerical solutions, can be computed [14, 15, 16, 17] . Most of them are based on the target-cell limited model to represent chronic/acute infections according to the infection resolution respect to the target cell production and natural death rates [18] . This way, the equilibrium states differ from isolated equilibrium points (i.e.: disease free and infected equilibria) for the former to a continuous of equilibrium points (i.e.: disease free equilibrium set) for the latter ones. Note that for acute infections, the only feasible equilibria is the disease free, since the pathogen particles at the end of infection will be cleared independently of the in-host reproduction number [18, 8, 19] . The existence of healthy equilibrium set implies that stability analysis can be performed considering equilibrium sets as a generalization of equilibrium points, which gives an environment to employ set-theoretic methods [20, 21] , widely used in the design of set-based controllers, although not fully employed for modeling characterization and control of acute infections. Some preliminary results, which will be discuss later in this chapter, can be found in [22] . The control of infection can be modelled considering immune response mechanisms, where the infection is self-controlled by a combination of a non-specific and specific reactions [23, 24, 25] , or by drug therapies. The inclusion of pharmacokinetic (PK) and pharamacodynamic (PD) models of drug therapies allows the inclusion of therapeutic effects on the pathogen evolution [7, 18] . Therefore, the models parameters can be changed exogenously by dose frequency and quantity, naturally limited by the inhibitory potential of the drug (expressed in terms of EC50, or drug concentration for inhibiting 50% of antigen particles) and its cytotoxic effect (expressed in terms of IC50, or drug concentration which causes death to 50% of susceptible cells) [26] . Moreover, since drugs are normally administrated by pills or intravenous injections, instantaneous jumps are observed in the concentration of the drug in some tissues. This is mathematically conceptualized as a discontinuity of the first kind and gives rise to the so-called impulsive control systems [27] . This model representation has been used for optimal control and state-feedback control with constraints for infectious diseases, such as: influenza [10, 28] and HIV [29, 27] . Even though optimal dosage can be computed for chronic and acute models, the unstable healthy equilibrium of the former (under certain circumstances; for details, see [18, 15] ) and the availability of target cells above a critical level for the latter (as it is discuss later), involve the duration of drug therapy, with the presence of viral rebounds when therapy is disrupted. This effect has been noticed for chronic [30] and acute [31] infections. Taking into account this scenario, in this work, we formalize the existence of an optimal single interval drug delivery such that viral rebounds are avoided. Even though, the presented analysis is valid for the target-cell limited model for acute infections, taking into account the current worldwide contextual situation (COVID-19 pandemic), we prove our results using an identified model of infected patients with SARS-CoV-2 virus [32, 22, 33] . After the introduction given in Section 1 the article is organized as follows. Section 2 presents the general "in the host" models used to represent infectious diseases. Section 4 studies the way the antivirals affect the dynamic of the model, emphasizing the fact that the stability analysis made in Section 3 remains unmodified and, so, any control strategy must be designed accounting for these details. In Section 5 control design able to exploit the stability model characterization is introduced, and its benefits are shown by simulating several cases, in Section 6. Finally, conclusions are given in Section 7. First let us introduce some basic notation. We consider R n as n-dimensional Euclidean space equipped with the euclidean distance between two points defined by d(x, y) := x − y = [(x − y) ′ (x − y)] 1/2 . The euclidean distance from a point x to a set Y is given by d(x, Y) := x Y = inf{y ∈ Y : x − y }. With X we will denote the constraint set of R 3 , given by We will consider X endowed with the inherit topology of R 3 , i.e. the open sets are intersections of open set of R 3 with X . Thus, a open ball in X with center in x and radius ε > 0 is given by B ε (x) := {y ∈ X : x − y < ε} and an ε-neighborhood of set Y ⊂ X is given by The interior of Y is the set of all interior points of Y and it is denoted by int(Y). Mathematical models of in-host virus dynamic have shown to be useful to understand of the interactions that govern infections and, more important, to allows external intervention to moderate their effects [24] . According to recent research in the area [18, 22] , the following ordinary differential equations (ODEs) are used in this work to describe the interaction between uninfected target or susceptible cells U [cell/mm 3 ], infected cells I [cell/mm 3 ], and virus V [copies/mL]: is positive, which means that U (t) ≥ 0, I(t) ≥ 0 and V (t) ≥ 0, for all t ≥ 0. We denote x(t) := (U (t), I(t), V (t)) the state vector, and X = R 3 ≥0 the state constraints set. The initial conditions of (2.1), which represent a healthy steady state before the infection, are assumed to be V (t) = 0, I(t) = 0, and U (t) = U 0 > 0, for t < 0. Then, at time t = 0, a small quantity of virions enters the host body and, so, a discontinuity occurs in V (t). Indeed, V (t) jumps from 0 to a small positive value V 0 at t = 0 (formally, V (t) has a discontinuity of the first kind at t 0 , i.e., lim t→0 − V (t) = 0 while Although the solution of (2.1) for t ≥ t 0 , being t 0 an arbitrary time, is unknown, we know that it depends on the basic reproduction number 1 R = βp cδ and the initial conditions (U (t 0 ), is a non increasing function of t (by 2.1.a). From [22] and [17] , if c >> δ (as it is always the case 2 ) it is known that for U (t 0 )R ≤ 1, then V (t) is a non increasing function of t, for all t ≥ t 0 , and goes asymptotically to zero for t → ∞. On the other hands, if U (t 0 )R > 1, V (t) reaches a maximumV and then goes asymptotically to zero, for t → ∞. In this latter case it is said that the virus spreads in the host, since 1 The reproduction number is usually defined as R = U (t 0 ) βp cδ > 0, for UIV-type models. However, for the sake of convenience, we remove the initial value of U in our definition. 2 If c >> δ. system (2.1) can be approximated byU Then, since U (t 0 ) > 0, conditions for V to increase or decrease at t 0 , are given by U (t 0 )R > 1 and U (t 0 )R < 1, respectively. there is at least one time instant for whichV > 0 [22] . The so called critical value of U , U * , is defined as where R is assumed to remain constant for all t ≥ t 0 . The critical value U * can be seen as the counterpart of the "herd immunity" in the epidemiological SIR-type models: i.e., U (t) reaches U * approximately at the same time as V (t) and I(t) reach their peaks or, in other words, V (t) and I(t) cannot increase anymore once U (t) is below U * . This way, conditions U (t 0 )R > 1 and U (t 0 )R ≤ 1 that determines if V (t) increases or decreases for t ≥ t 0 can be rewritten as U (t 0 ) > U * and U (t 0 ) ≤ U * , respectively. In what follows, we assume that U (0) > U * (or U (0)R > 1), which corresponds to the case of the outbreak of the infection (i.e., the virus does spread in the host), at time t = 0. Let us now define U ∞ := lim t→∞ U (t), V ∞ := lim t→∞ V (t) and I ∞ := lim t→∞ I(t), which are values that depend on R and the initial conditions U (t 0 ), V (t 0 ), I(t 0 ). According to [22] , V ∞ = I ∞ = 0, while U ∞ is a value in (0, U (t 0 )), which will be characterized in the next section. To find the equilibrium set of model (2.1), with initial conditions (U (t 0 ), V (t 0 ), I(t 0 )) ∈ X at an arbitrary time t 0 ≥ 0,U (t),İ(t) andV (t) need to be equaled to zero, in (2.1). According to [22, 17] there is only one equilibrium set in X , which is a healthy one, and it is defined by To examine the stability of the equilibrium points in X s , a first attempt consists in linearizing system (2.1) at some state x s := (U s , I s , V s ) ∈ X s , and analyzing the eigenvalues of the Jacobian matrix. As it is shown in [22] , this matrix has one eigenvalue at zero (λ 1 = 0), one always negative (λ 2 < 0) and a third one, λ 3 , that is negative, zero or positive depending on if U s is smaller, equal of greater than U * , respectively. Since the maximum eigenvalue λ 3 is the one determining the stability of the system, it is possible to separate set X s into two subsets, according to its behaviour. Then, a first intuition is that the equilibrium subset is stable, and that the equilibrium subset is unstable. However, this is not a conclusive analysis, given that one of the eigenvalues of the linearized system is null and so the linear approximation cannot be used to fully determine the stability of a nonlinear system (Theorem of Hartman-Grobman [34] ). Formal asymptotic stability of set X st s , together with its corresponding domain of attraction, is analyzed in the next subsection. A key point to properly analyze the asymptotic stability (AS) of system (2.1) is to consider the stability of the equilibrium sets X st s and X un s , instead of the single points inside them (as defined in Definitions 7.2, 7.3 and 7.5, in Appendix 1). Indeed, even when every equilibrium point in X st s is ǫ − δ stable, there is no single equilibrium point in such set that is locally attractive. As stated in Definition 7.5, in Appendix 1, the AS of X st s requires both, attractivity and ǫ − δ stability, which are stated in the next two subsections, respectively. Finally, in Subsection 3.4 the AS theorem is formally stated. According to Definition 7.2, any set containing an attractive set is also attractive. So, we are in fact interested in finding the smallest closed attractive set in X \X un s . Proof. The proof is divided into two parts. First it is proved that X st s is an attractive set, and then, that it is the smallest one. Attractivity of X st s : To prove the attractivity of X st s in X (and to show that X un s is not attractive) we needs to prove that U ∞ ∈ [0, U * ] for any initial conditions and values of R. U ∞ can be expressed as a function of R and initial conditions, as follows where W (·) is (the principal branch of) the Lambert function and (U (t 0 ), I(t 0 ), V (t 0 )) are arbitrary initial conditions at a given time t 0 ≥ 0. The minimum of U ∞ is given by U ∞ = 0, and it is reached when U (t 0 ) = 0 (for any value of R, I(t 0 ) and V (t 0 )). The maximum of U ∞ , on the other hand, is given by U ∞ = U * , and it is reached only when U (t 0 ) = U * and I(t 0 ) = V (t 0 ) = 0 (for any value of R), as it is shown in Lemma 7.7, in Appendix 2. Then, for any (U (t 0 ), I(t 0 ), V (t 0 )) ∈ X and R > 0, U ∞ (R, U (t 0 ), I(t 0 ), V (t 0 )) ∈ [0, U * ], which means that X st s is attractive, and the proof of attractivity is complete. Figure 1 shows how U ∞ behaves as function of U (t 0 ) and V (t 0 ), when I(t 0 ) = 0 and I(t 0 ) = 5e 5 . The first one is the scenario corresponding to t = 0, when a certain amount of virus enters the healthy host. X st s is the smallest attractive set: It is clear from the previous analysis, that any initial state x(t 0 ) = (U (t 0 ), This means that X un s is not attractive for any point in X \X s . However, to show that X st s is the smallest attractive set, we need to prove that every point x s ∈ X st s is necessary for the attractiveness. Let us consider a initial state of the form (U * , I(t 0 ), 0) with I(t 0 ) ≥ 0. Since W is a bijective function from (−1/e, 0) to (−1, 0), then U ∞ (R, U * , ·, 0) is bijective from (0, +∞) to (0, U * ). Hence for every point x s ∈ int(X st s ) there exists I(t 0 ) ≥ 0 such that the initial state (U * , I(t 0 ), 0) converges to x s . Since every interior point of X st s is necessary for the attractiveness, then the smallest closed attractive set is X st s , and the proof is concluded. (a) U∞ as function of U(t0) and V (t0), when I(t0) = 0. The orange plane represents U * = 1/R = 1.5e 8 . The maximum of U∞ is reached when U(t0) = U * and V (t0) = 0, and is given by U * . Patient 'A'. (b) U∞ as function of U(t0) and V (t0), when I(t0) = 5e 5 . The orange plane represents U * = 1/R = 1.5e 8 . The maximum of U∞ is reached when U(t0) = U * and V (t0) = 0, and is smaller than U * . Patient 'A'. The next theorem shows the formal Lyapunov (or ǫ − δ) stability of the equilibrium set X st s . Proof. We proceed by analysing the stability of single equilibrium pointsx : For eachx let us consider the following Lyapunov function candidate This function is continuous in X , is positive for all nonegative x =x and J(x) = 0. Furthermore, for x(t) ∈ X and t ≥ 0 we havė whereẋ(t) represents system (2.1). FunctionJ(x(t)) depends on x(t) only through V (t). So, independently of the value of the parameterŪ ,J(x(t)) = 0 for V (t) ≡ 0. This means that for any single x(0) ∈ X s , V (0) = I(0) = 0 and so, V (t) = 0, for all t ≥ 0. SoJ(x(t)) is null for any x(0) ∈ X s (i.e, it is not only null for x(0) =x but for any x(0) ∈ X s ). On the other hand, for x(0) / ∈ X s , functionJ(x(t)) is negative, zero or positive, depending on if the parameterŪ is smaller, equal or greater than U * = δc βp , respectively, and this holds for all x(0) ∈ X and t ≥ 0. So, for anyx ∈ X st s ,J(x(t)) ≤ 0 (particularly, forx = (Ū , 0, 0) = (U * , 0, 0),J(x(t)) = 0, for all x(0) ∈ X and t ≥ 0) which means that eachx ∈ X st s is locally ǫ−δ stable (see Theorem 7.6 in Appendix 1). Finally, whenŪ = 0, i.e.x = (0, 0, 0), we define the Lyapunov functional as J(x) = U − I + δ/pV and we proceed analogously as before to prove the local ǫ − δ stability of the origin. Therefore, since every state in X st s is locally ǫ − δ stable and X st s is compact, by Lemma 7.4, the whole set X st s is locally ǫ − δ stable. Finally, since X st s is attractive in X \X un s then it is impossible for any x ∈ X un s to be ǫ − δ stable, which implies that X st s is also the largest locally ǫ − δ stable set in X s , which completes the proof. . This means that it is not true thatJ(x(t)) < 0 for every x =x, and this is the reason why we cannot use the last part of Theorem 7.6 to ensure the asymptotic stability of particular equilibrium points (or subsets of X st s ). In fact, they are ǫ − δ stable, but not attractive. In the next Theorem, based on the previous results concerning the attractivity and ǫ − δ stability of X st s , the asymptotic stability is formally stated. Proof. The proof follows from Theorems 3.1, which states that X st s is the smallest attractive in X , and 3.2, which states that X st s is the largest locally ǫ − δ stable set in X . Figures 2 shows phase portrait plots of system (2.1), corresponding to different initial conditions. In this section some characteristics of system (2.1) concerning the value of U ∞ as a function of the reproduction number R and the initial conditions are analyzed. Consider the next Property. This means that the closer U (t 0 ) is to U * from above, the closer will be U ∞ to U * from below. iii. For U (t 0 ) < U * and fixed and ) reaches its maximum over X , and the maximum value is given by U * (see Lemma 7.7 in Appendix 2). The proof of the properties are omitted for brevity. However, Figures 1 and 11 show how U ∞ behaves for different values of initial conditions. All along this work we use a virtual patient, denoted as patient 'A', to demonstrate the results of each section. The parameters of patient 'A' were estimated by using viral load data of a RT-PCR COVID-19 positive patient -reported in [32] and used in [22, 33] -and are given by The initial conditions are given by: U 0 = 4 × 10 8 , I 0 = 0 and V 0 = 0.31. Furthermore, the reproduction number is R = 1.84 × 10 −8 , while the critical value for the susceptible cells is U * = 5.44 × 10 7 . The final value of U (if no antiviral treatment is applied) is given by U ∞ = 2.57 × 10 5 , which means that the area under the curve (AUC) of V is given by AU C V = 5.45 × 10 7 . The peak of V is given bŷ V = 1.98 × 10 7 . Figure 3 shows the time response corresponding to patient 'A'. As predicted, U ∞ is (significantly) smaller than U * , which means that antivirals reducing (even for a finite period of time) either p or β will increase U ∞ and, so, will reduce the AUC and, probably, the peak of V . Remark 3.6. Note that the area under the curve of V , between times t d1 and t d2 is given by . This way, if U ∞ is increased with respect to the value corresponding to the untreated case, the AUC of viral load decreases. Moreover, as it was shown in [17] , the viral load at time to peak is monotonically decreasing with antiviral therapy reducing β or p. The idea now is to formally incorporate the pharmacodynamic (PD) and pharmacokinetics (PK) of antivirals into system 2.1, to obtain a controlled system, i.e. a system with certain control actions -given by the antivirals -that allows us to (even partially) modify the whole system dynamic according to some control objectives. In contrast to vaccines that kill the virus, antiviral just inhibits the virus infection and replication rates, so reducing the advance of the infections in the respiratory tract. The PD is introduced in system (2.1) as follows: (4.1a) where η(t) ∈ [0, 1) represents the inhibition antiviral effects affecting the infection rate β (note that, according to [17] , the effect of antivirals on the replication rate p, is anal-ogous to the one on β, since both parameters affect in the same way the reproduction number R). On the other hand, the PK is modeled as a one compartment with an impulsive input action (to properly account for pills intakes or injections): where D is the amount of drug available (with D(0) = D 0 = 0), δ D is the drug elimination rate and the antiviral dose u k enters the system impulsively at times t k := kT , with T > 0 being a fix time interval and k ∈ I. Time t − k denotes the time just before t k , i.e., is a continuous-time system impulsively controlled, which shows discontinuities of the first kind (jumps) at times t k and free responses in t ∈ [t k , t k+1 ) (see [27] for details). Finally, the way the drug D enters system (4.1) is by means of η as follows: where EC 50 represents the drug concentration in the blood where the drug is halfmaximal. η(t) is assumed to be in [0, η max ), with η max < 1 (not full antiviral effect is considered, since this is an unrealistic scenario). Based on the PK and PD previous analysis, the complete Covid-19 infection model, taking into account an antiviral treatment (the controlled system) reads as follows: (4.4a) with initial conditions given by x 0 = (U 0 , I 0 , V 0 , D 0 ). Given that D(t) ≥ 0 for all t ≥ 0, the constraint set X is enlarged to beX := R 4 ≥0 . Also, a constraint for the input, u, is defined asŨ := {u ∈ R : 0 ≤ u ≤ u max }, where u max represent the maximal antiviral dosage (u max is usually determined by the drug side effects and maximal effectivity, 0 ≤ η max < 1), while sets X is enlarged by consideringX := X × R ≥0 . A detailed study of the stability of impulsive systems can be seen in [35] . We resume here the simulation of the virtual patient 'A', to demonstrate the impulsive control actions describing the effects of antiviral administration. It is assumed that antivirals affect the infection rate β, while the initial condition for D is D 0 = 0, δ D = 2 (days −1 ) and EC p 50 = 75 (mg). A scenario of 30 days was simulated, and a permanent dose of u k = 20 (mg) of antivirals is administered each T days, starting at t i = 4 days, with T = 1, T = 2 and T = 0.5. As shown in Figures 4a -4b , the system response is quite different for different sampling times. For T = 1 days, the antiviral treatment is able to decreases V from the beginning. On the other hand, for T = 2, the treatment is unable to stop the spread of virus (V continue increasing after the treatment is initiated), as it is shown in Figure 4b . Clearly, the effect of larger values of T is equivalent to smaller values of the dose u k . In what follows, for the sake of clarity, T will be fixed in 1 day and only the (constant) value of the doses u k -together with the initial and final time of the treatment -will be modified to analyze the different outcomes. Control objectives in 'in host' infections can be defined in several ways. The peak of the virus load uses to be a critical index to minimize, since it is directly related to the severity of the infection and the ineffective capacity of the host. However, other indexes -usually put in a second place -are also important. This is the case of the time the infection lasts in the host over significant levels [17] -including virus rebounds after reaching a pseudo steady state, and the total viral load or infected cells at the end of the infection (i.e., the AUC of V and I). These latter indexes also informs (in a different manner) about the severity of the infection and the time during which the host is able to infect other individuals, and are directly determined by the amount of susceptible cells at the end of the infection. So the twofold control objective is defined as follows: As it was said in the previous Section, antivirals affect the infection rate β, by the time-variant factor (1 − η(t)). Accordingly, the reproduction number R will be also time varying, following the formula: and the original reproduction number -i.e., the one corresponding to no treatmentwill be denoted as R(0) for clarity (R(0) is the reproduction number at the outbreak of the infection, when u 0 = 0 and η = 0). We will assume in the following a single interval antiviral treatment, consisting in a single fixed dose of antiviral, applied during a finite period of time. At the outbreak of the infection (t = 0), it is (U (0), I(0), V (0)) := (U 0 , 0, ǫ), with ǫ > 0 arbitrary small. Then, the single interval treatment is defined by the following input function: where t i t i , but finite. Note that after t f , u(t k ) = 0, which means that η(t) → 0 and R → R(0). The control problem we want to solve first reads as follows: for a given initial time, t i t i (as it is always the case), the system converges to an equilibrium state (U ∞ , 0, 0) with U ∞ ≤ U * , being U * the critical value for U corresponding to no antiviral treatment, i.e., U * = 1/R(0). We proceed by contradiction. Assume that U ∞ > U * . Consider system (4.1) for t ≥ t f . Since the antiviral treatment is interrupted at time t f , then η(t) ց 0, for t ≥ t f 3 . By eq. (5.1) we have that R(t) ր R(0), for t ≥ t f . If we denote U c (t) = 1/R(t) we have U c ց U * for t ≥ t f . Since U ∞ > U * then U ∞ > U c (t) for t > T , with T large enough. Hence (U (t), I(t), V (t)) is converging to an equilibrium point in the unstable part, which is a contradiction. Therefore U ∞ ≤ U * , which concludes the proof. V tot := ∞ t=0 V (t)dt ≈ ∞ t=0 p c I(t) In the search of such a value, the next definition is stated. Definition 5.4 (Goldilocks antiviral dose). The goldilocks antiviral dose (GAD), u g = u g (t i ), is the one that, if applied at t i 0 arbitrary small, and R(0) such that U (0) > U * . Consider also single interval antiviral treatment (as the one defined in (5.2)), with a given starting time t i ∈ (0,t(R(0))), and a finite final time t f . Define soft and strong treatments depending on if u i < u g or u i > u g , respectively. Define also long and short term treatments depending on if the system reaches or does not reach a QSS at t f . Then, the following scenarios can take place: i. Quasi optimal single interval antiviral treatment: if u i = u g , and t f is such that ii. Soft long-term antiviral treatment: if (U (t f ), i.e., U (t) will remain approximately constant for t ≥ t f . Furthermore, the softer a soft long term antiviral treatment is, the smaller will be iii. Strong long-term antiviral treatment: if (U (t f ), I(t f ), V (t f )) reaches a QSS, and U (t f ) > U * , a second outbreak wave will necessarily take place at some timet > t f and, finally, the system will converge to an Furthermore, the stronger a strong long term antiviral treatment is, the larger will be the second wave and the smaller will be U ∞ (R(0), i. Given that u i = u g is implemented for t ∈ [t i , t f ], t f is finite but large enough and U ∞ (R g , U (t i ), I(t i ), V (t i )) = U * , then U (t f ) approaches U * and V (t f ) approaches zero, from above, as t f increases. This means that at t f , when the treatment is interrupted, (U (t f ), I(t f ), V (t f )) is close to the unstable equilibrium set X un s . Then, by Property 3.5.(ii), function is to the equilibrium point (U * , 0, 0), with U (t f ) > U * , the closer will be U ∞ (R(0), U (t f ), I(t f ), V (t f )) to U * , with U ∞ < U * (see the 'pine' shape of U ∞ around U * , for I ≈ 0, in Figure 11 ii. Given that (U (t f ), ) is close to the stable equilibrium set X st s , when the treatment is interrupted. Then, the system will converge to an equilibrium with . Softer antiviral treatment produces smaller values of U (t f ) and, by Property 3.5.(iii), smaller values of U (t f ) produce smaller values of U ∞ (R(0), U (t f ), I(t f ), V (t f )). 5 Indeed, by the ǫ − δ stability of the equilibrium state (U * , 0, 0), for each (arbitrary small) ǫ > 0, it there exists δ > 0, such that, if the system starts in a ball of radius δ centered at (U * , 0, 0), it will keeps indeterminately in the ball of radius ǫ centered at (U * , 0, 0). Furthermore, it is possible to define invariant sets around (U * , 0, 0) by considering the level sets of the Lyapunov function (3.5), withŪ = U * , or even the level sets of function J(U, I, V ) := U * − U∞(R, U, I, V ), with a fixed R > 0. This way, once the system enters any arbitrary small level set of the latter functions, it cannot leaves the set anymore. See, Figure 5 iii. Given that (U (t f ), I(t f ), V (t f )) approaches a steady state with U (t f ) > U * , then (U (t f ), I(t f ), V (t f )) is close to the unstable equilibrium set, X un s , when the treatment is interrupted. Then, the system will converges to an equilibrium in the stable equilibrium set, X st s , with U ∞ (R(0), U (t f ), I(t f ), V (t f )) < U * . Stronger antiviral treatment produces greater values of U (t f ) and, by Property 3.5.(ii), values of U (t f ) farther from U * , from above, produce values of U ∞ (R(0), U (t f ), I(t f ), V (t f )) farther from U * , from below. When U (t f ) is significantly greater than U * , no matter how large is t f and how small is V (t f ) 6 , the system will evolve to an equilibrium in X st s , with for t > t f , the system significantly increase V (t), and this effect is known as a second outbreak wave. iv. Given that (U (t f ), I(t f ), V (t f )) is a transitory state, then it does not approach any equilibrium. This means that V (t f ) is significantly greater than 0, and according to Lemma 7.7, in Appendix 2, the maximum of ≥0 : I ≥ ε, V ≥ ε is given by −W (−RU * e −R(U * +ε+ δ p ε) )/R, which is a decreasing function of ε, and reaches U * only when ε = 0 (see Figure 11 ). Then, independently of the value of U (t f ), , V (t f )) will be (maybe significantly) smaller than the one obtained with quasi optimal single interval treatment, in which ε ≈ 0. The quasi optimal single interval treatment clearly accounts for a steady state condition, given that any realistic treatment needs to be interrupted at a finite time. Furthermore, given that only one antiviral dose, u i , is considered for the treatment, once the quasi optimal single interval treatment is determined (u i = u g ), also is the peak (maximum over t) of the virus,V : i.e., there is a uniqueV for each single interval control action, u i . However, if a more general control action is considered, in such a way that u k assume several values in the interval from t i to t f ,V can be arbitrarily reduced. Indeed, given that U ∞ depends only on the fact that U ≈ U * and V ≈ 0 at t f , then stronger antiviral doses can be used at the beginning of the treatment to lower the peak of V . If for instance two consecutive single interval control actions are implemented -the first one with a high dose, applied from t i to t 1 , and the second one with the quasi optimal antiviral, u g (t 1 ), applied from t 1 to a large enough t f -a lower peak of V will necessarily be obtained in contrast to one corresponding to the quasi optimal single interval control. Although this chapter is not devoted to analyze control strategies different from the single interval one, it is worth to remark this latter point since it states that: (1) both control objectives are independent, in the the sense that if a given upper bound for V is stated from the the beginning (to avoid complication and/or to reduce the infectivity of the host) it is in general possible to design control strategies that both, make V (t) not to overpasses the upper bound, and make U ∞ ≈ U * , and (2) the entire concept of a maximum or peak for V , for a given treatment, has sense only when U ∞ ≤ U * ; since otherwise, a rebounds of the virus will occurs once the treatment is interrupted, and a new peak for V may be reached. Figures 10a and 10b , in the Simulation section, show an example of a two-steps interval treatment that produces a peak of V smaller than the one corresponding to the quasi optimal single interval one. In this section each of the cases of Theorem 5.8, together with the case of two-steps interval control action of Subsection 5.2 are simulated for data coming from patient 'A', introduced in section 3.6 and 4.2. As it was already said, δ D = 2 (days −1 ), EC p 50 = 75 (mg) and the sampling time is selected to be T = 1 day. Initial conditions are given by (U 0 , I 0 , V 0 ) = (4 × 10 8 , 0, 0.31). Also, recall that U * = 5.44 × 10 7 and the untreated peak of V is given byV = 1.98 × 10 7 . Figure 6a shows the time evolution of U (logarithmic scale), V , and u k for patient 'A', when strong long-term antiviral treatment is implemented. The treatment starts at t i = 4 days and finished at t f = 30 days, while several strong doses are administered: u i = [21, 25, 35] mg. As it can be seen, at t f the value of U is greater than U * while V ≈ 0, so the viral load V rebounds after some time, producing a second (and larger) peak. More At time t f , when the treatment is interrupted, V (t f ) ≈ 0 and U (t f ) > U * , so the system is close to an unstable equilibrium point. So, for t > t f the state is attracted to an equilibrium in the AS equilibrium set X st s , following outer level curves of J. Outer level curves of J means both, a small U ∞ and a largeV . Figure 7a shows the time evolution of U (logarithmic scale), V , and u k , when soft long term antiviral treatment is implemented. The treatment starts at t i = 4 days and finished at t f = 30 days, while several soft doses are administered: u i = [4, 6, 8] mg. As it can be seen, at t f the value of U is smaller than U * , while V ≈ 0, so the viral load V decreases after the treatment is interrupted. The values of U ∞ andV corresponding to the three doses are given by U ∞ = [8.98×10 6 , 1.90×10 7 , 3.25×10 7 ], andV = [1.36 × 10 7 , 1.16 × 10 7 , 9.70e × 10 6 ]. Figure 7b shows the phase portrait in the space U, V , together with the level curves of the Lyapunov function J(U, I, V ). At time t f , when the treatment is interrupted, V (t f ) ≈ 0 and U (t f ) < U * , so the system is close to a stable equilibrium point. So, for t > t f the state remains almost unmodified. Figure 8a shows the time evolution of U (logarithmic scale), V , and u k , when the quasi optimal single interval antiviral treatment is administered. The treatment starts at t i = 4 days and finished at t f = 30 days, while the Goldilocks dose is given by u i = u g (t i ) = 10.5 mg. The values of U ∞ andV are given by U ∞ = 5.34 × 10 7 and V = 7.73 × 10 6 , respectively Figure 8b shows the phase portrait in the space U, V , together with the level curves of the Lyapunov function J(U, I, V ). As it can be seen, the system follows the only one trajectory that goes directly from (U 0 , V 0 ) to (U * , 0): any other path goes necessarily to an equilibrium with U ∞ < U * . Figure 9a shows the time evolution of U (logarithmic scale), V , and u k , when a short term treatment is implemented. The treatment starts at t i = 4 days and finished at t f = 15 days, while several doses -smaller and greater than u g (t i ) are administered: u i = [10, 15, 20, 25] mg. The values of U ∞ andV corresponding to the four doses are given by U ∞ = [2.98 × 10 6 7, 1.39 × 10 7 , 5.06 × 10 7 , 2.26 × 10 7 ], andV = [8.05 × 10 6 , 4.98 × 10 6 , 6.67 × 10 6 , 1.01 × 10 1 ]. Figure 9b shows the phase portrait in the space U, V . Given that trajectories go along the level curves of the Lyapunov function J(U, I, V ), any short term treatmenti.e., producing V (t f ) ≈ 0 -will make the system to surround the state (U * , 0) by an outer level curve, thus finishing at some U ∞ significantly smaller than U * . As before, outer level curves of J means both, a small U ∞ and a largeV . 6.5. Two-steps treatment, lowering the peak of V Finally, a scenario is simulated to show that always it is possible to lower the the peak of V -while maintaining U ∞ ≈ U * -if a control sequence more complex that the single interval one is implemented. Figure 10a shows the time evolution of U (solid blue line, logarithmic scale) and V (solid red line) corresponding to a two-steps interval control: the first step consisting in u i = 25 mg, from t i = 4 to t m = 30 days, and the second one consisting in u i = u g (t m ) = 5.6 mg, from t m = 30 to t f = 60 (solid line). Also, the quasi optimal single interval control of Subsection 6.3 is shown, to compare the performance (dashed line). As it can be seen, the peak of V is significantly reduced: fromV = 7.73 × 10 6 toV = 2.57 × 10 6 , while U ∞ is almost the same in both cases. Figure 10b shows the phase portraits of the two control strategies (solid line, two-steps control; dashed line, single interval control), where it can be seen also the reduction of the virus peak. This simple two-step strategy shows that with a more sophisticated control strategy (i.e., by means of a proper optimal control formulation) the virus peak can be arbitrarily reduced, maintaining the condition U ∞ ≈ U * . This is indeed, matter of future research. In this work, the stability and general long term behavior of UIV-type models have been fully analysed. A quasi optimal control action -consisting in the finite-time single interval antiviral treatment producing the minimal possible final amount of death cells -was found. The analysis shows also that more complex control strategies can account for both control objectives simultaneously: minimize the virus peak, while keeping the final amount of death cells at its maximum. A detailed analysis of subotimal scenarios permits to enumerate the following main results: i. To apply soft antiviral treatment during a long time (even no treatment at all), expecting the non-infected cells would evolve alone to the critical value U * , is not ii. To apply strong antiviral treatment for a long time, expecting the virus will die out alone is not an option. Strong antiviral produces an values of U at the end of the treatment larger than U * , but this final values are artificially stable steady states, since once the treatment is interrupted or reduced, a virus rebound will necessarily occurs at some future time, and U ∞ will be significantly smaller than U * . iii. To apply any antiviral treatment (soft or strong) for short period of time, such that the system is not able to reach a quasi steady state (i.e., when V at the end of the treatment is not close to zero) is not an option. If the treatment is interrupted at a transient state, the initial conditions for the next time period are such that U ∞ will be significantly smaller than U * . iv. According to the latter results, the best option is to apply an antiviral treatment such that the system reaches a quasi steady state with U ≈ U * and V ≈ 0 at the end of the treatment. This is what we call "the quasi optimal single interval antiviral treatment", since it makes the system to approach the maximal final value of uninfected cells (U ∞ ≈ U * ), without infection rebounds. v. An important point to be remarked is that the quasi optimal single interval antiviral treatment does not determine the peak of the virus. Quasi optimal conditions for U ∞ are stationary, while condition for minimizingV are transitory, so both objective can be accounted for simultaneously. Future works include the study of more complex control strategies (mainly model based control strategies as MPC and similar) and the explicit consideration of timevarying immune system. Lemma 7.4. Let X s be a compact equilibrium set. If every x s ∈ X s is ǫ − δ locally stable, then X s is ǫ − δ locally stable. Proof. Given ǫ > 0, there exists δ = δ(x s ) > 0 for each x s ∈ X s such that if x ∈ B δ(xs) (x s ) then φ(t; x) ∈ B ǫ (x s ) for t ≥ 0. The family of δ-balls form a open cover of X s . Let us denote the union of this cover V , i.e. V := {B δ(xs) (x s ) : x s ∈ X s }. Since X s is compact and the complement of V is closed, then the distance between them is strictly positive, i.e. δ * := d(X s , V c ) > 0. Therefore, the δ * neighborhood of the equilibrium set X s is contained in V . Thus if x ∈ B δ * (X s ) ⊂ V then φ(t; x) ∈ B ǫ (X s ) for t ≥ 0 .Therefore X s is ǫ − δ locally stable. Definition 7.5 (Asymptotic stability (AS) of an equilibrium set). Consider system 7.1 constrained by X and a set Y ⊆ X . A closed equilibrium set X s ⊂ X is asymptotically stable (AS) in Y if it is ǫ − δ locally stable and attractive in Y. Next, the theorem of Lyapunov, which refers to single equilibrium points and provides sufficient conditions for both, local ǫ − δ stability and assymptotic stability, is introduced. Theorem 7.6. (Lyapunov's stablity theorem [36, Theorem 4.1]) Consider system 7.1 constrained by X and an equilibrium state x s ∈ X s . Let Y ⊂ X be a neighborhood of x s and consider a function V (x) : Y → R such that V (x) > 0 for x = x s , V (x s ) = 0 andV (x(t)) ≤ 0, denoted as Lyapunov function. Then, the existence of such a function in a neighborhood of x s implies that x s ∈ X s is locally ǫ − δ stable in Y. If in additioṅ V (x(t)) < 0 for all x = x s , then x s is asymptotically stable in Y. The following Lemma describe the behavior of the maximum of U ∞ on each Ω(ε). Lemma 7.7 (Maximum of the function U ∞ ). Consider the function U ∞ given by (7.2) and for each ε ≥ 0 the domains Ω(ε) given by (7.3) . Then the maximum of U ∞ (U, I, V ) in Ω(ε) is reached in (U * , ε, ε). In particular, the maximum value of U ∞ over Ω(0) is reached in (U * , 0, 0) and is given by U ∞ (U * , 0, 0) = U * , where U * = 1/R. Proof. According to (7.2) , U ∞ can be written as Through the change of variables x = RU and y = R(I + δ/pV ), f can be studied as a function of the form g(x, y) = xe −(x+y) . Note that (U, I, V ) ∈ Ω(ε) if and only if x ≥ 0 and y ≥ η where η := R(1 + δ p )ε ≥ 0. Therefore to find extremes of f in Ω(ε) it is enough to study the extreme points of g over Ω ′ = {(x, y) ∈ R 2 ≥0 : y ≥ η}. Since ∇g = [(1 − x)e −(x+y) , −xe −(x+y) ] does not vanish and g → 0 when (x, y) → ∞, then the maximum is reached at the boundaries of Ω ′ . A simple analysis shows that g restricted to the boundary of Ω ′ achieves its maximum in (1, η) . This means that f (U, I, V ) achieves its maximum in U = 1/R = U * and I = V = ε. In particular, when ε = 0, f (U, I, V ) reaches its maximum in (U * , 0, 0). Furthermore, which concludes the proof. Dynamics of HIV infection of CD4+ T cells An in vivo pharmacokinetic/pharmacodynamic model for antiretroviral combination Modeling the within-host dynamics of HIV infection Modeling the mechanisms of acute hepatitis b virus infection Hepatitis c virus kinetics Hepatitis c viral dynamics in vivo and the antiviral efficacy of interferon-α therapy Viral kinetic modeling: state of the art Kinetics of influenza A virus infection in humans Influenza A virus infection kinetics: quantitative data and models Passivity-based inverse optimal impulsive control for Influenza treatment in the host The role of antibody in enhancing dengue virus infection Modelling original antigenic sin in dengue viral infection Ebola Virus Infection Modelling and Identifiability Problems Reproduction numbers of infectious disease models Stability analysis of pathogen-immune interaction dynamics Virus dynamics: a global analysis Dynamical characterization of antiviral effects in covid-19 In-host modeling The mechanisms for within-host influenza virus control affect model-based assessment and prediction of antiviral treatment Model predictive control: theory, computation, and design Set-theoretic methods in control Gonzalez, Characterization of SARS-CoV-2 Dynamics in the Host Mathematical models for immunology: current state of the art and future research directions Modeling and Control of Infectious Diseases in the Host: With MATLAB and R Assessing mathematical models of influenza infections using features of the immune response Assessing bioavailablility of drug delivery systems: mathematical modeling Control strategies for nonzero set-point regulation of linear impulsive systems Inverse Optimal Impulsive Control Based Treatment of Influenza Infection Impulsive control of single-input nonlinear systems with application to hiv dynamics Modeling hepatitis C virus dynamics: liver regeneration and critical drug efficacy Neuraminidase inhibitors for treatment of human and avian strain influenza: A comparative modeling study Virological assessment of hospitalized patients with COVID-2019 -host Modelling of COVID-19 Kinetics in Humans, medRxiv Differential equations and dynamical systems On stability of nonzero set-point for non linear impulsive control systems In this section some basic definitions and results are given concerning the asymptotic stability of sets and Lyapunov theory, in the context of non-linear continuous-time systems ([20] , Appendix B). All the following definitions are referred to systeṁwhere x is the system state constrained to be in X ⊆ R n , f is a Lipschitz continuous nonlinear function, and φ(t; x) is the solution for time t and initial condition x.Definition 7.2 (Attractivity of an equilibrium set).If Y is a ε-neighborhood of X s for some η > 0, we say that X s is locally attractive.We define the domain of attraction (DOA) of an attractive set X s for the system 7.1 to be the set of all initial states x such that φ(t; x) Xs → 0 as t → ∞. We use the term region of attraction to denote any set of initial states contained in the domain of attraction.A closed subset of an attractive set (for instance, a single equilibrium point) is not necessarily attractive. On the other hand, any set containing an attractive set is attractive, so the significant attractivity concept in a constrained system is given by the smallest one 7 . Definition 7.3 (Local ǫ − δ stability of an equilibrium set). Consider system 7.1 constrained by X . A closed equilibrium set X s ⊂ X is ǫ − δ locally stable if for all ǫ > 0 there exists δ > 0 such that if x Xs < δ then φ(t; x) Xs < ǫ, for all t ≥ 0.Unlike attractive sets, a set containing a locally ǫ − δ stable equilibrium set is not necessarily locally ǫ − δ stable. Even more, a closed subset of a locally ǫ − δ stable equilibrium set (for instance, a single equilibrium point) is not necessarily locally ǫ − δ stable. However, any (finite) union of equilibrium sets locally ǫ−δ stable is also locally ǫ − δ stable. So the significant stability concept in a constrained system is given by the largest one.Although a finite union of equilibrium set locally ǫ − δ stable is also locally ǫ − δ stable, in general we cannot extend this result to the case of arbitrary unions of points. Thus, even when every equilibrium point of an equilibrium set is locally ǫ − δ stable, we cannot assure that the whole set would be locally ǫ − δ stable. This is due to the fact that given a fixed ǫ the δ chosen for each point depend on the point and so the infimum of them could be zero. However, if in addition we also assume that the set is compact, then the stability of the set can be inherited from the stability of its points.