Dynamical Analysis and Design of Computational Methods for Nonlinear Stochastic Leprosy Epidemic Model

In this article, we present the dynamical analysis of the stochastic leprosy epidemic model. Positivity and boundedness are the criteria used in the deterministic model. A primary technique known as the Euler Maruyama is employed in the solution of the said model. Standard and non-standard computational methods are applied in evaluating the design stability and ef-ﬁciency based on the chosen criteria. The standard computational methods like the Stochastic Euler and the Stochastic Runge Kutta fail to restore the essential features of biological problems. However, our proposed approach, the stochastic non-standard ﬁnite diﬀerence (NSFD), is used and found to be eﬃcient, cost-eﬀective, and accommodates all the desired feasible properties. Our method achieves all-time convergence against the backdrop of other classical techniques that perform conditionally or fail over a long period. In the end, a comparison between this scheme and the existing ones reviews the novelty of our approach.


Introduction
In 2018, Giraldo et al. studied a mathematical model based on nonlinear ordinary differential equations. The model explains the transmission dynamics of paucibacillary leprosy (PB) and multibacillary leprosy (MB) [1]. In 2018, Varella et al. presented the mathematical model related to leprosy type infection: a case study of the spread of leprosy in Brazil [2]. In 2021, Marathi et al. developed a multi-dimensional mathematical model for the transmission of leprosy disease [3]. In 2016, Smith captured a compartmental model representing the present knowledge of susceptibility and pathogenesis to leprosy, adapted to different populations [4]. In 2007, Abubakar et al. implemented the Markov decision model to decide the best cost of drugs for the treatment of leprosy disease [5]. In 2012, Mushayabase et al. investigated the dynamics of leprosy on behalf of treatment and control measures [6]. In 2013, Ghayas et al. studied leprosy disease transmission dynamics: a case of Gilgit-Baltistan, Pakistan [7]. In 2013, Deps et al. linked the co-infection among leprosy and HIV [8]. In 2013, Chiyaka et al. developed numerical and theoretical analysis for transmission dynamics of leprosy disease [9]. Similarly, Arco et al., in 2016, diagnosed the leprosy symptoms and transmission modes [10], while Le et al., in 2018, investigated the control plans for the disease-like infections: a case study of monitoring and detection in China [11]. In 2017, Nobre et al. presented a high-risk analysis of leprosy infection like the multibacillary [12]. Ezenduka et al., in 2012, analyzed the cost-effectiveness of the three methods of detecting leprosy in Northern Nigeria to determine the most viable plan for locating the disease [13]. In 2009, Koba et al. analyzed the trends and patterns of native leprosy in Japan [14]. However, in the countries of the Amazon region, not until 2020 when Schaub et al. in [15] demonstrated the prevalence of distinct features of leprosy among humans and animal reservoirs. Kumar et al., in 2007, came out with a survey of leprosy disease among non-family contacts (NFC) and the family contacts (FC) of patients in Agra, India [16], while Fischer et al., in 2008, assessed the local distribution of leprosy cases in Bangladesh during 15-years of disease control program in Bangladesh [17]. The rest of this paper is organized based on the following sections: formulation of the deterministic leprosy epidemic model with fundamental properties is introduced in section 2. In section 3, the transition probabilities, positivity, boundedness and implementation of the well-known methods, including the NSFD, and their comparative analysis are exposed. Finally, section 4 offers the concluding remarks.

Formulation of the model
In this part, we studied the dynamics of a leprosy disease in humans using a mathematical model in which the human population is divided into four subpopulations, respectively. The total number is represented by N : [0, ∞) → R with the function of time t ≥ 0. Furthermore, each of the subcomponents of  population is denoted by a non-negative differentiable function x, y, z, w : [0, ∞) → R. The representation of each of this subdivision is described as follows: x(t) depicts the susceptible individuals, y(t) denotes the asymptomatic infected individuals, z(t) represents the infected individuals by multibacillary leprosy, w(t) shows the infected individuals by paucibacillary leprosy, and N (t) represents the total population. The systematic flow of leprosy disease is presented in Fig. 1.
The physical interpretation of the model is presented as follows: f represents a fraction of humans who developed multibacillary leprosy, 1-f portrays a fraction of humans who developed paucibacillary leprosy, ρ describes the natural growth rate of humans, β p delineates the contact rate of paucibacillary leprosy individuals, β m shows the contact rate of multibacillary leprosy individuals, θ presents the rate at which asymptomatic individuals may be symptomatic individuals, µ m depicts the death rate of the infected individuals by multibacillary leprosy, while µ is the natural death rate of individuals. The model of the differential equations is as given below: with non-negative initial conditions x (0) ≥ 0, y (0) ≥ 0, ω (0) ≥ 0, z (0) ≥ 0.

