Numerical investigation for the nonlinear model of hepatitis-B virus with the existence of optimal solution

1.   Introduction

Mathematical modeling is an art of translating a physical phenomenon into a tractable mathematical formulation whose theoretical and numerical analysis provides a deeper about the phenomenon insight ^[1,2,3,4]. Mathematical modeling provides the precision and direction for problem solving and furnishes a thorough understanding of the real world system. These continuous models of infectious diseases are quite helpful in investigating and eradicating the infection in the society for a specific period of time ^[5,6,7]. There are some viral infections that can disrupt the economy of the nation, seriously ^[8,9]. Transportation is the basic need of the human beings. Some diseases communicate in the society by the process of diffusion. So, it becomes necessary to impose some precautionary restraints on the immigrants. Mathematicians, scientists and policy makers adopt different strategies for the welfare of the society. Mathematicians serve the community by developing and modifying the various mathematical models of the concerned problems ^{[10,11,12,13,14]}.

By analyzing such models, the communication of the various infections can be predicted more correctly ^[15,16,17]. In the current work, communication of the hepatitis-B virus is studied along with its physical properties. Hepatitis-B is a life threatening liver disease, it spreads by the hepatitis-B virus in the human body. If proper care and treatment are not provided to the infected person, he may become chronic. Chronic stage of hepatitis-B arises, if the disease lasts for more than six months. This stage may lead to the fatal diseases of liver such as cancer or cirrhosis. Infected individuals may get recovered fully, even if the infection is severe. Generally, the signs and symptoms develop within four months after being infected. Some viral diseases may be controlled by observing the precautionary measures depending upon the nature of the infection ^{[18,19,20,21]}. Some of the signs and symptoms of the hepatitis-B are upper abdominal pain, dark coloured urine, low fever, severe joint pain, appetite deficiency and fatigue etc. The hepatitis-B virus propagates from individual to individual through blood transfusion or by the body secretions and fluids for instance, semen, saliva etc. The most common modes of the HBV communication include the sexual contacts, infected syringes, blades and contaminated instruments such as razors, needles and forceps etc. One may become contaminated by doing unsafe sex with the polluted partner. Furthermore, this infection may be transferred from one person to the other by the use of poisoned needles, syringes, razors etc. Generally, this disease has two stages, that is, acute and chronic. At acute stage, the infection occurs less that twenty four weeks. Immune system of the body can fight against the infection and the infected person can recover within a few weeks. After a period of six months, the stage becomes chronic, this stage may not be curable. It remains throughout the life. Keeping in view the problems relating this disease, some specific models are devised to elaborate the disease dynamics. Various researchers took into account different necessary factors according to the disease spread. No one in this real world, can fulfil his daily life needs from one place. For this, he has to move from one place to the other. Now, if an infected individual travels widely in search of the community, the infection will diffuse in that population. To overcome this type of factors, an appropriate mathematical model will be required. The infection occurs mostly in Asia, Southern Europe, Africa and Latin America. This disease is a good example of epidemic in China, Asia and Africa. This virus can transmit both horizontally and vertically. The horizontal transmission occurs through the transfusion of the infected blood to a susceptible individual by other polluted fluids of the body, by the unsafe sexual activities or by the use of polluted syringes. On the other hand, the vertical type of transmission happens from mothers to infants during pregnancy. In this study, advection-reaction-diffusion model is taken into account for the investigation purpose of the hepatitis-B virus (HBV).

2.   Compartmental model of hepatitis-B virus

In this section, a mathematical model of hepatitis-B virus is {studied} for the numerical investigation. Let $\omega = (0, L)\times (0, T)\subseteq R^2$ represents the domain of the problem. Also, suppose that $s = s(x, t)$, $i_A = i_A(x, t)$, $i_C = i_C(x, t)$, $r = r(x, t)$ be the sufficiently smooth real valued functions defined on the whole domain $\Omega$. Let $h_1(x), h_2(x)$ and $h_3(x)$ be the continuously differentiable positive real valued functions defined on $(0, L)$. The proposed mathematical model of hepatitis-B virus with advection and diffusion terms, is defined as:

Since, $r(t)$ is absent in the first three equations so, above system reduces to the new system

The initial and boundary conditions, we will impose, for all $x\in(0, L)$

and

