Analysis of Structure-Preserving Discrete Models for Predator-Prey Systems with Anomalous Diffusion

: In this work, we investigate numerically a system of partial differential equations that describes the interactions between populations of predators and preys. The system considers the effects of anomalous diffusion and generalized Michaelis–Menten-type reactions. For the sake of generality, we consider an extended form of that system in various spatial dimensions and propose two ﬁnite-difference methods to approximate its solutions. Both methodologies are presented in alternative forms to facilitate their analyses and computer implementations. We show that both schemes are structure-preserving techniques, in the sense that they can keep the positive and bounded character of the computational approximations. This is in agreement with the relevant solutions of the original population model. Moreover, we prove rigorously that the schemes are consistent discretizations of the generalized continuous model and that they are stable and convergent. The methodologies were implemented efﬁciently using MATLAB. Some computer simulations are provided for illustration purposes. In particular, we use our schemes in the investigation of complex patterns in some two- and three-dimensional predator–prey systems with anomalous diffusion.


Introduction
The investigation of the interactions between populations of predators and preys in nature is a highly transited topic of research in applied mathematics currently. Indeed, this area of research has proven to be extremely fruitful in view of the wide range of possible scenarios that merit investigation. As examples, we can mention studies that report on the modeling and analysis of predator-prey models with disease in the prey [1], the analysis of stochastic systems with modified Leslie-Gower and Holling-type schemes [2], the dynamic behaviors of Lotka-Volterra predator-prey models that incorporate predator cannibalism [3], the analysis of diffusive predator-prey systems with Michaelis-Menten-type predator harvesting [4], synthetic Escherichia coli predator-prey ecosystems [5], the analytical investigation of stage-structured predator-prey models depending on maturation delay and death rate [6], and non-autonomous ratio dependent models with Holling-type functional response with temporal delay [7], among other interesting topics [8].
It is worth pointing out that the investigation has focused mainly on the analytical aspects of the problem [9]. However, the literature reports on various numerical methods that have been designed explicitly to solve efficiently various predator-prey systems. For instance, there are studies establish the existence, uniqueness, positivity, and boundedness of the solutions for both methods. Section 5 is devoted to establishing the numerical properties of the schemes, including the consistency, stability, and convergence. Some numerical applications are provided in Section 6, and we close this manuscript with some concluding observations.

Preliminaries
Throughout, assume that a i and b i are real numbers such that a i < b i , for each i = 1, 2, 3. Let T > 0 represent a fixed time. We define Ω = (a 1 , b 1 ) × (a 2 , b 2 ) × (a 3 , b 3 ) and Ω T = Ω × (0, T) and use Ω and Ω T to represent the closures of Ω and Ω T , respectively. In this manuscript, u : Ω T → R and v : Ω T → R represent sufficiently smooth functions, and define x = (x 1 , x 2 , x 3 ) ∈ Ω.
Definition 1 (Podlubny [42]). Assume that f : R → R, and suppose that n ∈ N ∪ {0} and α ∈ R are such that n − 1 < α < n. If it exists, we introduce the fractional derivative in the sense of Riesz of the function f of order α at the point x ∈ R as: (1) Here, the Gamma function is given by: Definition 2. Suppose that u : Ω T → R; assume that α > −1; and let n ∈ Z satisfy n − 1 < α ≤ n. If they exist, the Riesz space-fractional derivatives of the function u of order α with respect to x 1 , x 2 , and x 3 at the point (x, t) ∈ Ω T are respectively defined by: For the remainder of this paper, we will let a, c, d, D 1 , and D 2 be positive; suppose that b ∈ R, and let α, β ∈ R be such that 1 < α ≤ 2 and 1 < β ≤ 2. Let φ u : Ω → R and φ v : Ω → R be two functions that physically describe the initial conditions for populations of prey and predator, respectively. Under these conventions, the problem under investigation is the predator-prey model with Allee effects and diffusion of fractional order, which is described by: Convey that u = u(x, t) and v = v(x, t) for simplicity. The model (6) is a Michaelis-Menten-type reaction-diffusion predator-prey system where the diffusion is anomalous. Here, u(x, t) and v(x, t) represent the normalized densities of the prey and the predator, respectively, at the point x ∈ Ω and time t ≥ 0. The relative constant a is the intrinsic rate of growth of the prey; b ∈ (−1, 1) is the Allee effect; c ∈ (0, 1] denotes the rate of the energy rate from the prey to the predator; and d is the relative rate of death of the predator population. Meanwhile, D 1 and D 2 are non-negative constants that represent the speed of individual movements of u and v, respectively [43]. Notice that we can rewrite the system (6) in generalized form as: where the function F depends on u, while G u and G v depend on both u and v. It is easy to see that the system (7) reduces to the population model (6) when F, G u , and G v have the following expressions with u, v ∈ R + ∪ {0}: (7) reduces to the well-known Lotka-Volterra system.
We recall the following definition from the literature. It will be an essential tool to provide consistent discretizations of the general fractional problem (7).
(a) The following iterative formulas hold: Lemma 2 (Wang et al. [44]). Let 0 < α ≤ 2 and α = 1, and suppose that f ∈ C 5 (R) is a function whose derivatives up to order five are all integrable. For almost all x ∈ R,

