Modeling different infectious phases of hepatitis B with generalized saturated incidence: An analysis and control

1.   Introduction

The liver is a vital organ in a living body. Several diseases caused due to the consequences of liver infections. Hepatitis B is the liver inflammation produced by hepatitis B virus (HBV) ^[1]. Hepatitis B damages the liver cells and, as a result, produces liver inflammation ^[2,3]. Hepatitis B is a multi-infection disease causing acute and chronic hepatitis infection, which leads to a long-term existence of infection and risk of liver cancer. The initial stage is up to six months of the infection known as the acute HB-infection in which, usually, the immune system is capable of surviving the infection, but there may also be a chance for the infection to become more severe and lead to chronic hepatitis B infection ^[4]. More than two hundreds and fifty million individuals have chronic hepatitis B worldwide as per WHO information. Moreover, every year the infection of hepatitis B causes approximately 1 million deaths worldwide due to cirrhosis of liver and liver cancer. Particularly in the Western Pacific region and Africa, the burden of hepatitis B infection is high because an estimated 6.2% and 6.1% of the population, respectively, live with chronic hepatitis B infection. Nevertheless, the infection is not limited to the Western Pacific region and Africa only, but is a global health problem, and the infection rates are also significant in other ares such as Southeast Asia and the Eastern Mediterranean.

HBV transmits in various ways, but the critical transmission routes are vaginal secretions, transfusion of semen as well as blood, sharing razors without care, sexual interaction, and drug equipment contaminated with infected blood of the virus of hepatitis B. Another primary source of this virus is vertical or parental transmission, i.e., the infected mother can transmit it to their baby at birth. Children of ages 1–6 years have a 90% chance of getting the infection of hepatitis B after exposure to the virus ^[5]. However, the virus cannot transmit due to casual contact like drinking water, eating, kissing, and hugging. Similarly, the virus cannot be transmitted through general gatherings like in universities, schools, colleges, or other places ^[6,7]. Symptoms of this virus may include nausea, fatigue, vomiting, muscle and joint aches, yellow skin and dark urine, diarrhea, easy bruising, tiredness, etc.

Modeling the epidemiology of infectious disease is one of the important and emerging areas of applied mathematics as well as applied science that will be used to forecast the long dynamical behavior of various epidemics ^{[8,9,10,11,12,13,14,15]}. Several biologists, mathematicians, and researchers used the concept of mathematical modeling to investigate the dynamics of contagious diseases. Daniel Bernoulli was the first one in mathematical epidemiology to present a model describing smallpox dynamics in 1766 ^[16]. Kermack and Mckandrik presented the susceptible-infected-recovered (SIR) epidemiological model to represent the dynamics of infectious disease among three groups of the population ^[17]. In order to control the spread of HIV, a model with protection awareness is investigated by Zhai et al. ^[18]. As a global issue, numerous models are investigated to study the dynamics of HB ^[19,20,21]. Furthermore, a control strategy has been investigated by Medley et al. ^[22] to eliminate the contagious viral infection of HBV. Likewise, an age-structured model reported for the dynamics of HBV in China by Zhao et al. ^[23]. Similarly, the temporal dynamics of HBV with controls analysis have been investigated in ^[24,25]. Motivated by the work reported in ^[25], we study the dynamics of hepatitis B under the effect of various infectious phases and generalized incidence rate.

