A Mathematical Model to Study the Potential Hepatitis B Virus Carriers and Effect of Vaccination Strategies in China

1. Introduction

Hepatitis B is a potentially life-threatening infection of the liver caused by the hepatitis B virus [1]. According to the World Health Organization data [2], about 2 billion people worldwide are infected with the hepatitis B virus, 250 million people carry the hepatitis B virus, and about 1 million people die each year from chronic active hepatitis, cirrhosis or primary liver cancer. The World Health Organization has listed HBV infection as one of the top 10 causes of death worldwide. China is one of the highest burdens of hepatitis B in the world, with annual data on new cases of hepatitis B from 2003 to 2021, is shown in Figure 1, with an average of 1,007,634 new cases per year, and the lowest number of new cases 719,011 in 2003 and the highest number of cases 1,169,946 in 2007; and the mortality of liver cancer caused by hepatitis B is the highest too, with an average of 614 deaths per year, with the lowest number of death cases 352 in 2015 and the highest number death cases 995 in 2006, as shown in Figure 1; the average annual mortality rate is 0.0452 (units of 1/100,000), with the lowest being 0.0258 in 2015 and the highest being 0.0761 in 2006, as shown in Figure 2.

Hepatitis B is a viral infection that damages the liver. The main modes of transmission of hepatitis B virus are sexual transmission, mother-to-child transmission, father-to-child transmission, medical transmission, and blood transmission. Hepatitis B virus does not transmit through the following channels: sharing utensils, breastfeeding, hugging, kissing, shaking hands, coughing, sneezing, or playing in public swimming pools or similar behaviors [2].

The incubation period of hepatitis B virus infection is usually 1-6 months, with an average of 3 months, and patients with incubation period are infectious. Vaccination will be good to prevent the infection of hepatitis B virus, but the vaccine may be ineffective after a period of time and there will be a risk of re-infection. After having body fluid with a patient, the person needs to be observed for 6 months and tested for fifth hepatitis B tests to determine whether he is infected [1]. Hepatitis B carriers are infectious as long as the viral load reaches a certain concentration.

From 2002 to 2007, the World Health Organization and the Chinese Ministry of Health jointly implemented the China-GAVI program of hepatitis B vaccine catch-up for adolescents in 12 provinces in western China, which aims to reduce the high incidence of hepatitis B in China and has achieved good results. The project covers all regions of 12 provinces, autonomous regions and municipalities directly under the Central Government, including Sichuan, Guizhou, Yunnan, Tibet, Chongqing, Shaanxi, Gansu, Qinghai, Ningxia, Xinjiang, Guangxi, Inner Mongolia, as well as the national poor counties in 10 provinces of Hunan, Hubei, Shanxi, Jiangxi, Anhui, Henan, Heilongjiang, Hainan, Hebei and Jilin. (10 non-western provinces region). The target population of the project is all newly born children in the above-mentioned areas during the project implementation period [3]. After 2008, the number of hepatitis B cases showed a significant decline, but the total number of cases continued to be high. After 2015, as the number of cases continued to increase. Now, the vaccination rate for hepatitis B in children is 100%. Normally adult hepatitis B carriers who have been diagnosed as infectious should consciously control their sexual behavior and not infect others, but the annual number of hepatitis B cases remains high and even continues to increase, probably from potentially infectious hepatitis B carriers. Therefore, it is important to correctly assess the number of potentially infectious hepatitis B virus carriers and to explore the implementation of booster vaccination strategies for susceptible individuals for the prevention and control of hepatitis B virus.

Many researchers in China and abroad have studied the process of hepatitis B virus transmission and the effect of prevention and control measures by building mathematical models. Su et al. developed a dynamic model of hepatitis B virus transmission with mechanisms of vaccination and vertical transmission [4]; Sun et al. analyzed the effect of physical examination and immunization on hepatitis B [5];Li et al. established a mathematical model for pharmacological blockade of HBV transmission between mothers and infants during maternal pregnancy [6];Liu et al. considered a SIRS infectious disease model with age structure and continuous vaccination [7]; Tao YJ et al considered a model of hepatitis B virus transmission with vertical transmission and period of infection when the incidence function is nonlinear [8]; Zhao T. et al. constructed a hepatitis B kinetic model and analyzed the data of new cases of hepatitis B in Xinjiang from 2005 to 2014[9]; O'Leary, C. et al. developed, including mathematical models of HBV drug treatment and immunization, to compare the effects of vaccination and treatment [10]; Jianhua Pang et al. developed a model to explore the impact of vaccination and other HBV infection control measures [11]; Habenom. H et al. developed a fractional order mathematical model of hepatitis B virus transmission and studied possible strategies for hepatitis B vaccination and isolation control [12]; Ramses Djidjou Demasse et al. developed a model of the age structure of HBV transmission that included both symptomatic and asymptomatic infections [13]. From 2002 to 2007, the China-GAVI program was implemented and chronic HBV infection in children decreased, but the number of new cases per year in China continues to increase, which we believe is mainly due to potentially infectious hepatitis B virus carriers.

