Mathematical model analysis and numerical simulation of intervention strategies to reduce transmission and re-activation of hepatitis B disease

Introduction

The hepatitis B virus causes hepatitis B, an infectious disease that affects the liver (HBV). It is spread through contact with infectious blood, sperm, and other bodily fluids.¹ Furthermore, it can be passed on from infected women to their infants at the moment of birth, or from family members to their children in early childhood.² Because it is asymptomatic by nature, it develops complicated and can lead to chronic liver disease. People in this chronic stage are at a significant risk of dying from liver cirrhosis and cancer. Hepatitis B, on the other hand, is a liver infection caused by the HBV that can be prevented with a vaccine.³

For numerous decades, experts have paid close attention to the study of the hepatitis B virus. One of the reasons for this close scrutiny is a strong desire to learn everything there is to know about this lifethreatening virus and how it spreads and spreads. Risk of HBV also goes to healthcare workers who sustain accidental needle-stick injuries while caring for HBV infected people.⁴ Although it is asymptomatic and have ten times mode of transmission than HIV/AIDS, asafe and effective vaccine is available to prevent HBV infection.⁵ Another feature of HBV disease is that it might reactivate after malignancy, autoimmune disease, or organ transplantation immunosuppressive medication.

In the study of virus dynamics, mathematical modeling has proven to be the most effective method for understanding biological mechanisms and interpreting experimental results. To better understand the dynamics of HIV, HBV, and other virus infections, an early mathematical model for the basic dynamics of virus in vivo was devised and examined.⁶ Hepatitis B is a serious liver illness caused by the HBV, which is a major worldwide medical problem and the most frequent kind of viral hepatitis. It is found only in Asia and Africa.⁷ Around the world, around two billion people are asymptomatic or symptomatic of the virus, and approximately 360 million have the most severe infection.⁸

In 1990, 154 million people, or 2.9 percent of the global population, were migrants; however, the figure for 2013 was 3.2 percent,⁹ and these figures do not include undocumented migrants.

The United States, Canada, and a few other EU countries were among the top ten destinations for international migrants in 2013.¹⁰

According to research on some African migrants in the United States, the proportion of people infected with HBV is ten times that of the host population. This is most likely due to a mix of factors, including a lack of understanding of diseases, risk factors, and symptoms, as well as a lack of access to healthcare and health information. According to studies, workers and displaced persons are 25 times more likely to die from liver disease than the general population.^11,12 proposed a mathematical model for hepatitis B with migration. According to the findings, rigorous immigration rules such as screening and limiting the number of immigrants allowed into a given population could help curb the spread of the disease.

Reference 13 also analyzed the state of the art in modeling and anazing data received from hepatitis B virus infected individuals treated with antiviral medications. With more understanding and quantitative methods, it will be easier to test new medicines for antiviral and immune-modulating effects, and viral kinetic studies may eventually be used to predict long-term patient responses.¹² Investigated the role of infectious immigrant epidemic models and immunization on disease patterns. This straightforward configuration takes into account the likelihood of acquired immunity with a focus on SIR and SIS models Class has been immunized.

One of the main goals of research into hepatitis B virus (HBV) infection is to enhance control and, eventually, eliminate the infection from the community. Mathematical models can be an effective tool in this strategy, allowing us to optimize the use of finite resources or simply to improve the effectiveness of infection control methods. To demonstrate the effects of climate change, Anderson and May employed a simple mathematical model.¹⁴ Reference 15 created mathematical models of HBV transmission dynamics and optimal vaccination and treatment control for both vertical and horizontal HBV transmission. However, they failed to account formigration effect and HBV reactivation following recovery and re-infection of the recovered group. In this study, we looked at re-infection, which occurs when previously healthy people get reinfected after undergoing liver transplantation and chemotherapy and are exposed to HBV and migration of population which is currently becoming harsh in increasing number of infected population of infectious disease in ones country.

Mathematical model formulation and assumptions

