Lyapunov function and global asymptotic stability for a new multiscale viral dynamics model incorporating the immune system response: Implemented upon HCV

In this paper, a new mathematical model is formulated that describes the interaction between uninfected cells, infected cells, viruses, intracellular viral RNA, Cytotoxic T-lymphocytes (CTLs), and antibodies. Hence, the model contains certain biological relations that are thought to be key factors driving this interaction which allow us to obtain precise logical conclusions. Therefore, it improves our perception, that would otherwise not be possible, to comprehend the pathogenesis, to interpret clinical data, to control treatment, and to suggest new relations. This model can be used to study viral dynamics in patients for a wide range of infectious diseases like HIV, HPV, HBV, HCV, and Covid-19. Though, analysis of a new multiscale HCV model incorporating the immune system response is considered in detail, the analysis and results can be applied for all other viruses. The model utilizes a transformed multiscale model in the form of ordinary differential equations (ODE) and incorporates into it the interaction of the immune system. The role of CTLs and the role of antibody responses are investigated. The positivity of the solutions is proven, the basic reproduction number is obtained, and the equilibrium points are specified. The stability at the equilibrium points is analyzed based on the Lyapunov invariance principle. By using appropriate Lyapunov functions, the uninfected equilibrium point is proven to be globally asymptotically stable when the reproduction number is less than one and unstable otherwise. Global stability of the infected equilibrium points is considered, and it has been found that each equilibrium point has a specific domain of stability. Stability regions could be overlapped and a bistable equilibria could be found, which means the coexistence of two stable equilibrium points. Hence, the solution converges to one of them depending on the initial conditions.


Introduction
Hepatitis C virus (HCV) is a bloodborne virus that has become one of the most serious infectious diseases that threaten human health [1]. It can lead to both acute and chronic hepatitis, ranging in severity from a mild to a serious lifelong illness. A significant number of those who are chronically infected will develop cirrhosis or liver cancer, where HCV is a major cause of liver cancer. An estimated 71 million people have chronic HCV infections worldwide. The world health organization (WHO) estimated that in 2016, approximately 399 000 people died from hepatitis C, mostly from cirrhosis and liver cancer [2]. Treatment for chronic hepatitis C infection began in the early 1990s with interferon-alfa [3]. This injectable drug worked by improving the immune system, rather than by specifically attacking the virus. In 1998, the oral drug ribavirin was added to interferon [4]. The development of the treatment occurred in 2002 with the approval of pegylated interferon-alfa, a process that makes interferon more durable and effective [5]. New treatment options include direct-acting antiviral agents (DAAs) targeting specific HCV-life cycle components [6]. The quick steps of HCV drug development where a cure rate of more than 95% was achieved [2], have led to the hopeful prediction that full eradication of HCV is theoretically possible in the absence of a vaccine for HCV. There remain many barriers that need to be overcome. Further studies for factors that increase the eradication rate are needed. Such barriers include factors related to awareness, linkage to care, the development and availability of simplified and highly effective drug regimens, improving the rates of detection of infection, and the availability of funds expertise [7,8].
Mathematical modeling is a useful tool to study and analyze many engineering and physical problems. It has also been used to describe some biological processes such as heartbeats [9], tumor growth and cancer treatment [10,11], and virus dynamics for many types of viruses such as HCV, HBV, HIV, and Covid-19 [12][13][14]. It provides a powerful tool in the study of virus dynamics because it helps to understand the biological mechanisms and interpret the experimental results. Mathematical models can be used to predict the virus behavior under certain conditions or to determine which parameters increase disease spread. They can also be used to predict the number of medications required to help eradicate a disease or control it [9]. Mathematical modeling is also useful in public health policy formulation addressing the control of infectious diseases [15]. The early mathematical model for HCV was developed and analyzed in [16,17] as a scheme consisting of a system of ordinary differential equations describing the basic dynamics of the hepatitis C virus in-vivo. Models for HCV treatment with DAAs therapy is considered in [18][19][20][21][22]. A novel approximate analytical solution for solving the standard viral dynamic model for HCV is presented in [23]. Local and global stability analysis for basic virus dynamics models is studied in [24][25][26].
Interactions between replicating virus, liver cells, and different types of immune responses (CTLs and antibodies) are highly complex and nonlinear, so these interactions between HCV and the immune system were studied with mathematical models in [27,28] and stability has been analyzed for these models. However, these models can only describe the intercellular viral dynamics and cannot describe the intracellular viral dynamics which is required to capture the different antiviral effects corresponding to the action mechanisms of drugs.
Authors in [29][30][31][32] constructed a multiscale model that accounts for the dynamics of intracellular viral replication, and which includes the major stages in the HCV life cycle that are targeted by DAAs. These multiscale models have been developed using partial differential equations (PDEs). Since numerical PDE solvers are time-consuming and often converge poorly, a new approach has been suggested by Kitagawa,et al. [33] that converts a standard PDE multiscale model of the HCV infection into an equivalent system of ordinary differential equations (ODEs) without any assumptions. This transformed model prevents time-consuming calculations and has become widely available for further mathematical analysis. Kitagawa, et al. [34] derived the basic reproduction number of the transformed ODE model and studied the global stability of the model using Lyapunov-LaSalle's invariance principle and investigated all possible steady states of the model. Local stability analysis of this model is considered in [35] using Routh-Hurwitz criterion. In that work, sensitivity analysis had been performed to specify the influence of each parameter on the basic reproduction number.
It is worth mentioning that the classical multiscale model ignores the responses of the immune system, which have a significant role in reducing viral load. From the point of view of the mathematical analysis, considering the immune system with a multiscale model in the PDE form is an undesirable task. In this paper, a new mathematical model, that deals with the interaction between the transformed multiscale model of the HCV infection in ODE form and immune responses, is proposed. Different antiviral effects of multidrug treatments are presented by defining three efficacies which are responsible for blocking intracellular viral production, blocking virion assembly/secretion, and enhancing the degradation rate of vRNA. The model contributes to improving our realization to the interactions between HCV, drug treatments, infected cells, and immune system. For instance, the analysis of the model reveals the existence of five equilibrium points: an uninfected point, an infected point with no immune responses, an infected point with dominant antibody responses without CTLs, an infected point with dominant CTLs responses without antibody, and an infected point with coexistence responses of both CTLs and antibody. In section 2, the ODE model extracted from the multiscale model for describing the dynamics of the HCV infection is described. Consequently, this model is extended to consider the impact of the immune system response. The proof of the positivity of the extended model and calculation of the reproduction number for the model is presented in section 3. The equilibrium points are determined in section 4 and the stability analysis is presented in section 5. Finally, section 6 concludes the paper.

