Dynamics of a stochastic hepatitis B virus transmission model with media coverage and a case study of China

1.   Introduction

Hepatitis B, a viral liver infection caused by the hepatitis B virus (HBV), is a major threat to public health and security around the whole world. HBV can cause both acute and chronic hepatitis. Acute hepatitis B refers to the virus infection less than half a year and part of the patients can recover and get lifetime immunity. The disease course of chronic hepatitis B can be quite long, and it can lead to severe liver disease such as cirrhosis and liver cancer ^[1]. Globally $257$ million people were living with chronic hepatitis B infection, which resulted in estimated $887, 000$ deaths in $2015$. According to the latest fact sheets released by the World Health Organization (WHO), $296$ million people were living with chronic hepatitis B infection in $2019$, and there are $1.5$ million new infections each year ^[2]. Worldwide, hepatitis B resulted in an estimated $820, 000$ deaths in $2019$, mostly from cirrhosis and hepatocellular carcinoma (primary liver cancer) ^[2]. Hence, hepatitis B epidemic is still a major global health problem.

In China, hepatitis B is one of the top three infectious diseases reported by the Chinese Center for Disease Control and Prevention (CDC) ^[3]. Based on the sero-epidemiological investigation of hepatitis B in $1992$, about $9.75\%$ of the general population in China was chronic HBV carriers. That is, around $130$ million people of mainland China were carriers of HBV in the $1990$s ^[1]. Then an effective and national vaccination program has been conducted by the government, especially the hepatitis B planned immunization for newborn babies and children. The sero-epidemiological survey in $2006$ showed about $7.18\%$ of the population was hepatitis B carriers, and the number of hepatitis B carriers was reduced by $30$ million ^[3]. Lately in $2014$, an epidemiology survey of people younger than $29$ reported that the HBV infection rates were $0.32\%$ among the group of children aged one to four, $0.94\%$ in the group aged $5$ to $14$, and $4.18\%$ in the group aged $15$ to $29$, respectively. The government of mainland China has achieved remarkable results in the prevention and treatment of hepatitis B. Nevertheless, the epidemic situation remains to be severe and complicated. At present there are $80$ million HBV carriers in mainland China, among whom $28$ million carriers need immediate medical treatment. During the past decade, the number of reported incidence of hepatitis B is around one million each year (Figure 1(a)). The incidence rates of hepatitis B in mainland China from $2005$ to $2021$ are presented in Figure 1(b), and the incidence rates varied from $89.00$ per $100, 000$ in $2007$ (the highest during $2005$–$2021$) to $64.29$ per $100, 000$ in $2020$ (the lowest during $2005$–$2021$).

The incidence of hepatitis B is different among regions, and it could be influenced by various factors such as vaccination and economy. Figure 2 displays the hepatitis B incidence rates in $31$ provinces and municipalities of mainland China, and the data in $2007$ and $2020$ are chosen and displayed. The figure shows that the incidence decreased for most provinces ($22$ of $31$), especially for Ningxia, Gansu, Henan, Xinjiang and Qinghai. This strongly reveals that the immunization program with hepatitis B vaccine is successful in most provinces, and the achievement is remarkable for western areas in China. However, it should be pointed out, in $2020$, the incidence varied from $6.15$ per $100, 000$ in Beijing to $150.99$ per $100, 000$ in Qinghai. There is still a gap between the impoverished west and the prosperous east.

In the last three decades, the research on hepatitis B has attracted more and more attention from the perspective of epidemiology and biomathematics. For instance, in $1991$, Anderson and May used a simple mathematical model to study the impacts of carriers on the transmission of HBV ^[4]. Zhao et al. ^[5] studied the dynamics of HBV transmission and proposed vaccination strategy to prevent the prevalence in the population. In $2010$, Zou et al. ^[1] formulated a high dimensional deterministic model to analyze the transmission dynamics of HBV infection in China. Then in $2015$ Zou et al. ^[6] investigated the sexual transmission dynamics of HBV in China. Zhang et al. ^[7] proposed a four dimensional HBV epidemic model involving an exponential birth rate and vertical transmission, and fitted the parameters to the public data in Xinjiang, China. Recently Din and Li ^[8] analyzed a stochastic hepatitis B epidemic model with vaccination effect and showed a case study for Pakistan. For more mathematical epidemic models on HBV infection, one can refer to ^[9,10,11,12] and the references therein.

