key: cord-0909537-gfdi8sik
authors: Abuin, Pablo; Anderson, Alejandro; Ferramosca, Antonio; Hernandez-Vargas, Esteban A.; Gonzalez, Alejandro H.
title: Characterization of SARS-CoV-2 Dynamics in the Host
date: 2020-10-06
journal: Annu Rev Control
DOI: 10.1016/j.arcontrol.2020.09.008
sha: a73e631da0a40731ed59adf778915c286f937259
doc_id: 909537
cord_uid: gfdi8sik

While many epidemiological models were proposed to understand and handle COVID-19, too little has been invested to understand human viral replication and the potential use of novel antivirals to tackle the infection. In this work, using a control theoretical approach, validated mathematical models of SARS-CoV-2 in humans are characterized. A complete analysis of the main dynamic characteristic is developed based on the reproduction number. The equilibrium regions of the system are fully characterized, and the stability of such regions are formally established. Mathematical analysis highlights critical conditions to decrease monotonically SARS-CoV-2 in the host, as such conditions are relevant to tailor future antiviral treatments. Simulation results validates the aforementioned system characterization.

By December 2019, an outbreak of cases with pneumonia of unknown etiology was reported in Wuhan, Hubei province, China [1] . On January 7, a novel betacoronavirus was identified as the etiological agent by the Chinese Center of Disease Control and Prevention (CCDC), and subsequently named as Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) [2] . On February 11, the World Health Organization (WHO) named the disease as Coronavirus disease 2019 (COVID-19) [3] . Although prevention and control measures were implemented rapidly, from the early stages in Wuhan and other key areas of Hubei [4] , the first reported cases outside of China showed that the virus was starting to spread around the world [5].

On March 11, with more that 111.800 cases in 114 countries, and 4921 fatalities cases, COVID-19 was declared a pandemic by the WHO [5] . So far, with more than 7.000.000 total cases confirmed in 213 countries and territories [6, 7] , and an estimated case-fatality rate (CFR) of 5.7% (H1N1 pandemic, CFR<1%) [8] , the potential health risks are evident.

The virus spreads mainly from person-to-person through respiratory droplets produced when an infected person coughs, sneezes or talks [9] . The nonexistence of vaccines or specific therapeutic treatments, preventive measures such as social and physical distancing, hand washing, cleaning and disinfection of surfaces and the use of face masks, among others, have been implemented in order to decrease the transmission of the virus.

Epidemiological mathematical models [10, 11, 12, 13] have been proposed to predict the spread of the disease and evaluate the potential impact of infection prevention and control measures in outbreak management [14] . However, mathematical models at intra-host level that could be useful to understand the SARS-CoV-2 replication cycle and interaction with immune system as well as the pharmacological effect of potential drug therapies [15, 16] are needed. So far, there are approximately 109 trials (including those not yet recruiting, active, or completed) to asses pharmacological therapy for the treatment of COVID-19 in adult patients [17] , including antiviral drugs (i.e. Hydroxychloroquine, Remdesivir, Favipiravir, Lopinavir/Ritonavir, Ribavirin), immunomodulatory agents (i.e. Tocilizumab) and immunoglobulin therapy, among others. Recently, Hernandez-Vargas et. al. [18] proposed different intra-host mathematical models (2 based on target cell-limited model, with and without latent phase, and another considering immune response) for 9 patients with COVID-19. Numerical results in [18] showed intra-host reproductive number values consistent to influenza infection (1.7-5.35) .

Although models in [18] have been fitted to COVID-19 patients data, a control theoretical approach is needed to characterize the model dynamics. Even when the equilibrium states are known, a formal stability analysis is needed to understand the model behavior and, mainly, to design efficient control strategies. Note that the target cell model has been employed previously taking into account pharmacodynamic (PD) and pharmacokinetic (PK) models of antiviral therapies [19, 20] , and this can be potentially done also for COVID- 19. In this context, the main contribution of this article is twofold. First, a full characterization of equilibrium and stability proprieties is performed for the COVID-19 target cell-limited model [18] . Then, formal properties concerning the state variables behavior before convergence -including an analysis of the virus peak times -are given. A key aspect in the target cell model for acute infections shows some particularities such as it has a minimal nontrivial stable equilibrium set, whose stability does not depend on the reproduction number. On the other side, assuming a basic reproduction number greater than 1, the virus would not be cleared before the target cells decreases below under a given critical value, which is independent of the initial conditions.