Extended model for the transformed multiscale ODE HCV model
A multiscale model in the form of PDEs, which describes the intracellular life cycle had been proposed and applied by many researchers [29][30][31][32] for analyzing clinical data under multidrug treatment. The model is as follows: @iða; tÞ @t þ @iða; tÞ @a ¼ À d i a; t ð Þ ð3Þ The variables T(t) and V(t) are the numbers of target cells and viruses, respectively, and the variable i(a,t) represents the age distribution of infected cells. Similarly, R(a,t) is the age and time distribution of intracellular viral RNA (vRNA) in a cell with infection age a. The target cells are assumed to be produced at rate s, infected by viruses at rate β, and naturally die at rate r. The infected cells die at rate δ, and virions are cleared at rate c. The parameters α and μ denote the production and degradation rates of the intracellular viral RNA, respectively. Viral RNA is assumed to assemble along with viral proteins and to be secreted from an infected cell as viral particles at rate ρ. The model recognizes the different antiviral effects of multidrug treatments by defining three efficacies ε α , ε s , and k�1, which are responsible for the actions of blocking intracellular viral production, blocking virion assembly and/or secretion, and increasing the degradation rate of vRNA, respectively.
Kitagawa et al. [33] transformed the previous multiscale PDE model into the following ODEs model: where I(t) denotes the total number of infected cells and is defined as IðtÞ ¼ R 1 0 iða; tÞda, and P(t) is the total amount of intracellular viral RNA pooled in all infected cells and defined as PðtÞ ¼ R 1 0 Rða; tÞiða; tÞda. The entry virus-derived RNA starts to replicate from z copies in a newly infected cell and is fixed to 1 [33,34].
In this work, an extension to the transformed multiscale ODE model is proposed. Two more variables are added to stand for the immune system response. The first one represents the CTLs number which is responsible for killing the infected cells accordingly inhibiting the reproduction of the virus and is denoted by Z(t). The second one represents the number of the antibodies generated which is responsible for neutralizing the virus in-vivo and is denoted by W(t). Hence, the proposed model is described by the following ODEs system: The term f I(t)Z(t) in Eq (10) represents the rate of killing the infected cells by the CTL response and the term q V(t)W(t) in Eq (12) represents the rate of neutralizing virus particles by the antibodies. CTLs become activated in response to viral antigen derived from infected cells, and once activated, they are divided, and their population grows (clonal expansion). So, in Eq (13), the CTLs increase at a rate of uI(t)Z(t). The CTLs decay at a rate of bZ(t) due to the lack of antigenic stimulation. Antibodies are produced by B cells and initially they are attached to them. They serve as the receptor that can specifically recognize the virus. When the B cells are exposed to a free virus, they divide and secrete the antibodies. Accordingly, antibodies progress at a rate gV(t)W(t) and decay at a rate hW(t) in Eq (14).
To avoid complexity in the mathematical analysis, the saturation effects of the concentrations of all variables are not contained in the proposed model. Yet, since the proliferation terms are not limited by saturation, the model could predict unlimited increases in the values of these variables which are certainly unrealistic. In general, however, the model can describe the dynamics of these variables and to gain important insights, as long as one is aware of the model limitations, and the results obtained do not depend on the unrealistic values of variables.
The proposed model presented by Eqs (9)- (14) can be used for numerical simulations to represent a variety of medical cases under treatment and can be a valuable tool to comprehend the pathogenesis and in controlling treatment of chronic HCV. Yet, the computation of the basic reproduction number, the determination of the equilibrium points, and stability analysis for this model have to be accomplished under no treatment. No treatment can be specified by assigning ε α = 0, ε s = 0, and k = 1. For clarity and effectiveness, we demonstrate the form of the proposed model under no treatment as:

Non-negativity of the solutions
To retain the biological fidelity of the model, the solutions to the mathematical model have to be non-negative. Proof. We know that all of the parameters used in the system are positive. Thus, we can place lower bounds on each of the Eqs (9)- (14). Thus, Through basic differential equations methods, we can resolve the inequalities and produce: PðtÞ � e À ðkmþð1À ε s ÞrþdÞt � 0; WðtÞ � e À h t � 0; Thus, for all t2[0,τ], T(t), I(t), P(t), V(t), Z(t) and W(t) will remain non-negative in R 6 .

Computation of the basic reproduction number (R 0 )
The basic reproduction number is defined as the expected total number of viral particles newly produced during the whole period of infection from one typical viral particle in a population consisting only of uninfected cells. Accordingly, the basic reproduction number R 0 is calculated under no treatment condition described by Eqs (15)- (20), and it is also computed at disease-free equilibrium E 0 . E 0 is the uninfected equilibrium point, which will be explained in Many methods can be used to obtain the basic reproduction number, see for example [36]. The chosen method is the next-generation method, which was introduced Diekmann et. al., [37]. There are two principal approaches to apply this method elaborated by Driessche and Watmough [38] and by Castillo-Chavez, et. al., [39]. In this work, the second approach is considered, and an outline of this approach is given, proofs and further details can be found in [36,37,39]. Variables T(t), I(t), P(t), V(t), Z(t) and W(t) can be discretized into three groups: the non-infected group ϕ, the infected but not infectious group ψ, and the infected and infectious group γ. Hence, we have ϕ = (T,Z,W), ψ = (I,P), and γ = (V). The model in Eqs (15)- (20) can be written as The uninfected equilibrium point is given by Solving these two equations for I and P in terms of V, gives: Substituting in h (ϕ 0 ,I,P,V) leads to: Let hence: G can be written as: Hence, the basic reproductive number is the spectral radius: that means that:

The equilibrium points
Equilibrium points are the values of the variables T�, I�, P�, V�, Z� and W�, under no treatment, at which the derivatives of these variables, i.e., the left-hand sides in Eqs (15)-(20), vanish. These equilibrium points represent the steady states after the cease of medication. In Fact, in stability analysis the interest is in specifying the behavior of the virus after the cease of medication. Hence, these equilibrium points satisfy the following algebraic equations: The commercial program Mathematica 12 program is used solve these algebraic equations to obtain the equilibrium points. The program gives six equilibrium points, however, one of them has negative coordinates that have no biological meaning. The five other points are: where ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi The first point, E 0 , is a virus-free equilibrium point, while the other four points are virusinfected. These four infected equilibrium points are: an infected state with no immune responses, an infected state with dominant antibody responses without CTLs, an infected state with dominant CTL responses without antibodies, and an infected state with coexistence responses of both CTLs and antibodies, respectively.
Since the equilibrium points should have non-negative coordinates only, the following conditions for existence can be obtained: It is also required that 2 c r u μ 1 �k 3 −k 1 , however, this is always true and would not add a new condition. When R 0 = A 2 then E 3 = E 1 .

E 4 exists only if
These conditions are necessary and sufficient conditions for the existence and non-existence of the equilibrium points.

Global stability analysis
Usually, the global stability analysis of a dynamical system is a very complex problem. One of the most efficient methods to solve this problem is Lyapunov's theory. To build the Lyapunov function, the technique used in [24][25][26]40], which had been suggested and utilized for other models, is adopted. In this section, the global asymptotic stability of the model for both the uninfected and the infected equilibrium points is investigated.
Assume the following general form of Lyapunov function l(t): where T� means the value of T at the equilibrium point, and whenever T� is zero the corresponding ln term does not exist, and this is applied to all other variables. � 1 , � 2 , � 3 , � 4 , � 5 , and � 6 are constants and will be specified through the proof of the stability of equilibrium points.

Basic properties of Lyapunov function l(t)
The following properties should be demonstrated in any Lyapunov function l(t): 1. It is a continuously differentiable function defined in domain D�R 6 ; 02D, and defined for all T(t)�0, I(t)�0, P(t)�0, V(t)�0, Z(t)�0 and W(t)�0. This property is already satisfied in the proposed function in Eq (45).
2. It is always nonnegative function in R 6 , but equal to 0 at the equilibrium points. In the following subsections, the conditions required for the equilibrium points to fulfill this property are obtained.
3. It satisfies the radial unboundedness condition i.e., if any dependent variable tends to infinity, l(t) also tends to infinity. To show this, let us assume that T(t) in Eq (45) tends to infinity, then: Using L'Hôpital's rule, lnðT=T � Þ T=T � tends to zero, hence, l(t) tends to infinity.

Global stability of uninfected equilibrium point
Theorem 5.1. E 0 is globally asymptotically stable if R 0 �1. Proof: Consider the following Lyapunov function for the uninfected equilibrium point E 0 (T 0 ,0,0,0,0,0): The time derivative of l (t) is: Eq (45) can be simplified by substituting s = r T 0 , which leads to: Since the arithmetical mean is greater than the geometrical mean, then 2 À T Consequently, dlðtÞ dt < 0 for any coordinate values (T,I,P,V,Z,W) and dlðtÞ dt ¼ 0 at the coordinates of the uninfected equilibrium point (T 0 ,0,0,0,0,0). Therefore, R 0 �1 ensures that dlðtÞ dt � 0 which verifies that all the trajectories of the model converge to E 0 , that is, the uninfected equilibrium point E 0 is globally asymptotically stable when R 0 �1. Theorem 5.2

E 3 is globally asymptotically stable if
4. E 4 is globally asymptotically stable if R 0 �max(A 3 , A 4 ).

Proof:
Consider the Lyapunov function for the infected equilibrium points: The time derivative of l(t) is: For the infected equilibrium points, Eq (46) can be simplified by the following substitutions from Eqs (15)- (20): The six terms of Eq (46) can be simplified as shown: The first term: By using Eq (47), the first term can be simplified as shown: The second term: By using Eq (48), the second term can be simplified as shown: The third term: By using Eq (49), the third term can be simplified as shown: The fourth term: By using Eq (50), the fourth term can be simplified as shown: The fifth term: The sixth term: Substituting the sixth terms in Eq (46) gives: Substituting Eq (48) into the previous equation gives Since the arithmetical mean is greater than the geometrical mean, i.e., the terms For the infected equilibrium point E 1 : Substituting Z� = 0 in Eq (48) and then substituting the result combined with Z� = 0 and W� = 0 into Eq (51) gives As per our derivation, dlðtÞ dt ¼ 0 only at the equilibrium point E 1 . Furthermore, for I � � b u and V � � h g it follows that dlðtÞ dt < 0. Substituting the coordinates of the equilibrium point leads to R 0 � u s u sÀ d b and R 0 � 1 þ b h r g . Hence, according to Lyapunov-LaSalle's invariant principle combined with remark 4.1, E 1 exists and it is globally asymptotically stable if 1�R 0 �min(A 1 , For the infected equilibrium point E 2 : Substituting Z� = 0 in Eq (48) and then substituting the result combined with Z� = 0 and V � ¼ g h into Eq (51) gives: It can be noticed that dlðtÞ dt ¼ 0 only at the equilibrium point E 3 . For V � � h g it follows that dlðtÞ dt < 0. Substituting the coordinates of this equilibrium point gives k 1 þk 3 2b c m 1 u � h g , which by simplification leads to R 0 � b s u h ðaþdÞ r d ða b gþb s u h k 2 Þ . Hence, according to Lyapunov-LaSalle's invariant principle combined with remark 4.1, E 3 exists and it is globally asymptotically stable if For the infected equilibrium point E 4 : Substituting Only at the equilibrium point E 4 the derivative dlðtÞ dt ¼ 0 while dlðtÞ dt < 0 at any other point. According to Lyapunov-LaSalle's invariance principle combined with remark 4.1, E 4 exists and it is globally asymptotically stable if R 0 �max (A 3 , A 4 ). This ends the prove of theorem 5.2.