Besides the classical compartmental models for hepatitis B epidemic, a vibrant field of research on the virus transmission dynamics has emerged and developed rapidly ^[12]. The modeling approaches usually incorporate the interaction among intracellular viral dynamics, multicellular infection process, and immune responses ^{[13,14,15,16]}. For instance, Wang et al. ^[13] formulated a multi-scale computational model of SARS-CoV-2 infection using a combination of differential equations and stochastic modeling, and revealed heterogeneity among COVID-19 patients. Rihan and Alsakaji ^[14] analyzed a stochastic delay HBV infection model with cell-to-cell transmission and cytotoxic T lymphocytes (CTLs) immune response. Recently Mohajerani et al. ^[17] proposed a multiscale modeling of HBV capsid assembly pathways, by constructing Markov state models and employing transition path theory.

When an infectious disease breaks out in one area, the center for disease control and prevention will release authoritative information through the mass media at the first moment, including the potential risk of the disease and prevention measures ^[18]. Media coverage helps to raise public health awareness, increase vaccine rates and reduce the spread of infections. Recently, a number of epidemic models have taken into account the influence of media coverage on the spread of infectious disease ^{[19,20,21,22,23]}. Their study suggests that media coverage is critical in disease eradication ^[19,20]. There are different forms of function to represent the effect of media coverage in epidemic models. Let $S$ and $I$ denote the number of the population that are susceptible and infectious respectively, and let $\beta$ be the transmission rate. Cui et al. ^[19] used the exponential function $\beta exp(-mI)SI$ to denote the incidence rate with media coverage, where $m > 0$. In another work ^[20], the incidence rate incorporating media coverage takes the form: $(\beta-\beta_2f(I))SI$, where $\beta > \beta_2$ and the function $f(I)$ satisfies $f(0) = 0$, $f'(I)\geq 0$, $\lim_{I\rightarrow +\infty} f(I) = 1$. In the present paper, similar to the approach in the previous literatures ^[21,24,25], we adopt the most common form $f(I) = \frac{I}{b+I}$, where $b$ is a half saturation constant.

Furthermore, the transmission of HBV is disturbed by various random factors in the environment, such as population mobility and unpredictable exposure to infections. An increasing number of researchers have formulated stochastic hepatitis B epidemic models considering environmental noise ^[8,10,11]. In fact, environment disturbances have an important effect on the evolution of infectious diseases ^{[26,27,28,29,30]}, and Gaussian white noise is usually selected as an appropriate representation of environmental fluctuations ^[31]. For instance, Wang et al. ^[28] and Meng et al. ^[29] proved that a large disturbance of white noise can lead infectious diseases to extinction. A large number of works also indicate that stochastic disturbance can suppress disease outbreak ^[18,26].

In the present paper, motivated by the above discussion, we formulate and investigate a stochastic hepatitis B model with media coverage and saturated incidence rate. We obtain a sufficient condition for the extinction of the disease, and prove that the stochastic system has a unique stationary distribution under certain conditions. Moreover, as a case study, we utilize our model for fitting the available data of mainland China from 2005 to 2021. Based on our simulation result, the incidence rate of hepatitis B in China will remain around 50–60 per $100, 000$ in the long term.

The paper is organized as follows. In Section $2$, we present the model formulation. In Section $3$, the existence and uniqueness of positive solution is proved for the stochastic model. In Section $4$, we obtain a sufficient condition for the extinction of the disease. In Section $5$, the sufficient condition on the existence of stationary distribution is obtained, which indicates that all the compartments will be persistent. Moreover, we provide an estimation of lower bound of the expectation for the number of infected cases. In Section $6$, numerical simulations are conducted to illustrate our theoretical results. As a case study, we fit our model to the available HBV data of mainland China from 2005 to 2021.

2.   Model formulation

