Mathematical analysis of a within-host model of SARS-CoV-2

In this paper, we have mathematically analyzed a within-host model of SARS-CoV-2 which is used by Li et al. in the paper “The within-host viral kinetics of SARS-CoV-2” published in (Math. Biosci. Eng. 17(4):2853–2861, 2020). Important properties of the model, like nonnegativity of solutions and their boundedness, are established. Also, we have calculated the basic reproduction number which is an important parameter in the infection models. From stability analysis of the model, it is found that stability of the biologically feasible steady states are determined by the basic reproduction number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$(\chi _{0})$\end{document}(χ0). Numerical simulations are done in order to substantiate analytical results. A biological implication from this study is that a COVID-19 patient with less than one basic reproduction ratio can automatically recover from the infection.


Introduction
A deadly virus, called SARS-CoV-2 and mostly known as coronavirus, has affected most of the countries from November 2019. The World Health Organization declared it as a pandemic, as the virus affected more than a hundred countries and killed lots of people worldwide, especially in Europe and the USA in a short period [1]. During this rampant situation, to stop the virus propagation, almost all the regions of the world have been locked-down by applying strict social distancing measures and banned social gathering in all circumstances. Research is going on all over the world to control this epidemic from different points of view, such as pathology [2,3], microbiology [4,5], mathematics [6,7], and so on [8][9][10]. Until now, several mathematical models have been suggested to study the population dynamics of coronavirus [8], but few of them could predict the outcomes accurately.
As we know, mathematical studies demonstrate the best strategies for controlling any type of disease or virus infection in humans [11,12]. The modeling of any disease helps to decide on public health and introduces the disease treatment. When the primary studies of any deadly disease cannot give us any concluding remarks to take significant decisions to control the pandemic situation, the key facts of mathematical modeling always help the human civilization. Mathematical models of viral dynamics with viral and antiviral drug effects are examined and assisted by their transparent structures and biologically feasible parameters. There is much work to be done now and in the future to assess the figure of confirmed cases of COVID-19. Also, a few studies focused on confirmed cases, risk history, and disease schedule characteristics [8,12,13]. Mathematical models of epidemiology have been developed to help policymakers in making the right decisions at the right time. These models highlight that isolation plays an important role in controlling the spread of the virus [6,13]. But this is not a solution to the long-term vision as a country's economic growth globally is declining due to the locked-downs to maintain social distancing [7]. Therefore, researchers and scientists are now developing dynamic viral within-host models to understand the mechanisms of SARS-CoV-2 in the human body. To date, only a few host-level models [4,5] have been developed to understand the SARS-CoV-2 replication cycle and its interactions with the immune system. Most of them are linked to HIV [14], hepatitis virus [15], Ebola [16], influenza [9,17], and other models.
At present, clinically there is no effective treatment developed to remove the virus from the human body. However, the research is going on. Researchers provided many treatments (like plasma therapy, monoclonal antibody therapy, etc.) which are effective for early diseases like SARS-CoV, MERS-CoV, Ebola, Influenza, HIV-like virus disease. Also we all know that our body immune system gives a good response to fight any viruses or diseases [2].
At the time of SARS-CoV-2 infection, macrophages are first targeted by SARS-CoV-2 and after that these SARS-CoV-2 propagate to T cells. At this stage of the virus route, the T cells activate and further they involve differentiation. Besides, the T cells produce cytokines (INF-α, IL-6, and IL-10) associated with the different types of a T cell. A large amount of cytokines provides a greater activation of the immune response to fight the virus. Particularly T cells, CD4 + T cells, and CD8 + T cells have played a significant antiviral role in a fight against pathogens. There is also a risk of mounting autoimmunity or devastating inflammation. CD4 + T cells help the immune system of the body by generating virus-specific antibodies with the activation of T-dependent B cells. However, CD8 + T cells can kill virally infected cells, as they are cytotoxic. In general, a large number of CD8 + T cells in the infected SARS-CoV-2 body are found in nearly 80% of the total infiltrative inflammatory cells in the interstitial pulmonary tract, which play a significant role in clearing CoVs. The loss of CD4 + T cells is correlated with reduced conscription of lymphocytes and neutralizing the production of antibodies and cytokines, resulting in severe immune-mediated interstitial pneumonitis and delayed SARS-CoV lung clearance [2,18].
Researchers have shown that there is a long-lasting and persistent response of T cells to the S and other structural proteins (including the proteins M and N), which provides sufficient knowledge to draft SARS vaccine by combining viral structural proteins. These types of vaccine may provide a strong, efficient, and long-term response to the virus by memory cells [3]. Also, clinical trials show that a monoclonal antibody therapy is an effective treatment tool which responds better to SARS-CoV-2 [10]. In our paper, we studied a model of in-host viral kinetics, that describes the response of SARS-CoV-2 to epithelial cells. In Sect. 2, we describe a previously proposed model by Li et al. [5]. We have shown in Sect. 3 that the solutions of the model are biologically feasible. The basic reproduction number is computed and steady state analysis is done in Sect. 4. Sections 5 and 6 deal with the local and global stability of the model, respectively. In Sect. 7, we fitted our theoretical results with MATLAB. A concluding discussion is presented in Sect. 8.