Numerical Models
The purpose of this section is to propose two different methods based on finite-differences to approximate the solutions of (7). For the sake of convenience, we consider only the two-dimensional form of (7), which reads as follows: It is worth pointing out that an analysis of the three-dimensional model is also feasible, though it would require additional nomenclature. We preferred to carry out the full description and analysis in the two-dimensional case for the sake of a better explanation.
Agree that I p = {1, 2, . . . , p} andĪ p = I p ∪ {0}, for all p ∈ N. Let M, N, K ∈ N, and introduce uniform partitions of [a 1 , b 1 ] and [a 2 , b 2 ], respectively, denoted by: Obviously, x 1,m = a 1 + h x 1 m and x 2,n = a 2 + h x 2 n for each m ∈Ī M and n ∈Ī N . In this case, the partition norms in the x 1 and In a similar fashion, we fix a (not necessarily uniform) partition for the interval [0, T], which will be represented by: For each k ∈Ī K−1 , we let τ k = t k+1 − t k . Numerically, we define u k m,n and v k m,n , respectively, as the approximations to the analytical solutions u and v of (14) at the point (x 1,m , x 2,n , t k ) for each m ∈Ī M , n ∈Ī N , and k ∈Ī K . Definition 4. Let α ∈ (0, 1) ∪ (1,2]. Define the discrete linear operators: By Lemma 2, the discrete operators introduced above provide approximations of second order to the Riesz spatial derivatives of the functions u and v, with respect to x 1 and x 2 at the point (x 1,m , x 2,n , t k ). For the remainder of this work and without loss of generality, we assume that the partition of the interval [0, T] is uniform, in which case τ k = τ ∈ R + , for each k ∈Ī K−1 . This assumption will be imposed only for the sake of convenience in the use of our notation.

Explicit Method
We present here an explicit scheme to approximate the solutions of (14). In the first stage, we introduce additional discrete operators to describe the scheme. The nomenclature presented in Section 2 will be observed throughout this section.

Definition 5.
Let w be any of u or v. For each m ∈Ī M , n ∈Ī N , and k ∈ I K−1 , we introduce the standard linear operators: Recall now that the operator (19) yields a first-order estimate of the partial derivative of w with respect to time at (x 1,m , x 2,n , t), for each m ∈Ī M , n ∈Ī N , and k ∈ I K−1 . Substituting the differential operators of (7) for their finite-difference approximations, we obtain the following discrete model to approximate the solutions of (14): Here, the difference equations are valid for each m ∈Ī M , n ∈Ī N , and k ∈ I K−1 .
Substituting then the expressions of the discrete operators into (20), we obtain an alternative representation of our finite-difference scheme. More precisely, let: for each i = 1, 2 and α, β ∈ (1, 2]. After some algebraic operations, it is easy to check that the recursive equations of (20) can be rewritten equivalently as: where: Let represent matrix transposition. Observe then that (22) can be represented in an equivalent vector form. Indeed, let k ∈Ī K and j ∈Ī M , and define the (N + 1)-dimensional vectors: With these conventions, we introduce the (M + 1) × (N + 1)-dimensional vectors: where ⊕ represents the vector operation of juxtaposition.
The following are all matrices of dimension (N + 1) × (N + 1), for each m ∈ I M and k ∈ I K−1 : For each k ∈ I K−1 , let A k u and A k v be block matrices of sizes [(M + 1) × (N + 1)] × [(M + 1) × (N + 1)], which are defined respectively by: With this notation, the vector representation of (22) is given by the iterative system:

Implicit Method
The purpose of this section is to introduce a Crank-Nicolson-type technique to approximate the solutions of (14). In the first stage, we define some discrete operators used to design our implicit finite-difference scheme. In the present section, we will observe the notation presented previously. The purpose is to provide various equivalent representations of the Crank-Nicolson scheme, which will be mathematically useful. Definition 6. Let w be any of u or v. For each m ∈Ī M , n ∈Ī N , and k ∈ I K−1 , we introduce the discrete linear operators: Remember that the discrete operators introduced in the previous definition yield second-order estimates of the temporal partial derivative of w at the point (x 1,m , x 2,n , t) and the exact value of w at that point, respectively. With this notation, a second finite-difference methodology to calculate the solutions of (14) is provided by the implicit system: for each m ∈Ī M , n ∈Ī N , and k ∈ I K−1 .
As in the case of the explicit method, an equivalent implicit representation of (42) is readily at hand. Indeed, after some algebraic simplifications and convenient manipulations, the difference equations of the system (42) can be rewritten as: where: Let k ∈Ī K and j ∈Ī M , and define the (N + 1)-dimensional vectors: Define then: Additionally, define the matrices B m,i , C k m , F m,i , and G k m as in Section 3.1, but using now the constants (44)- (49). Next, define the matrices A k u and A k v through the expressions of A k u and A k v in (39), using the new constants (44)- (49). On the other hand, we will agree that I is the identity matrix of dimension (N + 1) × (N + 1), and set H k u and E k v be the block matrices of dimension (M + 1) × (N + 1), given by: With this nomenclature, the matrix representation of (43) is given by: (61)

Structural Properties
The present section is devoted to establishing the main structural properties of the schemes (22) and (43). More precisely, we show the existence and uniqueness of the solutions of both methods. Moreover, we prove that the methods preserve the positive and bounded character of the solutions under appropriate conditions on the parameters.

Definition 7.
A real matrix A is nonnegative if all of its components are nonnegative. In such a case, we use the nomenclature A ≥ 0. If ρ ∈ R, then A is said to be bounded from above by ρ if every component of A is at most equal to ρ. This will be denoted by A ≤ ρ. In the case that ρ > 0, then we write 0 ≤ A ≤ ρ to denote that the conditions A ≥ 0 and A ≤ ρ are both satisfied.
In the following, we will suppose that the functions F, G u , and G v are bounded. As a consequence, there exist constants s 1 , s 2 , and s 3 ∈ R + with the properties that: for each m ∈Ī M , n ∈Ī N , and k ∈Ī K . Moreover, in the following, we will let s = max{s 1 , s 2 , s 3 }. Using these conventions, the following result establishes the main structural properties of the explicit method. It is worth pointing out that the existence and uniqueness of solutions are obviously inherent properties of this scheme in light of its explicit character.
Proof. To prove that u k+1 ≥ 0 and v k+1 ≥ 0, we only need to show that A k u ≥ 0 and A k v ≥ 0. Notice that the off-diagonal elements of A k u are equal to b m,i (for some m, i ∈Ī M and i = m) or equal to c n,i (for some n, i ∈Ī N and i = n) or zero.
Similarly, we can establish that c n,i > 0. On the other hand, the elements in the diagonal are of the form a k m,n (for some m ∈Ī M and n ∈Ī N ). Using (63), one obtains that: It follows that A k u ≥ 0, and the proof that A k v is established analogously. As a consequence, u k+1 ≥ 0 and v k+1 ≥ 0. To show that the approximations u k and v k are bounded, we define the vector e of dimension (M + 1)(N + 1) × (M + 1)(N + 1) with all elements equal to one. We will prove that z k+1 u = e − u k+1 > 0 and that z k+1 v = e − v k+1 > 0. Substituting z k+1 u into the first equation of (40), we readily obtain that z k+1 u = e − u k+1 = e − A k u u k , and we only need to show now that e − u k+1 > 0. Notice that e − u k+1 is a vector whose components are of the form: for m ∈Ī M and n ∈Ī N . By hypothesis and (65), we have: This implies that z k+1 u > 0 or, equivalently, that u k+1 ≤ 1. In a similar fashion, we can readily establish the inequality v k+1 ≤ 1. We conclude that 0 ≤ u k+1 ≤ 1 and 0 ≤ v k+1 ≤ 1 are satisfied.
We turn our attention to the structural properties of the scheme (42). To that end, the concept of the Minkowski matrix and its properties will be of utmost importance.

