Numerical solution of diffusive HBV model in a fractional medium.

Evolution systems containing fractional derivatives can result to suitable mathematical models for describing better and important physical phenomena. In this paper, we consider a multi-components nonlinear fractional-in-space reaction-diffusion equations consisting of an improved deterministic model which describe the spread of hepatitis B virus disease in areas of high endemic communities. The model is analyzed. We give some useful biological results to show that the disease-free equilibrium is both locally and globally asymptotically stable when the basic reproduction number is less than unity. Our findings of this paper strongly recommend a combination of effective treatment and vaccination as a good control measure, is important to record the success of HBV disease control through a careful choice of parameters. Some simulation results are presented to support the analytical findings.

Background

In this paper, we consider in one dimensional space the multicomponent fractional reaction–diffusion system. The fractional derivative equation is obtained by replacing the second-order derivatives in classical n-variable reaction–diffusion systems, with order α in the interval 0 < α < 2. A coupled system of n (for n ≥ 1, n integer) species (population densities or concentrations of chemicals) which interact in a nonlinear fashion and diffuse may be modelled by the n system of equationswith the initial conditionSubject to any of the boundary conditions:

In the case of an infinite system, x ∈ (−∞,∞), here R is a subset of (−∞,∞).

no-flux or Neumann boundary condition for a finite system, and

called the Dirichlet or fixed concentration boundary condition, also for a fixed system where .

The fractional derivative operator(with similar expressions for u, i = 1, 2, … , n) is the Riemann–Liouville fractional gradient of order α. The diffusion matrix is defined as a diagonal matrix D = diag(D1, D2, … , D), and the diffusion coefficients D which do not depend on u must be positive. u(x,t) is a community of species density. Function F accounts for the local kinetics of the system.

For the purpose of simplicity, we may write,where for . The reaction kinetics in (1) is expected to satisfy the following conditions:

In  = {(u1, u2, … , u)|u > 0|}, the vector field (F1(u), F2(u), … , F(u)) has unstable equilibrium state at 0 = (0, 0, … , 0) and asymptotically stable at A = (A1, A2, … , A) with A > 0 for i = 1, 2, … , n.

The coefficientsmust be finite.

Fractional differential equations known as the differentiation and integration of noninteger order of the form (1) are becoming increasingly used as a modelling tool for diffusive processes associated with sub-diffusion, diffusion and super-diffusion scenarios, and have been applied to an increasing number of situations in biochemistry (Yuste et al. 2004), control (Wang and Zhou 2011; Wang et al. 2012), medicine (Erturk et al. 2011; Hall and Barrick 2008; Otte et al. 2016) and mathematical biology or physics (Atangana 2015; Barkari et al. 2000; Tomovski et al. 2012; Wang and Du 2013). In literature, many definitions of fractional derivatives have been given among which are the definitions by Caputo, Grnwald–Letnikov, Hadamard, Marchaud and Riemann–Liouville are regarded as the most useful definitions (Haghighi et al. 2014; Podlubny et al. 2009; Polyanin and Zaitsev 2004; Yang et al. 2010; Zeng et al. 2014, 2015; Zheng et al. 2015; Zhou 2014). Quite Unfortunate that, similar usages of different definitions give rise to different results (Podlubny 1999).

Over the years, several numerical and analytical methods of solution have been adopted to solve both linear and nonlinear equations. Among which are the homotopy analysis method (Alomari et al. 2009), homotopy perturbation method (He 2005; Yildirim and Sezer 2010), Adomian decomposition method (Hassan 2008; Ray 2009), multistep differential transform method (Li Ding and Lin-Jiang 2013; Erturk et al. 2011; Jiang et al. 2012), Kansa method (Chen et al. 2010), finite difference and tau methods (Celik and Duman 2012; Pang and Sun 2012; Saadatmandi and Dehghan 2011; Sausa 2009; Su et al. 2010) to mention a few. In the context of this paper, the spatial complexity of the domain imposes geometric constraints on the transport processes on all length scales, which can be measured as temporal correlations on all the time scales. Hence, we propose the study of fractional hepatitis B Virus (HBV) reaction–diffusion system in sub-diffusive (0 < α < 1), diffusive (α = 1) and super-diffusive (1 < α < 2) scenarios, using the Fourier spectral method in space to remove the stiffness issue associated with the spatial fractional derivative in a finite but large domain size L. For the temporal discretization, we employ higher-order exponential time differencing scheme.

