Modelling homosexual and heterosexual transmissions of hepatitis B virus in China

Abstract

Studies have shown that sexual transmission, both heterosexually and homosexually, is one of the main ways of HBV infection. Based on this fact, we propose a mathematical model to study the sexual transmission of HBV among adults by classifying adults into men and women and considering both same-sex and opposite-sex transmissions of HBV in adults. Firstly, we calculate the basic reproduction number $R_0$ and the disease-free equilibrium point $E_0$. Secondly, by analysing the sensitivity of $R_0$ in terms of model parameters, we find that the infection rate among people who have same-sex partners, the frequency of homosexual contact and the immunity rate of adults play important roles in the transmission of HBV. Moreover, we use our model to fit the reported data in China and forecast the trend of hepatitis B. Our results demonstrate that popularizing the basic knowledge of HBV among residents, advocating healthy and reasonable sexual life style, reducing the number of adult carriers, and increasing the immunization rate of adults are effective measures to prevent and control hepatitis B.

Keywords:

• Hepatitis B virus
• heterosexual transmission
• homosexual transmission
• mathematical modelling
• basic reproduction number
• sensitivity analysis

1. Introduction

Hepatitis B is a liver disease that results from infection with the hepatitis B virus (HBV). It ranges from a mild illness lasting a few weeks to a serious, lifelong illness. HBV spreads by blood, semen, or other body fluids. Among adults, HBV transmission occurs primarily among unvaccinated adults with risk behaviours, including having multiple sex partners and sex partners of people with chronic infection. In fact, HBV is easily transmitted through sexual activity and is believed 50–100 times more infectious that HIV (see [14,25]). It is also more infectious and more stable in the environment than hepatitis C virus. By 2015, 2 billion people worldwide had been infected with HBV, with an average annual death toll of 887,000. Till now, there is no effective treatment available for chronic HBV carriers, and immunization with hepatitis B vaccine remains the most important prevention measure.

In China, hepatitis B is one of the top three infectious diseases and about 300,000 people die of HBV-related diseases annually [5]. It has been estimated that 120 million people carry HBsAg based on the national hepatitis seroepidemiological survey in 1992 [16]. Hepatitis B infection was the one of the two greatest attributable proportion of cancer deaths by risk factor in China [13]. In 2002, the China-GAVI Project (CGP) was initiated to provide hepatitis B vaccine to infants born in the less developed areas of China including the western provinces and poverty counties of middle provinces. During the CGP between 2003 and 2009, an estimated 3.8 million chronic HBV infections and 680,000 deaths were prevented in CGP areas [11]. In 2009, it was suggested to retroactively immunize young people (younger than 15 years old) in China. Now the incidence in children is much better controlled, but the incidence in adults keeps a high level. The reported HBV infections in adults aged $>20$ years increased by $22\mathrm{\%}$, from 740,000 in 2004 to 903,000 in 2014 [33,39].

The risk of developing to a carrier is dependent on the age at infection, and the transmission has different routes for adults and children [36,37]. It is known that sexual transmission is an important route of HBV spread in adults [39]. In previous work, heterosexual transmission was considered as the major route of transmission of hepatitis B among adults [1,22]. Mathematical models have focused on the heterosexual transmission of HBV [39]. There is evidence showing that homosexual contact is responsible for about $10\mathrm{\%}$ of all cases of hepatitis B [9] and HBV is very efficiently transmitted sexually during male homosexual contact [15,21]. WHO recommends all men who have sex with men (MSM) receive HBV testing since MSM are a high-risk group for HBV infection. There are some studies on transmission of HBV in homosexual men in United States and Europe [4,28]. In China, the health problems, especially the HIV infections of men who have sex with men have also attracted some attentions [23,24,31]. The prevalence of HBV also keeps a high level in China. For example, in Kunming, a city of Southwestern China, the prevalence of HBV were $7.7\mathrm{\%}$ among MSM [32].

Over the past 10 years, the health problems and needs of lesbians have attracted attentions, and more research has been conducted exploring sexual behaviour and health among the larger category of women who have sex with women (WSW) [26]. Although the sex distribution of HBsAg carriers indicate that the number of males is more than that of females [39], we believe that women-to-women transmission may also occur for HBV. Some high risk behaviours like sex with men, multiple bisexual partners, and contract marriages (marriages between gay men and lesbian women) exist among WSW [2,3,26]. A Chinese sociological study of WSW reported that oral-genital contact had a prevalence of $60-80\mathrm{\%}$ [26], which may cause HBV transmission.

Mathematical models have been used to describe various transmission routes and to design control strategies for different areas. Vaccination has been emphasized mostly in modelling, such as [19,21,27,36,37]. The effect of mother-to-child transmission of HBV has been considered in the models of [9,27,37–39]. The transmission of HBV and the development into chronic infection change with age have been described as either multi-age-group models [35,39] or age-structured models [27,36,38]. Sexual transmission of HBV between adults has also been discussed for different countries in [9,27,39].

In this paper, we construct a mathematical model to analyse the sexual transmission, both heterosexually and homosexually, of hepatitis B in China. Different from previous studies of HBV models for China [34,36–39], we focus on the sexual transmission of hepatitis B which only exists within adults. Therefore, we only consider the adult population in China and ignore the new born babies and children in the model. We analyse the disease-free equilibrium $E_0$, compute the basic reproduction number $R_0$ of the model, simulate the HBV data on infected male and female adults reported by Chinese Center for Disease Control and Prevention, and perform some sensitivity analysis of the basic reproduction number in terms of model parameters. Finally we provide some corrections to the analysis in reference [39] by some of us.

2. Sexual transmission model of hepatitis B

2.1. Model formulation

We only consider the adult population and divide the total population into susceptible individuals, infective individuals, carriers, and recovered individuals. The susceptible individuals are further divided into female susceptible individuals $S_f(t)$ and male susceptible individuals $S_m(t)$. Infective individuals are further divided into female infective individuals $I_f(t)$ and male infective individuals $I_m(t)$. Carriers are further divided into female carriers $C_f(t)$ and male carriers $C_m(t),$ and we uniformly denote the recover individuals as $R(t)$, thus we have the total population $N(t)=S_f(t)+I_f(t)+C_f(t)+S_m(t)+I_m(t)+C_m(t)+R(t).$ We ignore the loss of immunity after immunization, the immigration and emigration of the total population, and thus regard the newly born population as the only way of population input. We use $A_f$ and $A_m$ to represent female and male recruitment rate, respectively. There are two ways for susceptible individuals to become infected: non-sexual contact and sexual contact, and sexual contact is also divided into homosexual contact and heterosexual contact.

