Numerical Study for Time Delay Multistrain Tuberculosis Model of Fractional Order

A novel mathematical fractional model of multistrain tuberculosis with time delay memory is presented. The proposed model is governed by a system of fractional delay differential equations, where the fractional derivative is defined in the sense of the Grünwald–Letinkov definition. Modified parameters are introduced to account for the fractional order. The stability of the equilibrium points is investigated for any time delay. Nonstandard finite deferencemethod is proposed to solve the resulting system of fractional-order delay differential equations. Numerical simulations show that nonstandard finite difference method can be applied to solve such fractional delay differential equations simply and effectively.


Introduction
It is known that tuberculosis (TB) is one of the most important infectious diseases and is considered as the second largest cause of mortality by infectious diseases and a challenging disease to control [1]. Time delays required to treatment of active TB present a major obstacle to the control of a TB epidemic [2]; it worsens the disease, increases the risk of death, and enhances tuberculosis transmission to the community [3,4]. Both patient and the health system may be responsible for the treatment delay [3]. On the other hand, mathematical models are quite important and efficient tool to describe and investigate TB diseases; see [5][6][7][8][9]. In [10], Silva et al. presented TB model with time-delay memory. Herein, we consider a general model of multistrain TB diseases with time-delay memory. A discrete time delay is incorporated, in the variables of active TB infection of two and three strains, to represent the required time to commencement of treatment and diagnosis [11]. The multistrain TB model incorporates three strains: extensively drug-resistant (XDR), emerging multidrug resistant (MDR), and drug sensitive, and has been developed by Arino and Soliman [12] in 2015. Several factors of spreading TB such as the exogenous reinfection, the fast infection, and secondary infection are included in this model. Sweilam et al. introduced some numerical studies for this model in [13][14][15][16].
Fractional differential equations have been the focus of many studies due to their frequent appearance in various sciences [13][14][15][16][17][18][19][20]. The general theory of differential equations with delays (DDEs) is widely developed and discussed in the literature [21][22][23][24][25]. Delayed fractional differential equations (DFDEs) are also used to describe dynamical systems [26][27][28]. Recently, DFDEs begin to raise the attention of many researchers [29][30][31][32][33]. Relatable and efficient numerical techniques for DFDEs are very necessary and important [34]. Nonstandard finite difference method (NSFDM) was firstly proposed by Mickens [35] in 1980s to solve numerically the ordinary differential equations (ODEs) and partial differential 2 Complexity   equations (PDEs) with more accuracy than standard finite difference method (SFDM). It is considered as a powerful numerical scheme that preserves properties of exact solutions of the differential equation [36]. The main aim of work is to study numerically the solutions of fractional-order model of multistrain TB with time delay memory. The presence of fractional-order and time delays in the model can lead to a notable increase in the complexity of the observed behavior, and the solution continuously depends on all the previous states. An efficient numerical method, NSFDM, is used to numerically solve the fractional-order delay model. The rest of the paper is organized as follows: In Section 2, we present a fractional order model with time delay for multistrain TB. Stability of equilibrium points is presented in Section 3. NSFD for fractional-order delay differential equations is introduced in Section 4. Some numerical simulations are given in Section 5, and conclusion in Section 6. Some definitions on fractional calculus and some properties of nonstandard discretization are given in Appendix.

Fractional Multistrain TB Model with Time Delay
In this section, a multistrain TB model of fractional-order and time delay memory is presented. The population of interest is divided into eight compartments depending on their epidemiological stages as follows: susceptible ( ); latently infected with drug sensitive TB ( ); latently infected with MDR TB ( ); latently infected with XDR TB ( ); sensitive drug TB infectious ( ); MDR TB infectious ( ); XDR TB infectious ( ); recovered . One biological meaning of the given parameters is given in Table 1. One of the main assumptions of this model is that the total population ( ), with ( ) = ( ) + ( ) + ( ) + ( ) + ( ) + ( ) + ( ) + ( ), is variable of the time. We introduce a discrete time delay in the state variables and , denoted by , that represents the time required for diagnosis and commencement of treatment of active TB infection of two and three strains. The parameters in the modified the model are described in Table 2; see [37]. The modified system of multistrain TB model of fractional-order and time delay is The initial conditions for system (1) From biological meaning, we further assume that > 0 for = 1, . . . , 8. Throughout this paper, we focus on the dynamics of the solutions of (1) in the restricted region, Ω = {( , , , , , , , ) ∈ R 8 | + + + + + + + ≤ / }. We refer here to [24,31], for local existence, uniqueness, and continuation results.

