Modelling the dynamics of acute and chronic hepatitis B with optimal control

This article examines hepatitis B dynamics under distinct infection phases and multiple transmissions. We formulate the epidemic problem based on the characteristics of the disease. It is shown that the epidemiological model is mathematically and biologically meaningful of its well-posedness (positivity, boundedness, and biologically feasible region). The reproductive number is then calculated to find the equilibria and the stability analysis of the epidemic model is performed. A backward bifurcation is also investigated in the proposed epidemic problem. With the help of two control measures (treatment and vaccination), we develop control strategies to minimize the infected population (acute and chronic). To solve the proposed control problem, we utilize Pontryagin’s Maximum Principle. Some simulations are conducted to illustrate the investigation of the analytical work and the effect of control analysis.


Problem formulation
Based on disease transmission characteristics, an epidemic model is proposed to investigate hepatitis B virus transmission.There are four classes of host populations, symbolized by N(t): the susceptible class S(t), the acutely infected compartment A(t), the chronically ill class B(t), and the immunized/recovered class R(t).The following assumptions are made in our model.a 1 .Each parameter, as well as the variable used in the proposed epidemic problem, is nonnegative.a 2 .The vaccine for hepatitis B is very effective because it provides indefinite protection; therefore, the suscep- tible individuals, after being vaccinated successfully, lead to the recovered population.a 3 .Both the acutely infected and chronically infected individuals will cause the infection to be susceptible, and by successful interaction, the susceptible will lead to the acute class.a 4 .Natural death occurs in each model group, while death from disease only occurs in the chronic class.a 5 .The portion of newborns with maternal infection leads to B(t).
The schematic disease transmission process is demonstrated by Fig. 1.By grouping all of the above assumptions, a system of autonomous differential equations can be derived that describes the complete model  Investigation of the model (1) is subject to the initial sizes of compartments The parameter in (1) is the rate of newborns, v is the vaccination parameter, and η is the maternally infected rate.The symbol gamma denotes the reduced transmission rate, and µ 0 illustrates the proportion of natural death.Similarly, µ 1 is the portion of deaths that occurs due to the disease.We represent the contact parameter by α and the recovery rate from the acute class by γ 1 .Moreover, γ 2 symbolizes the recovery in chronically infected population, and β is the proportion of those who move from acute class to chronic one.First, we prove the well-posedness by illustrating the following results.Proof To determine that the model (1) possesses a unique solution, we follow the methodology given in 32 and define the vector field of the proposed model H : where The right-hand side of the Eq. ( 3) implies that the function H is continuous and therefore ensures the existence of solution (S, A, B, R) over an interval [0, ∞) .In addition, calculating the derivative of H with respect to the model state variables gives the Jacobian matrix as given by Since, DH is continuous over R 4 and thus H is locally Lipschitz continuous on (−∞, ∞) × R 4 , therefore, the model solution (S, A, B, R) is uniquely determined on the interval [0, ∞) .

Proposition 2.2 (Positivity of solution)
Let the solution to the problem (1)-(2) be symbolized by (S, A, B, R), whenever exists, then it is positive for all t greater than zero.
Proof Obviously, right-side functions in the system (1) satisfy the conditions of differentiability, implying the existence of a unique maximal solution for any associated Cauchy problem.Thus, the first equation of system (1) takes the form where ϕ = q 1 + ψ .The solution of (6) looks like Following the same steps, the model, second equation can be re-casted as which leads to Similarly, the last two equations of the epidemiological model (1) can be re-written as Integrating, we then obtain www.nature.com/scientificreports/Thus, from the above Eqs.( 7)-( 9), it could be observed that all the state variables of the proposed epidemiological model satisfying the initial conditions remain non-negative.

