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.


Introduction
It is very well known that the adaptive immune response has a significant impact on the progress of the early stage of HBV [1].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][3][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 [5].Understanding the dynamics of the two main arms of the adaptive immune response, CTL cells and Bcells [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 [8], while NAs' function is to elongate DNA which leads to the inhibition of HBV replication [8].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 [8].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 acuteon-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][12][13][14][15][16][17][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 [11].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 2 International Journal of Differential Equations consider the latently infected cells, which are established at an acute stage [11], 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 [16] 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][21][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.

Introducing the Model
We defined the dynamics of the early stage of the HBV by the following system: with where (), (), V(), () , 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 .First, we analyse the model without drug therapy (i.e.,  1 =  2 = 0).More precisely, we consider the following model: We assume that the initial functions of the system of delayed differential equations (3) verify ( () ,  () , V () ,  () ,  ()) ∈ .
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  ≥ 0.
In order to prove boundedness of the solutions of system (3), we consider the following function: if we set  = min(/  , /2, , ℎ, ), we will have Using Gronwall's Lemma, we have that () is bounded, which makes the variables (), (), V(), () , 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  ∈ [0, ).Hence, the boundedness of (), V(), () , and () on [0, ) imply that  = ∞.This completes the proof of the theorem.

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  0 = (  , 0, 0, 0, 0).The second equilibrium This equilibrium exits only if R 0 > 1 and R 1 > 1 with We also have the equilibrium  2 = ( 2 ,  2 , V 2 , 0,  2 ), 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 R ⋆ ≥ 1, then we define International Journal of Differential Equations If R 2 > 1 and R  > 1, then the coordinates of  2 are given by Remark 2. If R 0 > 1, we can easily prove the following: (1) The equilibrium  1 exists if and only if And if R 2 > 1 , this condition combined with the condition R  > 1 could be simplified to From the two previous assessments, the model could have two equilibria  1 and  2 at same time if Finally, it is easy to see that if The third type of equilibrium,  3 , is characterised by no CTL cells response; that is,  = 0.For this reason, we define the threshold by If R ⬦ ≥ 1, we define   and  ℎ (with   ≤  ℎ ) by Hence, the coordinates of   3 = (  3 ,   3 , V 3 ,   3 , 0), with  = , ℎ, are given by Notice that the virus coordinate does not depend on   .However,   3 ,   3 , and   3 are increase functions with respect to   .
Notice that   3 > 0 require that R 0 > 1 and with Φ given by Remark 3. We consider the threshold R  defined by which 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,  4 = ( 4 ,  4 , V 4 ,  4 ,  4 ), where all the coordinates are nonzero.Using the same condition as the previous case, if R ⬦ ≥ 1, then there are two distinct   4 = (  4 ,  4 , V 4 ,  4 ,   4 ), 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 R ⬦ ≥ 1 for   to be feasible, and R  > 1,   > / , and (ℎ − /  )(  / − 1) > 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 by and we have the following results.
Proposition 4. The free-equilibrium point  0 is locally asymptotically stable when R 0 < 1 and unstable when R 0 > 1.
Proof.The characteristic polynomial of the Jacobian matrix (31) at  0 is given by and 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 R 0 < 1.We conclude that the free-equilibrium point  0 is locally asymptotically stable when R 0 < 1 and unstable when R 0 > 1.
Next result will give the condition of stability of the no immune response equilibrium  1 = ( 1 ,  1 , V 1 , 0, 0) , where its coordinates are defined in (14).
Proof.Since the positivity of  2 and  2 depends on the positive sign of R 0 − 1, we conclude that Next, we investigate the case where 1 < R 0 < 1+ − /.Using the Jacobian matrix (31), the characteristic equation at  1 is as follows: International Journal of Differential Equations with Using the form of V 1 and  1 given in ( 14), the two first eigenvalues V 1 −ℎ and  1 − are negative (resp.) if and only if (/ℎ)H 0 < 1 and H 0 < 1 (resp.).
(ii) As the delay  increases, by the inequality 1 < R 0 < 1 +  − /, the quantity R 0 will be a bit bigger than one.
Next, we study the condition of local stability of the equilibrium  2 .
Proof.We can easily notice that R 2 ≤ 1 and then R  ≤ 0, and then  2 < 0 and  2 < 0, which means that  2 does not exist.
On the other hand, if R 2 > 1 and Assume that R 2 > 1 and condition (22) holds.From (31), the characteristic equation at  2 is given by where with It is clear that V 2 − ℎ = ℎ(Φ − 1) is an eigenvalue of   2 .The sign of this eigenvalue is negative if Φ < 1, positive if Φ > 1 , and zero when Φ = 1.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 Φ < 1 .Consequently, if R 2 > 1 and R 2 /(R 2 − 1) < R 0 <     /4 , then  2 is unstable when Φ > 1 and locally asymptotically stable when Φ < 1.
Now, we aim to find the condition of local stability of the equilibrium  ℎ 3 ; we have the following result.
On the other hand, from the Routh-Hurwitz theorem, the other eigenvalues of the above matrix have a negative real part when  ℎ < ( +  − ℎ)/ℎ − .
Consequently, if 3 is locally asymptotically stable.

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  1 and  2 are the weight for the treatment.It is legitimate to assume that two control functions,  1 () and  2 () , 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 ( * 1 ,  * 2 ) such that where  is the control set defined by First, we need to ensure the existence of the optimal control pair.Using the results in Fleming and Rishel [33] and Lukes [34], we have the following theorem.Theorem 10.There exists an optimal control pair ( * 1 ,  * 2 ) ∈  such that The proof of this result is omitted since it is similar to the one in Tridane et al. [16].
Next, via Pontryagin's Minimum Principle [35], we give the necessary conditions for an optimal control problem.We convert solving our optimization problem into maximizing the Hamiltonian  ≡ (,,,V, , ,   , V  ,  1 ,  2 ,   ) point-wisely with respect to  1 and  2 as follows: with And   ,  = 1, 2, 3, 4, 5, are the adjoint functions to be determined.By applying Pontryagin's Minimum Principle in the case system with delay [35], we have the following theorem.

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 [36], used similarly by [37] and enhanced for delay differential equation system [38][39][40].
For the simulation, we use the parameter values given in Table 1.
First, we start our simulation by showing the effect of the delay on the dynamics of the different cells' population as well as the free virions particles.Figure 1 presents the time series of the uninfected cells, the infected cells, the free viruses, and the antibodies.The dashed curves represent the case with delay, while the solid curves show the case without delay.The delay has a clear effect on the dynamics of the early HBV infection by slowing down the overall time series by expending the time between the phases of each curve.However, there is no difference between the two cases as the time passes, which means that the time delay could have an effect on the time scale in planning the treatment period.However, the delay does not lead to periodic dynamics of the model.Hence, the delay cannot cause periodic oscillations.
The next illustrative simulation of the model aims to help in comparing the uninfected cells, the infected cells, the viral load, and the immune response with and without therapy.
Figure 2 shows an increase of the healthy hepatocytes (a) in the first three days, but it is clear that the therapy gives a substantial increase of healthy hepatocytes, with more than 200,000 cells, compared with the case without therapy.
We notice also that, in the absence of the therapy, the number of the infected hepatocytes (b) increases rapidly in the first four days, decreases within twenty days, and increases after 25 days, whereas, in the presence of treatment, the number of infected hepatocytes decreases asymptotically to an undetectable level.More precisely, the number of infected cells with control stabilises at 2.5482, while the number of infected cells without control reaches 2.264 × 10 5 , which makes the drug therapy efficiency in blocking the new infections at 98.73%.
In Figure 2, we see that the number of free virions (a) decreases rapidly towards an undetectable level after introducing the therapy.In fact, with control, the virus stabilises at 1.2112 while without control it reaches 9.768 × 10 6 , which represents a perfect efficiency of the drug therapy in inhibiting the viral production (about 99.99%).
Figure 3(b) shows the antibodies immune response as a function of time.Without the therapy, the antibody level shows relapse in count 50 days after the infection, before it persists over time.We can see clearly that the relapse of the antibody synchronised with the virus peak.On the other hand, the early therapy reduces the burden on the antibody as immune response is barely measured.
The optimal therapy protocol is represented by Figure 4.Each curve presents the optimal drug dosage efficiency and the drug timing during the time of therapy.The optimal therapy requires having a full dosage efficiency for both drugs; the efficiency should be for about 4 days for INF and about 2 weeks for NAs.After 4 days, the INF administration should be stopped and again retaken until it reaches 32% efficiency.Later on, the efficiency can be dropped to less than 10%.For the NAs drugs, after two weeks, the efficiency can be reduced to 50% and eventually dropped to 15% for the rest of the treatment duration.

Conclusion
In this paper, we investigated a mathematical model of the adaptive immune response of the early stage of HBV.The early stage is characterised by a delay in the infection process and a logistic growth of the healthy hepatocyte cells.The aim is to study the role of the two arms of the adaptive immune response, represented by the antibodies and the  CTL cells, in the progress of the HBV infection as the virus gains ground and becomes widespread.Our study showed the possibility of several outcomes depending on many thresholds, which led us to find the conditions of existence of four possible equilibria and investigate their local stability.The stability analysis of these equilibria was very involving and required rigorous calculations.Our mathematical analysis and numerical simulations show that the delay has the effect of slowing down the progress of the disease but does not lead to oscillatory behavior of the dynamics.
As a result of this finding, our next goal was to find the possibility of introducing the actual therapy, which includes standard interferon- and nucleoside analogues.For this purpose, we investigated the optimal control of this therapy via the proposed model.The implementation of such therapy in the early stage instead of the acute stage of HBV infection could be helpful in reducing the burden of the disease.The optimal therapy aims to increase the efficacy of the drug while keeping the healthy hepatocyte cells at the normal level and enhancing the immune response.To solve this problem,  we used the standard techniques to prove the condition of existence of a solution and to find the optimality system.A well-known numerical method was used to solve the optimality system and to identify the best treatment strategy of HBV infection to block new infections and prevent viral production using drug therapy with minimum side effects on the immune response and the healthy hepatocyte cells.

International Journal of Differential Equations
Our numerical results show that the optimal treatment strategies should have high efficiency at the beginning of the therapy, about four days for INF and two weeks for NAs; the efficiency can be adjusted to 10% for INF and to 50% for NAs, and gradually to 15%.
Since there is no clear guideline for the combination therapy in general [41] and for the early infection of HBV in particular, this work should serve as an initial step to consider an early combined use of IFN and NAs in HBV infection.Of course, more pharmacokinetic studies are needed to investigate the long time use of this therapy and the possible risk of treatment failure [8].

Figure 3 :
Figure 3: The HBV as a function of time (a).The antibody response as a function of time (b).

Figure 4 :
Figure 4: The optimal control  1 (a) and the optimal control  2 (b) versus time.
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,  1 represents the efficiency of IFN in reducing the new infected cells; the infection rate in the presence of the drug is (1− 1 ), while  2 is the efficiency of NAs in blocking the reverse transcription, such that the virions production rate under this drug is (1 −  2 ).
The healthy hepatocytes become infected by the virus at rate (V/), 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 .