On the Distributed Order Fractional Multi-Strain Tuberculosis Model: A Numerical Study

In this paper, a novel mathematical distributed order fractional model of multi-strain Tuberculosis is presented. The proposed model is governed by a system of distributed order fractional differential equations, where the distributed order fractional derivative is defined in the sense of the Grünwald−Letinkov definition. A nonstandard finite difference method is proposed to study the resulting system . The stability analysis of the proposed model is discussed. Numerical simulations show that, the nonstandard finite difference method can be applied to solve such distributed order fractional differential equations simply and effectively.


Introduction
In fact, the fractional calculus has been acknowledged as a promising mathematical tool to efficiently characterize the historical memory and global correlation of complex dynamic systems, phenomena or structures. However, various literature indicated that the memory and/or nonlocality of the system may change with time, space or other conditions [1], [2]. The variable-order (VO) fractional operators depending on their non-stationary powerlaw kernel can describe the memory and hereditary properties of many physical phenomena and processes [3], [4].
It is known that the distributed order derivatives are fractional order derivatives that have been integrated over the order of the derivative within a given range. The idea of fractional derivative of distributed order is stated by [18] and later developed by Caputo himself in [19,20] and Bagley and Torvik [21]. Many researchers used this idea, they applied it to some interesting mathematical models of partial fractional differential equation of distributed order. Diethelm and Ford [22] used a numerical technique along with its error analysis to solve the distributed order differential equation and analyze the physical phenomena and engineering problems; see [22] and the references cited therein. Recently Saberi Najafi et al. [23,24] studied the stability analysis of distributed order fractional differential equations with respect to the nonnegative density function.
Mathematical models can project how infectious diseases progress, to show the likely outcome of an epidemic, and help inform public health interventions. In epidemiology, compartmental models serve as the base mathematical framework for understanding the complex dynamics of these systems [5].
where n is the first integer which is not less than α, that is, n − 1 < α < n and Γ(·) is a Gamma function. The Caputo fractional derivative of z(t) is defined as It is known that Grünwald−Letinkov's definition and Riemann-Liouville definition are equivalence, but the Grünwald−Letinkov's definition is very easily utilized for the numerical evaluations see [8,11]. Hence, in order to apply nonstandard finite difference schemes laters, we have chosen Grünwald−Letinkov's approximation formula for the distributed order fractional derivative [11]: where, m = [ x h ] denotes the integer part of x h and h is the step-size. Equation 3 can be discretized as follows: Where g(t m , z(x m )) = C D α t z(t), x m = mh , and ω α r , are the Grünwald−Letinkov coefficients define as follows:

NSFD Discretization
It is known that, the numerical scheme is called nonstandard method if at least one of the following conditions are satisfied [30]: 1. The nonlocal approximation is used.
2. The discretization of the derivative is not traditional and uses a nonnegative function ( [12], [14]).
For more details on the nonstandard method see [16] and the references cited therein.

Remark 2
Almost all of standard procedures yield schemes which are convergent with restriction of step size while NSFDM are convergent for any step size. Also, in addition to the usual properties of consistency, stability and hence convergence for nonstandard finite difference schemes, they produce numerical solutions which also exhibit essential properties of the solution, for more details see ( [12]- [14], [16]).

Distributed Order Fractional Multi-Strain TB Model
In this section, the fractional multi-strain TB model is modified by integrating over all possible orders of the fractional time derivative. This is called a distributed order derivative. The population of interest is divided into eight compartments depending on their epidemiological stages as follows: susceptible (S); latently infected with drug sensitive TB (L s ); latently infected with MDR TB (L m ); latently infected with XDR TB (L x ); sensitive drug TB infectious (I s ); MDR TB infectious (I m ); XDR TB infectious (I x ); recovered R. All interpretation and meaning of parameters for this model see [28]. One of the main assumptions of this model is that, the total population ON THE DISTRIBUTED ORDER FRACTIONAL MULTI-STRAIN TUBERCULOSIS MODEL system is described by distributed order fractional derivatives as follows:

Stability analysis of distributed order fractional systems
The linear distributed order fractional systems can be expressed as: 8 , the matrix B ∈ R n×n and q(α), the density function, 0 < α ≤ 1. Saberi Najafi et al. [24] have obtained the general solution of the distributed order fractional systems (18), as follows: , r = e iπ and A(s) = ρ cos(πγ) + iρ sin(πγ). Now, we will recall some theorems and definitions on linear distributed order fractional equations and then we will show that theorem for nonlinear distributed order fractional equations as well.

