Dynamic modeling and analysis of Hepatitis B epidemic with general incidence

1.   Introduction

Hepatitis B is a global health problem with a high incidence rate in developing countries. According to statistics, asymptomatic HBV carriers in the world exceed 280 million, and China accounts for about 130 million. HBV infects liver cells when it enters the body. Most Hepatitis infections are caused by viruses, infections, germs or addiction to ethyl alcohol and medicines. The transmission of HBV can take place in a variety of ways, for instance, transmission of blood, bodily fluid transmission and mother-to-child vertical transmission, and so on. Vaccination against Hepatitis B is the most basic measure for preventing and controlling the disease. Many mathematicians and biologists have studied the Hepatitis B epidemic. They built different mathematical models to analyze the dynamic behavior of HBV. Then stability theory, bifurcation phenomenon, analysis of sensitivity and optimal control strategies of the infectious disease models are studied. This not only helps to achieve a reduction in Hepatitis B transmission, but also helps to prevent Hepatitis B in daily life. For example, in ^[1,2], the authors established the epidemic models with bilinear incidence for Hepatitis B. Hepatitis B can be effectively controlled through optimal control strategies. In ^[3], a Hepatitis B transmission model is established as follows:

where ${\alpha _1}$ represents the transmission rate of HBV, ${\alpha _2}$ represents the saturation rate. The kinetic behavior of HBV model is discussed with stability theory. In ^[4,5], the authors established the fractional HBV model. By means of the fractional Routh-Hurwitz stability criterion, the global dynamic behavior of fractional Hepatitis B model is studied.

Epidemics are strongly influenced by environmental changes. In the case of human diseases, since one person's contact with another is unpredictable, the prevalence and spread of infectious diseases is random. Therefore, it is necessary to incorporate random effects into mathematical models ^[6,7]. By doing so, we are able to develop a more reasonable model. Modeling using stochastic differential equations is a very suitable method in the study of epidemic dynamics. Many researchers use stochastic infectious disease models to investigate different diseases. Here we focus on the stochastic HBV models. Khan et al.^[8] proposed a stochastic HBV epidemic model with bilinear incidence as follows:

where $\alpha$ represents the bilinear incidence rate, $\eta^{2} > 0$ represents the intensity of white noise, $B(t)$ represents the Brownian motion. Based on the theory of stochastic Lyapunov function, the dynamic behavior of Hepatitis B stochastic model is studied. Anwarud Din et al. ^[9,10,11] proposed the stochastic models of Hepatitis B with standard incidence. Wu et al. ^[12] established a stochastic delay model of Hepatitis B with bilinear incidence. Anwarud Din et al. ^[13] built a stochastic time-delay model of Hepatitis B with standard incidence. In ^[14], a new stochastic Hepatitis B epidemic model that includes white noise, Markov switching, and vaccination control was developed. The above studies concluded that high noise can guarantee the extinction of Hepatitis B.

To sum up, most of the current research focus on HBV models with bilinear incidence, standard incidence and saturated incidence, and few studies on HBV models with environmental noise disturbance and general incidence. The general incidence rate is more realistic than the bilinear incidence rate, standard incidence rate and saturated incidence rate. This is the research motivation of this paper. New stochastic and deterministic models of Hepatitis B epidemic with general incidence rate are established. The dynamic behavior of HBV model is studied, and the influence of environmental white noise on the epidemic of HBV is analyzed, and the optimal control strategy for eliminating HBV is developed. The research work in this paper is an extension of the work on ^{[3,8,9,10,11]}. The meanings of parameters in the models studied in this paper are as follows:

$\Lambda$: the birth rate;

$\beta$: infection rate from susceptible population to Hepatitis B;

${{\mu _0}}$: the natural mortality rate;

${{\mu _1}}$: the mortality rate from HBV;

${{\gamma }}$: the recovery rate of HBV;

$\nu$: the vaccination rate of HBV.

Here is a breakdown of this article's organization. In Section 2, the new models of Hepatitis B are established. In Section 3, we obtain the basic reproduction number and the equilibrium points for the deterministic epidemic model of Hepatitis B. Lyapunov function is used to prove the local asymptotic stability. In Section 4, we verify that the stochastic model has one and only one global positive solution. The extinction, persistence and moment exponential stability of stochastic model are studied by means of stochastic Lyapunov function theory. In section 5, the optimal control strategy to eliminate HBV is developed by using the optimal control theory. To reduce Hepatitis B infection rates and to promote vaccination rates, three control variables are used, for instance, isolation of patients, treatment of patients, and vaccine inoculation. In Section 6, Runge-Kutta method is used for numerical simulation to support our main theoretical conclusions. Section 7 provides a brief summary and outlook of the main findings.

2.   Mathematical model formulation

New deterministic and stochastic mathematical models of HBV transmission are established. We make the following assumptions about the models:

$(A_{1}).$ $N(t)$ represents the total population at time $t$, which is divided into three parts: susceptible persons $S(t)$, infected persons $I(t)$ and convalescent patients $R(t)$. Namely, $N(t) = S(t) + I(t)+ R(t)$.

$(A_{2}).$ All parameter values of the models are non-negative.

$(A_{3}).$ The incidence is set as nonlinear incidence rate.