Recently, Khan et al. ^[10] proposed a stochastic model for the transmission of HBV. In their work, the population is divided into four compartments: susceptible humans $S$; acute HBV infections $I_1$; chronical HBV infections $I_2$; and recovered population $R$. They further assumed: (i) the contact of susceptible individuals with acutely and chronically infected hepatitis B individuals primarily causes acutely infected species; (ii) the transmission coefficient $\beta$ is subject to random fluctuation, that is, $\beta_i\rightarrow \beta_i+\eta_i\dot{B}_i(t)$ for $i = 1, 2$, where $B_i(t)$ is standard Brownian motion with $B_i(0) = 0$ and with the noise intensity $\eta_i^2 > 0$. Then the stochastic hepatitis B epidemic model in ^[10] was presented as follows

where $\Lambda$ is the recruitment rate of the population, $\mu_0$ is the natural mortality rate, $\mu_1$ is the disease mortality rate, and $\upsilon$ represents the vaccination rate of hepatitis B. Moreover, $\beta_i$ ($i = 1, 2$) represent the transmission rate of hepatitis B, $\gamma$ is the moving rate of acutely infected humans to chronic stage, and $\gamma_i$ ($i = 1, 2$) denote the recovery rates of acutely and chronically infected hepatitis B individuals, respectively.

In the above model, the authors chose bilinear incidence rate, that is, $\beta SI$. As a matter of fact, the transmission of infectious diseases is complicated, and nonlinear incidence rates have been wildly utilized in epidemic model ^[26,29]. In $1978$, to describe the phenomenon that incidence is increasing and the population is saturated with the infective, Capasso and Serio ^[32] proposed a saturated incidence rate $\frac{\beta SI}{1+\alpha I}$. Then such saturated incidence has been extensively used in epidemic models. For instance, Khan and Zaman ^[33] and Liu et al. ^[11] have formulated hepatitis B epidemic models with saturated incidence rate. In addition, as discussed in the introduction part, we also incorporate the impact of media coverage by using the function $f(I) = \frac{I}{b+I}$. Consequently, based on the deterministic part of model (2.1), we obtain the hepatitis B model with media coverage and saturated incidence rate as follows

where the forms $\frac{SI_i}{1+\alpha_i I_i}$ ($i = 1, 2$) represent the saturated incidence rates of acute and chronical infections, and the forms $\frac{I_i}{b_i+I_i}$ ($i = 1, 2$) denote the media coverage functions. Moreover, $\beta_{i1}$ ($i = 1, 2$) represent the transmission rates for acute and chronical infections, and $\beta_{i2}$ ($i = 1, 2$) denote the maximal reduced contact rates by mass media alert for acute and chronical infections, respectively. Recall that $\beta_{i1} > \beta_{i2}$ ($i = 1, 2$).

Furthermore, the disease transmission and biological populations are inevitably affected by environment noises. There are different approaches to add random perturbations to biological systems. The transmission rates $\beta_i$ ($i = 1, 2$) are assumed to be disturbed by Gaussian white noise in model (2.1). In this article, following the approach in ^[8,11,34,35], we assume that the environmental white noise is proportional to each state variable $S(t), I_1(t), I_2(t)$ and $R(t)$. Finally we extend model (2.2) to the following stochastic model

where $B_i(t)$ $(i = 1, 2, 3, 4)$ are independent standard Brownian motions with $B_i(0) = 0$, and $\sigma_i^2 > 0$ $(i = 1, 2, 3, 4)$ denote the intensities of white noise.

Throughout the paper, let $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t\geq0}, \mathbb{P})$ be a complete probability space with a filtration $\{\mathcal{F}_t\}_{t\geq0}$ satisfying the usual conditions (i.e., it is increasing and right continuous while $\mathcal{F}_0$ contains all $\mathbb{P}$-null sets), then $B_i(t)$ $(i = 1, 2, 3, 4)$ are defined on this complete probability space. We also introduce the following notations: $\mathbb{R}^d_+ = \{(x_1, x_2, ..., x_d)\in \mathbb{R}^d: x_i > 0, i = 1, 2, ..., d\}$, $a\wedge b = \min\{a, b\}$, $a\vee b = \max\{a, b\}$, $\langle f\rangle = \frac{1}{t}\int_0^t f(r)dr$.

3.   Existence and uniqueness of the positive solution