Susceptible S(t), Vaccinated V(t), Exposed E(t), Acutely infected A(t), Chronically infected C(t), Treated T(t), and Recovered R are the seven compartments that the total population N(t) is divided into based on the epidemiological status of individuals.

Corresponding compartmental diagram

Figure 1. Corresponding flow chart of extended model SVEACTR model.

Corresponding dynamical system

Mathematical analysis of the model

Under this section, we will check the positivity, boundedness, and existence of the solution of the model.

Theorem 1 (Positivity)

Let the initial data for the model be $S_0>0,V_0>0,E_0>0,A_0>0,C_0>0,T_0>0,R_0>0.$ Then, the solutions $S\left(t\right),V\left(t\right),E\left(t\right),A\left(t\right),C\left(t\right),T\left(t\right)$ and $R\left(t\right)$ of the model will be remaining positive for all time $t>0.$

Proof: Let $S_0>0,V_0>0,E_0>0,A_0>0,C_0>0,T_0>0,R_0>0$. Also, we assume that all the parameters are positive. To show this, we take these differential equations of the dynamical system given above and show that their solutions are non-negative as follows.

1. Let us take the first differential equation

After simplification, we get:

Since any exponential function with positive coefficients and $S_0$ are both positive and $\left(\frac{\left(\left(1-\theta \right)\pi +\left(1-\psi \right)\Lambda +\mathit{\tau V}+\mathit{\epsilon R}\right)}{\mu +\mathit{c\omega }\left(\frac{\mathit{\phi C}+A}{N}\right)}\right)$ is also positive, we can conclude that $S\left(t\right)>0.$

2. Let us consider the second ordinary differential equation $\frac{\mathit{dV}}{\mathit{dt}}=\mathit{\theta \pi }+\psi \Lambda -\mathit{\tau V}-\mathit{\mu V}.$

This equation is in the form of $\frac{\mathit{dy}}{\mathit{dt}}+p\left(x\right)y=q\left(x\right)$ which is first-order linear differential equation with integrating factor $I.F=e^{\int \left(\tau +\mu \right)\mathit{dt}}=e^{\left(\tau +\mu \right)t}$

Solution of $\frac{\mathit{dy}}{\mathit{dt}}+p\left(x\right)y=q\left(x\right)$ is given by:$y=\frac{1}{I.F}\int \left(I.F\times q\left(x\right)\right)\mathit{dx}$

So, the solution of $\frac{\mathit{dV}}{\mathit{dt}}+\left(\tau +\mu \right)V=\mathit{\theta \pi }+\psi \Lambda$ is given by:

Since $V_0$ is positive, any exponential function with positive coefficient is positive and $\frac{\mathit{\theta \pi }+\psi \Lambda }{\left(\tau +\mu \right)}$ is also positive, then we can conclude that $V\left(t\right)>0.$

3. Let us consider the third ordinary differential equation

After computing and simplification, we get:

Since $E_0$ is positive, any exponential function has range of positive real number and their product is positive and $\frac{\mathit{c\omega }\left(\frac{\mathit{\phi C}+A}{N}\right)S+\mathit{\gamma c\omega }\left(\frac{\mathit{\phi C}+A}{N}\right)R}{\left(\mu +c_1+c_2\right)}$ is also positive, then we can conclude that $E\left(t\right)>0.$

4. Let us consider the fourth differential equation

Computation and simplification gives us:

Since $A_0$ is positive, any exponential function with positive coefficient is positive and $\frac{c_{2}E}{\left({\delta }_2\alpha +{\delta }_1\alpha +\mu \right)}$ is also positive, then we can conclude that $A\left(t\right)>0.$

5. Let us consider the fifth ordinary differential equation

After some steps, we get:

Since $C_0$ is positive, any exponential function with positive coefficient is positive and $\frac{c_{1}E+{\delta }_2\mathit{\alpha A}}{\left(\mu +d+\sigma \right)}$ is also positive, then we can conclude that $C\left(t\right)>0.$