The disease incidence parameter is a critical term in epidemiological models that help researchers and policymakers correctly understand the burden and spread of diseases. Notably, it can measure disease frequency, helps with disease surveillance and assessment of disease prevention, enables resource allocation, evaluates the effectiveness of interventions, etc. The simple incidence rate is bilinear, $\beta{S}{I}$, which has been frequently exercised in numerous epidemiological models ^[26,27,28], where $\beta$ is the contact rate, $S$ is the amount of susceptible, and $I$ is the amount of infected. The concept of saturated incidence rate, $\frac{\beta{I}{S}}{1+\gamma{I}}$ has been introduced by Capasso et al. ^[29], where $\beta$ is the disease transmission ratio, while $\gamma$ is the parameter which captures the effect of behavioral changes on the spread of the disease. Afterward, many authors have been used the concept of saturated incidence rate while formulating epidemic models (see, for instance, ^[30,31]). Usually, when a disease outbreak occurs, individuals may change their behavior in various ways, such as practicing better hygiene, avoiding crowded places, or wearing masks, which ultimately reduce the transmission of the disease and slow its spread. The parameter $\frac{1}{1+\gamma{I}}$ measures the degree to which susceptible individuals change their behavior in response to the outbreak. Overall, the saturated incidence rate is more beneficial for studying the spread of infections and investigating behavioral changes' impact on disease transmission ^[32]. However, the saturated incidence rate may not work well for a disease with multiple infection phases because this formula assumes only a single infected group of individuals, which may not be accurate for infection with various phases or stages. We will try to fill this gap by introducing a more general form of saturated incidence rate as $\frac{\beta\left\{I_1+I_2+\ldots, I_n\right\}{S}}{1+\gamma\left\{I_1+I_2+\ldots, I_n\right\}}$, where $I_1, I_2, \ldots, I_n$ represent the different infection phases of the disease.

In this work, we construct an epidemiological model with a rate of generalized incidence to represent the dynamics of multi-infection disease of HBV. We develop the model according to the multi-infection phases of HB and divide the whole population into five subgroups of the compartmental population. We then discuss boundedness and positivity to show that the proposed epidemiological model is a well-posed dynamical system. We also find the reproductive quantity and model equilibria to investigate the qualitative analysis of the epidemic model. We discuss the local dynamical properties of the model with the aid of the linear stability approach, while to discuss the global properties of the problem, we use the Lyapunov theory. Based on the dynamics of the model, we then define a suitable control mechanism to control the infection of HBV transmission with the aid of optimal control theory. Three control measures, precautionary control measure, treatment of infected individuals, and vaccination control measure, are suggested to reduce the infected individuals with multi-infection phases and maximize the non-infected individuals. Finally, to show the validation of the analytical work and the effect of the control measures, we present the detailed numerical simulation using the well-known numerical procedure of the Runge-Kutta method.

The article structure is as follows: The detailed model formulation and its well-posedness are discussed in the Section 2, then in Section 3, we discuss the qualitative analysis of the model to derive the stability conditions. We then develop the control mechanism to eradicate the infection by taking the extended version of the model in Section 4. The detailed existence analysis and the charecterization of the optimal control problem are given in Section 5. Next, we present the model simulation to verify our theoretical analysis and show the effect of control measures implemented in Section 6. We give a conclusion in Section 7.

2.   The mathematical structure of the model

We develop the model to investigate the transmission dynamics of HBV. Since obviously hepatitis B is a multi-infection phase disease, we use the generalized saturated incidence $\frac{\beta\left\{I_1+I_2+\ldots, I_n\right\}{S}}{1+\gamma\left\{I_1+I_2+\ldots, I_n\right\}}$, while formulating the proposed model. Keeping in view the complex nature of the disease, we classify the various classes of population into sub-groups of susceptible, acute, chronic, hospitalized and recovered individuals denoted by $\mathcal{S}(t)$, $\mathcal{I}_a(t)$, $\mathcal{I}_c(t)$, $\mathcal{H}(t)$, and $\mathcal{R}(t)$, respectively. Since, $\mathcal{I}_a(t)$ and $\mathcal{H}(t)$ respectively represent the acute and hospitalized individuals, and usually do not transmit the disease to others, the disease transmission co-efficient for $\mathcal{I}_a(t)$ and $\mathcal{H}(t)$ are assumed to be zero. So, the generalized saturated incidence rate in this case looks like $\frac{\tau{\mathcal{S}(t)\mathcal{I}_{c}(t)}}{1+\theta_{1}\mathcal{N}_p(t)}$. In addition, because of the population dynamics, we assume that all the model states variables and parameters are non-negative values. Moreover, moving of the infected individuals (acute & chronic) leads to the hospitalized, while two types of recoveries are assumed as per HBV characteristics, natural for acute and due to treatment for the chronic portion of the population. We also assume the recovery is due to hospitalization. Natural death occurs in every population group, so the outflow of natural death is assumed in all groups of the compartmental population. In contrast, disease-induced death occurs due to the chronic infection. In addition, it is clear that vaccines for hepatitis B are available and effective, therefore we assume that the vaccinated individuals of the susceptible enter the recovered epidemiological group of the model. Thus, we present the model as follows:

with the initial compartmental population

In Eq (2.1), $\Pi$ is the birth rate, $\tau$ is the transmission coefficient of HBV, while $\theta_2$ is the parameter that describes the rate at which acutely infected individuals enter the chronic group. $\nu$ is assumed to be the vaccination rate. The parameter $\tau_1$ is the rate at which the acutely infected population enters the hospitalized group of individuals, and $\rho_1$ represents the death rate produced from the HBV infection. The individuals who move from the chronic to the hospitalized group are symbolized by $\tau_2$ while $\rho_2$ is also the death caused by the HBV infection in a hospitalized group. Moreover, to define the recovery of hospitalized individuals, we demonstrate it by $\tau_3$, and the natural death of all groups of the compartmental population is denoted by $\alpha_0$. In addition, $\mathcal{N}_p(t)$ is the sum of all infected population, i.e., $\mathcal{N}_p(t) = \mathcal{I}_a(t)+\mathcal{I}_c(t)+\mathcal{H}(t)$.

2.1.   Validity of the model

In this subsection, we will provide the validity of the proposed model against regional data taken from local hospitals of district Swat, Khyber Pakhtunkhwa, Pakistan, as shown by Figure 1(a). Since, according to 2017, the total population of district Swat is 2.31 million and thus the demographic parameter is calculated as $\Pi/\alpha_0 = 2,310,000$, where $\alpha_0 = (1/67.7)$ is the average life span in Pakistan. In addition, the procedure of ordinary least square (OLS) will be applied to obtain the best fit and to minimize the error between the reported data and the proposed model solutions. For this purpose, we use the objective function is given by

where $x_i$ and $\hat{x}_i$ respectively represent the actual cumulative hepatitis B reported cases and the associated model solution, while $n$ is the number of actual data points. Using the optimization algorithm, while updating the parameters values to derive better agreement with the real data and to minimize the error. The process is repeated until reaching the best model fit, as in Figure 1(b), which shows the validity of the model in the case of application to a real scenario.

2.2.   Positivity and boundedness

The proposed epidemiological model (2.1) represents the dynamics of compartmental populations; therefore, we need to check whether the state of the considered problem is non-negative for all $t$ with the initial compartmental population given in Eq (2.2). We also investigate that the considered model is bounded. For this, we study the results as given below.

Theorem 2.1. All the solutions $\left(\mathcal{S}(t), \mathcal{I}_a(t), \mathcal{I}_c(t), \mathcal{H}(t), \mathcal{R}(t)\right)$ of the model (2.1) with the initial compartmental sizes (2.2) remain non-negative and are uniformly bounded in the positively invariant region $\Omega$ for all non-negative $t$.

Proof. The solution of the model (2.1), first equation, can be written as

where $\psi(t) = \left\{\frac{\tau \mathcal{I}_c(t)}{1+\theta_1\mathcal{N}_p(t)}-(\alpha_0+\nu)\right\}\mathcal{S}(t)$. From the epidemiological model (2.1), second equation, we have

which leads to

In a similar fashion, the third equation of the model takes the following form:

After that, integration gives

By following the same steps, we may write the fourth and fifth equations of the model as

which implies that

and

We observed from the above equations that the system (2.1) satisfying the conditions (2.2) remains non-negative for every $t\geq0$.