A schematic diagram of the model is given in Figure 1. We denote λ as the non-sexually transmitted infection rate, and ${\lambda }_{mf}$, ${\lambda }_{fm}$, ${\lambda }_{ff}$, ${\lambda }_{mm}$ as sexually transmitted infection rates, where ${\lambda }_{mf}$ is the infectious rate from male to female, ${\lambda }_{fm}$ is the infectious rate from female to male, ${\lambda }_{ff}$ is the infectious rate from female to female, and ${\lambda }_{mm}$ is the infectious rate from male to male. Following the definitions of infectious rates given by Mclean and Blumberg in [19], we have the expression of each infectious rate as follows: (1) $\{\begin{array}{l}{\displaystyle {\lambda }_{fm}=\frac{2{\beta }_{fm}{c}_{fm}(I_f+ϵC_f)}{N},{\lambda }_{mf}={\displaystyle \frac{2{\beta }_{mf}{c}_{mf}(I_m+ϵC_m)}{N},}}\\ {\displaystyle {\lambda }_{mm}=\frac{2{\beta }_{mm}{c}_{mm}(I_m+ϵC_m)}{N},{\lambda }_{ff}={\displaystyle \frac{2{\beta }_{ff}{c}_{ff}(I_f+ϵC_f)}{N},}}\\ {\displaystyle \lambda =\beta (I_m+I_f)+ϵ\beta (C_m+C_f).}\end{array}$(1)

Figure 1. Flowchart of HBV transmission in adults.

The parameters are given in Table 1.

Table 1. Model parameters.

Parameters	Value	Comments	Source
a	0.75	The proportion of heterosex in female susceptibles	[5]
b	0.75	The proportion of heterosex in male susceptibles	[5]
$A_f$	8,304,000 ${\mathrm{year}}^{-1}$	Recruitment rate of female	[39]
$A_m$	8,304,000 ${\mathrm{year}}^{-1}$	Recruitment rate of male	[39]
${\mu }_0$	$6.9\times {10}^{-3}$ ${\mathrm{year}}^{-1}$	Natural mortality rate	[20]
${\mu }_1$	$2\times {10}^{-3}$ ${\mathrm{year}}^{-1}$	HBV related mortality rate	[6]
ϵ	0.16	Reduced transmission rate	[9]
${\gamma }_1$	0.26 ${\mathrm{year}}^{-1}$	Progress rate from acute to carrier	Fitting
${\gamma }_2$	0.025 ${\mathrm{year}}^{-1}$	Progress rate from carrier to immune	[9]
${\gamma }_3$	0.05 ${\mathrm{year}}^{-1}$	Vaccination rate for adults	[39]
$q_f$	0.05	Average probability a female adult fails to clear	[19]
		an acute infection and develops to carrier state
$q_m$	0.07	Average probability a male adult fails to clear	[19]
		an acute infection and develops to carrier state
β	$3.5\times {10}^{-11}$ ${\mathrm{year}}^{-1}$	Transmission rate for non-sexual transmission	Fitting
${\beta }_{mf}$	$2.26\times {10}^{-7}$	The transmission probability for male per sex	Fitting
		contact with female
${\beta }_{fm}$	$3.06\times {10}^{-3}$	The transmission probability for female per sex	Fitting
		contact with male
${\beta }_{ff}$	0.0142	The transmission probability for female per sex	Fitting
		contact with female
${\beta }_{mm}$	$9.73\times {10}^{-3}$	The transmission probability for male per sex	Fitting
		contact with male
$c_{mf}$	21.8 ${\mathrm{year}}^{-1}$	The average number of sex contacts of females	[12]
		with males
$c_{fm}$	21.8 ${\mathrm{year}}^{-1}$	The average number of sex contacts of males	[12]
		with females
$c_{ff}$	21.8 ${\mathrm{year}}^{-1}$	The average number of sex contacts of females	[12,17]
		with females
$c_{mm}$	65.4 ${\mathrm{year}}^{-1}$	The average number of sex contacts of males	[12,31]
		with males

Based on the above discussion and the flowchart in Figure 1, the model is formulated as follows: (2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_f}{\mathrm{d}t}=A_f-(\lambda +{\mu }_0+{\gamma }_3)S_f-a{\lambda }_{mf}{S}_f-(1-a){\lambda }_{ff}{S}_f,}\\ {\displaystyle \frac{\mathrm{d}{I}_f}{\mathrm{d}t}=\lambda S_f+a{\lambda }_{mf}{S}_f+(1-a){\lambda }_{ff}{S}_f-({\mu }_0+{\gamma }_1)I_f,}\\ {\displaystyle \frac{\mathrm{d}{C}_f}{\mathrm{d}t}=q_f{\gamma }_1{I}_f-({\mu }_0+{\mu }_1+{\gamma }_2)C_f,}\\ {\displaystyle \frac{\mathrm{d}{S}_m}{\mathrm{d}t}=A_m-(\lambda +{\mu }_0+{\gamma }_3)S_m-b{\lambda }_{fm}{S}_m-(1-b){\lambda }_{mm}{S}_m,}\\ {\displaystyle \frac{\mathrm{d}{I}_m}{\mathrm{d}t}=\lambda S_m+b{\lambda }_{fm}{S}_m+(1-b){\lambda }_{mm}{S}_m-({\mu }_0+{\gamma }_1)I_m,}\\ {\displaystyle \frac{\mathrm{d}{C}_m}{\mathrm{d}t}=q_m{\gamma }_1{I}_m-({\mu }_0+{\mu }_1+{\gamma }_2)C_m,}\\ {\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}={\gamma }_3(S_f+S_m)+(1-q_f){\gamma }_1{I}_f+(1-q_m){\gamma }_1{I}_m+{\gamma }_2(C_f+C_m)-{\mu }_{0}R,}\end{array}$(2) where the parameters are described in Table 1 and λ, ${\lambda }_{fm},$ ${\lambda }_{mf},$ ${\lambda }_{mm},$ ${\lambda }_{ff}$ are given in (1(1) $\{\begin{array}{l}{\displaystyle {\lambda }_{fm}=\frac{2{\beta }_{fm}{c}_{fm}(I_f+ϵC_f)}{N},{\lambda }_{mf}={\displaystyle \frac{2{\beta }_{mf}{c}_{mf}(I_m+ϵC_m)}{N},}}\\ {\displaystyle {\lambda }_{mm}=\frac{2{\beta }_{mm}{c}_{mm}(I_m+ϵC_m)}{N},{\lambda }_{ff}={\displaystyle \frac{2{\beta }_{ff}{c}_{ff}(I_f+ϵC_f)}{N},}}\\ {\displaystyle \lambda =\beta (I_m+I_f)+ϵ\beta (C_m+C_f).}\end{array}$(1) ).

2.2. The disease-free equilibrium $E_0$

It is easy to see that the model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_f}{\mathrm{d}t}=A_f-(\lambda +{\mu }_0+{\gamma }_3)S_f-a{\lambda }_{mf}{S}_f-(1-a){\lambda }_{ff}{S}_f,}\\ {\displaystyle \frac{\mathrm{d}{I}_f}{\mathrm{d}t}=\lambda S_f+a{\lambda }_{mf}{S}_f+(1-a){\lambda }_{ff}{S}_f-({\mu }_0+{\gamma }_1)I_f,}\\ {\displaystyle \frac{\mathrm{d}{C}_f}{\mathrm{d}t}=q_f{\gamma }_1{I}_f-({\mu }_0+{\mu }_1+{\gamma }_2)C_f,}\\ {\displaystyle \frac{\mathrm{d}{S}_m}{\mathrm{d}t}=A_m-(\lambda +{\mu }_0+{\gamma }_3)S_m-b{\lambda }_{fm}{S}_m-(1-b){\lambda }_{mm}{S}_m,}\\ {\displaystyle \frac{\mathrm{d}{I}_m}{\mathrm{d}t}=\lambda S_m+b{\lambda }_{fm}{S}_m+(1-b){\lambda }_{mm}{S}_m-({\mu }_0+{\gamma }_1)I_m,}\\ {\displaystyle \frac{\mathrm{d}{C}_m}{\mathrm{d}t}=q_m{\gamma }_1{I}_m-({\mu }_0+{\mu }_1+{\gamma }_2)C_m,}\\ {\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}={\gamma }_3(S_f+S_m)+(1-q_f){\gamma }_1{I}_f+(1-q_m){\gamma }_1{I}_m+{\gamma }_2(C_f+C_m)-{\mu }_{0}R,}\end{array}$(2) ) has a disease-free equilibrium $E_0=(S_f^0,0,0,S_m^0,$ $0,0,R^0),$ where (3) $S_f^0=\frac{A_f}{{\mu }_0+{\gamma }_3},S_m^0=\frac{A_m}{{\mu }_0+{\gamma }_3},R^0=\frac{{\gamma }_3(S_f^0+S_m^0)}{{\mu }_0}=\frac{{\gamma }_3(A_f+A_m)}{{\mu }_0({\mu }_0+{\gamma }_3)}.$(3)