Since $S, I_1, I_2$ and $R$ in system (2.3) denote the number of individuals, they should be nonnegative from the viewpoint of biology. We first introduce some basic definitions that will be used in the reminder of the article ^[36]. In general, let $X(t)$ be a homogeneous Markov process in the $d$-dimension Euclidean space $\mathbb{R}^d$ described by the stochastic differential equation

then the diffusion matrix is defined as

Furthermore, the differential operator $L$ is defined by

where $V(x)$ is an arbitrary twice continuously differential real-value function.

In this section, we obtain the following theorem which guarantees the existence and uniqueness of positive solution for system (2.3).

Theorem 3.1. For any initial value $(S(0), I_1(0), I_2(0), R(0))\in \mathbb{R}^4_+$, there is a unique positive solution $(S(t), I_1(t), I_2(t), R(t))$ for system (2.3) on $t\geq0$ and the solution will remain in $\mathbb{R}^4_+$ with probability one, namely, $(S(t), I_1(t), I_2(t), R(t))\in \mathbb{R}^4_+$ for all $t\geq0$ almost surely.

Proof. Since the coefficients of system (2.3) are locally Lipschitz continuous in $\mathbb{R}^4_+$, then for any initial value $(S(0), I_1(0), I_2(0), R(0))\in \mathbb{R}^4_+$, there exists a unique positive solution $(S(t), I_1(t), I_2(t), R(t))$ on $t\in[0, \tau_e)$, where $\tau_e$ is the explosion time ^[37]. Thus, it suffices to verify $S(t)$, $I_1(t)$, $I_2(t)$ and $R(t)$ do not explode to infinity in a finite time, that is, $\tau_e = \infty$ a.s. Let $k_0 > 0$ be sufficiently large such that $S(0)$, $I_1(0)$, $I_2(0)$ and $R(0)$ all lie within the interval $[\frac{1}{k_0}, k_0]$. For each integer $k\geq k_0$, define the stopping time

and throughout this paper we set $\inf\varnothing = \infty$ (as usual $\varnothing$ is the empty set). Clearly, $\tau_k$ is increasing as $k\rightarrow \infty$. Let $\tau_\infty = \lim\limits_{k\rightarrow \infty}\tau_k$, then $\tau_\infty\leq\tau_e$ a.s. Thus, $\tau_\infty = \infty$ a.s. implies $\tau_e = \infty$ a.s. Now we state that $\tau_\infty = \infty$. If this assertion is false, then there is a pair of constants $T > 0$ and $\epsilon\in(0, 1)$ such that

Thus there exists an integer $k_1\geq k_0$ such that

Define a $C^2$-function $V:\mathbb{R}^4_+ \rightarrow \mathbb{R}_+$ by

where $a$ and $b$ are positive constants to be determined later. The nonnegativity of this function can be obtained from

Apply Itô's formula ^[37] to $V$, and we obtain

where

Choose $a = \min\{\frac{\mu_0}{\beta_{11}}, \frac{\mu_0+\mu_1}{\beta_{21}}\}$ and $b = \frac{\mu_0 b_1\alpha_1}{\beta_{12}}$, then

and

where $K$ is a positive number. Thus, according to Eq (3.3), one can get

Let $k\geq k_1$ and integrate the above inequality from $0$ to $\tau_k\wedge T$, then we have

Taking the expectation of both sides of Eq (3.4) yields

Set $\Omega_k = \{\omega\in\Omega:\tau_k = \tau_k(\omega)\leq T\}$ for $k\geq k_1$, then we have $\mathbb{P}(\Omega_k)\geq\epsilon$ by $(3.2)$. Note that for every $\omega\in\Omega_k$, at least one of $S(\tau_k, \omega), I_1(\tau_k, \omega), I_2(\tau_k, \omega)$ and $R(\tau_k, \omega)$ equals to either $k$ or $\frac{1}{k}$, So $V(S(\tau_k, \omega), I_1(\tau_k, \omega), I_2(\tau_k, \omega), R(\tau_k, \omega))$ is no less than either

Thus, one can obtain

It then follows from Eq (3.5) that