To proceed further, let $\mathcal{N}(t)$ demonstrate the size of whole population, then

which yields that

From the Eq (2.8), we note that $\mathcal{N}(t)$ becomes less than or equal to $\frac{\Pi}{\alpha_0}$ as time grows unboundedly. Thus it follows that the total population in the region $\{R^5_+\backslash\{0\}\}$ is bounded by $\frac{\Pi}{\alpha_0}$ with growing time ($t\rightarrow \infty$), therefore, the solution trajectories of the model satisfying the initial conditions are bounded.

Since the model state variables are non-negative and $\mathcal{N}(t)\leq\frac{\Pi}{\alpha_0}$, it is implied that the proposed problems (2.1) and (2.2) is well-posed. Further, we assume that $\mathcal{N}(0)\leq\frac{\Pi}{\alpha_0}$, then from Eq (2.8), we conclude that $\mathcal{N}(t)\leq\frac{\Pi}{\alpha_0}$, and thus every solution of the proposed epidemic problem with initial conditions in $R^5_{+}$ remains in $\Omega$ as

3.   Qualitative analysis

We discuss the dynamics of the proposed model by investigating the model equilibria. Clearly, the disease-free state (DFE) of the system (2.1) is represented by $\mathcal{E}_{df}$ and written as $\mathcal{E}_{df} = (\mathcal{S}^0, \mathcal{I}^0_a, \mathcal{I}^0_c, \mathcal{H}^0, \mathcal{R}^0) = \left(\frac{\Pi}{\alpha_0+\nu}, 0, 0, 0, \frac{\nu \Pi}{\alpha_0(\alpha_0+\nu)}\right)$. Before calculating the hepatitis B endemic state, we find the reproductive parameter, denoted by $\mathcal{R}_0$, defined to be the threshold quantity (basic reproductive number), while demonstrating the average of newly infected caused by an infected after introducing them into a susceptible population. In the case of classical epidemiological models, when $\mathcal{R}_0 < 1$, the disease dies out. If $\mathcal{R}_0 > 1$, then it will be expected that the disease is spreading. We use the next-generation method to investigate this quantity as reported in ^[33,34]. Upon using the same methodology as adopted by ^[34], we assume that $\mathcal{Y} = \left(\mathcal{I}_a, \mathcal{I}_c, \mathcal{H}\right)^T$. Then

where

The reproductive number, $\mathcal{R}_0$ is defined to be the spectral radius of the matrix $\mathcal{F}\mathcal{V}^{-1}$, thus it is given by

where $\psi_1 = \alpha_{0}+\nu$, $\psi_2 = \alpha_{0}+\theta_{2}+\kappa_{1}+\tau_{1}$, and $\psi_3 = \alpha_{0}+\kappa_{2}+\rho_{1}+\tau_{2}$. The reproductive number is a dimensionless rate representing the average of secondary hepatitis B cases produced whenever a hepatitis B-infected person is introduced into the susceptible population. Thus, it is clear that if $\mathcal{R}_0 < 1$ and the initial sizes of the population's compartments are in the hepatitis-free state, then the hepatitis B disease vanishes. For this, we prove the subsequent result.

Theorem 3.1. If $\mathcal{R}_{0} < 1$, then the epidemiological problem (2.1) is stable locally and globally at the hepatitis-free state, $\mathcal{E}_{df} = \left(\frac{\Pi}{\alpha_{0}+\nu}, 0, 0, 0, \frac{\nu \Pi}{\alpha_{0}(\alpha_{0}+\nu)}\right)$, otherwise $\mathcal{E}_{df}$ is unstable and is a saddle point.

Proof. We calculate the linearized matrix of the system (2.1) at hepatitis-free state $(\mathcal{E}_{df})$ as

Calculating the eigenvalues of $J(\mathcal{E}_{df})$ implies that three eigenvalues are negative, i.e., $\psi_1$, $\psi_4$ and $\alpha_{0}$ are negative. For the remaining two eigenvalues, we take a $2\times2$ matrix given as

It is enough to show that the trace of $\mathcal{B}$ is $\mathcal{B} < 0$ and the determinant of $\mathcal{B}$ is $\mathcal{B} > 0$ for the Routh-Hurwitz criteria, thus whenever $\mathcal{R}_{0} < 1$, then

and

imply that the above criteria hold subject to the condition of $\mathcal{R}_{0} < 1$. Clearly, $trace(\mathcal{B}) < 0$ and $determinant(\mathcal{B}) > 0$ when $\mathcal{R}_{0} < 1$. However, in case of $\mathcal{R}_{0} > 1$, the $trace(\mathcal{B}) < 0$ as well as $determinant(\mathcal{B}) < 0$, which implies that the eigenvalues of $B$ have the alternative sign, i.e., positive as well as negative, so the $\mathcal{E}_{df}$ is an unstable equilibrium point.

We calculate the global dynamics properties of the model (2.1) at $\mathcal{E}_{df}$, and therefore define a Lyapunov function, such that

where $\lambda_{i}$ for $i = 1, 2, 3, 4$ are constants assumed to be positive. The derivative of the function (3.2) with the use of values from the model (2.1), leads to

By assuming the positive constants in such a way that $\lambda_{1} = \lambda_{2} = \psi_1$ and $\lambda_{3} = \tau\Pi$ in Eq (3.3), we get the following equation:

Algebraic manipulation gives that

which implies that

Therefore, $\frac{d\varphi}{dt}$ is negative, if $\mathcal{R}_0\leq 1$. Moreover, $\frac{d\varphi}{dt} = 0$, whenever $\mathcal{S} = \mathcal{S}^0, \mathcal{I}^0_a = 0, \mathcal{I}^0_c = 0, \mathcal{H}^0 = 0$.

Next, we examine the properties of the hepatitis B endemic state for the considered epidemic problem. To shorten our calculation, we assume that $\psi_4 = \alpha_0+\rho_2+\tau_3$ and $\mathcal{E}_{ee} = \left(\mathcal{S}^\sigma, \mathcal{I}_a^\sigma, \mathcal{I}_c^\sigma, \mathcal{H}^\sigma, \mathcal{R}^\sigma\right)$ is the endemic state of the model, then

To present the dynamics of the hepatitis B endemic state, we prove the result as given below.

Theorem 3.2. The proposed problem (2.1) is stable at $\mathcal{E}_{ee} = (\mathcal{S}^\sigma, \mathcal{I}^\sigma_{a}, \mathcal{I}^\sigma_{c}, \mathcal{H}^\sigma, \mathcal{R}^\sigma)$, if $\mathcal{R}_{0} > 1$.

Proof. Calculating the Jacobian at the $\mathcal{E}_{ee}$ of the proposed model, we may obtain

It is clear from the above matrix that $-\alpha_{0}$ and $-\psi_4$ are the two negative eigenvalues. To find the eignvalue's nature, we may define

Calculating the characteristic polynomial of the above matrix $G$, we have

where

All $k_i > 0$ for $i$ = 1, 2, 3 and $k_1k_2 > k_3$, which ensure the criteria of Routh-Hurwitz because $k_1 > 0$, $k_3 > 0$ and $k_1k_2 > k_3$, whenever $\mathcal{R}_0 > 1$. Thus, $f(y)$ has negative roots whenever $\mathcal{R}_0 > 1$, and so $\mathcal{E}_{ee}$ is the stable state of the model that is under consideration.

We now discuss the properties of global analysis of the proposed problem (2.1) at the hepatitis B endemic state. To investigate the global properties of the problem at the hepatitis B endemic state, we define the following function:

The temporal derivative of Eq (3.5) with respect to $t$, and using model (2.1), gives

which implies that

or equivalently we can write

The re-arrangement with full simplification leads to the following assertion:

By simplifying and re-writing, we obtain