The article is organized as follows. Section 2 presents the general in-host target cell-limited model used to represent SARS-CoV-2 infection dynamic. Section 3 characterizes the equilibrium sets of the system, and establishes their formal asymptotic stability, by proving both, the attractivity of the equilibrium set in a given domain, and its − δ (Lyapunov) local stability. Then, in Section 4, some dynamical properties of the system are stated, concerning the values of the states at the infection time t = 0. In Section 5 the general model for the SARS-CoV-2 infection is described and the general characteristics of the infection are analyzed. Finally, Section 6 gives the conclusion of the work, while several mathematical formalisms -necessary to support the results of Sections 3 and 4 -are presented in the Appendices.

R and I denote the real and integer numbers, respectively. The real vector space of dimension n is denoted as R n . R n ≥0 represents the vectors of dimension n whose components are equal or greater than zero. The distance from a point x ∈ R n to a set X ⊂ R n is defined by x X := inf z∈X x − z 2 , where · 2 denotes the norm-2. The open ball of radius around a point x ∈ R n , with respect to set X , is defined as B (x) := {z ∈ X : x − z 2 < }. For the real function f (z) = ze z , the so-called Lambert function is defined as the inverse of f (·), i.e., W (z) := f −1 (z) in such a way that W (f (z)) = z.

Although incomplete by definition, mathematical models of in-host virus dynamic improve the understanding of the interactions that govern infections and, more importantly, permit the human intervention to moderate their effects [21] . Basic in-host infection dynamic models usually include the susceptible cells, infected cells, and the pathogen particles [22] . Among the most used mathematical models, the target cell-limited model has been employed to represent and control HIV infection [23, 24, 25] , influenza [26, 27, 28, 19] , Ebola [29] , dengue [30, 31] among others.

In this work, we consider the mathematical model proposed by Hernandez-Vargas [18] given by the following set of ordinary differential equations (ODEs) : [27] .

The parameter values of the target cell model were fitted by [18] using viral kinetics reported by [32] in patients with COVID-19. The Differential Evolution (DE) algorithm was shown to be more robust to initial guesses of parameters than other mentioned methods [33] . Akaike information criterion (AIC) was used to compare the goodness-of-fit for models that evaluate different hypotheses in [18] . The target cell model showed better fitting than exponential growth and logarithmic decay models as well as the target cell model with eclipse phase [18] .

The model (2.1) is non-negative, which means that U (t) ≥ 0, I(t) ≥ 0 and V (t) ≥ 0, for all t ≥ 0. If we denote x(t) := (U (t), I(t), V (t)), then the states are constrained to belong to the invariant set:

Another meaningful set is the one consisting in all the states in X with strictly positive amount of virus and susceptible cells, i.e.,

Note that the set X is an open set. The initial conditions of (2.1) are assumed such at a healthy steady state before the infection time t = 0, i.e., V (t) = 0, I(t) = 0, and U (t) = U 0 , for t < 0. At time t = 0, a small quantity of virions enters to 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 = 0 (formally, V (t) has a discontinuity of the first kind at t 0 , i.e., lim t→0 − V (t) = 0 while lim t→0 + V (t) = V 0 > 0). The same scenario arises, for instance, when an antiviral treatment affects either parameter p or β. The jump of p or β can be considered as a discontinuity of the first kind. In any case, for the time after the discontinuity, the virus may spread or be cleared in the body, depending on its infection effectiveness. The following (mathematical) definition is given Definition 1 (Spreadability of the virus in the host). Consider system (2.1), constrained by the positive set X, at some time t 0 , with U (t 0 ) > 0, I(t 0 ) ≥ 0 and V (t 0 ) > 0 (i.e., x(t 0 ) = (U (t 0 ), I(t 0 ), V (t 0 )) ∈ X ). Then, it is said that the virus spreads in the host for t > t 0 if there exists at least one t * > t 0 such thaṫ

The latter definition states that the virus spreads in the body host if V (t) has at least one local maximum. On the other hand, the virus does not spread if V (t) is strictly decreasing for all t > t 0 . As it will be stated later on (Property 1), lim t→∞ V (t) = 0 for system (2.1), independently of the fact that the virus reaches or not a maximum (this is a key difference between acute and chronic infection models [21, 22] ).