$(A_{4}).$ Once successfully vaccinated or cured by treatment, immunity is considered permanent.

The supposed conditions $(A_{1})$–$(A_{4})$ lead to Hepatitis B epidemic model as below:

with $S\left(0 \right) = S_{0} \ge 0, {I}\left(0 \right) = I_{0 }\ge 0, R\left(0 \right) = R_{0} \ge 0,$ where $\beta$ represents the transmission rate of Hepatitis B, $f\left({S, I} \right) = 1 + {a_1}S + {a_2}I + {a_3}SI$. For convenience, we define $g\left({S, I} \right) = \frac{{\beta SI}}{{f\left({S, I} \right)}} = \frac{{\beta SI}}{{1 + {a_1}S + {a_2}I + {a_3}SI}}$. Here $g\left({S, I} \right)$ is the incidence rate, where $a_{1}, \; a_{2}, \; a_{3}\geq 0$. Morbidity is the number of new infections per population in a given time period. To simulate the spread of disease, some scholars use bilinear incidence (see ^[15,16,17]), standard incidence (see ^[18,19,20]) and saturated incidence (see ^[21,22]). However, the nonlinear morbidity $g(S, I)$ takes many forms, each of which has its own advantages, as below:

(i) When ${a_1} = {a_2} = {a_3} = 0$, $g(S, I)$ represents the bilinear incidence;

(ii) When ${a_1} = {a_3} = 0$ or ${a_2} = {a_3} = 0$, $g(S, I)$ is the saturation incidence^[3];

(iii) When ${a_3} = 0$, $g(S, I)$ is the Beddington-DeAngelis functional response as shown in ^[23];

(iv) When ${a_3} = {a_1} {a_2}$, $g(S, I)$ is the Crowley-Martin functional response as shown in ^[24].

On the other hand, the following stochastic epidemic model of Hepatitis B is studied by incorporating environmental noise into the above model (2.1):

with $S\left(0 \right) = S_{0} \ge 0, {I}\left(0 \right) = I_{0 }\ge 0, R\left(0 \right) = R_{0} \ge 0,$ where ${B}(t)$ represents the Brownian motion, ${\sigma } > 0$ represents the white noise intensity.

3.   Asymptotic behavior of deterministic Hepatitis B epidemic model (2.1)

3.1.   Basic reproduction number and disease-free equilibrium point ${E^0}$

Let $\frac{{dS}}{{dt}} = 0,$ $\frac{{d{I}}}{{dt}} = 0,$ $\frac{{dR}}{{dt}} = 0,$ then ${E^0} = ({S^0}, I^0, {R^0}) = \left({\frac{\Lambda }{{{\mu _0} + \nu }}, 0, \frac{{\nu \Lambda }}{{\left({{\mu _0} + \nu } \right){\mu _0}}}} \right)$, for $I^0 = 0$. The basic reproduction number of Hepatitis B model (2.1) is given below. Let $X = (S, I, R)$, (2.1) is written as $\frac{{dX}}{{dt}} = \mathcal{F}(X) - \mathcal{V}(X),$ where

The Jacobian of equation (3.1) around ${E^0}$. ${\rm{R}}_0^d$ can be obtained by calculating the spectral radius of $\mathcal{F}\mathcal{V}^{-1}$. Then, ${\rm{R}}_0^d = \frac{{\beta \Lambda }}{{\left({{\mu _0} + \nu + {a_1}\Lambda } \right)\left({{\mu _0} + {\mu _1} + \gamma } \right)}}.$

Theorem 3.1 ${E^0}$ is locally asymptotically stable if ${\rm{R}}_0^d < 1$, otherwise ${E^0}$ is unstable.

Proof. The Jacobian matrix of Hepatitis B epidemic model (2.1) at ${E^0}$ can be calculated as follows:

After calculation, $J{|_{{E^0}}}$ has three eigenvalues:

Obviously, the sign of ${\lambda _3}$ depends on ${\rm{R}}_0^d$. If ${\rm{R}}_0^d < 1$, then the eigenvalues of $J{|_{{E^0}}}$ are all negative. Thus, HBV model (2.1) is locally asymptotically stable at ${E^0}$. On the contrary, if ${\rm{R}}_0^d > 1$, then ${\lambda _3}$ is positive, so ${E^0}$ is unstable.

3.2.   Endemic equilibrium point ${E^*}$ and stability

By simple calculation, we have ${E^*} = \left({{S^*}, I^*, {R^*}} \right)$, for $I^* \neq 0,$ where ${t_1} = {\mu _0} + {\mu _1} + \gamma, \; {t_2} = {\mu _0} + \nu,$

Lemma 3.1 If ${\rm{R}}_0^d > 1$, then ${E^*}$ exists, otherwise it does not exist.

Proof. If ${\rm{R}}_0^d > 1$, we have