Proof Let
Taking the temporal derivatives of this equation and exploiting values from model (1), one may obtain Since by assumption µ 0 is a positive parameter and B ≥ 0 .Consequently we may write dN dt + µ 0 N ≤ � .Solution of this equation subject to the initial conditions (2) gives It is obvious that whenever t → ∞ , the last equation yields 0 < N ≤ � µ 0 .
Proposition 2.4 (Positively invariant set) Let N be the total population as given in (10), then the feasible region represented by is invariant positively and attracting for the proposed epidemiological model (1).
, then clearly Eq. ( 12) implies that N(t) ≤ � µ 0 .But on the other hand, on a contrary basis, if N(0) ≥ � µ 0 , then either the total population N(t) converge to � µ 0 as t increases without bound or the solution trajectories enter the feasible region within finite time, which implies that all the state variables initiated in R 4  + enter or converge � µ 0 eventually.

Basic reproductive number
There are two possible non-trivial equilibrium points of model (1), namely, the endemic and disease-free states.The disease-free state is represented by E 1 and is calculated as E 1 = (S 1 , 0, 0, R 1 ) , such that We use this state to calculate the so-called basic reproductive quantity, R 0 , which describes the average number of secondary infectious created by an index case, i.e. when an infective is presented into a susceptible population so the secondary infections are produced during its total infection age 33 .The reproductive number R 0 is then conveniently used to characterize the endemic equilibrium.Let us assume that J represents the linearized matrix of the system (1).Direct calculations show that the matrix J has the form with q 1 , q 2 and q 3 as given in (4).We follow Watmough and Driessche 34 to determine the threshold number of the epidemiological model that is under consideration.By assuming X = (A(t), B(t)) T , one can write from (1) that where F and V , in the above equation, are defined as The Jacobian or the linearized matrices of the above-defined F and V at the infection-free state (14) are respec- tively calculated as (9) B(t) ≥ B(0) exp −q 3 t and R(t) ≥ R(0) exp (−µ 0 t).
Vol www.nature.com/scientificreports/ The threshold quantity, R 0 is given by the largest eigenvalue of the matrix FV −1 .That is, R 0 = ρ(FV −1 ) .We deduced that Now we find the endemic-state of the model using (1) and ( 16), we obtain The characterization of the occurrence of the no-infection state ( 14) and of the disease-endemic state ( 17) is investigated in "Existence of backward bifurcation" section 4, while in "Stability analysis" section 5 we investigate their global analysis.

Existence of backward bifurcation
In epidemiological models, one of the necessary conditions to control the infection is R 0 < 1 .In contrast, this condition may not always be sufficient, owing to backward bifurcation, i.e., a stable endemic state co-exists with a stable infection-free state whenever R 0 < 1 .It is a common phenomenon in epidemiological models 35 .In this case, disease control depends upon the various sub-populations sizes of the epidemic problem.To investigate the existence of bifurcation, we suppose that at least one of the infected groups in system (1) is nonzero.In this situation, the solution of our proposed model (1) around steady state yields For B * � = 0 , we insert A * and S * in system (1) around steady state, and utilizing (4), we obtain the following equations It could be noted from the last relation that whenever the condition of R 0 < 1 holds, then b and c are non- negative.Also if R 0 > 1 then b < 0 .Clearly, a > 0 , so a positive solution of equation ( 18) exists, which depends on the signs of b, proving that the equilibrium continuously depends on the threshold quantity.Moreover equation (18) implies that For distinct ranges of the parameters, we state the underlying result.

Theorem 4.1 The considered epidemic problem (1) has:
(i) a unique endemic state in the biologically reasonable region (13) whenever b is negative and R 0 > 1; (ii) a unique endemic state in (13) if b = 0; (iii) two endemic steady states in (13) whenever b > 0.
It could be noted in the epidemiological models that one of the classical requirements for disease elimination is R 0 < 1 , while this is not sufficient 35 .Thus the condition R 0 < 1 is necessary for the control of hepatitis B but is not sufficient.Moreover, the backward bifurcation's presence in the model states that elimination of hepatitis B in case of R 0 < 1 depends on sub-populations of the model, and whenever R 0 = 1 , we have described the following.