In “Basic definitions and numerical techniques for fractional diffusion problems” section, we provide some required basic information and definitions of fractional calculus, we also formulate Fourier spectral method for multicomponents fractional-in-space reaction-diffusion system. We present fractional HBV model, and examine the main system for both local and global stability in “Main model” section. Numerical experiment of the full nontrivial result is provided in “Numerical simulations” section. The final section concludes the paper.

Basic definitions and numerical techniques for fractional diffusion problems

Two important tasks are done in this section. The first one is to present some of the basic required information and definitions that guide the solution of fractional diffusion equations. The second task is centered on the introduction of Fourier spectral methods as an efficient and easy-to-implement for the integration of fractional-in-space reaction–diffusion systems.

Basics definitions of fractional calculus

Definition 1

If f(x) ∈ [a, b] and a < x < b, then the Riemann–Liouville fractional integral is given as:

Definition 2

For α ∈ (0, 1),is known as the Riemann–Liouville fractional integral of order α (Pindza and Owolabi 2016; Saadatmandi and Dehghan 2011; Zeng et al. 2014; Zheng et al. 2015) where

Definition 3

For −1 < α < n, n ∈ N, x > 0, and , is the Caputo fractional derivative of order α is defined as:

Definition 4

For n − 1 < α < n (Haghighi et al. 2014; Ortigueira 2006; Samko et al. 1993)Is the Riesz fractional derivative, where and

Primarily, the Caputo fractional differential equation computes an ordinary differential equation (ODE), followed by a fractional integral to obtain the desired order of fractional derivative, but in the reverse order we compute the Riemann–Liouville fractional differential system. The Caputo fractional differential equation permits the inclusion of traditional initial and boundary conditions in the problem formulation. These two operators coincide in the case of homogeneous initial condition. For details geometric and physical interpretation for fractional calculus for both the Riemann–Liouville and Caputo types, readers are referred to the classical books (Kilbas et al. 2005; Podlubny 1999).

Theorem 1

Consider a one-component fractional reaction–diffusion modelwhere μ, t > 0, x ∈ R, α, η,, ν are real and positive parameters subject to the constraints 0 < α ≤ 2, |ϑ| < min(α, 2 − α), andis the generalized Riemann–Liouville fractional derivative operator given bywith the conditionswhere 1 < η ≤ 2, 0 ≤ ν ≤ 1, ρ > 0 is a constant with the kinetic (or reaction) term.is called the Riesz–Feller space fractional derivative of order α, skewness ϑ is symmetry given in terms of the Fourier transformwhere, 0 < α ≤ 2, |ϑ|, min(α, 2 − α), μ is the diffusion coefficient and(x, t) is a term that belongs to the area of reaction–diffusion function. The solution of (4) govern by the above conditions is

Proof

For the time variable, we apply Laplace transform w.r.t. time t, also for the space variable, we find the Fourier transform w.r.t. x using Laplace parameter s and Fourier parameter k, subject to the initial conditions and Eq. (6) with the Laplace transform of (5) given asfor n − 1 < η ≤ n, n ∈ u, 0 ≤ ν ≤ 1, we then transform the given equation into the formwhere * is the Fourier transform w.r.t. x with Laplace and Fourier parameters s and k. By convention we use (.) to denote the Laplace transform of (). So, solving for we obtainWe take the inverse Laplace transform of Eq. (8) and following the result of Saxena et al. (2006), we havewhere the real part of the parameters are positive, bear in mind thatwhich by inverse Fourier transform, Eq. (9) becomesas the required result. □