And by calculation, we get $\Delta _1^2 - 4t_1^2{a_3}{\Delta _2} > 0$. Thus, ${I^*} > 0.$ Because $2\Lambda {t_1}{a_3} + {\Delta _1} > \sqrt {\Delta _1^2 - 4t_1^2{a_3}{\Delta _2}}$, so $\Lambda - \frac{{ - {\Delta _1} + \sqrt {\Delta _1^2 - 4t_1^2{a_3}{\Delta _2}} }}{{2{t_1}{a_3}}} > 0.$ Then, ${\Lambda - {t_1}{I^*}} > 0$. Hence, we have ${S^*} > 0$, thereby ${R^*} > 0$. Consequently, ${E^*}$ exists, if ${\rm{R}}_0^d > 1$.

Theorem 3.2 ${E^*}$ is locally asymptotically stable if ${\rm{R}}_0^d > 1$, otherwise ${E^*}$ is unstable.

Proof. The Jacobian matrix of Hepatitis B epidemic model (2.1) at ${E^*}$ can be calculated as follows:

By calculation, the first eigenvalue of $J{|_{{E^*}}}$ is ${\lambda _1} = - {\mu _0} < 0$. Take

By $\frac{{dI}}{{dt}} = 0$, one has $\frac{{\beta {S^*}}}{{f\left({{S^*}, {I^*}} \right)}} = \left({{\mu _0} + {\mu _1} + \gamma } \right).$ Thus, $\beta {S^*} = \left({{\mu _0} + {\mu _1} + \gamma } \right)f\left({{S^*}, {I^*}} \right)$. By the definition of $f$, we get $\beta {S^*}\left({1 + {a_1}{S^*}} \right) < \left({{\mu _0} + {\mu _1} + \gamma } \right){f^2}\left({{S^*}, {I^*}} \right)$. So,

Then, by (3.2), we obtain

and

Thus, ${\rm{trac}}\left(A \right) < 0, \; {\rm{det}}\left(A \right) > 0$ if ${\rm{R}}_0^d > 1$. According to Routh-Hurwitz criterion, the two eigenvalues of the matrix $A$ are negative. It can therefore be concluded that ${E^*}$ is locally asymptotically stable when ${\rm{R}}_0^d > 1$.

3.3.   Numerical simulation

Example 3.1 In model (2.1), let $\Lambda = 0.8, \; \beta = 0.01, \; {a_1} = 0.1, \; {a_2} = 0.2, \; {a_3} = 1.2, \; {\mu _0} = 0.02, \; {\mu _1} = 0.1, \; \nu = 0.01, \; \gamma = 0.01$, $S(0) = 200,$ ${I}\left(0 \right) = 100,$ $R(0) = 100.$ After the calculation, we get $E^{0} = (27, 0, 13)$ and ${\rm{R}}_0^d = 0.56 < 1$. By observing Figure 1(a), it can be concluded that $E^{0}$ is locally asymptotically stable, which verifies the rationality of Theorem 3.1.

Example 3.2 The values of the parameters except $\beta = 0.9$ are the same as those in Example 3.1. After calculation, we get $E^{*} = (5, 5, 10)$ and ${\rm{R}}_0^d = 50 > 1$. By observing Figure 1(b), it can be concluded that ${E^*}$ is locally asymptotically stable, which verifies the rationality of Theorem 3.2.

4.   Dynamical behavior of the stochastic HBV model (2.2)

Let $(\Omega, \mathcal{F}, \mathbb{P})$ be one complete probability space whose filtration ${\left\{ {{\mathcal{F}_t}} \right\}_{t \ge 0}}$ satisfies the usual conditions (i.e., ${\left\{ {{\mathcal{F}_t}} \right\}_{t \ge 0}}$ is monotonically increasing and right-continuous while ${\mathcal{F}_0}$ contains all $\mathbb{P}$-null sets). We think over the following $n$-dimensional stochastic differential equation:

with $x(0) = x_{0}\in \mathbb{R}^{n}$, where $f(x(t), t): {\mathbb{R}^n} \times \left[ {0, +\infty} \right] \to {\mathbb{R}^n}$ and $g(x(t), t):{\mathbb{R}^n} \times \left[ {0, +\infty} \right] \to {\mathbb{R}^{n \times m}}$ are locally Lipschitz functions in $x$. $B(t)$ represents the $n$-dimensional Brownian motion defined on $(\Omega, \mathcal{F}, {\left\{ {{\mathcal{F}_t}} \right\}_{t \ge 0}}, \mathbb{P})$. ${C^{2, 1}}({\mathbb{R}^n} \times [0, + \infty ], {\mathbb{R}_ + })$ is a family of all nonnegative functions $V(x, t)$ defined on ${\mathbb{R}^n} \times [0, + \infty ]$, making them continuously differentiable twice in $x$ and once in $t$. The differential operator $L$ associated with (4.1) is defined ^[25] by:

If $L$ acts on $V \in {C^{2, 1}}({\mathbb{R}^n} \times [{t_0}, \infty ], {\mathbb{R}_ + })$, then

If $x(t)\in \mathbb{R}^{n}$, by Itô's formula, one has

Definition 4.1 ^[26] The equilibrium point $x = 0$ of (4.1) is considered to be $p$th moment exponentially stable, if there are $C_{1}, \; C_{2} > 0$ so that

4.1.   The existence and uniqueness of positive solutions

A bounded set $\Delta$ is defined as below:

Theorem 4.1 For $\left({S\left(0 \right), {I}\left(0 \right), R\left(0 \right)} \right) \in \Delta$, model (2.2) possesses one unique positive solution $(S(t), {I}(t), R(t))$ on $t\geq 0$. The solution is still in $\mathbb{R}_ + ^3$ with probability 1, namely, $(S(t), {I}(t), R(t)) \in \mathbb{R}_ + ^3$ on $t\geq 0$ a.s.

Proof. See Appendix A.

Corollary 4.1 The set $\Delta$ is positively invariant; in other words, if $(S(0), I(0), R(0)) \in \Delta$, then $\mathbb{P} ((S(t), I(t), R(t)) \in \Delta) = 1$, $\forall t\geq 0$.

4.2.   Moment exponential stability

Lemma 4.1 ^[26] If there is the function $V(t, x) \in {C^{1, 2}}(\mathbb{R} \times {\mathbb{R}^n})$ so that

then the equilibrium point of the Eq (4.1) is $p$th moment exponentially stable. This implies that the number of infected people goes extinct at an exponential rate. When $p = 2$, it usually means that it is exponentially stable at the mean square, and global asymptotically stable at the equilibrium point $x = 0$.

Lemma 4.2 ^[27] Let $p\geq 2$, $\varepsilon, x, y > 0$, then the following two inequalities are true

Theorem 4.2 Let $p\geq 2$. If ${\rm{R}}_0^d < 1$ and

then $E^{0}$ is $p$th moment exponentially stable in $\Delta$.

Proof. See Appendix B.

Corollary 4.2 If ${\rm{R}}_0^d < 1$ and

then $E^{0}$ is globally asymptotically stable in $\Delta$.

4.3.   Extinction of the disease

In studying Hepatitis B model, we are interested in when the disease becomes extinct and when it becomes persistent in the population. This section establishes sufficient conditions for disease elimination in model (2.2). First, we need to give the basic symbols and lemmas associated with this problem. For the integrable function $x (t)$ defined on $(0, \infty)$, is denoted as $\langle {x\left(t \right)} \rangle = \frac{1}{t}\int_0^t {x\left(s \right)ds}.$ The threshold ${\rm{R}}_0^s$ of model (2.2) is as follows

Theorem 4.3 Let $\left({S\left(t \right), {I}\left(t\right), R\left(t \right)} \right)$ be the solution of model (2.2) with $\left({S\left(0 \right), {I}\left(0 \right), R\left(0 \right)} \right) \in \Delta$. Suppose (i) ${\sigma ^2} > \frac{{{\beta ^2}}}{{2\left({{\mu _0} + {\mu _1} + \gamma } \right)}}$; or (ii) ${\rm{R}}_0^s < 1$, ${\sigma ^2} \le \frac{{2\beta \left({{\mu _0} + {a_1}\Lambda } \right)}}{\Lambda }.$ Then

In other words, $I(t)$ approaches zero at an exponential rate ${\rm{a}}{\rm{.s}}{\rm{.}}$, namely, extinction is a certainty.

Proof. See Appendix C.

4.4.   Persistence

The conditions that allow the disease to persist are discussed in this section. Now, let's start with the relevant knowledge.

Definition 4.2 If $\mathop {\lim {\rm{inf}}}\limits_{t \to \infty } \langle {I\left(t \right)} \rangle > 0\quad {\rm{a}}.{\rm{s}}.$, then the disease is persistence in the mean.

Lemma 4.3 ^[28] Let $F \in C([0, \infty)\times \Omega, \mathbb{R})$, $f \in C([0, \infty)\times \Omega, (0, \infty))$ such that $\mathop {{\rm{lim}}}\limits_{t \to \infty } \frac{{F\left(t \right)}}{t} = 0\; {\rm{a}}.{\rm{s}}.$ Assume there are $\lambda_{0}, \lambda > 0$ so that for all $t\geq 0$,

Thus,

Theorem 4.4 If ${\rm{R}}_0^s > 1$, then $I(t)$ is persistence in the mean, i.e.,

Proof. See Appendix D.

4.5.   Numerical simulation

Example 4.1 In model (2.2), let $\beta = 0.01$, $\sigma = 0.4$, $S(0) = 50$, $I(0) = 20$, $R(0) = 20$. For other parameter values, see Example 3.1. Through calculation, we can get ${\sigma ^2} = 0.16 < \frac{{2\left[ {\left({{\mu _0} + {\mu _1} + \gamma } \right)\left({{\mu _0} + {a_1}\Lambda } \right) - \beta \Lambda } \right]\left({{\mu _0} + {a_1}\Lambda } \right)}}{{{\Lambda ^2}}} = 0.38125, \; {\rm{R}}_0^d = 0.56 < 1.$ Then the conditions of Corollary 4.2 are verified. Figure 2(a) illustrates that $E^{0}$ is globally asymptotically stable, which verifies Corollary 4.2.

