Approximate solutions for HBV infection with stability analysis using LHAM during antiviral therapy.

Hepatitis B virus (HBV) is a life-threatening virus that causes very serious liver-related diseases from the family of Hepadnaviridae having very rare qualities resembling retroviruses. In this paper, we analyze the effect of antiviral therapy through mathematical modeling by using Liao's homotopy analysis method (LHAM) that defines the connection between the target liver cells and the HBV. We also examine the basic nonlinear differential equation by LHAM to get a semi-analytical solution. This can be a very straight and direct method which provides the appropriate solution. Moreover, the local and global stability analysis of disease-free and endemic equilibrium is done using Lyapunov function. Mathematica 12 software is used to find out the solutions and graphical representations. We also discuss the numerical simulations up to sixth-order approximation and error analysis using the same software.

Introduction

Hepatitis B virus (HBV) infection is one among the very serious infections that threatens the people worldwide, less than 300 million people live with chronic hepatitis B infection worldwide. It directly or indirectly affects the liver, causing diseases affecting also infants and young children. This virus is genetically from the family of Hepadnaviridae having very rare qualities resembling retroviruses [1]. Hepatitis B is a viral infection that attacks the liver and can cause both acute and chronic diseases. The infection is a potentially life-threatening liver infection caused by the hepatitis B virus (HBV). In 2015, hepatitis B resulted in an estimated 887,000 deaths, mostly from cirrhosis and hepatocellular carcinoma. As of 2016, 27 million people (10.5% of all people estimated to be living with hepatitis B) were aware of their infection, while 4.5 million (16.7%) of the people diagnosed were on treatment. Hepatitis B can be prevented by vaccines that are safe, available, and effective. The modalities of transmission of the HVB are completely different from those of the HAV (hepatitis A virus) but can be compared to those of the HCV and HIV virus. HCV is most commonly transmitted through contact with blood or other body fluids as well as from a positive mother to a child during birth and delivery. This poses a major public health problem for the protection of those living with affected patients, healthcare professionals and requires serious and effective control of blood transfusions. Transmission of the virus may also occur through the reuse of needles and syringes either in healthcare settings or among persons who inject drugs [2-4]. In addition, infection can occur during medical, surgical, and dental procedures, through tattooing, or through the use of razors and similar objects that are contaminated with infected blood. HBV infection can be either acute or chronic. The illness can be asymptomatic, symptomatic and may progress to cirrhosis. There is no specific treatment for acute hepatitis B. Chronic HBV infection is defined as persistence of hepatitis B surface antigen for six months and more after acute infection with HBV. It is a major public health problem because the majority of people are unaware of their HBV infection; they are asymptomatic and can transmit the disease. Chronic hepatitis B infection can be treated with medicines, including oral antiviral agents. Only a proportion (estimates vary from 10% to 40% depending on the setting and eligibility criteria) of people with chronic hepatitis B infection will require treatment. So, to prevent this disease and infection, hepatitis B vaccine should be taken. The vaccine is effective to prevent the infection increasing the protective levels of antibodies [5, 6]. Mathematical model is an essential tool to understand the virus dynamics, and it intensifies our understanding of the virus. The first important mathematical models covering immune response were developed by Anderson and May on infectious diseases of humans, in particular HIV and other viral diseases [7]. Other authors proposed different mathematical models concerning the transmission of viral infection. Blessing proposed a mathematical model of hepatitis B virus transmission dynamics with considerations of different classes of individuals, namely immunized, susceptible, latent, infected, and recovered classes [8]. Kamyad et al. suggested a mathematical model considering the immune response of vaccinated infants, what treatments given to an infected person, and how it controls the transmission of HBV [9]. The humoral immune response is also important in HCV infection and is the key to vaccine development [10]. Su, B. et al. [11] proposed an HBV transmission model utilized to know the vaccination effectiveness and transmission control.

Application of the mathematical models to different parameters shows that some patients get cleared of the virus quickly due to drug therapy, whereas in case of other patients the drug therapy works much slower [12]. Therefore, in some cases treatment has to be stopped. Then, the viral degeneration is observed in most patients [13]. In this paper, we try to get a deeper understanding of the complex virus dynamics that is seen during drug therapy in HBV infections. Such knowledge helps to improve the treatment by informing what drugs to use, when to start, and when to stop [14-17]. The main motivation of the paper is that this will help the young researchers of medicine as well as science find a solution comparing to real world problems pertaining to HBV antiviral therapy. Finding the analytical solutions for this model is very challenging. But we found the analytical solution for Eq. (2.1) using the LHAM method. Our model is very useful and it is easy to find the analytical solutions. We also can find the numerical simulations by using MATLAB for the same equations.

