# Mathematical Modeling of the Adaptive Immune Responses in the Early Stage of the HBV Infection

Authors: Karam Allali; Adil Meskaf; Abdessamad Tridane
Journal: International Journal of Differential Equations
 (2018)

---

## Abstract

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.


---

## Body

## 1. 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–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 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 [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 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 [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 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–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 Section2, 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:(1)dxdt=rxt1-TtTm-β1-u1tvtxtTt,dydt=βe-kτ1-u1tvt-τxt-τTt-τ-ayt-pytzt,dvdt=1-u2taNyt-δvt-qvtwt,dwdt=gvtwt-hwt,dzdt=cytzt-bztwith(2)Tt=xt+yt,where x(t), y(t), v(t), w(t), and z(t) 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, Tm, at a maximum per capita proliferation rate r. The healthy hepatocytes become infected by the virus at rate βvx/T, where β is a constant. Infected cells y die at rate a and are killed by the CTLs response at rate p. The infected non-virus-producing cells have a death rate k; these cells start producing viruses after delay time τ, hence e-kτ is the probability of survival between time t-τ and t. The free virus particles are produced at rate aN, where N is the number of free virions produced by the infected cells during their lifespan, and decay at rate δ. Parameter c represents the rate of expansion of CTL cell z and b is its decay rate in the absence of antigenic stimulation. The antibodies are developed in response to free virus at rate g and decay at rate h. Finally, u1 represents the efficiency of IFN in reducing the new infected cells; the infection rate in the presence of the drug is (1-u1)β, while u2 is the efficiency of NAs in blocking the reverse transcription, such that the virions production rate under this drug is (1-u2)aN.First, we analyse the model without drug therapy (i.e.,u1=u2=0). More precisely, we consider the following model: (3)dxdt=rxt1-TtTm-βvtxtTt,dydt=βe-kτvt-τxt-τTt-τ-ayt-pytzt,dvdt=aNyt-δvt-qvtwt,dwdt=gvtwt-hwt,dzdt=cytzt-bzt.LetX=C-τ,0;R5 be the Banach space of continuous mapping from [-τ,0] to R5 with respect to the norm (4)φ=sup-τ≤t≤0⁡φt.We assume that the initial functions of the system of delayed differential equations (3) verify (5)xt,yt,vt,zt,wt∈X.Following the standard approach, we assume that (6)xt>0,yt≥0,vt≥0,zt≥0,wt≥0,fort∈-τ,0,(7)Tm≥Tt=xt+yt>0,fort∈-τ,0.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 t≥0.Proof.
Notice that system (3) is locally Lipschitzian at t=0, which implies that the solution of system (3), subject to (7), exists and is unique on 0,b, where b is the maximal existence time for the solution of system (3). Notice that if x(0)=0, then x(t)≡0 for all t>0. Hence, we can assume that x(0)>0. Notice also that if y(0)=0, then, from (6), we have y′(0)=βv(-τ)x(-τ)/T(-τ)≥0t, which implies that, for small t>0, we have y(t)>0. Similarly, if v(0)=0, then v′(0)=aNy(0)>0, which implies that, for small t>0, we have v(t)>0. Moreover, if w(0)=0, z(0)=0, then w(t)≡0, z(t)≡0 for all t>0. Hence, we assume below that w(0)>0, z(0)>0.
Assume first that there isb>t1>0 such that x(t1)=0 and x(t)>0, y(t)>0, v(t)>0, for t∈[0,t1]. Observe that (8)dxtdt=rxt1-TtTm-βvtxtTt;it is easy to show that 0<T(t)<Tm for t∈[0,t1]; we can see that dx(t)/dt≥-βvtxt/Tt, and clearly y(t)<T(t), for t∈[0,t1]; these observations imply that, for t∈[0,t1], we have dx(t)/dt≥-βvtxt/yt.
Hence,(9)xt1≥x0e-∫0t1βvs/ysds>0,
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:
(10)
F
t
=
N
c
g
e
-
k
τ
x
t
+
N
c
g
y
t
+
τ
+
c
g
2
v
t
+
τ
+
c
q
2
w
t
+
τ
+
N
g
p
z
t
+
τ
.From (3), we have (11)dFtdt=Ncge-kτrxt-rxtTtTm-βvtxtTt+Ncgβe-kτvtxtTt-ayt+τ-pyt+τzt+τ+cg2aNyt+τ-δvt+τ-qvt+τwt+τ+cq2gvt+τwt+τ-hwt+τ+Ngpcyt+τzt+τ-bzt+τ;since 0<T(t)<Tm, x(t)<Tm, -x(t)T(t)<-x(t) for t>0, it follows that (12)dFtdt≤Ncge-kτrTm-Ncge-kτrTmxt-aNcg2yt+τ-δcg2vt+τ-hcq2wt+τ-Ngpbzt+τ;if we set ϱ=min⁡(r/Tm,a/2,δ,h,b), we will have(13)dFtdt≤Ncge-kτ-ϱFt.Using Gronwall’s Lemma, we have that F(t) is bounded, which makes the variables x(t), y(t), v(t), w(t), and z(t) 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 t∈0,b. Hence, the boundedness of T(t),v(t),w(t), and z(t) on 0,b imply that b=∞. 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 E0=(Tm,0,0,0,0). The second equilibrium E1=(x1,y1,v1,0,0) represents the no immune response equilibrium with(14)x1=TmβaNR0δrR1-1,y1=TmβaNR0δrR0-1R1-1,v1=aN2βTmδ2R0rR0-1R1-1.This equilibrium exits only if R0>1 and R1>1 with(15)R0=βe-kτNδ,R1=δrR0+βaNβaNR0=re-kτa+11R0.We also have the equilibriumE2=(x2,y2,v2,0,z2), 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:(16)R⋆=Tmδcr4βaNb.If R⋆≥1, then we define (17)R2=cTm2b1+1-1R⋆,Rz=R01-1R2.IfR2>1 and Rz>1, then the coordinates of E2 are given by (18)x2=bcR2-1,y2=bc,v2=aNbδc,z2=apRz-1.Remark 2.
IfR0>1, we can easily prove the following: (1)
The equilibriumE1 exists if and only if (19)1<R0<re-kτa+1.(2)
R⋆>1 is equivalent to(20)1<R0<Tmcre-kτ4ab.And ifR2>1, this condition combined with the condition Rz>1 could be simplified to (21)1<R2R2-1<R0<Tmcre-kτ4ab. From the two previous assessments, the model could have two equilibria E1 and E2 at same time if(22)1<R2R2-1<R0<min⁡Tmcre-kτ4ab,re-kτa+1.
Finally, it is easy to see that ifb/c≤Tm/2 or b/c≥Tm, then R2>1.The third type of equilibrium,E3, is characterised by no CTL cells response; that is, z=0. For this reason, we define the threshold by (23)R⋄=Tmgr4βh.IfR⋄≥1, we define Tl and Th (with Tl≤Th) by(24)Th=Tm21+1-1R⋄,Tl=Tm21-1-1R⋄.Hence, the coordinates ofE3i=(x3i,y3i,v3,w3i,0), with i=l,h, are given by(25)x3i=agTi2agTi+βe-kτhy3i=βe-kτhTiagTi+βe-kτhv3=hgw3i=aNgqβe-kτTiagTi+βe-kτh-δq.Notice that the virus coordinate does not depend onTi. However, x3i, y3i, and w3i are increase functions with respect to Ti.Notice thatw3i>0 require that R0>1 and (26)Ti>1ΦR0R0-1with Φ given by(27)Φ=aNδ·gh.Remark 3.
We consider the thresholdRT defined by(28)RT=aNgbδhc=Φbc,which represents the survival rate of the virus, with ignoring the antibody effect, aN/δ times g/h the survival rate of the antibody, over the survival rate of the CTL cells c/b.
It is easy to see that(29)R∗>R⋄⟺RT<1resp.R∗<R⋄⟺RT>1.Finally, we have the endemic equilibria,E4=(x4,y4,v4,w4,z4), where all the coordinates are nonzero. Using the same condition as the previous case, if R⋄≥1, then there are two distinct E4i=(x4i,y4,v4,w4,z4i), with i=l,h, with the coordinate given by (30)x4i=Ti-bc,y4=bc,v4=hg,w4=δqRT-1,z4i=βe-kτhcpbgTiTi-bc-ap,and Tl and Th 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 requiresR⋄≥1 for Ti to be feasible, and RT>1, Ti>b/c, and βhe-kτ/agTiTic/b-1>1.