Example 4.2 In model (2.2), let $\beta = 0.4$, $\sigma = 0.8$. For other parameter values, see Example 3.1. Through calculation, we can get ${\sigma ^2} = 0.64 > \frac{{{\beta ^2}}}{{2\left({{\mu _0} + {\mu _1} + \gamma } \right)}}\doteq 0.6154$. Then the condition (i) of Theorem 4.3 is verified. By Theorem 4.3, one has $\mathop {{\rm{limsup}}}\limits_{t \to \infty } \frac{{\ln I\left(t \right)}}{t} \le \frac{{{\beta ^2}}}{{2{\sigma ^2}}} - \left({{\mu _0} + {\mu _1} + \gamma } \right) = -0.005 < 0\; {\rm{a}}{\rm{.s}}{\rm{.}}$ As a result, $I(t)$ approaches 0 at an exponential rate with probability 1. In other words, the disease has disappeared. As a confirmation of our findings, we present simulations based on the Euler-Maruyama (EM) method, as shown in Figure 2(b).

Example 4.3 In model (2.2), let $\beta = 0.1$, $\sigma = 0.15$. For other parameter values, see Example 3.1. Through calculation, we can get ${\rm{R}}_0^s \doteq 0.615 < 1$, ${\sigma ^2} = 0.0225 \le \frac{{2\beta \left({{\mu _0} + {a_1}\Lambda } \right)}}{\Lambda } = 0.025$. Then the condition (ii) of Theorem 4.3 is verified. By Theorem 4.3, one has $\mathop {{\rm{limsup}}}\limits_{t \to \infty } \frac{{\ln I\left(t \right)}}{t} \le \left({{\rm{R}}_0^s - 1} \right)\left({{\mu _0} + {\mu _1} + \gamma } \right) \doteq -0.05005 < 0\; {\rm{a}}{\rm{.s}}{\rm{.}}$ As a result, $I(t)$ approaches 0 at an exponential rate with probability 1. In other words, the disease has disappeared. The simulation result is shown in Figure 2(c).

Example 4.4 In model (2.2), let $\beta = 0.8$, $\sigma = 0.4$. For other parameter values, see Example 3.1. Through calculation, we can get ${\rm{R}}_0^s \doteq 9.8 > 1$. Then Theorem 4.4 is verified. By Theorem 4.4, one has

Thus, the disease persists. The simulation result is shown in Figure 2(d).

5.   Optimal control analysis

In order to control the spread of HBV, this section adopts the optimal control theory^[29,30,31]. Our aim is to seek an effective control strategy to reduce HBV infection in the population. In order to reduce Hepatitis B infection rates and to promote vaccination rates, three control variables are used, for instance, $u_{1}(t)$, $u_{2}(t)$ and $u_{3}(t)$. The specific explanation is as below:

1) $u_{1}(t)$ is the isolation rate. Through this control variable, infected persons are isolated to avoid contact between infected persons and susceptible persons;

2) $u_{2}(t)$ represents the cure rate. Through this control variable, the number of patients can be decreased by using effective drugs to treat the infected persons;

3) $u_{3}(t)$ represents vaccination rate. The spread of Hepatitis B can be reduced through vaccination.

5.1.   Optimal control of Hepatitis B model (2.1)

To design a control strategy to eliminate Hepatitis B, we will consider the optimal strategy of deterministic model (2.1). The control strategy is achieved by minimizing the target function as below

subjugate to the control model

with

In Eq (5.1), $\omega_{1}, \; \omega_{2}, \; \omega_{3}, \; \omega_{4} > 0$. In the target function, $\omega_{1}$ represents the weight constant of Hepatitis B infection $I(t)$. $\omega_{2}$, $\omega_{3}$ and $\omega_{4}$ represent the weight constants of quarantine of infected persons and susceptible persons, infected persons' treatment and vaccination, respectively. $\frac{1}{2}{\omega_2}u_1^2\left(t \right)$, $\frac{1}{2}{\omega_3}u_2^2\left(t \right)$ and $\frac{1}{2}{\omega_4}u_3^2\left(t \right)$ represent costs related to segregation, treatment and vaccine inoculation, respectively. The objective of this section is to seek $u_1^*, \; u_2^*, \; u_3^*$ so that

subordinate to problems (5.2) and (5.3), where

5.1.1.   Existence of solution

As a result of this part, it will be demonstrated that the control problems (5.2) and (5.3) have a solution. Let

where

It is obvious that Eq (5.6) is a nonlinear system with bounded coefficients. Let

which satisfies

where ${K_1} = \max \left\{ {{k_i}} \right\}, i = 1, 2, 3$ is not affected by state parameters in system (5.2). And we could write it the same way for this

where ${K_2} = \max \left\{ {{K_1}, \left\| {{L_1}} \right\|} \right\} < \infty$, which implies that $F$ is continuous and uniformly Lipschitz. It goes without saying that state variables and control variables can't be negative. Thus, the solution of model (5.2) exists. As a next step, we will determine the control variables for minimizing the objective function.

Theorem 5.1 There is an optimal control ${u^*} = \left({u_1^*, u_2^*, u_3^*} \right) \in U$, so that

subjugate to systems (5.2) and (5.3).