The proposed system $(2.1)$–$(2.4)$ describes the dynamics of the acute and the chronic hepatitis-B virus, where $s(x, t), i_A(x, t)$ and $i_C(x, t)$ are represented by the total number of susceptible individuals to the hepatitis-B virus, infected with the acute and the chronic hepatitis-B at any point $x$ and any time $t$ respectively. Also, $r(x, t)$ defined in Eq $(2.4)$, represents the number of recovered individuals. Due to some biological reasons, we assume that $s, i_A$ and $i_C$ be the positive functions ^{[11,12,22,23,24]}. The positive constants $d_1, d_2$ and $d_3$ serve as the spatial rate of diffusion in each of the population involved in the model $(2.5)$–$(2.7)$. In the prescribed model, $\delta_1$ represents the birth rate, $\alpha_c$ is the transfer rate of individuals from susceptible to infected class with the acute hepatitis-B, $r_1$ is the rate of moving from the acute stage to the chronic stage, $r_2$ is the recovery rate from the compartment $i_A$ to the compartment $r$, $r_3$ is the recovery rate from $i_C$ to the compartment $r$, $\mu_r$ is the natural death rate, $\delta_2$ is the death rate due to the hepatitis-B, $\nu$ represents the vaccinated rate of the hepatitis-B.

The disease free equilibrium of the system is

The endemic equilibrium of the model is:

and

where

The system $(2.1)$–$(2.4)$ is the modification of the system described in ^[13]. The local and global stability analysis is also described in ^[13]. Various authors considered the advection and diffusion in the nonlinear mathematical models which have their own physical importance in the dynamics of the system ^{[25,26,27,29]}. In the current article, advection and diffusive effects of the compartmental population of an epidemic model are addressed.

In the next section, the methodology of nonstandard finite difference scheme is described for the model $(2.1)$–$(2.4)$ with the initial and boundary conditions $(2.8)$–$(2.11)$.

2.1.   Optimal analysis of the model

The main purpose to work on this section is, to study the existence theory for the admissible solutions of the epidemic model for the chronic hepatitis-B. The solution tetroid, namely the susceptible, acute infected, chronic infected and recovered populations are given by the coupled system $(2.1)$–$(2.4)$, which involves the corresponding advection as well as the diffusion with constant rates. The above system $(2.1)$–$(2.4)$ is considered under the given physical conditions $(2.8)$–$(2.11)$. Physical conditions, here, are meant to be the nonnegativity. The classical solutions of the system of differential equations depend on the choice of the sufficiently smooth functions appearing in the underlying system. Under these assumptions, one may choose the fixed point results for analysis. The classical solutions generally, lie in the Banach space of continuous functions, say $C$. For the explicit results, we consider the following closed subsets in the space of the continuous functions $C([0, L]\times[0, T])$, as follows:

Considering the system $(2.1)$–$(2.4)$ with $(2.8)$–$(2.11)$ within the framework of fixed point theory with the following fixed point operators

Also, let

We apply the Schauder fixed point theorem to the above operator. We are interested in the second version of this theorem where we have to verify the following two mandatory conditions:

(i) $\aleph^i (X):B\rightarrow B$ (self map).

(ii) $\aleph^i(B)$ is relatively compact.

So,

This bound of $\digamma_i$ is obvious in case the advection and diffusion coefficients are continuous but in the current problem they are constants and hence continuous, which laid to the fact that $\digamma_i$ is continuous and must be locally bounded.

Now,

if

This implies

Assume that, the optimal value of $T$ be $T^*$, then

Obviously, it gives the maximum value of $T$ for which we can find the solution of the proposed problem and the largest (global) interval of the solution of the system $(2.1)$–$(2.4)$ is $[0, T^*]$.

$\textbf{(ii)}$ For the second condition, for all operators $T^iX$, i = 1, 2, 3, 4, we assume that the image sequence $T^i_jX$ for the pre-images of the components of $X$, that is,

where $p^i$ is a special number representing the desired solution vector $(S, I_A, I_c, R)$. This is the family of images $\aleph_j^i(x)$ of $p_j^i$ and we have to show that $\aleph_j^ix$ is equi-continuous. So, suppose, we have a specially chosen point $t^*$, such that

where $t^*\in [0, t]$

Since, $\digamma_i$ are the family of continuous functions hence, bounded and bounded by a number (say) $\mathcal{C}_i$, then

which shows that

This implies that $\aleph_j^i(x, t)$ are equicontinuous for each $i = 1, 2, 3, 4$ respectively. Consequently, all the sequences of the fixed point operators for each of the state variables turn out to be equicontinuous. Then by the Arzela-Ascoli, there exist uniformly convergent subsequences of each of $\aleph_j^i$, $i = 1, 2, 3, 4$, then each $\aleph^i$ is relatively compact, so, $\aleph^i(B)$ is relatively compact for a $\mathcal{C}_i(r)$, depending on the radius of the ball stated above. So, we formulate the following result.