Lemma 4.2
The existence of backward bifurcation for the model (1) depends on the value of R 0 and exists whenever R 0 = 1 , while experiences backward bifurcation in case condition (iii) of Theorem 4.1 holds.

Stability analysis
We now demonstrate the global analysis of the problem (1) at both the non-trivial equilibria.For the global properties around disease-free state ( 14), we use the classical Lyapunov function theory, while at disease-endemic state (17), we use the geometrical approach.
Proof Let h 1 > 0 , h 2 > 0 and h 3 > 0 be constants to be determined later.Consider a function of the form The temporal differentiation of this equation, along with values from system (1), gives Let us assume h 1 = h 2 = q 1 , and h 3 = αβ�γ q 3 .Further, from our previous calculations have S 1 = q 1 .Then from the last equation, we have Simplification of the above equation leads to Thus, when R 0 < 1 , we have 0 < R 01 < 1 and 0 < R 02 < 1 , then dF dt < 0 .Also, dF dt = 0 if S = S 1 , and B = A = 0 .Hence, the principle of LaSalle's 36,37 reveals that ( 14) is stable globally asymptotically.  1is globally stable.Endemic state of (1) is unstable whenever R 0 > 1 does not hold.
Proof Reducing model (1) by removing R(t) with the fact that it can be derived from the relation of the total populace, i.e., N = A + S + B + R , which implies that R = N − S − A − B .So without losing generality, it is enough to discuss the dynamics of the reduced model for the original.Thus, let J 2 is the Jacobian while J |2| 2 is the second order compound matrix of the proposed model (1), then where We consider a function in the form of P = P(S, A, B) = diag{S/A, S/A, S/A} , then taking its inverse and dif- ferentiating with respect to t, i.e.P f (χ) , we have Implying that , Let ℓ(B) is the Lozinski measure 38 with respect to the above equation, then it becomes where g i = �B ij + ℓ(B ii )� for i, j = 1, 2 and i = j , so then g 1 and g 2 are defined by and In the above last two equations �B 12 � = �η + γ αS and �B 21 � = max{β, 0} = β .Consequently g 1 and g 2 can be written as gives From here, we can write ℓ(B) ≤ Ṡ S − 2µ 0 .The application of integration for ℓ(B) in [0, t] with lim as t approaches ∞ gives is negative, which proves the conclusion.
According to the stability analysis of the model, the reproductive number is an important parameter controlling disease dynamics.In the following subsection, we will discuss how parameters affect reproductive numbers.
Influential parameters and its relative impact.For the purpose of creating an effective control mechanism for disease elimination, we conduct a sensitivity analysis of the model parameters in order to uncover the most influential parameters that highly affect the basic reproductive number.By using the formula given where φ is any epidemic parameter of the proposed epidemiological model associated with the reproductive number R 0 .By following the formula (20), we obtain the sensitivity indices of the proposed model parameters as Clearly, the normalized sensitivity indices of the parameters α , v, γ 1 and γ 2 show that α is positively correlated with the reproductive quantity, while v, γ 1 and γ 2 are negatively correlated.An increase in the values positively correlated results increase in the value of the reproductive number.But on, the value of the reproductive number will be reduced whenever the value of the negatively correlated parameters increases.Hence, control efforts should be formulated by taking suitable control measures to reduce the burden of hepatitis B virus transmission.