In this study, a dynamic model based on the pathogenesis of hepatitis B is developed to assess the number of potential hepatitis B virus carriers and to explore effective strategies to control hepatitis B virus transmission through vaccination.

The paper includes four parts. The transmission of hepatitis B virus in China is introduced in the first subsection; in the second subsection, a model is developed and the dynamic analysis is performed; in the third subsection, a simulation analysis is performed to explore the effect of the parameters on potential hepatitis B virus carriers K, the percentage of K, and$R_c$; the fourth part is a simple discussion.

2. Methods

2.1. Model building

The main objectives of this paper are to accurately estimate the number of potential hepatitis B virus carriers and we develop a dynamic model that include potential hepatitis B virus carriers based on the transmission mechanism of hepatitis B virus. Because acute infection individuals have an extremely low mortality rate with more than 90% being cured [14], and the mortality due to disease of acute infection individuals is not considered. About 88% of acute infection individuals will become chronic HBV carriers, and the recovery rate of chronic HBV carriers is very low and they can easily become hepatitis B patients again after recovery, therefore recovery of chronic HBV carriers is not considered. Acute infection individuals and chronic HBV carriers consciously control their sexual behavior because they know they are carriers of the virus, so assuming that it is no longer infectious, hepatitis B is contagious during the incubation period, and potential hepatitis B virus carriers are also contagious. Divided the population into 7 compartments: susceptible population S, exposed E, acute infection individuals A, chronic HBV carriers C, potential hepatitis B virus carriers K, recovered R, and vaccinated V. Based on the process of hepatitis B virus transmission the following flow chart is shown in Figure 3.

The meaning of each parameter symbol: $\mathsf{\Lambda }$ is the number of new births rate every year. $\alpha$ is the vaccination rate. $\omega$ is the vaccine failure rate. $\beta$ is the base transmission rate. $\epsilon$ is the transfer rate of A to C. $\rho$ is the transfer rate from A to R. $\mu$ is the natural mortality rate and ${\lambda }_1$ is the conversion rate of E to A, ${\lambda }_2$ is the conversion rate of E to C, ${\lambda }_3$ is the conversion rate of E to K. ${\gamma }_1$ is the mortality rate of chronic HBV carriers due to disease$\mathrm{a}\mathrm{n}\mathrm{d}{\gamma }_2$ is the mortality rate of patients in K due to disease. ${\eta }_1$ is the transfer rate from E to A, ${\eta }_2$ is the transfer rate from E to C and ${\eta }_3$ is the transfer rate from E to K respectively. q is the proportion of people with acute infection who become chronically infected. All the parameters are nonnegative.

Based on Figure 3, the following dynamic model of hepatitis B virus transmission is developed:

$$\left\{\begin{array}{c}\frac{dS}{dt}=\Lambda +\omega V-\alpha S-\beta S(E+K)-\mu S,\hfill \\ \frac{dE}{dt}=\beta S\left(E+K\right)-{\lambda }_1{\eta }_{1}E-{\lambda }_2{\eta }_{2}E-{\lambda }_3{\eta }_{3}E-\mu E,\hfill \\ \frac{dA}{dt}={\lambda }_1{\eta }_{1}E-q\epsilon A-\left(1-q\right)\rho A-\mu A,\hfill \\ \frac{dC}{dt}={\lambda }_2{\eta }_{2}E+q\epsilon A-{\gamma }_{1}C-\mu C,\hfill \\ \frac{dK}{dt}={\lambda }_3{\eta }_{3}E-{\gamma }_{2}K-\mu K, \hfill \\ \frac{dR}{dt}=(1-q)\rho A-\mu R, \hfill \\ \frac{dV}{dt}=\alpha S-\omega V-\mu V. \hfill \end{array}\right.$$

(1)

2.2. Control reproduction number

The control reproduction number refers to the number of secondary cases after an infected individual enters the susceptible population [15]. The next generation matrix method is used to calculate it.

The disease state of system (1) are E, A, C, K, to calculate the number of control reproduction of system (1), let

$$\mathcal{F}=\left[\begin{array}{c}\beta S(E+K)\\ 0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\right]，\mathcal{V}=\left[\begin{array}{c}{\lambda }_1{\eta }_{1}E+{\lambda }_2{\eta }_{2}E+{\lambda }_3{\eta }_{3}E+\mu E\\ q\epsilon A+\left(1-q\right)\rho A+\mu A-{\lambda }_1{\eta }_{1}E\\ \begin{array}{c}{\gamma }_{1}C+\mu C-{\lambda }_2{\eta }_{2}E-q\epsilon A \\ {\gamma }_{2}K+\mu K-{\lambda }_3{\eta }_{3}E\end{array}\end{array}\right].$$

The Jacobi matrices of $\mathcal{F}$、$\mathcal{V}$ at the disease-free equilibrium point $P_0$ are:

$$F=\left[\begin{array}{c} \beta S_0 \begin{array}{ccc}0& 0& \beta S_0\end{array}\\ 0\begin{array}{ccc} 0& 0& 0\end{array}\\ \begin{array}{c}0\begin{array}{ccc} 0& 0& 0\end{array}\\ 0\begin{array}{ccc} 0& 0& 0\end{array}\end{array}\end{array}\right],$$

$$V=\left[\begin{array}{ccc}{\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu & 0& \begin{array}{cc}0 & 0\end{array}\\ -{\lambda }_1{\eta }_1& q\epsilon +\left(1-q\right)\rho +\mu & \begin{array}{cc}0 & 0\end{array}\\ \begin{array}{c}-{\lambda }_2{\eta }_2\\ -{\lambda }_3{\eta }_3\end{array}& \begin{array}{c}-q\epsilon \\ 0\end{array}& \begin{array}{c}\begin{array}{cc}{\gamma }_1+\mu & 0 \end{array}\\ \begin{array}{cc} 0 & { \gamma }_2+\mu \end{array}\end{array}\end{array}\right].$$

At the disease-free equilibrium point, we can derive

$$F{V}^{-1}=\left[\begin{array}{cccc}\frac{\beta S_0}{{\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu }+\frac{\beta S_0{\lambda }_3{\eta }_3}{\left({\gamma }_2+\mu \right)\left({\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right)}& 0& 0& -\frac{\beta S_0}{{ \gamma }_2+\mu }\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$$

The value of the control reproduction number $R_c$ is the spectral radius of $F{V}^{-1}$, from which we can calculate

$$R_c=\rho \left(F{V}^{-1}\right)=\frac{\beta S_0\left({\gamma }_2+\mu +{\lambda }_3{\eta }_3\right)}{\left({\gamma }_2+\mu \right)\left({\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right)}.$$

2.3. Disease-free equilibrium point and its stability

Clearly, there is a disease-free equilibrium point $P_0=\left(S_0{,\mathrm{0,0},\mathrm{0,0},0,V}_0\right)$ in system (1). Where

$$S_0=\frac{\left(\omega +\mu \right)\mathsf{\Lambda }}{\mu \left(\alpha +\mu +\omega \right)}，V_0=\frac{\mathsf{\alpha }\mathsf{\Lambda }}{\alpha \mu +{\mu }^2+\mu \omega }.$$

Theorem 1. 

The disease-free equilibrium point$P_0$of the system (1) is locally asymptotically stable when $R_c<1$; it is unstable when$R_c>1$.

Proof: 

the Jacobi matrix of the system (1) at the disease-free equilibrium point $P_0$ is

$$\left[\begin{array}{ccc}-\alpha -\mu & -\beta S_0& 0\\ 0& \beta S_0-{\lambda }_1{\eta }_1-{\lambda }_2{\eta }_2-{\lambda }_3{\eta }_3-\mu & 0\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \alpha \end{array}\end{array}& \begin{array}{c}{\lambda }_1{\eta }_1\\ {\lambda }_2{\eta }_2\\ \begin{array}{c}\begin{array}{c}{\lambda }_3{\eta }_3\\ 0\end{array}\\ 0\end{array}\end{array}& \begin{array}{c}-q\epsilon -\left(1-q\right)\rho -\mu \\ q\epsilon \\ \begin{array}{c}0\\ \begin{array}{c}\left(1-q\right)\rho \\ 0\end{array}\end{array}\end{array}\end{array} \begin{array}{ccc}0& \begin{array}{cc}-\beta S_0& 0\end{array}& \omega \\ 0& \begin{array}{cc}\beta S_0 & 0\end{array}& 0\\ \begin{array}{c}0\\ -{\gamma }_1-\mu \\ \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ 0\end{array}\end{array}& \begin{array}{c} \begin{array}{cc} 0 & 0\end{array}\\ \begin{array}{cc}0& 0\end{array}\\ \begin{array}{c}\begin{array}{cc}-{\gamma }_2-\mu & 0\end{array}\\ \begin{array}{c}\begin{array}{cc}0& -\mu \end{array}\\ \begin{array}{cc}0& 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ -\omega -\mu \end{array}\end{array}\end{array}\right].$$

The characteristic polynomial of this matrix is:

Clearly

$${\varphi }_1=-\mu <0,$$

$${\varphi }_2=-{\gamma }_1-\mu <0,$$

$${\varphi }_3=-\mathrm{q}\mathsf{\epsilon }-\left(1-\mathrm{q}\right)\mathsf{\rho }-\mathsf{\mu }<0.$$

The other two eigenvalues ${\varphi }_4$、${\varphi }_5$ satisfy the equation

$$f_1\left(\varphi \right)={\varphi }^2+{\mathrm{A}}_1\varphi +A_2,$$

where

$$A_1=\omega +2\mu +\alpha , A_2=\mu \left(\alpha +\omega +\mu \right).$$

We obtain

$${\Delta }_1=m_0=A_1>0$$

$${∆}_2=\left|\begin{array}{cc}{m}_0& 1\\ 0& m_1\end{array}\right|=m_0{m}_1=A_1{·A}_2>0.$$

This means that ${\varphi }_4$、${\varphi }_5$ have negative real parts. The other two eigenvalues ${\varphi }_6$、${\varphi }_7$ are the zeros of the function

$$f_2\left(\varphi \right)={\varphi }^2+A_3\varphi +A_4,$$

Here

$${ A}_3=-\beta S_0+{\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu {+\gamma }_2+\mu ,$$

$$A_4=-\beta S_0{\gamma }_2+{\lambda }_1{\eta }_1{\gamma }_2+{\lambda }_2{\eta }_2{\gamma }_2+{\lambda }_3{\eta }_3{\gamma }_2+\mu {\gamma }_2-\beta S_0\mu +{\lambda }_1{\eta }_1\mu +{\lambda }_2{\eta }_2\mu +{\lambda }_3{\eta }_3 \mu +{\mu }^2-\beta S_0{\lambda }_3{\eta }_3.$$

When $R_c<1$, $A_3>0$，$A_4>0$, according to the Routh-Hurwitz discriminant, the characteristic roots ${\varphi }_6$、${\varphi }_7$ have negative real parts.

When $R_c>1$, $A_4<0$ , $f\left(\varphi \right)$ has at least one characteristic root with positive real part, the system (1) is unstable at the disease-free equilibrium point $P_0$.

In summary, when $R_c<1$, all eigenvalues of the Jacobi matrix of system (1) have negative real parts and system (1) is locally asymptotically stable at the disease-free equilibrium point $P_0$;when $R_c>1$,the system (1) is unstable at the disease-free equilibrium point $P_0$.

2.4. Endemic equilibrium point and its stability

When $A^{*}\ne 0$, the system (1) exists the endemic equilibrium point

$$P^{*}=\left(S^{*}{,E}^{*}{,A}^{*},{C,}^{*}{K}^{*},V^{*},R^{*}\right)$$

Here

$$S^{*}=\frac{\Lambda \left(\omega +\mu \right){\lambda }_1{\eta }_1\left({\gamma }_2+\mu \right)}{\mu {\lambda }_1{\eta }_1\left({\gamma }_2+\mu \right)\left(\alpha -\omega -\mu \right)-\beta \left(\omega +\mu \right)\left(\epsilon q{A}^{*}+\left(1-q\right)\rho A^{*}+\mu A^{*}\right)\left({\gamma }_2+\mu +{\lambda }_3{\eta }_3\right)},$$

$$E^{*}=\frac{\epsilon q{A}^{*}+\left(1-q\right)\rho A^{*}+\mu A^{*}}{{\lambda }_1{\eta }_1},$$

$$C^{*}=\frac{{\lambda }_2{\eta }_2\left(\epsilon q{A}^{*}+\left(1-q\right)\rho A^{*}+\mu A^{*}\right)+{\lambda }_1{\eta }_{1}q\epsilon A^{*}}{{\lambda }_1{\eta }_1\left({\gamma }_1+\mu \right)},$$

$$K^{*}=\frac{{\lambda }_3{\eta }_3\left(\epsilon q{A}^{*}+\left(1-q\right)\rho A^{*}+\mu A^{*}\right)}{{\lambda }_1{\eta }_1\left({\gamma }_2+\mu \right)},$$

$$V^{*}=\frac{\Lambda \alpha {{\lambda }_1\eta }_1\left({\gamma }_2+\mu \right)}{\mu {\lambda }_1{\eta }_1\left({\gamma }_2+\mu \right)\left(\alpha -\omega -\mu \right)-\beta \left(\omega +\mu \right)\left(\epsilon q{A}^{*}+\left(1-q\right)\rho A^{*}+\mu A^{*}\right)\left({\gamma }_2+\mu +{\lambda }_3{\eta }_3\right)},$$

$$R^{*}=\frac{(1-q)\rho A^{*}}{\mu }.$$

Theorem 2. 

When $R_c>1$ the system (1) has a unique endemic equilibrium point.

Proof. 

Using the equilibrium equation of the system (1), we can obtain

$$K^{*}=\frac{{\lambda }_3{\eta }_3{E}^{*}}{{\gamma }_2+\mu },$$

$$A^{*}=\frac{{\lambda }_1{\eta }_1{E}^{*}}{q\epsilon +\left(1-q\right)\rho +\mu },$$

$$C^{*}=\frac{{\lambda }_2{\eta }_2{E}^{*}}{{\gamma }_1+\mu }+\frac{{q\epsilon \lambda }_1{\eta }_1{E}^{*}}{\left({\gamma }_1+\mu \right)\left(q\epsilon +\left(1-q\right)\rho +\mu \right)},$$

$$R^{*}=\frac{\left(1-q\right)\rho {\lambda }_1{\eta }_1{E}^{*}}{\left(q\epsilon +\left(1-q\right)\rho +\mu \right)\mu },$$

$$S^{*}=\frac{\Lambda }{\beta \left(1+\frac{{\lambda }_3{\eta }_3}{{\gamma }_2+\mu }\right)E^{*}+\frac{\alpha \mu }{\omega +\mu }+\mu },$$

$$V^{*}=\frac{\alpha }{\omega +\mu }\bullet \frac{\Lambda }{\beta \left(1+\frac{{\lambda }_3{\eta }_3}{{\gamma }_2+\mu }\right)E^{*}+\frac{\alpha \mu }{\omega +\mu }+\mu },$$

Further, according to the equilibrium equation we get

$$\Lambda -\mu S^{*}-{\lambda }_1{\eta }_1{E}^{*}-{\lambda }_2{\eta }_2{E}^{*}-{\lambda }_3{\eta }_3{E}^{*}-\mu E^{*}+\mu V^{*}=0.$$

From this, we can get

$$E^{*}=\frac{\mathsf{\Lambda }\beta \left(1+\frac{{\lambda }_3{\eta }_3}{{\gamma }_2+\mu }\right)-\left(\frac{\alpha \mu }{\omega +\mu }+\mu \right)\left({\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right)}{-\beta \left(1+\frac{{\lambda }_3{\eta }_3}{{\gamma }_2+\mu }\right)\left({\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right)}.$$

(2)

When $R_c>1,$

$$\beta S_0\left({\gamma }_2+\mu +{\lambda }_3{\eta }_3\right)>\left({\gamma }_2+\mu \right)\left({\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right).$$

(3)

Put

$$S_0=\frac{\Lambda }{\mu \left(\frac{\alpha }{\omega +\mu }+1\right)}.$$

into equation (3) and simplify we get

$$\beta \Lambda \left(1+\frac{{\lambda }_3{\eta }_3}{{\gamma }_2+\mu }\right)>\mu \left({\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right)\left(\frac{\alpha }{\omega +\mu }+1\right)$$

$$1-R_c=\frac{\mu \left(\frac{\alpha }{\omega +\mu }+1\right)\left({\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right)-\beta \Lambda \left(1+\frac{{\lambda }_3{\eta }_3}{{\gamma }_2+\mu }\right)}{\mu \left(\frac{\alpha }{\omega +\mu }+1\right)\left({\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right)}.$$

Here

$$E^{*}=\frac{\mu \left({\gamma }_2+\mu \right)\left(\alpha +\omega +\mu \right)\left(R_c-1\right)}{\beta \left({\gamma }_2+\mu +{\lambda }_3{\eta }_3\right)\left(\omega +\mu \right)}.$$

when $R_c>1$，$E^{\mathrm{*}}>0.$

Therefore, when $R_c>1$, the system (1) has a unique endemic equilibrium point $P^{*}.$

Theorem 3. 

When $R_c>1,$ the system (1) is locally asymptotically stable at the endemic equilibrium point$P^{*}$.

Proof. 

Let the Jacobi matrix of the system (1) at the endemic equilibrium point $P^{*}$

be J.

$$\left[\begin{array}{ccc}-\alpha -\mu -\beta \left(E^{*}+K^{*}\right)& -\beta S^{*}& 0\\ \beta \left(E^{*}+K^{*}\right)& \beta S^{*}-{\lambda }_1{\eta }_1-{\lambda }_2{\eta }_2-{\lambda }_3{\eta }_3-\mu & 0\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \alpha \end{array}\end{array}& \begin{array}{c}{\lambda }_1{\eta }_1\\ {\lambda }_2{\eta }_2\\ \begin{array}{c}\begin{array}{c}{\lambda }_3{\eta }_3\\ 0\end{array}\\ 0\end{array}\end{array}& \begin{array}{c}-q\epsilon -\left(1-q\right)\rho -\mu \\ q\epsilon \\ \begin{array}{c}0\\ \begin{array}{c}\left(1-q\right)\rho \\ 0\end{array}\end{array}\end{array}\end{array} \begin{array}{ccc}0& \begin{array}{cc}-\beta S^{*}& 0\end{array}& \omega \\ 0& \begin{array}{cc}\beta S^{*} & 0\end{array}& 0\\ \begin{array}{c}0\\ -{\gamma }_1-\mu \\ \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ 0\end{array}\end{array}& \begin{array}{c} \begin{array}{cc} 0 & 0\end{array}\\ \begin{array}{cc}0& 0\end{array}\\ \begin{array}{c}\begin{array}{cc}-{\gamma }_2-\mu & 0\end{array}\\ \begin{array}{c}\begin{array}{cc}0& -\mu \end{array}\\ \begin{array}{cc}0& 0\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}0\\ 0\\ \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ -\omega -\mu \end{array}\end{array}\end{array}\right].$$

The characteristic polynomial is obtained by calculating

$$f\left(\varphi \right)=g_1\left(\varphi \right)·g_2\left(\varphi \right),$$

where

$$g_1\left(\varphi \right)=\left(\varphi +\mu \right)\left(\varphi +{\gamma }_1+\mu \right)\left(\varphi +q\epsilon +\left(1-q\right)\rho +\mu \right),$$

$$\begin{gathered} g_2\left(\varphi \right)=\left(\varphi +\omega +\mu \right)\left\{\left(\varphi +{\gamma }_2+\mu \right)\left[\left(\varphi +\alpha +\mu \right)\left(\varphi -\beta S^{*}+{\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right)+ \\ \beta S^{*}\left(E^{*}+K^{*}\right)\left(\varphi +{\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right)\right]-\beta S^{*}{\lambda }_3{\eta }_3\left(\varphi +\alpha +\mu \right)\right\} \end{gathered}$$

$$-\alpha \omega \left[\left(\varphi -\beta S^{*}+{\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right)\left(\varphi +{\gamma }_2+\mu \right)-\beta S^{*}{\lambda }_3{\eta }_3\right].$$

It is easy to know that all the eigenvalues of $g_1\left(\varphi \right)$ have negative real parts, and the other four eigenvalues are the zeros of the function $g_2\left(\varphi \right)$.

From $g_2\left(\varphi \right)=0$ we get

$$\begin{gathered} \beta S^{*}{\lambda }_3{\eta }_3\left(\varphi +\mu \right)\left(\varphi +\alpha +\omega +\mu \right) \\ =\left(\varphi +\mu \right)\left(\varphi +\alpha +\omega +\mu \right)\left(\varphi +{\gamma }_2+\mu \right)\left(\varphi -\beta S^{*}+{\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3 \\ +\mu \right)+\beta S^{*}\left(E^{*}+K^{*}\right)\left(\varphi +\omega +\mu \right)\left(\varphi +{\gamma }_2+\mu \right)\left(\varphi +{\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right). \end{gathered}$$

Here

$$\begin{gathered} 1=\frac{\left(\varphi +{\gamma }_2+\mu \right)\left(\varphi -\beta S^{*}+{\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right)}{\beta S^{*}{\lambda }_3{\eta }_3} \\ +\frac{\beta S^{*}\left(E^{*}+K^{*}\right)\left(\varphi +\omega +\mu \right)\left(\varphi +{\gamma }_2+\mu \right)\left(\varphi +{\lambda }_1{\eta }_1+{\lambda }_2{\eta }_2+{\lambda }_3{\eta }_3+\mu \right)}{\beta S^{*}{\lambda }_3{\eta }_3\left(\varphi +\mu \right)\left(\varphi +\alpha +\omega +\mu \right)}. \end{gathered}$$

(4)

Assume that $Re\left(\varphi \right)\ge 0$ when $g_2\left(\varphi \right)=0$. Taking the norm on the right side of the equal sign of equation (4) can get

(5)

We can find that the results of equation (4) and equation (5) contradict each other, so the roots of $g_2\left(\varphi \right)$ also have negative real parts.

When $R_c>1$ , all the characteristic roots of the Jacobi matrix of system (1) have negative real parts, so the system (1) is locally asymptotically stable at the endemic equilibrium point $P^{*}$.

3. Results

The average incubation period of hepatitis B virus infection is 105 days, then we get ${\eta }_1={\eta }_2={\eta }_3=3.4$. Because the GAVI program vaccinates all newborns[16], the vaccination rate $\alpha =$ 1. The parameters ${\lambda }_1$ 、${\lambda }_2$ are obtained by data fitting, and ${\lambda }_1+{\lambda }_2+{\lambda }_3=1$. ${\lambda }_3$ is obtained by calculation. ${\gamma }_1$、${\gamma }_2$ are obtained from the mortality rate of hepatitis B in China published by the Chinese Center for Disease Control and Prevention [17].

The initial value of S is taken from the total national population in 2003 published by the National Bureau of Statistics of China [18], the initial value of C is the total number of hepatitis B cases nationwide in 2003 published by the National Health Commission of China [19], The initial value of A is taken by the ratio of the number of acute hepatitis B incidence to the number of chronic hepatitis B incidence of 1:5. The initial value of V is taken from the total number of people vaccinated against hepatitis B nationwide in 2002 published [18], and the initial value of R is calculated by known parameters. So here S(0) =1292270000, E(0)=250000, A(0)=143802, C(0)=719011, K(0)=180000, V(0)=6540000, R(0)=68737.

We obtained the values of parameters β、${\lambda }_1$、${\lambda }_2$ by fitting nonlinear least squares to the data of new cases per year, and the results are shown in Table 1. It can be seen that the percentage of potentially infectious in the compartment is as high as 30.49%, the fitting results are shown in Figure 4. The actual cumulative number of cases and the curve of the fitted cumulative number of cases are shown in Figure 5. The trend of the fitted curve is consistent with the actual curve.

3.1. The number of A, C and potential virus carriers K with time

We substitute the fitted parameters into the control regeneration number formula and calculate the control reproduction number R_c = 1.741. The fitted parameters are used to calculate the number of acute infection individuals A, chronic HBV carriers C and potential hepatitis B virus carriers K, and shown in Figure 6. The number of people with chronic hepatitis B infection C peaked in 2006, while the number of people with chronic hepatitis B infection was 1,181,405, after which the number of people with chronic hepatitis B infection gradually decreased. The number of acute hepatitis B infections first dropped to 2471 and gradually stabilized at 2188 after 2006. Based on the fitted data, the average number of potential hepatitis B carriers is estimated to be 449535 (95% CI [415651, 483420]).The number of potentially infectious hepatitis B virus carriers K is in a continuous rise from 2003, reaching 466,907 in 2009 and stabilizing at 466,907 after 2009. This suggests that the GAVI program, which was implemented from 2002 to 2007 and only provided booster immunization to children under 15 years of age, did not have a major impact on the number of potentially infectious infections among adults. However, the number of new cases of acute infections in the latter 8 years had a good control effect, significantly reducing the number of cases of infection in adolescents entering adulthood.

3.2. The impact of parameter ${\mathit{\lambda }}_3$

on hepatitis B virus transmission

The conversion rate λ₃ from latently infected E to potential hepatitis B virus carrier K has an impact on the number of hepatitis B patients, when the conversion rate λ₂ from latently infected E to chronic hepatitis B infected C increases, the conversion rate λ₁ from latently infected E to acute hepatitis B patient A remains constant, and the conversion rate λ₃ from latently infected E to potential hepatitis B virus carrier K gradually decreases, as shown in Figure 7. This indicates that the number of new cases of hepatitis B per year decreases as λ₃ decreases, and if λ₃ decreases by 0.2, the number of new cases of hepatitis B per year decreases by an average of 257882. Therefore, increasing the level of hepatitis B in the testing population will reduce the percentage of potential hepatitis B virus carriers.

on hepatitis B virus transmission.

3.3. Effect of vaccine failure rate $\mathit{\omega }$

on potential hepatitis B virus carriers

We obtain that the vaccine failure rate ω has an effect on the change of the number of potential hepatitis B virus carriers K according to numerical simulation. As the vaccine failure rate ω decreases, the number of potential hepatitis B virus carriers K gradually decreases, but K decreases more and more slowly, see Figure 8, the lower the vaccine failure rate, the smaller the amount of potential hepatitis B virus carriers will be, reducing the vaccine failure rate ω by 1%, the average number of potential hepatitis B virus carriers K decreases by 17,454. Therefore, the number of potential hepatitis B virus carriers K can be reduced by decreasing the vaccine failure rate.

on potential virus carriers K.

The vaccine failure rate ω affects the proportion of potential virus carriers K, see Figure 9. If the vaccine failure rate ω is reduced by 2%, the number of potential hepatitis B virus carriers K decreases in proportion by 1.41% on average. The trend of all curves is consistent, and the time of occurrence of the lowest and highest points is also consistent.

When the vaccine failure rate changed, the ratio of the number of potential hepatitis B virus carriers K to the total number of cases (A+C+K) also changed, as shown in Table 2. The smaller the vaccine failure rate ω, the smaller the proportion of K and the smaller the number of potential hepatitis B virus carriers K. Therefore, reducing the vaccine failure rate can effectively reduce the number of potential hepatitis B virus carriers K.

on the percentage of potential virus carriers K.

3.4. Effect of vaccine failure rates $\mathit{\omega }$

on the number of C when vaccination rate is 100%

In China, vaccination rates for children under five years of age have reached more than 95% and almost 100% since 2014. When the vaccination rate is 100% ($\alpha$=1), the effect of different vaccine failure rates $\omega$ on the number of chronic HBV carriers , the rate of change in the number of chronic HBV carriers is determined by the vaccine failure rate $\omega$, and the number of chronic HBV carriers decreases as the vaccine failure rate $\omega$ decreases, shown in Figure 10, and we obtain that decreasing the vaccine failure rate $\omega$ by 2% can reduce the number of chronic HBV carriers by an average of 51993, which indicates that the lower the vaccine failure rate, the lower the number of chronic HBV carriers.

on the number of chronic HBV carriers.

3.5. Effect of vaccination rates $\mathit{\alpha }$

on the number of C when vaccine failure rate is 0.1

When the vaccine failure rate $\omega$ is 0.1, different vaccination rates $\alpha$ will have an impact on the number of chronic HBV carriers, shown in Figure 11. When the vaccine failure rate $\omega$=0.1, from which it can be seen that the vaccination rate $\alpha$ determines the rate of change in the number of chronic HBV carriers; the number of chronic HBV carriers decreases as the vaccination rate $\alpha$ increases, increasing $\alpha$ by 20%, the number of chronic HBV carriers decreases by 169,469 on average, indicating that the higher the vaccination rate, the lower the number of chronic hepatitis B cases. It shows that the peak number of chronic HBV carriers is higher when $\alpha$ is smaller, and the peak occurs later than when $\alpha$=1.

on the number of chronic HBV carriers C.

3.6. Effect on $R_c$

$R_c$ is the threshold value to determine whether the disease is extinct or not.When $R_c<1$，the disease gradually dies out; when $R_c>1$，the disease develops into an endemic disease[4].

The vaccine failure rate $\omega$ and the basal transmission rate $\beta$, and the vaccination rate $\alpha$ have an effect on the control reproduction number $R_c$, shown in Figure 12 and Figure 13. It shows the control reproduction number $R_c$ increases with the increase of vaccine failure rate $\omega$ and basal transmission rate $\beta$; when the vaccination rate $\alpha$ increases and the vaccine failure rate $\omega$ decreases, the control reproduction number $R_c$ decreases, indicating that the higher the vaccination rate, the slower the transmission of hepatitis B. To better control the spread of hepatitis B virus, the best prevention and control strategy is to strengthen hepatitis B vaccination, not only for newborns, but also for people of different ages with booster vaccination.

and the basal transmission rate $\beta$ on $R_c$.

and the vaccination rate $\alpha$ on $R_c$.

4. Conclusion

China is one of the countries with the highest burden of hepatitis B in the world, increasing incidence of hepatitis B in China, and a reasonable estimate of the number of potential HBV carriers is crucial for the prevention and control of hepatitis B. Therefore, in this paper, a model of hepatitis B virus transmission dynamics with vaccination is developed based on the transmission mechanism of hepatitis B virus to study the number of potential hepatitis B virus carriers. The number of potential hepatitis B virus carriers K is obtained as 449535(95%CI[415651, 483420]) by fitting the parameters in the model using nonlinear least squares method, K is calculated to account for 30.49% of those with hepatitis B. The control reproduction number $R_c$=1.741 for hepatitis B virus transmission in China is obtained according to the fitted parameters and the next generation matrix method, this indicates that the transmission of hepatitis B virus in China is still in a critical situation, the number of new hepatitis B cases continues to increase at a high level in the population, and additional control is needed to slow the spread of the virus.

In this study, the effect of different parameters on hepatitis B virus transmission is considered, and it is shown that the conversion rate ${\lambda }_3$ from exposed E to potential hepatitis B virus carrier K affects the transmission of hepatitis B virus. When the vaccination rate $\alpha$ increases and the vaccine failure rate $\omega$ decreases, the control reproduction number $R_c$ decreases. The analysis reveals that reducing the vaccine failure rate, increasing the vaccination rate, and increasing hepatitis B testing for the public are essential to reduce the number of hepatitis B cases and control the number of potential hepatitis B virus carriers K. Therefore, encouraging hepatitis B vaccination for all types of people and increasing the rate of hepatitis B vaccination can effectively control the spread of hepatitis B virus.

This paper estimates that the number of potentially infectious hepatitis B virus carriers is high, maybe this is the reason the number of hepatitis B patients in China to continue to rise. To better control the spread of hepatitis B virus, an optional prevention and control strategy is to increase vaccination of different age groups; at the same time, it is necessary to increase publicity on prevention of hepatitis B virus transmission, to help the public to correctly understand the transmission of hepatitis B and do a good job of protection, and to call on the public to strengthen their own exercise. All have a positive impact on the prevention and treatment of hepatitis B. Our study is based only on the number of annual incidences of hepatitis B. Because the official website of the Chinese Health and Welfare Bureau only has data on the annual incidence of hepatitis B, there are no data on drug users and gender, although these two data have an impact on the assessment of epidemic transmission. Annual data on new cases illustrate that there are significant differences in the number of hepatitis B cases by age group, and we will discuss the impact of age heterogeneity on hepatitis B virus transmission in a follow-up study to find a more optimal vaccination strategy.