The infection severity can be related with the virus spreadability established in Definition 1. Liu et.al. [34] have shown that patients with severe COVID-19 tend to have a high viral load and a long virus shedding period. The mean viral load of severe cases was around 60 times higher than that of mild cases, suggesting that higher viral loads might be associated with severe clinical outcomes. Furthermore, they found that the viral load of severe cases remained significantly higher for the first 12 days after the appearance of the symptoms than those of corresponding mild cases. Mild cases were also found to have an early viral clearance, with 90% of these patients repeatedly testing negative on reverse transcription polymerase chain reaction (RT-PCR) by day 10 post symptoms onset (pso). By contrast, all severe cases still tested positive at or beyond day 10 pso. In addition, Zheng et.al. [35] reported (from a study with 96 SARS-CoV-2 patients, 22 with mild and 74 with severe disease) a longer duration of SARS-CoV-2 in lower respiratory samples of severe patients. For patients with severe disease the virus permanence was significantly longer (21 days, 14-30 days) than in patients with mild disease (14 days, 10-21 days; p=0.04). Moreover, higher viral loads were detected in respiratory samples, although no differences were found in stool and serum samples. While these findings suggest that reducing the viral load through clinical means and strengthening management should help to prevent the spread of the virus, they are preliminary and it remains controversial whether virus persistence is necessary to drive the dysfunctional immune response characteristic of COVID-19 patients [36] .

To properly establish conditions under which the virus does not spread for t > 0 (i.e., after the infection time t = 0) the so-called intra-host basic reproduction number is defined next.

Particularly, for t = 0, this number describes the number of infected cells produced by one infected cell, when a small amount of virus, V 0 , is introduced into a healthy stationary population of uninfected target cells, U 0 ,

A discussion about the way this value is obtained is given in Appendix 2. The relation between the basic reproduction number at the infection time (R 0 ) and the virus spreadiblity is stated in the next theorem.

Theorem 2.1. Consider the system (2.1), constrained by the positive set X, at the beginning of the infection, i.e., U (0) = U 0 > 0, I(0) = 0 and V (0) = V 0 > 0 (i.e., x(0) = (U (0), I(0), V (0)) ∈ X ). Then, a sufficient condition (not necessary) for the virus not to spread in the host is given by R 0 < 1.

Proof: In Theorem 4.1, Section 4, it is shown that if the virus spreads, then R 0 > 1. This means that (contrapositive of the statement) if R 0 ≤ 1 (particularly, R 0 < 1), then the virus does not spread in the host body.

Before proceeding with a full dynamic analysis of system (2.1), let us define first the so-called critical value of the susceptible cells, which is a threshold to properly understand the spread of the virus.

which, for fixed system parameters β, p, δ and c, is a constant.

Note that U (t) < U c if and only if R(t) < 1, for every t ≥ 0.

By equatingU ,İ andV to zero in (2.1), it can be shown that the system only has healthy equilibria of the form x s = (U s , 0, 0), with U s being an arbitrary positive value, i.e., U s ∈ [0, ∞). Thus, there is only one equilibrium set, which is the disease-free one, and it is defined by

To examine the stability of the equilibrium points in X s , system (2.1) can be linearized at a general state

. Then, the Jacobian matrix is given by

which evaluated at any point x s ∈ X s reads

The first eigenvalue is trivially given by λ 1 = 0. The other two, are given by:

To analyze the eigenvalues qualitatively, note that for U s = U c

which means that λ 2 = 0 and λ 3 = −(c + δ) < 0 (given that c, δ > 0). Furthermore, λ 2 < 0 and λ 3 < 0 for U s < U c ; and λ 2 > 0 and λ 3 < 0 for U s > U c . Since the maximum eigenvalue is the one dominating the stability behavior of the equilibrium under consideration, it is possible to infer how the system behaves near some segments of X s . The first intuition is that the equilibrium set

is stable, and that the equilibrium set

is unstable. These are just intuitions, 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 the nonlinear system (Theorem of Hartman-Grobman [37, 38] ). To formally prove the asymptotic stability of X 1 s in a given domain, it is necessary to show its global attractivity (in such domain) and local -δ stability.

A key point to analyze the general asymptotic stability (AS) of system (2.1) is to consider stability of the complete equilibrium sets X 1 s and X 2 s , and not of the single points inside them (as defined in Definitions 5, 6 and 7, in Appendix 1). As it is shown in the next subsections, there is no single AS equilibrium points in this system, although there is an AS equilibrium set (i.e., X 1 s ). As stated in Definition 7, in Appendix 1, the AS of X 1 s requires both, attractivity and − δ stability, which are stated in the next two subsections, respectively. Then, in Subsection 3.3 the AS theorem is formally stated.