Numerical techniques for fractional reaction–diffusion equations

Spectral algorithms are considered to be extremely valuable for generating numerical schemes in almost every areas of applied mathematics, because it has no excuse to differentiate poorly but possesses the ability to yield spectrally accurate fractional or spatial derivatives. We can always guarantee of getting good and efficient working codes, when spectral methods is coupled with Fast Fourier transforms and implement with high-level language like Matlab. In one component, we use the integrating factor technique to Fourier transform of (1) to getwhere U is the Fourier transform of species u(x, t). This implies that,

We let to explicitly remove the issue of stiffness in the fractional derivative term, and set , so that

Now, in practice, we discretize spatial domain by considering N the equispaced points direction of x. We employ the discrete fast Fourier transform (DFFT; Owolabi and Patidar 2015) so that (12) becomes a system of ODEswhere u = u(x) and . The boundary conditions can be clamped at extremes of the domain. At this point, we have transformed the system to ODEs, more importantly, the spatial derivative and the associated stiffness are gone. We can employ any higher-order explicit solver.

By using the standard notation of order s Runge–Kutta (RK) schemes, with time step Δt, we can advance from t = nΔt to t = (n + 1)Δt for the ODEwhere

Here, we apply the general explicit Runge–Kutta scheme to (13) for . For brevity, we denote μ as k′ in and set the replacement variable as

Finally, we write the s-stage RK scheme aswith modified termand U values at intermediate step

This implies that we work entirely in the spectral domain and invert a transform to recover u. It should be mentioned that once the stiffness is removed, one can rapidly and accurately advance in time with any explicit higher-order time-solvers, see for instance (Kassam and Trefethen 2005; Owolabi and Patidar 2014) for details.

For the second time stepping solver, we formulate and utilize the modified Krogstad (2005) version of the Cox and Matthews scheme (Cox and Matthews 2002), which we denoted here for brevity as ETDRK4.with the stages μ given as

The functions φ, are given aswhere L =  is the fractional derivative operator with order α, and N is the Fourier transform that accounts for the nonlinear reaction functions. Readers are referred to Cox and Matthews (2002) and Owolabi and Patidar (2016a, b) for details derivation and stability of the family of exponential time differencing schemes.

At this point, we need to check the performance of the fractional Fourier transform technique in conjunction with both the classical fourth-order Runge–Kutta (Owolabi and Patidar 2014) and the fourth-order exponential time-differencing Runge–Kutta schemes. Here we consider (1) in one component and set the reaction term to zero, so that the given multicomponents system reduces to fractional diffusion equation.

In Fig. 1, we report the L-norm error defined aswhere and u are the approximate and exact solutions, and N represents the computational grids number, on x ∈ [0, 0.5] at t = 1 and initial condition u0(x) = exp[− 10x2/(1 − x2)]. The reference solution is taken by evaluating the fractional diffusion equation using 212 Fourier modes. It is proper to say that regardless of fractional power α, our approach (Fourier transform) is able to attain a spectral convergence up to machine precision.

Fig. 1

a Comparison of the numerical and analytical solutions of the fractional diffusion equation at α = 1.5, D = 0.5, obtained at instances of t = 2 (dots), t = 3 (circles) and t = 4 (crosses). The thick lines represent the analytical solutions. b Convergence of the fractional Fourier method at instances of α at final time t = 1

Timing results (independent of fractional power α) are given in Table 1 for comparison between both schemes. Indication here is that, the ETDRK4 displays a better performance when used in conjunction with the Fourier spectral method, especially when N is large. In terms of efficiency and accuracy, Fourier methods when compared to low-order methods have proven to be advantageous relative to memory requirement, computationally efficient and fast in execution times. The result obtained in Table 1 justifies the reason for abandoning the RK4 scheme in the main computations.

Table 1

The relative error for fractional diffusion equation for various values of discretization, at D = 0.5 and t = 1