Proof. Our method of demonstrating optimal control is the one proposed in ^[31,32]. Because the state parameters and the control variables are both positive. Therefore, in the minimization problem, the necessary convexity in $u_{1}(t)$, $u_{2}(t)$ and $u_{3}(t)$ of the objective functional defined in Eq (5.1) is satisfied. As defined, $u_{1}$, $u_{2}$, $u_{3} \in U$ are specified as enclosed and convex variables. In order for optimal control to exist, the optimal system (5.2) must have a bound, which ensures its compactness. In addition, the function is composed of control variables and state variables, which shows the convexity of the objective function. Therefore, the problem under consideration satisfies all assumptions, so $\left({u_1^*, u_2^*, u_3^*} \right)$ exists.

5.1.2.   Optimality conditions

Define $u = (u_{1}, u_{2}, u_{3})$, $x = \left({S, I, R} \right)$. Then the definition of Lagrangian $L$ is as below:

And the related Hamiltonian $H$ is defined as below:

where

with

Secondly, the main research tool for optimal solution of control problem is the standard Pontryagin maximum principle. Suppose that $\left({{x^*}, {u^*}} \right)$ is the optimal solution of (5.1)–(5.3), then there is one nontrivial vector function $\lambda$, so that

with

and the transversal condition

Theorem 5.2 For the optimal control problems (5.1) and (5.2), ${S^*}, I^*$ and ${R^*}$ are the optimal state solutions of the optimal control variables $\left({u_1^*, u_2^*, u_3^*} \right)$. Then there are adjoint variables ${\lambda _1}\left(t \right), {\lambda _2}\left(t \right)$ and ${\lambda _3}\left(t \right)$, so that

with transversality conditions

And the specific forms of ${u_1^*(t), u_2^*(t)}$ and ${u_3^*(t)}$ are as below

Proof. Taking the partial derivatives of $S(t)$, ${I}(t)$, $R(t)$ in (5.13) yield the adjoint variables (5.18). In addition, to calculate ${u_1^*, \; u_2^*, \; u_3^*}$, we take partial derivatives of $u_{1}$, $u_{2}$, $u_{3}$ in (5.13). Then we obtain the optimal control variables (5.20)–(5.22).

By giving the equation of the adjoint variables (5.18) and its related conditions (5.3) and (5.19) and the optimal control parameters (5.20)–(5.22), control variables and state variables of the control problem are solved. The next step is to investigate the control theory of stochastic systems when the same control measures are applied.

5.2.   Optimal control for HBV stochastic model (2.2)

As mentioned in Section 5.1, in this section, we study the stochastic optimal control of model (2.2), and the following stochastic control system can be obtained

with

For easy reading, the vectors are given as below:

and

with $x\left(0 \right) = {\left({S\left(0 \right), {I}\left(0 \right), R\left(0 \right)} \right)^\prime } = {x_0},$ where

The quadratic function is considered as below:

where ${C_i} > 0, \; i = \overline {1, 7}.$ The objective of this section is to obtain the control vector ${u^*}\left(t \right) = \left({u_1^*\left(t \right), u_2^*\left(t \right), u_3^*\left(t \right)} \right)$ so that

where the specific form of $U$ is as below:

where ${u_i^{\max }} > 0$. The next step is to apply the stochastic maximum principle to define the Hamiltonian $H\left({x, u, p, q} \right)$, which has the following form:

where $q = {\left({{q_1}, {q_2}, {q_3}} \right)^\prime }, \; p = {\left({{p_1}, {p_2}, {p_3}} \right)^\prime }$ are two different adjoint vectors, and $\langle { \cdot, \cdot } \rangle$ denotes Euclidean inner product space. As a result of applying the maximum principle, one has

In this case, $x^{*}(t)$ is one optimal path of $x(t)$. It is observed that, at both initial and terminal conditions of Eqs (5.32) and (5.33) are

As shown in Eq (5.34), the optimal control $x^{*}(t)$ is an operator of $q(t)$, $p(t)$ and $x^{*}(t)$. Thus, it means that

In this case, $\Phi$ can be computed by (5.34). Therefore, Equations (5.32) and (5.33) are written as the following equations

So,

According to random maximum principle,

We have

An auxiliary initial conditions and end conditions are granted as below

In this case, ${p_1}\left({{t_f}} \right) = - {k_1}S, \; {p_2}\left({{t_f}} \right) = - {k_2}{I}, \; {p_3}\left({{t_f}} \right) = - {k_3}R.$ In the Hamiltonian equation, by taking the derivatives of $u_{1}$, $u_{2}$, $u_{3}$, we can figure out that

In control theory, the desired objective is achieved by adjusting control variables. By substituting the control parameters into the system, the optimal goal of the power system can be obtained. The limits of control can be set according to Eq (5.27). We will establish the objective function by referring to the method of Eq (5.28). As a result of objective functions, there is a direct relationship between optimality and optimality. So special attention should be paid when selecting the objective function. Whenever the control objective function has multiple factors, the more important items should be given weight. Before applying the Pontryagin Maximum Rules^[33], the existence and compactness of an optimal control are tested. The optimal control makes the goal function reach a maximum or minimum value at a certain point. Differential equations can be optimized to Hamiltonian at a certain point. Hamiltonian is defined as below:

Optimal control involves finding the necessary point $u^{*}$ to maximize the Hamiltonian equation. We are able to obtain the adjoint system (5.43) by taking the derivative of H w.r.t. with respect to the state variable and substituting it with the final condition.

6.   Numerical simulations