s in X Before proceeding with the formal theorems of the attractivity of X 1 s , let us consider the following key property of system (2.1) concerning the attractivity of X s . Property 1 (Attractivity of X s ). Consider system (2.1) constrained by the positive set X, at some arbi-

Property 1 states that X s is an attractive set for system (2.1), in X , but not the smallest attractive set. Next, conditions are given to show that the smallest attractive set is given by X 1 s . Theorem 3.1 (Attractivity of X 1 s ). Consider system (2.1) constrained by the positive set X. Then, the set X 1 s defined in (2.8) is the smallest attractive set in X . Furthermore, X 2 s , defined in (2.9), is not attractive. Proof: The proof is divided into two parts. First it is proved that X 1 s is an attractive set, and then, that it is the smallest one.

Attractivity of X 1 s : The attractivity of X s in X was already proved in Property 1. So, to prove the attractivity of X 1 s in X (and to show that X 2 s is not attractive) it remains to demonstrate that U ∞ ∈ [0, U c ). From system (2.1), by replacing (2.1a) in (2.1b), it follows thaṫ

which implies that

Rearranging (2.1c) yields

This latter equation can be integrated, for general initial conditions U 0 , I 0 and V 0 , as follows:

Now, by defining U ∞ := lim t→∞ U (t), I ∞ := lim t→∞ I(t), V ∞ := lim t→∞ V (t), and recalling from Property 1 that I ∞ = V ∞ = 0, the latter equation for t → ∞, reads

where R 0 := βp cδ U 0 (as it was defined in (2.5)) and

Note that R 0 is a function of U 0 while K 0 is a function of I 0 and V 0 , and, furthermore, R 0 > 0 and

Now, by denoting

and

the latter equation can be written as

where W (·) is a Lambert function. Figure 1 shows the graph of such a function, where it can be seen that it has two branches, denoted as W p and W m . However, Figure 2 shows a plot of function z(R 0 , K 0 ) for negative values of K 0 and positive values of R 0 ), and W p maps (−1/e, 0] into (−1, 0], which implies that

for R 0 > 0 and K 0 < 0. Thus, by (3.13) , it follows that

which completes the proof. X 1 s is the smallest attractive set: It is clear from the previous analysis, that any initial state

. This means that X 2 s is not attractive in X . Let us consider now a state x s ∈ X 1 s and an arbitrary small ball of radius > 0, w.r.t. X , around it, B (x s ) ∈ X . Take two arbitrary initial states x 0,1 = (U 0,1 , I 0,1 , V 0,1 ) and x 0,2 = (U 0,2 , I 0,2 , V 0,2 ) in B (x s ), such that U 0,1 = U 0,2 and V 0,1 = V 0,2 . These two states converge, according to equation (3.15) , to x ∞,1 = (U ∞,1 , 0, 0) and x ∞,2 = (U ∞,2 , 0, 0), respectively, with U ∞,1 , U ∞,2 ∈ [0, U c ). Given that function z(R, K) is monotone (injective) in R 0 (and so in U 0 ) and W (z) is monotone (injective) in z, then U ∞,1 = U ∞,2 . This means that, although both initial states converge to some state in X 1 s , they necessarily converge to different points. Therefore neither single states x s ∈ X 1 s nor subsets of X 1 s are attractive in X . So, X 1 s is the smallest attractive set and the proof is concluded. Remark 2. Note that X 1 s and X 2 s are in the closure of the open set X , which is not in X . In other words, Theorem 3.1 shows that any initial state in X converges to a point onto the boundary of X that does not belong to X . Furthermore note that, an initial state of the form (U 0 , 0, 0), U 0 > U c , (i.e., a state in X 2 s ) cannot be attracted by any set since it is an equilibrium state (every state in X 2 s will remains unmodified). This is the reason why it is not possible to consider the attractivity of X 2 s in X . Proof: Let us consider a particular equilbrium point x s := (U s , 0, 0), with U s ∈ [0, U c ) (i.e., x s ∈ X 1 s ). Then a Lyapunov function candidate is given by (similar to one used in [39] for chronic infections)

This function is continuous in X, is positive for all nonegative x = x s and, furthermore, J(x s ) = 0. Function J evaluated at the solutions of system (2.1) reads:

Now, given U s ∈ [0, U c ), with U c = δc βp , it follows thatJ(x(t)) ≤ 0 for every x ∈ X (note that it is not true thatJ(x(t)) < 0 for x = x s , as shown next, in Remark 3). Then, J is a Lyapunov function for system (2.1), which means that each x s ∈ X 1 s is − δ stable (see Theorem 6.1 in Appendix 1). Therefore, it is easy to see that the equilibrium set X 1 s is also − δ stable, which completes the proof. Remark 3. Note that, in the latter proof, it is not true thatJ(x(t)) < 0 for every nonegative x = x s . If for instance, the functionJ(x(t)) is evaluated atx s = (Û , 0, 0), withÛ = U s , we have thatJ(x s (t)) = 0.

s are − δ stable, but not attractive. A schematic plot of such a behavior can be seen in Figure 3 . s is − δ stable but not attractive. Initial states x0 starting arbitrarily close to xs remain (for all t ≥ 0) arbitrarily close to xs, but do not converge to xs. As a consequence, set X 1 s is AS but the points inside it are not. [40, 41] :Ṡ = βSI,İ = βSI − δI, being S the susceptible and I the infected individuals. In this latter model, R 0 := (δ/β)S 0 and the critical value for S is S c = δ/β. The AS set is given by all the states of the form x s := (S s , 0), with S s ∈ [0, S c ). Furthermore, for this system, the maximum of I occurs when S = S c .

= Ax, when A = [0 − 1; 0 − 1 ], or the 2-state Kermack-McKendrick epidemic model

In the next Theorem, based on the previous results concerning the attractivity and − δ stability of X 1 s , the asymptotic stability is formally stated. Theorem 3.3. Consider system (2.1) constrained by the positive set X. Then, the set X 1 s defined in (2.8) is smallest asymptotically stable (AS) equilibrium set, with a domain of attraction given by X .

Proof: The proof follows from Theorems 3.1, which states that X 1 s is the smallest attractive in X , and 3.2, which states the local − δ stability of X 1 s .

A critical consequence of the latter Theorem is that no equilibrium point in X s (neither in X 1 s , nor in X 2 s ) can be used as setpoint in a control strategy design. The effect of antivirals (pharmocodynamic), for instance, is just to reduce the virus infectivity (by reducing the infection rate β) or the production of infectious virions (by reducing the replication rate p) [21] . So, the previous stability analysis is still valid for such controlled systems, since only a modification of some of the parameters defining U c is done. In such a context, only a controller able to consider the whole set X 1 s as a target (a set-based control strategy, as zone MPC [42, 43] ) will be fully successful in controlling system (2.1). Further details concerning antiviral treatments are given next, in Section 4.1.

In this section some further properties of system (2.1) concerning its dynamic are stated, based on the initial conditions at the infection time t = 0. The objective is to fully characterize the states behavior in a qualitative way, including the times at which the virus and the infected cells reach their peaks. First, Property 2 states some characteristics of U ∞ for different initial conditions. Then, Theorem 4.1 states a general relationship between the peak times of V and I and the time at which U reaches its critical value U c . 

iii. z(R 0 ) = −R 0 e −R0 is strictly decreasing for R 0 ∈ (0, 1) (note that R 01 := cδU0,1 βp and R 02 := cδU0,2 βp are in (0, 1), since they are smaller than U c ), while −W p (·) is strictly decreasing in (−1/e, 0) . Figure 4 shows U ∞ as a function of U 0 , taking V 0 as a parameter. Even more, by hypothesis the virus spreads, which means that V (t) reaches a local maximum at some timet V > 0. Therefore, V (t) must reach a local minimum at some 0 <ť V <t V . Now, by Lemma 1 in Appendix 3, it is R(ť V ) > 1 and R(t V ) < 1, respectively, and it is easy to see that R(t) is a decreasing function, so it follows that R 0 > R(ť V ) > 1. Then there exists α(0) > 0 such that R 0 > 1 + α(0) and, besides, 0 <ť V < t c <t V . From the minimum and maximum conditions of V , at timesť V and t V , we haveV (ť V ) = 0,V (ť V ) > 0 andV (t V ) = 0,V (t V ) < 0, respectively. After some algebraic computation, it is easy to see thatİ(ť V ) > 0 andİ(t V ) < 0, which means that I(t) must reach a maximum at some timet I , fulfillingť V <t I <t V . Moreover, it must bė

Given thatV (t) > 0 forť V < t <t V (it goes from its minimum to its maximum), then by (2.1.a),

Replacing this latter condition in (4.1), it follows that

which implies that R(t I ) = βpU (tI ) δc > 1 and, then,t I < t c . Therefore, t 0 <ť V <t I < t c <t V , which concludes the proof.