This paper is divided into eight parts. The first part of the paper forms an introduction which deals with the existing literature and proposed work. The second and third parts talk about mathematical modeling and analytical study of the HBV model. In these parts we analyze the local and global stability of HBV. The fourth part explains the LHAM method, and the fifth part is applications that have been used to find out the solutions using LHAM. The sixth part of the paper contains numerical experiments which deal with numerical simulations obtained up to sixth-order approximations using Mathematica 12 software. Error analysis forms the seventh part of the paper, the eighth part is discussion, and final part of the paper is conclusion.

Mathematical modeling

We consider the mathematical model for the basic virus dynamics [18]. This model represents an in-host model for the interaction between liver cells (uninfected and infected) and the virus, which is derived from the mathematical model. Here, X represents the target uninfected cells (uninfected hepatocytes), Y represents the infected cells (infected hepatocytes), and Z represents the HBV virus. This model represents the target cells that are infected at rate β, and infection occurs by contact between the target cells and the virus. It also considers that the infected cells die at a rate of δ and production of new virus happens at a rate of v. The rate of dead virus is c. The new target cells introduced in the liver may die before getting infected. The constant production rate is denoted by s, and the natural death rate is denoted by . Parameters and variables are presented in Table 1 [19]. The variables are considered either horizontal (through contact with infective individuals) or vertical transmission (directly from the parents). The horizontal transmission is less efficient, but the vertical transmission occurs with great frequency (70% to 100%) when the mother has acute hepatitis B during delivery.

Table 1

Values of the parameter model

Parameter	Explanation	Values
s	production of constant hepatocyte	1 cell day−1 mL−1
d_T	death rate of hepatocyte	0.01–9 × 10−4 day−1
β	rate of infectivity	1 × 10−10–6.6 × 10−8 mL virions−1 day−1
δ	infected hepatocytes killing rate	0.06–0.25 day−1
ρ	rate of recovery	0 day−1
v	virus production	1.4–164 virions cell−1 day−1
c	virus clearance	0.18–1 day−1
η	effect of therapy in blocking infection	0.2–0.5
ε	effect of drug in blocking new virus production	0.99934–1

These assumptions lead to the following model: The initial and boundary conditions of finding the solution of Eq. (2.1) are

We consider three cells as the target cells, the infected cells, and the hepatitis B virus cells [20]. This model defines the interaction between the target liver cells and the HBV [21]. The HBV gets cleared for strong immune response that cures cells and makes them not to get reinfected. All parameter estimations are obtained in Table 1. Also, it represents the efficacy of the treatment on blocking the production of virion [9], and compartmental diagrams are given in Fig. 1.

Figure 1

Diagram for basic virus dynamics model

Analytical study of the HBV model

Boundedness and positivity

Theorem 3.1

Let. Ifthen ∀tin, then the solution is.

To prove this for , , and are positive. Therefore, the solutions of model (2.1) are positive.

Theorem 3.2

Let, where. The solutions of system (2.1) are bounded in a compact set. Obviously, the solutions of (2.1) are bounded in the region A.

Equilibrium analysis

Disease-free equilibrium and endemic equilibrium

Basic reproduction number

Therefore, Hence, by [22],

When the virus vanishes, otherwise it continues. Since the antiviral therapy is given, will never be zero.

Stability analysis

Proposition 3.3

is locally asymptotically stable forand unstable otherwise.

Proof

Therefore the eigenvalues of the disease-free equilibrium are , , and −c.

Eigenvalues at disease-free equilibrium are all negative, and hence it is stable. Therefore if , is locally asymptotically stable. □

Proposition 3.4

is globally asymptotically stable forand unstable otherwise.

Let the Lyapunov function be defined as

Differentiating (3.1) with respect to “t”, we have since which implies . Therefore only when in (2.1) such that and as . Hence, by [23], is globally asymptotically stable when .

Stability analysis has been discussed in this study, which gives the stable equilibrium points obtained from the characteristic equation systems of differential equations. The type of an equilibrium point is determined by the eigenvalues of the Jacobian matrix. Stability analysis is obtained by determining the eigenvalues of the Jacobian matrix around equilibrium points. □

Proposition 3.5

is locally asymptotically stable forand unstable otherwise. The characteristic equation of equilibriumiswhereBy the Routh–Hurwitz criterion, if

Finally, is locally asymptotically stable.

Proposition 3.6

If, thenis globally asymptotically stable and unstable otherwise.

Consider the Lyapunov function Differentiating (3.2) with respect to “t”, we have From (3.3) if then . Also iff , , .

Hence, by [23], is globally asymptotically stable when . □

Liao’s homotopy analysis method (LHAM)

Let us consider the equation

We get the following zero-order deformation equation from [22-26]:

Here, ℑ is a supplementary linear operator such that for integral constants ().

When and , the zero-order deformation equation can be written as follows:

Using the Taylor series expansion of with respect to y, we get where

Differentiating the equation with respect to y by i times, then setting , and lastly dividing them by i!, we get the ith-order deformation equations: where and