Our analysis results are supported by approximate simulations of Hepatitis B models (2.1) and (2.2), respectively. Simulations can be performed from qualitative aspects. To test the rationality of the results, we use the stochastic Runge-Kutta method is adopted to simulate model (2.2), and the calculation model is obtained as below:

where ${\sigma } > 0$ is the white noise value, ${\xi _{k}}\left({k = \overline {1, n} } \right)$ is a standalone Gaussian stochastic variable with $N(0, 1)$ and the step length $\Delta t > 0$.

Next, the qualitative characteristics of deterministic and stochastic optimal controls are modeled. Firstly, Runge-Kutta iterative technique is used to simulate the deterministic model. The time interval is set as [0,100] in the positive direction, and the state system (5.2) and transverse condition (5.19) are solved by the prescribed method. Then the adjoint equation (5.18) with the same time interval is simulated by Runge-Kutta iterative method in the backward direction supported by transverse condition (5.19). Note that the values of other parameters except $\beta$ are shown in Example 3.1. The results are shown in Figure 3(a). Figure 3(a) shows the dynamic curves of susceptible population, infected person and recovered person in HBV deterministic model (2.1) with and without optimal control value. A clear difference can be observed between the two conditions with and without control. The simulation results show that with the implementation of control measures, the number of the infected and susceptible population tend to decrease, while the number of recovered patients tend to increase.

Next, we simulate optimal control techniques for stochastic models. The stochastic Runge-Kutta iterative technique is applied to simulate the optimal control system. Considering the transverse condition, the optimal control strategy is realized by approximating the state and adjoint model. First, we apply the stochastic Runge-Kutta iterative method to calculate the state system (5.42). Secondly, under the transverse condition (5.43), the corresponding adjoint equation (5.42) of the system is obtained by using the reverse technique and the iterative technique of the state equation. The control is then modified by applying the convex combination of the control and the values from the characterizations (5.45)–(5.47). The algorithm is repeated over and over again, and the iteration is performed until the difference between the values obtained in two successive iterations is very small. Note that the values of other parameters except $\beta, \; \sigma$ are shown in Example 3.1. The results are shown in Figure 3(b). Figure 3(b) respectively shows the dynamic curves of susceptible population, infected persons and recovered persons in the stochastic model (2.2) with and without optimal control values. Under the optimal control parameter values, the dynamic behavior before and after the optimal control is significantly different. The images show that through isolation, treatment and vaccination, the number of susceptible and infected persons is decreasing, but the number of recovered people is increasing.

7.   Conclusions

This paper investigates the transmission dynamics of HBV. The optimal control strategy is developed to control the transmission of HBV in the population. To this end, we first establish new HBV models with general incidence rate. We calculate the basic reproduction number, equilibrium points of deterministic Hepatitis B model to study the local asymptotic stability under certain conditions. Secondly, we calculate the random threshold. The random Lyapunov function theory is applied to verify that the model has one unique global positive solution. The extinction, persistence and stability of stochastic Hepatitis B model are given. These conditions are expressed as expressions containing the stochastic system parameters and the intensity of the noise term. It is clear that noise intensity has an important effect on disease transmission. To control the transmission of HBV, optimal control strategies are used to eliminate the transmission of HBV. In order to reduce Hepatitis B infection rates and to promote vaccination rates, three control variables are used, for instance, isolation of patients, treatment of patients, and vaccine inoculation. Runge-Kutta method is used for numerical simulations to support the theoretical results. It can be found that when the white noise is stronger, the extinction rate of the disease is higher. Disease is more persistent when white noise is lower in intensity. Virus dynamics-based stochastic epidemic models perform better in our study. A broad range of biomedical applications can be made from this theory, as it provides a solid foundation for studying similar diseases. An infection dynamics model based on stochastic delayed infection can, for example, be considered for studying the effects of incubation periods. Additionally, our research can be applied to analyze other epidemics, such as COIVD-19, tuberculosis, HIV and so on.

Acknowledgments

This research is funded by the Natural Science Foundation of Xinjiang Uygur Autonomous Region (Grant No. 2022D01A246, 2021D01B35, 2021D01A65), Natural Science Foundation of Colleges and Universities in Xinjiang Uygur Autonomous Region (Grant No. XJEDU2021Y048) and Doctoral Initiation Fund of Xinjiang Institute of Engineering (Grant No. 2020xgy012302).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix A.   Proof of Theorem 4.1.

Let $\left({S\left(0 \right), I\left(0 \right), R\left(0 \right)} \right) \in \Delta$, $N(t) = S(t)+I(t)+R(t)$. By model (2.2), one has

Then,

By integration, we get

So $N\left(\tau \right) < \frac{\Lambda }{{{\mu _0}}}$,

There is no doubt that model (2.2) meets the local Lipschitz condition. For $\left({S\left(0 \right), {I}\left(0 \right), R\left(0 \right)} \right) \in \mathbb{R}_ + ^3$, in this model (2.2), there is a unique local solution $(S(t), {I}(t), R(t))$, $t \in [0, {\tau _e})$. In this case, ${\tau _e}$ refers to the duration of the explosion. In order to prove that ${\tau _e} = \infty$ a.s., we must do the following. Make $k_{0} > 0$ large enough so that $\left({S\left(0 \right), {I}\left(0 \right), R\left(0 \right)} \right) > k_{0}$. For each integer $k \leq k_{0}$, the stopping time is defined as below