Remark 5. The value of α(0) is necessary to properly understand and characterize the system behavior according to the initial conditions and parameters. Opposite to what happens, for instance, in epidemiological models (SIR, etc.) , where R 0 > 1 is a necessary and sufficient condition for the disease to spread in a population, in our case R 0 > 1 is not a sufficient condition for the virus to spread in the host body. The only thing Theorem 4.1 ensures (by its contrapositive) is that a sufficient condition for the virus to not spread in the host body at time t > 0 is given by R 0 < 1 (or U (0) < U c ). See Figure 6 , lower plot, for an example. The value of α(0) can be computed numerically and it is usually small in comparison with R 0 (for all the patients simulated in Section 5, α(0) < 1 × 10 −4 ).

To clarify the results of this section, Figures 5 and 6 show a phase portrait and a state time evolution corresponding to system (2.1), when all parameters are equal to 1 (for simplicity), which means that U c = 1. The first plot ( Figure 5 ) depicts how every state trajectory -even those starting close to X 2 sconverges to X 1 s . As stated in Property 2, U ∞ approaches U c from below, as U (0) approaches U c from above. Also it can be seen how the virus load starts do decrease only once U (t) is smaller than U c , as stated in Theorem 4.1. On the other hand, the second plot ( Figure 6) shows the time evolution of U , I and V , for two different initial conditions. In the upper plot, initial conditions are selected such that 1 + α(0) < R 0 , while in the lower plot, the initial conditions produce 1 < R 0 < 1 + α(0). As it can be seen, only in the first case the virus spread in the host body (i.e.,V (t) > 0, for some t > 0), as stated in Theorem 4.1. 

Even though the analysis of potential antiviral treatments is out of the scope of this work, in this section some comments concerning the implications of Theorem 4.1 (and the system characterization) will be made. The antiviral effect can be modeled as a reduction of the virus infectivity in the presence of reverse transcriptase inhibitors (by reducing the infection rate β) and/or as a reduction in the production of infectious virions in the presence of protease inhibitors (by reducing the replication rate p). Let us assume that the antiviral pharmacodynamics (PD) corresponding to an antiviral is modeled as p(1 − η(t r )) (the analysis for β is almost the same), being η(t r ) ∈ (0, 1) the effectiveness of the antiviral and t r the time of treatment initiation. The antiviral pharmacokynetics (PK) is not considered, for simplicity, which means that the antivirals instantaneously modify η at time t r . Then, as the virus monotonically goes to zero only once U (t) is below U c , the antiviral will be effective (in the sense that the virus load starts decreasing as the treatment begins, and it does not increase again) only if the value of η(t r ) is such that U (t r ) < U c (t r ) := cδ p(1−η(tr))β (i.e., such that R(t r ) < 1 + α(t r ) ≈ 1). This condition defines a threshold for the antiviral effectiveness (say, a minimal critical value η c (t r )) that may explain, from a pure mathematical point of view, why some antiviral may not work for some patients.

From a control theory point of view, the assertions made in Theorem 4.1 means that a control strategy devoted to steers V (t) to zero at any time by administering a time-variant dose of antivirals (for instance by using η(t) < η c (t), for t > t r ), may be counterproductive. Indeed, to slow down V (t) by decreasing p or β, implies that U c = cδ pβ increases, but also soften the decreasing behavior of U (t). As a result, the time t c (and so, the virus peak timet V ) may be delayed, which means that V (t) is maintained in a high level for a longer time. According to preliminary simulations, the delay of the virus peak may be significantly long for antiviral with maximal effectiveness smaller than the critical value.

A formal (and detailed) mathematical analysis of the antiviral effect in COVID-19 is the main topic of a future work.