Theorem 1. The problem $(2.1)$–$(2.4)$ with the initial conditions $(2.8)$–$(2.11)$ has the solution provided that the condition $(2.12)$ holds.

Due to a number of parameters in epidemological models, it is hard, some times impossible, to find the exact solutions of these models. So, numerical solutions become the necessity for the study of multi-parameterized systems. In the next section, we use the nonstandard finite difference scheme to approximate the solutions of the proposed model.

2.2.   Numerical modeling

Let, for any two natural numbers $M$ and $N$, $h = \frac{L}{M}$, $k = \frac{T}{N}$ be the positive constants. Consider the partitions of the intervals $[0, L]$ and $[0, T]$ with partition norms $h$ and $k$ respectively. Define $x_i = ih$ and $t_n = qk$ where $i\in \{0, 1, 2..., M\}$ and $n\in \{0, 1, 2..., N\}$. Also, suppose that $S_i^n, {I_A}_i^n$ and ${I_C}_i^n$ be the numerical approximations of the exact values of $s(x_i, t_n), I_A(x_i, t_n)$ and $i_C(x_i, t_n)$ respectively at the grid point $(ih, nk)$, $i\in \{0, 1, 2, ..., M\}$ and $n\in \{0, 1, 2, ..., N\}$. Furthermore, if, J represents any of the values $S, I_A$ and $I_C$, then, we have

The continuous system of the partial differential equations $(2.5)$–$(2.7)$ can be discretized by using the linear discrete operators.

where $i\in \{0, 1, 2..., M\}$ and $n\in \{0, 1, 2..., N\}$.

The linear operator defined in Eq (2.13) approximates the first partal derivative of $J$ with respect to $t$ at points $(x_i, t_n)$ and $(x_i, t_{n+1})$. The operators $(2.14)$ and $(2.15)$ approximate the first and second partial derivatives with respect to $x$ respectively at the same points $(x_i, t_n)$ and $(x_i, t_{n+1})$.

After using these operators in the system $(2.5)$–$(2.7)$, we get a discrete system of algebraic equations.

After simplifications, $(2.16)$–$(2.18)$ give

where $P = \frac{k}{h}$, $Q_1 = \frac{kd_1}{h^2}$, $Q_2 = \frac{kd_2}{h^2}$, $Q_3 = \frac{kd_3}{h^2}$ and $i\in \{1, 2, ..., M\}$, $n \in \{0, 1, 2, ..., N-1\}$.

The discretization of initial and boundary conditions will be

and

To make the comparison, we choose two other well known schemes. The first one is the upwind implicit scheme. After some simple calculations, on $(2.5)$–$(2.7)$, we can get

Also, consider an other renowned Crank-Nicolson type numerical scheme for the system $(2.5)$–$(2.7)$

3.   Physical properties

In this section, some important physical and numerical traits for the system $(2.16)$–$(2.18)$ will be derived. These properties play a vital role to approximate the solution of continuous system from a discretized one.

Definition 3.1 [Z-matrix] A real matrix $A$ is said to be a $Z$-matrix if all of its off-diagonal entries are non-positive.

Definition 3.2 [M-matrix] A square matrix $A$ over $R$ is said to be an $M$-matrix, if

(a) $A$ is a $Z$-matrix

(b) all the principal diagonal elements of $A$ are positive

(c) $A$ is strictly diagonally dominent.

$M$-matrix theory ^[30] is quite useful to prove the positivity of the mathematical models concerning engineering, economics, autocatalytic chemical reactions, etc. In the next result, positivity of the discrete system $(2.19)$–$(2.21)$ will be proved with the help of $M$-matrix theory. Note that, if, a matrix is an $M$-matrix, so, it is invertible and its inverse is positive matrix ^[28].

Theorem 2. Let $h_1, h_2$ and $h_3$ be the positive functions defined on $(0, L)$ then the system $(2.16)$–$(2.19)$ with initial and boundary conditions, is solvable for all $k > 0$, $h > 0$ and the solutions are positive.

Proof. As we know that, the vector form of the system $(2.19)$–$(2.21)$ can be written as

where $F, G$ and $H$ are the square matrices of order $(M+1)$ then

and

where

and

Now, to prove the positivity of the system $(2.16)$–$(2.18)$, we adopt the technique of mathematical induction. Since, $S^0, I_A^0$ and $I_C^0$ are positive component vectors defined in the initial data, so, we suppose that $S^n, I_A^n$ and $I_C^n$, $(n\in 0, 1, 2, \cdot \cdot \cdot, N-1)$, are positive vectors with positive components. By the above calculations, it is shown that $F, G$ and $H$ are $M$-matrices, so, they are invertible and have positive inverses. Moreover, right hand side in the system $(3.1)$–$(3.3)$ is positive. Therefore,