N	N = 64	N = 128	N = 256	N = 512
RK4ETDRK4	1.960e−061.6391e−08	4.9539e−062.0899e−08	7.5637e−062.9017e−08	1.5528e−054.7984e−08
Ratio	72.9669	237.0400	260.6644	323.6079

Main model

Over the years, mathematical modelling has become an important tool in the application areas medical and life sciences to address some of the health challenging problems that are not approachable experimentally. Hepatitis B is commonly referred to as a life threatening infectious disease caused by the hepatitis B virus (HBV) which leads to inflammation or causing serious damage to the liver. According to World Health Organization’s (WHO) 2002 data report, over 2000 million people have been effected, more than 350 million individuals remain chronically infected and carriers of the virus. An estimated population of 4 million people are considered as acute clinical cases of the virus.

Research have shown that children and adolescents are most vulnerable to the disease than adults due to exposure which may show some clinical symptom and have a higher percentage chance of being acutely infected. Based on report, about a quarter of chronic infected individuals die of liver cancer annually. As a result, hepatitis B is known to be one of the most common viral source of cancer in the world nowadays. Hence, HBV infection is a disease of global health and its prevalence varies from one region to the other.

A lot of researchers have worked on HBV in the past, among which are the notable papers of Medley et al. (2001) where compartmentalized model was used to describe the spread of the disease. Almost a decade later, Zou et al. (2010) worked on the modified version of the model (Medley et al. 2001). They develop a model to explore the impact of vaccination and other controlling measures of HBV infection. Their model has simple dynamical behavior which has a globally asymptotically stable disease-free equilibrium when the basic reproduction number R0 < 1, and a globally asymptotically stable endemic equilibrium when R0 > 1. In the year 2014, Kimbir et al. (2014) give an extension to the earlier report in Zou et al. (2010) by including the treatment of chronically infected HBV carriers, it was also suggested in their report that the acute infected individuals are not subjected to antiviral treatment due to natural recovery. Wiah et al. (2015) employed a nonlinear extended deterministic model to address the impact of immigration on the population spread of HBV infection with acute and chronic infected carriers.

In this work, we present a deterministic model consisting of coupled nonlinear fractional partial differential equations of order α. This new model provides an extension of the models discussed earlier by Zou et al. (2010, 2015). The model population consists of five local kinetics broken into the susceptible individuals (U1), exposed class (U2) which are infected but yet to be infectious, acute infection individuals (U3), chronic HBV class (U4) and temporary protective immunity referred to as the recovered individuals (U5). The fractional reaction–diffusion system is given aswhere u = (u1, u2, … , u) is a vector of concentration or density for interacting species at position x and time t, and F, i = 1, … , 5, are the local reaction terms. The terms D > 0, i = 1, … , 5 are the diffusion coefficients. The initial densities are expected to be non-negative and the problem (17) is confined by imposing the appropriate choice of boundary conditions. The boundary conditions are taken as zero fluxwith similar relations for the U, i = 2, … , 5. The initial data are taken as in Eq. (2) in the form of some small perturbations from the uniform solution *:with analogous expressions for the remaining species. Subject to the following coupled reaction kineticswhere γ represents the rate of recruitment into a susceptible individuals, ω stands for the transmission rate of disease from one infection class to another. The HBV induced rate is given by β, while τ is the natural death rate. Parameters δ, σ are the progressive rate from the exposed (U2) to acute (U3) infection individuals and acute to chronic infection class (U4) respectively. The natural recovery rate for the exposed individuals for the latent HBV is denoted by ϵ, while ψ ≫ 1 is the transmission rate multiplier. Finally, ϕ and φ are the respective vaccination success rate for the class U1 and treatment success rate for (U3) infected class.

Stability analysis of the disease free equilibrium (DFE) point

In this section, we analyze the local stability of the disease-free equilibrium (DFE). It is the stability of at DFE that can guarantee a biologically meaningful results. Here we assume that the disease variables U2 = U3 = U4 = 0. If otherwise, the disease will persist and put the whole population of susceptible individuals into serious danger.

