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.

The fractional derivative operator (with similar expressions for u i , i = 1, 2, … , n) is the Riemann-Liouville fractional gradient of order α. The diffusion matrix is defined as a diagonal matrix D = diag(D 1 , D 2 , … , D n ), and the diffusion coefficients D i 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 γ im > 0 for i � = m. The reaction kinetics in (1) is expected to satisfy the following conditions: 1. In R n + = {(u 1 , u 2 , … , u n )|u i > 0|}, the vector field (F 1 (u), F 2 (u), … , F n (u)) has unstable equilibrium state at 0 = (0, 0, … , 0) and asymptotically stable at A = (A 1 , A 2 , … , A n ) with A i > 0 for i = 1, 2, … , n. 2. The coefficients must 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. 2014Zeng et al. , 2015Zheng 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 subdiffusive (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 f ∈ C n −1 , 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 α � = 1, C α = 1 2 cos ( πα 2 ) 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 model
where μ, t > 0, x ∈ R, α, η, σ, ν are real and positive parameters subject to the constraints 0 < α ≤ 2, |ϑ| < min(α, 2 − α), and D η,v t is the generalized Riemann-Liouville fractional derivative operator given by with the conditions where 1 < η ≤ 2, 0 ≤ ν ≤ 1, ρ > 0 is a constant with the kinetic (or reaction) term. x D α ϑ is called the Riesz-Feller space fractional derivative of order α, skewness ϑ is symmetry given in terms of the Fourier transform where ̟ α ϑ (k) = k α e i(signk) ϑπ 2 , 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 as for n − 1 < η ≤ n, n ∈ u, 0 ≤ ν ≤ 1, we then transform the given equation into the form where * 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 ū * (k, s) we obtain We take the inverse Laplace transform of Eq. (8) and following the result of Saxena et al. (2006), we have where the real part of the parameters are positive, bear in mind that which by inverse Fourier transform, Eq. (9) becomes as 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 get where U is the Fourier transform of species u(x, t). This implies that, We let Ω α = (X α x ) to explicitly remove the issue of stiffness in the fractional derivative term, and set U = e −Ω α tŪ , so that Now, in practice, we discretize spatial domain by considering N x 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 ODEs 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 = nΔt to t n+1 = (n + 1)Δt for the ODE where Here, we apply the general explicit Runge-Kutta scheme to (13) for Ū i . For brevity, we denote μ i as k′ in U i and set the replacement variable as Finally, we write the s-stage RK scheme as with modified term and 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 μ i given as The functions φ i , are given as where L = ∂ α u ∂x α 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 as where ū 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 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.
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.

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  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 R 0 < 1, and a globally asymptotically stable endemic equilibrium when R 0 > 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. (2010Zou et al. ( , 2015. The model population consists of five local kinetics broken into the susceptible individuals (U 1 ), exposed class (U 2 ) which are infected but yet to be infectious, acute infection individuals (U 3 ), chronic HBV class (U 4 ) and temporary protective immunity referred to as the recovered individuals (U 5 ). The fractional reaction-diffusion system is given as where u = (u 1 , u 2 , … , u n ) is a vector of concentration or density for interacting species at position x and time t, and F i , i = 1, … , 5, are the local reaction terms. The terms D i > 0, i = 1, … , 5 are the diffusion coefficients. The initial densities are expected to be nonnegative and the problem (17) is confined by imposing the appropriate choice of boundary conditions. The boundary conditions are taken as zero flux with similar relations for the U i , 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 kinetics where γ 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 (U 2 ) to acute (U 3 ) infection individuals and acute to chronic infection class (U 4 ) 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 U 1 and treatment success rate for (U 3 ) 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 U 2 = U 3 = U 4 = 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 R 0 , 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 (U 1 > 0, U 2 > 0, U 3 > 0, U 4 > 0, U 5 > 0) = (Û 1 ,Û 2 ,Û 3 ,Û 4 ,Û 5 ), 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 Ê = (γ /(τ + φ), 0, 0, 0, τ φ/(τ + φ)τ ), we define the reproduction number as R 0 = (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 point Ê = (γ /(τ + φ), 0, 0, 0, τ φ/(τ + φ)τ ) is locally asymptotically stable for the spatially homogeneous stationary solution of model (17) with kinetics (20) if R 0 < 1, and unstable if otherwise. (20) Proof The Jacobian matrix of the spatially homogeneous stationary solution of model (17) with kinetics (20) at point Ê is given by From the characteristics equation of (22), we have two negative eigen values 1 = −τ , 2 = −(τ + φ). After substitution, we have the rest of the characteristic polynomial given as where A ij , i, j = 1, 2, 3 is a 3 × 3 matrix, and I identity matrix. From (23), we obtain the characteristic equation X( ) = 3 − (a 11 a 22 a 33 ) 2 + (a 11 a 22 + a 11 a 33 + a 22 a 33 +a 12 a 21 ) − a 12 a 21 a 33 − a 12 a 22 a 33 − a 13 a 21 a 32 = 0.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 r 0 , r 1 , r 2 , r 3 are positive with negative real parts, where r 3 = 1, r 2 = a 11 a 22 + a 11 a 33 , r 1 = (a 11 a 22 + a 11 a 33 + a 22 a 33 + a 12 a 21 ), r 0 = a 12 a 21 a 33 − a 12 a 22 a 33 − a 13 a 21 a 32 (2) r 1 r 2 − r 0 a 3 > 0. We verified from the coefficients that r 0 > 0, r 1 > 0, r 2 > 0, r 3 > 0 and r 1 r 2 − r 0 a 3 = − R 0 > 1 which implies that R 0 < 1. Hence we complete the proof since R 0 < 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.

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 parameters with 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 U 1 , U 2 , U 5 could actually be free of HBV as time is progressed. Clearly, those in the classes U 3 and U 4 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 U 3 and U 4 . Species attractor at t = 3 and instances of ϵ is given in panel (f ).
Since it is not possible to completely wipe out both U 3 and U 4 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 (26) φ = 2.5, γ = 2.0, ω = 0.5, ψ = 2.5, δ = 0.5, ε = 1.5, ϕ = 0.75, β = 2.5  (26) of individuals U 4 respond to treatment over time. To checkmate the spread of HBV for the group U 3 , we need to increase parameter value δ from 1.5 to 10.5 and above. These assertion is evident in Fig. 4.

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 reactiondiffusion model (17)  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. 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 U 1 = 95.5 * ones(N, 1); U 2 = 2 * ones(N, 1); U 3 = 3 * ones(N, 1); U 4 = 1.5 * ones(N, 1); U 5 = 1 * ones(N, 1); (28) {φ = 2.5, τ = 1.5, γ = 2.0, ω = 0.5, ψ = 2.5, δ = 1.5, ε = 0.5, ϕ = 0.75, β = 2.5, D1 = 0.2, D2 = 0.5, D3 = 0.25, D4 = 0.1, D5 = 0.5} 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 δ.

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.