all are positive vectors.

Hence, by the principle of mathematical induction, we have concluded the result.

Definition 3.3. Let $M_h = \bigg\{x_i\in R: i\in \{0, 1, 2, \cdot\cdot \cdot, M \}\bigg\}$ be the set of grid points and $S_h$ be the vector space of real valued functions defined on $M_h$. Define a norm $\|.\|$ such that

and

An other important structural property of our proposed numerical scheme is, the consistency. For this purpose, define the differential operators

Also, the discrete differential operator be defined as

3.1.   Consistency of proposed scheme

By Taylor's series expansion, we can get the accuracy of proposed scheme. For this, let

After using Taylor's series expansions and some easy calculations, we get

and

Also,

Hence, the proposed numerical scheme is consistant with the system $(2.5)$–$(2.7)$.

Using definition $3.3$ and Eqs $(3.4)$–$(3.9)$, we can establish a result.

Theorem 3. If $s, i_A, i_C, r \in C^{2, 2}_{x, t}(\bar{\Omega})$, then there exists a non-negative constant $\xi$, independent of $h$ and $k$ such that $\max \{\|\theta-\theta'\|_\infty, \|\phi-\phi'\|_\infty, \|\psi-\psi'\|_\infty\}$.

3.2.   Stability

The main concern in the study of approximating the solutions of the system of differential equations is, the growth of round off errors in the numerical solutions. An other main thing is important to observe, is that, a small change in the initial conditions may cause a large deviation in the solution of the underlying system. In this scenario, if, the slight change in the initial data does not produce a huge variation in the approximate and exact solutions, we say, the numerical scheme which gives such approximate solutions, is stable. To discuss the stability analysis of the proposed scheme, we use the method of Von-Nuemann. For this purpose, we decompose the numerical error occurred in the numerical solutions in to the Fourier series. So, linearizing the Eqs $(2.19)$–$(2.21)$ and substitute

we get

By putting

we have

Similarly, from $(2.21)$, we have

Hence, the proposed scheme is Von-Neumann stable.

4.   Numerical experiment

4.1.   Example 1

The initial conditions are

The values of parameters involved in the model are $\delta_1 = 0.4$, $r_1 = 0.01$, $\mu_r = 0.03$, $\delta_2 = 0.002$, $r_2 = 0.05$, $r_3 = 0.06$, $\nu = 0.02$. For endemic equilibrium or steady state, we use $\alpha_C = 5$; and for disease free equilibrium or steady state $\alpha_C = 0.005$ and $d_1 = d_2 = d_3 = 0.02$.

Now, we give the numerical simulations to validate numerical results. In the infectious disease dynamics, the steady states of the epidemic model usually show the stable behavior with reference to the basic reproductive number $(R_0)$. As we mentioned before, the epidemic models of the hepatitis-B virus including advection and diffusion have two steady states, one is infection free when $(R_0 < 1)$ and second is infection existence state when $(R_0 > 1)$. Now, this steady state behavior must be shown graphically by the numerical methods.

In the Figure 1, the graphical behavior is elucidated at endemic state by using well-known classical techniques. The mesh solution graphs of susceptible individuals are depicted with the aid of upwind implicit technique and Crank Nicolson technique. It can be observed that both techniques demonstrate the negative solution of population dynamics which is not the part of the continuous advection diffusion reaction system. Moreover, the combined plots of $S(x, t)$, $I_A(x, t)$ and $I_C(x, t)$ describe that these classical schemes diverge for different step sizes.

The graphs of $s(x, t)$, $I_A(x, t)$ and $I_C(x, t)$ are plotted with the help of proposed scheme as shown in Figure 2. The parametric values and step sizes are the same as in the Figure 1. The upwind and Crank Nicolson methods fail to provide the positive solution. Moreover, these schemes diverge for slightly greater values of some parameters as described in the 1. On the other side, our scheme submit the positive solutions for every value of the parameters as mentioned earlier. So, it is apposite to say that our scheme is more reliable than that of the famous upwind and Crank Nicolson schemes.

Now, we present the numerical simulations of Hepatitis infection model by using the proposed technique in the Figure 3 and it is observed that this method not only sustains the positivity of the solution but also the stability of the endemic state as well. It is clear from the solution graphs that the proposed NSFD technique holds all the important structure of hepatitis-B continuous model.

Figure 4 exhibits the numerical simulations of all the sub-populations of the proposed epidemic advection diffusion system at disease free state with the assistance of the designed upwind NSFD technique. It is verified from the Figure 4 that NSFD technique reveals the positivity of the solution and also retains the stability of disease free state of continuous model. From all the figures presented in this manuscript, it is concluded that the proposed technique is dynamically consistent with the hepatitis-B continuous system under study.

4.2.   Example 2

In this example, we consider the continuous initial conditions for the underlying model and observed the behavior of state variable as describe in Figures 5 and 6. The values of parameters involved in for the model are $\delta_1 = 0.4$, $r_1 = 0.01$, $\mu_r = 0.03$, $\delta_2 = 0.002$, $r_2 = 0.05$, $r_3 = 0.06$, $\nu = 0.02$. For endemic state, we used $\alpha_C = 5$; and for disease free state $\alpha_C = 0.005$ and $d_1 = d_2 = d_3 = 0.02$.

The initial conditions are $S(x, 0) = 1+0.01\sin(\pi x)$, $I_A(x, 0) = 1-0.12\sin(\pi x)$ and $I_C(x, 0) = 1/3$ with homogeneous Neumann boundary conditions.

All the three graphs in the Figure 5 illustrate the dynamical behavior of the state variables $S(x, t), I_A (x, t)$ and $I_C (x, t)$ at disease free fixed point by using the proposed NSFD technique. It is evident from the numerical designs that the nonlinear system attains the stable state. Moreover, the disease dies out at this stage and whole the population become zero. These important biological facts are demonstrated by the graphs in the Figure 5.

Hence our projected numerical schemes converge towards the exact disease-free state. Numerical graphs in the Figure 6. Show the dynamical behavior of the state variables at endemic equilibrium position. Non-standard finite difference numerical scheme is applied to reach at the numerical solutions. When the disease exists in the domain population, the susceptible and infected population have some positive sizes. The Figure 6(a–c) ascertains this biological fact, by the converging towards the correct endemic steady state. From the Figure 6(a–c) it is evident that our numerical scheme converges towards the true fixed point. Hence, the prescribed NSFD device is a reliable numerical tool to solve no linear dynamical systems.

Remark 1. The underlying model is about the infectious disease chronic Hepatitis B virus. It has two steady states, namely the disease free steady state and endemic equilibrium. Basic reproductive number $R_{0}$ plays a decisive role in describing the steady states.The figures in the numerical experiment show that disease free steady state is attainable, when $R_{0} < 1$. Endemic equilibrium is attained, when $R_{0} > 1$. Moreover, the numerical graphs illustrates how the proposed numerical design help to attain the steady states for the plausible set of parametric values. Also the dynamical behavior of the state variables, to attain the equilibrium point can be observed by the graph.

5.   Conclusions

In the current article, we have studied, numerically, an advection-reaction-diffusion, nonlinear epidemic model of chronic hepatitis-B model. The concept of the optimal existence of the solutions for the described model is established. In the scenario of optimization, we have developed the results regarding the existence and uniqueness of the solutions of the system. The auxiliary data, initial and boundary conditions, are investigated. As the solutions of the governing equations lie in the spaces of functions but it is better to consider the subsets (closed ball) of a function space. A closed ball is chosen for this purpose and the explicit estimates are found. This closed ball is called the optimal ball. Under the certain conditions together with the explicit estimates, existence of the solutions of the physical model is proved via Shauder fixed point theorem. Addition of advection and diffusion in the considered set of equations, make it more general and realistic. After the analysis of the system, a numerical solution is computed by the nonstandard finite difference method. By using the approximations of the spatial and temporal derivatives, a corresponding discrete model is generated. The prominent feature of the numerical technique used for the proposed model, is that, it is structure preserving, that is, the discrete system obtained from the numerical scheme has the same properties that the corresponding continuous system possessed. We proved that the developed scheme is consistent with the proposed model. Stability of the numerical scheme is verified by Von Nuemann criteria. An other important feature is the positivity of the values of the state variables involving in the system under study, so, by using the $M$-matrix theory, positivity of the concerned model is proved. Finally a test problem is given to verify the results with the help of numerical simulations. We can extend our work to two and three dimensions. Hence, the proposed NSFD scheme is a value addition in the numerical solutions of the nonlinear infectious disease model.

Acknowledgment

The authors wish to thank the anonymous reviewers for their comments and criticisms. All of their comments were taken into account in the revised version of the paper, resulting in a substantial improvement with respect to the original submission.