6. Let us consider the sixth ordinary differential equation

Simplification after some step yields:

Since $T_0$ is positive and any exponential function with positive coefficient is positive and $\frac{\mathit{\sigma C}}{\left({\phi }_1+\mu \right)}$ is also positive, then we can conclude that $T\left(t\right)>0.$

7. Let us consider the seventh ordinary differential equation

After computation and simplification, we get:

Since $R_0$ is positive, any exponential function with positive coefficient is positive and $\frac{{\phi }_{1}T+{\delta }_1\mathit{\alpha A}}{\left(\mathit{\gamma c\omega }\left(\frac{\mathit{\phi C}+A}{N}\right)+\epsilon +\mu \right)}$ is also positive, then we can conclude that $R\left(t\right)>0.$

This completes the proof of the theorem. Therefore, the solution of the model is positive.

Theorem 2 (Boundedness)

To show the boundedness of the solution, we have to show lower bound and upper bound. But, initially, $N\left(0\right)=N_0>0,S\left(0\right)=S_0>0,V\left(0\right)=V_0>0,E\left(0\right)=E_0>0,A\left(0\right)=A_0>0,C\left(0\right)=C_0>0,T\left(0\right)=T_0>0,R\left(0\right)=R_0>0$

These initial conditions are considered as lower bounds. Now, we are going to show the upper bound. By taking the relation

and differentiating both sides of the equation with respect to time, we get:

In the absence of mortality due to Hepatitis B virus, that is at $d=0$,

$\Rightarrow \frac{\mathit{dN}}{\mathit{dt}}=V+\pi -\mathit{\mu N}$, Integrating both sides give us:

$\Rightarrow \frac{-1}{\mu }ln\left(\Lambda +\pi -\mathit{\mu N}\right)\le t+c$ which is simplified to

$\Rightarrow \Lambda +\pi -\mathit{\mu N}\ge A{e}^{-\mathit{\mu t}}$ where $e^{-\mathit{\mu c}}=A=\text{Constant}$

By applying the initial condition $N\left(0\right)=N_0$, we get:

$\Rightarrow \Lambda +\pi -\mathit{\mu N}\ge A$ which upon substitution in

Then, $-\mathit{\mu N}\ge -\Lambda -\pi +\left(\Lambda +\pi -\mu N_0\right)e^{-\mathit{\mu t}}$

As $t\to \infty ,$ the population size $N\to \frac{\Lambda +\pi }{\mu }$ since $\frac{\left(\Lambda +\pi -\mu N_0\right)e^{-\mathit{\mu t}}}{\mu }\to 0$, as $t\to \infty .$

This implies that $0\le N\le \frac{\Lambda +\pi }{\mu }.$ Thus, feasible solution set of the system equation of the model enters and remains in the region.

Therefore, the basic model is well posed epidemiologically and mathematically. Hence, it is sufficient to study the dynamics of the basic model in $\Omega$.

Disease-free equilibrium point

To find the disease-free equilibrium (DFE), we equated the right-hand side of the model (1)-(7) to zero, evaluating it at E=A=C=T=0 and solving for the noninfected and noncarrier state variables. Therefore, the disease-free equilibrium is given by:

Effective reproduction number

We calculate the basic reproduction number denoted by $R_{\mathit{Eff}}$ using the van den Driesch and Warmouth next-generation matrix approach from Ref. 16. The basic reproduction number is obtained by taking the largest (dominant) eigenvalue (spectral radius) of the matrix: $F{V}^{-1}=\left[\partial {\mathcal{F}}_i\left(\frac{E^0}{\partial x_j}\right)\right]{\left[\partial v_i\left(\frac{E^0}{\partial x_j}\right)\right]}^{-1}$, where ${\mathcal{F}}_i$ is the rate of appearance of new infection in compartment$i$, $v_i$ is the transfer of infections from one compartment $i$ to another, and $E^0$ is the disease-free equilibrium point. The reproduction number $R_{\mathit{Eff}}$ of our HBV model is obtained by rearranging the differential equation of the dynamical system (1)-(7) above in terms of $\frac{d{X}_i}{\mathit{dt}}={\mathcal{F}}_i-v_i={\mathcal{F}}_i-\left(v_i^{-}-v_i^{+}\right)$.

