A Fractional Approach to a Computational Eco-Epidemiological Model with Holling Type-II Functional Response

Abstract

Eco-epidemiological can be considered as a significant combination of two research fields of computational biology and epidemiology. These problems mainly take ecological systems into account of the impact of epidemiological factors. In this paper, we examine the chaotic nature of a computational system related to the spread of disease into a specific environment involving a novel differential operator called the Atangana–Baleanu fractional derivative. To approximate the solutions of this fractional system, an efficient numerical method is adopted. The numerical method is an implicit approximate method that can provide very suitable numerical approximations for fractional problems due to symmetry. Symmetry is one of the distinguishing features of this technique compared to other methods in the literature. Through considering different choices of parameters in the model, several meaningful numerical simulations are presented. It is clear that hiring a new derivative operator greatly increases the flexibility of the model in describing the different scenarios in the model. The results of this paper can be very useful help for decision-makers to describe the situation related to the problem, in a more efficient way, and control the epidemic.

1. Introduction

Nature is full of interactions between different species of living beings to provide food, shelter, and other essential needs. In many cases, it is necessary to gain a better understanding of these interactions to preserve different species of living organisms in nature. In recent years, the use of mathematical modeling in the study of problems in computational biology has attracted the attention of many experts in the fields of mathematics and biology [1,2,3,4,5,6,7,8,9,10,11,12,13]. As a result, and based on these important applications, many effective techniques have been proposed in solving mathematical models arising from such problems [14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39]. One of the most interesting aspects of such problems is examining cases where the environment has been affected by an infectious disease. Such models are referred to as eco-epidemiological models and have also been explored with an aim of disease control. Certainly, in this case, the spread of disease among the populations considered in the model can have a direct impact on the survival or extinction of species in the problem.

Nowadays, with the spread of infectious and contagious diseases throughout the natural world and their direct impact on food interactions in ecosystems, many models have been presented and studied. For example, taking the weak Allee effect and harvesting in prey population into account, the following nonlinear delay system of equations has been proposed [40]:

$$\begin{array}{c}\frac{d\mathsf{S}\left(t\right)}{dt}=\mathsf{S}\left(t\right)\left[(1-\mathsf{S}\left(t\right)-\mathsf{I}\left(t\right))\frac{\mathsf{S}\left(t\right)}{\mathsf{S}\left(t\right)+\theta }-\beta \mathsf{I}\left(t\right)-q_{1}E\right],\hfill \\ \frac{d\mathsf{I}\left(t\right)}{dt}=\mathsf{I}\left(t\right)[\beta \mathsf{S}\left(t\right)-a\mathsf{P}\left(t\right)-\mu -q_{2}E,\hfill \\ \frac{d\mathsf{P}\left(t\right)}{dt}=\alpha \mathsf{P}(t-\tau )\mathsf{I}(t-\tau )-d\mathsf{P}\left(t\right),\hfill \end{array}$$

(1)

where existing parameters in the model are introduced in Reference [40].

Furthermore, an eco-epidemiological model has been investigated in [41] under the combined influence of strong-Allee parameter and competition coefficients, given by

$$\begin{array}{c}\frac{d\mathsf{S}\left(t\right)}{dt}=\left[r-b\mathsf{S}\left(t\right)-c\mathsf{I}\left(t\right)\right]\mathsf{S}\left(t\right)\left(1-\frac{\theta +f}{\mathsf{S}\left(t\right)+f}\right)-\frac{\beta \mathsf{I}\left(t\right)\mathsf{S}\left(t\right)}{a+\mathsf{S}\left(t\right)},\hfill \\ \frac{d\mathsf{I}\left(t\right)}{dt}=\frac{\beta \mathsf{I}\left(t\right)\mathsf{S}\left(t\right)}{a+\mathsf{S}\left(t\right)}-\frac{{\alpha }_1\mathsf{I}\left(t\right)\mathsf{P}\left(t\right)}{d+\mathsf{I}\left(t\right)}-\mu \mathsf{I}\left(t\right),\hfill \\ \frac{d\mathsf{P}\left(t\right)}{dt}=\frac{{\gamma }_1\mathsf{I}\left(t\right)\mathsf{P}\left(t\right)}{d+\mathsf{I}\left(t\right)}-m\mathsf{P}\left(t\right).\hfill \end{array}$$

(2)

To get some more information about existing parameters in Model (2), please refer to the work in [41].

The authors of [42] have considered that an infectious disease has spread to the prey and predator populations in the environment. In their work, they have used the following nonlinear system to describe these interactions:

$$\begin{array}{c}\frac{d\mathsf{S}\left(t\right)}{dt}=\left[r-b\mathsf{S}\left(t\right)-c\mathsf{I}\left(t\right)-\frac{\beta \mathsf{I}\left(t\right)}{a+\mathsf{S}\left(t\right)}-\frac{{\alpha }_1\mathsf{P}\left(t\right)}{e+\mathsf{S}\left(t\right)}\right]\mathsf{S}\left(t\right),\hfill \\ \frac{d\mathsf{I}\left(t\right)}{dt}=\left[\frac{\beta \mathsf{S}\left(t\right)}{a+\mathsf{S}\left(t\right)}-\frac{{\alpha }_2\mathsf{P}\left(t\right)}{d+\mathsf{I}\left(t\right)}-\mu \right]\mathsf{I}\left(t\right),\hfill \\ \frac{d\mathsf{P}\left(t\right)}{dt}=\left[\left(\frac{c_1{\alpha }_1\mathsf{S}\left(t\right)}{e+\mathsf{S}\left(t\right)}+\frac{c_2{\alpha }_2\mathsf{I}\left(t\right)}{d+\mathsf{I}\left(t\right)}\right)\frac{\mathsf{P}\left(t\right)}{\theta +\mathsf{P}\left(t\right)}-m\right]\mathsf{P}\left(t\right).\hfill \end{array}$$

(3)

In all these models, $\mathsf{S}\left(t\right)$ and $\mathsf{I}\left(t\right)$ denote the susceptible and the diseased prey, respectively. Further, $\mathsf{P}\left(t\right)$ has been employed to describe the predator population. More explanations about existing parameters in Model (3) can be found in [42].

In recent years, the use of fractional derivative operators in the modeling various problems in mathematics, physics, and engineering has increased significantly. For example, In [43], a SIR epidemic model with Crowley–Martin type functional response and Holling type-II treatment rate is investigated. Moreover, a comparison based on newly defined fractional derivative operators which are called as Caputo–Fabrizio (CF) and Atangana–Baleanu (AB) has been presented in [44]. A model for HIV-1 primary infection with treatment in fractional order along with its Global dynamics is considered in [45]. The FitzHugh–Nagumo (F–N) model has been studied in [46] using the nerve impulses process. In [47], several chaotic systems, involving Atangana–Baleanu fractional derivative operator with interesting behaviors are presented. In [48], a fractal-fractional differentiation for the modeling and mathematical analysis of nonlinear diarrhea transmission dynamics under the use of real data has been investigated. In [49], hybrid mathematical models of new strains and co-infection in Caputo, Caputo–Fabrizio, and Atangana–Baleanu are studied. Further, a mathematical model in describing the transmission of Nipah virus within a targeted population has been proposed in [50], More examples can be found in [13,23,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65].

One of the main reasons for this increase in popularity is perhaps the fact that these operators are equipped with the concept of memory. In fact, to determine the value of the derivative in these fractional operators, it is necessary to collect all the previous information of those phenomena in the calculations from the beginning to that specific time situation. This valuable feature of operators plays a very vital and influential role in biological modeling, where the current behavior of the variables in the models is very much influenced by their overall behavior in the past. Perhaps it is for these reasons that many research activities in this field have been carried out with the help of these fractional operators [66,67,68,69]. One of the most prominent applications that has received much attention recently is the use of this tool in the mathematical modeling of problems related to Covid-19 disease [70,71,72,73,74,75,76]. As this disease has had a devastating effect on the lives of many of us around the world, we still need to study some more related models.

Motivated by the above contribution, this paper considers a newly proposed eco-epidemiological system. The key point in this model is to use the Atangana–Baleanu fractional derivative [77,78,79,80]. Similar to the works in [81,82,83,84,85,86,87,88,89], in some special situations, in this model chaotic behaviors [90,91,92,93,94,95,96,97,98,99] occur. This property indicates the high sensitivity of the model with respect to some values for its parameters. The study of these conditions has wide applications in the study of natural ecosystems.

The rest parts of this contribution are managed as follows. First, some necessary prerequisites, including definitions and properties related to fractional operators are presented in Section 2. The main proposed system of the article is proposed then investigated in Section 3. In Section 4, some theoretical aspects of the model, such as the calculation of equilibrium points of the model along with their stability, and the existence and uniqueness of the solution for the model are studied. Then, we present an efficient numerical algorithm to solve the problem in Section 5. Moreover, simulation results and discussion corresponding to the employed numerical technique are outlined in Section 6. Finally, some of the achievements and conclusions related to the results are drawn in Section 7.

2. An Overview of Fractional Calculus

In what follows, we outline some important and essential prerequisites in fractional calculus.

Definition 1. 

For a given function $\omega \left(t\right)$, calculus operators of the Caputo type are, respectively, defined as [100]

$${}^{\mathsf{C}}{\mathcal{D}}^{ϵ}\omega \left(t\right)=\frac{1}{\mathsf{\Gamma }(k-ϵ)}{\int }_0^t{(t-\rho )}^{m-ϵ-1}{\omega }^{\left(k\right)}\left(\rho \right)d\rho ,k-1<ϵ\le k,k\in \mathbb{N},$$

(4)

and

$$\begin{array}{cc}\hfill {}^{\mathsf{C}}{\mathcal{I}}^{ϵ}\omega \left(t\right)& =\frac{1}{\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\omega \left(\rho \right)d\rho .0<ϵ<1.\hfill \end{array}$$

(5)

Definition 2. 

For a given function $\omega \left(t\right)$, calculus operators of the Caputo–Fabrizio ($\mathsf{CF}$) type are, respectively, defined as [101]

$${}^{\mathsf{CF}}{\mathcal{D}}^{ϵ}\omega \left(t\right)=\frac{\mathsf{M}\left(ϵ\right)}{k-ϵ}{\int }_0^t\frac{d^k\omega \left(\rho \right)}{d{\rho }^k}exp\left[-\frac{ϵ}{n-ϵ}(t-\rho )\right]d\rho ,k-1<ϵ\le k$$

(6)

$${}^{\mathsf{CF}}{\mathcal{I}}^{ϵ}\omega \left(t\right)=\frac{2(1-ϵ)}{(2-ϵ)\mathsf{M}\left(ϵ\right)}\omega \left(t\right)+\frac{2ϵ}{(2-ϵ)\mathsf{M}\left(ϵ\right)}{\int }_0^t\omega \left(\rho \right)d\rho ,$$

(7)

where

$$\mathsf{M}\left(ϵ\right)=\frac{2}{2-ϵ}.$$

(8)

Definition 3. 

For a given function $\omega \left(t\right)$, calculus operators of the Atangana–Baleanu ($\mathsf{AB}$) type are, respectively, defined as [78]

$${}^{\mathsf{AB}}{\mathcal{D}}^{ϵ}\omega \left(t\right)=\frac{\mathsf{\Theta }\left(ϵ\right)}{1-ϵ}{\int }_0^t{\mathcal{E}}_{ϵ}\left[-\frac{ϵ}{1-ϵ}(t-\rho )\right]{\omega }^{\prime }\left(\rho \right)d\rho ,ϵ\in (0,1),$$

(9)

$${}^{\mathsf{AB}}{\mathcal{I}}^{ϵ}\omega \left(t\right)=\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\omega \left(t\right)+\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_{t_0}^t{(t-\mu )}^{ϵ-1}\omega \left(\rho \right)d\rho ,$$

(10)

where $\mathsf{\Theta }\left(ϵ\right)=1-ϵ+\frac{ϵ}{\mathsf{\Gamma }\left(ϵ\right)}$, and ${\mathcal{E}}_{ϵ}\left(\rho \right)$ is the Mittag–Leffler function given by [102]

$${\mathcal{E}}_{ϵ}\left(t\right)=\sum _{j=0}^{\infty }\frac{t^j}{\mathsf{\Gamma }(1+jϵ)}.$$

(11)

Moreover, the combination of the $\mathsf{AB}$ derivative and differential operators yields

$$\begin{array}{c}\hfill {}^{\mathsf{AB}}{\mathcal{I}}^{ϵ}\left({}^{\mathsf{AB}}{\mathcal{D}}^{ϵ}\omega \left(t\right)\right)=\omega \left(t\right)-\omega \left(0\right).\end{array}$$

(12)

3. The Model Formulation

In this paper, we consider a eco-epidemiological system that models the interaction between three species: ${\mathcal{X}}_1\left(t\right),{\mathcal{X}}_2\left(t\right)$, and ${\mathcal{X}}_3\left(t\right)$. The description of the interactions between the variables in this system is done using the following nonlinear system [103]:

$$\begin{array}{c}\hfill \begin{array}{cc}\frac{d{\mathcal{X}}_1\left(t\right)}{dt}\hfill & =r{\mathcal{X}}_1\left(t\right)-b{\mathcal{X}}_1^2\left(t\right)-c{\mathcal{X}}_1\left(t\right){\mathcal{X}}_2\left(t\right)-\frac{{\alpha }_1{\mathcal{X}}_1\left(t\right){\mathcal{X}}_3\left(t\right)}{e+{\mathcal{X}}_1\left(t\right)}-\frac{\beta {\mathcal{X}}_1\left(t\right){\mathcal{X}}_2\left(t\right)}{a+{\mathcal{X}}_1\left(t\right)},\hfill \\ \frac{d{\mathcal{X}}_2\left(t\right)}{dt}\hfill & =\frac{\beta {\mathcal{X}}_1\left(t\right){\mathcal{X}}_2\left(t\right)}{a+{\mathcal{X}}_1\left(t\right)}-\frac{{\alpha }_2{\mathcal{X}}_3\left(t\right){\mathcal{X}}_2\left(t\right)}{d+{\mathcal{X}}_2\left(t\right)}-\mu {\mathcal{X}}_2\left(t\right),\hfill \\ \frac{d{\mathcal{X}}_3\left(t\right)}{dt}\hfill & =\frac{{\alpha }_1{c}_1{\mathcal{X}}_1\left(t\right){\mathcal{X}}_3\left(t\right)}{e+{\mathcal{X}}_1\left(t\right)}+\frac{c_2{\alpha }_2{\mathcal{X}}_3\left(t\right){\mathcal{X}}_2\left(t\right)}{d+{\mathcal{X}}_2\left(t\right)}-m{\mathcal{X}}_3\left(t\right).\hfill \end{array}\end{array}$$

(13)

Three state variables exist in this model, including susceptible prey ${\mathcal{X}}_1\left(t\right)$ and infected prey ${\mathcal{X}}_2\left(t\right)$, and predator population ${\mathcal{X}}_2\left(t\right)$. These variables have considered the saturated incidence $\frac{\beta {\mathcal{X}}_1\left(t\right){\mathcal{X}}_2\left(t\right)}{a+{\mathcal{X}}_1\left(t\right)}$, where $\beta$ is the force of infection and a is saturation constant. Furthermore, r represents the growth rate of ${\mathcal{X}}_1\left(t\right)$, and b and c are the intra-class and inter-class competition coefficients, respectively. The predation rate of a susceptible prey and infected prey are respectively denoted by ${\alpha }_1$ and ${\alpha }_2$. Both $e,d$ are the half-saturation constants, and $c_1$ and $c_2$ are the conversion efficiency of the predator on susceptible and infected prey, respectively. Moreover, $\mu$ is the death rate of infected prey. Finally, m is used to explain the natural mortality rate of the predator population. More details about the model can be found in [103].

In order to benefit from the valuable features of fractional differential calculus in the studied model (13), let us replace the normal derivative in the model with the $\mathsf{AB}$ fractional derivative ${}^{\mathsf{AB}}{\mathcal{D}}^{ϵ}$. Subsequently, the new structure for the model (13) is considered as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}{}^{\mathsf{AB}}{\mathcal{D}}^{ϵ}{\mathcal{X}}_1\left(t\right)\hfill & =r{\mathcal{X}}_1\left(t\right)-b{{\mathcal{X}}_1}^2\left(t\right)-c{\mathcal{X}}_1\left(t\right){\mathcal{X}}_2\left(t\right)-\frac{{\alpha }_1{\mathcal{X}}_1\left(t\right){\mathcal{X}}_3\left(t\right)}{e+{\mathcal{X}}_1\left(t\right)}-\frac{\beta {\mathcal{X}}_1\left(t\right){\mathcal{X}}_2\left(t\right)}{a+{\mathcal{X}}_1\left(t\right)},\hfill \\ {}^{\mathsf{AB}}{\mathcal{D}}^{ϵ}{\mathcal{X}}_2\left(t\right)\hfill & =\frac{\beta {\mathcal{X}}_1\left(t\right){\mathcal{X}}_2\left(t\right)}{a+{\mathcal{X}}_1\left(t\right)}-\frac{{\alpha }_2{\mathcal{X}}_3\left(t\right){\mathcal{X}}_2\left(t\right)}{d+{\mathcal{X}}_2\left(t\right)}-\mu {\mathcal{X}}_2\left(t\right),\hfill \\ {}^{\mathsf{AB}}{\mathcal{D}}^{ϵ}{\mathcal{X}}_3\left(t\right)\hfill & =\frac{{\alpha }_1{c}_1{\mathcal{X}}_1\left(t\right){\mathcal{X}}_3\left(t\right)}{e+{\mathcal{X}}_1\left(t\right)}+\frac{c_2{\alpha }_2{\mathcal{X}}_3\left(t\right){\mathcal{X}}_2\left(t\right)}{d+{\mathcal{X}}_2\left(t\right)}-m{\mathcal{X}}_3\left(t\right),\hfill \end{array}\end{array}$$

(14)

subject to initial conditions $({\mathcal{X}}_1\left(t\right),{\mathcal{X}}_2\left(t\right),{\mathcal{X}}_3\left(t\right)){|}_{t=0}=({\mathcal{X}}_{1,0},{\mathcal{X}}_{2,0},{\mathcal{X}}_{3,0})\ge 0$.

4. Mathematical Analysis for Model

In this section, some theoretical features corresponding to the model (14) are explored.

4.1. The Equilibrium Points

The model considered in (14) includes the following equilibrium points:

• ${\mathcal{W}}_1=\left(0,0,0\right),$ that explains the trivial equilibrium point.
• ${\mathcal{W}}_2=\left(0,\frac{dm}{{\alpha }_2{c}_2-m},\frac{-\mu d{c}_2}{{\alpha }_2{c}_2-m}\right),$ which is meaningless from a biological point of view.
• ${\mathcal{W}}_3=\left(\frac{r}{b},0,0\right)$ that explains the existence of susceptible prey-only situation.
• ${\mathcal{W}}_4=\left(\frac{em}{{\alpha }_1{c}_1-m},0,\frac{e{c}_1({\alpha }_1{c}_{1}r-bem-mr)}{{\left({\alpha }_1{c}_1-m\right)}^2}\right),$ that explains the coexistence of suspected prey and predator populations situation. A sufficient condition for the existence of these points is that we have

  $${\alpha }_1{c}_1>m,\frac{{\alpha }_1{c}_{1}r}{be+r}>m.$$
• ${\mathcal{W}}_5=\left(\frac{a\mu }{\beta -\mu },-\frac{a(ab\mu -\beta r+\mu r)}{(\beta -\mu )(ac+\beta -\mu )},0\right),$ that implies the coexistence of suspected and infected prey populations in the model. A sufficient condition for the existence of these points is that we have

  $$\mu <\beta ,\frac{b\mu a}{\beta -\mu }<r.$$

The local stability for the equilibrium points of the model (14) at ${\mathcal{W}}_i=({\mathcal{X}}_1^{*},{\mathcal{X}}_2^{*},{\mathcal{X}}_3^{*})$ can be investigated using the following Jacobian matrix

$$\left[\begin{array}{ccc}-c{\mathcal{X}}_2^{*}-2b{\mathcal{X}}_1^{*}+r-\frac{\beta {\mathcal{X}}_2^{*}}{a+{\mathcal{X}}_1^{*}}+\frac{\beta {\mathcal{X}}_1^{*}{\mathcal{X}}_2^{*}}{{\left(a+{\mathcal{X}}_1^{*}\right)}^2}-\frac{{\alpha }_1{\mathcal{X}}_3^{*}}{e+{\mathcal{X}}_1^{*}}+\frac{{\alpha }_1{\mathcal{X}}_1^{*}{\mathcal{X}}_3^{*}}{{\left(e+{\mathcal{X}}_1^{*}\right)}^2}& -c{\mathcal{X}}_1^{*}-\frac{\beta {\mathcal{X}}_1^{*}}{a+{\mathcal{X}}_1^{*}}& -\frac{{\alpha }_1{\mathcal{X}}_1^{*}}{e+{\mathcal{X}}_1^{*}}\\ \frac{\beta {\mathcal{X}}_2^{*}}{a+{\mathcal{X}}_1^{*}}-\frac{\beta {\mathcal{X}}_1^{*}{\mathcal{X}}_2^{*}}{{\left(a+{\mathcal{X}}_1^{*}\right)}^2}& \frac{\beta {\mathcal{X}}_1^{*}}{a+{\mathcal{X}}_1^{*}}-\frac{{\alpha }_2{\mathcal{X}}_3^{*}}{d+{\mathcal{X}}_2^{*}}+\frac{{\alpha }_2{\mathcal{X}}_2^{*}{\mathcal{X}}_3^{*}}{{\left(d+{\mathcal{X}}_2^{*}\right)}^2}-\mu & -\frac{{\alpha }_2{\mathcal{X}}_2^{*}}{d+{\mathcal{X}}_2^{*}}\\ \frac{{\alpha }_1{c}_1{\mathcal{X}}_3^{*}}{e+{\mathcal{X}}_1^{*}}-\frac{c_1{\alpha }_1{\mathcal{X}}_1^{*}{\mathcal{X}}_3^{*}}{{\left(e+{\mathcal{X}}_1^{*}\right)}^2}& \frac{{\alpha }_2{c}_2{\mathcal{X}}_3^{*}}{d+{\mathcal{X}}_2^{*}}-\frac{{\alpha }_2{c}_2{\mathcal{X}}_2^{*}{\mathcal{X}}_3^{*}}{{\left(d+{\mathcal{X}}_2^{*}\right)}^2}& \frac{c_1{\alpha }_1{\mathcal{X}}_1^{*}}{e+{\mathcal{X}}_1^{*}}+\frac{{\alpha }_2{c}_2{\mathcal{X}}_2^{*}}{d+{\mathcal{X}}_2^{*}}-m\end{array}\right].$$

(15)

4.2. The Existence of the Solution

Here, we aim to confirm that the fractional model always possesses a solution. To this end, after incorporating the $\mathsf{AB}$ integral operator (10) on both sides of equations of the system (14) one gets

$$\begin{array}{c}\hfill {\mathcal{X}}_1\left(t\right)-{\mathcal{X}}_1\left(0\right)=\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}{\mathcal{Q}}_1\left(\mathbf{Y}\left(\mathbf{t}\right)\right)+\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}{\mathcal{Q}}_1\left(\mathbf{Y}\left(\rho \right)\right)d\rho ,\\ \hfill {\mathcal{X}}_2\left(t\right)-{\mathcal{X}}_2\left(0\right)=\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}{\mathcal{Q}}_2\left(\mathbf{Y}\left(\mathbf{t}\right)\right)+\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}{\mathcal{Q}}_2\left(\mathbf{Y}\left(\rho \right)\right)d\rho ,\\ \hfill {\mathcal{X}}_3\left(t\right)-{\mathcal{X}}_3\left(0\right)=\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}{\mathcal{Q}}_3\left(\mathbf{Y}\left(\mathbf{t}\right)\right)+\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}{\mathcal{Q}}_3\left(\mathbf{Y}\left(\rho \right)\right)d\rho ,\end{array}$$