The basic reproduction number, commonly denoted as R0, gives the total number of secondary infections that an average infectious class will induce given that the rest of the population is susceptible. By using the notation in Pang et al. (2010), we denote the emergence of new infection by F, and the transfer of individual from one class to another by V. The endemic equilibrium dynamics in the region , do not correspond to biologically meaningful results since it encourages the spread of HBV disease. Hence, we do not capture the endemic equilibrium results in the analysis. For all possible parameter values, the spatially homogeneous stationary solution of model (17) with kinetics (20) has a disease free equilibrium point , we define the reproduction number as R0 = (FV−1), where

After some algebraic manipulations, we obtain the average value of the expected number of secondary cases produced by a single infected individuals

Theorem 2

The disease free equilibrium pointis locally asymptotically stable for the spatially homogeneous stationary solution of model (17) with kinetics (20) if R0 < 1, and unstable if otherwise.

The Jacobian matrix of the spatially homogeneous stationary solution of model (17) with kinetics (20) at point is given byFrom the characteristics equation of (22), we have two negative eigen values . After substitution, we have the rest of the characteristic polynomial given aswhere A, i, j = 1, 2, 3 is a 3 × 3 matrix, and I identity matrix. From (23), we obtain the characteristic equation .By adopting the Routh–Hurwitz stability conditions of the linear differential equations, system (17) with kinetics (20) is stable for the disease free equilibrium point if: (1) the roots r0, r1, r2, r3 are positive with negative real parts, where r3 = 1, r2 = a11a22 + a11a33, r1 = (a11a22 + a11a33 + a22a33 + a12a21), r0 = a12a21a33 − a12a22a33 − a13a21a32 (2) r1r2 − r0a3 > 0. We verified from the coefficients that r0 > 0, r1 > 0, r2 > 0, r3 > 0 and r1r2 − r0a3 = − R0 > 1 which implies that R0 < 1. Hence we complete the proof since R0 < 1 which shows that the disease free equilibrium point is locally asymptotically stable. □

Global stability of the diseases free equilibrium point

Here, we are concerned with the global stability of DFE point. We adopt a similar technique to the proof of Theorem 2.

Theorem 3

The disease free equilibrium pointis globally asymptotically stable for the spatially homogeneous stationary solution of model (17) with kinetics (20), if R0 < 1.

Let G be a Lyapunov function given in the formOn differentiating (24) and substitute the spatially homogeneous version of model (17) with kinetics (20) into it at , we haveon substituting for and in (15) we getIt follows from the point that G = 0 since the derivatives of Ui = 0, i = 2, 3, 4. This implies that as U → 0, i = 2, 3, 4, also U1 → γ/(τ + ϕ) and U5 → τϕ/(τ + ϕ)τ as t → ∞. So the DFE point is globally asymptotically stable if the inequality is satisfies the condition R0 < 1, but unstable if R0 > 1. The proof is completed. □

Numerical simulations

In this section, we start to simulate numerically the solution of the spatially homogeneous system (17) with kinetics (20) using the numerical techniques formulated in “Basic definitions and numerical techniques for fractional diffusion problems” section above to substantiate our analytical findings.

Non-diffusive example

We consider the set of parameterswith reasonable initial populations (millions)

For the set of ecological parameters, we realized that the conditions given in Theorems 2 and 3 for non-diffusive system are satisfied for the disease free-equilibrium state.

In Fig. 2, we simulate the non-diffusive system numerically at different instances of tau to study the behaviour of the species with time (year) t = 1. It was observed only the class of individuals U1, U2, U5 could actually be free of HBV as time is progressed. Clearly, those in the classes U3 and U4 which corresponds to the number of acute and chronic individuals will continue to live and spread the virus as years roll-by until the entire population is endermic. Similarly in Fig. 3, we fixed all parameters and initial conditions as in (26) and (27), and simulate with different instances of ϵ. It is obvious, despite the necessary measures such as treatments and lots of sensitization programs put in place for those in the acute and chronic class to get recovered with increasing time t = 3, the number of casualties is increasing. This is evident on the axis containing the population profiles of U3 and U4. Species attractor at t = 3 and instances of ϵ is given in panel (f).