Model analysis
In this section, we discuss feasible region of the model. The variables x (t), y (t), ω (t) and z (t) must be non-negative. Therefore, the obtainable results of the epidemic model remain positive and bounded in the region M = (x, y, z, ω) ∈ R 4 + : N (t) ≤ κ µ , x ≥ 0, y ≥ 0, z ≥ 0, ω ≥ 0 at any time t ≥ 0.

Positivity of the model
The solution (x, y, ω, z) ∈ R 4 + of the system (1)- (4) is positive at any time t ≥ 0 given any non-negative initial conditions. Proof: It is clear from equations (1)-(4) that: as desired.

Theorem 2 The solution of the system (1)-(4) is bounded and lie in the feasible region M
Proof: Let us consider the population density function as follows: From Gronwall's inequality, So, it is clear that the solution of (1)-(4) is bounded and lies in M.

Local stability
In this section, we present two well-known theorems for local stability. Again, consider the system (1)-(4) as a function of F, G, H and K as follows: The partial derivates of the function regarding state variables are as follows: Thus, the Jacobian matrix of the model takes the following form Proof: The Jacobian matrix evaluated at L 1 is as follows Letting which is transformed into its counterpart polynomial form as If we examine the coefficients of the characteristic equation above as The Routh-Hurwitz stability criterion for the 3 rd -degree polynomial is satisfied, and hence the disease-free equilibrium is stable.
Proof: The Jacobian matrix of the system at L 2 is as follows: and proceeding in the same manner as before, we have: Consequently, the characteristic equation becomes and We observe that the Routh-Hurwitz stability criterion for the 4 th -degree polynomial is also satisfied, since Hence, the endemic equilibrium is locally stable.

Reproduction number
In this section, we investigate the next-generation matrix method to find the transmission and transition matrices (F and V respectively) of the system (1)-(4) after substituting the disease-free equilibrium as follows: The dominant eigenvalue of the matrix F V −1 is called the reproduction number and is denoted as

Stochastic model
Let us consider the vector C = [x, y, z, w] T for the system (1)-(4). Transition probabilities are presented in Table 1 for the expectations E * [∆C] and variance E * ∆C∆C T . Table 1 Possible changes in the process of the model.

Transition
Probabilities We have the following quantities: where: The following additional quantities are defined: Finally, the stochastic differential equation based on (10) and (11) takes the following form

Euler Maruyama scheme
In this section, we utilize the Euler Maruyama scheme to determine the numerical solution of the stochastic differential equation (12). The employed model parameters are reported in Table 2. The following computational procedure holds: β m x n z n N + β p x n w n N + θy n + µy n 0 −θy n 0 0 f θy n + µz n −f θy n 0 −θy n −f θy n θy n + f θy n + µw n where ∆t stands for the discretization parameter.

Parametric perturbation of the model
In what follows, we introduce parametric perturbation in the system (1)-(4) as follows and we get where σ i , i = 1, 2 are the randomness of the model and B(t) is the Brownian motion.

Positivity and boundedness of the stochastic model
Consider a complete probability space which is denoted by (Ω, I c , P ) and with filtration represented by {F t } t∈R . The initial set F 0 contains all F null sets satisfying the conditions -right continuous and increasing. Let U (t) = (x (t) , y (t) , z (t) , w (t)). Its norm follows In addition, let C 2,1 1 R 4 × (0, ∞) : R + be the families of all non-negative function V (U, t) defined on R 4 × (0, ∞), being twice differentiable in U and once in "t".