2.3. The basic reproduction number $R_0$

In this section we use the methods in [7,8] to calculate the basic reproduction number. Firstly, we let $X={\left[\begin{array}{c}{I}_f,I_m,C_f,C_m,S_f,S_m,R\end{array}\right]}^T$ and rewrite system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_f}{\mathrm{d}t}=A_f-(\lambda +{\mu }_0+{\gamma }_3)S_f-a{\lambda }_{mf}{S}_f-(1-a){\lambda }_{ff}{S}_f,}\\ {\displaystyle \frac{\mathrm{d}{I}_f}{\mathrm{d}t}=\lambda S_f+a{\lambda }_{mf}{S}_f+(1-a){\lambda }_{ff}{S}_f-({\mu }_0+{\gamma }_1)I_f,}\\ {\displaystyle \frac{\mathrm{d}{C}_f}{\mathrm{d}t}=q_f{\gamma }_1{I}_f-({\mu }_0+{\mu }_1+{\gamma }_2)C_f,}\\ {\displaystyle \frac{\mathrm{d}{S}_m}{\mathrm{d}t}=A_m-(\lambda +{\mu }_0+{\gamma }_3)S_m-b{\lambda }_{fm}{S}_m-(1-b){\lambda }_{mm}{S}_m,}\\ {\displaystyle \frac{\mathrm{d}{I}_m}{\mathrm{d}t}=\lambda S_m+b{\lambda }_{fm}{S}_m+(1-b){\lambda }_{mm}{S}_m-({\mu }_0+{\gamma }_1)I_m,}\\ {\displaystyle \frac{\mathrm{d}{C}_m}{\mathrm{d}t}=q_m{\gamma }_1{I}_m-({\mu }_0+{\mu }_1+{\gamma }_2)C_m,}\\ {\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}={\gamma }_3(S_f+S_m)+(1-q_f){\gamma }_1{I}_f+(1-q_m){\gamma }_1{I}_m+{\gamma }_2(C_f+C_m)-{\mu }_{0}R,}\end{array}$(2) ) as (4) $X^{\mathrm{\prime }}=\mathcal{F}-\mathcal{V},$(4) where $\begin{array}{rl}\mathcal{F}& =\left(\begin{array}{c}\lambda S_f+a{\lambda }_{mf}{S}_f+(1-a){\lambda }_{ff}{S}_f\\ \lambda S_m+b{\lambda }_{fm}{S}_m+(1-b){\lambda }_{mm}{S}_m\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right),\\ \mathcal{V}& =\left(\begin{array}{c}({\mu }_0+{\gamma }_1)I_f\\ ({\mu }_0+{\gamma }_1)I_m\\ ({\mu }_0+{\mu }_1+{\gamma }_2)C_f-q_f{\gamma }_1{I}_f\\ ({\mu }_0+{\mu }_1+{\gamma }_2)C_m-q_m{\gamma }_1{I}_m\\ -A_f+(\lambda +{\mu }_0+{\gamma }_3)S_f+a{\lambda }_{mf}{S}_f+(1-a){\lambda }_{ff}{S}_f\\ -A_m+(\lambda +{\mu }_0+{\gamma }_3)S_m+b{\lambda }_{fm}{S}_m+(1-b){\lambda }_{mm}{S}_m\\ -{\gamma }_3(S_f+S_m)-(1-q_f){\gamma }_1{I}_f-(1-q_m){\gamma }_1{I}_m-{\gamma }_2(C_f+C_m)+{\mu }_{0}R\end{array}\right).\end{array}$

Obviously, model (4(4) $X^{\mathrm{\prime }}=\mathcal{F}-\mathcal{V},$(4) ) has a disease-free equilibrium $\widetilde{E_0}=(0,0,0,0,S_f^0,S_m^0,R^0)$, and the linearization of model (4(4) $X^{\mathrm{\prime }}=\mathcal{F}-\mathcal{V},$(4) ) at the disease-free equilibrium $\widetilde{E_0}$ is given by: $D\mathcal{F}(\widetilde{E_0})-D\mathcal{V}(\widetilde{E_0})=\left(\begin{array}{cc}F-V& 0\\ -H& J\end{array}\right),$where H and J are $3\times 3$ matrices, and $F=\left(\begin{array}{cc}{F}_1& F_2\\ 0& 0\end{array}\right)$, $F_1=\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)$, $F_2=ϵF_1$, $V=\left(\begin{array}{cc}{D}_1& 0\\ -D_2& V_2\end{array}\right)$, $D_1=({\mu }_0+{\gamma }_1)I_{2\times 2}$, $D_2=\left(\begin{array}{cc}{q}_f{\gamma }_1& 0\\ 0& q_m{\gamma }_1\end{array}\right)$, $V_2=({\mu }_0+{\mu }_1+{\gamma }_2)I_{2\times 2}$, $\begin{array}{rl}A& =\beta S_f^0+\frac{2(1-a){\beta }_{ff}{c}_{ff}{S}_f^0}{S_f^0+S_m^0+R^0},B=\beta S_f^0+\frac{2a{\beta }_{mf}{c}_{mf}{S}_f^0}{S_f^0+S_m^0+R^0},\\ C& =\beta S_m^0+\frac{2b{\beta }_{fm}{}_{fm}{S}_m^0}{S_f^0+S_m^0+R^0},D=\beta S_m^0+\frac{2(1-b){\beta }_{mm}{c}_{mm}{S}_m^0}{S_f^0+S_m^0+R^0}.\end{array}$