The SARS-CoV-infection model
To study the within-host dynamics of SARS-CoV-2 infection, we consider the mathematical model used by Li et al. [5], which is represented by the following ordinary differential equations system: (2.1) with the initial conditions as below: Here, the model consists of three population compartments: virus-free pulmonary epithelial cells denoted by E p (t), virus infected pulmonary epithelial cells denoted by E * p (t), and the free virus denoted by v(t). Also, d E , d E * , and d v represent the death rates of virus-free epithelial cells, virus infected epithelial cells, and the free virus, respectively. Also β denotes the rate at which the virus-free epithelial cells are infected by free virus, π v is the production rate of free virus. In this model, the term d E E p (0) describes the constant regeneration rate of the virus-free pulmonary epithelial cells.

Basic properties of model (2.1)
Lemma 1 All analytic solutions of model (2.1) with initial conditions (2.2) are nonnegative for all t > 0.
Proof We can write the first equation of model (2.1) as Therefore, which implies From the above equation, it is clear that E p (t) > 0 for t > 0. It can be shown in the same way that E * p (t) ≥ 0, v(t) ≥ 0 for t > 0 [19].
Proof Adding the first two equations of the model (2.1), we get where δ =min{d E , d E * }. Therefore, the numbers of both the virus-free and infected pulmonary epithelial cells are always bounded. From the third equation of the model, we can easily show that the population of the free virus v is also bounded above. We consider K ≥ 0 as the maximum of the above bounds. Therefore, we get the bounded set where K ≥ 0. It is clear that this set is convex. Also using the nonnegativity criteria of solutions of the model (2.1) from Lemma 1 and the boundedness conditions from Lemma 2, it is clear that with the initial conditions (E p (0), E * p (0), v(0)) ∈ , we have solutions of the model (2.1) again in , i.e., (E p (t), E * p (t), v(t)) ∈ for all t ≥ 0. Hence, the set is positively invariant w.r.t. the model (2.1).

Basic reproduction number and steady states
It is straightforward to show that the system has a noninfected steady state S 1 (E p (0), 0, 0). Now, the basic reproduction ratio of the model (2.1) will be calculated with the help of the next generation matrix method. Considering Y = (E * p (t), v(t), E p (t)), we can write the system (2.1) as Here A(Y ) stands for the rate of new infections, T(Y ) stands for the transfer rate of individuals; they are given by The jacobian matrices DA(E p (0)) and DT(E p (0)) at the noninfected steady state are Therefore, the basic reproduction number (χ 0 ) of the system (2.1) is given by the spectral radius of the matrix at -1 and hence χ 0 = Also, from a simple calculation it is clear that the model (2.1) has a unique infected steady state

Local stability
In this section, we study the local stability behaviors at the different steady states of the model.
Theorem 1 If χ 0 < 1, then the noninfected steady state S 1 is locally asymptotically stable in while for χ 0 > 1 it becomes unstable.
Proof For the noninfected steady state S 1 , the Jacobian matrix J(S 1 ) is given by The characteristic equation of the matrix J(S 1 ) is where Therefore, if χ 0 < 1, then all the conditions of the Routh-Hurwitz criterion are satisfied. Hence, the noninfected steady state S 1 becomes locally asymptotically stable if χ 0 < 1 and unstable if χ 0 > 1.

Theorem 2 The virus-infected steady state S 2 is locally asymptotically stable whenever
Proof For the infected steady state S 2 , the Jacobian matrix J(S 2 ) is given by The characteristic equation of the above matrix becomes This implies that all conditions of the Routh-Hurwitz criterion are satisfied. Hence, the infected steady state S 2 is locally asymptotically stable whenever it exists.

Global stability
From local stability analysis, the behaviors of the system (2.1) near the steady states are obtained. For the behavior of the system (2.1) far away from the steady states, we have carried out a global stability analysis in this section.