A $C^{2}$-function $V:\mathbb{R}_ + ^3 \to {\mathbb{R}_ + }$ is defined as below

By applying the formula Itô, we can get

For $\tau \in \left[ {0, t \wedge {\tau _k}} \right]$, one has

Thus,

where $K = 3{\mu _0} + {\mu _1} + \gamma + \nu + \frac{{\beta \Lambda }}{{{\mu _0} + {a_1}\Lambda }} + \frac{{{\sigma ^2}{\Lambda ^2}}}{{{{\left({{\mu _0} + {a_1}\Lambda } \right)}^2}}}.$ Integrate the above inequality from $0$ to $\tau_{k} \wedge t$, and then taking the expectation, according to the properties of Brownian motion, we get

Because $V\left({S\left({{\tau _k} \wedge t} \right), {I}\left({{\tau _k} \wedge t} \right), R\left({{\tau _k} \wedge t} \right)} \right) > 0$, so

For ${\tau _k}$, some component of $S\left({{\tau _k} } \right), {I}\left({{\tau _k} } \right), R\left({{\tau _k} } \right)$ is equal to $k$. Thus, $V\left({S\left({{\tau _k}} \right), I\left({{\tau _k}} \right), R\left({{\tau _k}} \right)} \right) \ge - \ln \left({\frac{{k{\mu _0}}}{\Lambda }} \right).$ So,

By (7.8), we obtain

Extending $k$ to 0, one has $\mathbb{P} \left({{\tau _0} \le t} \right) = 0$, $t > 0$. Therefore, $\mathbb{P}({\tau _0} = \infty) = 1$. Consequently, ${\tau _0} = {\tau _e} = \infty\; \; {\rm{a}}{\rm{.s}}.$

Appendix B.   Proof of Theorem 4.2.

Let $p\geq 2$. The Lyapunov function is considered as below

Here ${\tau _1} > 0$ will be determined later. The first inequality in (4.6) is easily proved to be true. Then,

In $\Delta$, we get

According to Lemma 4.2, we can get

Then

Select $\varepsilon$ small enough to make the coefficient of ${\left({\frac{\Lambda }{{{\mu _0} + \nu }} - S} \right)^p}$ negative. By (4.8), we have ${ - \left({{\mu _0} + {\mu _1} + \gamma } \right) + \frac{{\beta \Lambda }}{{{\mu _0} + {a_1}\Lambda }} + \frac{{\left({p - 1} \right){\sigma ^2}{\Lambda ^2}}}{{2{{\left({{\mu _0} + {a_1}\Lambda } \right)}^2}}}} < 0$. When ${ - \left({{\mu _0} + {\mu _1} + \gamma } \right) + \frac{{\beta \Lambda }}{{{\mu _0} + {a_1}\Lambda }} + \frac{{\left({p - 1} \right){\sigma ^2}{\Lambda ^2}}}{{2{{\left({{\mu _0} + {a_1}\Lambda } \right)}^2}}}} < 0$, select $\tau_{1}$ is positive, so that the coefficient of $I^{p}$ is negative.

Appendix C.   Proof of Theorem 4.3.

Applying Itô formula for the second equation of model (2.2), one has

Integrating the above equation from 0 to $t$, and dividing both sides by $t$ at the same time, then we get

where $M(t): = \int_0^t {\frac{{\sigma S\left(\tau \right)}}{{f\left({S\left(\tau \right), I\left(\tau \right)} \right)}}dB\left(\tau \right)}.$ Then

The following formula can be obtained from the martingale theorem of large numbers

If condition (i) is met, then by (7.17), (7.18), we have

From the definition of $f(S, I)$, we get

If condition (ii) is met, then by (7.20), one has

If condition (ii) is met, then by (7.21), one has

The above inequality indicates that

Appendix D.   Proof of Theorem 4.4.

By Corollary 4.1, we have $S(t) < \frac{\Lambda }{{{\mu _0}}}\; {\rm{a}}{\rm{.s}}{\rm{.}}$ Then,

By model (2.2), one has

Integrating the above equation from 0 to $t$, and dividing both sides by $t$ at the same time, then

Then

here $\varphi \left(t \right) = {{\left[ {S(t) + I(t) + R(t) - \left({S(0) + I(0) + R(0)} \right)} \right]} / {\left({{\mu _0}t} \right)}}$. In this situation,

Applying Itô formula for the second equation of model (2.2) and combining with (7.20), (7.24), we get

Integrating the above equation from 0 to $t$, and combining with (7.27), then

here $\Phi \left(t \right) = \ln I\left(0 \right) + \int_0^t {\frac{{\sigma S(\tau)}}{{f\left({S(\tau), I(\tau)} \right)}}dB(\tau)} - \frac{{\beta {\mu _0}}}{{{\mu _0} + {a_1}\Lambda }}\left({\langle {R(t)} \rangle + \varphi \left(t \right)} \right)t.$ By combining the martingale theorem of large numbers with (7.28), we can get

Lemma 4.3 implies that