where $I_{\Omega_k(\omega)}$ is the indicator function of $\Omega_k$. Taking $k\rightarrow \infty$ will induce $\infty > V(S(0), I_1(0), I_2(0), R(0))+KT\geq+\infty$, and it is a contradiction. Hence, we must have $\tau_\infty = \infty$ a.s.

The conclusion is confirmed.

4.   Extinction

In this section, we mainly discuss the extinction of hepatitis B under some conditions. To begin with, we present the following theorem.

Theorem 4.1. If $\min\{\mu_0, \mu_1\} > \frac{(\sigma_1^2\vee\sigma_2^2\vee\sigma_3^2\vee\sigma_4^2)}{2}$, then the solution $(S(t), I_2(t), S_h(t), I_h(t))$ of system (2.3) has the following properties

and

The proof is similar to those in ^[38] and hence is omitted here.

It is easy to compute that the deterministic system (2.2) admits a disease-free equilibrium point $E_0(\frac{\Lambda}{\mu_0+\upsilon}, 0, 0, \frac{\Lambda\upsilon}{\mu_0(\mu_0+\upsilon)})$. For the stochastic system (2.3), we obtain a sufficient condition on the extinction of disease in the following theorem.

Theorem 4.2. Let $(S(t), I_1(t), I_2(t), R(t))$ be a solution of system (2.3) with initial value $(S(0)$, $I_1(0)$, $I_2(0)$, $R(0))\in \mathbb{R}_+^4$. If $R_0^s: = \frac{4\Lambda(2\beta_{11}-\beta_{12}+2\beta_{21}-\beta_{22})}{\mu_0(\sigma_2^2\wedge\sigma_3^2)} < 1$ and $\min\{\mu_0, \mu_1\} > \frac{(\sigma_1^2\vee\sigma_2^2\vee\sigma_3^2\vee\sigma_4^2)}{2}$ hold, then the disease will tend to extinction with probability one, and the solution of system (2.3) satisfies

Proof. Denote

and applying Itô's formula to $\ln Q$ yields

Integrate inequality (4.1) from $0$ to $t$ and divide by $t$ on both sides, then we get

According to system (2.3), we have

Integrate the above inequality form $0$ to $t$ and divide by $t$ on both sides, then

Furthermore, it follows from Theorem 4.1 that

Therefore,

Take the superior limit on both sides of inequality (4.2), then we have

Consequently,

which implies that

On the other hand, according to Eq (4.3), we have

According to Theorem 4.1 and Eq (4.7), we obtain

The conclusion is confirmed.

Remark 1. According to the expression of $R_0^s$, it is clear that the value of $R_0^s$ decreases with the increase of $\sigma_{2, 3}$, $\beta_{12}$ and $\beta_{22}$. Since $\beta_{i2}$ ($i = 1, 2$) represent the maximal reduced contact rates by mass media, the media coverage plays an important role in disease eradication.

5.   Stationary distribution and persistence

In this section, we focus on the existence of stationary distribution for system (2.3). From the biological point of view, stationary distribution can be interpreted as a weak stability of the system, and the disease will be persistent in the time mean sense. We first present a fundamental lemma.

Lemma 5.1. (^[39]) The Markov process X(t), the solution of system (3.1), has a unique ergodic stationary distribution $\pi(\cdot)$, if there exists a bounded domain $D\subset \mathbb{R}^d$ with regular boundary $\Gamma$ and

(C.1) there is a positive number M such that $\sum_{i, j = 1}^da_{ij}(x)\xi_i\xi_j\geq M\left|\xi\right|^2$, $x\in D$, $\xi\in \mathbb{R}^d$,

(C.2) there exists a nonnegative $C^2$-function $\text{V}$ such that $\text{LV}$ is negative for any $x\in\mathbb{R}^d\backslash D$. Then

for all $x\in \mathbb{R}^d$, where $f(\cdot)$ is a function integrable with respect to the measure $\pi$.

Define

Theorem 5.2. If $\hat{R}_0^s > 1$, then there exists a unique stationary distribution for system (2.3) and it has the ergodic property.

Proof. We will prove the theorem by verifying the conditions in Lemma 5.1. The diffusion matrix of model (2.3) is given by

