Characterization of SARS-CoV-2 dynamics in the host.

While many epidemiological models were proposed to understand and handle COVID-19 pandemic, 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 is 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 show the aforementioned system characterization.

Introduction

By December 2019, an outbreak of cases with pneumonia of unknown etiology was reported in Wuhan, Hubei province, China (Lu, Stratton, & Tang, 2020). 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) (Gorbalenya, 2020). On February 11, the World Health Organization (WHO) named the disease as Coronavirus disease 2019 (COVID-19) (Who, 2020). Although prevention and control measures were implemented rapidly, from the early stages in Wuhan and other key areas of Hubei Who, 2020, the first reported cases outside of China showed that the virus was starting to spread around the world (whotimeline, 2020).

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 (whotimeline, 2020). So far, with more than 7.000.000 total cases confirmed in 213 countries and territories (Coronavirus disease 2019, COVID-19), and an estimated case-fatality rate (CFR) of 5.7% (H1N1 pandemic, CFR < 1%) (Who, 2020), 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 (How covid-19 spreads). 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 (Acuna-Zegarra, Comas-Garcia, Hernandez-Vargas, Santana-Cibrian, Velasco-Hernandez, 2020, Alanis, Member, Hernandez-vargas, Nancy, Ríos-rivera, 2020, Giordano et al., 2020, Read, Bridgen, Cummings, Ho, Jewell, 2020) have been proposed to predict the spread of the disease and evaluate the potential impact of infection prevention and control measures in outbreak management (Anderson, Heesterbeek, Klinkenberg, & Hollingsworth, 2020). However, mathematical models at in-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 (Liu, Zhou, Li, Garner, Watkins, Carter, Smoot, Gregg, Daniels, Jervey, et al., Mitjà, Clotet, 2020) 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 (Sanders, Monogue, Jodlowski, & Cutrell, 2020), including antiviral drugs (i.e. Hydroxychloroquine, Remdesivir, Lopinavir/Ritonavir, Ribavirin), immunomodulatory agents (i.e. Tocilizumab) and immunoglobulin therapy, among others. Recently, Hernandez-Vargas & Velasco-Hernandez, 2020 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 Hernandez-Vargas & Velasco-Hernandez, 2020 showed intra-host reproductive number values consistent to influenza infection (1.7-5.35).

Although models in Hernandez-Vargas & Velasco-Hernandez, 2020 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 (Boianelli, Sharma-Chawla, Bruder, Hernandez-Vargas, 2016, Hernandez-Mejia, Alanis, Hernandez-Gonzalez, Findeisen, Hernandez-Vargas, 2019), 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 (Hernandez-Vargas & Velasco-Hernandez, 2020). 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 . 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.

Notation

and denote the real and integer numbers, respectively. The real vector space of dimension n is denoted as . represents the vectors of dimension n whose components are equal or greater than zero. The distance from a point to a set is defined by where ‖ · ‖2 denotes the norm-2. The open ball of radius ϵ around a point with respect to set is defined as . For the real function the so-called Lambert function is defined as the inverse of f( · ), i.e., in such a way that .

SARS-CoV-2 in-host mathematical model

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 (Hernandez-Vargas, 2019). Basic in-host infection dynamic models usually include the susceptible cells, infected cells, and the pathogen particles (Ciupe & Heffernan, 2017). Among the most used mathematical models, the target cell-limited model has been employed to represent and control HIV infection (Legrand, Comets, Aymard, Tubiana, Katlama, Diquet, 2003, Perelson, Kirschner, De Boer, 1993, Perelson, Ribeiro, 2013), influenza (Baccam, Beauchemin, Macken, Hayden, Perelson, 2006, Hernandez-Mejia, Alanis, Hernandez-Gonzalez, Findeisen, Hernandez-Vargas, 2019, Larson, Dominik, Rowberg, Higbee, 1976, Smith, Perelson, 2011), Ebola (Nguyen, Binder, Boianelli, Meyer-Hermann, & Hernandez-Vargas, 2015), dengue (Nikin-Beers, Ciupe, 2015, Nikin-Beers, Ciupe, 2018) among others.

In this work, we consider the mathematical model proposed by Hernandez-Vargas & Velasco-Hernandez, 2020 given by the following set of ordinary differential equations (ODEs) : where U [cells], I [cells] and V [copies/mL] represent the susceptible cells, the infected cells, and the virus load, respectively. The parameter β is the infection rate of susceptible cells by the virus. δ is the death rate of infected cells. p is the replication rate of free virus from the infected cells. c is the degradation (or clearance) rate of virus V. The effects of immune responses are not explicitly described in this model, but they are implicitly included in the death rate of infected cells (δ) and the clearance rate of virus (c) (Baccam et al., 2006).