Formulation of control problem
The theory of optimization is a prominent tool and is used frequently in the dynamics of infectious epidemiology.
With the help of this, we can formulate strategies for the minimization of various kinds of infections.We follow the approach as presented by Gul et al. and others (see, 39-43 ) to set up a control problem for the reduction of the hepatitis B infection.We propose a control mechanism to minimize the burden of hepatitis B virus transmission by utilizing two control measures (vaccination and treatment) because the normalized sensitivity indices of v, γ 1 , and γ 2 have an inverse relationship with reproduction quantity.The aim is to vaccinate susceptible (control u 1 (t) ), and the treatment of infected (control u 2 (t) ).This is in contrast with 31 , where an inconsistent control problem was developed and consequently solved.Indeed, in 31 , the authors vaccinate all the infected individuals at the same rate as they vaccinate the susceptible.This seems not to be coherent with medical practice.In our case, the control system is obtained from (1) by placing two control measures already mentioned u 1 (t) and u 2 (t) with the description to vaccinate S as well as treatment of A and B. This implies that system (1) becomes a specific case of the proposed control problem whenever u 1 (t) ≡ v (vaccination with constant rate) and u 2 (t) ≡ 0 (when there is no treatment).Thus, the control problem takes the form subject to the problem In (22), w 1 and w 2 describe the relative weight constants of acute and chronic individuals, respectively.Also, the constants w 3 , w 4 ≥ 0 measure the associated costs of vaccination and treatment, respectively.It could be illustrated from ( 22) that the control problem has a clear purpose, namely, to reduce the ratio of A and B by implementing the control measures costs u 1 (t) and u 2 (t) .However, it is not our objective to reduce or increase the number of susceptible, as inconsistently proposed in 31 .Indeed, our goal is to determine the control functions (u * 1 , u * 2 ) like under the control problem (23).The set U (control set) is such that In addition, to discuss the existence analysis, the control problem can be expressed as where where Q = max {µ 0 , µ 0 + µ 1 } , known as the Lipschitz constants, and hence the function F(Y) is Lipschitz con- tinuous, which ensures the existence of an optimal solution to the proposed control problem.Now to determine the existence of optimal controls, one has to prove their existence.
Existence analysis.We now perform the existence analysis of the optimal controls for the proposed control problem as stated by the equations ( 22)- (25).For this, we prove that the set of control and associated states of the model are not empty as well as the control set is closed and convex.The state system of the control problem is linear in the control variables, while the integrand of the objective functional is convex over the control set U. Thus, regarding the existence analysis, we state the underlying result.
Theorem 6.1 For the control problem ( 22)-( 25), there exists a pair of optimal values u * = (u * 1 , u * 2 ) ∈ U , such that Proof In order to discuss the existence of optimal controls, we use Theorem 9.2.1, p.182, given by Lukes in 44 , and followed by various authors 45,46 .Note that both the state and control variable of the model ( 23) have a nonnegative value and bounded co-efficient, which implies that the control set and associated state variables are non-empty.Also, solutions are bounded, and hence the control set is closed and convex.The state system of the control problem is linear in the control variables, implying that the optimal system is bounded.In addition, the integrand of the objective functional ( 22) is convex over U, therefore there exists a constants ζ > 1 and positive numbers ξ 1 and ξ 2 such that which is enough for the optimal control pair existence.

Optimality conditions.
To examine an optimal solution to ( 22)-( 25), we determine the Hamiltonian on the basis of Lagrangian by assuming that x = (S, A, B, R) and u = (u 1 , u 2 ) be the vectors of state and control measure, respectively.Denoting the Lagrangian and Hamiltonian by L and H respectively, then one may write and In the above equations, we have taken = ( 1 , 2 , 3 , 4 ) and g = g 1 , g 2 , g 3 , g 4 where We now find the solution (optimal) to the control model and exploit the conventional Maximum Principle 47 : Let (x * , u * ) represents the optimal solution, then one can find a non-trivial function satisfying with the maximality condition and the transversality condition www.nature.com/scientificreports/satisfied.The next theorem follows from ( 26)- (31).We note that the optimality conditions given in 31 are inconsistent.It is clear that in the case of minimization and a Hamiltonian, together with the associated multiplier as well as Lagrangian +1 , the Pontryagin principle emphasizes a condition of minimality rather than a condition maximal- ity.Here, condition (30) is the maximality condition for the proposed minimization problem as the multiplier −1 of the Hamiltonian associated with L. Theorem 6.2 Assume S * , A * , B * and R * respectively denote the optimal state variables with the accompanying optimal measures (u * 1 , u * 2 ) for the problem ( 22)-( 25), then i (t) , i = 1, . . ., 4 exists i.e., the adjoint variables exist which satisfy The terminal (transversality) conditions associated are The optimal measures u * 1 and u * 2 are as and Proof System (32) is derived from the Pontryagin Maximum Principle i.e. from the 2nd relation of (??) with the Hamiltonian as described in ( 26)-( 28), while conditions (33) follow from the transversality condition (31).To derive u * 1 and u * 2 , we differentiate the Hamiltonian partially and solve ∂H ∂u 1 = 0 and ∂H ∂u 2 = 0 for control measures.Finally, with the help of the maximality condition (30), we derive (34)- (35).