Apparently, the matrix $A$ is positive definite for any compact subset of $\mathbb{R}_+^4$, so the condition $(C.1)$ in Lemma 5.1 holds.

Define $V_1 = -\ln S,$ then

Define

where $a_1, a_2, c_1, c_2$ and $c_3$ are positive constants to be chosen later. By the use of Itô's formula, we obtain

Choose

and

It follows from inequality (5.1), that

Define

then by inequality (5.2) one can derive

where

and

Define

where $\theta$ is a positive constant that is less than one. Thus, we obtain

and

where

Define a $C^2$-function $\tilde{V}:\mathbb{R}_+^4\rightarrow \mathbb{R}$ by

Choose a positive constant $M$ such that

where

It is easy to check that

here, $D_n = (\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n)\times(\frac{1}{n}, n)$. Thus one can see that there exists at least one minimum point $(S^*, I_1^*, I_2^*, R^*)$ for the function $\tilde{V}(S, I_1, I_2, R)$. Define a nonnegative $C^2$-function $V:\mathbb{R}_+^4\rightarrow \mathbb{R}_+$ by

According to inequalities (5.3)–(5.7), we obtain

Now we denote a bounded closed set as follows

where $\epsilon > 0$, and we choose $\epsilon$ sufficiently small such that the following conditions hold in the set $\mathbb{R}_+^4\setminus D_\epsilon$

where $D, F, G, H, J$ are all positive constants which will be determined later. Hence, $\mathbb{R}_+^4\setminus D_\epsilon$ can be divided into the following ten domains,

Obviously, $\mathbb{R}_+^4\setminus D_\epsilon = D_\epsilon^1\bigcup D_\epsilon^2\bigcup D_\epsilon^3\bigcup D_\epsilon^4\bigcup D_\epsilon^5\bigcup D_\epsilon^6\bigcup D_\epsilon^7\bigcup D_\epsilon^8\bigcup D_\epsilon^9\bigcup D_\epsilon^{10}$. Our next task is to verify $LV\leq-1$ on $\mathbb{R}_+^4\setminus D_\epsilon$.

Case $1$. If $(S, I_1, I_2, R)\in D_\epsilon^1$, then

where

According to inequality (5.9), we obtain that $LV\leq-1$ for all $(S, I_1, I_2, R)\in D_\epsilon^1$.

Case $2$. If $(S, I_1, I_2, R)\in D_\epsilon^2$, then

It then follows from inequality (5.9) that $LV\leq-1$ for all $(S, I_1, I_2, R)\in D_\epsilon^2$.

Case $3$. If $(S, I_1, I_2, R)\in D_\epsilon^3$, then

Combining inequalities (5.8) and (5.10) we derive that $LV\leq-1$ for all $(S, I_1, I_2, R)\in D_\epsilon^3$.

Case $4$. If $(S, I_1, I_2, R)\in D_\epsilon^4$, then

By inequality (5.9), one can derive that $LV\leq-1$ for all $(S, I_1, I_2, R)\in D_\epsilon^4$.

Case $5$. If $(S, I_1, I_2, R)\in D_\epsilon^5$, then

Combining with inequality (5.9), we obtain that $LV\leq-1$ for all $(S, I_1, I_2, R)\in D_\epsilon^5$.

Case $6$. If $(S, I_1, I_2, R)\in D_\epsilon^6$, we have

According to inequality (5.9), we can deduce that $LV\leq-1$ for all $(S, I_1, I_2, R)\in D_\epsilon^6$.

Case $7$. If $(S, I_1, I_2, R)\in D_\epsilon^7$, we have

where

Combining with inequality (5.11), we can derive that $LV\leq-1$ for all $(S, I_1, I_2, R)\in D_\epsilon^7$.

Case $8$. If $(S, I_1, I_2, R)\in D_\epsilon^8$, then

where

By virtue of inequality (5.12), we obtain that $LV\leq-1$ for all $(S, I_1, I_2, R)\in D_\epsilon^8$.

Case $9$. If $(S, I_1, I_2, R)\in D_\epsilon^9$, then

where

By virtue of inequality (5.13), we can conclude that $LV\leq-1$ for all $(S, I_1, I_2, R)\in D_\epsilon^9$.

Case $10$. If $(S, I_1, I_2, R)\in D_\epsilon^{10}$, then

where

Together with inequality (5.14), one can obtain that $LV\leq-1$ for all $(S, I_1, I_2, R)\in D_\epsilon^{10}$.

Combining inequalities (5.15)–(5.24), we finally get a sufficiently small $\epsilon$ such that $LV\leq-1$ for all $(S, I_1, I_2, R)\in \mathbb{R}_+^4\backslash D_\epsilon$. Therefore the condition $(C.2)$ of Lemma 5.1 holds. According to Lemma 5.1, system (2.3) has a unique stationary distribution and it has the ergodic property.

The conclusion is confirmed.

Consequently, we provide an estimation of lower bound for the expectation of infective population in the following theorem.

Theorem 5.3. Let $(S(t), I_1(t), I_2(t), R(t))$ be a solution of system (2.3) for any initial value $(S(0)$, $I_1(0)$, $I_2(0)$, $R(0))\in \mathbb{R}_+^4$. If $\hat{R}_0^s > 1$, then

where $\lambda_2 = \max\{\frac{c_3\alpha_2\gamma}{\mu_0+\mu_1+\gamma_2}+\frac{a_2\alpha_1\gamma}{\mu_0+\gamma+\gamma_1}, \frac{a_2\alpha_1\gamma_2}{\mu_0+\gamma+\gamma_1}\}$.

Proof. Recall the function $V_2$ in the proof of Theorem 5.2

Apply Itô's formula to $V_2$, then by inequality (5.2) we have

Integrate the inequality (5.26) from $0$ to $t$ and divide by $t$ on both sides, then we obtain

where

According to the strong law of large numbers ^[37] and Theorem 4.1, it then follows that

Therefore,

namely, $\liminf\limits_{t\rightarrow \infty}\langle I_1+I_2\rangle\geq\frac{(\mu_0+\gamma+\gamma_1+\frac{\sigma_2^2}{2})(\hat{R}_0^s-1)}{\lambda_2}\quad a.s.$

Remark 2. From the viewpoint of biology, the existence of stationary distribution indicates that all the compartments will be persistent in the time mean sense. As is discussed in ^[40], the lower bound for the expectation of infective populations in Theorem 5.3 clearly shows that the disease will prevail if $\hat{R}_0^s > 1$.

6.   Numerical simulations and a case study of China

In this section, we provide some numerical simulations for system (2.3) to illustrate the feasibility of our theoretical results. Applying Milstein's higher order method ^[41] to system (2.3), we obtain the corresponding discretization equation as follows

where $\Delta t > 0$, and $\xi_k, \eta_k, \zeta_k, \varsigma_k$ ($k = 1, 2, ...n$) are independent Gaussian random variables $N(0, 1)$, and $\sigma_i^2 > 0$ ($i = 1, 2, 3, 4$) are the intensities of white noise.

First of all, take the data of hepatitis B of mainland China as a case study. We have provided the reported incidence rates (1/100, 000) of HBV in mainland China during $2005$–$2021$ in Figure 1(b). Besides, the incidence rates of hepatitis B in $31$ provinces are displayed in Figure 2. By comparing Figure 1(c) in ^[1] and Figure 2 (the data of $2007$), one can see that the reported incidence rates of HBV were taken as acute incidence rates therein. It is reasonable because the clinical difference between acute and chronic HBV infections depends on the length of infection, and acute hepatitis B usually refers to the virus infection less than six months. Therefore, inspired by the method in ^[1], we simulate the reported incidence rates in Figure 1(b) by computing the percentage of acute infection. In the process of model fitting, we first fix the initial values such that the percentage of acute infection is close to the incidence rate 75/100, 000. Then similar to the approaches in Table 3 of ^[8], we fix a relatively large value of the parameter $\Lambda$. The rest of the parameters are estimated and selected partially according to the values in Table 1 of ^[1]. More specifically, we choose the initial value $(S(0), I_1(0), I_2(0), R(0))$ as $(5000, 0.4, 20, 10)$, and the parameters in model (2.3) are taken as

