# Analysis of a New Delayed HBV Model with Exposed State and Immune Response to Infected Cells and Viruses

Authors: Deshun Sun; Fei Liu
Journal: BioMed Research International
 (2017)

---

## Abstract

We propose a comprehensive delayed HBV model, which not only considers the immune response to both infected cells and viruses and a time delay for the immune system to clear viruses but also incorporates an exposed state and the proliferation of hepatocytes. We prove the positivity and boundedness of solutions and analyze the global stability of two boundary equilibria and then study the local asymptotic stability and Hopf bifurcation of the positive (infection) equilibrium and also the stability of the bifurcating periodic solutions. Moreover, we illustrate how the factors such as the time delay, the immune response to infected cells and viruses, and the proliferation of hepatocytes affect the dynamics of the model by numerical simulation.


---

## Body

## 1. Introduction

Hepatitis B virus (HBV) has become one of the serious infectious diseases threatening global human health, which can cause chronic liver infection and further result in liver inflammation, fibrosis, cirrhosis, or even cancer [1]. Each year more than 1 million people die of end-stage liver diseases like cancer due to the HBV infection [2].Mathematical modeling and analysis of the dynamics of such infectious viruses as HBV play important roles in understanding the factors that govern the infectious disease progression and offering insights into developing treatment strategies and guiding antiviral drug therapies [3]. So far, there have been plenty of mathematical models proposed to describe and analyze virus infection, immune responses, and antiretroviral treatment [4–10].Among these works, the development of virus models with immune responses is gaining much attention [3, 11, 12]. The immune system is essential in controlling the level of virus reproduction in terms of the strength of the Cytotoxic T Lymphocyte (CTL) response. A small change of the CTL response may have a large effect on virus production and infected cells load. As to this aspect, typical work can be summarized as follows. Chen et al. [12] indicated that the immunity system can not only clear free viruses but also kill infected cells. Elaiw and AlShamrani [3] proposed a four-dimensional model with humoral immunity response and general function and analyzed the global asymptotic stability of all equilibria based on the general function. However, both models in [3, 12] do not consider time delay. In order to characterize the time of a body’s immune response after the virus infection of target cells, time delay has been taken into account [13–16]. For example, Zhu et al. [16] proposed an HIV infection model with CTL response delay and analyzed the effect of time delay on the stability of equilibria. Besides, a latent period would be necessary to be incorporated into a virus model because when viruses infect a healthy organ like liver, it will not be pathogenetic at once, as it takes about six weeks to six months from the infection to the incidence [17–20]. For example, Medley et al. [17] proposed an HBV model with an exposed state, namely, infected but not yet infectious. Moreover, it was modeled in [21, 22] that the liver can regenerate cells and compensate the lost infected hepatocytes by the proliferation of hepatocytes.In this paper, we will propose a more comprehensive model than those existing ones, which not only considers the immune response to both infected cells and viruses and a time delay for the immune system to clear viruses but also incorporates an exposed state and the proliferation of hepatocytes. We first discuss the existence of two boundary equilibria and one positive (infection) equilibrium. We then analyze the global stability of the two boundary equilibria, the local asymptotic stability and Hopf bifurcation of the positive equilibrium and also the stability of the bifurcating periodic solutions. Moreover, we perform numerical simulations to illustrate some of the theoretical results we obtain and also illustrate how the factors such as the immune response to infected cells and viruses and the proliferation of hepatocytes affect the dynamics of the model under time delay.The paper is structured as follows. In Section2, a delayed mathematical model is proposed, and the positivity and boundedness of solutions, existence of two boundary equilibria, and one positive equilibrium are discussed, followed by the global stability analysis of these two boundary equilibria and the local asymptotic stability and Hopf bifurcation of the positive equilibrium in Section 3. The stability of the bifurcating periodic solutions is studied in Section 4. In Section 5, some numerical simulations and discussions are given. Finally, a conclusion is given in Section 6.


## 2. Mathematical Model

Wang et al. [23] proposed a virus infection model of four-dimensional equations with delayed humoral immune response, which, however, does not involve an exposed state and consider the proliferation of hepatocytes. Although it considers the immune response to viruses, it does not involve the immune response to infected cells.Based on this model, we propose a new and comprehensive HBV model, which not only considers the immune response to both infected cells and viruses and a time delay for the immune system to clear viruses but also incorporates an exposed state and the proliferation of hepatocytes. To better understand our model, we illustrate its mechanism in Figure1.Figure 1
The mechanism of our model.The model is then given as follows:(1)x˙t=λ+rxt-dxt-βxtvt,e˙t=βxtvt-a1et-a2et,y˙t=a2et-a1yt-k1ytzt,v˙t=kyt-εvt-k2vtzt,z˙t=k3vt-τzt-τ-k4zt,w˙t=k5+k1k6ytzt-k7wt,where x,  e,  y,  v,  z, and w denote the number of uninfected cells, exposed cells, infected cells, free viruses, CTLs, and alanine aminotransferases (ALT), respectively. The parameter λ represents the natural production rate of uninfected cells. rx is a new term which is introduced to represent the proliferation of hepatocytes, where r is the proliferation rate. Parameters d, (and the following) a1,  ε,  k4, and k7 represent the natural death rate of uninfected cells, exposed cells, infected cells, free viruses, CTLs, and ALT, respectively. β represents the infection rate from uninfected cells to exposed cells and a2 the transfer rate from exposed cells to infected cells. The production rate of free viruses from infected cells is denoted by k, and the production rate of CTLs by k3. k5 represents the production rate of ALT from the extrahepatic tissue and k1k6 the production rate of ALT when the infected hepatocytes are killed by CTL. The immunity-induced clearance for infected cells is modeled by a term k1yz, where k1 represents the clearance rate of infected cells. Similarly, the immunity-induced clearance for free viruses is modeled by k2vz, where k2 represents the clearance rate of free viruses. τ is time delay. All the parameters in this paper are positive and d>r. For convenience, we define new parameter ρ=d-r.


### 2.1. Positivity and Boundedness of Solutions

In this subsection, we prove the positivity and the boundedness of solutions of system (1).We denotex(0)≥0,  e(0)≥0,  y(0)≥0,  w(0)≥0,  v(t)≥0,  z(t)≥0,  t∈[-τ,0]. From the first equation of system (1), we have x(t)=e-∫0t(ρ+βv(s))dsx(0)+λ∫0te-∫st[ρ+βv(ξ)]dξds; therefore, x(t)≥0 for ∀t>0 if x(0)≥0. Next, we consider the second, third, and fourth equation in system (1) as a nonautonomous system for e(t),y(t),v(t):(2)e˙=βxv-a1+a2et,y˙=a2e-a1+k1zyt,v˙=ky-ε+k2zvt.Based on Theorem2.1 in [24], we have e(t)≥0,  y(t)≥0,  v(t)≥0 if e(0)≥0,  y(0)≥0,  z(0)≥0.z
(
t
)
=
e
-
k
4
t
y
(
0
)
+
∫
0
t
k
3
v
(
t
-
γ
)
z
(
t
-
γ
)
e
-
k
4
(
t
-
γ
)
d
γ, we have z(t)≥0,  ∀t>0 if z(t)≥0,  t∈[-τ,0].w
(
t
)
=
e
-
k
7
t
w
(
0
)
+
e
-
k
7
t
∫
0
t
[
(
k
5
+
k
1
k
6
y
(
s
)
z
(
s
)
]
e
-
k
7
t
d
s, because y(t),z(t)≥0, so we have w(t)≥0,  ∀t>0 if w(0)≥0.Hence, the nonnegative is proved. In what follows, we will study the boundedness of solutions. We defineG(t) as a linear combination of x,  e,  y,  v,  z:(3)Gt=xt+et+yt+a12kvt+a1k22kk3zt+τ,δ=mind-r,a12,a1,ε,k4,dGtdt=λ-ρx-βxv+βxv-a1e-a2e+a2e-a1y-k1yz+a12kky-εv-k2vz+a1k22kk3k3vtzt-k4zt+τ=λ-d-rx-a12e-a1y-a1ε2kv-a1k2k42kk3zt+τ-k1yz≤λ-δGt.Therefore, we obtainlimt→∞⁡G(t)≤λ/δ, namely, x(t)+e(t)+y(t)+a1/2kv(t)+a1k2/2kk3z(t+τ)≤λ/δ. So we have 0≤x(t),e(t),y(t),v(t),z(t)≤λ/δ Because the boundedness of x(t),e(t),y(t),v(t),z(t),  limt→∞⁡w(t)≤(k5+k1k6)λ2/k7δ2.The boundedness is proved.


### 2.2. Equilibrium

In this subsection, we study the equilibria of system (1). The method to obtain equilibria is setting x˙=e˙=y˙=v˙=z˙=w˙=0 and computes the following:(4)λ-ρx-βxv=0,βxv-a1e-a2e=0,a2e-a1y-k1yz=0,ky-εv-k2vz=0,k3vt-τzt-τ-k4z=0,k5+k1k6yz-k7w=0.The system (1) has two boundary equilibria (an infection-free equilibrium E00 in which x≠0,  w≠0,  e=y=v=z=0 and an equilibrium without immune response E11 in which x≠0,  e≠0,  y≠0,  v≠0,  w≠0,  z=0) and a positive (infection) equilibrium E22 in which x≠0,  E≠0,  y≠0,  v≠0,  z≠0,  w≠0.The infection-free equilibrium isE00=(x0,0,0,0,0,w0), where x0=λ/ρ, and w0=k5/k7, and the basic reproductive number is obtained by the following method.Based on integral operator spectral radius, the basic reproductive number isR0=ρ(FV-1), where(5)F=00βx0a2000k0,V=a1+a2000a1+k1z0000ε.Hence, we have the basic reproductive number being R0=a2kx0β/a1ε(a1+a2).The equilibrium without immune response isE11=(x1,e1,y1,v1,0,w1), where x1=a1ε(a1+a2)/a2kβ,  e1=a1ε/a2kv1,  y1=ε/kv1,  v1=a2kλ/a1ε(a1+a2)-ρ/β, and w1=k5/k7. Similarly, we have the basic reproductive number is R1=k3v1/k4+kk3y1k1/a1k2k4 at E11=(x1,e1,y1,v1,0,w1).The infected positive equilibrium isE22=(x2,e2,y2,v2,z2,w2), where(6)v2=k4k3,x2=λρ+βv2,e2=βx2v2a1+a2,y2=A+A2+4B2,A=a1k2k4-k1k3εv2kk1k3,B=a2k2k4kk1k3>0,z2=k3ky2-εv2k2k4,w2=k5+k1k6y2z2k7.


## 2.1. Positivity and Boundedness of Solutions

In this subsection, we prove the positivity and the boundedness of solutions of system (1).We denotex(0)≥0,  e(0)≥0,  y(0)≥0,  w(0)≥0,  v(t)≥0,  z(t)≥0,  t∈[-τ,0]. From the first equation of system (1), we have x(t)=e-∫0t(ρ+βv(s))dsx(0)+λ∫0te-∫st[ρ+βv(ξ)]dξds; therefore, x(t)≥0 for ∀t>0 if x(0)≥0. Next, we consider the second, third, and fourth equation in system (1) as a nonautonomous system for e(t),y(t),v(t):(2)e˙=βxv-a1+a2et,y˙=a2e-a1+k1zyt,v˙=ky-ε+k2zvt.Based on Theorem2.1 in [24], we have e(t)≥0,  y(t)≥0,  v(t)≥0 if e(0)≥0,  y(0)≥0,  z(0)≥0.z
(
t
)
=
e
-
k
4
t
y
(
0
)
+
∫
0
t
k
3
v
(
t
-
γ
)
z
(
t
-
γ
)
e
-
k
4
(
t
-
γ
)
d
γ, we have z(t)≥0,  ∀t>0 if z(t)≥0,  t∈[-τ,0].w
(
t
)
=
e
-
k
7
t
w
(
0
)
+
e
-
k
7
t
∫
0
t
[
(
k
5
+
k
1
k
6
y
(
s
)
z
(
s
)
]
e
-
k
7
t
d
s, because y(t),z(t)≥0, so we have w(t)≥0,  ∀t>0 if w(0)≥0.Hence, the nonnegative is proved. In what follows, we will study the boundedness of solutions. We defineG(t) as a linear combination of x,  e,  y,  v,  z:(3)Gt=xt+et+yt+a12kvt+a1k22kk3zt+τ,δ=mind-r,a12,a1,ε,k4,dGtdt=λ-ρx-βxv+βxv-a1e-a2e+a2e-a1y-k1yz+a12kky-εv-k2vz+a1k22kk3k3vtzt-k4zt+τ=λ-d-rx-a12e-a1y-a1ε2kv-a1k2k42kk3zt+τ-k1yz≤λ-δGt.Therefore, we obtainlimt→∞⁡G(t)≤λ/δ, namely, x(t)+e(t)+y(t)+a1/2kv(t)+a1k2/2kk3z(t+τ)≤λ/δ. So we have 0≤x(t),e(t),y(t),v(t),z(t)≤λ/δ Because the boundedness of x(t),e(t),y(t),v(t),z(t),  limt→∞⁡w(t)≤(k5+k1k6)λ2/k7δ2.The boundedness is proved.


## 2.2. Equilibrium

In this subsection, we study the equilibria of system (1). The method to obtain equilibria is setting x˙=e˙=y˙=v˙=z˙=w˙=0 and computes the following:(4)λ-ρx-βxv=0,βxv-a1e-a2e=0,a2e-a1y-k1yz=0,ky-εv-k2vz=0,k3vt-τzt-τ-k4z=0,k5+k1k6yz-k7w=0.The system (1) has two boundary equilibria (an infection-free equilibrium E00 in which x≠0,  w≠0,  e=y=v=z=0 and an equilibrium without immune response E11 in which x≠0,  e≠0,  y≠0,  v≠0,  w≠0,  z=0) and a positive (infection) equilibrium E22 in which x≠0,  E≠0,  y≠0,  v≠0,  z≠0,  w≠0.The infection-free equilibrium isE00=(x0,0,0,0,0,w0), where x0=λ/ρ, and w0=k5/k7, and the basic reproductive number is obtained by the following method.Based on integral operator spectral radius, the basic reproductive number isR0=ρ(FV-1), where(5)F=00βx0a2000k0,V=a1+a2000a1+k1z0000ε.Hence, we have the basic reproductive number being R0=a2kx0β/a1ε(a1+a2).The equilibrium without immune response isE11=(x1,e1,y1,v1,0,w1), where x1=a1ε(a1+a2)/a2kβ,  e1=a1ε/a2kv1,  y1=ε/kv1,  v1=a2kλ/a1ε(a1+a2)-ρ/β, and w1=k5/k7. Similarly, we have the basic reproductive number is R1=k3v1/k4+kk3y1k1/a1k2k4 at E11=(x1,e1,y1,v1,0,w1).The infected positive equilibrium isE22=(x2,e2,y2,v2,z2,w2), where(6)v2=k4k3,x2=λρ+βv2,e2=βx2v2a1+a2,y2=A+A2+4B2,A=a1k2k4-k1k3εv2kk1k3,B=a2k2k4kk1k3>0,z2=k3ky2-εv2k2k4,w2=k5+k1k6y2z2k7.


## 3. Analysis

### 3.1. Global Stability Analysis of the Two Boundary Equilibria

In this section, we will employ the direct Lyapunov method and LaSalle’s invariance principle to establish the global asymptotic stability of the two boundary equilibria.Theorem 1.
The infection-free equilibriumE00 is globally asymptotically stable if and only if R0<1.
See AppendixA for proof.Theorem 2.
The equilibrium without immune responseE11 is globally asymptotically stable if and only if R1<1.
See AppendixB for proof.


### 3.2. Local Asymptotic Stability and Hopf Bifurcation of the Positive Equilibrium

In this section, we will discuss the local asymptotic stability and Hopf bifurcation of the positive equilibriumE22.The characteristic equation of system (1) at E22 is as follows:(7)Hλ;τ=λ+ρ+βv200βx200-βv2λ+a1+a20-βx2000-a2λ+a1+k1z20k1y2000-kλ+ε+k2z2k2v20000-k3z2e-λτλ-k3v2e-λτ+k4000-k1k6z20-k1k6y2λ+k7=0.Define(8)A1=2a1+a2+kz2,A2=a1+a2a1+kz2,A3=ε+k4-k2z2,A4=k4ε-k2z2,A5=a2kβx2,A6=a2kk4βx2,M1=-k3v2,M2=2k2z2-εk3v2,M3=kk1k3y2z2,M4=-a2kk3βx2v2+kk1k3y2z2a1+a2,s1=A1+A3,s2=A2+A4+A1A3,s3=A5+A1A4+A2A3,s4=A6+A2A4,s5=M1,s6=M2+A1M1,s7=A1M2+A2M1+M3,s8=A2M2+M4,s9=ρ+βv2,s10=β2x2v2a2k,s11=β2x2v2a2kk4,s12=-β2x2v22a2kk3,B1=s1+s9,B2=s2+s1s9,B3=s3+s2s9,B4=s4+s3s9+s10,B5=s11+s4s9,B6=s5,B7=s6+s5s9,B8=s7+s6s9,B9=s8+s7s9,B10=s12+s8s9,D1=B1+k7,D2=B2+B1k7,D3=B3+B2k7,D4=B4+B3k7,D5=B5+B4k7,D6=B6+B5k7,S1=B6,S2=B7+B6k7,S3=B8+B7k7,S4=B9+B8k7,S5=B10+B9k7,S6=B10k7.Then the characteristic equationH(λ;τ) above becomes(9)Hλ;τ=λ6+D1λ5+D2λ4+D3λ3+D4λ2+D5λ+D6+S1λ5e-λτ+S2λ4e-λτ+S3λ3e-λτ+S4λ2e-λτ+S5λe-λτ+S6e-λτ=0.When τ=0, (9) further becomes(10)Hλ;τ=λ6+n1λ5+n2λ4+n3λ3+n4λ2+n5λ+n6=0,where(11)n1=D1+S1,n2=D2+S2,n3=D3+S3,n4=D4+S4,n5=D5+S5,n6=D6+S6.Using the Routh-Hurwitz criterion [23], we obtain the following lemma.Lemma 3.
If (10) satisfies Δ1≡n1>0,  Δ2≡n11n3n2>0,  Δ3≡n110n3n2n1n5n4n3>0 and Δ4=n1100n3n2n11n5n4n3n20n6n5n4>0, then the positive equilibrium E22 is locally asymptotically stable when τ=0.Proof.
By the Routh-Hurwitz criterion, if the four conditions are satisfied, then all roots of (10) have negative real parts. Therefore, the positive equilibrium E22 is locally asymptotically stable when τ=0. For more details, we refer the readers to [25, 26].
From Lemma3, we know that all roots of H(λ;τ) lie to the left of the imaginary axis when τ=0. However, with τ increasing from zero, some of its roots may cross the imaginary axis to the right. In this case, there are some roots having positive real parts, and therefore the equilibrium E22 becomes unstable. Next, we will discuss the stability of system (1) at E22 when τ>0.
We first divide (9) into two parts and obtain(12)λ6+D1λ5+D2λ4+D3λ3+D4λ2+D5λ+D62=S1λ5+S2λ4+S3λ3+S4λ2+S5λ+S62e-λτ2.
Suppose (9) has a purely imaginary root λ=iω(ω>0). Substituting λ=iω into (12) yields(13)-ω6+D1ω5i+D2ω4-D3ω3i-D4ω2+D5ωi+D62=S1ω5i+S2ω4-S3ω3i-S4ω2+S5ωi+S62.By separating the real part and imaginary part, the following real part is obtained:(14)ω12+C1ω10+C2ω8+C3ω6+C4ω4+C5ω2+C6=0,where(15)C1=-D12-2D2+S12,C2=D22+2D4+2D1D3-S22-2S1S3,C3=-D32-2D6-2D1D5-2D2D4+S32+2S1S5+2S2S4,C4=D42+2D2D6+2D3D5-S42-2S2S6-2S3S5,C5=-D52-2D4D6+S52+2S4S6,C6=D62-S62.Let(16)Gx=x6+C1x5+C2x4+C3x3+C4x2+C5x+C6.Therefore, if (9) has a purely imaginary root iω, it is equivalent to the fact that G(x)=0 has a positive real root ω2.Theorem 4.
IfG(x)=0 has no positive real roots, then the positive equilibrium E22 is locally asymptotically stable for any τ>0.Proof.
IfG(x)=0 has no positive real roots, then obviously (9) has no positive real roots. Therefore, the positive equilibrium E22 is locally asymptotically stable for any τ>0.
Substitutingλ=iω into (22), we obtain the real part,(17)-ω6+D2ω4-D4ω2+D6+S1ω5-S3ω3+S5ωsin⁡ωτ+S2ω4-S4ω2+S6cos⁡ωτ=0,and imaginary part,(18)D1ω5-D3ω3+D5ω+-S2ω4+S4ω2-S6sin⁡ωτ+S1ω5-S3ω3+S5ωcos⁡ωτ=0.Assuming that G(x)=0 has n~(1≤n~≤6) positive real roots, denoted by xn(1≤n≤n~). As xn=ω, we then have(19)cos⁡xnτ=Qn=-S2xn2+S4xn-S6-xn3+D2xn2-D4xn+D6S2xn2-S4xn+S62+S1xn2xn-S3xnxn+S5xn2-S1xn2xn-S3xnxn+S5xnD1xn2xn-D3xnxn+D5xnS2xn2-S4xn+S62+S1xn2xn-S3xnxn+S5xn2,sin⁡xnτ=Pn=S2xn2-S4xn+S6D1xn2xn-D3xnxn+D5xnS2xn2-S4xn+S62+S1xn2xn-S3xnxn+S5xn2-S1xn2xn-S3xnxn+S5xn-xn3+D2xn2-D4xn+D6S2xn2-S4xn+S62+S1xn2xn-S3xnxn+S5xn2.Let(20)τnj=1xnarccos⁡Qn+2jπ,ifPn≥0,1xn2π-arccos⁡Qn+2jπ,ifPn<0,where 1≤n≤n~ and j=0,1,2,….
Therefore, the characteristic equationH(λ;τn(j))=0 has a pair of purely imaginary roots ±ixn. For every integer j and 1≤n≤n~, define λn(j)(τ)=αn(j)(τ)+iωn(j)(τ) as the root of (9) near τn(j), satisfying αn(j)(τn(j))=0 and ωn(j)(τn(j))=xn. Then the following theorem is obtained.Theorem 5.
IfG(x)=0 has some positive real roots, then E22 is locally asymptotically stable for τ∈[0,τn0(0)), where (21)τn00=min⁡τnj∣1≤n≤n~,j=0,1,2,….Proof.
Forτn0(0)=min⁡τn(j)∣1≤n≤n~,j=0,1,2,…,  G(x)=0 has no positive real roots when τ∈[0,τn0(0)), which means that all the roots of (9) have strictly negative real parts when τ∈[0,τn0(0)). Therefore, E22 is locally asymptotically stable for τ∈[0,τn0(0)).Theorem 6.
Ifxn0 is a simple root of G(x)=0, then there is a Hopf bifurcation for the system as τ increases past τn0(0).
See AppendixC for proof.


## 3.1. Global Stability Analysis of the Two Boundary Equilibria

In this section, we will employ the direct Lyapunov method and LaSalle’s invariance principle to establish the global asymptotic stability of the two boundary equilibria.Theorem 1.
The infection-free equilibriumE00 is globally asymptotically stable if and only if R0<1.
See AppendixA for proof.Theorem 2.
The equilibrium without immune responseE11 is globally asymptotically stable if and only if R1<1.
See AppendixB for proof.


## 3.2. Local Asymptotic Stability and Hopf Bifurcation of the Positive Equilibrium

In this section, we will discuss the local asymptotic stability and Hopf bifurcation of the positive equilibriumE22.The characteristic equation of system (1) at E22 is as follows:(7)Hλ;τ=λ+ρ+βv200βx200-βv2λ+a1+a20-βx2000-a2λ+a1+k1z20k1y2000-kλ+ε+k2z2k2v20000-k3z2e-λτλ-k3v2e-λτ+k4000-k1k6z20-k1k6y2λ+k7=0.Define(8)A1=2a1+a2+kz2,A2=a1+a2a1+kz2,A3=ε+k4-k2z2,A4=k4ε-k2z2,A5=a2kβx2,A6=a2kk4βx2,M1=-k3v2,M2=2k2z2-εk3v2,M3=kk1k3y2z2,M4=-a2kk3βx2v2+kk1k3y2z2a1+a2,s1=A1+A3,s2=A2+A4+A1A3,s3=A5+A1A4+A2A3,s4=A6+A2A4,s5=M1,s6=M2+A1M1,s7=A1M2+A2M1+M3,s8=A2M2+M4,s9=ρ+βv2,s10=β2x2v2a2k,s11=β2x2v2a2kk4,s12=-β2x2v22a2kk3,B1=s1+s9,B2=s2+s1s9,B3=s3+s2s9,B4=s4+s3s9+s10,B5=s11+s4s9,B6=s5,B7=s6+s5s9,B8=s7+s6s9,B9=s8+s7s9,B10=s12+s8s9,D1=B1+k7,D2=B2+B1k7,D3=B3+B2k7,D4=B4+B3k7,D5=B5+B4k7,D6=B6+B5k7,S1=B6,S2=B7+B6k7,S3=B8+B7k7,S4=B9+B8k7,S5=B10+B9k7,S6=B10k7.Then the characteristic equationH(λ;τ) above becomes(9)Hλ;τ=λ6+D1λ5+D2λ4+D3λ3+D4λ2+D5λ+D6+S1λ5e-λτ+S2λ4e-λτ+S3λ3e-λτ+S4λ2e-λτ+S5λe-λτ+S6e-λτ=0.When τ=0, (9) further becomes(10)Hλ;τ=λ6+n1λ5+n2λ4+n3λ3+n4λ2+n5λ+n6=0,where(11)n1=D1+S1,n2=D2+S2,n3=D3+S3,n4=D4+S4,n5=D5+S5,n6=D6+S6.Using the Routh-Hurwitz criterion [23], we obtain the following lemma.Lemma 3.
If (10) satisfies Δ1≡n1>0,  Δ2≡n11n3n2>0,  Δ3≡n110n3n2n1n5n4n3>0 and Δ4=n1100n3n2n11n5n4n3n20n6n5n4>0, then the positive equilibrium E22 is locally asymptotically stable when τ=0.Proof.
By the Routh-Hurwitz criterion, if the four conditions are satisfied, then all roots of (10) have negative real parts. Therefore, the positive equilibrium E22 is locally asymptotically stable when τ=0. For more details, we refer the readers to [25, 26].
From Lemma3, we know that all roots of H(λ;τ) lie to the left of the imaginary axis when τ=0. However, with τ increasing from zero, some of its roots may cross the imaginary axis to the right. In this case, there are some roots having positive real parts, and therefore the equilibrium E22 becomes unstable. Next, we will discuss the stability of system (1) at E22 when τ>0.
We first divide (9) into two parts and obtain(12)λ6+D1λ5+D2λ4+D3λ3+D4λ2+D5λ+D62=S1λ5+S2λ4+S3λ3+S4λ2+S5λ+S62e-λτ2.
Suppose (9) has a purely imaginary root λ=iω(ω>0). Substituting λ=iω into (12) yields(13)-ω6+D1ω5i+D2ω4-D3ω3i-D4ω2+D5ωi+D62=S1ω5i+S2ω4-S3ω3i-S4ω2+S5ωi+S62.By separating the real part and imaginary part, the following real part is obtained:(14)ω12+C1ω10+C2ω8+C3ω6+C4ω4+C5ω2+C6=0,where(15)C1=-D12-2D2+S12,C2=D22+2D4+2D1D3-S22-2S1S3,C3=-D32-2D6-2D1D5-2D2D4+S32+2S1S5+2S2S4,C4=D42+2D2D6+2D3D5-S42-2S2S6-2S3S5,C5=-D52-2D4D6+S52+2S4S6,C6=D62-S62.Let(16)Gx=x6+C1x5+C2x4+C3x3+C4x2+C5x+C6.Therefore, if (9) has a purely imaginary root iω, it is equivalent to the fact that G(x)=0 has a positive real root ω2.Theorem 4.
IfG(x)=0 has no positive real roots, then the positive equilibrium E22 is locally asymptotically stable for any τ>0.Proof.
IfG(x)=0 has no positive real roots, then obviously (9) has no positive real roots. Therefore, the positive equilibrium E22 is locally asymptotically stable for any τ>0.
Substitutingλ=iω into (22), we obtain the real part,(17)-ω6+D2ω4-D4ω2+D6+S1ω5-S3ω3+S5ωsin⁡ωτ+S2ω4-S4ω2+S6cos⁡ωτ=0,and imaginary part,(18)D1ω5-D3ω3+D5ω+-S2ω4+S4ω2-S6sin⁡ωτ+S1ω5-S3ω3+S5ωcos⁡ωτ=0.Assuming that G(x)=0 has n~(1≤n~≤6) positive real roots, denoted by xn(1≤n≤n~). As xn=ω, we then have(19)cos⁡xnτ=Qn=-S2xn2+S4xn-S6-xn3+D2xn2-D4xn+D6S2xn2-S4xn+S62+S1xn2xn-S3xnxn+S5xn2-S1xn2xn-S3xnxn+S5xnD1xn2xn-D3xnxn+D5xnS2xn2-S4xn+S62+S1xn2xn-S3xnxn+S5xn2,sin⁡xnτ=Pn=S2xn2-S4xn+S6D1xn2xn-D3xnxn+D5xnS2xn2-S4xn+S62+S1xn2xn-S3xnxn+S5xn2-S1xn2xn-S3xnxn+S5xn-xn3+D2xn2-D4xn+D6S2xn2-S4xn+S62+S1xn2xn-S3xnxn+S5xn2.Let(20)τnj=1xnarccos⁡Qn+2jπ,ifPn≥0,1xn2π-arccos⁡Qn+2jπ,ifPn<0,where 1≤n≤n~ and j=0,1,2,….
Therefore, the characteristic equationH(λ;τn(j))=0 has a pair of purely imaginary roots ±ixn. For every integer j and 1≤n≤n~, define λn(j)(τ)=αn(j)(τ)+iωn(j)(τ) as the root of (9) near τn(j), satisfying αn(j)(τn(j))=0 and ωn(j)(τn(j))=xn. Then the following theorem is obtained.Theorem 5.
IfG(x)=0 has some positive real roots, then E22 is locally asymptotically stable for τ∈[0,τn0(0)), where (21)τn00=min⁡τnj∣1≤n≤n~,j=0,1,2,….Proof.
Forτn0(0)=min⁡τn(j)∣1≤n≤n~,j=0,1,2,…,  G(x)=0 has no positive real roots when τ∈[0,τn0(0)), which means that all the roots of (9) have strictly negative real parts when τ∈[0,τn0(0)). Therefore, E22 is locally asymptotically stable for τ∈[0,τn0(0)).Theorem 6.
Ifxn0 is a simple root of G(x)=0, then there is a Hopf bifurcation for the system as τ increases past τn0(0).
See AppendixC for proof.


## 4. Stability of the Bifurcating Periodic Solutions

In this section, we will continue to derive the explicit formulas for determining the stability, direction, and other properties of the Hopf bifurcation at a critical valueτn0(0) by means of the normal form and the center manifold theory [27].First, we make the following hypotheses.(1) Equation (9) has a pair of purely imaginary roots ±iω0 at τ=τ0, where τ0∈τn(j)∣1≤n≤n~,j=0,1,2,….(2) The remaining roots of (22) have strictly negative real parts.(3)ω0 is a simple root of G(x)=0.We useu=τ-τ0 to represent a new bifurcation parameter. Let X(t)=(x-x2,e-e2,y-y2,v-v2,z-z2,w-w2)T, and Xt(θ)=X(t+θ), where θ∈-τ,0. Therefore, system (1) can be written as the following functional differential equation:(22)X˙t=LuXt+fXt·,u,where (23)Luϕ=F1ϕ0+F2ϕ-τ,F1=-ρ-βv200-βx200βv2-a1-a20βx2000a2-a1-k1z20-k1y2000k-ε-k2z2-k2v200000-k4000k1k6z20k1k6y2-k7,F2=000000000000000000000000000k3z2k3v20000000,fϕ,u=-βϕ10ϕ40βϕ10ϕ40-k1ϕ30ϕ50-k2ϕ40ϕ50k3ϕ4-τϕ5-τk1k6ϕ30ϕ50.By the Riesz representation theorem [28], there exists a 6×6 matrix-valued function such that(24)Luϕ=∫-τ0dηθ,uϕθ,where dη(θ,u)=F1δ(θ)dθ+F2δ(θ+τ)dθ.Forϕ∈C([-τ,0],R6), we further define(25)Auϕθ=dϕθdθ,ifθ∈-τ,0,∫-τ0dηξ,uϕξ≡Luϕ,ifθ=0,Ruϕθ=0,ifθ∈-τ,0,fϕ,u,ifθ=0.Then system (22) can be written as(26)X˙tθ=AuXtθ+RuXtθ.For φ∈C([-τ,0],R6), define(27)A∗0φs=dϕsds,ifs∈-τ,0,∫-τ0dηTξ,0φ-ξ,ifs=0,and an inner product of ϕ,φ(28)φ,ϕ=φ-T0ϕ0-∫θ=-τ0∫ξ=0θφ-Tξ-θdηθϕξdξ,where η(θ)=η(θ,0) and ϕ∈C([-τ,0],R6). Then A(0) and A∗(0) are adjoint operators.Leth(θ) and h∗(s) be the eigenvectors of A(0) and A∗(0) corresponding to the eigenvalues iω0 and -iω0, respectively. We choose h(θ) and h∗(s) as (29)hθ=1,h2,h3,h4,h5,h6Teiω0θ,h∗s=D1,h2∗,h3∗,h4∗,h5∗,h6∗Teiω0s,so that h∗(s),h(θ)=1 is satisfied. We give the detailed computation of (29) in Appendix D.In the following, we will compute the coefficients,g20,  g11,  g02, and g21, using the method given in [27]. The detailed computation of g20,  g11,  g02, and g21 is presented in Appendix E.Then the following values can be computed:(30)c10=i2ω0g11g20-2g112-g0223+g212,u2=-Re⁡c10Re⁡λ′τ0,β2=2Re⁡c10.The signs of u2,  β2 determine the direction of the Hopf bifurcation and the stability of bifurcating periodic solutions, respectively [27]. From (C.11) of Appendix C, we obtain signdαn(j)(τ)/dττ=τnj=sign(dG/dx)x=xn.Letu2∗=-Re(c1(0))/G′(ω02). We obtain the following theorem.Theorem 7.
Assume the hypotheses (1), (2), and (3) at the beginning of Section4 hold.
(1) Ifu2∗>0(u2∗<0), then the bifurcating periodic solutions exist for τ>τ0(τ<τ0) in a τ0-neighborhood.
(2) Ifβ2<0(β2>0), the bifurcating periodic solutions are orbitally asymptotically stable as t→+∞(t→-∞).Proof.
(1) Ifτ0=τn0(0), where τn0(0)=min⁡τn(j)∣1≤n≤n~,j=0,1,2,…, and the hypotheses (1), (2), and (3) hold, then, from Theorem 6, we draw the conclusion that the existence and stability of the bifurcating periodic solutions are only determined by Re(c1(0)).
(2) Ifβ2<0, namely, Re(c1(0))<0, then there exist stable periodic solutions for τ>τn0(0) in a τ0-neighborhood. So the bifurcating periodic solutions are orbitally asymptotically stable as t→+∞.


## 5. Simulation and Discussions

In this section, we will numerically illustrate the theoretical results obtained above and also discuss how the factors such as the immune response to infected cells and viruses and the proliferation of hepatocytes affect the dynamics of the model under time delay.For the following simulations, we choose the parameter values for system (1) as follows:(31)λ=4.0551,r=0.6933,d=4.4096,β=4.6178,a1=0.0638,a2=1.8858,k1=0.8391,k=2.7011,ε=0.5083,k2=0.1963,k3=4.6661,k4=4.8580,k5=1.8046,k6=3.2210,k7=0.3397.We set the initial values tox(t)=1,  e(t)=1,  y(t)=1,  v(t)=1,  z(t)=1, and w(t)=1, where t∈[-τ,0].


### 5.1. Hopf Bifurcation and the Stability of Periodic Solutions

With the parameter values given in (31), we have the positive equilibrium E22=(0.4757,1.1732,0.3281,1.0411,1.7470,9.8727) and the critical time value τn0(0)=0.041.Whenτ>0.041, we obtain stable bifurcating periodic solutions. For example, when τ=0.05, the simulation result is shown Figures 2 and 3. Figure 2 indicates that a stable limit cycle is obtained as expected and Figure 3 indicates the state dynamics of uninfected cells, exposed cells, infected cells, free viruses, CTLs, and ALT, which are periodically oscillating.Figure 2
The positive equilibriumE22 bifurcates into a periodic solution at τ=0.05 (a limit cycle).Figure 3
The positive equilibriumE22 bifurcates into a periodic solution at τ=0.05.
(a)
(b)
(c)
(d)
(e)
(f)Whenτ<0.041, the bifurcating periodic solutions are unstable. For example, when τ=0.03, the simulation result is shown in Figures 4 and 5. From Figures 4 and 5, we know that the positive equilibrium E22 is asymptotically stable and the system will converge to E22.Figure 4
The positive equilibriumE22 at τ=0.03.Figure 5
The positive equilibriumE22 remains stale at τ=0.03.
(a)
(b)
(c)
(d)
(e)
(f)With the increasing of time delay(τ), the radius of limit cycle will increase. The simulation result is shown Figure 6.Figure 6
The different periodic solutions atτ=0.05,τ=0.07 and τ=0.09 (three limit cycles).


### 5.2. The Immune Response to Infected Cells

Here we will investigate the effect of the immune response to infected cells on the model dynamics under time delay. When we changek1=0.8391 to k1=0.01 by fixing other values given in (31), we obtain a simulation result at τ=0.05, illustrated in Figure 7. Comparing Figure 2 when k1=0.8391 and Figure 7 when k1=0.01, we can see that, with the decrease of k1, the stable periodic solution becomes unstable, that is, asymptotically stable.Figure 7
The positive equilibriumE22 remains stable at τ=0.05 when k1=0.01 (blue line) and exists at stable periodic solutions when k1=0.8391 (red line).
(a)
(b)
(c)
(d)
(e)
(f)


### 5.3. The Immune Response to Viruses

We continue to investigate the effect of the immune response to viruses on the model dynamics under time delay. When we changek2=0.1963 to k2=0.001 by fixing other values given in (31), we obtain a simulation result at τ=0.05, illustrated in Figure 8. Similarly, comparing Figures 2 and 8, we can see that, with the decrease of k2, the stable periodic solution also becomes unstable, that is, asymptotically stable. We further can see that the effect of the immune response to infected cells on the model dynamics is similar to that of the immune response to viruses.Figure 8
The positive equilibriumE22 remains stable at τ=0.05 when k2=0.001 (blue line) and exists at stable periodic solutions when k2=0.1963 (red line).
(a)
(b)
(c)
(d)
(e)
(f)


### 5.4. Proliferation of Hepatocytes

Then we investigate the effect of proliferation of hepatocytes on the model dynamics under time delay. For this, we still keepτ=0.05 and change the value of parameter r. When we change r=0.6933 to r=0.001 by fixing other values given in (31), we obtain a simulation result at τ=0.05, illustrated in Figure 9. Figure 9 shows that when r=0.001, the bifurcating periodic solution is stable, compared with Figure 2 when r=0.6933. We also try other values of parameter rand obtaining similar results. Thus, we can see parameter r has a small effect on the model dynamics, which reflect in periodicity and the positive equilibrium E22.Figure 9
The positive equilibriumE22 remains stable at τ=0.05 when r=0.001 (blue line) and exists at stable periodic solutions when r=0.6933 (red line).
(a)
(b)
(c)
(d)
(e)
(f)


## 5.1. Hopf Bifurcation and the Stability of Periodic Solutions

With the parameter values given in (31), we have the positive equilibrium E22=(0.4757,1.1732,0.3281,1.0411,1.7470,9.8727) and the critical time value τn0(0)=0.041.Whenτ>0.041, we obtain stable bifurcating periodic solutions. For example, when τ=0.05, the simulation result is shown Figures 2 and 3. Figure 2 indicates that a stable limit cycle is obtained as expected and Figure 3 indicates the state dynamics of uninfected cells, exposed cells, infected cells, free viruses, CTLs, and ALT, which are periodically oscillating.Figure 2
The positive equilibriumE22 bifurcates into a periodic solution at τ=0.05 (a limit cycle).Figure 3
The positive equilibriumE22 bifurcates into a periodic solution at τ=0.05.
(a)
(b)
(c)
(d)
(e)
(f)Whenτ<0.041, the bifurcating periodic solutions are unstable. For example, when τ=0.03, the simulation result is shown in Figures 4 and 5. From Figures 4 and 5, we know that the positive equilibrium E22 is asymptotically stable and the system will converge to E22.Figure 4
The positive equilibriumE22 at τ=0.03.Figure 5
The positive equilibriumE22 remains stale at τ=0.03.
(a)
(b)
(c)
(d)
(e)
(f)With the increasing of time delay(τ), the radius of limit cycle will increase. The simulation result is shown Figure 6.Figure 6
The different periodic solutions atτ=0.05,τ=0.07 and τ=0.09 (three limit cycles).


## 5.2. The Immune Response to Infected Cells

Here we will investigate the effect of the immune response to infected cells on the model dynamics under time delay. When we changek1=0.8391 to k1=0.01 by fixing other values given in (31), we obtain a simulation result at τ=0.05, illustrated in Figure 7. Comparing Figure 2 when k1=0.8391 and Figure 7 when k1=0.01, we can see that, with the decrease of k1, the stable periodic solution becomes unstable, that is, asymptotically stable.Figure 7
The positive equilibriumE22 remains stable at τ=0.05 when k1=0.01 (blue line) and exists at stable periodic solutions when k1=0.8391 (red line).
(a)
(b)
(c)
(d)
(e)
(f)


## 5.3. The Immune Response to Viruses

We continue to investigate the effect of the immune response to viruses on the model dynamics under time delay. When we changek2=0.1963 to k2=0.001 by fixing other values given in (31), we obtain a simulation result at τ=0.05, illustrated in Figure 8. Similarly, comparing Figures 2 and 8, we can see that, with the decrease of k2, the stable periodic solution also becomes unstable, that is, asymptotically stable. We further can see that the effect of the immune response to infected cells on the model dynamics is similar to that of the immune response to viruses.Figure 8
The positive equilibriumE22 remains stable at τ=0.05 when k2=0.001 (blue line) and exists at stable periodic solutions when k2=0.1963 (red line).
(a)
(b)
(c)
(d)
(e)
(f)


## 5.4. Proliferation of Hepatocytes

Then we investigate the effect of proliferation of hepatocytes on the model dynamics under time delay. For this, we still keepτ=0.05 and change the value of parameter r. When we change r=0.6933 to r=0.001 by fixing other values given in (31), we obtain a simulation result at τ=0.05, illustrated in Figure 9. Figure 9 shows that when r=0.001, the bifurcating periodic solution is stable, compared with Figure 2 when r=0.6933. We also try other values of parameter rand obtaining similar results. Thus, we can see parameter r has a small effect on the model dynamics, which reflect in periodicity and the positive equilibrium E22.Figure 9
The positive equilibriumE22 remains stable at τ=0.05 when r=0.001 (blue line) and exists at stable periodic solutions when r=0.6933 (red line).
(a)
(b)
(c)
(d)
(e)
(f)


## 6. Conclusions

In this paper, we consider a comprehensive delayed HBV model. Different from other existing models, our model not only considers the immune response to both infected cells and viruses and a time delay for the immune system to clear viruses but also incorporates an exposed state and the proliferation of hepatocytes.We then prove the positivity and boundedness of solutions and analyze the global stability of two boundary equilibria and investigate the local asymptotic stability and Hopf bifurcation of the positive (infection) equilibrium and also the stability of the bifurcating periodic solutions. We also numerically illustrate the Hopf bifurcation and the stability of the bifurcating periodic solutions.Moreover, we numerically illustrate how the factors such as the time delay, the immune response to infected cells and viruses, and the proliferation of hepatocytes affect the dynamics of the model, which shows that the former two factors have a big effect on the model dynamics, while the latter one does not have a big effect.



---

*Source: 7805675-2017-11-16.xml*