By [7,8], we know that the basic reproduction number of system (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_f}{\mathrm{d}t}=A_f-(\lambda +{\mu }_0+{\gamma }_3)S_f-a{\lambda }_{mf}{S}_f-(1-a){\lambda }_{ff}{S}_f,}\\ {\displaystyle \frac{\mathrm{d}{I}_f}{\mathrm{d}t}=\lambda S_f+a{\lambda }_{mf}{S}_f+(1-a){\lambda }_{ff}{S}_f-({\mu }_0+{\gamma }_1)I_f,}\\ {\displaystyle \frac{\mathrm{d}{C}_f}{\mathrm{d}t}=q_f{\gamma }_1{I}_f-({\mu }_0+{\mu }_1+{\gamma }_2)C_f,}\\ {\displaystyle \frac{\mathrm{d}{S}_m}{\mathrm{d}t}=A_m-(\lambda +{\mu }_0+{\gamma }_3)S_m-b{\lambda }_{fm}{S}_m-(1-b){\lambda }_{mm}{S}_m,}\\ {\displaystyle \frac{\mathrm{d}{I}_m}{\mathrm{d}t}=\lambda S_m+b{\lambda }_{fm}{S}_m+(1-b){\lambda }_{mm}{S}_m-({\mu }_0+{\gamma }_1)I_m,}\\ {\displaystyle \frac{\mathrm{d}{C}_m}{\mathrm{d}t}=q_m{\gamma }_1{I}_m-({\mu }_0+{\mu }_1+{\gamma }_2)C_m,}\\ {\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}={\gamma }_3(S_f+S_m)+(1-q_f){\gamma }_1{I}_f+(1-q_m){\gamma }_1{I}_m+{\gamma }_2(C_f+C_m)-{\mu }_{0}R,}\end{array}$(2) ) is the spectral radius of the matrix $F{V}^{-1}$, denoted by $\rho (F{V}^{-1})$, which is the maximum of the module of the eigenvalues of $F{V}^{-1}$. By simple calculation, we have $F{V}^{-1}=(F_1+F_2{V}_2^{-1}{D}_2)D_1^{-1}=\frac{1}{{\mu }_0+{\mu }_1}\left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right),$where $\begin{array}{l}{a}_{11}={\displaystyle A(1+\frac{ϵq_f{\gamma }_1}{{\mu }_0+{\mu }_1+{\gamma }_2}),a_{12}=B(1+\frac{ϵq_m{\gamma }_1}{{\mu }_0+{\mu }_1+{\gamma }_2}),}\\ a_{21}={\displaystyle C(1+\frac{ϵq_f{\gamma }_1}{{\mu }_0+{\mu }_1+{\gamma }_2}),a_{22}=D(1+\frac{ϵq_m{\gamma }_1}{{\mu }_0+{\mu }_1+{\gamma }_2}).}\end{array}$The eigenvalues of $F{V}^{-1}$ satisfy the following characteristic equation (5) ${\lambda }^2-\frac{a_{11}+a_{22}}{{\mu }_0+{\gamma }_1}\lambda +\frac{a_{11}{a}_{22}-a_{12}{a}_{21}}{({\mu }_0+{\gamma }_1{)}^2}=0,$(5) which has eigenvalues ${\lambda }_{1,2}=\frac{1}{{\mu }_0+{\gamma }_1}\frac{a_{11}+a_{22}\pm \sqrt{(a_{11}+a_{22}{)}^2-4(a_{11}{a}_{22}-a_{12}{a}_{21})}}{2}=\frac{r_2\pm \sqrt{r_2^2-4r_1{r}_3}}{2r_1},$ where $\begin{array}{rl}{r}_1& ={({\gamma }_1+{\mu }_0)}^2{({\mu }_0+{\mu }_1+{\gamma }_2)}^2{({\mu }_0+{\gamma }_3)}^2,\\ r_2& =({\gamma }_1+{\mu }_0)({\mu }_0+{\mu }_1+{\gamma }_2)({\mu }_0+{\gamma }_3)\\ & \times [A_f(\beta +\frac{2(1-a){\beta }_{ff}{c}_{ff}}{N^0})({\mu }_0+{\mu }_1+{\gamma }_2+ϵq_f{\gamma }_1)\\ & +A_m(\beta +{\displaystyle \frac{2(1-b){\beta }_{mm}{c}_{mm}}{N^0}})({\mu }_0+{\mu }_1+{\gamma }_2+ϵq_m{\gamma }_1)],\\ r_3& =A_f{A}_m({\mu }_0+{\mu }_1+{\gamma }_2+ϵq_f{\gamma }_1)({\mu }_0+{\mu }_1+{\gamma }_2+ϵq_m{\gamma }_1)\\ & \{(\beta +{\displaystyle \frac{2(1-a){\beta }_{ff}{c}_{ff}}{N^0}})(\beta +{\displaystyle \frac{2(1-b){\beta }_{mm}{c}_{mm}}{N^0}})\\ & -(\beta +{\displaystyle \frac{2a{\beta }_{mf}{c}_{mf}}{N^0}})(\beta +{\displaystyle \frac{2b{\beta }_{fm}{c}_{fm}}{N^0}})\},\\ N^0& =S_f^0+S_m^0+R^0=\frac{A_f+A_m}{{\mu }_0}.\end{array}$

Thus, the basic reproduction number can be expressed as (6) $R_0={\displaystyle \frac{r_2+\sqrt{r_2^2-4r_1{r}_3}}{2r_1}.}$(6) Following from Theorem 2 of [8], we obtain that the disease-free equilibrium $E_0$ is asymptotically stable if $R_0<1$, while it is unstable if $R_0>1$.

3. Simulations and sensitivity analysis

3.1. Simulations

In order to simulate infected proportions of both female and male adults, we now let $x_1=\frac{S_f}{N},x_2=\frac{I_f}{N},x_3=\frac{C_f}{N},x_4=\frac{S_m}{N},x_5=\frac{I_m}{N},x_6=\frac{C_m}{N},x_7=\frac{R}{N}.$Then model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_f}{\mathrm{d}t}=A_f-(\lambda +{\mu }_0+{\gamma }_3)S_f-a{\lambda }_{mf}{S}_f-(1-a){\lambda }_{ff}{S}_f,}\\ {\displaystyle \frac{\mathrm{d}{I}_f}{\mathrm{d}t}=\lambda S_f+a{\lambda }_{mf}{S}_f+(1-a){\lambda }_{ff}{S}_f-({\mu }_0+{\gamma }_1)I_f,}\\ {\displaystyle \frac{\mathrm{d}{C}_f}{\mathrm{d}t}=q_f{\gamma }_1{I}_f-({\mu }_0+{\mu }_1+{\gamma }_2)C_f,}\\ {\displaystyle \frac{\mathrm{d}{S}_m}{\mathrm{d}t}=A_m-(\lambda +{\mu }_0+{\gamma }_3)S_m-b{\lambda }_{fm}{S}_m-(1-b){\lambda }_{mm}{S}_m,}\\ {\displaystyle \frac{\mathrm{d}{I}_m}{\mathrm{d}t}=\lambda S_m+b{\lambda }_{fm}{S}_m+(1-b){\lambda }_{mm}{S}_m-({\mu }_0+{\gamma }_1)I_m,}\\ {\displaystyle \frac{\mathrm{d}{C}_m}{\mathrm{d}t}=q_m{\gamma }_1{I}_m-({\mu }_0+{\mu }_1+{\gamma }_2)C_m,}\\ {\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}={\gamma }_3(S_f+S_m)+(1-q_f){\gamma }_1{I}_f+(1-q_m){\gamma }_1{I}_m+{\gamma }_2(C_f+C_m)-{\mu }_{0}R,}\end{array}$(2) ) becomes as the following “proportional” model (7) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=\frac{A_f}{N}-\frac{A_f+A_m}{N}{x}_1-[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N+{\gamma }_3]x_1}\\ -2a{\beta }_{mf}{c}_{mf}(x_5+ϵx_6)x_1\\ {\displaystyle -2(1-a){\beta }_{ff}{c}_{ff}(x_2+ϵx_3)x_1+{\mu }_1{x}_1(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N]x_1}\\ +2a{\beta }_{mf}{c}_{mf}(x_5+ϵx_6)x_1+2(1-a){\beta }_{ff}{c}_{ff}(x_2+ϵx_3)x_1\\ {\displaystyle -{\gamma }_1{x}_2-\frac{A_f+A_m}{N}{x}_2+{\mu }_1{x}_2(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_3}{\mathrm{d}t}=q_f{\gamma }_1{x}_2-({\mu }_1+{\gamma }_2)x_3-\frac{A_f+A_m}{N}{x}_3+{\mu }_1{x}_3(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_4}{\mathrm{d}t}=\frac{A_m}{N}-\frac{A_f+A_m}{N}{x}_4-[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N+{\gamma }_3]x_4}\\ -2b{\beta }_{fm}{c}_{fm}(x_2+ϵx_3)x_4\\ {\displaystyle -2(1-b){\beta }_{mm}{c}_{mm}(x_5+ϵx_6)x_4+{\mu }_1{x}_4(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_5}{\mathrm{d}t}=[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N]x_4+2b{\beta }_{fm}{c}_{fm}(x_2+ϵx_3)x_4}\\ +2(1-b){\beta }_{mm}{c}_{mm}(x_5+ϵx_6)x_4\\ {\displaystyle -{\gamma }_1{x}_5-\frac{A_f+A_m}{N}{x}_5+{\mu }_1{x}_5(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_6}{\mathrm{d}t}=q_m{\gamma }_1{x}_5-({\mu }_1+{\gamma }_2)x_6-\frac{A_f+A_m}{N}{x}_6+{\mu }_1{x}_6(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}N}{\mathrm{d}t}=A_f+A_m-{\mu }_{0}N-{\mu }_1(x_3+x_6)N,}\end{array}$(7) where $A_f$ in Table 1 represents the recruitment rate of the whole susceptible adult female population, $A_m$ represents the recruitment rate of the whole susceptible adult male population. The calculation method is as follows [39]: first, the average number of Chinese under-aged population from 2005 to 2016 was 332, 160, 000, growth rate from under-aged to adult was 1/20, and the proportion of men and women was 1/2, respectively. Then we can obtain $A_f$ and $A_m$ as 8, 304, 000 $(=332,160,000\times 1/20\times 1/2)$. From Chinese Center for Disease Control and Prevention, we obtain the data of hepatitis B among adult males and females from 2005 to 2016, and the total numbers of adult males and females are obtained from the National Bureau of Statistics. For the convenience of classification, any person over twenty years old is counted as an adult. In order to get the infected proportions of both female and male adults, we divide the numbers of female and male adults with hepatitis B by the total number of adults, respectively. The specific infected proportions for both female and male adults are shown in Table 2.

Table 2. Acute infected proportions of female and male adults in China.

Year	2005	2006	2007	2008	2009	2010	2011
Male	0.000575071	0.000645782	0.000683863	0.000685808	0.00068905	0.000618111	0.000629059
Female	0.000296771	0.000349889	0.000377479	0.000383057	0.000389468	0.000353505	0.000374578
Year	2012	2013	2014	2015	2016
Male	0.000622946	0.000550488	0.000534022	0.000536865	0.000539815
Female	0.000373096	0.00032709	0.000315454	0.000307392	0.000310499

We estimate the unknown parameters by calculating the minimum sum of: (8) $min\sum _{i=2005}^{2016}\frac{(Dat{a}_i-Simulatio{n}_i{)}^2}{Simulatio{n}_i},$(8) with MATLAB tool fminsearch. (8(8) $min\sum _{i=2005}^{2016}\frac{(Dat{a}_i-Simulatio{n}_i{)}^2}{Simulatio{n}_i},$(8) ) includes the sum over infected proportions of female and male adults in Table 2. First, the initial values of $x_1(0)=\frac{S_f(0)}{N(0)}=0.42$, $x_2(0)=\frac{I_f(0)}{N(0)}=0.000296771$, $x_4(0)=\frac{S_m(0)}{N(0)}=0.47$, $x_5(0)=\frac{I_m(0)}{N(0)}=0.000575071$ and $N(0)=931,400,000$ are calculated from the data of China CDC's website. Then $x_3(0)=\frac{C_f(0)}{N(0)}=8.09\times {10}^{-3}$, $x_6(0)=\frac{C_m(0)}{N(0)}=7.3\times {10}^{-4}$ are the optimal unknown initial values fitted by MATLAB. Moreover, we can also get some of the optimal unknown parameter values: ${\gamma }_1=0.26$, $\beta =3.5\times {10}^{-11}$, ${\beta }_{mf}=2.26\times {10}^{-7}$, ${\beta }_{fm}=3.06\times {10}^{-3}$, ${\beta }_{ff}=0.0142$, ${\beta }_{mm}=9.73\times {10}^{-3}$. The degree of freedom is 15, we get that the Chi-Square value at the 95$\mathrm{\%}$ significance level should be less than 24.996 by searching for the table in [10]. By looking for the parameter intervals whose Chi-Square is less than 24.996, we estimate the confidence interval of the unknown initial values and parameters one by one (see Table 3).

Table 3. The 95$\mathrm{\%}$ confidence intervals for unknown parameters and initial values.

Parameters	95$\mathrm{\%}$ Confidence interval	Initial values	95$\mathrm{\%}$ Confidence interval
${\gamma }_1$	$[0.154,0.388]$	$\frac{C_f(0)}{N(0)}$	$[2.55\times {10}^{-3},0.0178]$
β	$[0,1.17\times {10}^{-10}]$	$\frac{C_m(0)}{N(0)}$	$[0,0.00777]$
${\beta }_{mf}$	$[0,0.0147]$
${\beta }_{fm}$	$[0,0.0101]$
${\beta }_{ff}$	$[0.0024,0.0335]$
${\beta }_{mm}$	$[0,0.0191]$

In order to predict the $95\mathrm{\%}$ confidence interval of the model, we assume that the results of the model conform to the normal distribution, and the mean value of the normal distribution is the simulated value of the model. Therefore, we only need to know the variance of the distribution. To calculate the standard deviation, 10,000 points are randomly selected in the $95\mathrm{\%}$ confidence interval of all parameters, and the error of each point is $\sqrt{\frac{1}{22}\sum _{2005}^{2016}(Dat{a}_i-Simulatio{n}_i{)}^2}.$Here, 22 is chosen because there is no error when the first values of infected proportions of female and male adults are taken as the initial values. The mean value of 10,000 errors is taken as the standard deviation, which is $\sigma =1.36\times {10}^{-4}$. Then the $95\mathrm{\%}$ confidence interval of the prediction model is $[max\{0,Simulatio{n}_i-1.645\sigma \},Simulatio{n}_i+1.645\sigma ]$, $i=2004,2005,\dots ,2035$ [10].

Based on the parameter values given in Table 1, we calculate the basic reproduction number $R_0=0.107<1$, the disease-free equilibrium of model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_f}{\mathrm{d}t}=A_f-(\lambda +{\mu }_0+{\gamma }_3)S_f-a{\lambda }_{mf}{S}_f-(1-a){\lambda }_{ff}{S}_f,}\\ {\displaystyle \frac{\mathrm{d}{I}_f}{\mathrm{d}t}=\lambda S_f+a{\lambda }_{mf}{S}_f+(1-a){\lambda }_{ff}{S}_f-({\mu }_0+{\gamma }_1)I_f,}\\ {\displaystyle \frac{\mathrm{d}{C}_f}{\mathrm{d}t}=q_f{\gamma }_1{I}_f-({\mu }_0+{\mu }_1+{\gamma }_2)C_f,}\\ {\displaystyle \frac{\mathrm{d}{S}_m}{\mathrm{d}t}=A_m-(\lambda +{\mu }_0+{\gamma }_3)S_m-b{\lambda }_{fm}{S}_m-(1-b){\lambda }_{mm}{S}_m,}\\ {\displaystyle \frac{\mathrm{d}{I}_m}{\mathrm{d}t}=\lambda S_m+b{\lambda }_{fm}{S}_m+(1-b){\lambda }_{mm}{S}_m-({\mu }_0+{\gamma }_1)I_m,}\\ {\displaystyle \frac{\mathrm{d}{C}_m}{\mathrm{d}t}=q_m{\gamma }_1{I}_m-({\mu }_0+{\mu }_1+{\gamma }_2)C_m,}\\ {\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}={\gamma }_3(S_f+S_m)+(1-q_f){\gamma }_1{I}_f+(1-q_m){\gamma }_1{I}_m+{\gamma }_2(C_f+C_m)-{\mu }_{0}R,}\end{array}$(2) ) is $E_0(1.46\times {10}^8,0,0,1.46\times {10}^8,0,0,2.12\times {10}^9)$, and model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{S}_f}{\mathrm{d}t}=A_f-(\lambda +{\mu }_0+{\gamma }_3)S_f-a{\lambda }_{mf}{S}_f-(1-a){\lambda }_{ff}{S}_f,}\\ {\displaystyle \frac{\mathrm{d}{I}_f}{\mathrm{d}t}=\lambda S_f+a{\lambda }_{mf}{S}_f+(1-a){\lambda }_{ff}{S}_f-({\mu }_0+{\gamma }_1)I_f,}\\ {\displaystyle \frac{\mathrm{d}{C}_f}{\mathrm{d}t}=q_f{\gamma }_1{I}_f-({\mu }_0+{\mu }_1+{\gamma }_2)C_f,}\\ {\displaystyle \frac{\mathrm{d}{S}_m}{\mathrm{d}t}=A_m-(\lambda +{\mu }_0+{\gamma }_3)S_m-b{\lambda }_{fm}{S}_m-(1-b){\lambda }_{mm}{S}_m,}\\ {\displaystyle \frac{\mathrm{d}{I}_m}{\mathrm{d}t}=\lambda S_m+b{\lambda }_{fm}{S}_m+(1-b){\lambda }_{mm}{S}_m-({\mu }_0+{\gamma }_1)I_m,}\\ {\displaystyle \frac{\mathrm{d}{C}_m}{\mathrm{d}t}=q_m{\gamma }_1{I}_m-({\mu }_0+{\mu }_1+{\gamma }_2)C_m,}\\ {\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}={\gamma }_3(S_f+S_m)+(1-q_f){\gamma }_1{I}_f+(1-q_m){\gamma }_1{I}_m+{\gamma }_2(C_f+C_m)-{\mu }_{0}R,}\end{array}$(2) ) has no positive equilibrium under this set of parameters, which implies the disease-free equilibrium $E_0$ is globally asymptotically stable. This also means that the number of adults patients with hepatitis B in China will be gradually reduced. Next, we use model (7(7) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=\frac{A_f}{N}-\frac{A_f+A_m}{N}{x}_1-[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N+{\gamma }_3]x_1}\\ -2a{\beta }_{mf}{c}_{mf}(x_5+ϵx_6)x_1\\ {\displaystyle -2(1-a){\beta }_{ff}{c}_{ff}(x_2+ϵx_3)x_1+{\mu }_1{x}_1(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N]x_1}\\ +2a{\beta }_{mf}{c}_{mf}(x_5+ϵx_6)x_1+2(1-a){\beta }_{ff}{c}_{ff}(x_2+ϵx_3)x_1\\ {\displaystyle -{\gamma }_1{x}_2-\frac{A_f+A_m}{N}{x}_2+{\mu }_1{x}_2(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_3}{\mathrm{d}t}=q_f{\gamma }_1{x}_2-({\mu }_1+{\gamma }_2)x_3-\frac{A_f+A_m}{N}{x}_3+{\mu }_1{x}_3(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_4}{\mathrm{d}t}=\frac{A_m}{N}-\frac{A_f+A_m}{N}{x}_4-[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N+{\gamma }_3]x_4}\\ -2b{\beta }_{fm}{c}_{fm}(x_2+ϵx_3)x_4\\ {\displaystyle -2(1-b){\beta }_{mm}{c}_{mm}(x_5+ϵx_6)x_4+{\mu }_1{x}_4(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_5}{\mathrm{d}t}=[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N]x_4+2b{\beta }_{fm}{c}_{fm}(x_2+ϵx_3)x_4}\\ +2(1-b){\beta }_{mm}{c}_{mm}(x_5+ϵx_6)x_4\\ {\displaystyle -{\gamma }_1{x}_5-\frac{A_f+A_m}{N}{x}_5+{\mu }_1{x}_5(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_6}{\mathrm{d}t}=q_m{\gamma }_1{x}_5-({\mu }_1+{\gamma }_2)x_6-\frac{A_f+A_m}{N}{x}_6+{\mu }_1{x}_6(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}N}{\mathrm{d}t}=A_f+A_m-{\mu }_{0}N-{\mu }_1(x_3+x_6)N,}\end{array}$(7) ) to simulate the data from 2005 to 2016 and predict the trends of infected proportions of female and male adults in China.

In Figure 2(a ,b), the solid yellow curves with circles represent the solution curves of infected proportions of male and female adults of model (7(7) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=\frac{A_f}{N}-\frac{A_f+A_m}{N}{x}_1-[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N+{\gamma }_3]x_1}\\ -2a{\beta }_{mf}{c}_{mf}(x_5+ϵx_6)x_1\\ {\displaystyle -2(1-a){\beta }_{ff}{c}_{ff}(x_2+ϵx_3)x_1+{\mu }_1{x}_1(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N]x_1}\\ +2a{\beta }_{mf}{c}_{mf}(x_5+ϵx_6)x_1+2(1-a){\beta }_{ff}{c}_{ff}(x_2+ϵx_3)x_1\\ {\displaystyle -{\gamma }_1{x}_2-\frac{A_f+A_m}{N}{x}_2+{\mu }_1{x}_2(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_3}{\mathrm{d}t}=q_f{\gamma }_1{x}_2-({\mu }_1+{\gamma }_2)x_3-\frac{A_f+A_m}{N}{x}_3+{\mu }_1{x}_3(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_4}{\mathrm{d}t}=\frac{A_m}{N}-\frac{A_f+A_m}{N}{x}_4-[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N+{\gamma }_3]x_4}\\ -2b{\beta }_{fm}{c}_{fm}(x_2+ϵx_3)x_4\\ {\displaystyle -2(1-b){\beta }_{mm}{c}_{mm}(x_5+ϵx_6)x_4+{\mu }_1{x}_4(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_5}{\mathrm{d}t}=[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N]x_4+2b{\beta }_{fm}{c}_{fm}(x_2+ϵx_3)x_4}\\ +2(1-b){\beta }_{mm}{c}_{mm}(x_5+ϵx_6)x_4\\ {\displaystyle -{\gamma }_1{x}_5-\frac{A_f+A_m}{N}{x}_5+{\mu }_1{x}_5(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_6}{\mathrm{d}t}=q_m{\gamma }_1{x}_5-({\mu }_1+{\gamma }_2)x_6-\frac{A_f+A_m}{N}{x}_6+{\mu }_1{x}_6(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}N}{\mathrm{d}t}=A_f+A_m-{\mu }_{0}N-{\mu }_1(x_3+x_6)N,}\end{array}$(7) ), respectively. We can see that these two solution curves match the real data curves. The infected proportions of female and male adults exhibit a decreasing trend, and so do the trend of the development of hepatitis B in the next 20 years (see Figure 3(a ,b)). That is to say, by implementing the current control measures, such as strengthening the regulation of blood circulation, raising the awareness of the public, the number of adults patients with hepatitis B in China will be gradually reduced.

Figure 2. Simulations of reported data on the prevalence of hepatitis B in Chinese adults. The solid blue curves represent the reported data while the solid yellow curves with circles are simulated by using model (7(7) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1}{\mathrm{d}t}=\frac{A_f}{N}-\frac{A_f+A_m}{N}{x}_1-[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N+{\gamma }_3]x_1}\\ -2a{\beta }_{mf}{c}_{mf}(x_5+ϵx_6)x_1\\ {\displaystyle -2(1-a){\beta }_{ff}{c}_{ff}(x_2+ϵx_3)x_1+{\mu }_1{x}_1(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N]x_1}\\ +2a{\beta }_{mf}{c}_{mf}(x_5+ϵx_6)x_1+2(1-a){\beta }_{ff}{c}_{ff}(x_2+ϵx_3)x_1\\ {\displaystyle -{\gamma }_1{x}_2-\frac{A_f+A_m}{N}{x}_2+{\mu }_1{x}_2(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_3}{\mathrm{d}t}=q_f{\gamma }_1{x}_2-({\mu }_1+{\gamma }_2)x_3-\frac{A_f+A_m}{N}{x}_3+{\mu }_1{x}_3(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_4}{\mathrm{d}t}=\frac{A_m}{N}-\frac{A_f+A_m}{N}{x}_4-[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N+{\gamma }_3]x_4}\\ -2b{\beta }_{fm}{c}_{fm}(x_2+ϵx_3)x_4\\ {\displaystyle -2(1-b){\beta }_{mm}{c}_{mm}(x_5+ϵx_6)x_4+{\mu }_1{x}_4(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_5}{\mathrm{d}t}=[\beta (x_2+x_5)N+ϵ\beta (x_3+x_6)N]x_4+2b{\beta }_{fm}{c}_{fm}(x_2+ϵx_3)x_4}\\ +2(1-b){\beta }_{mm}{c}_{mm}(x_5+ϵx_6)x_4\\ {\displaystyle -{\gamma }_1{x}_5-\frac{A_f+A_m}{N}{x}_5+{\mu }_1{x}_5(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}{x}_6}{\mathrm{d}t}=q_m{\gamma }_1{x}_5-({\mu }_1+{\gamma }_2)x_6-\frac{A_f+A_m}{N}{x}_6+{\mu }_1{x}_6(x_3+x_6),}\\ {\displaystyle \frac{\mathrm{d}N}{\mathrm{d}t}=A_f+A_m-{\mu }_{0}N-{\mu }_1(x_3+x_6)N,}\end{array}$(7) ), the vertical segments show the 95$\mathrm{\%}$ confidence intervals of infected proportions of male and female adults.

Figure 3. Predictions of the infected proportions of male and female adults in China in the next 20 years.

The effect of each transmission probability or rate on infected proportions of female and male adults is shown in Figure 4. Specifically, in Figure 4(a ,b), we perturb an optimal value β to $\beta \pm 2\times {10}^{-11}$. In Figure 4(c,e ,f), we perturb optimal values ${\beta }_{fm}$, ${\beta }_{mm}$, ${\beta }_{ff}$ (i.e. ${\beta }_{ij}$) to ${\beta }_{ij}\pm 2\times {10}^{-3}$, respectively. In particular, we perturb an optimal value ${\beta }_{mf}$ to ${\beta }_{mf}+2\times {10}^{-3}$ and ${\beta }_{mf}-2.26\times {10}^{-7}=0$ since the optimal value ${\beta }_{mf}$ is less than the disturbance value $2\times {10}^{-3}$. From Figure 4 and the disturbance value, we can see that β (the non-sexual transmission rate) is more important than ${\beta }_{ij}$ (the sexual transmission rates); that is, if β increases slightly, then hepatitis B will break out. Therefore, under the current situation even though we cannot reduce the non-sexual transmission rate, in order to prevent the outbreak of disease we still need to carefully control the increase of the non-sexual transmission rate.

Figure 4. Simulations of the prevalence of hepatitis B in Chinese male and female adults in terms of different parameters: (a) (b) β (the non-sexual transmission rate); (c) ${\beta }_{fm}$ (the transmission rate from females to males); (d) ${\beta }_{mf}$ (the transmission rate from males to females); (e) ${\beta }_{mm}$ (the transmission rate from males to males); (f) ${\beta }_{ff}$ (the transmission rate from females to females).

Figure 4(c–e) shows that decreasing ${\beta }_{fm}$ (the transmission rate from female to male) and ${\beta }_{mm}$ (the transmission rate from male to male) can dramatically control hepatitis B in males. Although, ${\beta }_{fm}$ (the transmission rate from female to male) is more important than ${\beta }_{mm}$ (the transmission rate of male to male) to control hepatitis B in males, the influence of the homosexual transmission for male on the control of hepatitis B should not be ignored. Moreover, in Figure 4(d–f) we can see that ${\beta }_{mf}$ (the transmission rate from male to female) and ${\beta }_{ff}$ (the transmission rate from female to female) are equally important to control hepatitis B in females. Therefore, reducing sexual transmission of hepatitis B is an effective way to control hepatitis B.

The effect of changes in the initial values of hepatitis B carrier proportions on the prevalence of hepatitis B in males and females respectively is shown in Figure 5. We can see that changing the initial values of $\frac{C_m(0)}{N(0)}$ and $\frac{C_f(0)}{N(0)}$ has a significant effect on adult prevalence, especially male prevalence. Thus, great attention should be paid to adult carriers as well.

Figure 5. Simulations of the prevalence of hepatitis B in Chinese male and female adults in terms of different initial conditions: (a) (b) $\frac{C_m(0)}{N(0)}$ (male carrier proportions); (c) (d) $\frac{C_f(0)}{N(0)}$ (female carrier proportions).

3.2. The sensitivity analysis of $R_0$

Now we analyse the sensitivity of $R_0$ with respect to various model parameters. Figure 6 shows that β (the non-sexual transmission rate) is more important than ${\beta }_{ij}$ (the sexual transmission rate). In addition, Figure 6(a) shows that the influences of ${\beta }_{fm}$ and ${\beta }_{mm}$ on $R_0$ are similar. However, the situation changes as the values of the parameters increase, Figure 6(b) shows that controlling ${\beta }_{mm}$ is more effective than controlling ${\beta }_{fm}$ for male if they are large. It indicates that the control of homosexual transmission for males with hepatitis B virus is necessary. In Figure 6(c), we see that controlling ${\beta }_{mf}$ to female is as important as controlling ${\beta }_{ff}$ if they are very small. However, Figure 6(d) shows that the influence of ${\beta }_{ff}$ on $R_0$ is the strongest when the parameters are large. Thus, in order to control and prevent the infection of HBV in China, great attention should be paid to the homosexual transmission of the virus in adults.

Figure 6. Plots of the basic reproduction number $R_0$ in terms of (a) ${\beta }_{fm}$, ${\beta }_{mm}$ and β on a small scale, (b) ${\beta }_{fm}$, ${\beta }_{mm}$ and β on a large scale, (c) ${\beta }_{mf}$, ${\beta }_{ff}$ and β on a small scale, and (d) ${\beta }_{mf}$, ${\beta }_{ff}$ and β on a large scale.

Figure 7(a,b) shows that the basic reproduction number $R_0$ can be lowered by increasing either ${\gamma }_1$ (the progress rate from acute to carrier) or ${\gamma }_3$ (vaccination rate) for adults. In Figure 7(c,d), the basic reproduction number $R_0$ can be lowered by decreasing $q_m$ (average probability a male adult fails to clear an acute infection) and $q_f$ (average probability a female adult fails to clear an acute infection). However, the influence of $q_m$ and $q_f$ on $R_0$ is very small and it is very difficult to control.

Figure 7. Plots of the basic reproduction number $R_0$ in terms of (a) ${\gamma }_1$ (progress rate from acute to carrier), (b) ${\gamma }_3$ (vaccination rate for adults), (c) $q_m$ (average probability a male adult fails to clear an acute infection and develops to carrier state), and (d) $q_f$ (average probability a female adult fails to clear an acute infection and develops to carrier state).

4. Discussion

Hepatitis B has been one of the major health problems for decades in China. However, the prevalence and the routes of infections have been changing. In 1980s, more children were infected and the main transmission was from mothers to babies [37]. In the last 10 years, HBV incidence rate has decreased significantly and steadily in children because of the intensive and strong immunization campaign and programme by local and central governments. In 2016, the WHO set a goal to eliminate viral hepatitis as a public health threat by 2030 (reducing new infections by $90\mathrm{\%}$ and mortality by $65\mathrm{\%}$) [29]. Various studies have reported that sexual transmission, both homosexually and heterosexually, is an important route of HBV spread in adults [4,26,28,32,39]. Previous mathematical modelling studies on the transmission of HBV in China focused on the heterosexual transmission of the virus [39]. In this study, we proposed a mathematical model to study the sexual transmission of HBV in China in which both homosexual and heterosexual transmissions are included. In order to focus on the sexual spread of the virus, we ignored the loss of immunity after immunization and the immigration and emigration of the total population.

The data in China from 2005 to 2016 show that the incidence of HBV in men is still much higher than that of women (see Table 2). From the sensitivity analysis, we see that it is important to increase the immunization rate for high-risk adults. Different from most other sexually transmitted infections, HBV vaccine is included in the Chinese National Vaccination Programs. Hepatitis B vaccination is recommended for all unvaccinated, uninfected persons being evaluated or treated for an STD [25]. Especially, all MSM should be tested for HBsAg to detect chronic HBV infection. The sensitivity analysis also indicates that it is important and crucial to control the sexual transmission for both heterosexual and homosexual activities. The effective control measures include enhancing public education and awareness about hepatitis B virus, particularly for MSM and WSW. Sex partners of persons with HBsAg should be counselled to use latex condoms to protect themselves from sexual exposure to infectious body fluids (including semen and vaginal secretions), unless they have been demonstrated to be immune after vaccination (anti-HBs $\ge 10$ mIU/mL) or previously infected (anti-HBc positive) [28]. Several condoms for females are also available, which can provide protection from transmission of HBV as well as other sexual transmitted diseases [28].

The co-infection of HBV and HIV needs further investigation. On one hand, HBV and HIV share similar transmission routes. On the other hand, HIV infection can impair the immune response to hepatitis B vaccination. One more challenge with co-infection is that cross-resistance between HIV and HBV drugs increased liver injury [30]. More attention should be paid to prevent HBV infection within the HIV infections. For example, it is recommended that people with HIV infection be tested for anti-HBs 1–2 months after the third vaccine dose [28].

Moreover, HBV, HCV and HDV share the similar transmission route. It is estimated that $5\mathrm{\%}$ of HBsAg-positive carriers, i.e. approximately 15 million people, are co-infection with HDV worldwide. For HBV-infections, hepatitis C virus co-infection accelerates progression of liver disease and increases the risk of HCV. In China, HBV infected populations are also at risk of TB infection [30]. The dynamics of HBV/HDV co-infection, HBV/HCV co-infection, as well as the HBV/TB co-infection are more complicated and design further investigation.

Acknowledgments

The authors would like to thank the editor and anonymous reviewers for their helpful comments and suggestions which helped us to improve the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Research was partially supported by NSFC [grant numbers 11771168, 11831012, 11871235] and the Fundamental Research Funds for the Central Universities [grant number CCNU19TS030].