The parameter values of the target cell model were fitted by Hernandez-Vargas & Velasco-Hernandez, 2020 using viral kinetics reported by Wölfel et al. (2020) 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 (Torres-Cerna, Alanis, Poblete-Castro, Bermejo-Jambrina, & Hernandez-vargas, 2016). Akaike information criterion (AIC) was used to compare the goodness-of-fit for models that evaluate different hypotheses in Hernandez-Vargas & Velasco-Hernandez, 2020. The target cell model showed better fitting than exponential growth and logarithmic decay models as well as the target cell model with eclipse phase (Hernandez-Vargas & Velasco-Hernandez, 2020).

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 with strictly positive amount of virus and susceptible cells, i.e.,Note that the set is an open set.

The initial conditions of (2.1) are assumed such at a healthy steady state before the infection time i.e., and for t < 0. At time 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 (formally, V(t) has a discontinuity of the first kind at t 0, i.e., while ). 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

Spreadability of the virus in the host

Consider the system (2.1), constrained by the positive set at some time t 0, with U(t 0) > 0, I(t 0) ≥ 0 and V(t 0) > 0 (i.e., ). 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 that .

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), 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 (Ciupe, Heffernan, 2017, Hernandez-Vargas, 2019)).

The infection severity can be related with the virus spreadability established in Definition 1. Liu et al. (2020b) 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. (2020) 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 (Tay, Poh, Rénia, MacAry, & Ng, 2020).

Note that the virus spreadability may or may not cause a severe infection (a disease that eventually causes host death) which depends on how much time the virus is above a given value.

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

The intra-host basic reproduction number is defined as the number of infected cells (or virus particles) that are produced by one infected cell (or virus particle), at a given time. Its mathematical expression is given by:

Particularly, for 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 () and the virus spreadiblity is stated in the next theorem.

Consider the system (2.1) , constrained by the positive set at the beginning of the infection, i.e., and (i.e., ). Then, a sufficient condition (not necessary) for the virus not to spread is given by .

In Theorem 4.1, 4, it is shown that if the virus spreads, then . This means that (contrapositive of the statement) if (particularly, ), 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.

The critical value for U, is defined aswhich, for fixed system parameters β, p, δ and c, is a constant.

Note that if and only if for every t ≥ 0.

Equilibrium set characterization

By equating and to zero in (2.1), it can be shown that the system only has healthy equilibria of the form with U being an arbitrary positive value, i.e., U ∈ [0, ∞). Thus, there is only one equilibrium set, which is the disease-free one, and it is defined byTo examine the stability of the equilibrium points in system (2.1) can be linearized at a general state . From (2.1) we have . Then, the Jacobian matrix is given bywhich evaluated at any point readswith U ∈ [0, ∞). Then, the eigenvalues (λ 1, λ 2, λ 3) are given by the solution to i.e.,The first eigenvalue is trivially given by . The other two, are given by:To analyze the eigenvalues qualitatively, note that for which means that and (given that c, δ > 0). Furthermore, λ 2 < 0 and λ 3 < 0 for ; and λ 2 > 0 and λ 3 < 0 for . 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 . The first intuition is that the equilibrium setis stable, and that the equilibrium setis 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 (1982); Perko (2013)). To formally prove the asymptotic stability of in a given domain, it is necessary to show its global attractivity (in such domain) and local stability.

Asymptotic stability of the equilibrium sets

A key point to analyze the general asymptotic stability (AS) of system (2.1) is to consider stability of the complete equilibrium sets and and not of the single points inside them (as defined in Definition 5, Definition 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., ).

As stated in Definition 7, in 1, the AS of requires both, attractivity and stability, which are stated in the next two subsections, respectively. Then, in Section 3.3 the AS theorem is formally stated.

Attractivity of set in

Before proceeding with the formal theorems of the attractivity of let us consider the following key property of system (2.1) concerning the attractivity of .

(2.1) constrained by the positive set at some arbitrary time t 0, with U(t 0) > 0, I(t 0) ≥ 0 and V(t 0) > 0 (i.e., ). Then, U ∞ ≔ lim U(t) is a constant value smaller than U(t 0), and which means that tends to some state in .

Since for all t ≥ 0 and all by (2.1a) U(t) is a decreasing function (no oscillation can occur). Since U(t 0) > 0 and V(t 0) > 0, then is a constant value in [0, U(t 0)). Given that U(t) converges to a finite fixed value, then as t → ∞ by (2.1a). This implies - by the same Eq. (2.1a) - that as t → ∞ and, so, from Eq. (2.1b), that as t → ∞. Then . Finally, by Eq. (2.1c)), as t → ∞. Then which completes the proof. □