Fig. 2

Time series results for the spatially homogeneous version of system (17) with kinetics (20), obtained at different instants of τ and final time t = 1. Other parameters are fixed in (26)

Fig. 3

Time series results for the spatially homogeneous version of system (17) with kinetics (20), obtained at different instants of ϵ and final time t = 3. Other parameters are fixed in (26)

Since it is not possible to completely wipe out both U3 and U4 classes in model (or in the population). Then, the question of ‘what can be done?’ sets in. In the context of this paper, we came out with the opinion of varying some of the parameters, the correct choice of parameters that will put control to the spread of HBV is attained by reducing the values of τ from 1.5 to 0.5 and that of β from 2.5 → 0.05 to enable the class of individuals U4 respond to treatment over time. To checkmate the spread of HBV for the group U3, we need to increase parameter value δ from 1.5 to 10.5 and above. These assertion is evident in Fig. 4.

Fig. 4

Behaviour of the acute and chronic individuals to treatment at instances δ and time t. Parameters are: τ = 0.5, β = 0.05 at t = 5 and t = 10 for species U 3 and U 4. Both species respond to treatment at the control parameters level with respect to time. Other parameters are fixed in (26)

Fractional reaction–diffusion example

Here, we simulate the whole system in the presence of diffusion with fractional power index α. We now illustrate through numerical experiments with Neumann boundary conditions, a multiple coexistence steady states in the fractional-in-space reaction–diffusion model (17) with kinetics (20). The boundary conditions are clamped the the extreme points of spatial domain of size [0, L], L = 20. We compute the initial conditions as:where N is the number of discretization.

We fixed parameters as:to obtain the surface plots displayed in Fig. 5. Obviously, the behaviour of the species correspond to those in Figs. 2 and 3, though the amplitudes differ.

Fig. 5

Surface plots for fractional reaction–diffusion system (17) with kinetics (20) at sub-diffusive scenario when α = 0.5 and t = 3. The dynamics obtained are similar to those presented in Figs. 2 and 3. Parameters are given in (26)

We repeat the parameter values used in Fig. 4 to obtain the surface plots displayed in Fig. 6 for the super-diffusive case, when α = 1.5. Clearly, all the species are disease free when necessary vaccinations and treatments (drugs) are administered over time. The risen amplitudes indicating high cases of HBV drop-down or tend to zero to show that the spread of HBV has significantly reduced, as shown in Fig. 6. By the computer simulation of the fractional reaction–diffusion system, we have given evidence that hepatitis B virus model and control can be studied in fractional scenarios. Our findings in this paper strongly recommend a combination of effective treatment and vaccination as a good control measure is important to record the success of HBV disease control as done via parameters τ, β and δ.

Fig. 6

Surface plots showing super-diffusive results for fractional reaction–diffusion system (17) with kinetics (20) obtained at fractional derivative index α = 1.5 and t = 5 to mimic results obtained in Fig. 4

Conclusion

In this paper, a mathematical model for investigating the hepatitis B virus disease in fractional medium is derived. The model disease free equilibrium state is analyzed. We established via theorems that the model disease-free equilibrium is both locally and globally asymptotically stable, if the basic reproduction number is less than unity. Our aim is to examine the behaviour of diffusive fractional reaction–diffusion model in sub-diffusive and super-diffusive scenarios, derive efficient and reliable numerical techniques. By the computer experiment of the fractional reaction–diffusion system we have given enough evidence that numerical solution in the diffusive (fractional) scenario, at 0 < α < 2 is practicably the same as in the case of non-diffusive case when applied to model Hepatitis B virus system. Our findings in this work strongly recommend a combination of effective treatment and vaccination as a good control measure is important to record the success of HBV disease control. It should be noted that the methodology presented in this paper can be applied to model other physical phenomena in higher dimensions.