Remark 4
The value of det(A(s)I − B) = 0, is the characteristic function of the matrix B with respect to the distributed The inertia of the system (18) is the triple: Theorem 3.2 [24] The linear distributed order fractional system (18) is asymptotically stable if and only if any of the following equivalent conditions holds: Next, we will mainly discuss the stability of a nonlinear autonomous distributed order fractional system, which can be described by .
System (21) can be written as with the initial value ϱ(0) = Y 0 − Y * . The analytical procedure of linearization is based on the fact that if the matrix J has no purely imaginary eigenvalues, then the trajectories of the nonlinear system in the neighborhood of the equilibrium point have the same form as the trajectories of the linear system (23). Hence, by applying Theorem 3.2, the linear system (23) is asymptotically stable if and only if all roots s of the characteristic function of J with respect to A(s) = (A 1 (s), A 2 (s), ..., A n (s)) T satisfy | arg(s)| > π 2 , which implies that the equilibrium Y * of the nonlinear distributed order fractional system (20) is as asymptotically stable.

Remark 5
The nonlinear distributed order fractional system (20) in the point Y * is asymptotically stable if and only if Θ nA(s) B = Φ nA(s) B = 0.

NSFD for Distributed Order Fractional Differential Equations
NSFD schemes were introduced by Mickens in the 1980s as a powerful numerical method that preserves significant properties of exact solutions of the involved differential equation [7]. We defined distributed order fractional derivative using Grünwald−Letinkov's sense and applying NSFD. Then by using the midpoint quadrature rule. The system 10 − 17 can be discretized as follows: The discretizations for N (t) is given as: N n = S n + L n s + L n m + L n x + I n s + I n m + I n x + R n . Where, ω αi 0 = (φ s (h)) −αi , s = 1, 2, · · ·, 8 and 0 < α i < 1 . The nonlocal approximations are used for the nonlinear terms and the following denominator functions are used: Then, we can obtain

Numerical Results and Discussions
The purpose of this section is to show that NSFD designed in this paper provides good approximations for distributed order fractional differential equations. In all simulations we use the initial conditions (S(0), L s (0), L m (0), L x (0), I s (0), I m (0), I x (0), R(0)) = (5000,50, 50,50,30,30,30,60), with the parameters in Table 1. The approximate solutions of the proposed system are given in figures 1-5 at different values of q(α). In Fig. 1, the effect of different q(α) on the behavior the approximate solutions for state variables I s , I m , I x and R with comparered the obtianed results with integer case when q(α) = 1. This figure is given to demonstrate how the distributed order fractional model is a generalization of the integer order model. Fig. 2, shows the behavior of the approximate solutions for state variables S(t) and I x (t) by using NSFDM and SFDM at h = 2, we can clearly see, NSFDM converge to correct endemic equilibrium but the obtain soluations of SFDM is unstable. Fig. 3, shows the behavior of the approximate solutions of the state variables I m , I x , R at different Dirac delta function by using NSFDM. Fig. 4 show that, the obtianed solutions when q(α) = δ(α − 0.75) are excellent agreement with q(α) = 0.75. Fig.  5, shows the relationship between S(t) with I x (t) and R(t) with I m (t) when q(α) = δ(α − 0.5) using NSFDM. For analyzing the stability of the proposed system, we compute Υ n A(s) (J) in the case that the density varies. The results are shown in Table 2. Numerical results are provided to show that NSFD is computationally efficient. From the numerical experiments, we observed that NSFD schemes allow to replicate quite well the behavior of the solution in the classical problems with integer-order. Hence, we can conclude that the developed NSFD schemes are useful for solving distributed order fractional equations.

Conclusions
In this paper, a distributed order fractional multi-strain TB model has been presented as a general model. A nonstandard numerical scheme has been introduced to numerically study the approximate solution of proposed model problem. Some figures are given to demonstrate how the distributed order fractional model is a generalization of the integer and fractional order models where this dynamical model is more suitable to describe the biological phenomena with memory than the integer and fractional order model. The stability of the equilibrium points is investigated. It is concluded that NSFDM can be applied to solve such distributed order differential equations simply and given the large stability regions.