In this section, the model parameters in (2.1) will be associated to the patients labeled as A, B, C, D, E, F, G, H and I -reported in [32] . The initial number of target cells U 0 is estimated as approximately 10 7 cells [18] . I 0 is assumed to be 0 while V 0 is determined by interpolation considering an incubation period of 7 days (note, that V 0 ranges from 0.02 to 5.01 copies/mL which is below the detectable level of about 100 copies/mL). Moreover, the onset of the symptoms is assumed to occurs 4 to 7 days after the infection time (day 0, Figure 7 and 8) . The parameters and the initial conditions (U 0 , I 0 and V 0 , with t 0 = 0 the infection time) of each patient are collected in Table 1 . According to the system analysis of the previous sections, some relevant dynamical values are shown in Table 2 . Constant α(0) (defined in Theorem 4.1) is smaller than 10 × 10 −4 for all the patients, so it is not taken into account for the study. Figures 7 and 8 show the dynamics of V and U . As expected, the states converge to X 1 s , although significantly different behaviors can be observed for the different patients. From Figure 8 it can be seen that the healthy cells final value U ∞ is reduced in cases of patients with large values of R 0 , in spite all simulations have the same initial U 0 . This can be explained from the fact that W (R 0 e −R0 e K0 ) is monotonically decreasing for R 0 > 1 (see Figures 1 and 2) , and therefore, 0 < U ∞ (R 01 ) < U ∞ (R 02 ) for R 01 > R 02 > 1 (see Property 2, above). Note that the susceptible cells of patient C converges to U ∞ equals to 4.810 × 10 −10 [cell], which can be explained by the fact that this patient has a reproduction number (R 0 ) of 37.57, which is 5.2 times above the cohort mean value of 7.21. Figure 7 and Table 2 show that the viral load of patient C reaches the peak at 1.69 days post infection (dpi) (40.56 hours post infection, hpi).

Furthermore, from Figure 7 , it can be seen that for all the cases the viral load spreads (i.e.: the virus presents a peak) although R V (0) < 0 for all patients (i.e., I 0 = 0). This can be justified since U 0 U c and, therefore, R 0 will be greater than 1 + α(0) for all patients (note that, α(0) < 10 × 10 −4 ). Moreover, from Table 2 , we can corroborate thatt I > t c >t V which is in accordance to what is stated in Theorem 4.1.

Concerning the immune response, this model makes the assumption that it is constant and independent on viral load as well as infected cells. Furthermore, neither innate or adaptive response are modeled, being the viral load dynamic mainly limited by target cells availability. Since recent studies have shown a dysfunctional immune response (i.e.: lymphogenia, desregulated secretion of pro-inflammatory cytokines, excessive infiltration of monocytes, macrophages and T cells, among others) [36, 44] , this effect should be added in the proposed model, in order to have a more reliable representation (and, eventually, a more realistic control objective). In addition, a more reliable standard to measure the severity of disease could be related with the viral spreadability as well as the deregulated inflammatory response. [18] . The patient labeling is as presented in [32] . Vclear denotes a value of 50 [copies/ml] under which the virus is not detectable and it is considered is cleared. [18] . The patient labeling is as presented in [32] . Simulation for the patient C shows a very low value of U∞ (practically zero), which suggests that the selected value of U0 = 1.0e 7 may be large.

In this work a full dynamical characterization of a COVID-19 in-host target-cell model is performed. It is shown that there exists a minimal stable equilibrium set depending only on the system parameters. Furthermore, it is shown that there exists a parameter-depending threshold for the susceptible cells that fully characterizes the virus and infected cells qualitative behavior. Simulations demonstrate the potential utility of such system dynamic characterization to tailor the most valuable pipeline drugs against SARS-CoV-2.

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. All the following definitions are referred to systemẋ

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 4 (Equilibrium set). Consider system 6.1 constrained by X. The set X s ⊂ X is an equilibrium set if each point x ∈ X s is such that f (x) = 0 (this implying that φ(t; x) = x for all t ≥ 0).

Definition 5 (Attractivity of an equilibrium set). Consider system 6.1 constrained by X. A closed equilibrium set X s ⊂ X is attractive in X ⊂ X if lim t→∞ φ(t; x) Xs = 0 for all x ∈ X .

Any set containing an attractive set is attractive, so the significant attractivity concept in a constrained system is given by the smallest one.

Definition 6 ( − δ local stability of an equilibrium set). Consider system 6.1 constrained by X. A closed equilibrium set X s ⊂ X is − δ locally stable if for all > 0 it there exists δ > 0 such that in a given boundary of X s , x Xs < δ, it follows that φ(t; x) Xs < , for all t ≥ 0. Definition 7 (Asymptotic stability (AS) of an equilibrium set). Consider system 6.1 constrained by X. A closed equilibrium set X s ∈ X is asymptotically stable (AS) in X ⊂ X if it is − δ locally stable and attractive in X . Theorem 6.1 (Lyapunov theorem [45] ). Consider system 6.1 constrained by X and an equilibrium state

) ≤ 0, denoted as Lyapunov function. Then, the existence of such a function implies that x s ∈ X s is − δ locally stable. If in additionV (x(t)) < 0 for all x = x s andV (x s ) = 0, then x s ∈ X s is asymptotically stable.