Then.

Then,

Eigenvalues of $F{V}^{-1}$ is given by:

Then, the spectral radius (Effective reproduction number, $R_{\mathit{Eff}}$) of $F{V}^{-1}$ of our HBV model is:

Local stability analysis of the disease-free equilibrium point

Theorem 3: The disease-free equilibrium point of the model (1)-(7) above is locally asymptotically stable if $R_{\mathit{Eff}}<1$ and unstable if$R_{\mathit{Eff}}<1$.

Proof: The local stability of the disease-free equilibrium point of the system (1)-(7) can be studied from its Jacobian matrix at the disease-free equilibrium point $E^1=\left(S^0,V^0,E^0,A^0,C^0,T^0,R^0\right)=\left(\left(\frac{\Lambda +\pi }{\mu }-\frac{\mathit{\Lambda \psi }+\mathit{\theta \pi }}{\left(\tau +\mu \right)}\right),\frac{\mathit{\Lambda \psi }+\mathit{\theta \pi }}{\tau +\mu },0,0,0,0,0\right)$ and Routh-Hurwitz stability criteria. Then, the Jacobian matrix of the dynamical system (1)-(7) at $E^1=\left(\left(\frac{\Lambda +\pi }{\mu }-\frac{\mathit{\Lambda \psi }+\mathit{\theta \pi }}{\left(\tau +\mu \right)}\right),\frac{\mathit{\Lambda \psi }+\mathit{\theta \pi }}{\tau +\mu },0,0,0,0,0\right)$ is given by:

Then, the characteristic equation of the above Jacobian matrix is given by:

Let $m=\tau +\mu ,f=\frac{\mathit{c\omega S}}{N},l=\mu +c_1+c_2,g={\delta }_1\alpha +{\delta }_2\alpha +\mu ,h=\mu +d+\sigma ,i={\phi }_1+\mu ,j=\epsilon +\mu ,k=\frac{\mathit{\text{cωφS}}}{N},p={\delta }_2\alpha ,q={\delta }_1\alpha$