Theorem 3
If χ 0 ≤ 1, then the noninfected steady state S 1 approaches a globally asymptotically stable state in and it becomes unstable if χ 0 > 1.
Proof Define Lyapunov functional of the model (2.1) as It is clear from (6.2) that when χ 0 ≤ 1, dL dt ≤ 0. Define S to represent the set of solutions of model (2.1) where dL dt = 0. We have dL dt = 0 when v = 0 or χ 0 = 1 and E p ≤ E p (0). Using Lyapunov-Lasalle theorem [20], we have all curves in approach the set S which is also positively invariant. Again, on the boundary of where v = 0, we have E * p = 0 and . Hence, E p → E p (0) as t → ∞. Therefore, when χ 0 ≤ 1, all solution curves in the domain go to the virus-free steady state S 1 and hence S 1 is globally asymptotically stable.
From J(S 1 ), it can be easily seen that for χ 0 > 1, one root of the characteristic equation will be positive. Hence, the noninfected steady state S 1 is unstable when χ 0 > 1. Now, the global behavior of the virus-infected steady state S 2 whenever it exists will be studied by applying Li and Muldowney criterion [21]. It is easy to visualize the simple connectedness of the interior of set . Also, for χ 0 > 1, the system (2.1) has the unique steady state S 2 in int( ). Theorem 3 verifies that H ⊂ exist where H is an absorbing compact set for model (2.1). Hence, all the assumptions for Li and Muldowney global stability criterion [21] are satisfied, and we have the following theorem.

Theorem 4
The infected steady state S 2 is globally asymptotically stable whenever R 0 > 1.
Proof The Jacobian matrix of the model (2.1) is The associated second compound matrix J [2] [21,22] is given by to compute the matrix W f , each entry w ij of matrix W has to be replaced by its derivative in the direction of f . Then, , and Now, define Lozinskii measure as follows: where g 1 = μ(B 11 ) + |B 12 | and g 2 = |B 21 | + μ(B 22 ). Here, |B 12 | and |B 21 | are with respect to vector norm. Now, Using model equations, Hence, Therefore, we get , v(t)) be an arbitrary solution such that (E p (0), E * p (0), v(0)) ∈ H and consider any solution (E p (t), E * p (t), v(t)) ∈ H for all t ≥t wheret is a sufficiently large time. Then, for t >t, Consequently,

Numerical simulation
This section is devoted to numerical simulations in order to substantiate the analytic results. Using Matlab software, analytic results are fitted with the parameters from biologically feasible range for SARS-CoV-2. Based on chest radiograph score data, Li et al. [5] have estimated the parameters of the model (2.1) using the Monte Carlo Markov Chain (MCMC) method. Using that parameter set, more specifically, E p (0) = 22.41, E p * (0) = 2.59, v(0) = 0.061, π v = 0.24, β = 0.55, d E = 10 -3 , d E * = 0.11, and d v = 5.36, we have χ 0 = 5.01716 > 1. Theorem 4 implies that the infected equilibrium point S 2 (4.467, 0.163, 0.007) is globally asymptotically stable. Figure 1(a) verifies this analytic result. Also, Fig. 1(b) describes the dynamics the numbers of virus-free epithelial cells, virus infected epithelial cells, and free virus.
In order to verify the analytic results of the noninfected steady point S 1 , we consider lowering the infection rate from β to 0.1β by implementing a drug [5]. Hence with this assumption, considering the other parameters the same as in the above simulation, we get χ 0 = 0.501716 < 1. As stated by Theorem 3, the noninfected steady state S 1 (22.41, 0, 0) is globally asymptotically stable. This analytic result can be verified by Fig. 2(a). Also, the dynamics of the numbers of virus-free epithelial cells, virus-infected epithelial cells, and virus in this scenario are depicted by Fig. 2(b).
In the above simulations, we have considered that infection rate β can be lowered by implementing a drug, which leads to the removal of the virus from the human body. But with greater infection rate and other parameters in such a way that χ > 1, the viral load cannot be controlled and infection persists in the human body. This situation can be handled with a proper treatment strategy, otherwise the number of the virus-free epithelial cells will decrease and those of the virus-infected epithelial cells and the free virus will increase in a patient's body with time. As a result, the patient's condition will deteriorate with time, which will lead to extreme conditions like death of the patient.

Discussion
In this paper, we have considered a basic model for within-host dynamics of SARS-CoV-2 used by Li et al. [5] and mathematically analyzed it. First of all, we have proved that all analytic solutions of model (2.1) are nonnegative and uniformly bounded, conditions which are necessary for a biologically feasible model. For the model, we have found two biologically feasible steady states, noninfected and infected steady states. Local stability of both steady states was discussed. Also, we have found the global stability criteria for both steady states. Biologically, it follows that for the basic reproduction number χ 0 < 1, infection will be cleared from a human body without any treatment; otherwise we have to implement some treatment in order to reduce and to remove the infection from the body. It is our future work to apply different treatment policies in this model in order to clear the virus from an infected human body. It is also found that if basic reproduction number is greater than one then infection will persist for any amount of viral load in the host's body.