Then we compute the percentage of acute HBV infection by Eq (6.1), and compare it with the incidence rates of HBV in mainland China during $2005$–$2021$. The simulation is displayed in Figure 3(a), which shows that stochastic model fits the data well by selecting appropriate parameter values. The long-term solution is given in Figure 3(b), and one can see that the incidence rate of HBV in China will remain around 50–60 per $100, 000$ in the long term.

More simulations are conducted to illustrate our theoretical results. Firstly, let the initial value $(S(0), I_1(0), I_2(0), R(0)) = (0.9, 0.4, 0.2, 0.1)$, and we choose the parameter values in the stochastic model (2.3) as follows

It is easy to compute that

Thus, the condition of Theorem 4.2 is satisfied and the disease will tend to extinction. The numerical simulation is depicted in Figure 4.

Then we fix the same initial value $(S(0), I_1(0), I_2(0), R(0)) = (0.9, 0.4, 0.2, 0.1)$, and the parameter values in the stochastic model (2.3) are taken as

One can compute that

Therefore, the condition of Theorem 5.2 is satisfied, and system (2.3) has a stationary distribution. The numerical simulation is shown in Figures 5 and 6, which indicates that the system will be persistent in mean and the disease will prevail.

7.   Conclusions and discussion

Since most realistic systems are disturbed by various stochastic factors, in this paper, we have investigated a stochastic HBV transmission model with media coverage and saturated incidence rate. To begin with, we prove the existence and uniqueness of global positive solution of system (2.3). Then we prove that hepatitis B will tend to extinction if $R_0^s < 1$ and $\min\{\mu_0, \mu_1\} > \frac{(\sigma_1^2\vee\sigma_2^2\vee\sigma_3^2\vee\sigma_4^2)}{2}$, while system (2.3) has a unique ergodic stationary distribution if $\hat{R}_0^s > 1$. According to the expression of $R_0^s$ and $\hat{R}_0^s$, the noise intensities and mass media alert are crucial factors in the disease control. We also obtain an estimation of lower bound of the expectation for the number of infectious cases. Moreover, as a case study, we fit our stochastic model to the data of reported incidence rates of HBV in mainland China, and it is anticipated that the incidence rate of HBV in China will remain around 50–60 per $100, 000$ for a long time to come.

It should be mentioned that the present HBV transmission model is formulated from a standard SIR epidemic model. This modeling method often divides the total population into several compartments under the assumption that individuals are homogeneous ^[42]. In other words, our modeling approach is conducted on the single spatial scale, without the consideration of human behavior, contact heterogeneity, and population spatial or social structure ^[12]. In recent years, multiscale models have been introduced to improve the modeling of disease transmission ^[43]. Guo et al. ^[42] developed a heterogeneous graph modeling approach to describe the dynamic process of influenza virus transmission. Since the outbreak of coronavirus disease 2019 (COVID-19) throughout the world, the modeling of SARS-CoV-2 dynamics has further motivated the trends on multiscale modeling ^[44,45], from the small scale of the virus itself and cells to the large scale of individuals and further up to the collective behavior of populations ^[46]. For instance, Hayden et al. ^[47] extended a classical SIR model to a SIRC model by considering the coronavirus concentration in the air (denoted by C). The researchers proposed multi-scale epidemic models by linking the disease transmission to information dissemination dynamics ^[48] and to the behavior change dynamics ^[49]. Such multiscale modeling approaches provide important insights into HBV transmission dynamics.

Some interesting research topics deserve further consideration. In stochastic epidemic modeling, Gaussian white noise has been usually adopted to represent environment disturbances and to reflect the fluctuations of disease transmission. Meanwhile different types of noise has also been investigated in the literatures ^{[50,51,52,53,54]}. For instance, Lévy noise is introduced to represent some abrupt environmental shocks and disasters ^[53], and telephone noise (also known as telegraph noise or burst noise) can be regarded as instantaneous transitions between different regimes ^[54]. We hope to formulate more realistic models considering different types of noise in future research.

Acknowledgement

The authors are grateful to the anonymous referees for their valuable comments and suggestions which led to the improvement of this paper.

Conflict of interest

The authors declare that there is no conflict of interest regarding the publication of this manuscript.