Here, it is obvious that, ${\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5$ are all negative.

To check the remaining eigenvalues are negatives for the quadratic equation $L_2{\lambda }^2+L_1\lambda +L_0=0$

Let us consider ${\lambda }^2+\left(l+\mu +k{c}_1\right)\lambda +\mathit{ehg}-f\left(\mathit{g\phi }{c}_1+c_{2}h+\mathit{\phi \alpha }+{\delta }_2\right)+\mathit{f\epsilon \sigma }{c}_1{\phi }_1=0$

Assuming that $L_2=1$, $,L_1=l+\mu +k{c}_1,$ $L_0=-f\left(\mathit{g\phi }{c}_1+c_{2}h+\mathit{\phi \alpha }{\delta }_2\right)+\mathit{f\epsilon \sigma }{c}_1{\phi }_1$

we can apply Routh Hurwitz stability criteria since $L_2=1>0$ both, $L_1$ and $L_0$ should be positive.

Now, $L_1=l+\mu +k{c}_1>0$ and

Hence, the last two eigenvalues are negatives if $R_{\mathit{Eff}}<1.$

Thus, since all the eigenvalues are negatives, the disease-free equilibrium point of the HBV dynamical system (1)-(7) above is locally asymptotically stable if $R_{\mathit{Eff}}<1.$

Global stability of disease-free equilibrium point (DFEP)

Theorem 3: The disease-free equilibrium point of the dynamical system (1)-(7) which is given by:

$E^1=\left(S^0,V^0,E^0,A^0,C^0,T^0,R^0\right)=\left(\left(\frac{\Lambda +\pi }{\mu }-\frac{\mathit{\Lambda \psi }+\mathit{\theta \pi }}{\left(\tau +\mu \right)}\right),\frac{\mathit{\Lambda \psi }+\mathit{\theta \pi }}{\tau +\mu },0,0,0,0,0\right)$ of is globally asymptotically stable if $R_{\mathit{Eff}}\le 1$ otherwise unstable.

Proof: To show global stability of the DFEP, we applied Lyapunov function method as.^6,17

Let the Lyapunov function $L:R_{+}^7\to R_{+}$ is defined by:

$L=\left(S,V,E,A,C,T,R\right)={\alpha }_{1}S+{\alpha }_{2}E+{\alpha }_{3}A+{\alpha }_{4}C$. Then,

Taking the coefficients of $S,E,A,$ and$C$ are equal to zero, we get:

Partial derivative in terms of $S$ gives us: ${\alpha }_1\left(\left(1-\theta \right)\pi -\left(1-\psi \right)\Lambda -\mu \right)=0$.

This gives

Partial derivative in terms of $E$ gives us: $-\left(\mu +c_1+c_2\right){\alpha }_2+c_2{\alpha }_3+c_1{\alpha }_4=0$.

Partial derivative in terms of $A$ gives us: $-\frac{\mathit{c\omega }{S}^0}{N}{\alpha }_1+{\alpha }_2\frac{\mathit{c\omega }{S}^0}{N}-\left({\delta }_2\alpha +{\delta }_1\alpha +\mu \right){\alpha }_3+\alpha {\delta }_2{\alpha }_4=0.$

Partial derivative in terms of $C$ gives us: $\frac{\mathit{c\omega }{S}^0}{N}{\alpha }_1-{\alpha }_2\frac{\mathit{c\omega }{S}^0}{N}-\left(\mu +d+\sigma \right){\alpha }_4=0.$

Since $S\le S^0,$ and from equations (8)-(11) solving for ${\alpha }_2,{\alpha }_3$ and ${\alpha }_4$ and after simplification, we get:

Since $\frac{S^0}{N}=\frac{\tau +\mu -\mathit{\mu \theta }}{\tau +\mu }$, and $\frac{\mathit{c\omega }{S}^0}{N}=\mathit{c\omega }\left(\frac{\tau +\mu -\mathit{\mu \theta }}{\tau +\mu }\right)$,

$\Rightarrow \frac{\mathit{dL}}{\mathit{dt}}={\alpha }_2{c}_2\left(R_{\mathit{Eff}}-1\right)$, since ${\alpha }_2{c}_2>0,$

Then, $\frac{\mathit{dL}}{\mathit{dt}}={\alpha }_2{c}_2\left(R_{\mathit{Eff}}-1\right)\le 0$ where $R_{\mathit{Eff}}\le 1$ and $\frac{\mathit{dL}}{\mathit{dt}}=0$ if and only if $S=S^1,E=E^1,A=A^1$ and $C=C^1$. Therefore, the largest compact invariant set in $\left\{S,E,A,C\right\}\in {\Omega }_1:\frac{\mathit{dL}}{\mathit{dt}}=0$ is the singleton $\left\{E^1\right\}$ where $\left\{E^1\right\}$ is the disease-free equilibrium point of the model (1)-(7).

Thus, by LaSalle’s invariance principle,¹⁸ it implies that the disease-free equilibrium point $E^1=\left(\left(\frac{\Lambda +\pi }{\mu }-\frac{\mathit{\Lambda \psi }+\mathit{\theta \pi }}{\left(\tau +\mu \right)}\right),\frac{\mathit{\Lambda \psi }+\mathit{\theta \pi }}{\tau +\mu },0,0,0,0,0\right)$ is globally asymptotically stable in ${\Omega }_1$ if $R_{\mathit{Eff}}\le 1$ otherwise it is unstable.

Endemic equilibrium points (EEP), E₁

Let the endemic equilibrium of our HBV model system (1)-(7) be denoted by $E_1$ = $(V^{\ast },S^{\ast },E^{\ast },A^{\ast },C^{\ast },T^{\ast },R^{\ast })$. It is obtained by setting the right-hand side of each equation of our model (1)-(7) equal to zero and solving for the state variables in terms of the force of $\lambda =\mathit{c\omega }\left(\frac{\mathit{\phi C}+A}{N}\right)$. That is:

$E_1$ = $(V^{\ast },S^{\ast },E^{\ast },A^{\ast },C^{\ast },T^{\ast },R^{\ast }$)

Thus, after some calculation, we get the endemic equilibrium point is given by:$E_1=\left(V^{\ast },S^{\ast },E^{\ast },A^{\ast },C^{\ast },T^{\ast },R^{\ast }\right)$ where $V^{\ast }=\frac{\mathit{\theta \pi }}{\tau +\mu }$;

$R^{\ast }=\left(\frac{1}{1-a\left(\epsilon +\sigma \right)}\right)\left(a\left(1-\theta \right)\pi +\frac{\tau \left(\mathit{\theta \pi }+\psi \Lambda \right)}{\left(\tau +\mu \right)}\right)$, where:

Local stability analysis of the Endemic Equilibrium Point (EEP)

The disease-free equilibrium point of the model (1)-(7) above is locally asymptotically stable if $R_{\mathit{Eff}}<1$ and unstable if$R_{\mathit{Eff}}<1$.

Theorem 4: The endemic equilibrium point $E_1=\left(S^{\ast },V^{\ast },E^{\ast },A^{\ast },C^{\ast },T^{\ast },R^{\ast }\right)$ of the HBV model, (1)-(7) is locally asymptotically stable (LAS) if $R_{\mathit{eff}}>1.$

Proof: To show the local stability of the endemic equilibrium point, we use the method of the Jacobian matrix and Routh Hurwitz stability criteria. Then, the Jacobian matrix of the dynamical system (1)-(7) at the endemic equilibrium point $E_1.$ The Jacobean matrix, $J$($E_1$) of model (1)-(7) with respect to $E_1=\left(S^{\ast },V^{\ast },E^{\ast },A^{\ast },C^{\ast },T^{\ast },R^{\ast }\right)$ at the endemic equilibrium point is:

$a=\left(\tau +\mu \right),b=\left(\mu +c\omega \left(\frac{\phi C^{\ast }+A^{\ast }}{N}\right)\right),c=\left(\mu +c_1+c_2\right),d=\left({\delta }_1\alpha +{\delta }_2\alpha +\mu \right),f=\left(\mu +d+\sigma \right),g=\left({\phi }_1+\mu \right),h=\left(\gamma c\omega \left(\frac{\phi C^{\ast }+A^{\ast }}{N}\right)+\epsilon +\mu \right),j=\frac{c\omega S^{\ast }}{N},k=\frac{\mathrm{c\omega \phi }{S}^{\ast }}{N},l=\gamma c\omega \left(\frac{\phi C^{\ast }+A^{\ast }}{N}\right),m=\frac{\mathrm{c\omega \phi }{S}^{\ast }}{N},n=\frac{c\omega S^{\ast }}{N},p={\delta }_1\alpha ,r=\left(\gamma +\epsilon +\mu \right),u=c_1{E}^{\ast }$

Then, determinant of this Jacobian matrix is given by:

After simplification and adjustments, we get:

Then, we have: $\left(-a-\lambda \right)\left(-g-\lambda \right)\left(-r-\lambda \right)\left(-b-\lambda \right)\left(A_3{\lambda }^3+A_2{\lambda }^2+A_1\lambda +A_0\right)=0.$

Here, we get ${\lambda }_1=-a<0$,${\lambda }_2=-g<0$,${\lambda }_3=-r<0$ or ${\lambda }_4=-b<0$and$A_3{\lambda }^3+A_2{\lambda }^2+A_1\lambda +A_0=0.$

Where: