Mathematical Modeling in Virology by Differential EquationsView this Special Issue
Research Article | Open Access
Mathematical Modeling of the Adaptive Immune Responses in the Early Stage of the HBV Infection
The aim of this paper is to study the early stage of HBV infection and impact delay in the infection process on the adaptive immune response, which includes cytotoxic T-lymphocytes and antibodies. In this stage, the growth of the healthy hepatocyte cells is logistic while the growth of the infected ones is linear. To investigate the role of the treatment at this stage, we also consider two types of treatment: interferon- (IFN) and nucleoside analogues (NAs). To find the best strategy to use this treatment, an optimal control approach is developed to find the possibility of having a functional cure to HBV.
It is very well known that the adaptive immune response has a significant impact on the progress of the early stage of HBV . This response can either lead to complete cure from the infection, and it is characterised by the production of neutralizing antibodies against HBV surface antigen (HBsAg) and adequate cytotoxic lymphocyte T-cell (CTL) responses [2–4], or it could result in chronic infection that leads to liver cancer (HHC), cirrhosis, or liver failure.
During the incubation period of HBV, which is 30 to 180 days, the dynamics of the adaptive immune response are not fully understood, since the majority of the cases are clinically known after the infection is established and the patient is in the acute stage . Understanding the dynamics of the two main arms of the adaptive immune response, CTL cells and B-cells [6, 7], will help grasp how the virus escapes the adaptive immunity and improve the ability of the immune system to control the virus in early infection.
Moreover, it is known that the actual therapy, which includes the standard interferon- (INF) and the nucleoside analogues (NAs), is also initiated during the acute stage of HBV infection. INF helps eliminate the infected cells by reducing cccDNA , while NAs’ function is to elongate DNA which leads to the inhibition of HBV replication . As monotherapy, the NAs come in different types of drugs (Entecavir, Adefovir, and Lamivudine), which are also known for enhancing the functions of natural killers . The question is, what if we could initiate therapy even earlier? Is it possible to eradicate the virus within this period with the therapy? In fact, recent studies [9, 10] showed that early Lamivudine treatment could lead to better outcome in acute-on-chronic stages and with less liver damage. Therefore, our goal is to understand the dynamics of the adaptive immune response, via the CTL cells and B-cells, in the early stage of HBV, and investigate the impact of early HBV treatment therapy on disease progress via a mathematical model of virus-immune response.
Mathematical modeling of the immune response to HBV is a subject that has been heavily investigated over the years by many authors [11–18], just to name a few. To our knowledge, there is no study that investigates the adaptive immune response in the early stage of the infection and effect of the early treatment on the progress of the disease. In this work, we are aiming to investigate this issue by considering an augmented model of our recent works [18, 19], and we consider the logistic growth only for the healthy hepatocyte cells and the infected hepatocyte cells . This assumption is made to reflect the nature of the growth of these two types of cells in the early stage of the infection. We also did not consider the latently infected cells, which are established at an acute stage , and we did not consider noncytolytic carrying processes since no data support such assumption in this stage. Moreover, we have also considered a more generalized incident function  and the delay in this incident function to reflect the time between the infection and the cells becoming productively infected (infected and producing new viruses). The optimal control of the HBV therapy aims to find the optimal strategy of the drugs that allow blocking the virus production and infection. Several papers studied the optimal control of the HBV therapy [16, 20–22]. In our case, the therapy will have an antiviral effect, and we ignore its immunomodulatory effect since we do not know what impact the use of the therapy could have on the immune system in the early stage of HBV infection.
The paper is organized as follows. In Section 2, we introduce our model, and we investigate the basic properties of the model without therapy, which includes positivity and boundedness of solutions. In Section 3, we focus on the stability analysis of the different types of steady states. Next, we will investigate optimal control of the treatment therapy, and we will numerically solve the optimality conditions. Finally, we will give a discussion and a conclusion to our work.
2. Introducing the Model
We defined the dynamics of the early stage of the HBV by the following system:withwhere , , , , and denote the concentrations of uninfected cells, infected cells, viruses, antibodies, and cytotoxic T-lymphocytes (CTLs), respectively. The uninfected hepatocytes grow at a rate that depends on the liver size, , at a maximum per capita proliferation rate . The healthy hepatocytes become infected by the virus at rate , where is a constant. Infected cells die at rate and are killed by the CTLs response at rate . The infected non-virus-producing cells have a death rate ; these cells start producing viruses after delay time , hence is the probability of survival between time and . The free virus particles are produced at rate , where is the number of free virions produced by the infected cells during their lifespan, and decay at rate . Parameter represents the rate of expansion of CTL cell and is its decay rate in the absence of antigenic stimulation. The antibodies are developed in response to free virus at rate and decay at rate . Finally, represents the efficiency of IFN in reducing the new infected cells; the infection rate in the presence of the drug is , while is the efficiency of NAs in blocking the reverse transcription, such that the virions production rate under this drug is .
First, we analyse the model without drug therapy (i.e., ). More precisely, we consider the following model:
Let be the Banach space of continuous mapping from to with respect to the norm We assume that the initial functions of the system of delayed differential equations (3) verify Following the standard approach, we assume that
Under these initial conditions, the solutions of (3) satisfy the following theorem.
Theorem 1. System (3) has a unique solution; in addition, this solution is nonnegative and bounded for all .
Proof. Notice that system (3) is locally Lipschitzian at , which implies that the solution of system (3), subject to (7), exists and is unique on , where is the maximal existence time for the solution of system (3). Notice that if , then for all . Hence, we can assume that . Notice also that if , then, from (6), we have t, which implies that, for small , we have . Similarly, if , then , which implies that, for small , we have . Moreover, if , , then , for all . Hence, we assume below that , .
Assume first that there is such that and , , , for . Observe that it is easy to show that for ; we can see that , and clearly , for ; these observations imply that, for , we have
Hence,which contradicts our assumption.
Following a similar approach, we can prove that all the variables of system (3) are positive.
In order to prove boundedness of the solutions of system (3), we consider the following function:From (3), we have since , , for , it follows that if we set , we will haveUsing Gronwall’s Lemma, we have that is bounded, which makes the variables , , , , and bounded, which makes us unsure that the solution exists globally.
The above results show that the components of the solution of system (3) are nonnegative for all . Hence, the boundedness of , and on imply that . This completes the proof of the theorem.
3. Equilibrium Points and Their Stability
First, we aim to find all possible equilibria points. It is clear that our model (3) has disease-free equilibrium . The second equilibrium represents the no immune response equilibrium withThis equilibrium exits only if and with
We also have the equilibrium , which represents an equilibrium where CTL cells are the only adaptive immune response and B-cells are zero. To define such equilibrium, we introduce the following thresholds:If , then we define
If and , then the coordinates of are given by
Remark 2. If , we can easily prove the following: (1)The equilibrium exists if and only if (2) is equivalent to And if , this condition combined with the condition could be simplified to From the two previous assessments, the model could have two equilibria and at same time ifFinally, it is easy to see that if or , then .
The third type of equilibrium, , is characterised by no CTL cells response; that is, . For this reason, we define the threshold by
If , we define and (with ) by
Hence, the coordinates of , with , are given by
Notice that the virus coordinate does not depend on . However, , , and are increase functions with respect to .
Notice that require that and with given by
Remark 3. We consider the threshold defined bywhich represents the survival rate of the virus, with ignoring the antibody effect, times the survival rate of the antibody, over the survival rate of the CTL cells .
It is easy to see that
Finally, we have the endemic equilibria, , where all the coordinates are nonzero. Using the same condition as the previous case, if , then there are two distinct , with , with the coordinate given by and and are defined in (24).
The endemic equilibria are characterised by two possible levels of the healthy cells and corresponding CTL cells. On the other hand, coordinates of the endemic equilibria are constant with respect to the rest of the variables. It is important to mention that the existence of these two endemic equilibria requires for to be feasible, and , , and .
3.1. The Stability Analysis
In this section, we investigate the condition of stability of each possible equilibria point. First, the Jacobian matrix of system (3) is given byand we have the following results.
Proposition 4. The free-equilibrium point is locally asymptotically stable when and unstable when .
Proof. The characteristic polynomial of the Jacobian matrix (31) at is given byand then the eigenvalues of the Jacobian matrix at are It is clear that the first four eigenvalues are negative. The fifth one is negative when . We conclude that the free-equilibrium point is locally asymptotically stable when and unstable when .
Next result will give the condition of stability of the no immune response equilibrium , where its coordinates are defined in (14).
Theorem 5. (1) If , then the point does not exist.
(2) If , then .
(3) If , then is locally asymptotically stable if ; it is unstable for , with
Proof. Since the positivity of and depends on the positive sign of , we conclude that does not exist if . Moreover, if , it is easy to say that .
Next, we investigate the case where . Using the Jacobian matrix (31), the characteristic equation at is as follows: withUsing the form of and given in (14), the two first eigenvalues and are negative (resp.) if and only if and (resp.).
On the other hand, from the Routh-Hurwitz theorem, the other eigenvalues of the above matrix have a negative real part when .
Remark 6. (i) If , the condition can be replaced by .
(ii) As the delay increases, by the inequality , the quantity will be a bit bigger than one.
Next, we study the condition of local stability of the equilibrium .
Theorem 7. Assume that ; then the following applies: (1)If , then does not exist.(2)If (a)If , then does not exist.(b)If , then .(c)If (i)If , then is locally asymptotically stable.(ii)If , then is unstable, where is defined in (27).
Proof. We can easily notice that and then , and then and , which means that does not exist.
On the other hand, if and , then , and if and , then and then and .
Assume that and condition (22) holds. From (31), the characteristic equation at is given bywhere with It is clear that is an eigenvalue of . The sign of this eigenvalue is negative if , positive if , and zero when . On the other hand, from the Routh-Hurwitz theorem applied to the fourth-order polynomial in the characteristic equation, the other eigenvalues of the above matrix have negative real parts when . Consequently, if and , then is unstable when and locally asymptotically stable when .
Now, we aim to find the condition of local stability of the equilibrium ; we have the following result.
Theorem 8. Assume that and :(1)If , then equilibria for do not exist and when .(2)If , then are locally asymptotically stable.(3)If , then is unstable.
Proof. It is easy to see that if , then the equilibrium does not exist and if the two points and coincide.
If , using the Jacobian matrix (31), we get the following characteristic equation at :where withIt is clear that is an eigenvalue of . The sign of this eigenvalue is negative if , which is equivalent to . The sign of this eigenvalue is positive if , which will give, with , the condition of instability of the theorem.
On the other hand, from the Routh-Hurwitz theorem, the other eigenvalues of the above matrix have a negative real part when
Consequently, if , then is locally asymptotically stable.
Theorem 9. (1) If or , then the point with does not exist. Moreover, when and when .
(2) If and , then is locally asymptotically stable.
4. The Optimal Control Therapy Analysis
In this section, we consider the optimal control of the HBV drug therapy; as we mentioned previously, the therapy has an antiviral effect by reducing the viral production rate and blocking the shedding and bending of the virus to the uninfected cells. For this purpose, we consider the controlled version of system (3) defined as follows:
The optimization problem that we consider is to maximize the following objective functional: where stands for the time period of treatment. The two positive constants and are the weight for the treatment. It is legitimate to assume that two control functions, and , are bounded and Lebesgue integrable. These assumptions align with the fact that the drug has a limited dosage and time to use.
The goal is to decrease the viral load while increasing the number of the uninfected cells and maximizing the immune responses. This should be done with minimizing the cost of treatment. We can achieve this goal by maximizing the objective functional defined in (45), which means finding the optimal control pair such thatwhere is the control set defined by
Theorem 10. There exists an optimal control pair such that
The proof of this result is omitted since it is similar to the one in Tridane et al. .
Next, via Pontryagin’s Minimum Principle , we give the necessary conditions for an optimal control problem. We convert solving our optimization problem into maximizing the Hamiltonian point-wisely with respect to and as follows:withAnd , are the adjoint functions to be determined. By applying Pontryagin’s Minimum Principle in the case system with delay , we have the following theorem.
Theorem 11. Given optimal controls and solutions , , , , and of the corresponding state system (3), there exist adjoint variables, , , , , and satisfying the equationswith being an indicator function and also the transversality conditions Moreover, the optimal control is given by
5. Numerical Simulations
In order to solve our optimization system, we use a numerical schema based on the forward and backward finite difference approximation. This schema was originally presented in the case of ODE system in , used similarly by  and enhanced for delay differential equation system [38–40].
We consider the step size and with and . We take knots to left of and right of , to get the following partition: which gives . The state and the adjoint variables are , , , , , , , , , and and the controls are , in terms of nodal points , , , , , , , ,