Nonstandard theta Milstein method for solving stochastic multi-strain tuberculosis model

In this article, a novel stochastic multi-strain tuberculosis model is presented. Numerical simulations for this model are the main aim of this work. A non-standard theta Milstein method is constructed to study the proposed model, where the proposed method is based on choosing the weight factor theta. The main advantage of this method is it can be explicit or implicit with large stability regions using the idea of the weighed step introduced by R.E. Mickens. Mean-square stability of nonstandard theta Milstein method is studied. The new scheme shows a greater behavior compared to the theta Milstein method. It is concluded that the proposed scheme preserves the positivity of the solution and numerical stability in larger region than the standard method.


Introduction
It is well known that the solution of stochastic differential equations is either difficult in general or we do not have explicit solutions. Numerical schemes provide an easy way to integrate these equations, see [1][2][3][4][5][6][7][8]. Moreover, numerical simulations are considered the only way to solve these mathematical models in general or to derive the desired information. Therefore, the accuracy of these numerical solutions could be a major factor in choosing the appropriate numerical method and solving mathematical models. The usual numerical schemes even in the deterministic case such as Euler, Runge-Kutta, and Euler-Maruyama in the stochastic case do not preserve dynamical properties without conditions on the time step of the numerical integration; see [9] and references therein. The non-standard finite difference (NSFD) schemes were firstly proposed by Mickens