Numerical simulations
To support our theoretical results, we present the numerical investigations using the numerical procedure of Runge-Kutta method of the 4th order.Parameters for a disease-free state are taken as follows: We have taken some parameter values from the literature.In addition, some are assumed with feasible values based on sufficient analysis and calculation of the conditions that satisfy the stability results.In particular, , µ 0 , β , γ and γ 2 are taken from 3 , while all other parameters are assumed.Clearly, in this case, the proposed prob- lem, as stated by equation (1), has only the infection-free state and is stable globally asymptotically (see Fig. 2).Note that for these values, the calculation of the basic reproductive number gives that R 0 = 0.203 implies that R 0 < 1 , so the stability results at the disease-free states hold.In addition, the theoretical interpretation states that if R 0 < 1 , each solution curve of S approximately taking five months to reach its associated equilibrium position as depicted in Fig. 2a.Similarly, the dynamics of the acute and chronic population are demonstrated in Fig. 2b,c, which describe that the solution curves of A and B take approximately ten months to reach the stable equilibrium position.On the other hand, the dynamics of the recovered population reveal that there will always be recovered individuals, as presented by Fig. 2d.Thus biologically, the results state that eliminating the hepatitis B virus from the community is subject to the threshold parameter's value.Whenever it is less than unity, the disease will be easily eradicated.
For investigating the stability of the proposed epidemiological model, we assume the same parameter values as in Eq. ( 36), except for α. .If α = 0.95 , then R 0 = 2.03 > 1 implies that the endemic state exists, as illustrated in Fig. 3, which clearly shows the results of the analytical analysis.. In this case, the biological investigations reveal that if no control measures are adopted appropriately, the disease will reach its associated endemic position.It could be noted that the susceptible portion of the population decreases from the initial and leads to zero over (31) = 0; whenever T = 0, (36) � = 0.0121, η = 0.8, µ 0 = 0.00693, v = 0.002, α = 0.95, γ = 0.16, γ 1 = 0.004, β = 0.33, γ 2 = 0.002, µ 1 = 0.8.time, as shown in Fig. 3a.However, there will always be an infected population, i.e., chronic and acute individuals, as shown in Fig. 3b,c, respectively.Similarly, we simulate the problem to study the dynamics of the recovered population as illustrated in Fig. 3d.The time dynamics of the recovered population state that the amount of recovered individuals decreases as time grows while leading to its associated endemic position.
We now perform the numerical investigation of the control problem to verify and support the theoretical examinations for optimal control analysis.We again use the 4th-order Runge-Kutta (RK) technique to perform the numerical simulations of the control problem.More precisely, we solve the system (23) via the 4th order RK scheme by taking the time unit from 0 to 50.We then solve the adjoint variables system as given by Eq. ( 32) with the help of backward RK procedure of the 4th order at the same interval of time along with the use of the transversality conditions stated by Eq. ( 33) as well as with the solution of the state system.To investigate the endemic state of the model, we use the same parametric values.The weight constants and initial conditions are, however, as follows: We then execute the above procedure with the aid of Matlab and obtain the graphical visualization as presented in Fig. 4, demonstrating the time dynamics of epidemiological groups of susceptible, acute, chronic, and recovered individuals.Our numerical results illustrate clearly the effect of applying the controls: to minimize acute and chronically infected populations while maximizing the recovered individuals.The illustration of the susceptible population with and without control is described in Fig. 4a.Moreover, Fig. 4b depicts the acute population with and without control.In a similar fashion, the dynamics of a chronic population with and without controls are shown in Fig. 4c.Finally, Fig. 4d visualizes the simulation of the recovered population with the application and without the application of controls.Further, Fig. 5 presents the control profiles.These analyses clearly, reflect the importance and the effectiveness of the implantation of the proposed control mechanism.We note that the numerical simulations presented in 31 , aiming to illustrate the usefulness of optimal control theory, are inconsistent.Indeed, in 31 , the authors solved an inconsistent control problem using two control measures: isolation and treatment.However, in the case of hepatitis B, where more than two billion people are infected, the control isolation is inconsistent according to the WHO guideline and is never exercised.
Additionally, we observe that the reproductive number is a key parameter, and when this quantity is greater than unity, the disease persists, while if it is less than one, it becomes extinct.Using the basic reproduction numbers as a reference, we conduct a sensitivity analysis of the model parameters against the control measures.Some parameters are directly proportional to the reproductive number, while others are negatively proportional.Vaccination and treatment controls are negatively correlated with reproductive quantity, such that whenever their values increase, the reproductive number decreases significantly, as shown in Fig. 6.Based on the parametric values given in Eq. ( 36), the sensitivity indexes of the vaccination and treatment control measures are calculated.It can be seen in Fig. 6a,b that if the control measures were increased by 10 percent, the reproductive number would decrease by 7.14 percent.Hence, we conclude that implementing the proposed control strategies in a true sense will result in the eradication of contagious hepatitis B virus infections.