(16)

where $\mathbf{Y}\left(\mathbf{t}\right)=\left[{\mathcal{X}}_1\left(t\right),{\mathcal{X}}_2\left(t\right),{\mathcal{X}}_3\left(t\right)\right]$, and also we define

$$\begin{array}{c}\hfill \begin{array}{cc}{\mathcal{Q}}_1\left(\mathbf{Y}\left(\mathbf{t}\right)\right)\hfill & =r{\mathcal{X}}_1\left(t\right)-b{{\mathcal{X}}_1}^2\left(t\right)-c{\mathcal{X}}_1\left(t\right){\mathcal{X}}_2\left(t\right)-\frac{{\alpha }_1{\mathcal{X}}_1\left(t\right){\mathcal{X}}_3\left(t\right)}{e+{\mathcal{X}}_1\left(t\right)}-\frac{\beta {\mathcal{X}}_1\left(t\right){\mathcal{X}}_2\left(t\right)}{a+{\mathcal{X}}_1\left(t\right)},\hfill \\ {\mathcal{Q}}_2\left(\mathbf{Y}\left(\mathbf{t}\right)\right)\hfill & =\frac{\beta {\mathcal{X}}_1\left(t\right){\mathcal{X}}_2\left(t\right)}{a+{\mathcal{X}}_1\left(t\right)}-\frac{{\alpha }_2{\mathcal{X}}_3\left(t\right){\mathcal{X}}_2\left(t\right)}{d+{\mathcal{X}}_2\left(t\right)}-\mu {\mathcal{X}}_2\left(t\right),\hfill \\ {\mathcal{Q}}_3\left(\mathbf{Y}\left(\mathbf{t}\right)\right)\hfill & =\frac{{\alpha }_1{c}_1{\mathcal{X}}_1\left(t\right){\mathcal{X}}_3\left(t\right)}{e+{\mathcal{X}}_1\left(t\right)}+\frac{c_2{\alpha }_2{\mathcal{X}}_3\left(t\right){\mathcal{X}}_2\left(t\right)}{d+{\mathcal{X}}_2\left(t\right)}-m{\mathcal{X}}_3\left(t\right).\hfill \end{array}\end{array}$$

