Mathematical Model of Hepatitis B Virus Treatment with Support of Immune System

Abstract

In the current paper, the classification of the equilibrium points of an HBV mathematical model with combined therapy is presented. The influence of right-hand side changes on solution behavior is estimated, and regulation with delays in upper- and lower-bound integral limits that presents a time period with IL-2 support therapy are researched.

1. Introduction

Hepatitis B (HBV) is one of the most serious chronic diseases of the 21st century, affecting 240 million individuals worldwide [1] and resulting in approximately 350,000 deaths per year [2]. In addition, HBV has a significant influence on human health due to the highly increased probability of cirrhosis and cancer [3]. A member of the Herpesviridae class of viruses, chronic HBV increases death cases of carcinoma by at least 50% [4] and has a high mortality rate with respect to other common known diseases. The common approach for reduction of mortality cases is the development of vaccines. Even so, the immune response caused by antibody cells rapidly decreases with time. The commonly used cure for chronic HBV cases is the fusion of interferon and analogous nucleoside methods [5]. The five analogues of interferon strategies have been approved in the USA [1]. Unfortunately, the previously described approach generally does not lead to full recovery. Alternative methods can be based on a combination of IL-2 (interleukin-2), NCT02360592 and NCT00451984. The validation of fusion between described methods reveals the great potential for decreasing the viral growth rate. The addition of IL-2 therapy was defined by a mathematical model and proposed in this article. This method was validated to be highly effective. Previous models have predicted the success of an antiviral remedy, but the patient’s immune response was not taken into account [6,7,8], thereby leading to non-realistic results in prediction [9]. Therefore, we propose a robust model that will take into account the response of the immune system of the patient and will use immune therapy as a support.

2. Introduction of the Model

In [10], a mathematical model of hepatitis C was proposed. The model was defined by the system of non-linear equations described below:

$$\left\{\begin{array}{c}{X}^{\prime }=r-dX-\beta VX\hfill \\ Y^{\prime }=\beta VX-aY-pYZ\hfill \\ V^{\prime }=kY-uV-qVW\hfill \\ W^{\prime }=-hW+gVW\hfill \\ Z^{\prime }=cYZ-bZ\hfill \end{array}\right..$$

(1)

The status which describes the relationship between virus and host at each given time can be defined by the mathematical model based on the system of five differential equations with five variables. These variables are: uninfected cells, infected cells, free virus numbers, antibody response, and cytotoxic T lymphocyte (CTL) response, denoted, respectively, by X(cells/mL), Y(cells/mL), V(IU/mL), W(IU/mL), and Z(cells/mL). X cells are produced at rate r, die at rate $dX$ and are infected by the virus at rate $\beta VX$. The infected cells, Y, grow at rate $\beta VX$, die (naturally) at rate $aY$, and are killed by the CTLs at rate $pYZ$. Free viruses V are produced by infected cells Y at rate $kY$, decay at rate $uV$ and are neutralized by antibodies at rate $qVW$. Antibodies W develop in response to free viruses at rate $gvW$ and decay at rate $hW$. The number of CTLs, Z, expands in response to infected cells at rate $cYZ$ and decays in the absence of infection at rate $bZ$. The number of CTLs, Z, expands in response to infected cells at rate $cYZ$ and decays in the absence of infection at rate $bZ$.

We introduce a cure that involves two methods: interferons and nucleosides, denoted by $0<\eta <1$ and $0<\epsilon <1$, respectively The efficiency of the currency is based on two strategies: that interferons and nucleosides can be estimated, denoted by $\eta ,\epsilon$ and $0<\epsilon ,\eta <1$, respectively.

Mathematical models play a very important role in medicine (see, for example, [11,12,13]). One of the best methods for the study of viral infections and related cures as well as predictions of patient condition after treatment and immune response can be based on mathematical models described by the system of differential equations.

In this paper, the influence of an IL-2 treatment is studied, and the distributed control function in the integral form is developed. The exponential stability results with regulation are proposed.

The following modified model is here researched [14]:

$$\left\{\begin{array}{c}{X}^{\prime }=r-dX-(1-\eta )\beta VX\hfill \\ Y^{\prime }=(1-\eta )\beta VX-aY-pYZ\hfill \\ V^{\prime }=(1-\epsilon )kY-uV-qVW\hfill \\ W^{\prime }=-hW+gVW\hfill \\ Z^{\prime }=cYZ-bZ+DG\hfill \end{array}\right.,$$

(2)

where

$$G\left(t\right)={\int }_0^t{e}^{-\alpha (t-s)}Z\left(s\right)ds.$$

(3)

The modified model (2) and (3) takes into account the response of the patient immune system and uses immune therapy as a support. In the current paper, classification of equilibrium points of the system (2) is presented, with those not relevant from the medical point of view explained. In additional, the influence of right-hand side changes on solution behavior is estimated, and immune therapy as a support in integral terms with delays in upper- and lower-bound integral limits is researched (3), because drug assimilation in different patients can have varying rates of influence, a factor not taken into account in previous research.

3. Classification of Equilibrium Points

Using Maple, we found five equilibrium points of the system (2).

$$P_1=\{X,Y,V,W,Z,G\}=\left(\frac{r}{d},0,0,0,0,0\right),$$

(4)

$$P_2=\left\{\begin{array}{c}X=\frac{\alpha u}{\beta k(\epsilon -1)(\eta -1)}\hfill \\ Y=\frac{kr(\epsilon -1)(\eta -1)\beta -d\alpha u}{\beta \alpha k(\epsilon -1)(\eta -1)}\hfill \\ V=\frac{-kr(\epsilon -1)(\eta -1)\beta +d\alpha u}{\alpha \beta u(\eta -1)}\hfill \\ W=0\hfill \\ Z=0\hfill \\ G=0\hfill \end{array}\right.,$$

$$P_3=\left\{\begin{array}{c}X=\frac{\alpha cru}{k(\eta -1)(\epsilon -1)(\alpha b-D)\beta +\alpha cdu}\hfill \\ Y=\frac{\alpha b-D}{\alpha c}\hfill \\ V=\frac{k(\epsilon -1)(-\alpha b+D)}{\alpha cu}\hfill \\ W=0\hfill \\ Z=-\frac{\left(k\right(\eta -1\left)\right(\epsilon -1\left)\right(\alpha b-cr-D)\beta +\alpha cdu)\alpha }{\left(k\right(\eta -1\left)\right(\epsilon -1\left)\right(\alpha b-D)\beta +\alpha cdu)p}\hfill \\ G=\frac{-k(\eta -1)(\epsilon -1)(\alpha b-cr-D)\beta -\alpha cdu}{\left(k\right(\eta -1\left)\right(\epsilon -1\left)\right(\alpha b-D)\beta +\alpha cdu)p}\hfill \end{array}\right.,$$

$$P_4=\left\{\begin{array}{c}X=\frac{rg}{-h(\eta -1)\beta +dg}\hfill \\ Y=-\frac{\beta hr(\eta -1)}{\alpha (-h(\eta -1)\beta +dg)}\hfill \\ V=\frac{h}{g}\hfill \\ W=\frac{(\eta -1)\left(kr\right(\epsilon -1)g+\alpha hu)\beta -\alpha dug}{(-h(\eta -1)\beta +dg)\alpha q}\hfill \\ Z=0\hfill \\ G=0\hfill \end{array}\right.,$$

$$P_5=\left\{\begin{array}{c}X=\frac{rg}{-h(\eta -1)\beta +dg}\hfill \\ Y=\frac{\alpha b-D}{\alpha c}\hfill \\ V=\frac{h}{g}\hfill \\ W=\frac{(-kb(\epsilon -1)g-chu)\alpha +gkD(\epsilon -1)}{\alpha chq}\hfill \\ Z=-\frac{\alpha \left(h\right(\eta -1\left)\right(\alpha b-cr-D)\beta -dg(\alpha b-D\left)\right)}{(\alpha b-D)\left(h\right(\eta -1)\beta -dg)p}\hfill \\ G=\frac{-h(\eta -1)(\alpha b-cr-D)\beta +dg(\alpha b-D)}{(\alpha b-D)\left(h\right(\eta -1)\beta -dg)p}\hfill \end{array}\right..$$

$P_1$ is the disease-free equilibrium, so it is the state of a healthy subject. In our system, the variables W and Z represents a patient’s immune system. Therefore, $P_2$ represents a patient with immune system failure. Even if, from a modeling perspective, we consider variables having zero values, in biology, this rarely happens. Therefore, $P_2$ is an extreme case that can approximate the state of HBV infection in immuno-compromised patients (for example, those affected by common variable immunodeficiency), where there are very low values of both cytotoxic and humoral response (see, for example, [15]).

Additionally, it was proven that the HBV infection state causes an immuno-depressive state (see, for example, [16]).

Similarly, even $P_3$ and $P_4$ can be considered as extreme cases of real clinical conditions, such as X-linked agammaglobulinemia or HIV infection (see, for example, [17,18]).

Finally, $P_5$ represents the most meaningful state to describe a chronic patient with HBV infection, because it represents the condition of a patient affected by the virus, but with a still existent and functional immune system.

In this article, the disease-free equilibrium point $P_1$ is studied since it represents a patient who is either a healthy subject or a recovered one. From this stationary state, it is possible to deduce the others, since we solve the system in the neighborhood of the equilibrium-free point.

4. Influence of Right-Hand Side Changes on Solution Behavior

Reducing the integro-differential system (2) to the system of ordinary differential equations obtains:

$$\left\{\begin{array}{c}{X}^{\prime }=r-dX-(1-\eta )\beta VX\hfill \\ Y^{\prime }=(1-\eta )\beta VX-aY-pYZ\hfill \\ V^{\prime }=(1-\epsilon )kY-uV-qVW\hfill \\ W^{\prime }=-hW+gVW\hfill \\ Z^{\prime }=cYZ-bZ+DG\hfill \\ G^{\prime }=Z-\alpha G\hfill \end{array}\right..$$

(5)

Linearizing system (5) in the disease-free equilibrium (4), a linear system is obtained:

$$\left\{\begin{array}{c}{x}_1^{\prime }=r-d{x}_1-\frac{(1-\eta )\beta r}{d}{x}_3\hfill \\ x_2^{\prime }=-a{x}_2+\frac{(1-\eta )\beta r}{d}{x}_3\hfill \\ x_3^{\prime }=(1-\epsilon )k{x}_2-u{x}_3\hfill \\ x_4^{\prime }=-h{x}_4\hfill \\ x_5^{\prime }=D{x}_6-b{x}_5\hfill \\ x_6^{\prime }=x_5-\alpha x_6\hfill \end{array}\right.,$$

(6)

where

$$x_1=X-\frac{r}{d},x_2=Y,x_3=V,x_4=W,x_5=Z,x_6=G,$$

and corresponding homogeneous systems are as follows:

$$\left\{\begin{array}{c}{x}_1^{\prime }=-d{x}_1-\frac{(1-\eta )\beta r}{d}{x}_3\hfill \\ x_2^{\prime }=-a{x}_2+\frac{(1-\eta )\beta r}{d}{x}_3\hfill \\ x_3^{\prime }=(1-\epsilon )k{x}_2-u{x}_3\hfill \\ x_4^{\prime }=-h{x}_4\hfill \\ x_5^{\prime }=D{x}_6-b{x}_5\hfill \\ x_6^{\prime }=x_5-\alpha x_6\hfill \end{array}\right..$$

(7)

Denote the coefficients matrix (7) as

$$A=\left(\begin{array}{cccccc}-d& 0& -\frac{(1-\eta )\beta r}{d}& 0& 0& 0\\ 0& -a& \frac{(1-\eta )\beta r}{d}& 0& 0& 0\\ 0& (1-\epsilon )k& -u& 0& 0& 0\\ 0& 0& 0& -h& 0& 0\\ 0& 0& 0& 0& -b& D\\ 0& 0& 0& 0& 1& -\alpha \end{array}\right) .$$

(8)

Let us denote

$$(\epsilon -1)(\eta -1)<\frac{aud}{\beta kr},$$

(9)

$$D<\alpha b.$$

(10)

Theorem 1.

([14]). If all the coefficients of system (2) are positive, and η and ε are parameters defined between 0 to 1, and if inequalities (9) and (10) are fulfilled, then system (2) is exponentially stable (see [19]) in the neighborhood of disease-free equilibrium.

Constructing a mathematical model, we neglect the influences of different factors that seem nonessential. We also cannot know exactly the values of the coefficients describing the model.

Consider the system

$$Y^{\prime }\left(t\right)-AY\left(t\right)=F\left(t\right)+\Delta F\left(t\right),$$

(11)

where the matrix A is the coefficient matrix of system (6) defined by (8), and $\Delta F\left(t\right)\in L_{\infty }^6$ describes the change of the right-hand side. In the following assertion, we estimate the difference between the solution-vector.

Theorem 2.

Under the assumption of Theorem 1, system (5) is exponentially stable in the neighborhood of disease-free equilibrium, and the following inequality is true:

$$\left|\right|Y\left(t\right)-X\left(t\right)\left|\right|\le \left|\right|C\left|\right|\left|\right|\Delta F\left|\right|,$$

where $C(t,s)$ is the Cauchy matrix of system (6) and

$$∥C∥=\underset{1\le i\le 6}{max}\left(\underset{t\ge 0}{sup}{\int }_0^t\sum _{j=0}^6\left|c_{ij}\left(t,s\right)\right|\right)ds,∥△F\left(t\right)∥=\underset{1\le i\le 6}{max}ess\underset{t\ge 0}{sup}\left|△F_i\left(t\right)\right|,$$

$$\left|\right|Y\left(t\right)-X\left(t\right)\left|\right|=ma{x}_{1\le i\le 6}su{p}_{t\ge 0}|y_i\left(t\right)-x_i\left(t\right)|.$$

Proof. 

The proof follows from the solution representation of system (6), and the estimates $\left|\right|C\left|\right|$ can be obtained from estimates of Cauchy matrix elements of system (6) obtained in [14]. □

5. Exponential Stability of the System with Delays in Upper and Lower Limits of Control Function

Consider system (2), where

$$G\left(t\right)={\int }_{{\tau }_1\left(t\right)}^{t-{\tau }_2\left(t\right)}{e}^{-\alpha (t-s)}Z\left(s\right)ds,$$

(12)

where ${\tau }_1\left(t\right)<t-{\tau }_2\left(t\right),{\tau }_1\left(t\right)>0,{\tau }_2\left(t\right)<t$.

Let us denote

$$\tilde{G}\left(t\right)={\int }_0^t{e}^{-\alpha (t-s)}Z\left(s\right)ds.$$

(13)

We can write (12) in the following form:

$$\begin{array}{c}\hfill G\left(t\right)=e^{-\alpha {\tau }_2\left(t\right)}{\int }_0^{t-{\tau }_2\left(t\right)}{e}^{-\alpha \left((t-{\tau }_2\left(t\right))-s\right)}Z\left(s\right)ds-e^{-\alpha (t-{\tau }_1\left(t\right))}{\int }_0^{{\tau }_1\left(t\right)}{e}^{-\alpha \left({\tau }_1\left(t\right)-s\right)}Z\left(s\right)ds=\\ \hfill =e^{-\alpha {\tau }_2\left(t\right)}\tilde{G}(t-{\tau }_2\left(t\right))-e^{-\alpha (t-{\tau }_1\left(t\right))}\tilde{G}\left({\tau }_1\left(t\right)\right).\end{array}$$

(14)

Reducing the integro-differential system (2) $G\left(t\right)$ defined by (12) to the system of ordinary differential equations, we obtain:

$$\left\{\begin{array}{c}{X}^{\prime }\left(t\right)=r-dX\left(t\right)-(1-\eta )\beta V\left(t\right)X\left(t\right)\hfill \\ Y^{\prime }\left(t\right)=(1-\eta )\beta V\left(t\right)X\left(t\right)-aY\left(t\right)-pY\left(t\right)Z\left(t\right)\hfill \\ V^{\prime }\left(t\right)=(1-\epsilon )kY\left(t\right)-uV\left(t\right)-qV\left(t\right)W\left(t\right)\hfill \\ W^{\prime }\left(t\right)=-hW\left(t\right)+gV\left(t\right)W\left(t\right)\hfill \\ Z^{\prime }\left(t\right)=cY\left(t\right)Z\left(t\right)-bZ\left(t\right)+D\tilde{G}\left(t\right)+e^{-\alpha {\tau }_2\left(t\right)}\tilde{G}(t-{\tau }_2\left(t\right))-e^{-\alpha (t-{\tau }_1\left(t\right))}\tilde{G}\left({\tau }_1\left(t\right)\right)-D\tilde{G}\left(t\right)\hfill \\ {\tilde{G}}^{\prime }\left(t\right)=Z\left(t\right)-\alpha \tilde{G}\left(t\right)\hfill \end{array}\right..$$

(15)

Linearizing system (15) in the neighborhood of disease-free equilibriun (4), the following homogeneous system is obtained:

$$\left\{\begin{array}{c}{x}_1^{\prime }\left(t\right)=-d{x}_1\left(t\right)-\frac{(1-\eta )\beta r}{d}{x}_3\left(t\right)\hfill \\ x_2^{\prime }\left(t\right)=-a{x}_2\left(t\right)+\frac{(1-\eta )\beta r}{d}{x}_3\left(t\right)\hfill \\ x_3^{\prime }\left(t\right)=(1-\epsilon )k{x}_2\left(t\right)-u{x}_3\left(t\right)\hfill \\ x_4^{\prime }\left(t\right)=-h{x}_4\left(t\right)\hfill \\ x_5^{\prime }\left(t\right)=D{x}_6\left(t\right)-b{x}_5\left(t\right)+e^{-\alpha {\tau }_2\left(t\right)}{x}_6(t-{\tau }_2\left(t\right))-e^{-\alpha (t-{\tau }_1\left(t\right))}{x}_6\left({\tau }_1\left(t\right)\right)-D{x}_6\left(t\right)\hfill \\ x_6^{\prime }\left(t\right)=x_5\left(t\right)-\alpha x_6\left(t\right)\hfill \end{array}\right.,$$

(16)

where

$$x_1=X-\frac{r}{d},x_2=Y,x_3=V,x_4=W,x_5=Z,x_6=\tilde{G}.$$

Let us denote

$$\begin{array}{c}{c}_{11}=\frac{d(a+u)}{2d},c_{12}=\frac{\sqrt{4k\beta rd(\epsilon -1)(\eta -1)+d^2{(a-u)}^2}}{2d},c_{13}=\frac{\alpha +b}{2},\hfill \\ c_{14}=\frac{\sqrt{{(\alpha -b)}^2+4D}}{2},c_{15}=\frac{\beta r(\eta -1)}{d(-c_{11}+c_{12}+d)},c_{16}=\frac{\beta r(\eta -1)}{d(-c_{11}+c_{12}+a)},\\ c_{17}=\frac{\beta r(\eta -1)}{d(-c_{11}-c_{12}+d)},c_{18}=\frac{\beta r(\eta -1)}{d(-c_{11}-c_{12}+a)},c_{19}=\frac{D}{-\frac{\alpha }{2}+\frac{b}{2}+c_{14}},\\ c_{21}=\frac{D}{-\frac{\alpha }{2}+\frac{b}{2}-c_{14}},c_{22}=\frac{c_{15}-c_{17}}{c_{16}-c_{18}},c_{23}=\frac{1}{c_{16}-c_{18}},c_{25}=\frac{c_{15}·c_{18}-c_{16}·c_{17}}{c_{16}-c_{18}},\\ \hfill c_{26}=\frac{c_{18}}{c_{16}-c_{18}},c_{27}=\frac{c_{16}}{c_{16}-c_{18}},c_{28}=\frac{1}{c_{19}-c_{21}},c_{31}=\frac{c_{21}}{c_{19}-c_{21}},c_{32}=\frac{c_{19}}{c_{19}-c_{21}}.\end{array}$$

It can be assumed that the denominators are not zero.

The columns of the Cauchy matrix of system (6) (see the building of the Cauchy matrix in [14]):

$$C_1(t,s)=\left(\begin{array}{c}1\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right)e^{-d(t-s)},$$

$$C_2(t,s)=\left(\begin{array}{c}{c}_{22}\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right)e^{-d(t-s)}+\left(\begin{array}{c}-c_{15}·c_{23}\\ c_{16}·c_{23}\\ -c_{23}\\ 0\\ 0\\ 0\end{array}\right)e^{(-c_{11}+c_{12})(t-s)}+\left(\begin{array}{c}{c}_{17}·c_{23}\\ -c_{18}·c_{23}\\ c_{23}\\ 0\\ 0\\ 0\end{array}\right)e^{(-c_{11}-c_{12})(t-s)},$$

$$C_3(t,s)=\left(\begin{array}{c}{c}_{25}\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right)e^{-d(t-s)}+\left(\begin{array}{c}-c_{15}·c_{26}\\ c_{16}·c_{26}\\ -c_{26}\\ 0\\ 0\\ 0\end{array}\right)e^{(-c_{11}+c_{12})(t-s)}+\left(\begin{array}{c}{c}_{17}·c_{27}\\ -c_{18}·c_{27}\\ c_{27}\\ 0\\ 0\\ 0\end{array}\right)e^{(-c_{11}-c_{12})(t-s)},$$

$$C_4(t,s)=\left(\begin{array}{c}0\\ 0\\ 0\\ 1\\ 0\\ 0\end{array}\right)e^{-h(t-s)},$$

$$C_5(t,s)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ c_{19}·c_{28}\\ c_{28}\end{array}\right)e^{(-c_{13}+c_{14})(t-s)}+\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ -c_{21}·c_{28}\\ -c_{28}\end{array}\right)e^{(-c_{13}-c_{14})(t-s)},$$

$$C_6(t,s)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ -c_{19}·c_{31}\\ -c_{31}\end{array}\right)e^{(-c_{13}+c_{14})(t-s)}+\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ c_{21}·c_{32}\\ c_{32}\end{array}\right)e^{(-c_{13}-c_{14})(t-s)}.$$

(17)

We can write system (7) in the following form:

$$X^{\prime }\left(t\right)=AX\left(t\right),$$

(18)

where

$$X\left(t\right)=col\{x_1\left(t\right),x_2\left(t\right),x_3\left(t\right),x_4\left(t\right),x_5\left(t\right),x_6\left(t\right)\}.$$

It is known that the general solution of the system

$$X^{\prime }\left(t\right)-AX\left(t\right)=F\left(t\right)$$

(19)

can be written in the following form (see [19]):

$$X\left(t\right)={\int }_0^{t}C(t,s)F\left(s\right)ds+C(t,0)X\left(0\right),$$

(20)

where $C(t,s)$ is a Cauchy matrix of system (7). We can rewrite (16) in the following form:

$$\left\{\begin{array}{c}{x}_1^{\prime }\left(t\right)+d{x}_1\left(t\right)+\frac{(1-\eta )\beta r}{d}{x}_3\left(t\right)=0\hfill \\ x_2^{\prime }\left(t\right)+a{x}_2\left(t\right)-\frac{(1-\eta )\beta r}{d}{x}_3\left(t\right)=0\hfill \\ x_3^{\prime }\left(t\right)-(1-\epsilon )k{x}_2\left(t\right)+u{x}_3\left(t\right)=0\hfill \\ x_4^{\prime }\left(t\right)+h{x}_4\left(t\right)=0\hfill \\ x_5^{\prime }\left(t\right)-D{x}_6\left(t\right)+b{x}_5\left(t\right)=D\left(e^{-\alpha {\tau }_2\left(t\right)}{x}_6(t-{\tau }_2\left(t\right))-e^{-\alpha (t-{\tau }_1\left(t\right))}{x}_6\left({\tau }_1\left(t\right)\right)-x_6\left(t\right)\right)\hfill \\ x_6^{\prime }\left(t\right)-x_5\left(t\right)+\alpha x_6\left(t\right)=0\hfill \end{array}\right..$$

(21)

Without loss of generality, we can assume that $X\left(0\right)=0$.

Substituting $X\left(t\right)={\int }_0^{t}C(t,s)W\left(s\right)ds$ into system (19), we obtain

$$x_6(t-\tau \left(t\right))=\sum _{i=1}^6{\int }_0^{t-\tau \left(t\right)}{c}_{6i}(t-\tau \left(t\right),s)w_i\left(s\right)ds,$$

and for $F\left(t\right)$,

$$F\left(t\right)=(\Omega F)\left(t\right),$$

where $F\left(t\right)=col\{f_1\left(t\right),f_2\left(t\right),f_3\left(t\right),f_4\left(t\right),f_5\left(t\right),f_6\left(t\right)\}$, and $\Omega :L_{\infty }^5\to L_{\infty }^5$ is the operator that is defined by

$$\left(\Omega F\right)\left(t\right)=col\{0,0,0,0,\omega F\left(t\right),0\},$$

where

$$\begin{array}{c}\omega F\left(t\right)=D{e}^{-\alpha {\tau }_2\left(t\right)}\sum _{i=1}^6{\int }_0^{t-{\tau }_2\left(t\right)}{c}_{6i}(t-{\tau }_2\left(t\right),s)w_i\left(s\right)ds\\ \\ -D{e}^{-\alpha (t-{\tau }_1\left(t\right))}\sum _{i=1}^6{\int }_0^{{\tau }_1\left(t\right)}{c}_{6i}({\tau }_1\left(t\right),s)w_i\left(s\right)ds\\ \\ -D\sum _{i=1}^6{\int }_0^t{c}_{6i}(t,s)w_i\left(s\right)ds.\end{array}$$

Let us denote ${\tau }_1^{*}=ess{sup}_{t\ge 0}\left|{\tau }_1\left(t\right)\right|$ and ${\tau }_2^{*}=ess{sup}_{t\ge 0}\left|{\tau }_2\left(t\right)\right|$, such that the following inequalities are fulfilled ${\tau }_1\left(t\right)<t-{\tau }_2\left(t\right),{\tau }_1\left(t\right)>0,{\tau }_2\left(t\right)<t$ and

$$R=\frac{|c_{28}|+|c_{31}|}{|c_{14}-c_{13}|}\left(e^{\alpha {\tau }_1^{*}}+e^{-{\tau }_2^{*}(c_{14}-c_{13})}+1\right)+\frac{|c_{28}|+|c_{32}|}{|c_{13}+c_{14}|}\left(e^{\alpha {\tau }_1^{*}}+e^{{\tau }_2^{*}(c_{14}+c_{13})}+1\right).$$

Theorem 3.

If assumptions of Theorem 1 are fulfilled, and $R<1$, then the disease-free equilibrium is exponentially stable.

Proof. 

Estimating the norm of operator $\Omega$, we obtain the assertion of this theorem. □

6. Simulations

Remark 1.

From [20], the following values of the parameters in system (1) were used:

$$d=0.00333,g=q=5,b=0.112,\beta =7,h=2,a=0.56,u=0.67,p=c=5.14,$$

$$k=20,r=6.17\ast {10}^{-4}.$$

Using Theorem 1, we have to choose parameters $\epsilon ,\eta ,\alpha$ and D, such that the following inequalities are fulfilled:

$$(\epsilon -1)(\eta -1)<0.014,\frac{D}{\alpha }<0.112.$$

Figure 1: initial conditions of $X\left(0\right)=0.01,Y\left(0\right)=1,V\left(0\right)=0.5,W\left(0\right)=0.1,Z\left(0\right)=0.2,$ and $U\left(0\right)=0.5$ and drug dose of $\epsilon =0.9,\eta =0.9,\alpha =10,$ and $D=1.1$.

Figure 2: initial conditions of $X\left(0\right)=0.1,Y\left(0\right)=0.5,V\left(0\right)=0.5,W\left(0\right)=0.1,Z\left(0\right)=0.1,$ and $U\left(0\right)=0$ and drug dose of $\epsilon =0.8,\eta =0.95,\alpha =20,$ and $D=2$.

Figure 3: initial conditions of $X\left(0\right)=0.5,Y\left(0\right)=0.1,V\left(0\right)=0.1,W\left(0\right)=0.1,Z\left(0\right)=0.1,$ and $U\left(0\right)=0.1$ and drug dose of $\epsilon =0.5,\eta =0.99,\alpha =50,$ and $D=5.5$.

The results of the simulations are presented in Figure 1, Figure 2 and Figure 3 accordingly.

From the figures, we see the immune system response to the virus. The immune system reacts to a large number of infected cells. This reaction occurs with a delay, and when the number of infected cells is small enough, the immune system stabilizes.

7. Conclusions

In this paper, a mathematical model for hepatitis B virus combination treatment is presented, considering a standard of care and IL-2 in the integral form.

In Section 3, the classification of the equilibrium points of a hepatitis B virus combination treatment mathematical model with combined therapy is given.

$P_1$ is the disease free equilibrium, so it is the state of a healthy subject. $P_2$ represents a patient with a failure in its immune system. $P_2$ is an extreme case that can approximate the state of HBV infection in immuno-compromised patients (for example, those affected by common variable immunodeficiency), where there are very low values of both cytotoxic and humoral response. Similarly, even $P_3$ and $P_4$ can be considered as extreme cases of real clinical conditions, such as X-linked agammaglobulinemia or HIV infection. Finally, $P_5$ represents the most meaningful state to describe a chronic patient with HBV infection. Thus, from a medical point of view, equilibrium points $P_2$, $P_3$ and $P_4$ are not relevant.

As noted, this is because the assimilation of a drug in the body of different patients can have different rates of influence—a factor overlooked in the previous works. Section 4 and Section 5 offer estimations of the influence of right-hand side changes on solution behavior. In addition, regulation with delays in upper- and lower-bound integral limits that present a time period with IL-2 support therapy is evaluated.

In Section 6, in order to validate the proposed model (5), a validation set was chosen, one defined by the initial conditions and dose of the drugs.