Theorem 5 For system
Proof: By Ito's formula, (13)-(16) admits a positive solution in the sense of unique local on [0, τ e ] where explosion time is denoted τ e due to the local Lipschitz coefficients of the model.
Next, we shall prove that the equations (13)- (16) assume this solution in a global sense -meaning that τ e = ∞.
Let m 0 = 0 be sufficiently large for x (0) , y (0) , z (0) and w (0) lying with the interval 1 m0 , m 0 . For each integer m ≥ 0, define a sequence that is so-called stopping time as where we set ϕ = ∞ (ϕ is an empty set).

Computational methods
This section discusses well-known methods like stochastic Euler, stochastic Runge Kutta, and stochastic non-standard finite difference method with nonnegative initial conditions.

Stochastic Euler
This approach (also known as the forward Euler method) is a first-order numerical approximation for solving stochastic differential equations (SDEs). The discretization of the system (13)- (16) under the rules of this scheme is as follows: where h is any discretization parameter and n ≥ 0. Using the values of the parameters as presented in Table 2, the graphs of stochastic Euler's scheme for disease-free equilibrium (DFE) and endemic equilibrium (EE) plotted with MATLAB software are shown in Figures 2-5.

Stochastic Runge-Kutta method
This technique is also known as Runge Kutta of order 4th method. In numerical analysis, it is of the family of an implicit and explicit iterative process. In 1900 two German mathematicians Carl Runge and Wilhelm Kutta, developed this method. The stochastic Runge-Kutta scheme for equations (13)-(16) can be established according to the following steps.

First
Step

Second
Step

Third
Step

Fourth
Step Final Step where h is any discretization parameter and n ≥ 0. The graphs of the stochastic Runge Kutta scheme for disease-free equilibrium (DFE) and endemic equilibrium (EE) plotted with MATLAB software are presented in Figures 6-9 for the parameters from Table 2.

Stochastic NSFD
The NSFD method is also known as the implicitly driven explicit method. From equations (13)-(16), we established the computational scheme in the following manner. Equation (13) yields x n+1 = x n + hρ Equation (14) gives y n+1 = y n + h βmx n z n N + βpx n w n N + σ1βmx n z n N ∆B n + σ2βpx n w n N ∆B n 1 + hθ + hµ .

Equation (15) leads to
Equation (16) produces where h is any discretization parameter and n ≥ 0.

Stability analysis
Theorem 6 Theorem: For any n ≥ 0, the proposed SNSFD method is stable if the eigenvalues of the system (33)-(37) lie in the unit circle for R 0 < 1.
Proof: Consider the right-hand sides of the system of equations (33)-(37) as functions A, B, C and D with the assumption that ∆B n = 0, as follows: The elements of the Jacobian matrix are and we have At the disease-free equilibrium (DFE-L 1 ) = ρ µ , 0, 0, 0 , we get Hence, using the Mathematica software, all the eigenvalues of the above Jacobian matrix were found to lie within the unit circle for R 0 < 1. Thus, the system (33)-(37) is stable.

Results comparison
Here, a comparison of the stochastic NSFD method with other stochastic numerical methods was carried out. We can easily conclude that some of the examined stochastic numerical approaches conditionally converge or diverge with larger time step values by looking at the numerical experimentations. However, the NSFD method remains convergent for all time step sizes.

Concluding remarks
Stochastic modelling is a reliable and efficient technique to handle highly nonlinear problems close to nature, and it may be understood as the extension of deterministic modelling. In this paper, we applied the parametric perturbation technique to the formulated model. We focused on the NSFD scheme that is dynamically consistent and positively bounded. Contrary to the assertions in some literature that most epidemiologists used Euler Maruyama, stochastic Euler, and stochastic Runge Kutta algorithms, we have demonstrated that the novelty of the stochastic NSFD scheme is also satisfactory. We have shown that the other techniques are conditionally convergent and even fail in the long run. However, the investigated system unveils the NSFD method in a stochastic form as always being convergent. Furthermore, in the future, we will extend this idea to other types of modelling, including spatial-temporal, fractional, fractal fractional and delay problems of dynamical systems.