Simulations
In this section, the proposed model and the transformed model are numerically simulated. For each of them, the corresponding system of ODE's is numerically solved using Mathematica 12 program. The parameter values, which were estimated from clinical datasets in [33], are summarized in (Table 1), and the immune system parameters, which were proposed in [41], are given in (Table 2). These parameters are used in the simulations. Some of the values of the parameters are not mentioned in the Tables 1 and 2 and the used values will be mentioned explicitly. To demonstrate the mutual relations between the basic reproduction number, the equilibrium points, and the stability analysis, the form of the proposed model under no treatment represented in Eqs (15)-(20) is simulated first. The simulation demonstrates the variation of all variables T(t), I(t), P(t), V(t), Z(t), W(t) with time. Each figure represents a case with a specific value for the basic reproduction number R 0 which leads to a corresponding equilibrium point that is stable according to theorems 5.1 and 5.2. Figs 1-5 illustrate that the variables converge to that equilibrium point. T 0 , I 0 , P 0 , V 0 , Z 0 , and W 0 are the initial values of the variables.
To demonstrate the effect of medical treatments, the form of the proposed model under treatment represented in Eqs (9)-(14) is considered now. The same cases considered in the first stage is reconsidered in this stage, however, with medical treatment. Figs 6-10 reveal the effect of medication. Whenever the line for the antibodies is not seen it coincides with the line for CTLs. The simulation proves the practicality and the effectiveness of the proposed model.
For the comparison purpose, the simulation for the transformed model is performed. The transformed model has two equilibrium points only. Hence, two cases are considered, one for the noninfected case and the other for infected case. No treatment illustrations are considered by substituting ε α = 0, ε s = 0, and k = 1 in Eqs 5-8 and for the medical treatment illustrations, the values in Table 1 for these parameters are used. Comparing Figs 11-14 with Figs 1, 2, 6 and 7 reveal some important outcomes for including the CTLs and antibodies in the model. The

Discussion
It can be notice from Eq (33) that the parameters related to immune system do not affect the basic reproduction number. To explain that, let us go back to the definition of R 0 which

PLOS ONE
implies that we assume that a virus particle enters the system, when t = 0 and the system is at the disease-free equilibrium E 0 . Accordingly, when t = 0, Eqs (15)- (18) implies that the rates of variables T,I,P,V have nonzero values since V(0) = 1, T(0) = T 0 . Hence, these variables will vary with time. Opposite to that, when t = 0, Eqs (19) and (20) imply that the rates of variables Z and W have zero values since Z(0) = 0, W(0) = 0. Hence, Z and W will not vary with time i.e., Z(t) = 0 and W(t) = 0. Biologically, when one typical viral particle placed in a population  consisting only of uninfected cells T, ultimately, some of these cells would be infected and became infected I cells and consequently P intracellular viral RNA cells. However, since initially there are no CTLs and antibodies, they could not be generated. Therefore, the parameters multiplied by Z and W would not appear in the formula of R 0 . Including the immune system, raised the number of the equilibrium points to five. The first point is a virus-free equilibrium point, and the second point is an infected point with no immune responses. These two equilibrium points are the same as those obtained using transformed model without the immune system in [34]. The last three equilibrium points give an insight about immune response on the stability of the system and notably they could not be Though role of CTL and antibodies for the resolution of HCV infection is debated in the literature, many medical evidence supports the notion of competition between the two branches of the immune system [42,43]. The infected equilibrium point with dominant antibody means that the antibody response is strong and diminishes virus load to a level that is too low to stimulate the CTL response. The infected equilibrium point with dominant CTLs means that the CTL response is strong and diminishes virus load to a level that is too low to stimulate the antibody response. In these two equilibrium points, the competition between the CTLs and

