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

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.

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 Section 2, 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 Figure 1.

Figure 1

The mechanism of our model.

The model is then given as follows: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) a 1,  ε,  k 4, and k 7 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 a 2 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 k 3. k 5 represents the production rate of ALT from the extrahepatic tissue and k 1 k 6 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 k 1 yz, where k 1 represents the clearance rate of infected cells. Similarly, the immunity-induced clearance for free viruses is modeled by k 2 vz, where k 2 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 denote x(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 −∫ x(0) + λ∫0 e −∫ 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):

Based on Theorem 2.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 − y(0) + ∫0 k 3 v(t − γ)z(t − γ)e − dγ, we have z(t) ≥ 0,  ∀t > 0 if z(t) ≥ 0,  t ∈ [−τ, 0].

w(t) = e − w(0) + e −∫0 [(k 5 + k 1 k 6 y(s)z(s)]e − ds, 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 define G(t) as a linear combination of x,  e,  y,  v,  z:

Therefore, we obtain lim⁡G(t) ≤ λ/δ, namely, x(t) + e(t) + y(t) + (a 1/2k)v(t) + (a 1 k 2/2kk 3)z(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),  lim⁡w(t) ≤ (k 5 + k 1 k 6)λ 2/k 7 δ 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 and computes the following:

The system (1) has two boundary equilibria (an infection-free equilibrium E 00 in which x ≠ 0,  w ≠ 0,  e = y = v = z = 0 and an equilibrium without immune response E 11 in which x ≠ 0,  e ≠ 0,  y ≠ 0,  v ≠ 0,  w ≠ 0,  z = 0) and a positive (infection) equilibrium E 22 in which x ≠ 0,  E ≠ 0,  y ≠ 0,  v ≠ 0,  z ≠ 0,  w ≠ 0.

The infection-free equilibrium is E 00 = (x 0, 0,0, 0,0, w 0), where x 0 = λ/ρ, and w 0 = k 5/k 7, and the basic reproductive number is obtained by the following method.

Based on integral operator spectral radius, the basic reproductive number is R 0 = ρ(FV −1), whereHence, we have the basic reproductive number being R 0 = a 2 kx 0 β/a 1 ε(a 1 + a 2).

The equilibrium without immune response is E 11 = (x 1, e 1, y 1, v 1, 0, w 1), where x 1 = a 1 ε(a 1 + a 2)/a 2 kβ,  e 1 = (a 1 ε/a 2 k)v 1,  y 1 = (ε/k)v 1,  v 1 = a 2 kλ/a 1 ε(a 1 + a 2) − ρ/β, and w 1 = k 5/k 7. Similarly, we have the basic reproductive number is R 1 = k 3 v 1/k 4 + kk 3 y 1 k 1/a 1 k 2 k 4 at E 11 = (x 1, e 1, y 1, v 1, 0, w 1).

The infected positive equilibrium is E 22 = (x 2, e 2, y 2, v 2, z 2, w 2), where

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 equilibrium E 00 is globally asymptotically stable if and only if R 0 < 1.

See Appendix A for proof.

Theorem 2 .

The equilibrium without immune response E 11 is globally asymptotically stable if and only if R 1 < 1.

See Appendix B 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 equilibrium E 22.

The characteristic equation of system (1) at E 22 is as follows:

Define

Then the characteristic equation H(λ; τ) above becomesWhen τ = 0, (9) further becomeswhereUsing the Routh-Hurwitz criterion [23], we obtain the following lemma.

Lemma 3 .

If (10) satisfies Δ1 ≡ n 1 > 0,  ,   and , then the positive equilibrium E 22 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 E 22 is locally asymptotically stable when τ = 0. For more details, we refer the readers to [25, 26].

From Lemma 3, 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 E 22 becomes unstable. Next, we will discuss the stability of system (1) at E 22 when τ > 0.

We first divide (9) into two parts and obtain

Suppose (9) has a purely imaginary root λ = iω  (ω > 0). Substituting λ = iω into (12) yieldsBy separating the real part and imaginary part, the following real part is obtained:whereLetTherefore, 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 .

If G(x) = 0 has no positive real roots, then the positive equilibrium E 22 is locally asymptotically stable for any τ > 0.

If G(x) = 0 has no positive real roots, then obviously (9) has no positive real roots. Therefore, the positive equilibrium E 22 is locally asymptotically stable for any τ > 0.

Substituting λ = iω into (22), we obtain the real part,and imaginary part,Assuming that G(x) = 0 has positive real roots, denoted by . As , we then haveLetwhere and j = 0,1, 2,….

Therefore, the characteristic equation H(λ; τ () = 0 has a pair of purely imaginary roots . For every integer j and , define λ ((τ) = α ((τ) + iω ((τ) as the root of (9) near τ (, satisfying α ((τ () = 0 and . Then the following theorem is obtained.

Theorem 5 .

If G(x) = 0 has some positive real roots, then E 22 is locally asymptotically stable for τ ∈ [0, τ (0)), where

For ,  G(x) = 0 has no positive real roots when τ ∈ [0, τ (0)), which means that all the roots of (9) have strictly negative real parts when τ ∈ [0, τ (0)). Therefore, E 22 is locally asymptotically stable for τ ∈ [0, τ (0)).

Theorem 6 .

If x is a simple root of G(x) = 0, then there is a Hopf bifurcation for the system as τ increases past τ (0).

See Appendix C 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 τ (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 .

(2) The remaining roots of (22) have strictly negative real parts.

(3) ω 0 is a simple root of G(x) = 0.

We use u = τ − τ 0 to represent a new bifurcation parameter. Let X(t) = (x − x 2, e − e 2, y − y 2, v − v 2, z − z 2, w − w 2), and X (θ) = X(t + θ), where θ ∈ [−τ, 0]. Therefore, system (1) can be written as the following functional differential equation:where

By the Riesz representation theorem [28], there exists a 6 × 6 matrix-valued function such thatwhere dη(θ, u) = F 1 δ(θ)dθ + F 2 δ(θ + τ)dθ.

For ϕ ∈ C([−τ, 0], R 6), we further defineThen system (22) can be written asFor φ ∈ C([−τ, 0], R 6), defineand an inner product of ϕ, φ where η(θ) = η(θ, 0) and ϕ ∈ C([−τ, 0], R 6). Then A(0) and A (0) are adjoint operators.

Let h(θ) 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 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, g 20,  g 11,  g 02, and g 21, using the method given in [27]. The detailed computation of g 20,  g 11,  g 02, and g 21 is presented in Appendix E.

Then the following values can be computed:The signs of u 2,  β 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 sign[(dα ((τ)/dτ)|] = sign[(dG/dx)|].

Let u 2 = −Re(c 1(0))/G′(ω 0 2). We obtain the following theorem.

Theorem 7 .

Assume the hypotheses (1), (2), and (3) at the beginning of Section 4 hold.

(1) If u 2 > 0  (u 2 < 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 → −∞).

(1) If τ 0 = τ (0), where , 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(c 1(0)).

(2) If β 2 < 0, namely, Re(c 1(0)) < 0, then there exist stable periodic solutions for τ > τ (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:

We set the initial values to x(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 E 22 = (0.4757,1.1732,0.3281,1.0411,1.7470,9.8727) and the critical time value τ (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 equilibrium E 22 bifurcates into a periodic solution at τ = 0.05 (a limit cycle).

Figure 3

The positive equilibrium E 22 bifurcates into a periodic solution at τ = 0.05.

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 E 22 is asymptotically stable and the system will converge to E 22.

Figure 4

The positive equilibrium E 22 at τ = 0.03.

Figure 5

The positive equilibrium E 22 remains stale at τ = 0.03.

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 change k 1 = 0.8391 to k 1 = 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 k 1 = 0.8391 and Figure 7 when k 1 = 0.01, we can see that, with the decrease of k 1, the stable periodic solution becomes unstable, that is, asymptotically stable.

Figure 7

The positive equilibrium E 22 remains stable at τ = 0.05 when k 1 = 0.01 (blue line) and exists at stable periodic solutions when k 1 = 0.8391 (red line).

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 change k 2 = 0.1963 to k 2 = 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 k 2, 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 equilibrium E 22 remains stable at τ = 0.05 when k 2 = 0.001 (blue line) and exists at stable periodic solutions when k 2 = 0.1963 (red line).

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 E 22.

Figure 9

The positive equilibrium E 22 remains stable at τ = 0.05 when r = 0.001 (blue line) and exists at stable periodic solutions when r = 0.6933 (red line).

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.