Property 1 states that is an attractive set for system (2.1), in but not the smallest attractive set. Next, conditions are given to show that the smallest attractive set is given by .

(2.1) constrained by the positive set . Then, the set defined in (2.8) is the smallest attractive set in . Furthermore, defined in (2.9) , is not attractive.

The proof is divided into two parts. First it is proved that is an attractive set, and then, that it is the smallest one.

Attractivity of :

The attractivity of in was already proved in Property 1. So, to prove the attractivity of in (and to show that is not attractive) it remains to demonstrate that . From system (2.1), by replacing (2.1a) in (2.1b), it follows that which implies that

Rearranging (2.1c) yields

Then, replacing (3.2) in (3.3), we have

Finally, by substituting (3.4) in (2.1a), and multiplying by 1/U(t) both sides of the equation (without loss of generality we can assume that ), it follows that

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

Now, by defining U ∞ ≔ lim U(t), I ∞ ≔ lim I(t), V ∞ ≔ lim V(t), and recalling from Property 1 that the latter equation for t → ∞, reads where (as it was defined in (2.5)) and

Note that is a function of U 0 while is a function of I 0 and V 0 and, furthermore, and for every . Then, after some manipulation, (3.7) reads

Now, by denotingand the latter equation can be written as or where W( · ) is a Lambert function. Fig. 1 shows the graph of such a function, where it can be seen that it has two branches, denoted as W and W. However, in this case, since for which has not biological sense. Note that U ∞ is a finite value in [0, U 0). Also for and (Fig. 2 shows a plot of function for negative values of and positive values of ), and W maps into which implies that for and . Thus, by (3.13), it follows that which completes the proof.

Fig. 1

Lambert function. W(z) has two branches, denoted as W (in blue) and W (in red). Both branches are defined for ; however while which means that only the branch W will be used in our analysis, as it is shown in the proof of Theorem 3.1.

Fig. 2

Function for and . Note that for all values of and .

is the smallest attractive set:

It is clear from the previous analysis, that any initial state in converges to a state with . This means that is not attractive in . Let us consider now a state and an arbitrary small ball of radius ϵ > 0, w.r.t. around it, . Take two arbitrary initial states and in such that U 0,1 ≠ U 0,2 and V 0,1 ≠ V 0,2. These two states converge, according to Eq. (3.15), to and respectively, with . Given that function z(R, K) is monotone (injective) in (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 they necessarily converge to different points. Therefore neither single states nor subsets of are attractive in . So, is the smallest attractive set and the proof is concluded. □

Note that and are in the closure of the open set which is not in . In other words, Theorem 3.1 shows that any initial state in converges to a point onto the boundary of that does not belong to . Furthermore note that, an initial state of the form (U 0, 0, 0), (i.e., a state in ) cannot be attracted by any set since it is an equilibrium state (every state in will remain unmodified). This is the reason why it is not possible to consider the attractivity of in .

Local stability of

The next theorem shows the formal Lyapunov (or ) stability of the equilibrium set .

Consider system (2.1) constrained by the positive set . Then, the equilibrium set defined in (2.8) is locally stable.

Let us consider a particular equilbrium point x ≔ (U, 0, 0), with (i.e., ). Then a Lyapunov function candidate is given by (similar to one used in Nangue (2019) for chronic infections)

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

Now, given with it follows that for every (note that it is not true that for x ≠ x, as shown next, in Remark 3). Then, J is a Lyapunov function for system (2.1), which means that each is stable (see Theorem A.1 in 1). Therefore, it is easy to see that the equilibrium set is also stable, which completes the proof. □

Note that, in the latter proof, it is not true that for every nonegative x ≠ x. If for instance, the function is evaluated at with we have that . In fact, is null along the whole U axis, given that this axis is an equilibrium set. This means that the (individual) states in are stable, but not attractive.

A schematic plot of such a behavior can be seen in Fig. 3 .

Fig. 3

Every point in is stable but not attractive. Initial states x0 starting arbitrarily close to x remain (for all t ≥ 0) arbitrarily close to x, but do not converge to x. As a consequence, set is AS but the points inside it are not.

A similar behavior can be seen in system when or the 2-state Kermack–McKendrick epidemic model (Brauer, 2005, Brauer, Castillo-Chavez, Castillo-Chavez, 2012): being S the susceptible and I the infected individuals. In this latter model, and the critical value for S is . The AS set is given by all the states of the form x ≔ (S, 0), with S ∈ [0, S). Furthermore, for this system, the maximum of I occurs when .

Asymptotic stability of

In the next Theorem, based on the previous results concerning the attractivity and stability of the asymptotic stability is formally stated.

Consider system (2.1) constrained by the positive set . Then, the set defined in (2.8) is smallest asymptotically stable (AS) equilibrium set, with a domain of attraction given by .

The proof follows from Theorems 3.1, which states that is the smallest attractive in and 3.2, which states the local stability of . □

A critical consequence of the latter Theorem is that no equilibrium point in (neither in nor in ) 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) (Hernandez-Vargas, 2019). So, the previous stability analysis is still valid for such controlled systems, since only a modification of some of the parameters defining is done. In such a context, only a controller able to consider the whole set as a target (a set-based control strategy, as zone MPC (Ferramosca, Limon, González, Odloak, Camacho, 2010, González, Rivadeneira, Ferramosca, Magdelaine, Moog, 2020)) will be fully successful in controlling system (2.1). Further details concerning antiviral treatments are given next, in Section 4.1.