Conclusions
Our study examined the dynamics of hepatitis B epidemics under acute and chronic transmission scenarios, both horizontally and vertically.After mathematically deriving and analyzing the proposed system, we determined its reproduction number to find the model equilibrium and stability.In the proposed problem, there are two equilibrium states: infection-free and disease-endemic.There is a detailed description of both steady states.Under certain conditions, the equilibria are stable.The global properties of the proposed epidemiological model were analyzed using the Lyapunov function theory and the geometrical approach.Additionally, we demonstrated that the proposed problem exhibits backward bifurcation.To develop an optimal control mechanism, the proposed problem incorporates two time-dependent control measures, vaccination, and treatment.According to Pontryagin's necessary conditions, there exists an optimal control mechanism for minimizing the infected (acute and chronic).We concluded by presenting numerical justifications and examining whether the derived analytical www.nature.com/scientificreports/findings are robust.Treatment and vaccination, along with the application of optimization theory, were found to be very effective in controlling hepatitis B virus infection.Accordingly, we do not recommend the isolation of individuals, as opposed to the results reported in 31 .Based on different perspectives of the results, which were investigated analytically and numerically, we concluded that one way to eradicate hepatitis B is to minimize the threshold quantity by keeping it below unity.The model also suggests that if hepatitis B persists, it will reach its 4. The temporal dynamics of the model with optimal control verse without control ( v = 0.02 ) against the parametric values (37).

Figure 1 .
Figure 1.The schematic diagram for the transmission of the disease.

Proposition 2 . 1 (
Existence and uniqueness) The proposed epidemiological model (1) with initial conditions described by Eq. (2) possesses a unique solution.

Figure 2 .
Figure 2. Solution curves of the system (1) around the disease free equilibrium against the parameters value given in Eq. (36) and for different initial sizes of population, where the value of threshold quantity (basic reproductive number), R 0 = 0.203 < 1.

Figure 3 .
Figure 3. Solution curves of system around the endemic equilibrium against the parameters values given in Eq. (36), except α = 0.95 and different initial sizes of population, which implies that R 0 = 2.03 > 1.

Figure 5 .
Figure 5.The plot represents the influence of the control measures on the basic reproductive number.