In this paper, a deterministic and stochastic fractional order model for lesser date moth (LDM) using mating disruption and natural enemies is proposed and analysed. The interaction between LDM larvae, fertilized LDM female, unfertilized LDM female, LDM male and the natural enemy is investigated. In order to clarify the characteristics of the proposed deterministic fractional order model, the analysis of existence, uniqueness, non-negativity and boundedness of the solutions of the proposed fractional-order model are examined. In addition, some sufficient conditions are obtained to ensure the local and global stability of equilibrium points. The occurrence of local bifurcation near the equilibrium points is investigated with the help of Sotomayor’s theorem. Numerical simulations are conducted to illustrate the properties of the proposed fractional order model with respect to the intrinsic growth rate of the LDM larvae, natural enemy’s mortality rate, predation rate, sex pheromone trap parameter, fractional order and environmental noise. The impact of mating disruption on lesser date moth is demonstrated. Also, a numerical approximation method is developed for the proposed stochastic fractional-order model.
Lesser date mothstochasticstabilitynatural enemiesSotomayor’s theoremmating disruptionIntroduction
The palm tree is considered one of the oldest major and basic crops in Southwest Asia, North Africa and many other places of the world. Palm trees are affected by many agricultural pests that cause significant losses to the palm trees and their fruits, as well as affecting the age and growth of the palm tree if it is left without control [1]. The LDM is one of the most dangerous pests of young palm trees and immature palm fruits. The damage is mainly due to the way in which the LDM larvae feed, as soon as they leave the egg until the date of their entry the virgin dwelling develops feeding and causing tunnels in the affected part of the palm in all directions and depths without any early signs showing the infection [2–5]. One of the most promising strategies for controlling LDM is the use of a mating disruption using the sex pheromone traps [5]. Natural enemies’ use to stop pest infestation has long been recommended as a clean and environmentally friendly way to protect crops. The main natural enemies that are used in agricultural pest control are larval predators. Goniozus swirskiana can be considered as one of the most important natural enemies that can attack the LDM [6,7]. In the real world, plant diseases models are always affected by the environmental noise. Thus, the stochastic models may be a more appropriate way of modelling agricultural pests in many circumstances [8]. Recently, fractional calculus has been applied to describe different mathematical models, and it has been shown to be more accurate in some cases compared to the classical models [9]. The main objective of this paper is to propose and analyse a deterministic and stochastic fractional order LDM model with mating disruption and sex pheromone traps taking into consideration the effect of the natural enemy on LDM.
The paper is arranged as follows: In Section 2, the deterministic mathematical model is described as well as the existence, uniqueness, non-negativity, and boundedness of the solutions of LDM system are verified. The local and global stability of equilibrium points of the LDM system is analyzed in Section 3. Using Sotomayor’s theorem, the local bifurcation conditions are derived in Section 4. In Section 5, we extend the deterministic fractional-order LDM model to the stochastic case and a numerical approximation method developed for the proposed stochastic fractional-order case. Some numerical examples are given in Section 6 to illustrate the theoretical findings. Finally, the discussion and conclusion are given in Section 7.
Dynamic of the Deterministic Fractional Order Model
Following [10–14], the model of lesser date moth with mating disruption and sex pheromone
where CDα is the Caputo fractional derivative of order α and 0<α<1. The detailed explanation of the variables and parameters for system (1) are listed in Tab. 1.
Biological description of parameters used in system (1)
Variables & parameters
Description
x
population density of LDM larvae
y
population density of unfertilized LDM female
f
population density of fertilized LDM female
m
population density of LDM male
z
population density of natural enemy
r
the LDM larvae intrinsic growth rate
k
the carrying capacity
α1
transfer rate from larvae to unfertilized LDM females
α2
transfer rate from unfertilized to fertilized LDM females
β
predation intensity between larvae and natural enemy
a
half saturation constant
μ1
mortality rate of LDM larvae
μ2
mortality rate of unfertilized LDM female
μ3
mortality rate of fertilized LDM female
μ4
mortality rate of LDM male
μ5
mortality rate of natural enemy
ϵ
proportion of larvae that emerge to unfertilized females
δ
multiple mating rate of fertilized LDM female
p
the maximal death rate by sex-pheromone trap
γ
the capture rate for male by sex pheromone trap
e
conversion coefficient for predation between LDM larvae and natural enemy.
Existence and Uniqueness
In this section, we investigate the existence and uniqueness of the solutions of the fractional order system (1) in the region Ω×(0,T] where
Ω={(x,y,f,m,z)∈R+5:max(|x|,|y|,|f|,|m|,|z|)≤φ},
Theorem 1.For eachX0=(x0,y0,f0,m0,z0)∈Ω, there exists a unique solutionX(t)∈Ωof the fractional order system (1), which is defined for allt≥0.
Proof. Define a mapping F(X)=(F1(X),F2(X),F3(X),F4(X),F5(X)), in which
Hence, F(X) satisfies the Lipschitz condition with respect to X. According to Cresson et al. [15], as F(X) locally Lipschitz. Then there exists unique local solution to the fractional order system (1).
Non-Negativity and Boundedness
The following results show the non-negativity of the solutions of the fractional order system (1). According to Cresson et al. [15], a model of the form dXdt=F(X) satisfies the positivity property if and only if for all i=1,….,5, Fi(X)≥0 for all X∈R+5 such that Xi=0. Thus, the solution of the integer-order model (1), with nonnegative initial conditions remains nonnegative. Also, the solution satisfies the Lipschitz condition, as stated in Theorem 1. By Theorem 5 and Theorem 6 in Cresson et al. [15], the solution of the fractional-order model (1) also satisfies the non-negativity. The boundedness of the solutions of model (1) are given in the following theorem.
Theorem 2.All the solutions of the fractional-order LDM Model (1) starting in R+5are uniformly bounded.
Proof. The approach of [16,17] is utilized. Let (x(t),y(t),f(t),m(t),z(t)) to be any solution of the system (1) with non-negative initial conditions. Let M(t)=x(t)+y(t)+f(t)+m(t)+z(t), then
where, ν<min{μ1,μ2,μ3,μ4,μ5}, thus, cαDαM(t)+νM≤rk4. In accordance with Lemma 9 in Choi et al. [18], it follows that, 0≤M(t)≤M(0)Eα(−νtα)+rk4νtαEα,α+1(−νtα),
where Eα is the Mittag-Leffler function. According to Lemma 5 and Corollary 6 in Choi et al. [18], it follows
0≤M(t)≤rk4ν,ast→∞.
Hence all the solutions of fractional-order LDM model (1) that start in R+5 are uniformly bounded in the region, H={(x,y,f,m,z)∈R+5:M(t)≤rk4ν+ξ,for\ anyξ>0}.
Equilibria and Stability
The LDM model (1) has the following three equilibrium points:
1) E0=(0,0,0,0,0), which always exists.
2) The free natural enemy equilibrium point E1=(x1,y1,f1,m1,0), where
The free natural enemy equilibrium point exists positively if R0>1, where, R0=rα1+μ1 is the basic offspring number obtained by using the next generation method [19].
3) The coexistence equilibrium point E2=(x2,y2,f2,m2,z2), where
4) The coexistence equilibrium point E2 exists if eβ>μ5 and R0>1+raμ5k(α1+μ1)(eβ−μ5).
The locally and globally asymptotically stable of equilibrium points of LDM system (1) are now investigated. The stability analysis of the equilibrium point E0=(0,0,0,0,0) is not considered because in this case all the population will go to extinction.
The stability of free natural enemy equilibrium point E1=(x1,y1,f1,m1,0) is investigated as follows.
Theorem 3.IfR0<1+raμ5k(α1+μ1)(eβ−μ5),then the equilibrium pointE1is locally asymptotically stable.
Proof. The first three eigenvalues of J(E1) are λ1=r−2rx1k−(α1+μ1)=(α1+μ1)(1−R0), λ2=−(γpy1+p+μ4) and λ3=eβx1a+x1−μ5 . The other roots are determined by
λ2+(α2+δ+μ2+μ3)λ+δμ2+μ3(α2+μ2)=0.
The roots of Eq. (3) have negative real parts. It can be observed that λ3<0, when eβx1a+x1<μ5 which equivalent R0<1+raμ5k(α1+μ1)(eβ−μ5). So, E1 is locally asymptotically stable if 1<R0<1+raμ5k(α1+μ1)(eβ−μ5).
Theorem 4.IfR0<1+raμ5k(α1+μ1)(eβ−μ5), then the natural enemy extinction equilibrium pointE1is globally asymptotically stable.
Proof. The following positive definite Lyapunov function is considered.
L1=eaa+x1(x−x1−x1lnxx1)+z.
By calculating the time derivative of L1 along the solution of system (1), one obtains,
Then cαDαL1≤0, when eβx1a+x1<μ5 which equivalent R0<1+raμ5k(α1+μ1)(eβ−μ5). According to generalized Lyapunov–Lasalle’s invariance principle [20], the natural enemy extinction equilibrium point E1 is globally asymptotically stable when R0<1+raμ5k(α1+μ1)(eβ−μ5).
The stability of the coexistence equilibrium point E2=(x2,y2,f2,m2,z2) is investigated as follows.
The first eigenvalue of J(E2) is λ1=−γpp+y2−μ4. The other eigenvalues are determined by
where B11=rx2k−βz2x2(a+x2)2 and B44=γpp+y2+μ4. Then, the proposition proposed in Matouk [21] can be used to determine the stability conditions of the equilibrium point E2. When B11>0, then θi>0,i=1,2,3,4. Also, θ1θ2θ3−θ32−θ12θ4=(B11(a+x2)3Ω4+Ω2Ω3)(B11δμ2(a+x2)3+Ω1(α2+μ2)+eaβ2w2x2(δ+μ3))2(a+x2)9, where
when Ω2Ω3>0, then θ1θ2θ3−θ32−θ12θ4>0, therefore all the eigenvalues of the Jacobian matrix J(E2) near the equilibrium point E2 have negative real parts. Thus, due to the Routh-Hurwitz criterion the equilibrium point E2 is locally asymptotically stable. The local stability of the coexistence equilibrium point E2 is given in the following theorem.
Theorem 5.IfΩ2Ω3>0andB11>0, then the coexistence equilibrium pointE2is locally asymptotically stable.
The global stability of the coexistence equilibrium point is investigated in the following theorem.
Theorem 6.Ifβz2a(a+x2)<rk, then the coexistence equilibrium pointE2is globally asymptotically stable.
Proof. The following positive definite Lyapunov function is considered.
L2=H(x−x2−x2lnxx2)+z−z2−z2lnzz2.
By calculating the time derivative of L2 along the solution of system (1), one obtains,
Choosing H=eaa+x2, then cαDαL2<0. According to generalized Lyapunov–Lasalle’s invariance principle [20], the coexistence equilibrium point E2 is globally asymptotically stable when βz2a(a+x2)<rk.
Bifurcation Analysis
In this section the local bifurcations near the equilibrium points of LDM model (1) are inves-tigated with the help of Sotomayor’s theorem [22]. The Hopf bifurcation theorem given in Liu [23] is also presented to discuss the bifurcation analysis of the underlying system. One can compute
where ν1 is any non zero real number. Similarly, suppose V2=(τ1,τ2,τ3,τ4,τ5)T be the eigenvector corresponding to the zero eigenvalue of the matrices J(E0), thus J(E0)TV2=0 gives V2=(1,0,0,0,0)T. Consider, ∂F∂r=Fr(X,r)=(x(1−xk),0,0,0,0)T, thus, V2TFr(E0,r∗)=0. Therefore, according to Sotomayor’s theorem for local bifurcation, the LDM model (1) has no saddle-node bifurcation near E0 at r∗=α1+μ1. One can note that r=α1+μ1 is equivalent to R0=1. Now,
DFr(E0,r∗)=(1000000000000000000000000),
then V2TDFr(E0,r∗)V1=ν1≠0. Using (5), one obtains
V2TD2F(X,r)(V1,V1)=−2rkν12−2βaν1ν5≠0.
Thus, according to Sotomayor’s theorem, the LDM system (1) has a transcritical bifurcation at r∗=α1+μ1 as the parameter r passes through the value r∗, thus the proof is complete.
Theorem 8.The LDM system (1) undergoes a transcritical bifurcation with respect to the bifurcation parameterμ5aroundE1=(x1,y1,f1,m1,0)ifμ5=μ5∗=eβx1a+x1.
Proof. The Jacobian matrix of the LDM system (1) at the free enemy equilibrium point E1 with μ5=μ5∗ has zero eigenvalue takes the form
where A11=r−2rx1k−(α1+μ1),A22=(α2+μ2),A33=(δ+μ3) and A44=γpy1+p+μ4.J(E1)V3=0, gives the eigenvector corresponding to the zero eigenvalue of the matrix J(E1), hence
where ν5 is any non zero real number. Similarly, J(E1)TV4=0 gives the eigenvector corresponding to the zero eigenvalue of the matrix J(E1)TV4=0, hence V4=(0,0,0,0,1)T. Consider, ∂F∂μ5=Fμ5(X,μ5)=(0,0,0,0,−z)T, thus, V4TFμ5(E1,μ5∗)=0. Therefore, according to Sotomayor’s theorem for local bifurcation, the LDM model (1) has no saddle-node bifurcation near E1 at μ5=μ5∗. Now,
DFμ5(E1,μ5∗)=(000000000000000000000000−1),
then, V4TDFμ5(E1,μ5∗)V3=−ν5≠0. Using (5), one obtains
V4TD2F(X,μ5)(V3,V3)=eaβν1ν5(a+x1)2≠0.
Thus, according to Sotomayor’s theorem, the LDM system (1) has a transcritical bifurcation at μ5∗=eβx1a+x1 as the parameter μ5 passes through the value μ5∗, thus the proof is complete.
In this part, we shall show that as the coexistence equilibrium loses stability, periodic solutions can bifurcate from the positive equilibrium. We first give the following lemma.
Lemma 9 .The characteristic Eq.(4) has a pair of purely imaginary roots and the remaining roots have negative real parts if and only ifz2=r(a+x2)2kβandθ1θ2θ3−θ32−θ12θ4=0.
Suppose (4) has two eigenvalues which have negative real parts and two complex conjugates eigenvalues (call them λ=m(φ)±in(φ)) such that m(φ∗)=0,n(φ∗)>0,dmdφ|φ=φ∗≠0. Substituting λ=m(φ)±in(φ) into (4), and separating the real and imaginary, we get
m4+θ1m3+θ2m3+θ3m+θ4−(6m2+3θ1m+θ2)n2+n4=0,
4m3+3θ1m2+2θ2m+θ3−(4m+θ1)n2=0,
Substituting (6) into (7), differentiating with respect to φ and utilizing m(φ∗)=0 and n(φ∗)≠0, we have
Theorem 10.For the coexistence equilibrium point E2of the integer order LDM system (1), the system aroundE2enters into the Hopf bifurcation whenφpassesφ∗if the coefficientsθj(φ)(j=1,2,3,4)atφ=φ∗satisfying the following condition:
where Wi(i=1,2,3,4,5) are independent standard Brownian motions with Wi(0)=0 and σi>0 denote the intensities of the white noise. The stochastic fractional-order LDM model (8) can be written in the general form:
CαDαX(t)=F(X)+g(X)dWdt,
where F(x) is given in (2), g(x)=(σ1x,σ2y,σ3f,σ4m,σ5z) and dWdt=(dW1dt,dW2dt,dW3dt,dW4dt,dW5dt)T. Applying Riemann–Liouville integral to both sides of (9), one can obtain the following stochastic Volterra integral equation.
According to Wang et al. [24,25], under some conditions on the coefficient functions, the global existence and uniqueness of solutions for the stochastic fractional-order system (8) can be investigated. Because Grunwald-Letnikov’s definition is the most straightforward from the point of view of numerical implementation, so we will use it to solve the LDM system of fractional order stochastic differential equations. Grunwald-Letnikov (GLDα) fractional derivative of order α defined by Aminikhah et al. [26,27]
GLαDαf(t)=Limh→0h−α∑j=0[t−ah](−1)j(αj)f(t−jh),
where [t−ah] means the integer part of t−ah. This formula (11) can be reduced to
GLαDαf(tn)≈h−α∑j=0nwjαf(tn−j),
where h is the time step, tn=nh and wjα are the Grunwald-Letnikov coefficients satisfy the following recurrence relationship
w0α=1,wjα=(1−1+αj)wj−1α,j=1,2,3,…
If f(t) is continuous function and f′(t) is integrable function in the interval [0,T], then the relation between Caputo and Grunwald-Letnikov fractional derivative takes the form [28,29]
To show the effect of the intrinsic growth rate of the LDM larvae, we draw the bifurcation diagram concerning r as a bifurcation parameter. It can be seen that a transcritical bifurcation occurs at r = 0.51 as shown in Fig. 1 and stated in Theorem 7. It can also be observed that when r > 0.51, the coexistence equilibrium point E2 = (1.6667, 25.2475, 24.7525, 25.36, 157.407r − 85) is locally asymptotically stable. According to Theorem 10, it can be seen that the supercritical Hopf bifurcation value localized at r = 1.31143. For r > 1.31143 the LDM system (1) undergoes limit cycle behaviour.
To show the effect of the natural enemy’s mortality rate around the coexistence equilibrium points, we draw the bifurcation diagram for µ5 as a bifurcation parameter. It can be seen that the Hopf bifurcation value localized at µ5 = 0.0284692 as shown in Fig. 2. It can also be observed that when µ5> 0.0284692, the coexistence equilibrium point E2 is locally asymptotically stable. For µ5< 0.3235, the system undergoes limit cycle behaviour. According to Theorem 7, for µ5> 0.0642346 the natural enemy goes extinct from the system and a transcritical bifurcation occurs at µ5 = 0.0642346 and the equilibrium point E1 = (26.94, 408.101, 400.009, 528.446, 0) is locally asymptotically stable.
In order to show the effect of the predation rate, we draw the bifurcation diagram with respect to β as a bifurcation parameter. It can be seen that a transcritical bifurcation occurs at β = 0.0501342 as shown in Fig. 3. It can also be observed that when β < 0.0501342, the free enemy equilibrium point E1 = (22.35, 338.569, 331.93, 436.688, 0) is locally asymptotically stable as indicated in Theorem 8. For 0.0501342 < β < 0.152449, the coexistence equilibrium point E2 is locally asymptotically stable. It can be seen that the supercritical Hopf bifurcation value localized at β = 0.152449. For β > 0.152449 the LDM system (1) undergoes limit cycle behaviour.
From Fig. 4a, it can be seen that the sex pheromone trap parameter p is important in that it affects the population density of the LDM male. One can observe from Fig. 4b. that the population density of LDM male decrease with increasing p. We conclude that the dynamics of LDM can be controlled by sex pheromone trap parameters p.
Fig. 5 shows that the deterministic fractional-order remains stable for different values of fractional-order α though solutions reach to equilibrium point E2(7.5, 113.614, 111.386, 140.149, 189) more slowly for a smaller value of fractional-order α. It is important to notice that when α = 1 the fractional order model for lesser date moth (1) reduces to the classical integer-order model [10–14].
When the strength of environmental noise is very close to zero, the fractional-order stochastic LDM system (8) behaves like the deterministic model. Fig. 6a indicates the dynamical behavior of stochastic LDM (8) without noise (i.e., σi = 0), which gives the deterministic model (1). In Fig. 6a, the dynamical behavior of the LDM system was unstable, which corresponds with the results of Theorem 5. Fig. 6b shows that the LDM smooth oscillations’ dynamical behavior when the strength of the noise was low (σi = 0.02). However, with an increase in the strength of noise, such as the medium-noise situation shown in Fig. 6c, the dynamical behavior of the LDM became more complex, and they tended to extinction. Fig. 6c represents the dynamical behavior of the model (1) when the noise strength was high (σi = 0.05). The natural enemy can die out due to the white noise stochastic disturbance. By comparing Figs. 6a and 6b, one can realize that if the noise is not strong, the stochastic perturbation does not cause sharp changes in the LDM model (8). However, when the environmental noise σi is sufficiently large (see Fig. 6c), the noise can force the natural enemy to become extinct.
Bifurcation diagram of LDM model (1) with respect to r
Bifurcation diagram of LDM model (1) with respect to µ5
Bifurcation diagram of LDM model (1) with respect to β
(a) Plot of the LDM male versus time with p = 0.1, 0.3, 0.8, (b) Plot of the LDM male versus time with 0 ≤ p ≤ 1
Time series of the LDM female with α = 0.94, 0.97, 1
Time series of the fractional order stochastic model (10) with σ = 0, 0.02, 0.05 (a) σi = 0 (b) σi = 0.02 (c) σi = 0.05
Discussion and Conclusion
In this paper, we consider a deterministic and stochastic fractional order model for LDM using mating disruption and natural enemies. We obtain some sufficient conditions that ensure the local and global stability of equilibrium points. We conclude that sex pheromone trap parameters can control dynamics of LDM. The occurrence of local bifurcation near the equilibrium point is performed using Sotomayor’s theorem. Numerical simulations are performed to support and illustrate the theoretical findings. From the numerical results, one can realize that if the noise is not strong, the stochastic perturbation does not cause sharp changes in the dynamics of the LDM model. However, when the environmental noise is sufficiently large, the noise can force the population to become extinct.
Funding Statement: The authors gratefully acknowledge Qassim University, represented by the Deanship of Scientific Research, on the financial support under the number (cosao-bs-2019-2-2-I-5469) during the academic year 1440 AH/2019 AD.
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.
ReferencesL. I.El-Juhany, “
Degradation of date palm trees and date production in Arab Countries: Causes and potential rehabilitation,”
, vol.
4, no.
8, pp.
3998–
4010,
2010.J.Al-Khayri, “
Date palm phoenix dactylifera micropropagation,” in Protocols for Micropropagation of Woody Trees and Fruits,
Springer, pp.
509–
526,
2007. A.Levi-Zada, D.Fefer, L.Anshelevitch, A.Litovsky, M.Bengtssonet al., “
Identification of the sex pheromone of the lesser date moth, Batrachedra amydraula, using sequential SPME auto-sampling,”
, vol.
52, no.
35, pp.
4550–
4553,
2011.M. A.Al-Deeb and H. A.Al-Dhaheri, “
Use of a pheromone-baited trap to monitor the population of the lesser date moth Batrachedra amydraula (Lepidoptera: Batrachedridae) in the UAE,”
, vol.
5, no.
6, pp.
2572–
2575,
2017.M.Hoddle, A. H.Al-Abbad, H.El-Shafie, J.Faleiro, A.Sallamet al., “
Assessing the impact of areawide pheromone trapping, pesticide applications, and eradication of infested date palms for Rhynchophorus ferrugineus (Coleoptera: Curculionidae) management in Al Ghowaybah, Saudi Arabia,”
, vol.
53, pp.
152–
160,
2013.E.Sadeghi, V.Baniameri and A.Marouf, “
Oviposition behaviour of Goniozus swirskiana (Hymenoptera: Bethylidae: Bethylinae) a parasitoid of Batrachedra amydraula Meyrick from the warmest desert of Iran,”
, vol.
20, no.
11, pp.
1493–
1498,
2012.M.Abbas, S.Al-Khatry, R.Shidi and N. A.Al-Ajmi, “
Natural enemies of the lesser date moth, Batrachedra amydraula Meyrick (Lepidoptera: Batrachedridae) with special reference to its parasitoid Goniozus sp,”
, vol.
24, no.
2, pp.
293–
296,
2014.L.Addison, “
Analysis of a predator-prey model: A deterministic and stochastic approach,”
, vol.
8, no.
4, pp.
359–
368,
2017.B.Ahmad, S. K.Ntouyas, A.Alsaedi and A. F.Albideewi, “
A study of a coupled system of Hadamard fractional differential equations with nonlocal coupled initial-multipoint conditions,”
, vol.
2021, no.
1, pp.
1–
16,
2021.R.Anguelov, C.Dufourd and Y.Dumont, “
Mathematical model for pest-insect control using mating disruption and trapping,”
, vol.
52, pp.
437–
457,
2017.J. P.Ntahomvukiye, A.Temgoua and S.Bowong, “
Study of the population dynamics of Busseola fusca, maize pest,”
, vol.
66, no.
4, pp.
379–
397,
2018.S.Xiang, Y.Pei and X.Liang, “
Analysis and optimization based on a sex pheromone and pesticide pest model with gestation delay,”
, vol.
12, no.
5, pp.
1950054,
2019.M. D.Tapi, L.Bagny-Beilhe and Y.Dumont, “
Miridae control using sex-pheromone traps. Modeling, analysis and simulations,”
, vol.
54, pp.
103082,
2020.M.El-Shahed and A.Al-Dubiban, “
Mathematical modelling of lesser date moth using sex pheromone traps and natural enemies,”
, vol.
2021, pp.
1–
14,
2021.J.Cresson and A.Szafran´ska, “
Discrete and continuous fractional persistence problems the positivity property and applications,”
, vol.
44, pp.
424–
448,
2017.H.Li, L.Zhang, C.Hu, Y.Jiang and Z.Teng, “
Dynamical analysis of a fractional-order predator- prey model incorporating a prey refuge,”
, vol.
54, pp.
435–
449,
2016.M.Sambath, P.Ramesh and K.Balachandran, “
Asymptotic behavior of the fractional order three species prey-predator model,”
, vol.
19, pp.
721–
733,
2018.S. K.Choi, B.Kang and N.Koo, “
Stability for Caputo fractional differential systems,”
, vol.
2014, pp.
1–
6,
2014.P. V.Driessche and J.Watmough, “
Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission,”
, vol.
180, pp.
29–
48,
2002.J.Huo, H.Zhao and L.Zhu, “
The effect of vaccines on backward bifurcation in a fractional order HIV model,”
, vol.
26, pp.
289–
305,
2015.A.Matouk, “
Stability conditions, hyperchaos and control in a novel fractional order hyperchaotic system,”
, vol.
373, no.
25, pp.
2166–
2173,
2009.L.Perko,
. Vol.
7,
Springer Science & Business Media,
2013.W. M.Liu, “
Criterion of Hopf bifurcations without using eigenvalues,”
, vol.
182, no.
1, pp.
250–
256,
1994.W.Wang, S.Cheng, Z.Guo and X.Yan, “
A note on the continuity for Caputo fractional stochastic differential equations,”
, vol.
30, no.
7, pp.
073106,
2020.T.Doan, P.Huong, P.Kloeden and A.Vu, “
Euler-Maruyama scheme for Caputo stochastic fractional differential equations,”
, vol.
380, pp.
112989,
2020.H.Aminikhah, A. H. R.Sheikhani, T.Houlari and H.Rezazadeh, “
Numerical solution of the distributed-order fractional Bagley-Torvik equation,”
, vol.
6, no.
3, pp.
760–
765,
2017.M.Hou, J.Yang, S.Shi and H.Liu, “
Logical stochastic resonance in a nonlinear fractional-order system,”
, vol.
135, no.
9, pp.
747,
2020.I.Podlubny,
.
New York, NY, USA:
Academic Press,
1999.C.Huang and M.Stynes, “
Error analysis of a finite element method with GMMP temporal discretisation for a time-fractional diffusion equation,”
, vol.
79, no.
9, pp.
2784–
2794,
2020.