Hence, $\frac{d\xi}{dt} < 0$ for all $\mathcal{S}$, $\mathcal{I}_a$, $\mathcal{I}_c$, $\mathcal{H}$, $\mathcal{R}$, and $\frac{d\xi}{dt} = 0$ at the endemic state, so the hepatitis B endemic state $\mathcal{E}_{ee}$ is the positively invariant set only containing $\left\{\left(\mathcal{S}, \mathcal{I}_a, \mathcal{I}_c, \mathcal{H}, \mathcal{R}\right):\mathcal{S} = \mathcal{S}^\sigma\right.$, $\left.\mathcal{I}_a = \mathcal{I}_a^\sigma, \mathcal{I}_c = \mathcal{I}_c^\sigma, \mathcal{H} = \mathcal{H}^\sigma, \mathcal{R} = \mathcal{R}^\sigma \right\}$, which implies that $\mathcal{E}_{ee}$ is a stable state.

4.   Optimal control

Optimal control theory is a powerful mathematical technique through which we can develop control strategies to control different infectious diseases, i.e., hepatitis B virus. We make a control mechanism for the infection of hepatitis B with the objective to reduce the number of infective by taking into account the maximization of $\mathcal{S}(t)$ as well as $\mathcal{R}(t)$ populations while to minimizing $\mathcal{I}_a(t)$ and $\mathcal{I}_c(t)$. For this purpose, three control measures dependent over time will be used, i.e., $u_i(t)$, $u_t(t)$, and $u_v(t)$, which physically represent the preventive measures of hepatitis B, treatment of infected individuals, as well as vaccination, respectively. Clearly, there are five state variables, i.e., $\mathcal{S}(t)$, $\mathcal{I}_a(t)$, $\mathcal{I}_c(t)$, $\mathcal{H}(t)$, and $\mathcal{R}(t)$, therefore, we then assume the above three control measures to design the control problem as

subject to

with

In the earlier section, i.e., in Section (2.1), the description of parameters is discussed in detail, while in the objective functional (4.1), $C_1$, $C_2$, and $C_3$ illustrate the weight constants for the proposed control strategies. The weight constants, $D_1$, $D_2$, and $D_3$ are the constants relating to the control measures of preventive measures, treatment of infected individuals, and vaccination, respectively. Moreover, the terms $\frac{1}{2}D_1u_{i}^2(t)$, $\frac{1}{2}D_2u_t^2(t)$, and $\frac{1}{2}D_3u_v^2(t)$ demonstrate the associated cost with the control measures. Thus, we wish to find the control measures that minimize the objective function, such that

subject to the system (4.2). The control set is described by

5.   Existence analysis

We discuss the existence analysis, and therefore assume the control system as stated by Eq (4.2) with initial conditions. It is very clear that for Lebesgue measurable and bounded control measures, initial conditions (positive) and bounded solutions (positive) to the proposed system exist. Moreover, we go back to the control problem as stated by Eq (4.2), as well as Eq (4.4), to figure out the optimal solution. We define the Lagrangian first and then the Hamiltonian for the said purposes, i.e., for the optimal problem (4.2) and (4.4). Consequently, the Lagrangian takes the form

In addition, for the minimal value of $\mathcal{L}_{opt}$, we define the Hamiltonian, such that

Thus, for the existence of such controls, first we prove the existence; therefore, regarding the existence, we illustrate the following result.

Theorem 5.1. There exists an optimal control $u_{opt}^\sigma = (u_{i}^\sigma, u_t^\sigma, u_v^\sigma)\in M$, such that

subject to the control system as reported in Eq (4.2).

Proof. Following the same methodology as used in ^[35,36], we explore the existence analysis of the optimal control functions. Clearly, the state, as well as control, are non-negative values. So, the necessary condition of convexity for the objective functional (4.1) over $u_{i}(t)$, $u_t(t)$, and $u_v(t)$ holds. The control variable set $u_{i}, u_t, u_v \in M$ is obviously closed as well as convex and so implies the optimal system's boundedness. This grants compactness. Also, the integrand $C_1\mathcal{S}(t)+C_2(t)\mathcal{I}_a(t)+C_3\mathcal{I}_c(t)+\frac{1}{2}\left\{D_1u_i^2(t)+D_2u_t^2(t)+D_3u_v^2(t)\right\}$ is convex over $M$.

Now the optimal solution will be determined to the proposed control problem. For this, the Pontryagin Maximum Principle will be utilized. If

where $y = (\mathcal{S}(t), \mathcal{I}_a(t), \mathcal{I}_c(t), \mathcal{H}(t), \mathcal{R}(t))$ and $u = (u_i, u_t, u_v)$, while $f$ is the Lagrangian of the objective function (4.1) and $g = (g_1, g_2, g_3, g_4, g_5)$ represent the right-hand side of the first, second, third, fourth, and fifth equations of the control system (4.2). Thus, if $\left(y^\sigma, u^\sigma\right)$ is an optimal solution, then a non-trivial vector $\chi(t)$ (a set of adjoint variables) exists, such that $\chi(t) = \left(\chi_1(t), \chi_2(t), \ldots, \chi_5(t)\right)$, satisfying

Moreover, the necessary condition will be applied to the Hamiltonian in terms of the following results.

Theorem 5.2. Let $\mathcal{S}^\sigma$, $\mathcal{I}_a^\sigma$, $\mathcal{I}_c^\sigma$, $\mathcal{H}^\sigma$, and $\mathcal{R}^\sigma$ be the optimal solution with associated optimal control measures $\left(u_{i}^\sigma, u_t^\sigma, u_v^\sigma\right)$ for the problem (4.2)–(4.4), then the adjoint variables $\chi_1(t)$, $\chi_2(t)$, $\chi_3(t)$, $\chi_4(t)$, and $\chi_5(t)$ exist and satisfy

with boundary (transversality) conditions

Further, the optimal control measures $u_{i}^\sigma(t)$, $u_t^\sigma(t)$, and $u_v^\sigma(t)$ are defined as

Proof. We use the Hamiltonian (5.3) for finding the adjoint system (5.5) as well as the transversality condition (5.6). We set $\mathcal{S}(t) = \mathcal{S}^\sigma$, $\mathcal{I}_a(t) = \mathcal{I}_a^\sigma$, $\mathcal{I}_c(t) = \mathcal{I}_c^\sigma$, $\mathcal{H}(t) = \mathcal{H}^\sigma$, and $\mathcal{R}(t) = \mathcal{R}^\sigma$, while the differentiation of $\mathcal{H}_{pmp}$ with respect to $\mathcal{S}(t)$, $\mathcal{I}_a(t)$, $\mathcal{I}_c(t)$, $\mathcal{H}(t)$, and $\mathcal{R}(t)$ leads to the system (5.5). Moreover, to get $u_{i}^\sigma$, $u_t^\sigma$, and $u_v^\sigma$, $\mathcal{H}_{pmp}$ will be differentiated respectively with respect to $u_{i}$, $u_t$ and $u_v$, and then the solution of $\frac{\partial \mathcal{H}_{pmp}}{\partial u_{i}} = 0$, $\frac{\partial \mathcal{H}_{pmp}}{\partial u_t} = 0$, and $\frac{\partial \mathcal{H}_{pmp}}{\partial u_v} = 0$, in the interior of control set with the application of optimality condition. In the end, the use of the control property $M$ gives the optimal value of the control variables.

We recall the formula $u^\sigma = (u_{i}^\sigma, u_t^\sigma, u_v^\sigma)$ to characterize the optimal control problem, which consists of the state system with the initial sizes of compartmental populations, the adjoint system with terminal conditions and the optimal measures. We further use the iterative procedure to solve the proposed optimal control problem.

6.   Numerical simulations