Characterization for different initial conditions

In this section some further properties of system (2.1) concerning its dynamic are stated, based on the initial conditions at the infection time . 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 .

Consider system (2.1) , constrained by the positive set at the beginning of the infection, i.e., and (i.e., ). Consider also that V 0 is small enough to describe the beginning of the infection. Then,

U ∞ → 0 when U 0 → ∞ or U 0 → 0.

when .

for initial conditions .

for initial conditions .

If and V 0 ≈ 0 then . Therefore and by (3.13).

when which means that either or . This implies that U 0 → 0 or U 0 → ∞, respectively.

when which is true if or, the same, when .

is strictly decreasing for (note that and are in (0,1), since they are smaller than ), while is strictly decreasing in . So, which implies that .

is strictly increasing for while is strictly decreasing in . So, which implies that . Fig. 4 shows U ∞ as a function of U 0, taking V 0 as a parameter. □

Fig. 4

According to Eq. (3.13), U∞(U0) is plotted for different values of V0. All parameters are equal to 1 for simplicity, which means that .

Virus behavior from the infection time

Consider system (2.1) , constrained by the positive set at the beginning of the infection, i.e., and . If the virus spreads (according to Definition 1 ), then for some α(0) > 0 (or, the same, ) and there exist positive times t such that where and are the times at which V(t) reaches a local minimum and a local maximum, respectively, is the time at which I(t) reaches a local maximum, and t(t) reaches . Furthermore, for all .

First, note that since the initial conditions are and V(0) > 0. Even more, assuming the virus spreads, which means that V(t) reaches a local maximum at some time . Therefore, V(t) must reach a local minimum at some .

Now, by Lemma 1 in Appendix 3, it is and respectively, and it is easy to see that is a decreasing function, so it follows that . Then there exists α(0) > 0 such that and, besides, .

From the minimum and maximum conditions of V, at times and we have and respectively. After some algebraic computation, it is easy to see that and which means that I(t) must reach a maximum at some time fulfilling . Moreover, it must be

Given that for (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 and, then, . Therefore, which concludes the proof. □

The value of α(0) is necessary to properly understand and characterize the system behavior according to the initial conditions and parameters. In epidemiological models (SIR, etc.), where is a necessary and sufficient condition for the disease to spread in a population, in our case 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 (or ). See Fig. 6, lower plot, for an example. The value of α(0) can be computed numerically and it is usually small in comparison with (for all the patients simulated in Section 5, ).

Fig. 6

Time evolution of U, I and V, with unitary parameters β, δ, p, c, for initial conditions (upper plot) and (lower plot).

To clarify the results of this section, Figs. 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 . The first plot (Fig. 5) depicts how every state trajectory - even those starting close to - converges to . As stated in Property 2, U ∞ approaches from below, as U(0) approaches from above. Also it can be seen how the virus load starts to decrease only once U(t) is smaller than as stated in Theorem 4.1. On the other hand, the second plot (Fig. 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 while in the lower plot, the initial conditions produce . As it can be seen, only in the first case the virus spread in the host body (i.e., for some t > 0), as stated in Theorem 4.1.

Fig. 5

Phase portrait of system (2.1), with unitary parameters. Empty circles represent the initial states, while solid circles represent final states. Note that only the initial states with corresponds to scenarios with .

Remarks concerning antivirals treatments

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 (the analysis for β is almost the same), being η(t) ∈ (0, 1) the effectiveness of the antiviral and t the time of treatment initiation. The antiviral pharmacokynetics (PK) is not considered for simplicity, which means that the antivirals instantaneously modify η at time t. Then, as the virus monotonically goes to zero only once U(t) is below U, 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) is such that (i.e., such that ). This condition defines a threshold for the antiviral effectiveness (say, a minimal critical value η(t)) 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) < η(t), for t > t), may be counterproductive. Indeed, to slow down V(t) by decreasing p or β, implies that increases, but also soften the decreasing behavior of U(t). As a result, the time t (and so, the virus peak time ) 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.