Definition 8.
If all the off-diagonal entries of a real matrix A are non-positive, then A is called a Z-matrix.

Definition 9. A real square matrix A is a Minkowski matrix if:
(i) A is a Z-matrix, (ii) all the diagonal components of A are positive, and (iii) A is strictly diagonally dominant.
Minkowski matrices are invertible, and their inverses are positive matrices [45]. This property will be exploited in our results.
Similarly, it is easy to see that c n,i < 0. On the other hand, the entries in the diagonal of A k u take the form a k m,n (for m ∈Ī M and n ∈Ī N ). Lemma 1(a) and the hypotheses show that a k m,n > 1 − aτF(u k m,n ) > 0. The strict diagonal dominance of A k u follows from Lemma 1, the hypotheses, and the fact that the following inequalities hold for each m ∈Ī M and each n ∈Ī N : It follows that A k u is a Minkowski matrix. That A k v is also a Minkowski matrix is proven similarly.
Theorem 2 (Existence and uniqueness). Let k ∈ I K−1 , and suppose that u k > 0 and v k > 0. If aτF(u k m,n ) < 1, cτG v (u k m,n , v k m,n ) < 1, and G u (u k m,n , v k m,n ) ≥ 0 for each m ∈Ī M and n ∈Ī N , then the recursive system (42) has a unique solution.
Proof. By hypothesis and Lemma 3, the matrices A k u and A k v are Minkowski matrices, so nonsingular. It follows that the equations A k u u k+1 = E k u u k−1 and A k v v k+1 = E k v v k−1 have unique solutions.
Theorem 3 (Positivity and boundedness). Let k ∈ I K−1 . Suppose that F is a bounded function over [0, 1] and that G u and G v are positive and bounded on the hold, then the solution of (42) is such that 0 ≤ u k+1 ≤ 1 and 0 ≤ v k+1 ≤ 1.
Proof. By hypothesis and Lemma 3, the matrices A k u and A k v are Minkowski matrices. To show that u k+1 ≥ 0 and v k+1 ≥ 0, use the identities of (61) to obtain: where (A k u ) −1 > 0 and (A k v ) −1 > 0. We need to show now that E k u and E k v are nonnegative. Note that the off-diagonal elements of E k u are either of the form −b m,i (with m, i ∈Ī M and i = m), or −c n,i (with n, i ∈Ī N and i = n), or zeros. As in the proof of Lemma 3, it follows that b m,i < 0 and c n,i < 0. Thus, −b m,i > 0 and −c n,i > 0. In turn, the elements in the diagonal are of the form 2 − a k m,n (for some m ∈Ī M and n ∈Ī N ). Observe then that (71) implies that: This shows that E k u > 0, and we can prove similarly that E k v > 0. It follows that the approximations u k+1 and v k+1 are nonnegative, and it only remains to establish the boundedness. Equivalently, we will prove that z k+1 u = e − u k+1 > 0 and z k+1 v = e − v k+1 > 0. Substituting z k+1 u into (61), we obtain: where z k+1 . It suffices to prove that A k u e − E k u u k−1 > 0, but the components of this vector are of the form: for m ∈Ī M and n ∈Ī N . By hypothesis and (73), it follows that: This implies that A k u e − E k u u k−1 > 0, which means that z k+1 u > 0. The proof that z k+1 v > 0 can be provided in an entirely analogous way. Finally, we conclude that 0 ≤ u k+1 ≤ 1 and 0 ≤ v k+1 ≤ 1.

Numerical Properties
The aim in this section is to prove the most important numerical features of the schemes (20) and (42). For the remainder of this section, we will use the following continuous functionals: Throughout, we employ the symbol h to represent the vector (h x 1 , h x 2 ).

Definition 10.
We define the infinity norm · ∞ : R q → R as: If E is a real matrix of size q × q, then its infinity norm is given by: For the remainder, if A is a real square matrix, then ρ(A) will represent its spectral radius.
Lemma 4 (Tian and Huang [46]). Let M = (m ij ) be an M-matrix, and let N = (n ij ) be a nonnegative matrix of the same size as M. If M is strictly diagonally dominant by rows, then ρ(M −1 N) satisfies:

Explicit Method
We will establish now the main numerical properties of the scheme (20). To that end, let us introduce the following discrete functionals, for each m ∈Ī M , n ∈Ī N , and k ∈ I K−1 : Theorem 4 (Consistency). Let u, v ∈ C 5 (Ω T ). If τ < 1, then there are constants C > 0 and C > 0 that are independent of τ, h x 1 , and h x 2 , with the property that for all m ∈Ī M , n ∈Ī N , and k ∈ I K−1 , |Lu k m,n − Lu k m,n | ≤ C(τ + h 2 ) and |Lv k m,n − Lv k m,n | ≤ C (τ + h 2 ).
Proof. We will apply standard arguments using Taylor's theorem and the formula (13) to prove the first inequality of (88). Since u ∈ C 5 (Ω T ), there exist constants C 1 , C 2,i > 0 for i ∈ {1, 2} that are independent of τ, h x 1 , and h x 2 , such that: for each m ∈Ī M , n ∈Ī N and k ∈ I K−1 . Defining C = max {C 1 , C 2,1 , C 2,2 } and applying the triangle inequality, we readily check that the first inequality of (88) holds. In a similar fashion, we can also establish the validity of the second inequality.
Theorem 5 (Nonlinear stability). Let (u k ) K k=0 , (v k ) K k=0 and (u k ) K k=0 , (v k ) K k=0 be two sets of positive and bounded solutions of (20) corresponding to the initial conditions (φ 0 u , φ 0 v ) and (φ 0 u , φ 0 v ), respectively. If the matrices A k u and A k u are identical to some constant matrix A k 1 and the matrices A k v and A k v are identical to some constant matrix A k 2 for each k ∈Ī K , then there exist constants C 1 , C 2 > 0 such that: Proof. We will only establish the first inequality from (91). Define k = u k − u k , for each k ∈Ī K−1 . By hypothesis, we have that k+1 = u k+1 − u k+1 = A k 1 (u k − u k ) = A k 1 k . Notice that recursion readily shows that k+1 = (A k 1 A k−1 1 · · · A 1 1 A 0 1 ) 0 . Taking the infinity norm in both sides, we obtain: The conclusion of this result readily follows when we take The proof for the second inequality is analogous.
Proof. Using the hypotheses, it is easy to see that for any m ∈Ī M and n ∈Ī N , By Lemma 4, we conclude that ρ(A k u ) < 1, and the inequality ρ(A k v ) < 1 is proven in similar fashion. As a consequence, we conclude that the scheme (20) is linearly stable, as desired.

Implicit Method
We turn our attention to the theoretical analysis of the scheme (42). Firstly, we establish the accuracy properties of that method. To this end, for each m ∈Ī M , n ∈Ī N , and k ∈ I K−1 , we define the discrete functionals Lu k m,n and Lv k m,n as: Lv k m,n = δ (1) t v k m,n − cG v (u k m,n , v k m,n )µ (1) t v k m,n + dµ Theorem 7 (Consistency). Let u, v ∈ C 5 (Ω T ). If τ < 1, then there exist constants C > 0 and C > 0 that are independent of τ, h x 1 and h x 2 such that for all m ∈Ī M , n ∈Ī N , and k ∈ I K−1 , |Lu k m,n − Lu k m,n | ≤ C(τ 2 + h 2 ) and |Lv k m,n − Lv k m,n | ≤ C (τ 2 + h 2 ).
Proof. By the regularity of u, there are constants C 1 , C 2 , C 3 , C 4,i > 0 for i = {1, 2}, such that: for each m ∈Ī M , n ∈Ī N , and k ∈ I K−1 . Let C = max {C 1 , C 2 , C 3 , C 4,1 , C 4,2 }, and use the triangle inequality to establish the first relation of (96). The proof of the second inequality is analogous.
Lemma 5 (Chen et al. [47]). Suppose that A is a matrix of size m × m and real components, which satisfies: Lemma 6. Let k ∈ I k−1 , and suppose that 0 ≤ u k ≤ 1 and 0 ≤ v k ≤ 1. If (73) and (74) are satisfied, Proof. According to Lemma 5, it suffices to show that A k u and A k v satisfy the inequality (101). Notice that the off-diagonal elements of A k u are equal to b m,i (for some m, i ∈Ī M and i = m), or equal to c n,i (for some n, i ∈Ī N and i = n), or zero. On the other hand, the elements in the diagonal are of the form a k m,n (for some m ∈Ī M and n ∈Ī N ). Using the inequality (73), we observe that: for each m ∈Ī M and each n ∈Ī N . Thus, the inequality (101) is satisfied for each row of A k u . In a similar fashion, we can prove that the inequality (101) holds for each row of A k v .

Proof.
We have already proven that E k u ≥ 0 in Theorem 3, so we only need to show that E k u ∞ < 1. Let m ∈Ī M and n ∈Ī N , and observe that the inequality (73) yields: Then, the sum of the all elements of each row of the matrix E k u is less than one. We conclude that E k u ∞ < 1. In a similar way, we can readily establish that E k v ≥ 0 and E k v ∞ < 1.
In our next results, we will establish the nonlinear and the linear stability of the method (42). It is worth pointing out that the study of convergence will be tackled in Theorem 10. An easy variation in the proof of that theorem readily establishes the linear convergence of the explicit scheme.
Theorem 8 (Nonlinear stability). Let (u k ) K k=0 , (v k ) K k=0 , and (u k ) K k=0 , (v k ) K k=0 be positive and bounded solutions of (42) corresponding to the initial conditions (φ 0 Suppose that the inequalities (71)-(74) hold for each k ∈Ī K . If the matrices A k Proof. Observe that (104) obviously holds if k ∈ {0, 1}, so assume that it is true for k ∈ {1, . . . , K − 1}. Using Lemmas 6 and 7, we have that: The conclusion of this result follows now by induction. In an analogous fashion, we may readily show that (105) holds.

Proof. We will use Lemma 4 to prove that ρ((A
Using the hypotheses, if m ∈Ī M and n ∈Ī N , then: Notice that ρ((A k u ) −1 E k u ) < 1 holds since every quotient is less than one. In a similar fashion, we can readily prove that ρ((A k v ) −1 E k v ) < 1. As a consequence, we conclude that (42) is linearly stable.
Theorem 10 (Convergence). Suppose that u, v ∈ C 5 (Ω T ) are positive and bounded solutions of (14). Let τ < 1, and suppose that (u k ) K k=0 and (v k ) K k=0 are positive and bounded solutions of (43). If (71)-(74) are satisfied for all k ∈ I K , then there are constants κ u , κ v ∈ R + that are independent of τ and h, such that: Proof. Define k = u k − u k , for each k ∈ I K . Since the exact and the numerical solutions coincide on the initial data, then 0 ∞ = 1 ∞ = 0. Using Theorem 7 and Lemmas 6 and 7, we obtain: which yields that k+1 ∞ − k−1 ∞ ≤ τC(τ 2 + h 2 ) for all k ∈ I K−1 . The conclusion follows now if we let κ u = TC. In a similar way, we can prove that there exists a constant κ v satisfying (109). As a conclusion, the scheme (43) is quadratically convergent.

Applications
The present section will be devoted to providing some computer simulations to illustrate the applicability of the finite difference schemes proposed in this work. In the first stage, we must point out that the explicit scheme (20) is obviously more suitable than the implicit model (42) when investigating two-and three-dimensional regimes [48,49]. To implement it efficiently, we notice firstly that the first equation of (20) can be rewritten as: where: for each m ∈ I M , n ∈ I N , and k ∈ I K−1 . Let W be the real matrix of size (M + 1) × (N + 1) whose entry at the row m ∈ I M and column n ∈ I N is W m,n = W(u k m,n , v k m,n ). Similarly, U k will represent the matrix of the same size as W such that U k m,n = u k m,n . Notice now that: Here, * represents the usual operation of matrix multiplication, and for each q ∈ N, H (α) q represents the real symmetric matrix of size (q + 1) × (q + 1) defined by: From this, it is easy to see now that the set of equations (111) can be rewritten equivalently in N . Summarizing, the explicit scheme (20) can be expressed alternatively in matrix form as: In this formula, we convey that Φ u m,n = φ 0 u,m,n and Φ v m,n = φ 0 v,m,n , for each (m, n) ∈ I M × I N . Moreover, we let Z be the real matrix of size (M + 1) × (N + 1) whose entry at the position (m, n) ∈ I M × I N is defined by Z(u k m,n , v k m,n ), where: Likewise, V k denotes the matrix of the same size as X with the property that V k m,n = v k m,n , for each (m, n) ∈ I M × I N and k ∈ I K . For the sake of convenience, we have included a MATLAB implementation of the scheme (116) in the Appendix. It is worth pointing out that the scheme is actually very general, not only in the sense that it accounts for different fractional differentiation orders, but also because the functions W and X therein are arbitrary.
The following examples will make use of variants of the MATLAB code provided in Appendix A. For the sake of illustration, the models considered will be in two spatial dimensions. Example 1. In this example, we let Ω = (0, 100) × (0, 100) and T = 3000. Throughout, we will consider the following diffusion-reaction system, defined for each (x, t) ∈ Ω: Here, the constants a, c, d, D 1 , and D 2 are positive, and we will define: It is easy to see that this point is a stationary solution of (118), which is a model that describes the spatio-temporal dynamics of a predator-prey system without Allee effects. Obviously, this system is a particular form of the more general model (7). To approximate solutions of the present system of equations, we will let φ v (x) be any sample from a normally distributed random variable with the mean equal to v * and the standard deviation equal to 0.01. Meanwhile, φ u will be equal to zero on all B, except on a central square at the middle of B, on which it will be constantly equal to 0.2. Let a = 0.8, c = 0.3, d = 0.1, D 1 = 0.01, and D 2 = 0.6. Figure 1 provides snapshots of the solution u in (118) at (a) t = 0, (b) t = 160, (c) t = 290, (d) t = 500, (e) t = 1010, and (f) t = 3000, using α = β = 2. The graphs exhibit the presence of complex patterns, in agreement with the results obtained in [43]. In those results, we agreed that x = x 1 and y = x 2 . For illustration purposes, Figures 2 and 3 provide similar results for the cases when α = β = 1.6 and α = β = 1.2, respectively. Turing patterns appear in those instances also, in spite of the fractional nature of those cases. In our simulations, we used our implementation of (20) shown in Appendix A, with τ = 0.02 and h x 1 = h x 2 = 1/3. Moreover, we imposed homogeneous Neumann conditions on the boundary of B.
It is worth pointing out that the literature lacks an analytical apparatus that justifies the presence of the Turing patterns in the fractional scenarios of Figures 2 and 3. In that sense, the explicit numerical methodology reported in the present paper may be a reliable tool to confirm any analytical results on the fully fractional form of (118).
In our last example, we will tackle the three-dimensional scenario. To that end, the code of Appendix A had to be adapted to the three-dimensional scenario, and parallel programming was needed in order to speed up the computer time.

Example 2.
We considered the three-dimensional form of problem (118) with the same model parameters, spatial domain Ω = (0, 100) × (0, 100) × (0, 100) ⊆ R 3 , τ = 0.02, and h x 1 = h x 2 = h x 3 = 1. Under these circumstances, Figure 4 shows snapshots of the solution u of the three-dimensional problem (118) at (a) t = 0, (b) t = 160, (c) t = 290, (d) t = 500, (e) t = 1010, and (f) t = 3000, letting α = β = 1.6. The graphs exhibit the presence of complex three-dimensional patterns that have not been investigated in the literature. For convenience, Figure 5 shows x-, y-and z-cross sections of the approximate solution u at the same times. In this case, the graphs exhibit the presence of two-dimensional patterns on the sides of the cube B. In turn, Figures 6 and 7 show similar results for the case when α = β = 1.2. Again, the presence of three-dimensional and two-dimensional patterns is obvious from the graphs. We performed more simulations (not included in this work to avoid redundancy) using various values of α and β. The results have shown the presence of three-dimensional complex patterns in all the cases considered. , and (f) t = 3000. We let φ v (x) be a random sample from a normally distributed random variable with the mean equal to v * and the standard deviation equal to 0.01, and φ u is the function depicted in (a). The approximations were calculated using our implementation of (20) shown in Appendix A, with τ = 0.02 and h x 1 = h x 2 = 1/3. , and (f) t = 3000. We let φ v (x) be a random sample from a normally distributed random variable with the mean equal to v * and the standard deviation equal to 0.01, and φ u is the function depicted in (a). The approximations were calculated using our implementation of (20) shown in Appendix A, with τ = 0.02 and h x 1 = h x 2 = 1/3. , and (f) t = 3000. We let φ v (x) be a random sample from a normally distributed random variable with the mean equal to v * and the standard deviation equal to 0.01, and φ u is the function depicted in (a). The approximations were calculated using our implementation of (20) shown in Appendix A, with τ = 0.02 and h x 1 = h x 2 = 1/3.  Our last example provides numerical evidence that the emerging Turing patterns preserve their share independent of the discretization step.  Figure 8 shows the results of our simulations considering various values of the spatial partition norms. Throughout, we let h = h x 1 = h x 2 .
The graphs correspond to the values (a) h = 1/3, (b) h = 1/4, h = 1/5, and h = 1/8. A similar pattern shape was obtained in all cases. These results provide numerical evidence that the type of Turing pattern is independent of the discretization spatial step-size. We performed similar experiments considering different temporal steps. The simulations are not reproduced to avoid redundancy, but they prove that the emerging patterns are also independent of τ.

Conclusions
In this manuscript, we studied computationally a system of two diffusive partial differential equations with coupled nonlinear reactions in generalized forms. Our system considered the presence of anomalous diffusion in multiple spatial dimensions along with suitable initial data. The model generalized various particular systems from the physical sciences, including the diffusive systems, which describe the interaction between populations of predators and preys with the Michaelis-Menten-type reaction. The system was discretized following a finite-difference approach, and two schemes were proposed to approximate the solutions. The two numerical models were rigorously analyzed to elucidate their structural and numerical properties.
As the main structural results, we established the existence and uniqueness of the numerical solutions. Moreover, we proved that the schemes were both capable of preserving the positivity and boundedness of approximations [50,51]. As in various other examples available in the literature [31,[52][53][54], this was in perfect agreement with the fact that the relevant solutions of normalized population models were positive and bounded. Numerically, we proved rigorously that the schemes were consistent, stable, and convergent. Some simulations were provided in this work to illustrate the performance of the schemes. In particular, we showed the capability of one of the schemes to be applied on the investigation of Turing patterns in anomalously diffusive systems describing predator-prey interactions. To that end, a fast computational implementation in MATLAB of one of the schemes was employed. Of course, the present approach may be applied to other different scenarios [55,56].
At the closure of this manuscript, it is interesting to point out that the two discretizations of the system (7) were capable of preserving the anomalous diffusion rate. Indeed, note that those rates were equal to D 1 and D 2 for each of the partial differential equations of the continuous model (7). On the other hand, in light of Lemma 2, it follows that: Similarly, it is easy to check that: It follows that our spatial discretizations are capable of preserving the anomalous diffusion rate. Funding: The first author would like to acknowledge the financial support of the National Council for Science and Technology of Mexico (CONACYT). The second author acknowledges financial support from CONACYT through Grant A1-S-45928.

Acknowledgments:
The authors wish to thank the guest editors for their kind invitation to submit a paper to the Special Issue of Mathematics MDPI on "Computational Mathematics and Neural Systems". They also wish to thank the anonymous reviewers for their comments and criticisms. All of their comments were taken into account in the revised version of the paper, resulting in a substantial improvement with respect to the original submission. Finally, the authors also wish to thank the Benemérica Universidad Autónoma de Aguascalientes.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. MATLAB Code
The following is a basic implementation in MATLAB of the computational model (116). This code was employed to solve Problem (118). Suitable variations of the code were considered in order to produce the simulations of that example. The parallelization of the scheme is not provided, though it is worth pointing out that the task is straightforward.