Basic Reproduction Number.
The basic reproduction number, 0 , is defined as the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual [32]. Herein, we apply the method in [32] to drive 0 . The order of the infected variables is I fl ( , , , , , ) . ( The vector representing new infections into the infected classes is given by ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) . ( The vector representing other flows within and out of the infected classes I is given by ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) .
The matrix of new infections and the matrix of transfers between compartments are the Jacobian matrices obtained by differentiating and with respect to the infected variables I and evaluating at the disease-free equilibrium. They take the form Then the basic reproduction number 0 for system (1) is the spectral radius of the next generation matrix and is given by where

Equilibrium Points and Their Asymptotic Stability
To discuss the local asymptotic stability for evaluating the equilibrium points, let us consider the following [38]: Let us consider the coordinate transformation: The characteristic equation associated with above matrix is [38] Δ ( ) = ( 0 ) − = 0, Lemma 1. If 0 < 1, then the disease-free equilibrium 0 is locally asymptotically stable for = 0.
Consider now the case > 0; we noted that the roots of the characteristic equation (15) have no pure imaginary roots for any value of the delay , if 0 > 1. Hence all the roots of the characteristic equation have negative real parts and all eigenvalues satisfy Matignon's conditions [39]. Therefore, the endemic equilibrium is locally asymptotically.

Numerical Results and Simulations
In this section, we show the effectiveness of the numerical technique for delay fractional differential equations. Throughout this section, all simulations are performed with initial conditions ( (0), (0), (0), (0), (0), (0), (0), (0)) = (5000, 50, 50, 50, 30, 30, 30, 60), with the parameters in Table 3. The approximate solutions of the proposed system are given in Figures 1-12 at different values of and . Figure 1 shows the behavior of the approximate solutions of ( ) in two cases with and without delay using dde23 at = 1 and = 0.3. In Figure 2, we use the same data in Figure 1 and use NSFDM; we noted that the number of individuals ( ) increases in the case of nondelay, that is, the delays in diagnosis and commencement of treatment to the individuals of and causing a shortage in the number of individuals of ( ). Figure 3 shows the relationship between ( ) and ( − ) and chaotic attractors at = 0.1 and = 1.   Figures 11 and 12 show the behavior of the approximate solutions with different value of , which are given to demonstrate how the fractional model is a generalization of the integer order model.

Conclusion
Fractional models have the potential to describe more complex dynamics than the integer models and can include easily the memory effect present in many real world phenomena. In this paper, multistrain TB model of fractional order derivatives with time delay memory is presented. A nonstandard numerical scheme is introduced to numerically study the approximate solution of proposed model problem.
The obtained results show that the delays in diagnosis and commencement of treatment to the individuals of and cause a shortage in the number of individuals of ( ). The approximate solution of proposed model changes when and   take different values. Some figures are given to demonstrate how the fractional delay model is a generalization of the integer order model. It is concluded that NSFDM can be applied to solve such fractional delay differential equations simply and effectively. All results are obtained by using MATLAB programming.

A. Preliminaries and Notations
In this section, some basic definitions and properties in the theory of the fractional calculus are presented. Moreover, we introduce the main aspects concerning nonstandard discretization methods.
A.1. Grünwald-Letinkov Fractional Derivatives (GLFDs). We will begin with the signal fractional differential equation (see [17,40,41]) A.2. NSFD Discretization. It is known that the numerical scheme is called nonstandard method if at least one of the following conditions is satisfied [36]: (1) the discretization of derivatives is not traditional and uses a nonnegative function [35,43], (2) nonlocal approximations are used.

Complexity 13
To construct the numerical scheme for system (1) using NSFDM, the approximations of temporal derivatives are made based on generalized forward scheme of first order. Hence, if ( ) ∈ 1 (R), we define its derivative as follows: ( ) = ( + △ ) − ( ) (△ ) + ( (△ )) , as △ → 0, where (△ ) is a real-valued function on R and △ = ℎ. In the following, the denominator functions are little complex functions of the step-size of time than the classical one [44].