# Application and Optimal Control for an HBV Model with Vaccination and Treatment

Authors: Jing Zhang; Suxia Zhang
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## 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 numberR0. 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.


---

## Body

## 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 Section2, 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 byS, I, C, T, and R, 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 bω or unsuccessfully at rate b(1-ω), with b being the birth rate and ω being the successfully immunized proportion. When unsuccessfully immunized, the population of b(1-ω)υC is assumed to enter into the chronic carriers, and the rest b(1-ω)(1-υC) 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 I, C, and T, and their different transmission rates are supposed to be β1, β2, and β3, respectively, with the assumption of β1>β2 and β1>β3, 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γ1, splitting in ηγ1 to chronic class and (1-η)γ1 to the immunized class due to different immune response in host. In chronic class, the carriers can clear virus and become immunized at rate γ2. We take consideration of treatment for both acute and chronic infection and assume that there are rates of θ1 and θ2 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.,b=μ.(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 Figure1 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:(1)S′t=b1-ω1-υC-β1I+β2C+β3TS-μS,I′t=β1I+β2C+β3TS-μ+γ1+θ1I,C′t=ηγ1I+b1-ωυC+εT-μ+γ2+θ2C,T′t=θ1I+θ2C-μ+σ+εT,R′t=bω+1-ηγ1I+γ2C+σT-μR.Notice that the variable R does not appear in other equations in (1); thus the equation of R 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 R, its dynamical behavior remains in the following positively invariant subset of R4: (2)Ω=S,I,C,T∣S,I,C,T≥0,S+I+C+T≤b1-ωμ.Figure 1
Flow chart of the epidemiological classes for the model.For convenience, let(3)K1=μ+γ1+θ1,K2=μ+γ2+θ2,K3=μ+σ+ε.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 (4)R0=b1-ωμβ1K1+β2ηγ1K3+θ1εK1K3K2-b1-ωυ-θ2ε+β3ηγ1θ2+θ1K2-b1-ωυK1K3K2-b1-ωυ-θ2ε.In the above expression of R0, 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 R0 is also illustrated numerically.There always exists a disease-free equilibriumP0=(S0,0,0,0) for system (1), and there is a unique endemic equilibrium P∗=(S∗,I∗,C∗,T∗) if and only if R0>1, where (5)S0=b1-ωμ,S∗=S0R0,I∗=b1-ω1-1R0K3K2-b1-ωυ-θ2εb1-ωυηγ1K3+θ1ε+K1K3K2-b1-ωυ-θ2ε,C∗=b1-ω1-1R0ηγ1K3+θ1εb1-ωυηγ1K3+θ1ε+K1K3K2-b1-ωυ-θ2ε,T∗=b1-ω1-1R0θ1K2-b1-ωυ+θ2ηγ1b1-ωυηγ1K3+θ1ε+K1K3K2-b1-ωυ-θ2ε.


### 3.1. Stability and Persistence