PLOS ONE
antibodies results in exclusion of one of them. The competition could result in co-existence of both as in the fifth equilibrium point.
Comparing the values of the variables T�, I�, P�, and V� for the equilibrium point with no immune responses E 1 and the equilibrium point with dominant antibody responses E 2 , it can be simply proven that T 2 �T 1 , I 2 �I 1 , P 2 �P 1 , and V 2 �V 1 when E 2 exists i.e., R 0 �A 1 . The same is applicable to equilibrium point with dominant CTL E 3 responses and equilibrium point with coexistence responses of both CTLs and antibodies E 4 . Though, the antibody activation and the CTL activation could not eliminate the viral load, they notably increase the uninfected cells, decrease the infected cells and intracellular viral RNA, and reduce the viral load. It is worth to notice that theorems 5.1 and 5.2 prove that each equilibrium point has a specific domain of stability. These domains of stability could be overlapped. For example, if A 1 �A 2 and A 3 �A 4 , the domains of global stability of E 2 and E 3 will be intersecting except if A 3 �A 2 . Hence, a bistable equilibrium could be found, which means the coexistence of two stable equilibrium points. A similar situation had been reported in many biological circumstances, like in multistrain disease dynamics discussed in [36], due to the low capacity for treatment of infective in epidemic models [44], and in investigating bifurcations and stability of an HIV model that incorporates the immune responses [45]. In the presence of bistable equilibria, the solution converges to one of the two stable equilibrium points depending on the initial conditions. Therefore, it is called bistable dominance since the species in the better position originally dominates [36].
Since the proposed model is a multiscale model that incorporate the immune system response, it considers the intracellular viral RNA with the introduction of age-dependency in addition to time-dependency. Hence, the model can explore the dynamics of HCV infection under therapy with DAAs by including both the intracellular viral RNA replication/degradation and the extracellular viral infection with age-dependency in addition to time-dependency. The parameters of the intracellular viral RNA, P, appears in both the basic reproduction number and in the coordinates of the equilibrium points. Therefore, the stability of the proposed model is considerably much more difficult to consider and to analyze compared to the corresponding classical model which could not describe the intracellular viral dynamics [27,28,46].

Conclusions
This work has utilized the transformed multiscale model for HCV in the form of ODE, which is direct and easier to analyze and modify. The immune system, which has a significant role in reducing the virus load, has been incorporated into the multiscale model. The proposed model could represent a valuable tool to comprehend the pathogenesis and controlling treatment of chronic HCV. One of the main advantages of the proposed model over classical multiscale model is its ability to obtain equilibrium points after the cease of medication while in the presence the immune effects. The basic reproduction number of the infection R 0 has a crucial rule in dealing with the stability of the spread of the HCV infection, hence; it has been identified. The parameters related to immune system do not affect the basic reproduction number. The disease-free equilibrium point and the endemic equilibrium points are specified. Conditions for the existence of these points are derived. At any state of the system, only a maximum of five total equilibrium points including the uninfected point can be available. The four infected equilibrium points are: a point with no immune responses, a point with dominant antibody responses without CTLs, a point with dominant CTL responses without antibodies, and a point with coexistence responses of both CTLs and antibodies, respectively. It has been revealed that the four infected equilibrium points are dependent upon the immune system parameters.
Global stability of the equilibrium points has been considered, the Lyapunov principle has been utilized. A new appropriate Lyapunov function has been suggested, hence, sufficient conditions have been derived for the global stability of the five equilibrium points. It has been proven that the uninfected equilibrium point is asymptotically stable if R 0 � 1 and unstable if R 0 > 1. The stability of the four infected equilibrium points depends upon the basic reproduction number and upon the parameters defined by the CTL response and antibody response. Therefore, these parameters play an important role to characterize the stability of the equilibrium points. The activation of antibodies and the activation of CTLs will not eliminate the viral load but they, remarkably, reduce the viral load.
For successful treatment, if R 0 >1 the treatment should be directed to improve the body parameters to ensure that R 0 �1 and then the treatment for reducing the virus could be conducted until the state of the body comes to the attracting zone for the stable uninfected point. Consequently, the immune system will lead the state of the body to a stable uninfected state. Otherwise, if an unstable uninfected equilibrium point exists, the virus could not be eradicated even if this uninfected point is approached. Also, a successful treatment ensures that the infected equilibrium points do not exist, so the system would not be attracted by any one of them if it exists.