(17)

Defining $\mathbf{N}\left(\mathbf{Y}\left(\mathbf{t}\right)\right)=\left[{\mathcal{Q}}_1\left(\mathbf{Y}\left(\mathbf{t}\right)\right),{\mathcal{Q}}_2\left(\mathbf{Y}\left(\mathbf{t}\right)\right),{\mathcal{Q}}_3\left(\mathbf{Y}\left(\mathbf{t}\right)\right)\right]$, and moreover ${\mathbf{Y}}_0=$$\left[{\mathcal{X}}_1\left(0\right),{\mathcal{X}}_2\left(0\right),{\mathcal{X}}_3\left(0\right)\right].$ Using these assumptions, Equation (16) can be rewritten as

$$\mathbf{Y}\left(\mathbf{t}\right)-{\mathbf{Y}}_0=\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{N}\left(\mathbf{Y}\left(\mathbf{t}\right)\right)+\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\mathbf{N}\left(\mathbf{Y}\left(\rho \right)\right)d\rho .$$

(18)

Now, inspired by the (18) and starting from ${\mathbf{Y}}_0\left(t\right)={\mathbf{Y}}_0$, we define the following iterative scheme

$${\mathbf{Y}}_n\left(t\right)-{\mathbf{Y}}_0=\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{N}\left({\mathbf{Y}}_{n-1}\left(t\right)\right)+\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\mathbf{N}\left({\mathbf{Y}}_{n-1}\left(\rho \right)\right)d\rho .$$

(19)

Considering Equation (19), we have

$${\mathbf{Y}}_n\left(t\right)-{\mathbf{Y}}_{n-1}\left(t\right)=\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\left[\mathbf{N}\left({\mathbf{Y}}_{n-1}\left(t\right)\right)-\mathbf{N}\left({\mathbf{Y}}_{n-2}\left(t\right)\right)\right]+\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\left[\mathbf{N}\left({\mathbf{Y}}_{n-1}\left(\rho \right)\right)-\mathbf{N}\left({\mathbf{Y}}_{n-2}\left(\rho \right)\right)\right]d\rho .$$

(20)

In this position, we define ${\varsigma }_n\left(t\right)={\mathbf{Y}}_n\left(t\right)-{\mathbf{Y}}_{n-1}\left(t\right)$. Then, it follows that

$${\mathbf{Y}}_n\left(t\right)=\sum _{i=0}^n{\varsigma }_i\left(t\right).$$

(21)

As a result, one gets

$$\begin{array}{ccc}\hfill \parallel {\varsigma }_n\left(t\right)\parallel & =& \parallel {\mathbf{Y}}_n\left(t\right)-{\mathbf{Y}}_{n-1}\left(t\right)\parallel \hfill \\ \hfill \parallel {\varsigma }_n\left(t\right)\parallel & =& \parallel \frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\left[\mathbf{N}\left({\mathbf{Y}}_{n-1}\left(t\right)\right)-\mathbf{N}\left({\mathbf{Y}}_{n-2}\left(t\right)\right)\right]+\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\left[\mathbf{N}\left({\mathbf{Y}}_{n-1}\left(\rho \right)\right)-\mathbf{N}\left({\mathbf{Y}}_{n-2}\left(\rho \right)\right)\right]d\rho \parallel .\hfill \end{array}$$

Therefore, we have

$$\begin{array}{c}\hfill \parallel {\varsigma }_n\left(t\right)\parallel \le \frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\parallel \mathbf{N}\left({\mathbf{Y}}_{n-1}\left(t\right)\right)-\mathbf{N}\left({\mathbf{Y}}_{n-2}\left(t\right)\right)\parallel +\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\parallel \mathbf{N}\left({\mathbf{Y}}_{n-1}\left(\rho \right)\right)-\mathbf{N}\left({\mathbf{Y}}_{n-2}\left(\rho \right)\right)\parallel d\rho .\end{array}$$

Now, if nonlinear operator $\mathbf{N}$ satisfies the Lipshitz condition, then one obtains

$$\begin{array}{c}\hfill \parallel {\varsigma }_n\left(t\right)\parallel \le \frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{L}\parallel {\mathbf{Y}}_{n-1}\left(t\right)-{\mathbf{Y}}_{n-2}\left(t\right)\parallel +\frac{ϵ\mathbf{L}}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\parallel {\mathbf{Y}}_{n-1}\left(t\right)-{\mathbf{Y}}_{n-2}\left(t\right)\parallel d\rho .\end{array}$$

Consequently, we derive the following inequality:

$$\begin{array}{c}\hfill \parallel {\varsigma }_n\left(t\right)\parallel \le \frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{L}\parallel {\varsigma }_{n-1}\left(t\right)\parallel +\frac{ϵ\mathbf{L}}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\parallel {\varsigma }_{n-1}\left(t\right)\parallel d\rho .\end{array}$$

Further, replacing $\parallel {\varsigma }_{n-1}\left(t\right)\parallel$ by its value, it reads

$$\begin{array}{c}\hfill \parallel {\varsigma }_n\left(t\right)\parallel \le {\left(\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{L}+\frac{ϵ\mathbf{L}{t}^{ϵ}}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }(ϵ+1)}\right)}^2\parallel {\varsigma }_{n-2}\left(t\right)\parallel .\end{array}$$

Furthermore, it reads

$$\begin{array}{c}\hfill \parallel {\varsigma }_n\left(t\right)\parallel \le {\left(\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{L}+\frac{ϵ\mathbf{L}{t}^{ϵ}}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }(ϵ+1)}\right)}^3\parallel {\varsigma }_{n-3}\left(t\right)\parallel .\end{array}$$

Moreover, finally, we obtain

$$\begin{array}{ccc}\hfill \parallel {\varsigma }_n\left(t\right)\parallel & \le & {\left(\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{L}+\frac{ϵ\mathbf{L}{t}^{ϵ}}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }(ϵ+1)}\right)}^n\parallel {\varsigma }_0\left(t\right)\parallel ,\hfill \\ & \le & {\left(\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}+\frac{ϵt^{ϵ}}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }(ϵ+1)}\right)}^n{\mathbf{L}}^n{max}_{t\in [0,T]}{\mathbf{Y}}_0\left(t\right).\hfill \end{array}$$

(22)

In this case, we examine the following definition:

$$\begin{array}{c}\hfill \mathbf{Y}\left(t\right)=\sum _{i=0}^n{\varsigma }_i\left(t\right).\end{array}$$

(23)

The general structure for $\mathbf{Y}\left(t\right)$ is then suggested as

$$\begin{array}{c}\hfill \mathbf{Y}\left(t\right)={\mathbf{Y}}_n\left(t\right)+{\mu }_n\left(t\right),\end{array}$$

(24)

where ${\mu }_n\left(t\right)\to 0$ when $n\left(t\right)\to \infty$. Thus,

$$\begin{array}{c}\hfill \mathbf{Y}\left(t\right)-{\mathbf{Y}}_n\left(t\right)=\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{N}(\mathbf{Y}\left(t\right)-{\mu }_n\left(t\right))+\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\mathbf{N}(\mathbf{Y}\left(\rho \right)-{\mu }_n\left(\rho \right))d\rho .\end{array}$$

(25)

Now, we can write

$$\begin{array}{ccc}\hfill & & \mathbf{Y}\left(t\right)-{\mathbf{Y}}_0-\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{N}(\mathbf{Y}\left(t\right)-{\mu }_n\left(t\right))-\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\mathbf{N}(\mathbf{Y}\left(\rho \right)-{\mu }_n\left(\rho \right))d\rho \hfill \\ & =& {\mu }_n\left(t\right)+\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\left[\mathbf{N}(\mathbf{Y}\left(t\right)-{\mu }_n\left(t\right))-\mathbf{N}\left(\mathbf{Y}\left(t\right)\right)\right]-\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\left[\mathbf{N}(\mathbf{Y}\left(\rho \right)-{\mu }_n\left(\rho \right))-\mathbf{N}\left(\mathbf{Y}\left(\rho \right)\right)\right]d\rho .\hfill \end{array}$$

Taking norm on both sides of the latter result, one gets

$$\begin{array}{ccc}\hfill & \parallel & \mathbf{Y}\left(t\right)-{\mathbf{Y}}_0\left(t\right)-\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{N}\left(\mathbf{Y}\left(t\right)\right)+\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\mathbf{N}(\mathbf{Y}\left(\rho \right)d\rho \parallel \hfill \\ \hfill & \le & \parallel {\mu }_n\left(t\right)\parallel +\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\parallel \mathbf{N}(\mathbf{Y}\left(\rho \right)-{\mu }_n\left(\rho \right))-\mathbf{N}\left(\mathbf{Y}\left(\rho \right)\right)\parallel \hfill \\ \hfill & & +\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\parallel \mathbf{N}(\mathbf{Y}\left(\rho \right)-{\mu }_n\left(\rho \right))-\mathbf{N}\left(\mathbf{Y}\left(\rho \right)\right)\parallel d\rho ,\hfill \\ \hfill & \le & \parallel {\mu }_n\left(t\right)\parallel +\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{L}\parallel {\mu }_{n-1}\left(t\right)\parallel ,+\frac{ϵt^{ϵ}}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }(ϵ+1)}\mathbf{L}\parallel {\mu }_{n-1}\left(t\right)\parallel .\hfill \end{array}$$

For large values n, the right side of the equation becomes zero, so one gets

$$\mathbf{Y}\left(\mathbf{t}\right)-{\mathbf{Y}}_0=\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{N}\left(\mathbf{Y}\left(\mathbf{t}\right)\right)+\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_0^t{(t-\rho )}^{ϵ-1}\mathbf{N}\left(\mathbf{Y}\left(\rho \right)\right)d\rho .$$

(26)

Moreover, this result provides clear proof of the uniqueness for the solution for the system.

4.3. The Uniqueness of the Solution

At this point in the article, we assume that the system has two solutions ${\mathbf{Y}}_1\left(t\right)$ and ${\mathbf{Y}}_2\left(t\right)$. Now, we obtain

$$\begin{array}{ccc}\hfill \parallel {\mathbf{Y}}_1\left(t\right)-{\mathbf{Y}}_2\left(t\right)\parallel & \le & \frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{L}\parallel {\mathbf{Y}}_1\left(t\right)-{\mathbf{Y}}_2\left(t\right)\parallel +\frac{ϵ\mathbf{L}{t}^{ϵ}}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }(ϵ+1)}\parallel {\mathbf{Y}}_1\left(t\right)-{\mathbf{Y}}_2\left(t\right)\parallel ,\hfill \\ & \le & \left(\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{L}+\frac{ϵ\mathbf{L}{t}^{ϵ}}{\mathsf{\Theta }(1+ϵ)\mathsf{\Gamma }\left(ϵ\right)}\right)\parallel {\mathbf{Y}}_1\left(t\right)-{\mathbf{Y}}_2\left(t\right)\parallel ,\hfill \\ & ⋮& \\ & \le & {\left(\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{L}+\frac{ϵ\mathbf{L}{t}^{ϵ}}{\mathsf{\Theta }(1+ϵ)\mathsf{\Gamma }\left(ϵ\right)}\right)}^n\parallel {\mathbf{Y}}_1\left(t\right)-{\mathbf{Y}}_2\left(t\right)\parallel .\hfill \end{array}$$

Then, if $\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{L}+\frac{ϵ\mathbf{L}{t}^{ϵ}}{\mathsf{\Theta }(1+ϵ)\mathsf{\Gamma }\left(ϵ\right)}<1$ holds then for $n\to \infty$, one obtains

$${\left(\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathbf{L}+\frac{ϵ\mathbf{L}{t}^{ϵ}}{\mathsf{\Theta }(1+ϵ)\mathsf{\Gamma }\left(ϵ\right)}\right)}^n\to 0.$$

Consequently, $\parallel {\mathbf{Y}}_1\left(t\right)-{\mathbf{Y}}_2\left(t\right)\parallel =0$ holds. Therefore, ${\mathbf{Y}}_1\left(t\right)={\mathbf{Y}}_2\left(t\right)$ is resulted.

5. An Approximate Approach to the Solution

The introduction of numerical methods has always been one of the consequences of presenting new definitions in the field of fractional differential calculus. In other words, each new definition for operators requires its numerical method. Each of these numerical methods has its own advantages, limitations, and requirements. In this article, we follow the numerical idea of the product integration (PI) rule [104].

To develop the numerical method, first, consider the following fractional system:

$${}^{\mathsf{AB}}{\mathcal{D}}^{ϵ}\mathsf{\Omega }\left(t\right)=\mathcal{N}(t,\mathsf{\Omega }).$$

(27)

The use of the integral operator introduced in (10) on both sides of (27) results in the following Volterra integral equation:

$$\mathsf{\Omega }\left(t\right)-\mathsf{\Omega }\left(t_0\right)=\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathcal{N}(t,\mathsf{\Omega }\left(t\right))+\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}{\int }_{t_0}^t{(t-\rho )}^{ϵ-1}\mathcal{N}(\rho ,\mathsf{\Omega }\left(\rho \right))d\rho .$$

(28)

Considering the time discretization of $t=t_n=t_0+n\hslash$ in (28) suggests

$$\mathsf{\Omega }\left(t_n\right)=\mathsf{\Omega }\left(t_0\right)+\frac{1-ϵ}{\mathsf{\Theta }\left(ϵ\right)}\mathcal{N}(t_n,\mathsf{\Omega }\left(t_n\right))+\frac{ϵ}{\mathsf{\Theta }\left(ϵ\right)\mathsf{\Gamma }\left(ϵ\right)}\sum _{i=0}^{n-1}{\int }_{t_i}^{t_{i+1}}{(t_n-\rho )}^{ϵ-1}\mathcal{N}(\rho ,\mathsf{\Omega }\left(\rho \right))d\rho .$$

(29)

Using the idea of linear interpolation, the function $\mathcal{N}(\rho ,\mathsf{\Omega }(\rho \left)\right)$ can be expanded as follows:

$$\mathcal{N}(\rho ,\mathsf{\Omega }\left(\rho \right))\approx \mathcal{N}(t_{i+1},{\mathsf{\Omega }}_{i+1})+\frac{\rho -t_{i+1}}{\hslash }\left(\mathcal{N}(t_{i+1},{\mathsf{\Omega }}_{i+1})-\mathcal{N}(t_i,{\mathsf{\Omega }}_i)\right),\rho \in [t_i,t_{i+1}],$$

(30)

where ${\mathsf{\Omega }}_i=\mathsf{\Omega }\left(t_i\right)$.

Taking the linear function (30) into account in the integrand in (29), and also by performing some required algebraic calculations, the following iterative scheme is obtained to approximate the $\mathsf{AB}$ fractional problem (27) as [79,105,106]

$$\begin{array}{cc}\hfill {\mathsf{\Omega }}_n& ={\mathsf{\Omega }}_0+\frac{ϵ\hslash }{\mathsf{\Theta }\left(ϵ\right)}\left({\nu }_n\mathcal{N}(t_0,{\mathsf{\Omega }}_0)+\sum _{i=1}^n{\zeta }_{n-i}\mathcal{N}(t_i,{\mathsf{\Omega }}_i)\right),n\ge 1,\hfill \end{array}$$

(31)

where

$$\begin{array}{cc}\hfill {\nu }_n& =\frac{{(-1+n)}^{ϵ+1}-n(-1+n-ϵ)}{\mathsf{\Gamma }(ϵ+2)},\hfill \\ \hfill {\zeta }_i& =\{\begin{array}{c}\frac{1-ϵ}{ϵ\hslash }+\frac{1}{\mathsf{\Gamma }(ϵ+2)},i=0,\hfill \\ \frac{{(i-1)}^{ϵ+1}-2i^{ϵ+1}+{(i+1)}^{ϵ+1}}{\mathsf{\Gamma }(ϵ+2)},i=1,2,\dots ,n-1.\hfill \end{array}\hfill \end{array}$$

(32)

The approximate method obtained in (31) and (32) can be efficiently employed to characterize the approximate solution to the model (14). In this case, it is obtained as

$$\begin{array}{cc}\hfill {{\mathcal{X}}_1}_n={{\mathcal{X}}_1}_0& +\frac{ϵ\hslash }{\mathsf{\Theta }\left(ϵ\right)}[{\nu }_n\left(r{\mathcal{X}}_{1,0}-b{\mathcal{X}}_{1,0}^2-c{\mathcal{X}}_{1,0}{\mathcal{X}}_{2,0}-\frac{{\alpha }_1{\mathcal{X}}_{1,0}{\mathcal{X}}_{3,0}}{e+{\mathcal{X}}_{1,0}}-\frac{\beta {\mathcal{X}}_{1,0}{\mathcal{X}}_{2,0}}{a+{\mathcal{X}}_{1,0}}\right)\hfill \\ & +\sum _{i=0}^n{\zeta }_{n-i}\left(r{\mathcal{X}}_{1,i}-b{\mathcal{X}}_{1,i}^2-c{\mathcal{X}}_{1,i}{\mathcal{X}}_{2,i}-\frac{{\alpha }_1{\mathcal{X}}_{1,i}{\mathcal{X}}_{3,i}}{e+{\mathcal{X}}_{1,i}}-\frac{\beta {\mathcal{X}}_{1,i}{\mathcal{X}}_{2,i}}{a+{\mathcal{X}}_{1,i}}\right)],\hfill \\ \hfill {{\mathcal{X}}_2}_n={{\mathcal{X}}_2}_0& +\frac{ϵ\hslash }{\mathsf{\Theta }\left(ϵ\right)}[{\nu }_n\left(\frac{\beta {\mathcal{X}}_{1,0}{\mathcal{X}}_{2,0}}{a+{\mathcal{X}}_{1,0}}-\frac{{\alpha }_2{\mathcal{X}}_{3,0}{\mathcal{X}}_{2,0}}{d+{\mathcal{X}}_{2,0}}-\mu {\mathcal{X}}_{2,0}\right)\hfill \\ & +\sum _{i=0}^n{\zeta }_{n-i}\left(\frac{\beta {\mathcal{X}}_{1,i}{\mathcal{X}}_{2,i}}{a+{\mathcal{X}}_{1,i}}-\frac{{\alpha }_2{\mathcal{X}}_{3,i}{\mathcal{X}}_{2,i}}{d+{\mathcal{X}}_{2,i}}-\mu {\mathcal{X}}_{2,i}\right)],\hfill \\ \hfill {{\mathcal{X}}_3}_n={{\mathcal{X}}_3}_0+& \frac{ϵ\hslash }{\mathsf{\Theta }\left(ϵ\right)}[{\nu }_n\left(\frac{{\alpha }_1{c}_1{\mathcal{X}}_{1,0}{\mathcal{X}}_{3,0}}{e+{\mathcal{X}}_{1,0}}+\frac{c_2{\alpha }_2{\mathcal{X}}_{3,0}{\mathcal{X}}_{2,0}}{d+{\mathcal{X}}_{2,0}}-m{\mathcal{X}}_{3,0}\right)\hfill \\ & +\sum _{i=0}^n{\zeta }_{n-i}\left(\frac{{\alpha }_1{c}_1{\mathcal{X}}_{1,i}{\mathcal{X}}_{3i}}{e+{\mathcal{X}}_{1,i}}+\frac{c_2{\alpha }_2{\mathcal{X}}_{3,i}{\mathcal{X}}_{2,i}}{d+{\mathcal{X}}_{2,i}}-m{\mathcal{X}}_{3,i}\right)].\hfill \end{array}$$

(33)

Applying these iterative forms will result in approximate solutions to the fractional problem we are considering in this paper. It is clear that these schemes are implicit forms and can be solved via existing efficient techniques, such as Newton’s method.

6. Discussion of Simulation Results

To solve the fractional system fractional, the iterative scheme introduced in (33) is used. The values for $ϵ$ that have been taken in the model are $0.85,0.874,0.898,0.922,0.946$, and $0.97$. Further, corresponding numerical simulations in the model are generated by taking

$$\begin{array}{c}\hfill c=0.01,b=1,r=1,\beta =0.6,a=0.363,{\alpha }_1=0.01,e=15,{\alpha }_2=0.05,d=0.5,\mu =0.4,c1=2,c2=1,m=0.01.\end{array}$$

(34)

Note that in all simulations, we have used the initial condition $\left({\mathcal{X}}_1\left(0\right),{\mathcal{X}}_2\left(0\right),{\mathcal{X}}_3\left(0\right)\right)=\left(0.93,0.25,1.2\right)$.

In Figure 1 and Figure 2, we have selected the value of $ϵ$ as $0.85,0.874,0.898,0.922,0.946$, and $0.97$ in the system (14). In these plots, it is clear that any value for the $ϵ$ causes a certain effect on the system, but eventually, the system tends to a single equilibrium point.

Now, we investigate the sensitivity of the model with respect to some of the parameters in the formulation presented in the system (14). First, we consider the values defined in (34) with $ϵ=0.95$ and examine the effect of $\beta$ on results via taking $0.57,0.58,0.59,0.6,0.61$, and $0.62$. Approximate solutions corresponding to these assumptions are shown in Figure 3 and Figure 4. In these plots, it is clear that as the amount of $\beta$ increased, the system tends to exhibit more malicious behavior. However, for smaller values, the system is stable and tends towards a certain equilibrium point.

In this section, we look at the effect of parameter $m_1$ on model results. To this end, the values of $m_1$ are taken as $0.01,0.018,0.026,0.034,0.042$, and $0.05$. In Figure 5 and Figure 6, we have depicted the acquired approximate results by taking these values. The results show that a change in this value for the parameter will change the type of system equilibrium point.

Moreover, we have examined the role of a through taking $0.3,0.34,$$0.38,0.42,0.46$, and $0.5$. Figure 7 and Figure 8 display corresponding to the obtained solutions of the system (14). Behaviors related to smaller values of a’s are more stable whilst for larger values of the parameter; a very severe oscillation behavior is evident in the system.

Then, we have taken the role of ${\alpha }_1$ on results via taking $0.01,0.02,0.03,0.04,0.05$, and $0.06$ in Figure 9 and Figure 10. In these graphs, it is clear that by increasing the value for ${\alpha }_1$, the behavior will be very chaotic and unstable.

Finally, we study the role of ${\alpha }_2$ on results through taking values $0.01,0.02,0.03,0.04,0.05$, and $0.06$ in Figure 11 and Figure 12.

The results obtained in these plots confirm that the system solution for each case converges to its different corresponding equilibrium points.

7. Conclusions

Eco-epidemiology can be considered as a meaningful combination of two research fields of ecology and epidemiology. These problems mainly take ecological systems into account epidemiological factors. By utilizing new modern mathematical tools, significant progress can be made in the modeling eco-epidemiological scenarios. By fully incorporating these elements, eco-epidemiology would be rooted in the investigation of the pathways by which biological and social experiences generate health and disease. In this contribution, the $\mathsf{AB}$ fractional derivative is employed to study some computational aspects of three species of the prey–predator model in mathematical biology. The concept of memory and symmetry are the main reasons for using these useful tools. According to the results obtained in this study, which is very compatible with the expected conditions, the method used in the article can be used to solve other problems in epidemiology. This can be considered as a direction for future research in describing other eco-epidemiological problems.