### 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 by(31)r1-2x+yTm-βvyx+y2-rxTm+βvxx+y2-βxx+y00βe-kτvyx+y2-βe-kτvxx+y2-a-pzβe-kτxx+y0-py0aN-δ-qw-qv000gwgv-h00cz00cy-band we have the following results.Proposition 4.
The free-equilibrium pointE0 is locally asymptotically stable when R0<1 and unstable when R0>1.Proof.
The characteristic polynomial of the Jacobian matrix (31) at E0 is given by(32)PE0λ=λ+rλ+b-λ-hλ2+a+δλ+aδ1-R0,and then the eigenvalues of the Jacobian matrix at Ef are (33)-r,-b,-h,-12a+δ+a+δ2+4aδR0-1,-12a+δ-a+δ2+4aδR0-1.
It is clear that the first four eigenvalues are negative. The fifth one is negative whenR0<1. We conclude that the free-equilibrium point E0 is locally asymptotically stable when R0<1 and unstable when R0>1.Next result will give the condition of stability of the no immune response equilibriumE1=(x1,y1,v1,0,0), where its coordinates are defined in (14).Theorem 5.
(1) IfR0<1, then the point E1 does not exist.
(2) IfR0=1, then E1=E0.
(3) If1<R0<1+re-kτ/a, then E1 is locally asymptotically stable if min⁡H0,gaN/hδH0<1; it is unstable for min(H0,gaN/hδH0)>1, with(34)H0=TmβaNR0δrR0-1R1-1.Proof.
Since the positivity ofy2 and z2 depends on the positive sign of R0-1, we conclude that E1 does not exist if R0<1. Moreover, if R0=1, it is easy to say that E1=Ef.Next, we investigate the case where1<R0<1+re-kτ/a. Using the Jacobian matrix (31), the characteristic equation at E1 is as follows: (35)PE1λ=cy1-b-λgv1-h-λλ3+a1λ2+a2λ+a3,a1=βe-kτTvx1+δ+βTvy1+2x1+y1Tmr-r,=βe-kτTvx1+δ+βTvy1+βaNR0δrre-kτa+11R0-1R0+1.a2=βe-kτTvx1δ+2x1+y1Tmr-r+βTvy1+2x1+y1Tmr-r+βTvy1a+δ-aNβe-kτTx,a3=βaNR0δrre-kτa+11R0-1R0+1+βTvy1δβe-kτTvx1+aδ+aNβe-kτTx+aNβ2e-kτTvy1Tx,with(36)Tv=v1x1+y12,Tx=x1x1+y1.Using the form of v1 and y1 given in (14), the two first eigenvalues gv1-h and cy1-b are negative (resp.) if and only if gaN/hδH0<1 and H0<1 (resp.).On the other hand, from the Routh-Hurwitz theorem, the other eigenvalues of the above matrix have a negative real part when1<R0<1+re-kτ/a.Remark 6.
(i) Ifhδ/aNg>1, the condition min⁡H0,gaN/hδH0<1 can be replaced by H0<1<hδ/aNg.
(ii) As the delayτ increases, by the inequality 1<R0<1+re-kτ/a, the quantity R0 will be a bit bigger than one.Next, we study the condition of local stability of the equilibriumE2.Theorem 7.
Assume thatR0>1; then the following applies: (1)
IfR2≤1, then E2 does not exist.(2)
IfR2>1
(a)
IfR2/R2-1>R0, then E2 does not exist.(b)
IfR2/R2-1=R0, then E2=E1.(c)
IfR2/R2-1<R0<Tmcrekτ/4ab
(i)
IfΦ<1, then E2 is locally asymptotically stable.(ii)
IfΦ>1, then E2 is unstable, where Φ is defined in (27).Proof.
We can easily notice thatR2≤1 and then Rz≤0, and then x2<0 and z2<0, which means that E2 does not exist.
On the other hand, ifR2>1 and 1<R0<R2/R2-1, then z2<0, and if R2>1 and 1<R0=R2/R2-1, then Rz=1 and then z2=w2=0 and E2=E1.
Assume thatR2>1 and condition (22) holds. From (31), the characteristic equation at E2 is given by(37)PE2λ=gv2-h-λλ4+b1λ3+b2λ2+b3λ+b4,where (38)b1=δ+βe-kτTvx2+a+pz2+βTvy2+2x2+y2Tmr-r,b2=δβe-kτTvx2+aδ+pδz2-aNβe-kτTx+py2+βTvy2+2x2+y2Tmr-rδ+βe-kτTvx2+a+pz2,b3=b-cy2δβe-kτTvx2+aδ+pδz2-aNβe-kτTx+py2cz2δ+βTvy2+2x2+y2Tmr-rδβe-kτTvx2+aδ+pδz2-aNβe-kτTx+py2,b4=βTvy2+2x2+y2Tmr-rcy2-bδβe-kτTvx2+aδ+pδz2-aNβe-kτTx-βTvy22x2+y2Tmr-rpy2cz2δ,with (39)Tv=v1x1+y12,Tx=x1x1+y1.It is clear that gv2-h=h(Φ-1) is an eigenvalue of JE2. 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 R2>1 and R2/R2-1<R0<Tmcrekτ/4ab, then E2 is unstable when Φ>1 and locally asymptotically stable when Φ<1.Now, we aim to find the condition of local stability of the equilibriumE3h; we have the following result.Theorem 8.
Assume thatR0>1 and R⋄>1:(1)
IfTh<βe-τkh/ag1/R0-1, then equilibria E3i for i=l,h do not exist and E3h=E1 when Th=βe-τkh/ag1/R0-1.(2)
Ifβe-τkh/ag1/R0-1<Th<1/R0(gN/δch+R0/c), then E3h are locally asymptotically stable.(3)
IfTh>maxβe-τkh/ag1/R0-1,1/R0gN/δch+R0/c, then E3 is unstable.Proof.
It is easy to see that ifTh<1/R0g+βhe-τk/ag, then the equilibrium E3h does not exist and if Th=βe-τkh/ag1/R0-1 the two points E3h and E1 coincide.
IfTh>βe-τkh/ag1/R0-1, using the Jacobian matrix (31), we get the following characteristic equation at E3i:(40)PE3iλ=cy3i-b-λλ4+c1λ3+c2λ2+c3λ+c4,where (41)c1=βe-kτTvx3i+a+δ+2x3i+y3iTmr-r+βTvy3i+qw3i,c2=qv3gw3i+δ+qw3ih-gv3+βe-kτTvx3i+aδ+qw3i-gv3+h+aNβe-kτTx+2x3i+y3iTmr-r+βTvy3iβe-kτTvx3i+a-qw3i+gv3-δ-h-βe-kτTvy3i-rx1iTm+βTvx3i,c3=-βe-kτTvx3i-aδ+qw3igv3-h+βe-kτTvx3i+aqv3gw3i+aNβe-kτTxgv3-h+-2x3i+y3iTmr+r-βTvy3i-βe-kτTvx3i-aδ+qw3i-gv3+h+aNβe-kτTx+-2x3i+y3iTmr+r-βTvy3iδ+qw3igv3-h+2x3i+y3iTmr-r+βTvy3iqv3gw3i+βe-kτTvy3iβTxaN+gv3-h-rx1iTm+βTvx3i+-rx1Tm+βTvx3iδ+qw3i,c4=-2x3i+y3iTmr+r-βTvy3i-βe-kτTvx3i-aδ+qw3ih-gv3+qv3gw3i+-rx3iTm+βTvx3iδ+qw3igv3βe-kτTvy3i-βe-kτTvy3iqv3gw3i-rx3iTm+βTvx3i+gv3-hβTxaN, with(42)Tv=v3x3i+y3i2,Tx=x3ix3i+y3i.It is clear that cy3h-b=bcβhe-kτTh/bg+βe-kτh-1 is an eigenvalue of JE3h. The sign of this eigenvalue is negative if Th<b(g+βe-kτh)/cβhe-kτ, which is equivalent to Th<1/R0gN/δch+R0/c. The sign of this eigenvalue is positive if Th>1/R0gN/δch+R0/c, which will give, with Th>1/R0g+βhe-τk/ag, 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 whenTh<b(g+βe-kτh)/cβhe-kτ.
Consequently, ifβe-τkh/ag1/R0-1<Th<1/R0gN/δch+R0/c, then E3h is locally asymptotically stable.Theorem 9.
(1) IfΦ<1 or Hiw,z<1, then the point E4i with i=l,h does not exist. Moreover, E4i=E2 when Φ=1 and E4i=E2 when Hiw,z=1.
(2) IfΦ>1 and Hiw,z>1, then E4i is locally asymptotically stable.
Here(43)Hiw,z=cβe-kτhTibg+βe-kτh;i=l,h.


## 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 by(31)r1-2x+yTm-βvyx+y2-rxTm+βvxx+y2-βxx+y00βe-kτvyx+y2-βe-kτvxx+y2-a-pzβe-kτxx+y0-py0aN-δ-qw-qv000gwgv-h00cz00cy-band we have the following results.Proposition 4.
The free-equilibrium pointE0 is locally asymptotically stable when R0<1 and unstable when R0>1.Proof.
The characteristic polynomial of the Jacobian matrix (31) at E0 is given by(32)PE0λ=λ+rλ+b-λ-hλ2+a+δλ+aδ1-R0,and then the eigenvalues of the Jacobian matrix at Ef are (33)-r,-b,-h,-12a+δ+a+δ2+4aδR0-1,-12a+δ-a+δ2+4aδR0-1.
It is clear that the first four eigenvalues are negative. The fifth one is negative whenR0<1. We conclude that the free-equilibrium point E0 is locally asymptotically stable when R0<1 and unstable when R0>1.Next result will give the condition of stability of the no immune response equilibriumE1=(x1,y1,v1,0,0), where its coordinates are defined in (14).Theorem 5.
(1) IfR0<1, then the point E1 does not exist.
(2) IfR0=1, then E1=E0.
(3) If1<R0<1+re-kτ/a, then E1 is locally asymptotically stable if min⁡H0,gaN/hδH0<1; it is unstable for min(H0,gaN/hδH0)>1, with(34)H0=TmβaNR0δrR0-1R1-1.Proof.
Since the positivity ofy2 and z2 depends on the positive sign of R0-1, we conclude that E1 does not exist if R0<1. Moreover, if R0=1, it is easy to say that E1=Ef.Next, we investigate the case where1<R0<1+re-kτ/a. Using the Jacobian matrix (31), the characteristic equation at E1 is as follows: (35)PE1λ=cy1-b-λgv1-h-λλ3+a1λ2+a2λ+a3,a1=βe-kτTvx1+δ+βTvy1+2x1+y1Tmr-r,=βe-kτTvx1+δ+βTvy1+βaNR0δrre-kτa+11R0-1R0+1.a2=βe-kτTvx1δ+2x1+y1Tmr-r+βTvy1+2x1+y1Tmr-r+βTvy1a+δ-aNβe-kτTx,a3=βaNR0δrre-kτa+11R0-1R0+1+βTvy1δβe-kτTvx1+aδ+aNβe-kτTx+aNβ2e-kτTvy1Tx,with(36)Tv=v1x1+y12,Tx=x1x1+y1.Using the form of v1 and y1 given in (14), the two first eigenvalues gv1-h and cy1-b are negative (resp.) if and only if gaN/hδH0<1 and H0<1 (resp.).On the other hand, from the Routh-Hurwitz theorem, the other eigenvalues of the above matrix have a negative real part when1<R0<1+re-kτ/a.Remark 6.
(i) Ifhδ/aNg>1, the condition min⁡H0,gaN/hδH0<1 can be replaced by H0<1<hδ/aNg.
(ii) As the delayτ increases, by the inequality 1<R0<1+re-kτ/a, the quantity R0 will be a bit bigger than one.Next, we study the condition of local stability of the equilibriumE2.Theorem 7.
Assume thatR0>1; then the following applies: (1)
IfR2≤1, then E2 does not exist.(2)
IfR2>1
(a)
IfR2/R2-1>R0, then E2 does not exist.(b)
IfR2/R2-1=R0, then E2=E1.(c)
IfR2/R2-1<R0<Tmcrekτ/4ab
(i)
IfΦ<1, then E2 is locally asymptotically stable.(ii)
IfΦ>1, then E2 is unstable, where Φ is defined in (27).Proof.
We can easily notice thatR2≤1 and then Rz≤0, and then x2<0 and z2<0, which means that E2 does not exist.
On the other hand, ifR2>1 and 1<R0<R2/R2-1, then z2<0, and if R2>1 and 1<R0=R2/R2-1, then Rz=1 and then z2=w2=0 and E2=E1.
Assume thatR2>1 and condition (22) holds. From (31), the characteristic equation at E2 is given by(37)PE2λ=gv2-h-λλ4+b1λ3+b2λ2+b3λ+b4,where (38)b1=δ+βe-kτTvx2+a+pz2+βTvy2+2x2+y2Tmr-r,b2=δβe-kτTvx2+aδ+pδz2-aNβe-kτTx+py2+βTvy2+2x2+y2Tmr-rδ+βe-kτTvx2+a+pz2,b3=b-cy2δβe-kτTvx2+aδ+pδz2-aNβe-kτTx+py2cz2δ+βTvy2+2x2+y2Tmr-rδβe-kτTvx2+aδ+pδz2-aNβe-kτTx+py2,b4=βTvy2+2x2+y2Tmr-rcy2-bδβe-kτTvx2+aδ+pδz2-aNβe-kτTx-βTvy22x2+y2Tmr-rpy2cz2δ,with (39)Tv=v1x1+y12,Tx=x1x1+y1.It is clear that gv2-h=h(Φ-1) is an eigenvalue of JE2. 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 R2>1 and R2/R2-1<R0<Tmcrekτ/4ab, then E2 is unstable when Φ>1 and locally asymptotically stable when Φ<1.Now, we aim to find the condition of local stability of the equilibriumE3h; we have the following result.Theorem 8.
Assume thatR0>1 and R⋄>1:(1)
IfTh<βe-τkh/ag1/R0-1, then equilibria E3i for i=l,h do not exist and E3h=E1 when Th=βe-τkh/ag1/R0-1.(2)
Ifβe-τkh/ag1/R0-1<Th<1/R0(gN/δch+R0/c), then E3h are locally asymptotically stable.(3)
IfTh>maxβe-τkh/ag1/R0-1,1/R0gN/δch+R0/c, then E3 is unstable.Proof.
It is easy to see that ifTh<1/R0g+βhe-τk/ag, then the equilibrium E3h does not exist and if Th=βe-τkh/ag1/R0-1 the two points E3h and E1 coincide.
IfTh>βe-τkh/ag1/R0-1, using the Jacobian matrix (31), we get the following characteristic equation at E3i:(40)PE3iλ=cy3i-b-λλ4+c1λ3+c2λ2+c3λ+c4,where (41)c1=βe-kτTvx3i+a+δ+2x3i+y3iTmr-r+βTvy3i+qw3i,c2=qv3gw3i+δ+qw3ih-gv3+βe-kτTvx3i+aδ+qw3i-gv3+h+aNβe-kτTx+2x3i+y3iTmr-r+βTvy3iβe-kτTvx3i+a-qw3i+gv3-δ-h-βe-kτTvy3i-rx1iTm+βTvx3i,c3=-βe-kτTvx3i-aδ+qw3igv3-h+βe-kτTvx3i+aqv3gw3i+aNβe-kτTxgv3-h+-2x3i+y3iTmr+r-βTvy3i-βe-kτTvx3i-aδ+qw3i-gv3+h+aNβe-kτTx+-2x3i+y3iTmr+r-βTvy3iδ+qw3igv3-h+2x3i+y3iTmr-r+βTvy3iqv3gw3i+βe-kτTvy3iβTxaN+gv3-h-rx1iTm+βTvx3i+-rx1Tm+βTvx3iδ+qw3i,c4=-2x3i+y3iTmr+r-βTvy3i-βe-kτTvx3i-aδ+qw3ih-gv3+qv3gw3i+-rx3iTm+βTvx3iδ+qw3igv3βe-kτTvy3i-βe-kτTvy3iqv3gw3i-rx3iTm+βTvx3i+gv3-hβTxaN, with(42)Tv=v3x3i+y3i2,Tx=x3ix3i+y3i.It is clear that cy3h-b=bcβhe-kτTh/bg+βe-kτh-1 is an eigenvalue of JE3h. The sign of this eigenvalue is negative if Th<b(g+βe-kτh)/cβhe-kτ, which is equivalent to Th<1/R0gN/δch+R0/c. The sign of this eigenvalue is positive if Th>1/R0gN/δch+R0/c, which will give, with Th>1/R0g+βhe-τk/ag, 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 whenTh<b(g+βe-kτh)/cβhe-kτ.
Consequently, ifβe-τkh/ag1/R0-1<Th<1/R0gN/δch+R0/c, then E3h is locally asymptotically stable.Theorem 9.
(1) IfΦ<1 or Hiw,z<1, then the point E4i with i=l,h does not exist. Moreover, E4i=E2 when Φ=1 and E4i=E2 when Hiw,z=1.
(2) IfΦ>1 and Hiw,z>1, then E4i is locally asymptotically stable.
Here(43)Hiw,z=cβe-kτhTibg+βe-kτh;i=l,h.


## 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:(44)dxdt=rxt1-TtTm-β1-u1tvtxtTt,dydt=βe-kτ1-u1tvt-τxt-τTt-τ-ayt-pytzt,dvdt=1-u2taNyt-δvt-qvtwt,dwdt=gvtwt-hwt,dzdt=cytzt-bzt.The optimization problem that we consider is to maximize the following objective functional:(45)Ju1,u2=∫0tfxt+zt+wt-A12u12t+A22u22tdt,where tf stands for the time period of treatment. The two positive constants A1 and A2 are the weight for the treatment. It is legitimate to assume that two control functions, u1(t) and u2(t), 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 (u1∗,u2∗) such that(46)Ju1∗,u2∗=max⁡Ju1,u2:u1,u2∈U,where U is the control set defined by(47)U=u1t,u2t:uitmeasurable,0≤uit≤1,t∈0,tf,i=1,2.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(u1∗,u2∗)∈U such that(48)Ju1∗,u2∗=maxu1,u2∈U⁡Ju1,u2.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 H≡H(t,x,y,v,z,w,xτ,vτ,u1,u2,λi) point-wisely with respect to u1 and u2 as follows:(49)H=A12u1t2+A22u2t2-xt-zt-wt+∑i=05λifiwith(50)f1=rxt1-TtTm-β1-u1tvtxtTt,f2=βe-kτ1-u1tvt-τxt-τTt-τ-ayt-pytzt,f3=1-u2taNyt-δvt-qvtwt,f4=gvtwt-hwt,f5=cytzt-bzt.And λi,i=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.Theorem 11.
Given optimal controlsu1∗,u2∗ and solutions x∗, y∗, v∗, z∗, and w∗ of the corresponding state system (3), there exist adjoint variables, λ1, λ2, λ3, λ4, and λ5 satisfying the equations(51)dλ1tdt=1-λ1tr1-T∗tTm-rx∗tTm-1-u1∗tβv∗ty∗tT∗2-χ0,tf-τtλ2t+τu1∗t+τ-1βe-kτv∗ty∗tT∗2t,dλ2tdt=λ1trx∗tTm-1-u1∗tβv∗tx∗tT∗2+λ2ta+pz-λ3t1-u2∗taN-cz∗tλ5t-χ0,tf-τtλ2t+τu1∗t+τ-1βe-kτv∗tx∗tT∗2t,dλ3tdt=λ1tβ1-u1∗tx∗tT∗t+λ3tδ+qwt-λ4tgw∗t+χ0,tf-τtλ2t+τβe-kτu1∗t+τ-1x∗tT∗t,dλ4tdt=1+λ3tqv∗t+λ4th-gv∗t,dλ5tdt=1+λ2tpy∗t+λ5tb-cy∗t,with χ being an indicator function and T∗(t)=x∗(t)+y∗(t) also the transversality conditions (52)λitf=0,i=1,…,5.Moreover, the optimal control is given by(53)u1∗=min⁡1,max⁡0,βA1λ2te-kτv∗t-τx∗t-τT∗t-τ-λ1tv∗tx∗tT∗u2∗=min⁡1,max⁡0,1A2λ3taNy∗t.


## 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 [36], used similarly by [37] and enhanced for delay differential equation system [38–40].We consider the step sizeh>0 and n,m∈N2 with τ=mh and tf-t0=nh. We take m knots to left of t0 and right of tf, to get the following partition: (54)Δ=t-m=-τ<⋯<t-1<t0=0<t1<⋯<tn=tf<tn+1<⋯<tn+m,which gives ti=t0+ih(-m≤i≤n+m). The state and the adjoint variables are x(t), y(t), v(t), w(t), z(t), λ1(t), λ2(t), λ3(t), λ4(t), and λ5(t) and the controls are u1(t), u2(t) in terms of nodal points xi, yi, vi, wi, zi, λ1i, λ2i, λ3i, λ4i, λ5i, u1i, and u2i. By combining the forward and backward difference approximation, we get Algorithm 1.Algorithm 1
Step  1:
fori=-m,…,0, do:
x
i
=
x
0, yi=y0, Ti=x0+y0, vi=v0, wi=w0, zi=z0, u1i=0, u2i=0.
end for
fori=n,…,n+m, do:
λ
1
i
=
0, λ2i=0, λ3i=0, λ4i=0, λ5i=0.
end for
Step  2:
fori=0,…,n-1, do:
x
i
+
1
=
x
i
+
h
r
x
i
1
-
T
i
T
m
-
β
1
-
u
1
i
v
i
x
i
T
i,
y
i
+
1
=
y
i
+
h
β
e
-
k
τ
1
-
u
1
i
v
i
-
m
x
i
-
m
T
i
-
m
-
a
y
i
-
p
y
i
z
i,
v
i
+
1
=
v
i
+
h
1
-
u
2
i
a
N
y
i
-
δ
v
i
-
q
v
i
w
i,
w
i
+
1
=
w
i
+
h
g
v
i
w
i
-
h
w
i,
z
i
+
1
=
z
i
+
h
c
y
i
z
i
-
b
z
i,
T
i
+
1
=
x
i
+
1
+
y
i
+
1,
λ
1
n
-
i
-
1
=
λ
1
n
-
i
-
h
1
-
λ
1
n
-
i
r
1
-
T
i
+
1
T
m
-
r
x
i
+
1
T
m
-
1
-
u
1
i
β
v
i
+
1
y
i
+
1
T
i
+
1
-
χ
0
,
t
f
-
τ
t
n
-
i
λ
2
n
-
i
+
m
u
1
i
+
m
-
1
β
e
-
k
τ
v
i
+
1
y
i
+
1
T
i
+
1
2,
λ
2
n
-
i
-
1
=
λ
2
n
-
i
-
h
λ
1
n
-
i
r
x
i
+
1
T
m
-
1
-
u
1
i
β
v
i
+
1
x
i
+
1
T
i
+
1
2
+
λ
2
n
-
i
a
+
p
z
i
+
1
-
λ
3
n
-
i
1
-
u
2
i
a
N
-
c
z
i
+
1
λ
5
n
-
i
-
χ
0
,
t
f
-
τ
t
n
-
i
λ
2
n
-
i
+
m
u
1
i
+
m
-
1
β
e
-
k
τ
v
i
+
1
x
i
+
1
T
i
+
1
2,
λ
3
n
-
i
-
1
=
λ
3
n
-
i
-
h
λ
1
n
-
i
1
-
u
1
i
β
x
i
+
1
T
i
+
1
+
λ
3
n
-
i
δ
+
q
w
i
+
1
-
λ
4
n
-
i
g
w
i
+
1
+
χ
0
,
t
f
-
τ
t
n
-
i
λ
2
n
-
i
+
m
u
1
i
+
m
-
1
β
e
-
k
τ
x
i
+
1
T
i
+
1,
λ
4
n
-
i
-
1
=
λ
4
n
-
i
-
h
1
+
q
λ
3
n
-
i
v
i
+
1
+
λ
4
n
-
i
h
-
g
v
i
+
1,
λ
5
n
-
i
-
1
=
λ
4
n
-
i
-
h
1
+
p
λ
2
n
-
i
y
i
+
1
+
λ
5
n
-
i
b
-
c
y
i
+
1,
R
1
i
+
1
=
β
A
1
λ
2
n
-
i
-
1
e
-
k
τ
v
i
-
m
+
1
x
i
-
m
+
1
T
i
-
m
+
1
-
λ
1
n
-
i
-
1
v
i
+
1
x
i
+
1
T
i
+
1
R
2
i
+
1
=
1
A
2
λ
3
n
-
i
-
1
a
N
y
i
+
1,
u
1
i
+
1
=
min
⁡
1
,
max
⁡
R
1
i
+
1
,
0,
u
2
i
+
1
=
min
⁡
1
,
max
⁡
R
2
i
+
1
,
0,
end for
Step  3:
fori=1,…,n, write
x
∗
t
i
=
x
i
,
y
∗
t
i
=
y
i
,
T
∗
t
i
=
T
i
,
v
∗
t
i
=
v
i
,
z
∗
t
i
=
z
i
,
w
∗
t
i
=
w
i
,
u
1
∗
t
i
=
u
1
i
,
u
2
∗
t
i
=
u
2
i
.
end forFor the simulation, we use the parameter values given in Table1.Table 1
Parameters, their symbols, and default values used in Model (3).
Parameters
Meaning
Value
References
r
Maximum hepatocyte growth rate
≤1.0 day−1
[23, 24]
Tm
Hepatocyte carrying capacity
2 × 1011 cells
[25]
β
Rate of virion infection of hepatocytes
3.6 × 10−5–1.8 × 10−3 cells virion−1 day−1
[26]
τ
Time delay
1 day
[27, 28]
k
Normal death rate for hepatocytes
.0039 day−1
[27, 28]
a
Infected hepatocyte death rate
0.0693–0.00693 day−1
[25]
p
Clearance rate of infection
7 ± 1.7 × 10−4 ml/cell day−1
[29]
N
Number of free viruses produced by infected cells
480
[12, 13]
δ
Free virion half-life
0.67 day−1
[30]
q
Neutralization rate of virion and antibodies
10−10–10−12 ml day−1
[31]
g
Activation rate of B-cells
1.38 × 10−2–10−4 day−1
[31]
h
Death rate of B-cells
0.03–0.1 day−1
[31]
c
Activation rate of CTL cells
4.4 ± 1.5 × 10−7 ml cell−1 day−1
[29]
b
Death rate of CTL cells
0.5 day−1
[32]As the parameters having been chosen from different references (see Table1), we use in our numerical simulations a set of parameters that are within the range of the estimation of these references; that is, r=1, Tm=2×1011, β=0.0018, k=1.1×10-2, τ=1, a=0.0693, N=480, δ=0.693, q=0.01, p=0.001, c=4.4×10-8, b=0.5, q=10-10, g=10-4, h=0.1, A1=250, and A2=2500.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. Figure1 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.Figure 1
The uninfected cells (a). The infected cells (b). The HBV (c). The antibody response (d).
(a)
(b)
(c)
(d)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.Figure2 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.Figure 2
The uninfected cells as a function of time (a). The infected cells as a function of time (b).
(a)
(b)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 at2.5482, while the number of infected cells without control reaches 2.264 × 105, which makes the drug therapy efficiency in blocking the new infections at 98.73%.In Figure2, 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 × 106, which represents a perfect efficiency of the drug therapy in inhibiting the viral production (about 99.99%).Figure3(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.Figure 3
The HBV as a function of time (a). The antibody response as a function of time (b).
(a)
(b)The optimal therapy protocol is represented by Figure4. 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.Figure 4
The optimal controlu1(a) and the optimal controlu2(b) versus time.
(a)
(b)


## 6. 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.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 to10%for INF and to50% for NAs, and gradually to15%.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].



---

*Source: 6710575-2018-02-01.xml*