Printer Friendly

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

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 ([approximately equal to]90%) than that in adults ([approximately equal to]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 S, 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[omega] or unsuccessfully at rate b(1 - [omega]), with b being the birth rate and [omega] being the successfully immunized proportion. When unsuccessfully immunized, the population of b(1-[omega])vC is assumed to enter into the chronic carriers, and the rest b(1 - [omega])(1 - vC) stays in the susceptible state, with v 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 [[beta].sub.1], [[beta].sub.2], and [[beta].sub.3], respectively, with the assumption of [[beta].sub.1] > [[beta].sub.2] and [[beta].sub.1] > [[beta].sub.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 [[gamma].sub.1], splitting in [eta][[gamma].sub.1] to chronic class and (1 - [eta])[[gamma].sub.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 [[gamma].sub.2]. We take consideration of treatment for both acute and chronic infection and assume that there are rates of [[theta].sub.1] and [[theta].sub.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 [sigma]. However, some people may experience treatment failure and move back to the chronic state at rate [epsilon]. In each class, natural death occurs at rate [mu].

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 = [mu].

(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:

[mathematical expression not reproducible] (1)

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 [R.sup.4]:

[OMEGA] = {(S, I, C, T) | S, I,C, T [greater than or equal to] 0, S + I + C + T [less than or equal to] b (1 - [omega])/[mu]}. (2)

For convenience, let

[K.sub.1] = [mu]+ [[gamma].sub.1] + [[theta].sub.1], [K.sub.2] = [mu]+ [[gamma].sub.2] + [[theta].sub.2], [K.sub.3] = [mu] + [sigma] + [epsilon]. (3)

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

[mathematical expression not reproducible]. (4)

In the above expression of [R.sub.0], 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 [R.sub.0] is also illustrated numerically.

There always exists a disease-free equilibrium [P.sub.0] = ([S.sub.0], 0, 0, 0) for system (1), and there is a unique endemic equilibrium [P.sup.*] = ([S.sup.*], [I.sup.*], [C.sup.*], [T.sup.*]) if and only if [R.sub.0] > 1, where

[mathematical expression not reproducible]. (5)

3.1. Stability and Persistence

Theorem 1. For system (1),

(i) when [R.sub.0] < 1, there is only a disease-free equilibrium [P.sub.0] and it is globally asymptotically stable;

(ii) when [R.sub.0] > 1, there exists a unique endemic equilibrium [P.sup.*] and it is locally stable.

Proof. (i) The global stability of [P.sub.0] can be proved by a Lyapunov function and LaSalle's invariance principle. Specifically, set [V.sub.1](I, C, T) = [A.sub.1]I + [A.sub.2]C + [A.sub.3]T with

[mathematical expression not reproducible]. (6)

Computing the derivative of V1 along system (1) leads to [V'.sub.1] [less than or equal to] [K.sub.1]([R.sub.0] -1)I, which implies that [P.sub.0] is the largest invariant set for all solutions of system (1) if [R.sub.0] < 1; thus we get the desired result of [P.sub.0].

(ii) At endemic equilibrium [P.sup.*], we have the corresponding characteristic equation as

[[lambda].sup.4] + [a.sub.1][[lambda].sup.3] + [a.sub.2][[lambda].sup.2] + [a.sub.3][lambda] + [a.sub.4] = 0, (7)


[mathematical expression not reproducible], (8)

and here

[mathematical expression not reproducible]. (9)

By the equations in (1), we have the fact that

[mathematical expression not reproducible], (10)

and then it is easy to get [a.sub.i] > 0 (i = 1, 2, 3, 4), and we can calculate to find that [a.sub.1][a.sub.2]-[a.sub.3] > 0 and [a.sub.3]([a.sub.1][a.sub.2]-[a.sub.3])- [a.sup.2.sub.1][a.sub.4] > 0. Therefore, by Hurwitz criterion, the characteristic equation (7) has only roots with negative real part, which implies that [P.sup.*] is locally stable.

Next we will focus on the uniform persistence for model (1) when [R.sub.0] > 1. To do this, we begin by letting

[mathematical expression not reproducible], (11)

and here [[PHI].sub.t] denotes the semiflow from X to X defined by system (1).

Theorem 2. System (1) is uniformly persistent with respect to ([X.sub.0], [partial derivative] [X.sub.0]) when [R.sub.0] > 1; i.e., for [delta] > 0, there is

[mathematical expression not reproducible]. (12)

Proof. We conduct this proof by verifying the following three results as in [26].

(a) [M.sub.[partial derivative]] = {(S, 0, 0, 0) | S [greater than or equal to] 0} is the maximal compact invariant set in [partial derivative][X.sub.0].

(b) When [R.sub.0] > 1, for any [epsilon] > 0, the solution of model (1) with initial value in [X.sub.0] satisfies [lim.sub.t[right arrow][infinity]] sup max {I(t), C(t), T(t)} > [epsilon].

(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.sub.[partial derivative]] 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

[mathematical expression not reproducible], (13)

and {(S, 0, 0, 0) | S [greater than or equal to] 0} [subset] [M.sub.[partial derivative]] is obvious, which implies that the result is valid.

Secondly, we deduce (b) still by contradiction. Assume that, for any [epsilon] > 0, there exists a time [t.sub.1] such that I(t) < [epsilon], C(t) < [epsilon], and T(t) < [epsilon] when t > [t.sub.1]. Then for the equation

[??]' = b(1 - [omega]) (1 - v[epsilon]) - ([[beta].sub.1][epsilon] + [[beta].sub.2][epsilon] + [[beta].sub.3][epsilon]) [??] - [mu][??], (14)

is known that the equilibrium of (14) is globally stable and it converges to S0 as [epsilon] [right arrow] 0. By comparison principle, we have S(t) [greater than or equal to] [??](t) [greater than or equal to] [S.sub.0] - 2[[epsilon].sub.1] for [[epsilon].sub.1] small enough.

To get the desired result, compared with model (1), we utilize the following system about variables [??], [??] and [??]:

[mathematical expression not reproducible]. (15)

For this system, the local stability of its solution is determined by the eigenvalues [lambda] of Jacobian matrix J of (15) with J = [J.sub.1] - 2[[epsilon].sub.1][J.sub.2], where

[mathematical expression not reproducible]. (16)

Denote the set of eigenvalues of matrix J as [sigma](J), s(J) = max{R[lambda]|[lambda][member of][sigma](J)} with R[lambda] being the real part of [lambda]. It can be easily shown that s([J.sub.1]) > 0 when [R.sub.0] > 1; thus s(J) > 0 for [[epsilon].sub.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) [member of] [sigma](J). Then there exists a corresponding vector [xi] > 0 such that J[xi] = s(J)[xi], leading to [mathematical expression not reproducible], 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 [X.sub.0] are positively invariant, and the semiflow [[PHI].sub.t] defined by (1) is compact and point dissipative, which ensures the existence of a global attractor for [[PHI].sub.t]. In addition, the fact that [M.sub.[partial derivative]] is the maximal compact invariant set in [partial derivative][X.sub.0] is indicated by (a), and we can choose {[P.sub.0]} as the Morse decomposition of [M.sub.[partial derivative]]; then [mathematical expression not reproducible], implying {[P.sub.0]} has the property of isolation. Furthermore, it has been shown in (b) that any solution from [X.sub.0] always remains in it, so [W.sup.s]([P.sub.0]) [intersection] [X.sub.0] = 0. Therefore, we have the uniform persistence of system (1) when [R.sub.0] > 1 with L = [X.sub.0].

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 [omega] and rate of transfer from treated to immunized [sigma], 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 [R.sub.0] with these two parameters, we have the following derivatives:

[mathematical expression not reproducible]. (17)

Obviously, [R.sub.0] is negatively related to [omega], implying that increasing the successful vaccination at birth will lead to the declining of [R.sub.0], which is desired in the prevention of the disease. From the second derivative, we can see that [R.sub.0] is also inversely related to [sigma] and [R.sub.0] will decrease with the increase of it. In fact, [sigma] 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 [R.sub.0] 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 [R.sub.0], 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 [R.sub.0]. 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 [R.sub.0] is most sensitive to parameter [omega] with P=0.8544, followed by [sigma] (P=-0.2277). Compared with them, the remaining six parameters contribute not so much to [R.sub.0], with P values only under 0.2. Furthermore, [[beta].sub.i] (i = 1, 2, 3) and [epsilon] have positive impact, while the parameters [omega], [[theta].sub.1], [[theta].sub.2], and [sigma] have negative impact. Thus we know that large [omega] and [sigma] are beneficial to reduce [R.sub.0]. 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 [R.sub.0] induced by the pairs of parameter [omega] and [sigma], [epsilon], and [[theta].sub.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 [omega] and [sigma], we know that for fixed [sigma], when [omega] is increased to a certain value, then [R.sub.0] can be smaller than one, which means that it is possible to control the disease by measures reflected by [omega] and [sigma]. At the same time, for some [omega], [R.sub.0] may also be reduced to less than 1 with increasing of [sigma]. For the other pair of [[theta].sub.2] and [epsilon], their joint effect on [R.sub.0] is estimated in the right plot, showing that remaining [[theta].sub.2] and [epsilon] to a certain range can also make the value of [R.sub.0] 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 1.3 x 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 [omega] 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, [mu] = 0.0072 by [29]).All these values are listed in Table 3. With these parameters, the basic reproductive number is calculated as [R.sub.0] = 1.0764.

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., [omega], [sigma], [epsilon], and [[theta].sub.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 [omega] = 0.65, 0.75, 0.85, we can see the different results with the increasing of [omega], 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, [sigma], [epsilon], and [[theta].sub.2], are related to treatment. If [sigma] = 0.06 is taken as a baseline, then the prevalence can be reduced under 5% before 2030 when [sigma] increases fourfold to 0.24, as shown in the plot for varying [sigma]. The influence of [epsilon] and [[theta].sub.2] on chronic carriers is illustrated in the second row of Figure 5, which gives the possible reduction in prevalence when [epsilon] is decreased from 0.2323 to 0.1323 and [[theta].sub.2] 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 [omega](t) and [sigma](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:

[mathematical expression not reproducible], (18)

Similarly, the equation of R(t) is omitted during the analysis as before. As the control functions [omega](t) and [sigma](t) are changed, the solution to the differential system will change. We assume that [omega](t) and [sigma](t) are measurable and have imposed bounds of 0 and 1; i.e., 0 [less than or equal to] [omega](t), [sigma](t) [less than or equal to] 1. Our optimal control problem consists of finding [omega](t), [sigma](t), and the associated state variables to minimize the given objective functional at a time interval [[t.sub.0], [t.sub.f]], which is given by

[mathematical expression not reproducible]. (19)

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 [A.sub.i] (i = 1, 2, 3, 4) and [B.sub.i] (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.sup.*] = ([[omega].sup.*], [[sigma].sup.*]), which is used to achieve the goal of minimizing the objective functional; i.e.,

[mathematical expression not reproducible], (20)

with [OMEGA] = {([omega](t), [sigma](t)) [member of] [L.sup.1]([t.sub.0], [t.sub.f]) | 0 [less than or equal to] [omega](t),[sigma](t) [less than or equal to] 1} being the control space. The existence of optimal control [u.sup.*] can be established since the integrand in (19) and the right hand side of system (18), denoted by f = ([f.sub.1], [f.sub.2], [f.sub.3], [f.sub.4]), are continuously differentiable functions and concave in both (S, I, C, T) and u = ([omega], [sigma]).

Theorem 3. Let the control function u = ([omega], [sigma]) [member of] [OMEGA] be measurable on [t.sub.0] [less than or equal to] t [less than or equal to] [t.sub.f] with value in [0, 1].Then there exists an optimal control [u.sup.*] = ([[omega].sup.*], [[sigma].sup.*]) minimizing the objective functional J([omega], [sigma]) of (19) with

[[omega].sup.*] = max {min {[b/[B.sub.1]] [[[lambda].sub.1] (1 - vC) + [[lambda].sub.3]vC], 1}, 0}, [[sigma].sup.*] = max {min {[[[lambda].sub.4]T/[B.sub.2]], 1}, 0}, (21)

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 Hamiltonian H, which is defined as follows:

H = [A.sub.1]S (t) + [A.sub.2]I (t) + [A.sub.3]C(t) + [A.sub.4]T (t) + [[B.sub.1]/2] [[omega].sup.2] (t) + [[B.sub.2]/2] [[sigma].sup.2] (t) + [4.summation over (i=1)] [[lambda].sub.i][f.sub.i], (22)

With [lambda] = ([[lambda].sub.1], [[lambda].sub.2], [[lambda].sub.3], [[lambda].sub.4]) being the associated adjoint variable, which can be written in terms of the Hamiltonian

[mathematical expression not reproducible], (23)

and the transversality conditions at the final time give

[[lambda].sub.1] ([t.sub.f]) = 0, [[lambda].sub.2] ([t.sub.f]) = 0, [[lambda].sub.3] ([t.sub.f]) = 1, [[lambda].sub.4] ([t.sub.f]) = 0. (24)

The information of differential equations (18) is attached onto the minimization of the objective functional by adjoint variable [lambda]. Therefore, the problem of finding [u.sup.*] 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

[partial derivative]H/[partial derivative][omega] = [B.sub.1][omega] - [[lambda].sub.1]b (1 - vC) - [[lambda].sub.3]bvC = 0, [partial derivative]H/[partial derivative][sigma] = [B.sub.2][sigma] - [[lambda].sub.4]T = 0, (25)

which can be solved in terms of state and adjoint variables to give [bar.u] = ([bar.[omega]], [bar.[sigma]]) for u = ([omega],[sigma]) as

[bar.[omega]] = [b/[B.sub.1]] [[[lambda].sub.1] (1 - vC) + [[lambda].sub.3]vC], [bar.[sigma]] = [[lambda].sub.4]T/[B.sub.2]. (26)

For the optimal control [u.sup.*], which requires considering the constrains on the control and the sign of [partial derivative]H/[partial derivative]u, we have

[mathematical expression not reproducible]. (27)

It is similar for [[sigma].sup.*], so that the optimal control can be written in a compact way as

[[omega].sup.*] = max {min {[bar.[omega]], 1}, 0}, [[sigma].sup.*] = max {min {[bar.[sigma]], 1}, 0} . (28)

To find the optimal state, we can deal with system (18) when substituting [u.sup.*] 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., [R.sub.0] = 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 [A.sub.i] = 1 (i = 1,2,3,4), implying that the minimization of population number in different classes has the same importance, with [B.sub.i] (i = 1, 2) being varied. Note that the value of [B.sub.i] 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 [B.sub.i] 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 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 [B.sub.1] = [10.sup.7], [B.sub.2] = [10.sup.12].

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 [omega](t) and treatment [sigma](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 [B.sub.i], i.e., [B.sub.1] = [10.sup.7], [B.sub.2] = [10.sup.12] and [B.sub.1] = [10.sup.10], [B.sub.2] = [10.sup.10]. 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 [omega] and [sigma]. 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 [omega] and a higher parameter for [sigma], meaning increasing the relative importance between them, the controls change as expected (left plot). [omega] is forced to stay at its upper bound for almost all the interval and reduced to zero at the end, but [sigma] stays at zero for a much longer time compared with [omega].

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 [R.sub.0], we conduct the theoretical analysis of the model which is determined by [R.sub.0] 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.


This work was supported by the National Natural Science Foundation of China (nos. 11501443 and 11571275).


[1] A. Suk-Fong Lok, "Hepatitis B: 50 years after the discovery of Australia antigen," Journal of Viral Hepatitis, vol. 23, no. 1, pp. 5-14, 2016.

[2] World Health Organization, "Hepatitis B. World Health Organization Fact Sheet No. 204," 2017, mediacentre/factsheets/fs204/en/index.html.

[3] J.H. Hoofnagle, E. Doo, T. J. Liang, R. Fleischer, and A. S. F. Lok, "Management of hepatitis B: summary of a clinical research workshop," Hepatology, vol. 45, no. 4, pp. 1056-1075, 2007.

[4] E. Mast and J. Ward, "Hepatitis B Vaccine," in Vaccines, S. Plotkin, W. Orenstein, and P. Offit, Eds., pp. 205-242, WB Saunders Company, Philadelphia, PA, USA, 5th edition, 2008.

[5] E. Franco, B. Bagnato, M. G. Marino, C. Meleleo, L. Serino, and L. Zaratti, "Hepatitis B: epidemiology and prevention in developing countries," World Journal of Hepatology, vol. 4, no. 3, pp. 74-80, 2012.

[6] T. M. Block, S. Rawat, and C. L. Brosgart, "Chronic hepatitis B: a wave of new therapies on the horizon," Antiviral Research, vol. 121, pp. 69-81, 2015.

[7] A. S. Lok, "Does antiviral therapy for hepatitis B and C prevent hepatocellular carcinoma?" Journal of Gastroenterology and Hepatology, vol. 26, no. 2, pp. 221-227, 2011.

[8] P. Bedossa, K. Patel, and L. Castera, "Histologic and noninvasive estimates of liver fibrosis," Clinical Liver Disease, vol. 6, no. 1, pp. 5-8, 2015.

[9] G. F. Medley, N. A. Lindop, W. J. Edmunds, and D. J. Nokes, "Hepatitis-B virus endemicity: Heterogeneity, catastrophic dynamics and control," Nature Medicine, vol. 7, no. 5, pp. 619-624, 2001.

[10] S. Thornley, C. Bullen, and M. Roberts, "Hepatitis B in a high prevalence New Zealand population: a mathematical model applied to infection control policy," Journal of Theoretical Biology, vol. 254, no. 3, pp. 599-603, 2008.

[11] T. Zhang, K. Wang, and X. Zhang, "Modelling and analyzing the transmission dynamics of HBV epidemic in Xinjiang, China," Plos One, 2015.

[12] L. Zou, W. Zhang, and S. Ruan, "Modeling the transmission dynamics and control of hepatitis B virus in China," Journal of Theoretical Biology, vol. 262, no. 2, pp. 330-338, 2010.

[13] S. Zhang and Y. Zhou, "The analysis and application of an HBV model," Applied Mathematical Modelling: Simulation and Computation for Engineering and Environmental Systems, vol. 36, no. 3, pp. 1302-1312, 2012.

[14] J. Mann and M. Roberts, "Modelling the epidemiology of hepatitis B in New Zealand," Journal of Theoretical Biology, vol. 269, pp. 266-272, 2011.

[15] S. Zhao, Z. Xu, and Y. Lu, "A mathematical model of hepatitis B virus transmission and its application for vaccination strategy in China," International Journal of Epidemiology, vol. 29, no. 4, pp. 744-752, 2000.

[16] L. Zou, S. Ruan, and W. Zhang, "An age-structured model for the transmission dynamics of hepatitis B," SIAM Journal on Applied Mathematics, vol. 70, no. 8, pp. 3121-3139, 2010.

[17] S. Zhang and X. Xu, "A mathematical model for hepatitis B with infection-age structure," Discrete and Continuous Dynamical Systems--Series B, vol. 21, no. 4, pp. 1329-1346, 2016.

[18] P. Liang, J. Zu, J. Yin et al., "The independent impact of newborn hepatitis B vaccination on reducing HBV prevalence in China, 1992-2006: a mathematical model analysis," Journal of Theoretical Biology, vol. 386, pp. 115-121, 2015.

[19] A. Kamyad, R. Akbari, A. Heydari et al., "Mathematical modeling of transmission dynamics and optimal control of vaccination and treatment for hepatitis B virus," Computational and Mathematical Methods in Medicine, vol. 2014, Article ID 475451, 15 pages, 2014.

[20] T. Khan, G. Zaman, and M. Ikhlaq Chohan, "The transmission dynamic and optimal control of acute and chronic hepatitis B," Journal of Biological Dynamics, vol. 11, no. 1, pp. 172-189, 2017.

[21] J. E. Forde, S. M. Ciupe, A. Cintron-Arias, and S. Lenhart, "Optimal control of drug therapy in a hepatitis B model," Applied Sciences, vol. 6, no. 8, article 219, 2016.

[22] A. Goyal and J.M. Murray, "Roadmap to control HBV and HDV epidemics in China," Journal of Theoretical Biology, vol. 423, pp. 41-52, 2017.

[23] National Health and Family Planning Commission of the People's Republic of China, "National report of notifiable diseases, 2004-2016,"

[24] J. Pang, J. A. Cui, and X. Zhou, "Dynamical behavior of a hepatitis B virus transmission model with vaccination," Journal of Theoretical Biology, vol. 265, no. 4, pp. 572-578, 2010.

[25] P. van den Driessche and J. Watmough, "Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission," Mathematical Biosciences, vol. 180, pp. 29-48, 2002.

[26] S. Zhang and X. Xu, "Dynamic analysis and optimal control for a model of hepatitis C with treatment," Communications in Nonlinear Science and Numerical Simulation, vol. 46, pp. 14-25, 2017.

[27] X. Q. Zhao, Dynamical Systems in Population Biology, Springer, New York, NY, USA, 2003.

[28] H. L. Smith, Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems, vol. 41 of Mathematical Surveys and Monographs, American Mathematical Society, Providence, RI, USA, 1995.

[29] National Bureau of Statistics of China, "China Statistical Yearbook 2016, Birth rate, Death rate and Natural growth rate of population,"

[30] W. J. Edmunds, G. F. Medley, and D. J. Nokes, "The transmission dynamics and control of hepatitis B virus in the Gambia," Statistics in Medicine, vol. 15, no. 20, pp. 2215-2233, 1996.

[31] L. Pontryagin, The Mathematical Theory of Optimal Processed, Taylor & Francis, London, UK, 1987.

[32] S. M. Lenhart and J. T. Workman, Optimal Control Applied to Biological Models, Taylor & Francis, Boca Raton, FL, USA, 2007.

Jing Zhang and Suxia Zhang (iD)

School of Science, Xi'an University of Technology, Xi'an 710048, China

Correspondence should be addressed to Suxia Zhang;

Received 26 April 2018; Accepted 13 August 2018; Published 2 September 2018

Academic Editor: Yong Zhou

Caption: Figure 1: Flow chart of the epidemiological classes for the model.

Caption: Figure 2: Sensitivity analysis for eight input parameters and outcome variable (the basic reproductive number).

Caption: Figure 3: Contour plots of the basic reproductive number [R.sub.0] versus the parameter pairs ([sigma], [omega]) (left) and ([[theta].sub.2], [epsilon]) (right).

Caption: Figure 4: Data fitting for the annually reported cases (pentagram) from 2004 to 2016 by model (1) (solid curve).

Caption: Figure 5: Proportion of chronic carriers against time with the same parameters and initial values as in Figure 4 (solid), compared with the curves of varying [omega], [sigma], [epsilon], and [[theta].sub.2] in each plot, respectively (dashed).

Caption: Figure 6: Optimal states for (18) plotted as time evolves. Dashed lines present the simulations with optimal control [u.sup.*] and solid lines present that with current interventions in China, with the same parameters as in Figure 4 and [R.sub.0] = 1.0764.

Caption: Figure 7: Simulations of new infected cases (left) and prevalence of chronic carriers (right) under optimal control and current interventions from 2017 to 2030.

Caption: Figure 8: Optimal control plots with different [B.sub.i] when [R.sub.0] = 1.0764. Left: [B.sub.1] = [10.sup.7], [B.sub.2] = [10.sup.12]. Right: [B.sub.1] = [10.sup.10], [B.sub.2] = [10.sup.10].
Table 1: Sensitivity values of the basic reproductive number on the

Input parameters     The basic reproductive
                        number [R.sub.0]

                     P value    P value

[[beta].sup.1]        0.0566     0.0113
[[beta].sup.2]        0.1557       0
[[beta].sup.3]        0.1881       0
[epsilon]             0.1956       0
[[theta].sub.1]      -0.0318     0.1554
[[theta].sub.2]      -0.1925       0
[sigma]              -0.2277       0
[omega]               0.8544       0

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,268

Table 3: List of parameters and initial values used in simulation.

Parameter                   Value             Source

b                    0.0096 [year.sup.-1]      [29]
v                     0.11 [year.sup.-1]     [15, 30]
[[gamma].sub.2]      0.025 [year.sup.-1]     [15, 30]
[[beta].sub.1]             12.4929              LSM
[[beta].sub.3]              1.5372              LSM
[[theta].sub.1]      0.0578 [year.sup.-1]       LSM
[sigma]               0.06 [year.sup.-1]        LSM


S(0)                        0.0745              LSM
T(0)                    .706423936e-3           LSM

Parameter                   Value             Source

[mu]                 0.0096 [year.sup.-1]      [29]
[[gamma].sub.1]        4 [year.sup.-1]       [15, 30]
[eta]                0.885 [year.sup.-1]     [15, 30]
[[beta].sub.2]              0.4002              LSM
[epsilon]            0.2323 [year.sup.-1]       LSM
[[theta].sub.2]      0.0936 [year.sup.-1]       LSM
[omega]               0.65 [year.sup.-1]      Assumed


I(0)                    .761409263e-3           LSM
C(0)                        0.076             Assumed
COPYRIGHT 2018 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2018 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Title Annotation:Research Article
Author:Zhang, Jing; Zhang, Suxia
Publication:Discrete Dynamics in Nature and Society
Geographic Code:9CHIN
Date:Jan 1, 2018
Previous Article:Bargaining Power Choices with Moral Hazard in a Supply Chain.
Next Article:Proactive Hedging European Call Option Pricing with Linear Position Strategy.

Terms of use | Privacy policy | Copyright © 2019 Farlex, Inc. | Feedback | For webmasters