Theorem 1.
For system (1),
(i) whenR0<1, there is only a disease-free equilibrium P0 and it is globally asymptotically stable;
(ii) whenR0>1, there exists a unique endemic equilibrium P∗ and it is locally stable.Proof.
(i) The global stability ofP0 can be proved by a Lyapunov function and LaSalle’s invariance principle. Specifically, set V1(I,C,T)=A1I+A2C+A3T with (6)A1=1,A2=S0β2K3+β3θ2K3K2-b1-ωυ-θ2ε,A3=S0β3K2-b1-ωυ+β2εK3K2-b1-ωυ-θ2ε.Computing the derivative of V1 along system (1) leads to V1′≤K1(R0-1)I, which implies that P0 is the largest invariant set for all solutions of system (1) if R0<1; thus we get the desired result of P0.
(ii) At endemic equilibriumP∗, we have the corresponding characteristic equation as(7)λ4+a1λ3+a2λ2+a3λ+a4=0,where (8)a1=K3+K4+K5+K7,a2=K3K4+K5+K7+K4K5+K4K7+K5K7-θ2ε-β12I∗+θ1β3+ηγ1β2S∗,a3=K4K7K3+K5+K3K5K4+K7+θ1β1I∗β3S∗-θ2εK4+K5-K4+K7θ1β3+β2εθ1+β12I∗K3+K7+β3ηγ1θ2+ηγ1β2K3+K4S∗,a4=K3K4K5K7+β1I∗ηγ1θ2β3S∗+ηγ1β3S∗β1I∗K6+K6β1I∗εθ1+θ1β1I∗β3S∗K7-β2εθ1S∗K4-β12S∗I∗K3K7-θ2εK4K5-β3S∗ηγ1θ2K4-ηγ1β2S∗K3K4-θ1β3S∗K4K7-θ1εβ12I∗S∗,and here (9)K4=μ+β1I∗+β2C∗+β3T∗,K5=K1-β1S∗,K6=b1-ωυ+β2S∗,K7=K2-b1-ωυ.By the equations in (1), we have the fact that (10)β12I∗+θ1β3+ηγ1β2S∗+θ2ε<K3K7,K4+K7θ1β3+β2εθ1+β12I∗K3+K7+β3ηγ1θ2+ηγ1β2K3+K4S∗<K3K5K4+K7+K3K4K7,and then it is easy to get ai>0(i=1,2,3,4), and we can calculate to find that a1a2-a3>0 and a3(a1a2-a3)-a12a4>0. Therefore, by Hurwitz criterion, the characteristic equation (7) has only roots with negative real part, which implies that P∗ is locally stable.Next we will focus on the uniform persistence for model (1) when R0>1. To do this, we begin by letting (11)X=S,I,C,T∣S≥0,I≥0,C≥0,T≥0,X0=S,I,C,T∈X∣S≥0,I>0,C>0,T>0,∂X0=X∖X0,M∂=S0,I0,C0,T0∈∂X0∣ΦtS0,I0,C0,T0∈∂X0,t≥0,and here Φt denotes the semiflow from X to X defined by system (1).Theorem 2.
System (1) is uniformly persistent with respect to (X0,∂X0) when R0>1; i.e., for δ>0, there is (12)limt→∞⁡infIt>δ,limt→∞infCt>δ,limt→∞infTt>δ.Proof.
We conduct this proof by verifying the following three results as in [26].
(a)M∂={(S,0,0,0)∣S≥0} is the maximal compact invariant set in ∂X0.
(b) WhenR0>1, for any ϵ>0, the solution of model (1) with initial value in X0 satisfies limt→∞supmax{I(t),C(t),T(t)}>ϵ.
(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 M∂ remains in it if and only if I(0)=C(0)=T(0)=0. Contradictorily, without loss of generality, if I(0)>0, we will examine the values of I(t),C(t), and T(t) in time interval [0,1]. In fact, from the equations of I, C, and T in (1), we find (13)It≥I0e-K1≜M1>0,Ct≥ηγ1M1K2-b1-ωυ1-e-K2-b1-ωυ≜M2>0,Tt>θ1M1+θ2M2K31-e-K3≜M3>0,and {(S,0,0,0)∣S≥0}⊂M∂ is obvious, which implies that the result is valid.
Secondly, we deduce (b) still by contradiction. Assume that, for anyϵ>0, there exists a time t1 such that I(t)<ϵ, C(t)<ϵ, and T(t)<ϵ when t>t1. Then for the equation(14)S~′=b1-ω1-υϵ-β1ϵ+β2ϵ+β3ϵS~-μS~,it is known that the equilibrium of (14) is globally stable and it converges to S0 as ϵ→0. By comparison principle, we have S(t)≥S~(t)≥S0-2ϵ1 for ϵ1 small enough.
To get the desired result, compared with model (1), we utilize the following system about variables I~,C~, and T~:(15)I~′=β1I~+β2C~+β3T~S0-2ϵ1-K1I~,C~′=ηγ1I~+b1-ωυC~+εT~-K2C~,T~′=θ1I~+θ2C~-K3T~.For this system, the local stability of its solution is determined by the eigenvalues λ of Jacobian matrix J of (15) with J=J1-2ϵ1J2, where(16)J1=β1S0-K1β2S0β3S0ηγ1b1-ωυ-K2εθ1θ2-K3,J2=β1β2β3000000.Denote the set of eigenvalues of matrix J as σ(J), s(J)=max{Rλ∣λ∈σ(J)} with Rλ being the real part of λ. It can be easily shown that s(J1)>0 when R0>1; thus s(J)>0 for ϵ1 small enough. Noticing that J is a quasi-positive matrix with nonnegative off-diagonal elements, by Corollary 4.3.2 in [28], s(J) is one of eigenvalues of J; i.e., s(J)∈σ(J). Then there exists a corresponding vector ξ>0 such that Jξ=s(J)ξ, leading to (I~(t),C~(t),T~(t))→∞ as t→∞. Note the fact that I(t)≥I~(t), C(t)≥C~(t), and T(t)≥T~(t), a contradiction.
Finally, two conditions (C1) and (C2) in Theorem 1.3.1 [27] are to be verified. From system (1), we know that X and X0 are positively invariant, and the semiflow Φt defined by (1) is compact and point dissipative, which ensures the existence of a global attractor for Φt. In addition, the fact that M∂ is the maximal compact invariant set in ∂X0 is indicated by (a), and we can choose {P0} as the Morse decomposition of M∂; then ⋃x∈M∂ω(x)={P0}, implying {P0} has the property of isolation. Furthermore, it has been shown in (b) that any solution from X0 always remains in it, so Ws(P0)∩X0=∅. Therefore, we have the uniform persistence of system (1) when R0>1 with L=X0.


### 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 R0 with these two parameters, we have the following derivatives: (17)∂R0∂ω=-bυK1K3ηγ1+θ1εβ2K3+β3θ2K12K3K-θ2ε2<0,∂R0∂σ=-K1ηγ1θ2+Kθ1β2ε+β3KK12K3K-θ2ε2<0.Obviously,R0 is negatively related to ω, implying that increasing the successful vaccination at birth will lead to the declining of R0, which is desired in the prevention of the disease. From the second derivative, we can see that R0 is also inversely related to σ and R0 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 ofR0 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 R0, 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 R0. 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 R0 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 R0, with P values only under 0.2. Furthermore, βi(i=1,2,3) and ε have positive impact, while the parameters ω, θ1, θ2, and σ have negative impact. Thus we know that large ω and σ are beneficial to reduce R0. In other words, newborn hepatitis B vaccination and effective treatment are considered as the most important measures to control the disease.Table 1
Sensitivity values of the basic reproductive number on the parameters.
Input parameters The basic reproductive numberR0 P value P value β1 0.0566 0.0113 β2 0.1557 0 β3 0.1881 0 ε 0.1956 0 θ1 -0.0318 0.1554 θ2 -0.1925 0 σ -0.2277 0 ω 0.8544 0Figure 2
Sensitivity analysis for eight input parameters and outcome variable (the basic reproductive number).Contour plots are presented in Figure3, which illustrate the variation of R0 induced by the pairs of parameter ω and σ, ε, and θ2, 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 R0 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 ω, R0 may also be reduced to less than 1 with increasing of σ. For the other pair of θ2 and ε, their joint effect on R0 is estimated in the right plot, showing that remaining θ2 and ε to a certain range can also make the value of R0 less than one.Figure 3
Contour plots of the basic reproductive numberR0 versus the parameter pairs (σ,ω) (left) and (θ2,ε) (right).


## 3.1. Stability and Persistence

Theorem 1.
For system (1),
(i) whenR0<1, there is only a disease-free equilibrium P0 and it is globally asymptotically stable;
(ii) whenR0>1, there exists a unique endemic equilibrium P∗ and it is locally stable.Proof.
(i) The global stability ofP0 can be proved by a Lyapunov function and LaSalle’s invariance principle. Specifically, set V1(I,C,T)=A1I+A2C+A3T with (6)A1=1,A2=S0β2K3+β3θ2K3K2-b1-ωυ-θ2ε,A3=S0β3K2-b1-ωυ+β2εK3K2-b1-ωυ-θ2ε.Computing the derivative of V1 along system (1) leads to V1′≤K1(R0-1)I, which implies that P0 is the largest invariant set for all solutions of system (1) if R0<1; thus we get the desired result of P0.
(ii) At endemic equilibriumP∗, we have the corresponding characteristic equation as(7)λ4+a1λ3+a2λ2+a3λ+a4=0,where (8)a1=K3+K4+K5+K7,a2=K3K4+K5+K7+K4K5+K4K7+K5K7-θ2ε-β12I∗+θ1β3+ηγ1β2S∗,a3=K4K7K3+K5+K3K5K4+K7+θ1β1I∗β3S∗-θ2εK4+K5-K4+K7θ1β3+β2εθ1+β12I∗K3+K7+β3ηγ1θ2+ηγ1β2K3+K4S∗,a4=K3K4K5K7+β1I∗ηγ1θ2β3S∗+ηγ1β3S∗β1I∗K6+K6β1I∗εθ1+θ1β1I∗β3S∗K7-β2εθ1S∗K4-β12S∗I∗K3K7-θ2εK4K5-β3S∗ηγ1θ2K4-ηγ1β2S∗K3K4-θ1β3S∗K4K7-θ1εβ12I∗S∗,and here (9)K4=μ+β1I∗+β2C∗+β3T∗,K5=K1-β1S∗,K6=b1-ωυ+β2S∗,K7=K2-b1-ωυ.By the equations in (1), we have the fact that (10)β12I∗+θ1β3+ηγ1β2S∗+θ2ε<K3K7,K4+K7θ1β3+β2εθ1+β12I∗K3+K7+β3ηγ1θ2+ηγ1β2K3+K4S∗<K3K5K4+K7+K3K4K7,and then it is easy to get ai>0(i=1,2,3,4), and we can calculate to find that a1a2-a3>0 and a3(a1a2-a3)-a12a4>0. Therefore, by Hurwitz criterion, the characteristic equation (7) has only roots with negative real part, which implies that P∗ is locally stable.Next we will focus on the uniform persistence for model (1) when R0>1. To do this, we begin by letting (11)X=S,I,C,T∣S≥0,I≥0,C≥0,T≥0,X0=S,I,C,T∈X∣S≥0,I>0,C>0,T>0,∂X0=X∖X0,M∂=S0,I0,C0,T0∈∂X0∣ΦtS0,I0,C0,T0∈∂X0,t≥0,and here Φt denotes the semiflow from X to X defined by system (1).Theorem 2.
System (1) is uniformly persistent with respect to (X0,∂X0) when R0>1; i.e., for δ>0, there is (12)limt→∞⁡infIt>δ,limt→∞infCt>δ,limt→∞infTt>δ.Proof.
We conduct this proof by verifying the following three results as in [26].
(a)M∂={(S,0,0,0)∣S≥0} is the maximal compact invariant set in ∂X0.
(b) WhenR0>1, for any ϵ>0, the solution of model (1) with initial value in X0 satisfies limt→∞supmax{I(t),C(t),T(t)}>ϵ.
(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 M∂ remains in it if and only if I(0)=C(0)=T(0)=0. Contradictorily, without loss of generality, if I(0)>0, we will examine the values of I(t),C(t), and T(t) in time interval [0,1]. In fact, from the equations of I, C, and T in (1), we find (13)It≥I0e-K1≜M1>0,Ct≥ηγ1M1K2-b1-ωυ1-e-K2-b1-ωυ≜M2>0,Tt>θ1M1+θ2M2K31-e-K3≜M3>0,and {(S,0,0,0)∣S≥0}⊂M∂ is obvious, which implies that the result is valid.
Secondly, we deduce (b) still by contradiction. Assume that, for anyϵ>0, there exists a time t1 such that I(t)<ϵ, C(t)<ϵ, and T(t)<ϵ when t>t1. Then for the equation(14)S~′=b1-ω1-υϵ-β1ϵ+β2ϵ+β3ϵS~-μS~,it is known that the equilibrium of (14) is globally stable and it converges to S0 as ϵ→0. By comparison principle, we have S(t)≥S~(t)≥S0-2ϵ1 for ϵ1 small enough.
To get the desired result, compared with model (1), we utilize the following system about variables I~,C~, and T~:(15)I~′=β1I~+β2C~+β3T~S0-2ϵ1-K1I~,C~′=ηγ1I~+b1-ωυC~+εT~-K2C~,T~′=θ1I~+θ2C~-K3T~.For this system, the local stability of its solution is determined by the eigenvalues λ of Jacobian matrix J of (15) with J=J1-2ϵ1J2, where(16)J1=β1S0-K1β2S0β3S0ηγ1b1-ωυ-K2εθ1θ2-K3,J2=β1β2β3000000.Denote the set of eigenvalues of matrix J as σ(J), s(J)=max{Rλ∣λ∈σ(J)} with Rλ being the real part of λ. It can be easily shown that s(J1)>0 when R0>1; thus s(J)>0 for ϵ1 small enough. Noticing that J is a quasi-positive matrix with nonnegative off-diagonal elements, by Corollary 4.3.2 in [28], s(J) is one of eigenvalues of J; i.e., s(J)∈σ(J). Then there exists a corresponding vector ξ>0 such that Jξ=s(J)ξ, leading to (I~(t),C~(t),T~(t))→∞ as t→∞. Note the fact that I(t)≥I~(t), C(t)≥C~(t), and T(t)≥T~(t), a contradiction.
Finally, two conditions (C1) and (C2) in Theorem 1.3.1 [27] are to be verified. From system (1), we know that X and X0 are positively invariant, and the semiflow Φt defined by (1) is compact and point dissipative, which ensures the existence of a global attractor for Φt. In addition, the fact that M∂ is the maximal compact invariant set in ∂X0 is indicated by (a), and we can choose {P0} as the Morse decomposition of M∂; then ⋃x∈M∂ω(x)={P0}, implying {P0} has the property of isolation. Furthermore, it has been shown in (b) that any solution from X0 always remains in it, so Ws(P0)∩X0=∅. Therefore, we have the uniform persistence of system (1) when R0>1 with L=X0.


## 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 R0 with these two parameters, we have the following derivatives: (17)∂R0∂ω=-bυK1K3ηγ1+θ1εβ2K3+β3θ2K12K3K-θ2ε2<0,∂R0∂σ=-K1ηγ1θ2+Kθ1β2ε+β3KK12K3K-θ2ε2<0.Obviously,R0 is negatively related to ω, implying that increasing the successful vaccination at birth will lead to the declining of R0, which is desired in the prevention of the disease. From the second derivative, we can see that R0 is also inversely related to σ and R0 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 ofR0 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 R0, 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 R0. 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 R0 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 R0, with P values only under 0.2. Furthermore, βi(i=1,2,3) and ε have positive impact, while the parameters ω, θ1, θ2, and σ have negative impact. Thus we know that large ω and σ are beneficial to reduce R0. In other words, newborn hepatitis B vaccination and effective treatment are considered as the most important measures to control the disease.Table 1
Sensitivity values of the basic reproductive number on the parameters.
Input parameters The basic reproductive numberR0 P value P value β1 0.0566 0.0113 β2 0.1557 0 β3 0.1881 0 ε 0.1956 0 θ1 -0.0318 0.1554 θ2 -0.1925 0 σ -0.2277 0 ω 0.8544 0Figure 2
Sensitivity analysis for eight input parameters and outcome variable (the basic reproductive number).Contour plots are presented in Figure3, which illustrate the variation of R0 induced by the pairs of parameter ω and σ, ε, and θ2, 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 R0 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 ω, R0 may also be reduced to less than 1 with increasing of σ. For the other pair of θ2 and ε, their joint effect on R0 is estimated in the right plot, showing that remaining θ2 and ε to a certain range can also make the value of R0 less than one.Figure 3
Contour plots of the basic reproductive numberR0 versus the parameter pairs (σ,ω) (left) and (θ2,ε) (right).


## 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 Table2.Table 2
Annually new reported HBV cases [23].
Year 2004 2005 2006 2007 2008  Cases 916,426 982,297 1,109,130 1,169,946 1,169,569  Year 2009 2010 2011 2012 2013  Cases 1,179,607 1,060,582 1,093,335 1,087,086 962,974  Year 2014 2015 2016  Cases 935,702 934,215 942,268Based on these epidemiological data in Table2, 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 1.3×109 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 C(0). Note that the birth rate and death rate are assumed to be identical in (H1), so we take the average of them in simulation (b=0.0121,μ=0.0072 by [29]). All these values are listed in Table 3. With these parameters, the basic reproductive number is calculated as R0=1.0764.Table 3
List of parameters and initial values used in simulation.
Parameter Value Source Parameter Value Source b 0.0096year-1 [29] μ 0.0096year-1 [29] υ 0.11year-1 [15, 30] γ 1 4year-1 [15, 30] γ2 0.025year-1 [15, 30] η 0.885year-1 [15, 30] β1 12.4929 LSM β 2 0.4002 LSM β3 1.5372 LSM ϵ 0.2323year-1 LSM θ1 0.0578year-1 LSM θ 2 0.0936year-1 LSM σ 0.06year-1 LSM ω 0.65year-1 Assumed  Variable   Variable  S(0) 0.0745 LSM I ( 0 ) .761409263e-3 LSM T(0) .706423936e-3 LSM C ( 0 ) 0.076 AssumedThe simulation result and the annually new reported cases are shown in Figure4. 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.Figure 4
Data fitting for the annually reported cases (pentagram) from 2004 to 2016 by model (1) (solid curve).In order to examine the impact of vaccination and treatment on HBV prevalence in China, we focus on four sensitive parameters, i.e.,ω, σ, ε, and θ2, 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 ω=0.65,0.75,0.85, 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 θ2, are related to treatment. If σ=0.06 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 θ2 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 θ2 is increased from 0.0936 to 0.1936, respectively.Figure 5
Proportion of chronic carriers against time with the same parameters and initial values as in Figure4 (solid), compared with the curves of varying ω, σ, ε, and θ2 in each plot, respectively (dashed).


## 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 ω(t) and σ(t), 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 (S(t), I(t), C(t), T(t), and R(t)) satisfy the following differential equations which depend on the control variables:(18)S′t=b1-ωt1-υC-β1I+β2C+β3TS-μS,I′t=β1I+β2C+β3TS-μ+γ1+θ1I,C′t=ηγ1I+b1-ωtυC+εT-μ+γ2+θ2C,T′t=θ1I+θ2C-μ+σt+εT,R′t=bωt+1-ηγ1I+γ2C+σtT-μR.Similarly, the equation ofR(t) is omitted during the analysis as before. As the control functions ω(t) and σ(t) are changed, the solution to the differential system will change. We assume that ω(t) and σ(t) are measurable and have imposed bounds of 0 and 1; i.e., 0≤ω(t),σ(t)≤1. Our optimal control problem consists of finding ω(t), σ(t), and the associated state variables to minimize the given objective functional at a time interval [t0,tf], which is given by(19)Jω,σ=Ctf+∫t0tfA1St+A2It+A3Ct+A4Tt+B12ω2t+B22σ2tdt.Note, 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 Ai(i=1,2,3,4) and Bi(i=1,2) 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 u∗=(ω∗,σ∗), which is used to achieve the goal of minimizing the objective functional; i.e., (20)Jω∗,σ∗=minΩ⁡Jωt,σt,with Ω={(ω(t),σ(t))∈L1(t0,tf)∣0≤ω(t),σ(t)≤1} being the control space. The existence of optimal control u∗ can be established since the integrand in (19) and the right hand side of system (18), denoted by f=(f1,f2,f3,f4), are continuously differentiable functions and concave in both (S,I,C,T) and u=(ω,σ).Theorem 3.
Let the control functionu=(ω,σ)∈Ω be measurable on t0≤t≤tf with value in [0,1]. Then there exists an optimal control u∗=(ω∗,σ∗) minimizing the objective functional J(ω,σ) of (19) with (21)ω∗=max⁡min⁡bB1λ11-υC+λ3υC,1,0,σ∗=max⁡min⁡λ4TB2,1,0,where (S,I,C,T) is the corresponding solution of (18).Proof.
In fact, we can generate the necessary conditions of this optimal problem by forming the HamiltonianH, which is defined as follows: (22)H=A1St+A2It+A3Ct+A4Tt+B12ω2t+B22σ2t+∑i=14λifi,with λ=(λ1,λ2,λ3,λ4) being the associated adjoint variable, which can be written in terms of the Hamiltonian (23)dλ1tdt=-∂H∂S=-A1+λ1μ+β1I+β2C+β3T-λ2β1I+β2C+β3T,dλ2tdt=-∂H∂I=-A2+λ1β1S+λ2μ+γ1+θ1-λ2β1S-λ3ηγ1-λ4θ1,dλ3tdt=-∂H∂C=-A3+λ1bυ1-ωt+λ1β2S-λ2β2S-λ3bυ1-ωt+λ3μ+γ2+θ2-λ4θ2,dλ4tdt=-∂H∂T=-A4+λ1β3S-λ2β3S-λ3ε+λ4μ+σt+ε,and the transversality conditions at the final time give (24)λ1tf=0,λ2tf=0,λ3tf=1,λ4tf=0.The information of differential equations (18) is attached onto the minimization of the objective functional by adjoint variable λ. Therefore, the problem of finding u∗ 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 to(25)∂H∂ω=B1ω-λ1b1-υC-λ3bυC=0,∂H∂σ=B2σ-λ4T=0,which can be solved in terms of state and adjoint variables to give u¯=(ω¯,σ¯) for u=(ω,σ) as (26)ω¯=bB1λ11-υC+λ3υC,σ¯=λ4TB2.For the optimal control u∗, which requires considering the constrains on the control and the sign of ∂H/∂u, we have (27)
ω
∗
=
0
,
i
f
∂
H
∂
ω
<
0
,
ω
¯
,
i
f
∂
H
∂
ω
=
0
,
1
,
i
f
∂
H
∂
ω
>
0
.It is similar for σ∗, so that the optimal control can be written in a compact way as(28)ω∗=maxminω¯,1,0,σ∗=max⁡min⁡σ¯,1,0.
To find the optimal state, we can deal with system (18) when substituting u∗ 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 Figure4, we take the case in China, for example, where the basic reproductive number is the same as in Section 4, i.e., R0=1.0764, 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 Ai=1(i=1,2,3,4), implying that the minimization of population number in different classes has the same importance, with Bi(i=1,2) being varied. Note that the value of Bi 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 Bi is too small compared with the population of state variables.The simulations of state variables are compared in Figure6, 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 S, I, and T 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 B1=107, B2=1012.Figure 6
Optimal states for (18) plotted as time evolves. Dashed lines present the simulations with optimal control u∗ and solid lines present that with current interventions in China, with the same parameters as in Figure 4 and R0=1.0764.For comparison, the simulation curves of new infected cases and the prevalence of chronic carriers are shown in Figure7, 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.Figure 7
Simulations of new infected cases (left) and prevalence of chronic carriers (right) under optimal control and current interventions from 2017 to 2030.Figure 8
Optimal control plots with differentBi when R0=1.0764. Left: B1=107, B2=1012. Right: B1=1010, B2=1010.With the aim of minimizing the objective functional subject, optimal control, including newborn vaccinationω(t) and treatment σ(t), 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 Bi, i.e., B1=107, B2=1012 and B1=1010, B2=1010. 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 ω.


## 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 ω(t) and σ(t), 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 (S(t), I(t), C(t), T(t), and R(t)) satisfy the following differential equations which depend on the control variables:(18)S′t=b1-ωt1-υC-β1I+β2C+β3TS-μS,I′t=β1I+β2C+β3TS-μ+γ1+θ1I,C′t=ηγ1I+b1-ωtυC+εT-μ+γ2+θ2C,T′t=θ1I+θ2C-μ+σt+εT,R′t=bωt+1-ηγ1I+γ2C+σtT-μR.Similarly, the equation ofR(t) is omitted during the analysis as before. As the control functions ω(t) and σ(t) are changed, the solution to the differential system will change. We assume that ω(t) and σ(t) are measurable and have imposed bounds of 0 and 1; i.e., 0≤ω(t),σ(t)≤1. Our optimal control problem consists of finding ω(t), σ(t), and the associated state variables to minimize the given objective functional at a time interval [t0,tf], which is given by(19)Jω,σ=Ctf+∫t0tfA1St+A2It+A3Ct+A4Tt+B12ω2t+B22σ2tdt.Note, 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 Ai(i=1,2,3,4) and Bi(i=1,2) 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 u∗=(ω∗,σ∗), which is used to achieve the goal of minimizing the objective functional; i.e., (20)Jω∗,σ∗=minΩ⁡Jωt,σt,with Ω={(ω(t),σ(t))∈L1(t0,tf)∣0≤ω(t),σ(t)≤1} being the control space. The existence of optimal control u∗ can be established since the integrand in (19) and the right hand side of system (18), denoted by f=(f1,f2,f3,f4), are continuously differentiable functions and concave in both (S,I,C,T) and u=(ω,σ).Theorem 3.
Let the control functionu=(ω,σ)∈Ω be measurable on t0≤t≤tf with value in [0,1]. Then there exists an optimal control u∗=(ω∗,σ∗) minimizing the objective functional J(ω,σ) of (19) with (21)ω∗=max⁡min⁡bB1λ11-υC+λ3υC,1,0,σ∗=max⁡min⁡λ4TB2,1,0,where (S,I,C,T) is the corresponding solution of (18).Proof.
In fact, we can generate the necessary conditions of this optimal problem by forming the HamiltonianH, which is defined as follows: (22)H=A1St+A2It+A3Ct+A4Tt+B12ω2t+B22σ2t+∑i=14λifi,with λ=(λ1,λ2,λ3,λ4) being the associated adjoint variable, which can be written in terms of the Hamiltonian (23)dλ1tdt=-∂H∂S=-A1+λ1μ+β1I+β2C+β3T-λ2β1I+β2C+β3T,dλ2tdt=-∂H∂I=-A2+λ1β1S+λ2μ+γ1+θ1-λ2β1S-λ3ηγ1-λ4θ1,dλ3tdt=-∂H∂C=-A3+λ1bυ1-ωt+λ1β2S-λ2β2S-λ3bυ1-ωt+λ3μ+γ2+θ2-λ4θ2,dλ4tdt=-∂H∂T=-A4+λ1β3S-λ2β3S-λ3ε+λ4μ+σt+ε,and the transversality conditions at the final time give (24)λ1tf=0,λ2tf=0,λ3tf=1,λ4tf=0.The information of differential equations (18) is attached onto the minimization of the objective functional by adjoint variable λ. Therefore, the problem of finding u∗ 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 to(25)∂H∂ω=B1ω-λ1b1-υC-λ3bυC=0,∂H∂σ=B2σ-λ4T=0,which can be solved in terms of state and adjoint variables to give u¯=(ω¯,σ¯) for u=(ω,σ) as (26)ω¯=bB1λ11-υC+λ3υC,σ¯=λ4TB2.For the optimal control u∗, which requires considering the constrains on the control and the sign of ∂H/∂u, we have (27)
ω
∗
=
0
,
i
f
∂
H
∂
ω
<
0
,
ω
¯
,
i
f
∂
H
∂
ω
=
0
,
1
,
i
f
∂
H
∂
ω
>
0
.It is similar for σ∗, so that the optimal control can be written in a compact way as(28)ω∗=maxminω¯,1,0,σ∗=max⁡min⁡σ¯,1,0.
To find the optimal state, we can deal with system (18) when substituting u∗ 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 Figure4, we take the case in China, for example, where the basic reproductive number is the same as in Section 4, i.e., R0=1.0764, 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 Ai=1(i=1,2,3,4), implying that the minimization of population number in different classes has the same importance, with Bi(i=1,2) being varied. Note that the value of Bi 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 Bi is too small compared with the population of state variables.The simulations of state variables are compared in Figure6, 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 S, I, and T 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 B1=107, B2=1012.Figure 6
Optimal states for (18) plotted as time evolves. Dashed lines present the simulations with optimal control u∗ and solid lines present that with current interventions in China, with the same parameters as in Figure 4 and R0=1.0764.For comparison, the simulation curves of new infected cases and the prevalence of chronic carriers are shown in Figure7, 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.Figure 7
Simulations of new infected cases (left) and prevalence of chronic carriers (right) under optimal control and current interventions from 2017 to 2030.Figure 8
Optimal control plots with differentBi when R0=1.0764. Left: B1=107, B2=1012. Right: B1=1010, B2=1010.With the aim of minimizing the objective functional subject, optimal control, including newborn vaccinationω(t) and treatment σ(t), 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 Bi, i.e., B1=107, B2=1012 and B1=1010, B2=1010. 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 R0, we conduct the theoretical analysis of the model which is determined by R0 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.



---

*Source: 2076983-2018-09-02.xml*