Abstract
In this study, we formulate a model for hepatitis B virus with control strategies of newborn vaccination and treatment. Mathematical analysis is done theoretically and numerically. The results indicate that the stability of equilibria and persistence of the disease are determined by the basic reproductive number . Using the least squares method, the model is applied to simulate yearly new infected cases of hepatitis B in China from 2004 to 2016. Moreover, optimal control problem with newborn vaccination and treatment appearing as functions of time is analyzed by classical optimal theory. The existence of the solution to optimality system is proved, and the simulations are conducted to show the results when optimal control or current intervention is used.
1. Introduction
It has been over 50 years since the identification of hepatitis B virus (HBV) surface antigen in 1967 [1]. An estimated 257 million people worldwide are living with HBV infection, and in 2015, hepatitis B resulted in 887,000 deaths, mostly from complications, including cirrhosis and hepatocellular carcinoma [2]. HBV infection exhibits an acute infection stage and a chronic liver infection, which is determined by the degree of virus replication and the intensity of host immune response. Infection in newborns, who are incapable of constructing a defective immune response, is more likely to result in chronic infection (ā90%) than that in adults (ā5%), in whom most primary infections are self-limited [3].
Safe and effective vaccine is available for all age groups to prevent HBV infection and development of chronic disease. More than 150 countries have vaccine immunization programs, with routine infant vaccination designated as a high priority in all countries [4]. However, about 4.5 million new infections occur each year, of which a quarter progress to liver disease, and vaccine coverage in developing countries with high endemicity is limited due to high cost and social hurdles [5]. When viral replication is observed in a patient, efficient therapy is needed to control the risk of disease progression. It is estimated that, if left untreated, approximately 15-25% of chronically infected individuals would develop liver cirrhosis and HCC after decades of infection [6]. And there is abundant evidence that antiviral therapy, in patients with long-term virological response, can improve liver histology by providing indirect support and possibly even reversing liver damage [7, 8].
Mathematical models for HBV infection have been developed to study the dynamical properties and control strategies analytically and numerically, with the forms of ordinary differential equations (ODE), partial differential equations (PDE), and discrete equations. Medley et al. [9] observed that the prevalence of HBV is largely determined by a feedback mechanism that relates the rate of transmission, average age at infection, and age-related probability of progressing to chronic infection. Several models were constructed to study the transmission dynamics of this disease and applied to simulate statistical data in some specific regions and countries [10, 11]. An SLICRV compartmental model was proposed to investigate the prevalence of HBV infection in China from 2003 to 2008 [12]. Assuming the immunity acquired from vaccination and recovery being both lifelong, an SEICR model was used for data fitting based on available HBV epidemic data in China [13]. Mann and Robers [14] divided the whole population into five age classes to model the epidemiology of hepatitis B in New Zealand. Considering the important consequences of age for HBV infection, PDE models were developed to capture the transmission characteristics, evaluating the long-term effectiveness of vaccination programmes and analyzing the transmission dynamics [15, 16]. A model with age structure was formulated to study the possible effects of variable infectivity, partitioning the entire infectious period into two stages corresponding to acute and chronic infection, respectively [17]. By a discrete dynamic model, Liang et al. [18] evaluated the independent impact of newborn vaccination on reducing HBV prevalence in China.
In this study, we propose a model for HBV with three infectious compartments in the human population, i.e., two infected classes and a treatment class, which is different from the models in [19ā22] and mentioned above. The model we consider is formulated based on HBV transmission among population, compared with the in-host model of delay differential equations in [21]. Due to different motivation in [22], we focus on the transmission of HBV monoinfection and take no consideration of its superinfection with HDV, because the prevalence and new infected cases of HDV are obviously much lower than those of HBV in China [23]. Model application is done to simulate HBV data released by the National Health Commission of China from 2004 to 2016, extending the numerical results in [12, 13, 24] by more epidemiological data. Furthermore, in order to examine the effect of newborn vaccination and treatment on its prevention, the corresponding two parameters are regarded as functions of time, resulting in an optimality system, different from the control strategies with vaccination for the susceptibles in [19, 20]. The existence of the solution to the optimality system is discussed by classical optimal theory, converting the problem of minimizing the objective functional to minimizing the Hamiltonian with respect to the control. Based on the parameter values confirmed in model application, simulations of optimal states, new infected cases, and prevalence of chronic carriers are compared to show the different results with optimal control and with current interventions from initial year 2017 till 2030. Optimal control is also plotted to provide information of implementation as time evolves.
This paper is organized as follows. In Section 2, the descriptions of the parameters and variables, and model formulation are done. In Section 3, we study the stability of steady states and uniform persistence of the disease briefly and conduct sensitivity analysis for input parameters and outcome variable. Model application is done to simulate the epidemiological data in China in Section 4. In Section 5, the optimality system is analyzed theoretically and numerically. We conclude in Section 6 with discussions.
2. Formulation of the Model
When modeling HBV transmission, we divide the whole population into five classes, i.e., susceptible individuals, acute infections, chronic carriers, treated patients, and immunized individuals, which are denoted by , , , , and , respectively. In our model, vertical transmission, i.e., transmission from mother to child, is incorporated, since it has great impact on HBV prevalence, especially in highly endemic areas. The assumptions about vertical transmission can be referred to in [13]. For the newborns, they are assumed to be immunized successfully at rate or unsuccessfully at rate , with being the birth rate and being the successfully immunized proportion. When unsuccessfully immunized, the population of is assumed to enter into the chronic carriers, and the rest stays in the susceptible state, with denoting the probability of children developing to chronic state born to carrier mothers. The virus can also be horizontally transmitted by patients in the classes of , , and , and their different transmission rates are supposed to be , , and , respectively, with the assumption of and , reflecting the fact that the acute patients are the most infectious in the three infected states.
The acute infections are assumed to progress at rate , splitting in to chronic class and to the immunized class due to different immune response in host. In chronic class, the carriers can clear virus and become immunized at rate . We take consideration of treatment for both acute and chronic infection and assume that there are rates of and in the acute and chronic class, respectively, receiving treatment and then moving to the treated class. Once being treated, some patients can be cured and then move to immunized class at rate . However, some people may experience treatment failure and move back to the chronic state at rate . In each class, natural death occurs at rate .
There are some other assumptions in modeling listed in the following.
(1) The recruitment into susceptible population is simplified as new birth, and the birth rate and death rate are assumed to be identical; i.e., .
(2) The infected people who experience treatment failure move back only to chronic state due to the short period of acute infection.
(3) We take no consideration of HBV induced death in view of the efforts of improving liver histology by treatment.
(4) The immunity acquired by transient infection and vaccination is lifelong.
A flow diagram is presented in Figure 1 to show the transitions from one state to another. Then our mathematical model for HBV transmission is given by the following system of ordinary differential equations:Notice that the variable does not appear in other equations in (1); thus the equation of can be ignored, and the reduced system, which is studied in the following, has the same dynamical behavior as the original system. For system (1) without equation , its dynamical behavior remains in the following positively invariant subset of :
For convenience, let For epidemiological model, the basic reproductive number is a crucial threshold which determines the dynamics of the model. It gives the expected number of secondary cases produced in a completely susceptible population by a typical infective patient and can be computed by the approach in [25]. Then the basic reproductive number for model (1) can be defined as In the above expression of , there are three separate components, which give the average number of new infections produced by three infectious classes, respectively.
3. Dynamic Properties
In the following, we mainly analyze the stabilities of two equilibria, i.e., disease-free equilibrium and endemic equilibrium. Uniform persistence of system (1) is explored theoretically. In addition to the evaluation of parameters related to vaccination and treatment, the sensitivity of parameter variations on the basic reproductive number is also illustrated numerically.
There always exists a disease-free equilibrium for system (1), and there is a unique endemic equilibrium if and only if , where
3.1. Stability and Persistence
Theorem 1. For system (1),
(i) when , there is only a disease-free equilibrium and it is globally asymptotically stable;
(ii) when , there exists a unique endemic equilibrium and it is locally stable.
Proof. (i) The global stability of can be proved by a Lyapunov function and LaSalleās invariance principle. Specifically, set with Computing the derivative of along system (1) leads to , which implies that is the largest invariant set for all solutions of system (1) if ; thus we get the desired result of .
(ii) At endemic equilibrium , we have the corresponding characteristic equation aswhere and here By the equations in (1), we have the fact that and then it is easy to get , and we can calculate to find that and . Therefore, by Hurwitz criterion, the characteristic equation (7) has only roots with negative real part, which implies that is locally stable.
Next we will focus on the uniform persistence for model (1) when . To do this, we begin by letting and here denotes the semiflow from to defined by system (1).
Theorem 2. System (1) is uniformly persistent with respect to when ; i.e., for , there is
Proof. We conduct this proof by verifying the following three results as in [26].
(a) is the maximal compact invariant set in .
(b) When , for any , the solution of model (1) with initial value in satisfies .
(c) Conditions in Theorem 1.3.1 in [27] are satisfied.
Firstly, for (a), it is suffice to show that the solution of model (1) from remains in it if and only if . Contradictorily, without loss of generality, if , we will examine the values of , and in time interval . In fact, from the equations of , , and in (1), we find and is obvious, which implies that the result is valid.
Secondly, we deduce (b) still by contradiction. Assume that, for any , there exists a time such that , , and when . Then for the equationit is known that the equilibrium of (14) is globally stable and it converges to as . By comparison principle, we have for small enough.
To get the desired result, compared with model (1), we utilize the following system about variables , and :For this system, the local stability of its solution is determined by the eigenvalues of Jacobian matrix of (15) with , whereDenote the set of eigenvalues of matrix as , with being the real part of . It can be easily shown that when ; thus for small enough. Noticing that is a quasi-positive matrix with nonnegative off-diagonal elements, by Corollary 4.3.2 in [28], is one of eigenvalues of ; i.e., . Then there exists a corresponding vector such that , leading to as . Note the fact that , , and , a contradiction.
Finally, two conditions (C1) and (C2) in Theorem 1.3.1 [27] are to be verified. From system (1), we know that and are positively invariant, and the semiflow defined by (1) is compact and point dissipative, which ensures the existence of a global attractor for . In addition, the fact that is the maximal compact invariant set in is indicated by (a), and we can choose as the Morse decomposition of ; then , implying has the property of isolation. Furthermore, it has been shown in (b) that any solution from always remains in it, so . Therefore, we have the uniform persistence of system (1) when with .
3.2. Numerical Simulations
In this section, numerical simulations of model (1) are provided to illustrate the variability of the basic reproductive number due to some interested parameters. Here parameters related to newborn vaccination and treatment are highlighted, especially the rate of birth that successfully immunized and rate of transfer from treated to immunized , which are evaluated analytically and numerically in the following to assess their possible impact on the disease transmission. To consider the variation in the basic reproductive number with these two parameters, we have the following derivatives:
Obviously, is negatively related to , implying that increasing the successful vaccination at birth will lead to the declining of , which is desired in the prevention of the disease. From the second derivative, we can see that is also inversely related to and will decrease with the increase of it. In fact, denotes the successful rate of treatment for HBV patients, so with the coverage and efficiency of newborn hepatitis B vaccination being improved and in combination with the promising future in the therapy of preventing complications from HBV infection, the prevalence of this disease may be reduced to some extent in highly HBV-endemic areas.
Based on the Latin Hypercube Sampling scheme (LHS), we calculate the partial rank correlation coefficients (PRCC, denoted as P) to explore the sensitivity of parameter variations on the basic reproductive number due to the uncertainty in estimating input parameters. In the calculation, eight parameters are chosen as the input variables and the value of as the output variable. For each input parameter, it is sampled 2000 times and assumed to be normal in its distribution. By the values of P calculated for the eight parameters, we can identify the significance of any parameter contribution to the variability of , which are presented in Table 1 as well as P values (it is denoted as zero when being smaller than 0.0001). The P values can indicate the correlation between interested parameters and the basic reproductive number . Specifically, the absolute values of P reveal the correlation being important, moderate, or not significantly different from zero, and plus or minus sign means the influence is positive or negative, respectively. From Figure 2, we can see that is most sensitive to parameter with P=0.8544, followed by (P=-0.2277). Compared with them, the remaining six parameters contribute not so much to , with P values only under 0.2. Furthermore, and have positive impact, while the parameters , , , and have negative impact. Thus we know that large and are beneficial to reduce . In other words, newborn hepatitis B vaccination and effective treatment are considered as the most important measures to control the disease.
Contour plots are presented in Figure 3, which illustrate the variation of induced by the pairs of parameter and , , and , respectively, since these parameters are relatively more sensitive compared to the other ones as shown by sensitivity analysis in Figure 2, coinciding with the control measures of birth vaccination and treatment. From the left plot of and , we know that for fixed , when is increased to a certain value, then can be smaller than one, which means that it is possible to control the disease by measures reflected by and . At the same time, for some , may also be reduced to less than 1 with increasing of . For the other pair of and , their joint effect on is estimated in the right plot, showing that remaining and to a certain range can also make the value of less than one.
4. Model Application
HBV infection has long been a major health problem in China, and it is still among the areas of high endemicity of HBV, largely due to the great source of chronic carriers. To control this disease, great efforts have been made to implement immunization programmes and to reduce the rate of chronic infection as much as possible. In 2002, routine vaccination for newborns was fully integrated into the National Children Immunization Program, so that the coverage of newborn vaccination was increased greatly from then on. In 2006, the third national serosurvey showed that the HBsAg prevalence in the overall people had dropped from 9.75% to 7.18%, which revealed that there are about 93 million people all over the country being infected by this virus. In the past few years, the incidence of hepatitis B has lowered slowly according to the National Health Commission of China. The reported number of newly infected was approximately 1,180,000 in 2009; after that the number reduced to 934,215 in 2015, although it went up a little in 2016. The annually new reported HBV cases are shown in Table 2.
Based on these epidemiological data in Table 2, model (1) is employed to fit the cases from 2004 to 2016 by the least squares method (LSM), with the aim of minimizing the difference between model simulation and yearly new reported data. In simulation, the total population is and then the variables are normalized. The values of most parameters and variables are determined according to the demographic data in [29], the typical features of HBV transmission [15, 30], and LSM for better accuracy, and the others are taken by estimation, such as the successful rate of newborn vaccination and initial value of chronic population . Note that the birth rate and death rate are assumed to be identical in (H1), so we take the average of them in simulation ( by [29]). All these values are listed in Table 3. With these parameters, the basic reproductive number is calculated as .
The simulation result and the annually new reported cases are shown in Figure 4. From Figure 4 we can see the comparison between the simulation curve and the real data. The prediction of HBV trends until 2020 is also presented under the same parameters, which shows that the total new infected cases would decrease gradually year by year if the current interventions remain. Obviously, there exists some error in the simulation result for the possible reasons of imperfect model, imprecise parameters, and so on, but it still gives us some hint about the incidence of HBV in China in the next few years.
In order to examine the impact of vaccination and treatment on HBV prevalence in China, we focus on four sensitive parameters, i.e., , , , and , to investigate the variation in prevalence of chronic carriers due to the change of them, which are shown in Figure 5. From the plot with , we can see the different results with the increasing of , showing that the higher the successful rate of newborn vaccination, the higher the contribution to reduction of HBV prevalence; thus the full 3-dose and timely birth dose vaccination coverage should be improved as much as possible. The other three parameters, , , and , are related to treatment. If is taken as a baseline, then the prevalence can be reduced under 5% before 2030 when increases fourfold to 0.24, as shown in the plot for varying . The influence of and on chronic carriers is illustrated in the second row of Figure 5, which gives the possible reduction in prevalence when is decreased from 0.2323 to 0.1323 and is increased from 0.0936 to 0.1936, respectively.
5. Optimal Control Problem
In this section, by classical optimal theory in [31, 32], combined strategies of newborn vaccination and treatment are used to deal with optimal control problem of hepatitis B infection. The desired goal is to minimize the number of susceptible and infected population and the cost of implementing the control strategy, which involves tradeoff between these two competing factors.
5.1. Analysis of Optimal Control
Based on model (1), the two control measures, i.e., newborn vaccination and treatment, enter the system as control functions, denoted by and , respectively, appearing as function of time in the optimality system, different from the corresponding form of constant in model (1). The analysis of this nonautonomous system can make contributions to control management. Specifically, in order to achieve the optimal control, what are the percentages of the newborns that should be vaccinated and the infected population that should be treated as time evolves in a given period? The state variables (, , , , and ) satisfy the following differential equations which depend on the control variables:
Similarly, the equation of is omitted during the analysis as before. As the control functions and are changed, the solution to the differential system will change. We assume that and are measurable and have imposed bounds of 0 and 1; i.e., . Our optimal control problem consists of finding , , and the associated state variables to minimize the given objective functional at a time interval [], which is given byNote, in the above objective functional, in addition to an integral of the state and control variables, there is a term of the chronic population at the final time, which is called a payoff term or referred to as the salvage term. In the integrand, the coefficients and are weight constants used to balance the size of concerned population and the importance of control cost. The form of quadratic costs is introduced to model decreasing returns of treatment and vaccination efforts and to make the problem more tractable. By adjusting the control functions, the optimality problem can be manipulated to find the optimal control , which is used to achieve the goal of minimizing the objective functional; i.e., with being the control space. The existence of optimal control can be established since the integrand in (19) and the right hand side of system (18), denoted by , are continuously differentiable functions and concave in both and .
Theorem 3. Let the control function be measurable on with value in . Then there exists an optimal control minimizing the objective functional of (19) with where is the corresponding solution of (18).
Proof. In fact, we can generate the necessary conditions of this optimal problem by forming the Hamiltonian , which is defined as follows: with being the associated adjoint variable, which can be written in terms of the Hamiltonian and the transversality conditions at the final time give The information of differential equations (18) is attached onto the minimization of the objective functional by adjoint variable . Therefore, the problem of finding that minimizes the objective functional subjected to (18), is converted to minimizing the Hamiltonian with respect to the control. Then by Pontryagin Maximum Principle, the optimality condition leads towhich can be solved in terms of state and adjoint variables to give for as For the optimal control , which requires considering the constrains on the control and the sign of , we have It is similar for , so that the optimal control can be written in a compact way asTo find the optimal state, we can deal with system (18) when substituting into the equations. The proof is completed.
5.2. Numerical Simulations
For the optimality system consisting of state equations (18) and adjoint equations, numerical solutions and optimal control measures are illustrated in this subsection to show the consequence of introducing optimal control in our model. The differential equations of state system are solved by standard ODE solver when giving an initial condition for the state and control variables. Simulations of adjoint variable begin from the transversality condition at the final time according to its equations in the optimality system by a fourth-order Runge-Kutta scheme. Then the characterization of optimal control is updated by using the new values of state and adjoint variables.
In the simulations, we mainly focus on the results of implementing optimal control to reveal the effectiveness of vaccination and antiviral treatment. Based on the parameter values used in Figure 4, we take the case in China, for example, where the basic reproductive number is the same as in Section 4, i.e., , and the initial year is taken as 2017, which means from that year the optimal control is implemented in the model to suppress the prevalence of the disease. We take the weight coefficients , implying that the minimization of population number in different classes has the same importance, with being varied. Note that the value of is taken according to the form of integrand in (19); otherwise the nonlinear term of optimal control has no impact on the objective functional if is too small compared with the population of state variables.
The simulations of state variables are compared in Figure 6, which show the results with optimal control, versus the corresponding population where current interventions are implemented, plotted together from 2017 to 2030, with the time unit being year. As depicted in Figure 6, the numbers of the susceptible and infected individuals, which are components in the objective function (19), are held much lower when the optimal control is used than that when current measures are used. Notice that the population of , , and increases slightly at the end of the time interval when the level of control implementation decreases, but it is still much less overall than that when current interventions is used. These comparisons begin from the initial year 2017, and the predictions of the subsequent years till 2030 are simulated with the same parameters as those in model application in Section 4, where the weight coefficients , .
For comparison, the simulation curves of new infected cases and the prevalence of chronic carriers are shown in Figure 7, which are two quantities that measure the severity of the disease. As desired, the new infected cases decrease from its initial amount, and the prevalence of chronic HBV behaves in the same way. Specially, the prevalence can be forced to less than 3%; otherwise, under present intervention, it will remain above 5% until 2030. However, it is plotted that the lowered control strength also allows for a slight increase in the number of new infected cases at the end of the interval, which is corresponding to the optimal control in Figure 8 when less significance is placed on control strategy in the last few years.
With the aim of minimizing the objective functional subject, optimal control, including newborn vaccination and treatment , must be manipulated properly and economically as time evolves, which is plotted in Figure 8. In the simulations, we examine the optimal control for two sets of , i.e., , and , . It is shown that the optimal control will be affected by these two weight parameters, which balance the relative importance of the two terms of and . The two controls have almost the same performance when they are equal (right plot), remaining at their upper bound for the majority of the interval and then decreasing rapidly to zero eventually. If we use a lower weight parameter for and a higher parameter for , meaning increasing the relative importance between them, the controls change as expected (left plot). is forced to stay at its upper bound for almost all the interval and reduced to zero at the end, but stays at zero for a much longer time compared with .
6. Discussion
The model studied in this paper is a different version of the existing ones investigated in [12, 13, 17ā19], with treated population as a distinguishing compartment. Treatment of HBV infection in an attempt to improve liver histology and then reduce the risk of progression is a desirable control, in addition to newborn vaccination designated as a high priority to prevent the disease. Thus, a deterministic mathematical model incorporating control measures of newborn vaccination and treatment is constructed to analyze the dynamic behaviors and its application of real situation in China. By calculating the basic reproductive number , we conduct the theoretical analysis of the model which is determined by and provide numerical simulations to explore the sensitivity of model parameters on it. Meanwhile, using least squares method, we use the model to simulate the yearly new reported data released in China from 2004 to 2016, comparing the results of data fitting by the model with the epidemiological data and providing the prediction in the next few years. Due to more available data, the application in this paper is comparatively more accurate than that in [13].
Optimal control techniques are of great use in developing optimal strategy for complex biological situations. Since control combination of newborn vaccination and treatment is highlighted in the disease prevention, it is studied by classical optimal theory when these two parameters appear as functions of time in an optimality system. The aim of the control problem is to minimize the given objective functional by adjusting the control at a time interval, which confirmed that there exists an optimal control by forming the Hamiltonian. Simulations suggest that the level of state population can be held much lower if the optimal control is implemented completely at the given interval than that with current interventions used in model application.
In addition to newborn vaccination and treatment, other measures may also exert considerable impact on HBV control, such as the prevention of horizontal transmission by vaccination in other age populations and safe injection practice. Meanwhile, the actual progression of HBV is complex, and the interventions in preventing this disease are changed with time according to their performances and the administration policies. For the sake of simplicity, we make some assumptions for the model formulation. However, to have more accurate description of HBV characteristics and better simulation result of practical application, more factors which influence the transmission of it should be included in the model, and it is valuable to discuss hepatitis B infection in some high risk group or at smaller spatial scales, which will be explored in our future research.
Data Availability
All the data used to support the findings of this study are included in our manuscript and can be accessed freely from the references.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (nos. 11501443 and 11571275).