A stochastic multi-strain TB model
In this section, we introduce the multi-strain TB model which is given in [33]; this model incorporates three strains: drug-sensitive, emerging multi-drug-resistant (MDR), and extensively drug-resistant (XDR). The population of interest is divided into eight compartments, see Table 1. Let us assume that ( , , { t } t≥t 0 , P) is a complete probability space with a filtration { t } t≥t 0 . Let W (t) be an eight-dimensional Wiener process defined on this probability space. We assume that the eight coordinates W 1 (t), W 2 (t), ..., W 8 (t) are mutually independent and ξ 1 , ξ 2 , , ..., ξ 8 represent the intensities of the white noises. The stochastic perturbation in our model is a white noise type that is directly proportional to the all variables of the model. Therefore, the stochastic system can be described by the Itô system as follows: dR =P 1 t 1s L s + P 2 t 2s I s + To investigate the dynamical behavior of the proposed model, the first thing is whether we have a unique global solution (i.e., no explosion in a finite time) for any given initial value. Consider the following d-dimensional stochastic system: where F(t, X(t)) is a function in R d defined in [ t 0 , ∞[ ×R and G(t, X(t)) is a d × m matrix, F and G are locally Lipschitz functions in X, and W = W (t) t≥0 is a d-dimensional Wiener process. We assume that X = 0 is a solution of system (1)- (8). The coefficients of Eq. (9) are generally required to satisfy linear growth conditions and local Lipschitz conditions [34][35][36]. If the coefficients F and G of (9) do not satisfy linear growth conditions, the solution of system (9) may explode at a finite time. From [36], the coefficients are locally Lipschitz. Consequently, the system has a unique local solution for any feasible initial state.
In order to study stability of disease-free equilibrium point E * , we assume that there are global solutions which are almost surely non-negative.

Stability of the SDE model for multi-strain tuberculosis model
In this section, the theorem below can be interpreted as follows; at least, the stochastic perturbations do not destabilize the system. Let (R 0s , R 0m , R 0x ) the basic reproduction number of drug-sensitive strain, MDR strain, and XDR strain. Let us define R * as follows: Theorem 1 If R * < 1, then disease-free equilibrium is a.s exponentially stable.
Proof Assume that [33]: Using the fact that S(t)+σ s R(t) N < 1 and that I s (t) ≤ I ∞ s at any t, it follows Using assumption (11), and for simplicity, let us define a 1 : Similarly, using assumptions (12) and (13), we can prove the following inequalities involving I m and I x . Define a C 8 - Thus, we can define Z = ln V by the Itô formula; we compute, Then, Therefore, where each G i (t) is a martingale defined as: Regarding the quadratic variations of the stochastic integral G i (t), we have t 0 By the strong law of large numbers for martingales [36], we therefore have It finally follows from (16) by dividing t on both sides and then letting t → ∞ that LVdt (a.s).
We note that This finally proves the (a.s.) exponential stability.

Basic properties of the solution
In this section, we remind classical results about continuous and discrete stochastic differential equation systems.

Stochastic boundedness
Stochastic boundedness is one of the most important topics because boundedness of a system guarantees its validity in a population system. We first present the definition of stochastically ultimate boundedness. ∈ (0, 1), there is a positive constant δ = δ( ) such that for any initial value X 0 ∈ R 8 + , the solution X(t) to (1)- (8) has the property that:

A nonstandard theta Milstein method
One of the merits of NFDM which is given by Mickens in 1980 [10, 22-24, 28, 32] is used to study numerically the behavior of the ordinary differential equations (ODE) and PDE efficiently. The NFDM is able to maintain the properties of the exact solution of the original ODE or PDE. The numerical scheme is called NFDM if at least one of the following conditions are satisfied [10][11][12][13][14][15][22][23][24][25][26][27][28]: 1. The nonlocal approximation is used. 2. The discretization of the derivative is not traditional and uses a nonnegative function.
In this paper, we consider numerical methods for the strong solution of Itô SDEs: where z(t) is a random variable with value in R m , f : R m −→ R m is called the drift function, g : R m −→ R m is called the diffusion function, and dW (t) is a one-dimensional Wiener process, whose increment W (t) = W (t + t) − W (t) is a Gaussian random variable N(0, t).

Definition 2 A general one-step theta Milstein method (TMM) is given as follows:
Definition 3 A general one-step nonstandard theta Milstein method (NTMM) is given as follows: where

Remark 1 In (19) if ϕ(h)
= h, then the scheme will be called theta Milstein method and its two special cases for θ are as follows: • If θ = 1, we obtain explicit Milstein method (EMM).

Mean-square stability of NSTMM
The mean-square stability is a stochastic version of absolute stability, and it is a very important concept in the numerical simulation of SDEs. A suitable way to find numerical schemes for SDEs is analysis of MS stability.
Definition 4 [38] The equilibrium position, z(t) ≡ 0, is said to be mean-square stable if for every ε > 0 there exists δ 1 > 0 such that where If, in addition to (20), there exists a δ 2 > 0 such that then the equilibrium position is said to be asymptotically mean-square stable.
Definition 5 [38] Suppose that the equilibrium position of Itô's SDE as (17). The Gaussian random variable N(0, t) is asymptotically mean-square stable. Then, a numerical scheme that produces the iterations z n to approximate the solution z(t) of (21) is said to be asymptotically mean-square stable if lim n−→∞ z n = 0.

years Assumed
with known solution z(t) = z 0 e ( a−b 2 2 ) , which is represented by where J is the standard Gaussian random variable J = W n √ h  ∼ N(0, 1). Saito and Mitsui [39] introduced the following definition of mean-square (MS) stability.

Definition 6 [39] The numerical method is said to be MS stable for a, b, h if
whereR(a, b, h) is called MS stability function of the numerical method. For θ = 1 of (18), EMM is given as follows, for θ = 0 of (18) by introducing implicitness in f (t n , z n )h, we arrive at IDMM. If in (23) we replace h by ψ(h), we have a new method called NSIDMM; this special case from (19) when theta = 0 and the function ψ(h) satisfies the following conditions: For more details about nonstandard method, see [10, 22-24, 28, 32] and the references cited therein. Applying the NSIDMM which given as: to the linear test Eq. (21), we obtain The MS stability function of the NSDIM method is given by The NSIDM will be MS stable if R 1 (p, q) < 1. The MS stability property of the NSIDMM is better than that of IDMM.

Conclusions
In this paper, we constructed NTMM to introduce numerically the approximate solution of a stochastic multi-strain TB model. The proposed method is based on choosing the weight factor θ. The main advantage of this method is that it can be explicit or implicit with large stability regions as we see in our results. Special attention is given to study mean-square stability. Some numerical results are used to show the accuracy of the NTMM, and some figures are used to demonstrate how the solutions change when θ take different values. All computations in this paper are performed using MATLAB programming.