Characterization of the SARS-CoV-2 target cell model

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 Wölfel et al. (2020). The initial number of target cells U 0 is estimated as approximately 107 cells (Hernandez-Vargas & Velasco-Hernandez, 2020). 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, Figs. 7 and 8 ). The parameters and the initial conditions (U 0, I 0 and V 0, with the infection time) of each patient are collected in Table 1 .

Fig. 7

SARS-CoV-2 Dynamics. The continuous blue line is the simulation with parameter values presented in Hernandez-Vargas & Velasco-Hernandez, 2020. The patient labeling is as presented in Wölfel et al. (2020). V denotes a value of 50 [copies/ml] under which the virus is not detectable.

Fig. 8

Susceptible cells dynamics. The continuous blue line is the simulation with parameter values presented in Hernandez-Vargas & Velasco-Hernandez, 2020. The patient labeling is as presented in Wölfel et al. (2020). Simulation for the patient C shows a very low value of U∞ (practically zero), which suggests that the selected value of may be large.

Table 1

Target limited cell model parameter values for different patients with COVID-19 (Hernandez-Vargas & Velasco-Hernandez, 2020).

Patient	β	δ	p	c
A	9.98×10^−8	0.61	9.3	2.3
B	1.77×10^−7	14.11	20.2	0.8
C	8.89×10^−7	79.51	134.4	0.4
D	3.15×10^−8	45.51	620.2	2.0
E	5.61×10^−8	7.51	96.4	5.0
F	1.41×10^−8	37.61	995.0	0.6
G	1.77×10^−8	8.21	338.4	5.0
H	1.58×10^−8	21.11	927.8	1.8
I	4.46×10^−9	4.21	994.6	4.3

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 for all the patients, so it is not taken into account for the study.

Table 2

Characterization Parameters of patients with COVID-19.

Patient	U_c	U∞	R_0	K_0	t^_I	tc	t^_V	Vmax
A	1.51 × 106	1.36 × 104	6.61	−2.17×10^−7	10.16	10.24	10.58	1.73 × 107
B	3.15 × 106	4.88 × 105	3.18	−6.87×10^−8	11.54	12.26	12.32	4.35 × 106
C	2.66 × 105	4.81×10^−10	37.57	−6.89×10^−7	1.43	1.67	1.69	1.47 × 107
D	4.65 × 106	1.67 × 106	2.15	−4.89×10^−9	9.04	9.42	9.44	2.33 × 107
E	6.94 × 106	4.58 × 106	1.44	−3.48×10^−9	15.02	15.16	15.24	4.03 × 106
F	1.61 × 106	2.03 × 104	6.21	−7.28×10^−9	7.12	7.76	7.78	1.42 × 108
G	6.84 × 106	4.43 × 106	1.46	−1.1×10^−9	14.80	14.92	15.00	1.44 × 107
H	2.59 × 106	2.3 × 105	3.86	−2.72×10^−9	5.16	5.44	5.48	1.577 × 108
I	4.08 × 106	1.14 × 106	2.45	−3.21×10^−10	9.28	9.38	9.50	2.60 × 108

Figs. 7 and 8 show the dynamics of V and U. As expected, the states converge to although significantly different behaviors can be observed for the different patients. From Fig. 8 it can be seen that the healthy cells final value U ∞ is reduced in cases of patients with large values of in spite all simulations have the same initial U 0. This can be explained from the fact that is monotonically decreasing for (see Figs. 1 and 2), and therefore, for R 01 > R 02 > 1 (see Property 2, above). Note that the susceptible cells of patient C converges to U ∞ equals to which can be explained by the fact that this patient has a reproduction number () of 37.57, which is 5.2 times above the cohort mean value of 7.21. Fig. 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 Fig. 7, it can be seen that for all the cases the viral load spreads (i.e.: the virus presents a peak) although for all patients (i.e., ). This can be justified since and, therefore, will be greater than for all patients (note that, ). Moreover, from Table 2, we can corroborate that 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) (Diao, Wang, Tan, Chen, Liu, Ning, Chen, Li, Liu, Wang, et al., 2020, Tay, Poh, Rénia, MacAry, Ng, 2020), 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.

Conclusions

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.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.