We illustrate the numerical findings of the analytical results to demonstrate the feasibility of the derived results graphically. We solve the proposed problem via the Runge-Kutta method of the 4th order, along with various initial sizes of population and different sets of parameter values as given in the captions of the figures. Moreover, the time unit is taken to be 0–100. Moreover, the parameters values are taken in accordance with the theoretical results that have been carried out in Theorems 3.1 and 3.2, and are given at the captions of the figures. We generate the following graphs as shown in Figures 2 and 3, which are respectively illustrating the verification of analytical results around disease–free and endemic states. More specifically, Figure 2 represents the dynamics of the compartmental populations of the model around the hepatitis free state, while the disease endemic state dynamics are represented by Figure 3. In Figure 2, the graphs (a)–(e) respectively describe the dynamics of susceptible, acute, chronic, hospitalized, and recovered populations, which show that all other individuals vanish except susceptible and the recovered populations, whenever $\mathcal{R}_0 < 1$. However, if $\mathcal{R}_0 > 1$, the disease may reach the endemic state as shown in Figure 3, implies that the infected individuals will always persist and there is a need for interventions strategies to control the transmission of the contagious disease of HBV.

In addition to showing the effect of the proposed control mechanism by presenting the simulation of the control analysis, we want to differentiate between the control and without control implementation. Therefore, we use the forward Runge-Kutta method of the 4th order to solve the state system (4.2), while the adjoint system (5.5) will be solved with the help of backward Runge-Kutta method of the 4th order. All the values for parameters are given in the captions of the figures representing the control implementation analysis. Thus, as a result, we obtain the graphs as represented by Figure 4, which describe the dynamics of susceptible (Figure 4(a)), acute (Figure 4(b)), chronic (Figure 4(c)), hospitalized (Figure 4(d)), and recovered (Figure 4(e)) individuals with and without control measures. The graphs clearly illustrate the effect of control strategies: to reduce the infected and to increase the recovered populations. The difference between the two cases is clearly visible. We observed that the collective implementation of the proposed control measures optimally leads to eradicating the contagious disease of HBV transmission.

7.   Conclusions

In this work, we presented a more generalized epidemiological model for HBV transmission by including the new features according to the characteristic of the disease. The incidence parameter plays an essential role in the dynamics and control of biological models; therefore, we use the generalized saturated incidence rate $\frac{\beta\{I_1+I_2+I_3+\ldots+I_n\}S}{1+\gamma\{I_1+I_2+I_3+\ldots+I_n\}}$ to study the temporal dynamics of hepatitis B, which is more suitable as compared to traditional incidence rates. Because hepatitis B is a multi-infection disease, the traditional saturated incidence rate $\frac{\beta{SI}}{1+\gamma{I}}$ is not appropriate while investigating the dynamics of HBV transmission. Thus, considering the disease's characteristics, we formulated the model and discuss the feasibility of the problem. We then calculate the disease-free equilibria and consequently the basic reproductive number $\mathcal{R}_0$ with the help of a well-known technique of the next-generation matrix approach. In a similar fashion, the endemic state is also derived using the reproductive number and then the detailed stability analysis is discussed via various approaches to derive the stability conditions. For this, the linearization, as well as the Lyapunov theory, are retrieved to discuss the local and global properties of the newly constructed model. In addition, very importantly, we then use three control measures and design a control mechanism with the aid of optimal control theory that, how to eliminate the infection of hepatitis B. At last, all the theoretical as well as analytical findings are supported via graphical representations with the help of numerical experiments to show the validity of the model and the effects of control implementation.

In the future, we will consider the protection awareness to separate the susceptible individuals into two groups, i.e., susceptible with protection awareness and susceptible without protection awareness.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by Research Supporting Project Number (RSP2024R421), King Saud University, Riyadh, Saudi Arabia. The work has been also supported by the UAE University, fund No. 12S107. Further, this work is supported by the University Innovation Foundation of China (Grant No. 2022IT101) and the Basic Education Quality Improvement Research Center project Foundation of Xinjiang province (Grant No. WKJDJSZD23002).

Conflict of interest

The authors declare there is no conflict of interest.