The derivation of the basic reproduction number R 0 will be given by means of the concept of nextgeneration matrix [46] . Consider system (2.1) and the healthy equilibrium x 0 = (U 0 , 0, 0), which is stable in the absence of virus. Of the complete state of system (2.1), x = (U, I, V ), only two states depend on infected cells, that is I and V . Let us rewrite the ODEs for this two states in the forṁ

where 

then matrix F G −1 , represents the so-called next-generation matrix. Each (i, j) entry of such a matrix represents the expected number of secondary infections in compartment i produced by an infected cell introduced in compartment j. The spectral radius of this matrix, that is, the maximum absolute value of its eigenvalues, defines the basic reproduction number R 0 . For the specific case of system (2.1), the next-generation matrix is given by

Therefore, the basic reproduction number R 0 is given by

Notice that R 0 coincides with the entry (1, 1) of matrix F G −1 , thus meaning that R 0 represents the expected number of secondary infections produced in compartment I by an infected cell originally in I.

The next Lemma characterizes the virus minimum and maximum times, for system (2.1), in terms of the value of the reproduction number R(t). Given that I(t * V ) > 0 (note that I(t) is positive for all t > 0), then βp c U (t * V ) − δ = 0, or

This way if an inflection point does occurs at t * V , then t * V = t c , where t c is the time at which R = 1. This proves item (iii).

Furthermore, if V reaches a local minimum at t * V , thenV (t * V ) > 0 (instead ofV (t * V ) = 0, as it is in (6.4), which by (6.3) implies that

This proves item (i).

On the other hand, if V reaches a local maximum at t * V , thenV (t * V ) < 0 (instead ofV (t * V ) = 0, as it is in (6.4)), which by (6.3) implies that

This proves item (ii).

by substituting (3.4) in (2.1a), and multiplying by 1/U (t) both sides of the equation

Outbreak of pneumonia of unknown etiology in wuhan china: the mystery and the miracle

Severe acute respiratory syndrome-related coronavirus-the species and its viruses, a statement of the coronavirus study group

Who director-general's remarks at the media briefing on

Report of the who-china joint mission on coronavirus disease

Coronavirus disease 2019 (covid-19) situation report -86

Covid-19 dashboard by the center for systems science and engineering (csse) at johns hopkins university

Who director-general's opening remarks at the media briefing on covid

A sidarthe model of covid-19 epidemic in italy

The sars-cov-2 epidemic outbreak: a review of plausible scenarios of containment and mitigation for mexico

Novel coronavirus 2019-ncov: early estimation of epidemiological parameters and epidemic predictions

Neural Control for Epidemic Model of Covid-19 with a Complex Network Approach

How will country-based mitigation measures influence the course of the covid-19 epidemic?

Research and development on therapeutic agents and vaccines for covid-19 and related human coronavirus diseases

Use of antiviral drugs to reduce covid-19 transmission, The Lancet Global Health

-host modelling of COVID-19 kinetics in humans

Passivity-based inverse optimal impulsive control for influenza treatment in the host

Oseltamivir pk/pd modeling and simulation to evaluate treatment strategies against influenza-pneumococcus coinfection, Frontiers in cellular and infection microbiology

Modeling and Control of Infectious Diseases in the Host: With MATLAB and R

In-host modeling

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

Influenza virus population dynamics in the respiratory tract of experimentally infected mice

Kinetics of influenza a virus infection in humans

Influenza a virus infection kinetics: quantitative data and models

Ebola virus infection modeling and identifiability problems

The role of antibody in enhancing dengue virus infection

Modelling original antigenic sin in dengue viral infection

Virological assessment of hospitalized patients with COVID-2019

A Comparative study of Differential Evolution Algorithms for Parameter Fitting Procedures, IEEE World Congress on Computational Intelligence (WCCI)

Viral dynamics in mild and severe cases of covid-19

Viral load dynamics and disease severity in patients infected with sars-cov-2 in zhejiang province, china

The trinity of covid-19: immunity, inflammation and intervention

Differential equations and dynamical systems

Global stability analysis of the original cellular model of hepatitis c virus infection under therapy

Mathematical models in population biology and epidemiology

The kermack-mckendrick epidemic model revisited

MPC for tracking zone regions

Stable impulsive zone mpc for type 1 diabetic patients based on a long-term model

Reduction and functional exhaustion of t cells in patients with coronavirus disease 2019 (covid-19)

Reproduction numbers of infectious disease models

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests