Mathematical modeling of hepatitis B infection with vaccination and optimal control interventions

Abstract

In this paper, we formulate a mathematical model to analyze the impact of hospitalization and vaccination on the transmission dynamics of hepatitis B infection. Initially, we evaluate the model equilibria and basic reproduction number. Further, using Lyapunov function we prove that the model has a globally asymptotically stable infection free equilibrium when the basic reproduction number is less than 1, and using geometric approach we show global asymptotic stability of endemic equilibrium, when $\mathcal{R}_{0}>$ 1. We perform sensitivity analysis of $\mathcal{R}_{0}$ to identify the dominant parameters that seriously affect the hepatitis B infection. After that, using optimal control theory and Pontryagin’s maximum principle, we develop the control problem and illustrate the necessary optimality conditions. Finally, we perform the numerical simulations in order to point out the effectiveness of control interventions.

Keywords

Hepatitis B infection model stability simulations control analysis hospitalization

1. Introduction

Hepatitis B infection caused by hepatitis B virus (HBV) is one the major and life threatening health problems in the world. This viral disease also known as blood-bone infection mainly attack on liver in human body, leading to cirrhosis and liver cancer. The infection of HBV can be either acute or chronic. In acute stage many people may not show any symptoms up to 6 months. However, few develop a rapid onset of sickness with dark urine, joint pain, fatigue, jaundice, vomit and severe pain. In the chronic stage this infection may either be asymptomatic or may be related with a chronic liver infection which may eventually develop cirrhosis over a period of several years [1]. This infection cannot spread among people by sharing eating utensils, holding hands, kissing, coughing, hugging, breastfeeding or sneezing. Most commonly, HBV infection can transmit among people due to reuse of syringes and needles, transfusion of blood and other body fluids, unprotected sex with infected partner etc. This type of HBV transmission is known as horizontal. It can also be transmitted vertically from infected mother to her child during childbirth. A mother infected with HBV infection has a twenty percent risk of shifting this infection to her offspring during delivery. According to WHO reports, HBV infection is mostly notified in Africa, Southern Europe, Asia and Latin America. Round about 240 million people have chronic HBV infection and more than 0.78 million die each year [1].

Mathematical modeling of transmittable diseases are one of the prominent tools to investigate the complex dynamics and provide a clue for possible controlling of the infection from the community. Using these models one can predict about the current and long term behavior of the concern disease which are beneficial for health authorities to take further steps for the disease eradication. A number of transmission models have been developed to explore the dynamics of HBV in different countries in recent years. Thornley et al. [3], proposed an HBV transmission model to estimate the effectiveness of health control strategies including infant vaccination in the New Zealand population. Pang et al. [4], constructed a mathematical model of HBV incorporating the vaccinated class in order to predict the impact of the vaccination intervention in the population of China. Khan et al. [5] developed an HBV model with migrant class and studied the impact of immigrants on the prevalence of HBV infection. Unprotected sexual contacts is one of the main causes of HBV transmission, therefore, Zou et al. [6] studied the impact of sexual transmission in China through a mathematical model.

Optimal control theory has been derived from the calculus of variation and is a strong mathematical tool in order to quantify the controlling strategies for a disease [7]. The dynamics of the epidemiological models are expressed by a state variable(s). The main goal of using this theory is to attain the maximum benefits with a lowest cast. A number of powerful techniques have been developed in the existing work to evaluate the optimal solution for various epidemiological models including HBV [8, 9]. Khan et al. [10], constructed an optimal control model for HBV infection and used the Pontryagin’s Principle to seek the desired optimal problem of the proposed mathematical models with the given constraints. They have considered three control measures that are isolation, vaccination and educational awareness. Recently, to carry out the impact of media coverage and other controlling parameters in the elimination of HBV in the community, a transmission model has been explored in [11].

Motivated by the existing literature we develop mathematical model for the dynamics of HBV infection with vaccinated and hospitalized classes. In early stage, vaccination is effective but in the patients who develop chronic infection need to visit hospital or health care centers to take appropriate treatment including antiviral medication such as tenofovir or interferon. Further, we will modify the proposed HBV model by including three different control functions (or variables) depending on time. The detail of remanning sections of the manuscript are: In the next section, we develop the HBV model. In Section 3, we investigate the basic mathematical analysis of the HBV model. In Section 4, we perform the threshold sensitivity analysis. The global stability results of the model equilibria are presented in Section 5. The control analysis of the system is discussed in Section 6. Numerical simulations with brief discussion of four different control strategies are presented in Section 7. Finally a brief concluding remarks are drawn in Section 8.

2. Formulation of the model

In the present section, we construct the mathematical model consisting of seven non-linear differential equations (DEs), in order to analyze the transmission dynamics of HBV infection. To develop the model, the total population $N(t)$ , is sub-divided into seven epidemiological states: susceptible $S(t)$ , exposed $E(t)$ , acute $A(t)$ , carrier $C(t)$ , hospitalized class $H(t)$ , recovered $R(t)$ and vaccinated class $V(t)$ such that

$\displaystyle N(t)=S(t)+E(t)+C(t)+A(t)+H(t)+V(t)+R(t).$

The hospitalized compartment contain those people who need proper treatment in the acute and chronic stages. The proposed model is expressed via the following system of nonlinear DEs describing the dynamics of HBV infection:

$\displaystyle\left\{\begin{array}[]{l}\frac{dS(t)}{dt}=\pi\omega(1-\xi C(t))-% \beta(A(t)+\eta C(t))S(t)-(\mu_{0}+\nu)S(t)+\varphi V(t)\\ \frac{dE(t)}{dt}=\beta(A(t)+\eta C(t))S(t)-(\delta+\mu_{0})E(t)+\pi\omega\xi C% (t),\\ \frac{dA(t)}{dt}=\delta E(t)-(h_{A}+\mu_{1}+\mu_{0}+\gamma)A(t),\\ \frac{dC(t)}{dt}=\gamma A(t)-(\mu_{2}+\mu_{0}+h_{C})C(t),\\ \frac{dH(t)}{dt}=h_{A}A(t)+h_{C}C(t)-(\mu_{0}+\kappa)H(t),\\ \frac{dR(t)}{dt}=\kappa H(t)-\mu_{0}R(t),\\ \frac{dV(t)}{dt}=\pi(1-\omega)+\nu S-(\varphi+\mu_{0})V(t).\end{array}\right.$ (1)

with the corresponding non-negative initial conditions given as

$\displaystyle S(0)=S_{0},E(0)=E_{0},A(0)=A_{0},C(0)=C_{0},H(0)=H_{0},R(0)=R_{0% },V(0)=V_{0}.$

In the model Eq. (1) we consider both vertical and horizontal transmission. The birth rate is denoted by parameter $\pi$ , the natural death rate is $\mu_{0}$ in all epidemiological stages and the disease induced rate in acute and chronic stages is $\mu_{1}$ and $\mu_{2}$ respectively. The proportion of unimmunized new born babies is denoted by $\omega$ , the parameter $\xi$ denotes parentally infected individuals rate, $\varphi$ is the immunity waning rate of vaccinated individuals and they rejoin the suspectable class. For the newborns, the factor $\pi(1-\omega)$ is assumed to be immunized successfully and enter to vaccinated class, whereas in the case of unsuccessfully immunized, they individuals $\pi\omega\xi C$ enter to exposed class and the rest $\pi\omega(1-\xi C)$ stay in the suspectable class. The effective contact rate between suspectable and HBV infected people is $\beta$ , the infectiousness of carriers corresponding to acute infection is denoted by $\eta$ and the vaccination rate of suspectable people is $\nu$ . The exposed individuals enter to acute stage at the rate $\delta$ , the acutely infected individuals further develop the infection and enter to chronic stage with the rate denoted by $\gamma$ . The population in the $A$ and $C$ compartments move to hospitalized class with the transmission rates $h_{A}$ and $h_{C}$ respectively whose further recovered after proper treatment and enter to recovered class at the rate $\kappa$ . From the system Eq. (1) it is noted that variable $R$ does not occurs in the reaming equations, so we can ignore the 6th equation of the model. For convenience in our latter calculations let us denote

$\displaystyle k_{1}=(\mu_{0}+\nu),k_{2}=(\delta+\mu_{0}),k_{3}=(h_{A}+\mu_{1}+% \mu_{0}+\gamma),$ $\displaystyle k_{4}=(\mu_{2}+\mu_{0}+h_{C}),k_{5}=(\mu_{0}+\kappa),k_{6}=(% \varphi+\mu_{0}).$

3. Basic mathematical properties

In this section we investigate some basic analysis of the proposed HBV model.

3.1 Invariant region

The HBV model given in Eq. (1) will be studied in biologically feasible region as given below:

$\displaystyle\Theta=\{(S(t),E(t),A(t),C(t),H(t),V(t))\in\mathbb{R}^{6}_{+}:N(t% )\leqslant\pi/\mu_{0}\}.$

From the addition of all equations in the above model Eq. (1) we get

$\displaystyle\frac{dN(t)}{dt}=\pi-\mu_{0}N(t)-\mu_{1}A-\mu_{2}C\leqslant\pi-% \mu_{0}N(t).$

Clearly, $\frac{dN(t)}{dt}\leqslant 0$ if $N\geqslant\pi/\mu_{0}$ . Further, as

$\displaystyle\frac{dN(t)}{dt}+\mu_{0}N(t)\leqslant\pi.$ (2)

Integrating Eq. (2) and implementing the standard comparison theorem [12] we obtained

$\displaystyle 0\leqslant N\leqslant\pi/\mu_{0}(1-e^{-\mu_{0}t})+N(0)e^{-\mu_{0% }t}.$

If $N(0)\leqslant\pi/\mu_{0}$ , then $N(t)\leqslant\pi/\mu_{0}$ . All the initial solutions in $\Theta$ remain in $\Theta$ for all positive values of $t$ . Thus, the region $\Theta$ is positively invariant and attracts all solutions in $\mathbb{R}^{6}_{+}$ .

3.2 Model equilibria

The model Eq. (1) has two biologically meaningful equilibria which are disease-free equilibrium (DFE) and endemic equilibrium (EE). The DFE $X_{0}$ is obtained by putting all derivatives in the model equation to zero and solving the resulting system without infection. Thus we get

$\displaystyle X_{0}=(S^{0},0,0,0,0,V^{0}),$

such that,

$\displaystyle S^{0}=\frac{\pi(\varphi+\mu_{0}\omega)}{\mu_{0}(\mu_{0}+\varphi+% \nu)},\quad V^{0}=\frac{\pi(\mu_{0}+\nu-\mu_{0}\omega)}{\mu_{0}(\mu_{0}+% \varphi+\nu)}.$

The EE $X_{**}=(S^{**},E^{**},A^{**},C^{**},H^{**},R^{**},V^{**})$ is derived by linearizing the model Eq. (1). Thus we have

$\displaystyle\left\{\begin{array}[]{l}S^{**}=\frac{S^{0}}{\mathcal{R}_{0}},\\ E^{**}=\frac{k_{3}}{\delta}A^{*},\\ C^{**}=\frac{\gamma}{k_{4}}A^{*},\\ H^{**}=\frac{h_{A}A^{*}+h_{C}C^{*}}{k_{5}},\\ R^{**}=\frac{\kappa}{\mu_{0}}H^{*},\\ V^{**}=\frac{\pi(1-\omega)+\nu S^{*}}{k_{6}}.\end{array}\right.$ (3)

3.3 Reproduction number

\mathcal{R}_{0}

The biologically threshold quantity called basic reproduction number $\mathcal{R}_{0}$ denotes the secondary infections averagely obtained by an infectious individual in a entire susceptible class. If $\mathcal{R}_{0}$ less than one, then the infection disappears, and if its value is greater than unity, then the infection will persist and can invade the population. To be more specific, it works as threshold dynamics for system Eq. (1), and is obtained by the next generation technique developed in [13]. The relevant matrices for the computations of $\mathcal{R}_{0}$ are given by

$\displaystyle\mathcal{F}=\begin{pmatrix}\beta(A+\eta C)S\\ 0\\ 0\\ \end{pmatrix},\quad\mathcal{V}=\begin{pmatrix}k_{2}E-\pi\omega\xi C\\ k_{3}A-\delta E\\ k_{4}C-\gamma A\\ \end{pmatrix}.$

The corresponding Jacobian matrices are give by

$\displaystyle F=\begin{pmatrix}0&\beta S^{0}&\eta\beta S^{0}\\ 0&0&0\\ 0&0&0\\ \end{pmatrix},\quad V=\begin{pmatrix}k_{2}&0&-\pi\xi\omega\\ -\delta&k_{3}&0\\ 0&-\gamma&k_{4}\\ \end{pmatrix}.$

Thus using [13], the expression for $\mathcal{R}_{0}$ is obtained as follows:

$\displaystyle\mathcal{R}_{0}=\frac{\beta\delta S^{o}(k_{4}+\gamma\eta)}{k_{2}k% _{3}k_{4}-\pi\gamma\delta\xi\omega}.$

4. Threshold’s sensitivity analysis

Sensitivity analysis of $\mathcal{R}_{0}$ for an epidemic models is carried out to point out the dominant factors that mainely affects the prevalence of the infection. In order to perform this analysis, the standard combination of Latin Hypercube Sampling(LHS) and partial rank correlation coefficient(PRCC) approaches are applied [14]. PRCC is one of the most efficient and effective approaches for measuring monotonic and nonlinear relationship between the model inputs and outcomes. This analysis provides PRCC and relevant $p$ values using which we can calculate the uncertainty level in a disease model. More generally, the parameters with high PRCC and small $p$ values are considered to be the highly dominant parameters in the model. Here, we accomplished PRCC analysis in order to calculate the parameters which are essential in contributing variability in the outcome of the basic reproduction number. The graphical PRCC results are depicted in Fig. 1 while the numerical PRCC and corresponding $p$ values of the model parameter taken in the analysis are given in Table 1. The PRCC values can indicate the correlation between model important parameters and $\mathcal{R}_{0}$ , whereas the plus or minus sign in Table 1 reflects that the influence is positive or negative, respectively. Specifically, the absolute values of PRCC reveal the correlation being important, moderate, or not significantly different from zero. From the Fig. 1 and Table 1, it is clear that the parameters $\beta$ , $\omega$ , $\eta$ , $\xi$ , $\delta$ and $\phi$ have positive influence, which means that $\mathcal{R}_{0}$ can be reduced by minimizing values of these parameters. On the other hand $h_{A}$ , $h_{C}$ , $\kappa$ and $\nu$ have negative impact thus the larger values of these parameters are beneficial to reduce $\mathcal{R}_{0}$ .

Table 1
PRCC and corresponding $p$ values of $\mathcal{R}_{0}$

Parameter	PRCC values	$p$ values
$\beta$	0.6676	0.0000
$\eta$	0.1181	0.0002
$\omega$	0.2958	0.000
$\xi$	0.0679	0.0331
$\delta$	0.3786	0.0000
$h_{A}$	$-$ 0.1592	0.0000
$h_{C}$	$-$ 0.0916	0.0040
$\kappa$	$-$ 0.0774	0.0150
$\nu$	$-$ 0.2964	0.0000
$\phi$	0.3394	0.0000

Figure 1.

Sensitivity analysis of model parameters and PRCC results for $\mathcal{R}_{0}$ .

5. Global stability analysis

.

For $\mathcal{R}_{0}<$ 1, the DFE, $X_{0}$ of the model Eq. (1), is globally asymptotically stable (GAS) on $\Theta$ and for $\mathcal{R}_{0}>$ 1 it is unstable.

Proof..

To confirm the GAS of the DFE, we consider the following Lyapunov function:

$\displaystyle W(t)=aE(t)+bA(t)+cC(t),$

where $a$ , $b$ , and $c$ are arbitrary positive constants to be chosen latter. Taking derivative of $W(t)$ we obtain

$\displaystyle\frac{dW(t)}{dt}=aE^{\prime}+bA^{\prime}+cC^{\prime}.$

Using system Eq. (1), we get

$\displaystyle\frac{dW(t)}{dt}=a\{\beta(A+\eta C)S+\pi\omega\xi C-k_{2}E\}+b\{% \delta E-k_{3}A\}+c\{\gamma A-k_{4}C\}\leqslant a\{\beta(A+\eta C)S^{0}+\pi% \omega\xi C-k_{2}E\}+b\{\delta E-k_{3}A\}+c\{\gamma A-k_{4}C\}=\{a\beta S^{0}+% \gamma c-k_{3}b\}A+\{\delta b-k_{2}a\}E+\{a\eta\beta S^{0}+\pi\omega\xi a-k_{4% }c\}C.$ (4)

Now choosing, $a=$ 1, $b=\frac{k_{2}}{\delta}$ , $c=\frac{(\eta\beta\delta S^{0}+\pi\omega\xi)}{k_{4}}$ and simplifying, we obtained

$\displaystyle\frac{dW(t)}{dt}\leqslant\frac{k_{2}k_{3}}{\delta}\{\mathcal{R}_{% 0}-1\}A.$ (5)

Clearly, if $\mathcal{R}_{0}<$ 1 then $\frac{dW(t)}{dt}$ is negative. Hence, LaSalle’s invariant principle [20] implies that $X_{0}$ is the only singleton largest compact invariant set in $\Theta$ , hence it is GAS in $\Theta$ . ∎

We apply the result proved by Li and Muldowney in [17, 18], to show the GAS of EE which states:

.

If $\mathcal{K}$ is a compact subset in the interior of $\mathcal{D}$ , and their is $\vartheta>$ 0 and the Lozinskii measure $\bar{\mu}(B)<-\vartheta$ for all $x\in\mathcal{K}$ , then each omega limit point of system Eq. (1) in the interior of $\mathcal{D}$ is an equilibrium in $\mathcal{K}$ [21].

Where the Lozinskii measure can be calculated as [22]:

$\displaystyle\bar{\mu}(B)=\inf\{c:D_{+}c\|\mathbf{z}\|,\text{ for all solution% of }\acute{\mathbf{z}}=B\mathbf{z}\}.$

$D_{+}$ denotes the right-hand derivative. In order to show the GAS, first we need to investigate the norm $\|.\|$ for which $\bar{\mu}(B)<-\vartheta$ , for all $x\in\mathcal{K}$ .

Initially, we take the below sub-model of Eq. (1):

$\displaystyle\frac{d}{dt}S(t)=\pi\omega(1-\xi C(t))-\beta(A(t)+\eta C(t))S(t)-% (\mu_{0}+\nu)S(t)+\varphi V(t),$ $\displaystyle\frac{d}{dt}E(t)=\beta(A(t)+\eta C(t))S(t)-(\delta+\mu_{0})E(t)+% \pi\omega\xi C(t),$ (6) $\displaystyle\frac{d}{dt}A(t)=\delta E(t)-(h_{A}+\mu_{1}+\mu_{0}+\gamma)A(t),$ $\displaystyle\frac{d}{dt}C(t)=\gamma A(t)-(\mu_{2}+\mu_{0}+h_{C})C(t).$

The Jacobian of Eq. (5) is

$\displaystyle J=\begin{pmatrix}-(\beta(A+\eta C)+k_{1})&0&-\beta S&-(\pi\omega% \xi+\beta\eta S)\\ \beta(A+\eta C)&-k_{2}&\beta S&\pi\omega\xi+\beta\eta S\\ 0&\delta&-k_{3}&0\\ 0&0&\gamma&-k_{4}\\ \end{pmatrix}.$ (7)

The corresponding second additive compound matrix $M$ of $J$ is calculated as:

$\displaystyle M=\begin{pmatrix}M_{11}&\beta S&\pi\omega\xi+\beta\eta S&\beta S% &\pi\omega\xi+\beta\eta S&0\\ \delta&M_{22}&0&0&0&\pi\omega\xi+\beta\eta S\\ 0&\gamma&M_{33}&0&0&-\beta S\\ 0&\beta(A+\eta C)&0&M_{44}&0&-(\pi\omega\xi+\beta\eta S)\\ 0&0&\beta(A+\kappa C)&\gamma&M_{55}&\beta S\\ 0&0&0&0&\delta&M_{66}\\ \end{pmatrix}$ (8) $\displaystyle M_{11}=-\beta(A+\eta C)-2\mu_{0}-\delta-\nu,\quad M_{22}=-\beta(% A+\eta C)-2\mu_{0}-\nu-h_{A}-\mu_{1}-\gamma,$ $\displaystyle M_{33}=-\beta(A+\eta C)-2\mu_{0}-\nu-\mu_{2}-h_{C},\quad M_{44}=% -\delta-2\mu_{0}-h_{A}-\mu_{1}-\gamma,$ $\displaystyle M_{55}=-\delta-2\mu_{0}-h_{C}-\mu_{2},\quad M_{66}=-2\mu_{0}-h_{% A}-h_{C}-\gamma-\mu_{1}-\mu_{2}.$

Consider $K$ as follows:

$\displaystyle K=\begin{pmatrix}\frac{1}{A}&0&0&0&0&0\\ 0&\frac{1}{A}&0&0&0&0\\ 0&0&0&\frac{1}{A}&0&0\\ 0&0&\frac{1}{C}&0&0&0\\ 0&0&0&0&\frac{1}{C}&0\\ 0&0&0&0&0&\frac{1}{C}\\ \end{pmatrix},$

then,

$\displaystyle K^{-1}=\begin{pmatrix}A&0&0&0&0&0\\ 0&A&0&0&0&0\\ 0&0&0&C&0&0\\ 0&0&A&0&0&0\\ 0&0&0&0&C&0\\ 0&0&0&0&0&C\\ \end{pmatrix}.$

Calculating $K_{h}$ , which represents the derivative of $K$ in the direction of the vector field $h$

$\displaystyle K_{h}=\begin{pmatrix}-\frac{\dot{A}}{A^{2}}&0&0&0&0&0\\ 0&-\frac{\dot{A}}{A^{2}}&0&0&0&0\\ 0&0&0&-\frac{\dot{A}}{A^{2}}&0&0\\ 0&0&-\frac{\dot{C}}{C^{2}}&0&0&0\\ 0&0&0&0&-\frac{\dot{C}}{C^{2}}&0\\ 0&0&0&0&0&-\frac{\dot{C}}{C^{2}}\\ \end{pmatrix}.$

Further,

$\displaystyle K_{h}K^{-1}=\begin{pmatrix}-\frac{\dot{A}}{A}&0&0&0&0&0\\ 0&-\frac{\dot{A}}{A}&0&0&0&0\\ 0&0&-\frac{\dot{A}}{A}&0&0&0\\ 0&0&0&-\frac{\dot{C}}{C}&0&0\\ 0&0&0&0&-\frac{\dot{C}}{C}&0\\ 0&0&0&0&0&-\frac{\dot{C}}{C}\\ \end{pmatrix},$

and,

$\displaystyle KMK^{-1}=\begin{pmatrix}M_{11}&\beta S&\beta S&\frac{C(\pi\omega% \xi+\beta\eta S)}{A}&\frac{C(\pi\omega\xi+\beta\eta S)}{A}&0\\ \delta&M_{22}&0&0&0&\frac{C(\pi\omega\xi+\beta\eta S)}{A}\\ 0&\beta(A+\eta C)&M_{44}&0&0&-\frac{C(\pi\omega\xi+\beta\eta S)}{A}\\ 0&\frac{A\gamma}{C}&0&M_{33}&0&-\beta S\\ 0&0&\frac{A\gamma}{C}&\beta(A+\eta C)&M_{55}&\beta S\\ 0&0&0&0&\delta&M_{66}\\ \end{pmatrix}.$

Hence using

$\displaystyle\frac{\dot{A}}{A}=\frac{\delta E}{A}-(h_{A}+\mu_{1}+\mu_{0}+% \gamma),\quad\frac{\dot{C}}{C}=\frac{\gamma A}{C}-(\mu_{2}+\mu_{0}+h_{C}),$

We obtain the matrix $B$ as follows:

$\displaystyle B=\begin{pmatrix}B_{11}&\beta S&\beta S&\frac{C(\pi\omega\xi+% \beta\eta S)}{A}&\frac{C(\pi\omega\xi+\beta\eta S)}{A}&0\\ \delta&B_{22}&0&0&0&\frac{C(\pi\omega\xi+\beta\eta S)}{A}\\ 0&\beta(A+\eta C)&B_{44}&0&0&-\frac{C(\pi\omega\xi+\beta\eta S)}{A}\\ 0&\frac{A\gamma}{C}&0&B_{33}&0&-\beta S\\ 0&0&\frac{A\gamma}{C}&\beta(A+\eta C)&B_{55}&\beta S\\ 0&0&0&0&\delta&B_{66}\\ \end{pmatrix},$ (9)

where,

$\displaystyle B_{11}=h_{A}+\mu_{1}+\gamma-\mu_{0}-\delta-\nu-\beta(A+\eta C)-% \frac{\delta E}{A},$ $\displaystyle B_{22}=-\mu_{0}-\nu-\beta(A+\eta C)-\frac{\delta E}{A},$ $\displaystyle B_{33}=-\mu_{0}-\delta-\frac{\delta E}{A},$ $\displaystyle B_{44}=-\mu_{0}-\nu-\beta(A+\eta C)-\frac{\gamma A}{C},$ $\displaystyle B_{55}=-\mu_{0}-\delta-\frac{\gamma A}{C},$ $\displaystyle B_{66}=-\mu_{0}-h_{A}-\gamma-\mu_{1}-\frac{\phi A}{C}.$

Following [19, 21], we state the following norm on $\mathbb{R}^{6}$ :

$\displaystyle\|\mathbf{z}\|=\max\{U_{1},U_{2}\},$ (10)

where $\mathbf{z}\in\mathbb{R}^{6}$ , with $w_{i}$ , $i=1,\ldots,6$ , and $U_{1}(z_{1},z_{2},z_{2})$ are defined in similar way as in [19, 21]. Moreover, we will make use of the below inequalities:

$\displaystyle|z_{1}|,|z_{2}|,|z_{3}|,|z_{2}+z_{3}|\leqslant U_{1},$

and

$\displaystyle|z_{i}+z_{j}|,|z_{i}|,|z_{4}+z_{5}+z_{6}|\leqslant U_{2}\quad i=4% ,5,6:i\neq j.$

.

For $\mathcal{R}_{0}>$ 1, model Eq. (5) has a unique EE and is GAS provided that

$\displaystyle h_{A}+\mu_{1}+\gamma+\chi_{1}<\delta<\mu_{0}.$ (11)

Proof..

To prove the required result we will use the right derivative $D_{+}\|\mathbf{z}\|$ of the norm Eq. (10). Corresponding to different orthants and the definition of the norm Eq. (10) regarding each orthant, there are sixteen different cases.

Case 1: Let $U_{1}$ is greater than $U_{2}$ , $z_{1},z_{2},z_{3}$ are positive and $|z_{1}|>|z_{2}|+|z_{3}|$ . Then we have:

$\displaystyle\|\mathbf{z}\|=|z_{1}|=z_{1},$ (12)

therefore,

$\displaystyle D_{+}\|\mathbf{z}\|=\acute{z}_{1}=B_{11}z_{1}+B_{12}z_{2}+B_{13}% z_{3}+B_{14}z_{4}+B_{15}z_{5}\leqslant\left(h_{A}+\mu_{1}+\gamma-\mu_{0}-% \delta-\nu-\beta(A+\eta C)-\frac{\delta E}{A}\right)|z_{1}|+\beta S(|z_{2}|+|z% _{3}|)+\frac{C(\pi\omega\xi+\beta\eta S)}{A}(|z_{4}|+|z_{5}|).$

After using of $|z_{2}|+|z_{3}|<|z_{1}|$ , $|z_{4}+z_{5}|<U_{2}<|z_{1}|$ and Eq. (12) we get

$\displaystyle D_{+}\|\mathbf{z}\|\leqslant\left(h_{A}+\mu_{1}+\gamma-\mu_{0}-% \nu-\delta-\beta(A+\eta C)-\frac{\delta E}{A}+\beta S+\frac{C(\pi\omega\xi+% \beta\eta S)}{A}\right)\|\mathbf{z}\|\leqslant((h_{A}+\mu_{1}+\gamma+\chi_{1})% -(\mu_{0}+\nu+\delta)\|\mathbf{z}\|,$

where $\chi_{1}=\max\left(-\beta(A+\eta C)-\frac{\delta E}{A}+\beta S+\frac{C(\pi% \omega\xi+\beta\eta S)}{A}\right)$ . Using condition Eq. (11), we get

$\displaystyle D_{+}\|\mathbf{z}\|\leqslant-(\mu_{0}+\nu)\|\mathbf{z}\|,$

Case 2: Consider $U_{1}>U_{2}$ , $z_{1},z_{2},z_{3}>$ 0 and $|z_{1}|<|z_{2}|+|z_{3}|$ , then

$\displaystyle\|\mathbf{z}\|=|z_{2}|+|z_{3}|=z_{2}+z_{3},$

therefore,

$\displaystyle D_{+}\|\mathbf{z}\|=\acute{z}_{2}+\acute{z}_{3}=B_{21}|z_{1}|+B_% {22}|z_{2}|+B_{26}|z_{6}|+B_{32}|z_{2}|+B_{33}|z_{3}|+B_{36}|z_{6}|=\delta|z_{% 1}|-\left(\mu_{0}+\nu+\frac{\delta E}{A}\right)|z_{2}|-\left(\mu_{0}+\delta+% \frac{\delta E}{A}\right)|z_{3}|.$

Using $|z_{1}|<|z_{2}|+|z_{3}|$ , $-\delta|z_{3}|\leqslant$ 0, $-\nu|z_{2}|\leqslant$ 0, $-\frac{\delta E}{A}(|z_{2}|+|z_{3}|)\leqslant$ 0 and simplifying we get

$\displaystyle D_{+}\|\mathbf{z}\|\leqslant\delta(|z_{2}|+|z_{3}|)-\mu_{0}(|z_{% 2}|+|z_{3}|)\leqslant-(\mu_{0}-\delta)\|\mathbf{z}\|.$

Applying Eq. (11), we have $\zeta=(\mu_{0}-\delta)>$ 0, so we obtained

$\displaystyle D_{+}\|\mathbf{z}\|\leqslant-\zeta\|\mathbf{z}\|.$

The rest of fourteen cases can be proved easily. Hence, according to the result (5.2) the system Eq. (5) has a unique EE which is GAS. Now we take the last two equations of model Eq. (1) that is

$\displaystyle\frac{dH(t)}{dt}=h_{A}A(t)+h_{C}C(t)-(\mu_{0}+\kappa)H(t),\quad% \frac{dV(t)}{dt}=\pi(1-\omega)+\nu S-(\varphi+\mu_{0})V(t),$ (13)

The corresponding limit system is

$\displaystyle\frac{dH(t)}{dt}=h_{A}A^{*}+h_{C}C^{*}-(\mu_{0}+\kappa)H,\quad% \frac{dV(t)}{dt}=\pi(1-\omega)+\nu S^{*}-(\varphi+\mu_{0})V.$ (14)

Solving Eq. (14) we get

$\displaystyle H(t)=e^{-(\mu_{0}+\kappa)t}\left\{H(0)+(h_{A}A^{*}+h_{C}C^{*})% \int_{0}^{t}e^{(\mu_{0}+\kappa)s}ds\right\},$ $\displaystyle V(t)=e^{-(\varphi+\mu_{0})t}\left\{V(0)+(\pi(1-\omega)+\nu S^{*}% )\int_{0}^{t}e^{(\varphi+\mu_{0})}ds\right\}.$

Therefore, we deduce that

$\displaystyle H(t)\rightarrow H^{*},\quad V(t)\rightarrow V^{*},\text{ as}\ t% \rightarrow\infty.$

Hence, the system Eq. (1) has a GAS at EE $X_{**}$ . ∎

6. Formulation of optimal control problem

The present section aims to construct an optimal control problem for the HBV model Eq. (1). The Pontryagin’s Maximum Principle is used in order to obtain the corresponding necessary and sufficient conditions for the HBV control model. Our main objective is to provide a best controlling mechanism to eradicate the HBV from the community with a lowest cost of implementation. In order to do this, we use three control variables in the HBV model Eq. (1). The suggested controls are: $u_{1}(t)$ denotes the isolation of HBV infected individuals, $u_{2}(t)$ represents an effective public interventions (including awareness, public education etc) in order to aware the people suffering from this infection to visit hospital and the third variable $u_{3}(t)$ is used for the possible treatment of infected individuals including vaccination of suspectable people and some other medicine of hospitalized individuals. In the result of these controls, the objective functional taken in the study is given as follows:

$\displaystyle J(u_{1},u_{2},u_{3})=\min\int_{0}^{t}\left(\mathcal{P}_{1}E+% \mathcal{P}_{2}A+\mathcal{P}_{3}C+\mathcal{P}_{4}H+\frac{\mathcal{M}_{1}}{2}u_% {1}^{2}+\frac{\mathcal{M}_{2}}{2}u_{2}^{2}+\frac{\mathcal{M}_{3}}{2}u_{3}^{2}% \right)dt,$ (15)

subject to the optimal control problem for the HBV transmission Eq. (1).

$\displaystyle\frac{dS(t)}{dt}=\pi\omega(1-\xi C(t))-\beta(A(t)+\eta C(t))(1-u_% {1}(t))S(t)-\mu_{0}S(t)-\nu u_{3}S(t)+\varphi V(t),$ $\displaystyle\frac{dE(t)}{dt}=\beta(A(t)+\eta C(t))(1-u_{1}(t))S(t)-(\delta+% \mu_{0})E(t)+\pi\omega\xi C(t),$ $\displaystyle\frac{dA(t)}{dt}=\delta E(t)-h_{A}u_{2}A(t)-(\mu_{1}+\mu_{0}+% \gamma)A(t),$ $\displaystyle\frac{dC(t)}{dt}=\gamma A(t)-h_{C}u_{2}C(t)-(\mu_{2}+\mu_{0})C(t),$ (16) $\displaystyle\frac{dH(t)}{dt}=h_{A}u_{2}A(t)+h_{C}u_{2}C(t)-(\mu_{0}+\kappa(1+% u_{3}(t)))H(t),$ $\displaystyle\frac{dR(t)}{dt}=\kappa(1+u_{3}(t))H(t)-dR(t),$ $\displaystyle\frac{dV(t)}{dt}=\pi(1-\omega)+\nu u_{3}S-(\varphi+\mu_{0})V(t),$

subject to some non negetaive initial conditions. In functional Eq. (15), $\mathcal{P}_{1},\mathcal{P}_{2},\mathcal{P}_{3}$ and $\mathcal{P}_{4}$ denote respectively the weighting constants for the indivisuals $E(t)$ , $A(t)$ , $C(t)$ and $H(t)$ classes, while $\mathcal{M}_{1}$ , $\mathcal{M}_{2}$ and $\mathcal{M}_{3}$ are the weight constants used for the isolation, human intervention and treatment control respetively.

Hence, we evaluate an optimal controls $u^{*}_{1},u^{*}_{2}$ and $u^{*}_{3}$ such that,

$\displaystyle J(u^{*}_{1},u^{*}_{2},u^{*}_{3})=\min J(u_{1},u_{2},u_{3}),% \Upsilon=\{(u_{1},u_{2},u_{3})|0\leqslant u_{i}\leqslant 1,i=1,2,3,t\in[0,T_{f% }]\}.$

Further, we need to find the Lagrangian $L$ and Hamiltonian $\mathcal{H}$ for the control model, which are necessary to find the optimal solution. We proceed as follows

$\displaystyle L(E,A,C,H,u_{1},u_{2},u_{3})=\mathcal{P}_{1}E+\mathcal{P}_{2}A+% \mathcal{P}_{3}C+\mathcal{P}_{4}H+\frac{\mathcal{M}_{1}}{2}u_{1}^{2}+\frac{% \mathcal{M}_{2}}{2}u_{2}^{2}+\frac{\mathcal{M}_{3}}{2}u_{3}^{2},$ (17)

and,

$\displaystyle\!\!\!\!\!\!\mathcal{H}=\mathcal{P}_{1}E+\mathcal{P}_{2}A+% \mathcal{P}_{3}C+\mathcal{P}_{4}H+\frac{1}{2}\left(\mathcal{M}_{1}u_{1}^{2}+% \mathcal{M}_{2}u_{2}^{2}+\mathcal{M}_{3}u_{3}^{2}\right)+\lambda_{1}\{\pi% \omega(1-\xi C(t))-\beta(A(t)+\eta C(t))(1-u_{1}(t))S(t)-\mu_{0}S(t)-\nu u_{3}% S(t)+\varphi V(t)\}+\lambda_{2}\{\beta(A(t)+\eta C(t))(1-u_{1}(t))S(t)-(\delta% +\mu_{0})E(t)+\pi\omega\xi C(t)\}+\lambda_{3}\{\delta E(t)-h_{A}u_{2}A(t)-(\mu% _{1}+\mu_{0}+\gamma)A(t)\}+\lambda_{4}\{\gamma A(t)-h_{C}u_{2}C(t)-(\mu_{2}+% \mu_{0})C(t)\}+\lambda_{5}\{h_{A}u_{2}A(t)+h_{C}u_{2}C(t)-(\mu_{0}+\kappa(1+u_% {3}(t)))H(t)\}+\lambda_{6}\{\kappa(1+u_{3}(t))H(t)-dR(t)\}+\lambda_{7}\{\pi(1-% \omega)+\nu u_{3}S-(\varphi+\mu_{0})V(t)\},$ (18)

where $\lambda_{i}$ , $i=$ 1, 2, …, 7, are the adjoint variables. To obtain the optimal solution of the control problem we follow the procedure presented in [25]. The following result is established for this purpose.

.

For the optimal controls $u_{i}^{\ast}$ , $i=$ 1, 2, 3, be the solution of the above control system Eqs (15) and (6). Then there exists adjoint variables $\lambda_{i}$ , satisfying

$\displaystyle\frac{-d\lambda_{i}}{dt}=\frac{\partial\mathcal{H}}{\partial i},$ (19)

where $i=$ 1(1)7, and with transversality conditions

$\displaystyle\lambda_{i}(t_{f})=0,$ (20)

and,

$\displaystyle u_{1}^{*}=\min\left\{1,\max\left(0,\frac{\beta(A+\eta C)(\lambda% _{2}-\lambda_{1})}{\mathcal{M}_{1}}\right)\right\},$ $\displaystyle u_{2}^{*}=\min\left\{1,\max\left(0,\frac{{h_{A}A(\lambda_{3}-% \lambda_{5})}+h_{C}C(\lambda_{4}-\lambda_{5})}{{\mathcal{M}_{2}}}\right)\right\},$ (21) $\displaystyle u_{3}^{*}=\min\left\{1,\max\left(0,\frac{\nu S(\lambda_{1}-% \lambda_{7})+\kappa H(\lambda_{5}-\lambda_{6})}{\mathcal{M}_{3}}\right)\right\}.$

Proof..

The corresponding adjoint system is obtained by using condition Eq. (19). Hence, we derived the following adjoint system.

$\displaystyle\frac{d\lambda_{1}}{dt}=\beta(A+\eta C)(1-u_{1})(\lambda_{1}-% \lambda_{2})+(\mu_{0}+\nu u_{3})\lambda_{1}-\nu u_{3}\lambda_{7},$ $\displaystyle\frac{d\lambda_{2}}{dt}=-\mathcal{P}_{1}+(\mu_{0}+\delta)\lambda_% {2}-\delta\lambda_{3},$ $\displaystyle\frac{d\lambda_{3}}{dt}=-\mathcal{P}_{2}+\beta(1-u_{1})S(\lambda_% {1}-\lambda_{2})+(\mu_{1}+\mu_{0}+\gamma)\lambda_{3}-\gamma\lambda_{4}+h_{A}u_% {2}(\lambda_{3}-\lambda_{5}),$ $\displaystyle\frac{d\lambda_{4}}{dt}=-\mathcal{P}_{3}+\pi\omega\xi(\lambda_{1}% -\lambda_{2})+\beta\eta(1-u_{1})S(\lambda_{1}-\lambda_{2})+\lambda_{4}(\mu_{2}% +\mu_{0})+h_{C}u_{2}(\lambda_{4}-\lambda_{5}),$ (22) $\displaystyle\frac{d\lambda_{5}}{dt}=-\mathcal{P}_{4}+\kappa(1+u_{3})(\lambda_% {5}-\lambda_{6})+\mu_{0}\lambda_{5},$ $\displaystyle\frac{d\lambda_{6}}{dt}=\mu_{0}\lambda_{6},$ $\displaystyle\frac{d\lambda_{7}}{dt}=\varphi(\lambda_{7}-\lambda_{1})+\mu_{0}% \lambda_{7}.$

Furthermore, for the optimal control functions Eq. (6), we make use of the condition $\frac{\partial H}{\partial u_{j}}=$ 0, where $j=$ 1, 2, 3. ∎

Table 2

Parameters values with description used in numerical simulations

Notation	Description	Values	Source
$\pi$	Recruitment rate	0.0143	[3]
$\mu_{o}$	Natural mortality rate	0.0143	[3]
$\omega$	Failure of immunity	0–1	[15]
$\xi$	Un-immunized children rate	0.7	[16]
$h_{A}$	Hospitalization rate of $A(t)$	0.005	Assumed
$h_{C}$	Hospitalization rate of $C(t)$	0.005	Assumed
$\beta$	Disease transmission rate	0.8	[23]
$\delta$	Flow rate from $E(t)$ to $A(t)$	6	[2]
$\eta$	Carriers infectiousness rate due to acute infection	0–1	[24]
$\varphi$	Loss of immunity	0.06–0.03	[16]
$\nu$	Vaccinated rate of $S(t)$	0.3	[10]
$\mu_{1}$	Mortality rate due to acute infection	0.2	Assumed
$\mu_{2}$	Carrier individuals death rate	0.2	Assumed
$\gamma$	Flow rate from $A(t)$ to $C(t)$	4	[2]
$\kappa$	Recovery rate of hospitalized class	0.1	Assumed

7. Simulations

In the present part of the paper, we perform and illustrate the simulations of both without control model Eq. (1) and with control model Eq. (6) together with the control measures. Both models are solved numerically using Runge-Kutta order four scheme. The model parameters used in simulation results are presented in Table 2 most of which are taken from literature. We have taken the time level upto 50 units (years). The weight constants used in our graphical results are chosen as $\mathcal{P}_{1}=$ 0.600, $\mathcal{P}_{2}=$ 0.600, $\mathcal{P}_{3}=$ 0.200, $\mathcal{P}_{4}=$ 0.10, $\mathcal{M}_{1}=$ 0.004, $\mathcal{M}_{2}=$ 0.003 and $\mathcal{M}_{3}=$ 0.001. In the comparison Figs 2–5, the control model behavior is plotted by blue dashed line and without control model behavior is shown by a red bold line. The following strategies case-1 to 4 are adopted and their numerical results are discussed.

Strategy 1: In the present strategy, we ignore the isolation control ( $u_{1}=$ 0) and activate the rest of two controls i.e., $u_{2}=u_{3}\neq$ 0. The effectiveness of this strategy in the form of graphical results are depicted in sub-plots (a–h) of Fig. 2. We observed that implementing this strategy, the population in the suspectable, acute and chronic classes decreased significantly in contrast of hospitalized population. Whereas the individuals in the recovered and vaccinated classes are increased. The corresponding control variable profile is shown in Fig. 2(h).

Strategy 2: In this strategy, the control $u_{2}$ (awareness and public education) is non-active whereas the rest of two control variables (i.e., $u_{1}=u_{3}\neq$ 0) are active. The simulations of this case are shown in Fig. 3 with sub-plots (a–h). From the graphical behavior of this strategy we can see that the population in the exposed, acute, chronic and hospitalized compartments are decreased more rapidly in compression to case 1 while the recovered individuals increased little. The susceptible and vaccinated classes show same behavior as in the previous strategy.

Strategy 3: In the third strategy, we ignore the treatment (or vaccination) control $u_{3}=0$ and consider the remaining two control variables i.e., $u_{1}=u_{2}\neq$ 0. The graphical impact of this case is shown in Fig. 4. In contrast of the previous cases, it is observed that without treatment control the individuals in suspectable and hospitalized compartments increased but there is no effect on the population of vaccinated individuals. Also the population of exposed, acute and chronic classes are decreased to zero.

Figure 2.

Simulations of the control model Eq. (6) with $u_{2}$ and $u_{3}$ only.

Figure 3.

Simulations of the control model Eq. (6) with $u_{1}$ and $u_{3}$ only.

Figure 4.

Simulations of the control model Eq. (6) with $u_{1}$ and $u_{2}$ only.

Figure 5.

Simulations of the control model Eq. (6) with $u_{2}$ , $u_{2}$ and $u_{3}$ .

Strategy 4: We conclude from the interpretation of previous strategies that each strategy is beneficial for some classes but not collectively well for all compartment. Therefore, in the last case, we activate all the proposed control measures at a time i.e., $u_{1}\neq$ 0, $u_{2}\neq$ 0, $u_{3}\neq$ 0 and perform the numerical simulations. The graphical interpretation of this strategy is shown in Fig. 5. From Fig. 5, we can see this strategy is more prominent and necessary in order to the control of HBV infection from community.

8. Conclusion

In this paper, we have successfully formulated a mathematical model to study the dynamics of HBV infection. We incorporated both the hospitalized and vaccinated classes in the existing HBV models. Initially, the proposed model is developed without control variables. The basic properties of the model such as feasible region, basic reproduction number and equilibria are investigated. The Lypanov stability technique and geometric approach are used to prove the global stability of disease free and endemic equilibria respectively. Further, using optimal control analysis we modified the proposed model by incorporating three time dependent control variables which are: $u_{1}(t)$ (stands for isolation), $u_{2}(t)$ (used for public interventions including awareness of the people via education and media etc.,) and $u_{3}(t)$ (used for treatment along with vaccination at early stage). The Pontryagin’s maximum principle is applied in order to calculate the necessary conditions of the adjoint system. To solve both control and without control model numerically, we used RK4 technique. Furthermore, we have analyzed four different strategies for HBV eradication. The graphical results of each strategy is depicted and discussed in detail. We observed that first three strategies are beneficial for some individuals and not collectively for all classes. Finally we conclude that the fourth strategy i.e., activating all control function at the same time are more prominent and provide better results for the possible minimization of HBV infection.

References

World Health Organization Media Centre. Available: http://www.who.int/mediacentre/factsheets/fs204/en/. Accessed 2016 July 1st, (2014).

Medley

G.F.

Lindop

N.A.

Edmunds

W.J.

and Nokes

D.J.

, Hepatitis-b virus endemicity: Heterogeneity, catastrophic dynamics and control, Nat Med 7 (2001), 619–624.

Thornley

Bullen

and Roberts

, Hepatitis B in a high prevalence New Zealand population: A mathematical model applied to infection control policy, J Theo Biol 254 (2008), 599–603.

Pang

Cui

and Zhou

, Dynamical behavior of a hepatitis B virus transmission model with vaccination, J Theo Bio 265 (2010), 572–578.

Khan

M.A.

Islam

Arif

and Haq

, Transmission model of hepatitis B virus with the migration effect, Bio Res Int 2013 (2013), 1–10.

Zou

Ruan

and Zhang

, On the sexual transmission dynamics of hepatitis B virus in China, J Theo Bio 369 (2015), 1–12.

Pontryagin

L.S.

Boltyanskii

V.G.

Gamkrelidze

R.V.

Mishchenko

E.F.

, The Mathematical Theory of Optimal Processes, Wiley, New York, 1962.

Okosun

K.O.

Khan

M.A.

Bonyah

and Ogunlade

S.T.

, On the dynamics of HIV-AIDS and cryptosporidiosis, Eur Phys J Plus 132 (2017), 1–5.

Bonyah

Khan

M.A.

Okosun

K.O.

and Islam

, A theoretical model for Zika virus transmission, PLoS ONE 12 (2017), 1–26.

10.

Khan

M.A.

Islam

Valverde

J.C.

and Khan

S.A.

, Control strategies of hepatitis B with three control variables, J Bio Sys 26 (2018), 1–21.

11.

Khan

M.A.

Islam

and Zaman

, Media coverage campaign in Hepatitis B transmission model, Applied Mathematics and Computation 331 (2018), 378–393.

12.

Birkhoff

and Rota

G.C.

, Ordinary Differential Equations, Ginn, Boston, 1982.

13.

Driessche

P.V.D.

and Watmough

, Reproduction number and sub-threshold endemic equilibria for compartmental models of disease transmission, Math Bios 180 (2002), 29–48.

14.

Marino

Hogue

I.B.

Ray

C.J.

and Kirschner

, A methodology for performing global uncertainty and sensitivity analysis in systems biology, Journal of Theoretical Biology 254 (2008), 178–196. doi: 10.1016/j.jtbi.2008.04.011.

15.

World health organization, WHO/CDS/CSR/LYO/2002.2: Hepatitis B, http://apps.who.int/iris/bitstream/10665/67746/1/WHO_CDS_CSR_LYO_2002.2_HEPATITIS_B.pdf.

16.

Shepardand

C.W.

and Simard

E.P.

, Hepatitis B virus infection: Epidemiology and vaccination, Epid Rev 28 (2006), 112–125.

17.

M.Y.

, and Muldowney

J.S.

, On RA Smiths autonomous convergence theorem, R M J Math 25 (1995), 365–379.

18.

M.Y.

and Muldowney

J.S.

, A geometric approach to global-stability problems, SIAM J Math Ana 27 (1996), 1070–1083.

19.

Jana

Haldar

Nandi

S.K.

and Kar

T.K.

20.

Lasalle

J.P.

, Stability theroy for difference equations, in: Studies in Ordinary Differential Equations Hale

J.K.

, Ed. Math Assoc of America, Washington DC, 1977.

21.

Buonomo

and Vargas-De-Leon

, Global stability for an HIV-1 infection model including an eclipse stage of infected cells, J Math Anal Appl 385 (2012), 709–720.

22.

Martin

R.H.

, Logarithmic norms and projections applied to linear differential systems, J Math Anal Appl 45 (1974), 432–454.

23.

Edmunds

W.J.

Medley

G.F.

and Nokes

D.J.

, The transmission dynamics and control of hepatitis B virus in the Gambia, Stat Med 15 (2002), 2215–2233.

24.

Pang

Cui

and Zhou

, Dynamical behavior of a hepatitis B virus transmission model with vaccination, J Theo Bio 265 (2010), 572–578.

25.

Lenhart

Workman

J.T.

, Optimal Control Applied to Biological Models, CRC Press, 2007.

Mathematical modeling of hepatitis B infection with vaccination and optimal control interventions

Abstract

Keywords

1. Introduction

2. Formulation of the model

3.1 Invariant region

4. Threshold’s sensitivity analysis

Table 1 PRCC and corresponding p values of ℛ 0

.

Proof..

.

.

Proof..

.

Proof..

References

Table 1
PRCC and corresponding $p$ values of $\mathcal{R}_{0}$