Applications

The solution of Eq. (2.1) is defined by using the LHAM method as follows:

To obtain the analytical solution, the homotopy is

Equating terms, we get

From Eq. (5.7) ⇒ , where and .

From Eq. (5.8) ⇒ .

From Eq. (5.9) ⇒ , where and .

Again, equating terms, we get From Eq. (5.10) ⇒ where ; and .

From Eq. (5.11) ⇒ where ; and .

From Eq. (5.12) ⇒ The analytical solution of this model using LHAM is where , , , , . where , , . where , .

Numerical results

Let us consider the values for numerical results as follows:

Let us use Mathematica software to obtain the sixth-order expansions for , , and :

Error analysis

In this paper, an error analysis has been done to obtain the optimal values of h. It is shown in Figs. 8 to 10. The optimum and minimum values of h are found from Figs. 11 to 13. For that, we substitute Eqs. (5.1) to (5.3) in (2.1) and get the residual functions, which are shown in what follows. The h value ranges are given in Table 2, and the minimum values are shown in Table 3. Also, the residual errors are calculated in Table 4.

Figure 8

The residual error function of Eq. (7.1)

Figure 10

The residual error function of Eq. (7.3)

Figure 11

The optimum and minimum values of

Figure 13

The optimum and minimum values of

Table 2

The h value is

X(t)	−1.3 ≤ h ≤ −0.5
Y(t)	−1.5 ≤ h ≤ −0.9
Z(t)	−1.5 ≤ h ≤ −0.9

Table 3

The minimum values of , ,

	h^∗	Minimum value
RX(h_1)	−0.877677	3.332434 × 10−6
RY(h_2)	−0.649126	2.425675 × 10−8
RZ(h_3)	−0.514657	1.57936 × 10−12

Table 4

The residual errors for , , and for

t	ER_1(X,Y,Z;h1∗)	ER_2(X,Y,Z;h2∗)	ER_3(X,Y,Z;h3∗)
0.0	7.677888 × 10−6	5.536737 × 10−5	3.646747 × 10−8
0.1	3.647548 × 10−3	9.745774 × 10−6	8.747839 × 10−9
0.2	5.566774 × 10−2	2.747899 × 10−4	4.889889 × 10−6
0.3	1.647689 × 10−5	6.738893 × 10−3	9.374839 × 10−7
0.4	9.673487 × 10−9	1.988787 × 10−4	7.838439 × 10−6
0.5	6.674639 × 10−8	7.889987 × 10−6	2.473748 × 10−5
0.6	9.782357 × 10−7	3.774898 × 10−5	1.838493 × 10−7
0.7	6.646388 × 10−1	8.738783 × 10−4	5.737888 × 10−6
0.8	2.663879 × 10−6	1.789888 × 10−6	6.737882 × 10−4
0.9	8.653467 × 10−4	7.893289 × 10−7	7.374378 × 10−9
1	4.778898 × 10−5	4.789897 × 10−6	3.747889 × 10−5

Let us consider the square residual error for sixth-order approximation [27, 28]:

The minimal values of , , are as follows:

We consider the optimal values of , , and for all of the cases:

Discussion

We found the solution from (5.13)–(5.15), which contain ‘h’ that shows an easy technique to control and adjust curves to confirm series solution to converge, which is recommended by Liao [25, 26]. Figures 2–7 show the plots of sixth and seventh term approximation of , , and . Through these curves, it is clear that the valid region of ‘h’ is parallel to the horizontal axis. The valid region is listed in Table 2. Figures 8–13 show the residual error function of Eqs. (7.1), (7.2), and (7.3) using the sixth-order approximate solution for the different values of , , and .

Figure 2

The h-curves of sixth-order approximations for

Figure 7

The h-curves of seventh-order approximations for

The residual error function of Eq. (7.2)

Conclusion

Hepatitis B virus is proved to be a dangerous disease affecting people enormously. Though there has been vaccine to cure this disease, the antiviral therapy is recognized as the best method to eradicate this disease utmost at the root level. In this paper, the reproduction number is used to find the global dynamics. If , the disease-free equilibrium is globally asymptotically stable. Furthermore, if , the endemic equilibrium is globally asymptotically stable. The potentiality of HAM depicts the convergence of sequence solution for nonlinear differential equations, which we proved in this paper confirming that HAM is a very effective and powerful technique to find the approximate semi-analytical solutions. The numerical simulations have been obtained up to sixth-order approximations, and error analysis has been done with the help of Mathematica software. The study of mathematical models of disease development will allow better knowledge of disease evolution to reduce the incidence of accidental infections among healthcare professionals and to improve the quality of life of patients who may be given therapies already experienced in other hepatitis [29]. This research paper can be a framework for the young researchers to do a further research and design an effective antiviral therapy and drug design.