# A Fractional-Order Eco-Epidemiological Leslie–Gower Model with Double Allee Effect and Disease in Predator

Authors: Emli Rahmi; Isnani Darti; Agus Suryanto; Trisilowati Trisilowati
Journal: International Journal of Differential Equations
 (2023)

---

## Abstract

In this paper, a fractional order of a modified Leslie–Gower predator-prey model with disease and the double Allee effect in predator population is proposed. Then, we analyze the important mathematical features of the proposed model such as the existence and uniqueness as well as the non-negativity and boundedness of solutions to the fractional-order system. Moreover, the local and global asymptotic stability conditions of all possible equilibrium points are investigated using Matignon’s condition and by constructing a suitable Lyapunov function, respectively. Finally, numerical simulations are presented to verify the theoretical results. We show numerically the occurrence of two limit cycles simultaneously driven by the order of the derivative, the bistability phenomenon for both the weak and strong Allee effect cases, and more dynamic behaviors such as the forward, backward, and saddle-node bifurcations which are driven by the transmission rate. We have found that the risk of extinction for the predator with a strong Allee effect is much higher when the spread of disease is relatively high.


---

## Body

## 1. Introduction

The study of dynamical behavior of the predator-prey model has been an important theme in mathematical biology; see, for instance [1–5], and references therein. The results of these studies could give us substantial elucidation about the qualitative behavior of the density of each population in the future time, such as stability, bifurcations, and chaos, without having any experimental laboratory. In the last few decades, special attention has been paid to the merger of the predator-prey model with the model for transmissible disease as a new branch of mathematical biology known as the eco-epidemiological model. It becomes an important tool in analyzing the effect of infectious diseases on the existence (extinction) of prey or/and predator populations. Several researchers have discussed the eco-epidemiological model incorporating ecological factors, such as harvesting [6], prey refuge [7, 8], social behavior [9], and the Alee effect [10].One interesting phenomenon was introduced by [11] which demonstrated a condition where, at low densities of population, the presence of conspecifics could increase the per capita growth rate of the population. This phenomenon is called the Allee effect. Some population models involving the Allee effect have been developed to exhibit the important impact of this phenomenon on many aspects of population ecology, such as conservation of threatened species [12, 13], controlling pest species [14], and harvesting management [15, 16]. In general, there are two types of the Allee effects: the strong Allee effect and the weak Allee effect. The strong Allee effect has a threshold size population which is called the Allee threshold. When the density of population is under this threshold then the per capita growth rate of the population becomes negative. On the other hand, the weak Allee effect gives a reduction in the per capita growth rate when the population is low. Furthermore, several researchers investigate the dynamical behavior of a system incorporating any situation in which two or more components of the Allee effect can work simultaneously on a single population which is called the double Allee effect [1, 17]. These components could arrive from the reproduction or/and survival mechanisms such as the difficulty of finding a mate, cooperative breeding dependencies, cooperative antipredator behavior, and environmental conditioning, as shown in Table 1 [18].Table 1
Parameter description in system (2).
ParametersEcological descriptionsr^1Intrinsic growth rate of preyr^2Intrinsic growth rate of susceptible predatorK^Carrying capacity of preyθ^Capture rate of susceptible predatorb^Prey saturation constanta^Another saturation constantc^The strength of interference between susceptible predators in predation processp^The Allee thresholdq^The auxiliary Allee effect constantm^Intraspecific competition rate among predatorsk^Alternative food for the predatorβ^Transmission rate of the diseaseδ^Death rate of infected predatorNowadays, fractional-order derivative modelling has become popular in many fields of science, such as physics, thermodynamics, biology, control theory, and many others (see, for example, [19–22]). Modeling using a fractional-order derivative has a major advantage in that it involves memory which comes from the fact that the fractional-order operators are nonlocal. More precisely, all of the previous conditions are captured in defining the fractional-order derivative. It has some advantages over the integer order derivative, especially the ability to describe the memory and hereditary properties which are inherent in various processes [23–25]. Some scholars had proven that systems with memory are more consistent and adequate with the real phenomena [26, 27].Recently, [28] studied a fractional-order predator-prey model by assuming that the predation follows the Beddington-DeAngelis functional response and the double Allee effect on the predator population. By the last assumption, there is always a bistability condition for the strong Allee effect case. In other words, there is an extinction or existing condition for the predator depending on the choice of the initial condition. In this paper, we consider a predator-prey model in [28] by adding the assumption that there is a transmission of disease in the predator population. Biologically, the growth rates of a population must depend on the history of its previous conditions not only on the local conditions, and for this reason, the proposed model will use the fractional-order derivative which has a memory effect, to make it more accurate in predicting the future condition of the population. To the best of our knowledge, no investigation has been carried out on the dynamics of the proposed system. An important objective of this paper is to study how the spread of disease may affect changes in predator densities over time, where the intrinsic growth rate of the predator is affected by the double Allee effect.We have arranged this paper in the following manner. First, we propose a model and give some useful preliminaries which are used in our analysis briefly. Next, we prove the basic properties of solutions. We also study the existence of the equilibrium points, their local stability, and the Hopf bifurcation. The global stability conditions for the equilibrium points are also examined. Moreover, some numerical simulations with a set of hypothetical parameters are demonstrated to validate the theoretical results and further explore the role of the disease in the predator-prey interaction. Finally, conclusions are given in the last section.


## 2. Formulation Model and Useful Preliminaries

In [28], the author’s considered a predator-prey model with memory effect, which is a modification of the Leslie–Gower model using the Caputo fractional derivative. The model abandoned the mass conservation principle held by the Lotka–Volterra model because of the assumption that both prey and predator obey the logistic law. They also assume that the predation process is governed by the Beddington-DeAngelis type functional response [29] which is a nonlinear prey-predator-dependent functional response, and that the intrinsic growth rate of the predator population is influenced by the double Allee effect, as given in the following form:(1)DtαCN=r^1N1−NK^−θ^NP1+c^N+q^P,DtαCP=Ps^P−p^P+q^−m^Pk^+N.In this paper, we assume that the predator population is split into two subpopulations, namely the susceptible predatorS and infected predator I to include the disease transmission among predator. It is assumed that only susceptible predators can reproduce according to the logistic law, where the intrinsic growth rate of the susceptible predator is induced by the double Allee effect. The prey is only captured by the susceptible predator. This is reasonable from an ecological view as the infected predators are physically weak and lose their skill to catch food [30, 31]. It is also assumed that the susceptible predator has intraspecific competitive interactions with both the susceptible and infected predators. The transmission of the nonrecoverable disease in the predator population in the manner of the simple law of mass action and is vertically transferred from the susceptible predator to the infected predator compartment with a constant transmission rate β. Here, we also assume that the disease does not spread from the predator population to the prey population.Considering the above assumptions, we construct a fractional-order model with the double Allee effect and disease in predator as follows:(2)DtαCN=r^1N1−NK^−θ^NSa^+b^N+c^S,DtαCDS=Sr^2S−p^S+q^−m^S+I^k^+N−β^SI,DtαCI=β^SI−δ^I,where DtαC represents the Caputo fractional derivative with 0<α<1. The variables N,S, and I denote the density of the prey, susceptible predator, and infected predator populations, respectively, and all parameters are shown in Table 1.Since the right hand side of system (2) have time dimension time−1 which is inconsistent with time dimension on the left hand side time−α, we modify the above system (2) as follows:(3)DtαCN=r^1αN1−NK^−θ^αNSa^+b^N+c^S,DtαCS=Sr^2αS−p^S+q^−m^αS+Ik^+N−β^αSI,DtαCI=β^αSI−δ^αI.For the sake of convenience, we redefine parameters of the system (3) as r1=r^1α,K=K^,θ=θ^α,a=a^,b=b^,c=c^,r2=r^2α,p=p^,q=q^,m=m^α,k=k^,β=β^α, and δ=δ^α. Therefore, the system (3) with the initial conditions can be written in the following form:(4)DtαCN=r1N1−NK−θNSa+bN+cS,DtαCS=Sr2S−pS+q−mS+Ik+N−βSI,DtαCI=βSI−δI,where N0>0,S0>0, and I0>0. Next, we briefly introduce definition and some basic properties of the Caputo fractional derivative.Definition 1.
(see [32]). Let α∈0,1. The Caputo fractional derivative for the real valued function ht is defined as follows:(5)DtαCht=1Γ1−α∫0th′st−sαds,where Γ⋅ is the Gamma function, t≥0.
To verify the non-negativity and uniform boundedness of the solutions of the system (4), we require the following lemma for the fractional derivative.Lemma 2 (see [33]).
Suppose thatht is continuous in a,b and its derivatives DtαCht is continuous in a,b, for 0<α<1, then we have(6)ht=ha+1ΓαDtαChτt−aα,a≤τ≤t,∀t∈a,b.Lemma 3 (see [33]).
Suppose thatht is continuous in a,b and DtαCht is continuous in a,b for 0<α<1. If DtαCht≥0,∀t∈a,b, then ht is a nondecreasing function for all t∈a,b. If DtαCht≤0,∀t∈a,b, then ht is a nonincreasing function for all t∈a,b.Lemma 4 (comparison lemma [34]).
Supposeht is continuous in 0,∞. If ht satisfies(7)DtαCht≤−φht+μ,h0=h0∈ℝ,where 0<α<1,φ,μ∈ℝ2,andφ≠0, then(8)ht≤h0−μφEα−φtα+μφ.Next, we use the following theorem to determine the local stability behavior of the equilibrium points of the system (4).Theorem 5 (Matignon condition [35]).
Consider the following Caputo fractional-order system with initial value(9)DtαCh⟶t=f⟶t,h⟶,h⟶0=h⟶0,where h⟶∈ℝn and 0<α<1. If h⟶∗satisfies f⟶t,h⟶∗=0⟶, then h⟶∗ is an equilibrium point. If all eigenvalues λi, i=1,…,n of the Jacobian matrix Jh⟶∗ satisfy argλi>απ/2 then the equilibrium points h⟶∗ is locally asymptotically stable.To establish the conditions for the global stability of the equilibrium points of the system (4), we provide the following lemma.Lemma 6 (Volterra-type Lyapunov function [36]).
Letht∈ℝ+ be a continuous and derivable function. Then, for t≥0,(10)DtαCht−h∗−h∗lnhth∗≤1−h∗htDtαCht,h∗∈ℝ+,∀α∈0,1.Lemma 7 (generalized LaSalle invariance principle [37]).
SupposeΩ is a bounded closed set. Every solution of DtαCht=fht which starts from a point in Ω remains in Ω for all t≥0. If there exists Vh:Ω⟶ℝ with continuous first partial derivatives satisfying(11)DtαCV|DtαCht=fht≤0,then every solution ht starting in Ω tends to M as t⟶∞ where M is the largest invariant set of E with(12)E≔hDtαCVDtαCht=fht=0.


## 3. Basic Properties of Solutions

In this section, we present some basic properties of solutions of the system (4) such as the non-negativity, boundedness, existence, and uniqueness of solutions.


### 3.1. Non-Negativity and Boundedness of Solutions

From the biological point of view, we only concern to system (4) in the non-negative and bounded solutions. We consider ℝ+3≔N,S,I∈ℝ3:N≥0,S≥0,I≥0 as the set of all non-negative solutions. The following theorem shows the non-negativity and boundedness of the solutions of system (4).Theorem 8.
All solutions of the system (4) which initiate in ℝ+3 are non-negative and uniformly bounded.Proof.
LetN0>0,S0>0, and I0>0 and assume that Nt≥0 for t∈0,∞. Suppose the assumption is not true, then there exists t∗ which satisfies 0≤t<t∗ such that Nt>0 for t∈0,t∗, Nt∗=0, and Nt<0 for t>t∗. From the first equation of system (4), we have(13)DtαCNt∗|Nt∗=0=0.
According to Lemma3, we get Nt=0 for t>t∗ which contradicts to Nt<0,∀t>t∗. Therefore, we have Nt≥0 for t∈0,∞. Using the similar way, it can be shown that St≥0 and It≥0 for t∈0,∞. Hence, all solutions of the system (4) are non-negative.
Next, we show the boundedness of solutions of the system (4) by using the fractional comparison lemma. From the first equation of system (4), we obtain(14)DtαCNt+Nt=r1N1−NK−θNSa+bN+cS+N≤−r1N2K+r1+1N=−r1KN−Kr1+12r12+Kr1+124r1≤Kr1+124r1.
Based on Lemma4, we have(15)Nt≤N0−Kr1+124r1Eα−tα+Kr1+124r1,where Eα is the Mittag–Leffler function. Lemma 5 and Corollary 6 in [38] implies that Eα−tα⟶0 as t⟶∞, then it follows(16)0≤Nt≤Kr1+124r1=η1,ast⟶∞.
Thus, for anyε>0, Nt with initial condition in ℝ+3 is uniformly bounded to the region(17)Θ1=Nt≤η1+ε,ε>0.
DefineWt=St+It and adding the last two equations of system (4). Then, we get(18)DtαCWt≤r2S−pr2SS+q−mS2k+N−δI.
Forξ>0 and Nt≤η1+ε, we have(19)DtαCWt+ξWt≤−mk+η1+εS−r2+ξk+η1+ε2m2+r2+ξ2k+η1+ε4m+ξ−δI−pr2SS+q.
If we chooseξ<δ, then(20)DtαCWt+ξWt≤−mk+η1+εS−r2+ξk+η1+ε2m2+r2+ξ2k+η1+ε4m−pr2SS+q.≤r2+ξ2k+η1+ε4m−pr2SS+q.
Using the same argument as before and lettingε⟶0, then for the weak Alee effect p<0 we have the following equation:(21)0≤Wt≤r2+ξ2k+η14ξm−pr2ξq=η2,ast⟶∞.
On the other hand, for the strong Alee effectp>0 we have the following equation:(22)0≤Wt≤r2+ξ2k+η14ξm=η3,ast⟶∞.
Thus, for anyε>0, the combination of the solutions St and It with initial conditions started in ℝ+3 are uniformly bounded in the region(23)Θ2=St+It≤η2+ε,ε>0,p<0,St+It≤η3+ε,ε>0,p>0.


### 3.2. Existence and Uniqueness of the Solution

To investigate the existence and uniqueness of solution of system (4), we apply the locally Lipschitz condition in the region 0,T×ϒ, where ϒ=N,S,I∈ℝ3:maxN,S,I≤Φ, for sufficient large Φ and T<∞.Theorem 9.
For eachZ0=N0,S0,I0∈ϒ, there exists a unique solution Zt=Nt,St,It∈ϒ of system (4) with initial condition Z0, which is defined for all t≥0.Proof.
For anyZ,Z¯∈ϒ, we consider a mapping of a continuous function FZ=F1Z,F2Z,F3Z where(24)F1Z=r1N−r1N2K−θNSa+bN+cS,F2Z=r2S2S+q−pr2SS+q−mS2k+N−mSIk+N−βSI,F3Z=βSI−δI.
It follows that(25)FZ−FZ¯=F1Z−F1Z¯+F2Z−F2Z¯+F3Z−F3Z¯=|r1N−r1N2K−θNSa+bN+cS−r1N¯−r1N¯2K−θN¯S¯a+bN¯+cS¯|+|r2S2S+q−pr2SS+q−mS2k+N−mSIk+N−βSI−r2S¯2S¯+q−pr2S¯S¯+q−mS¯2k+N¯−mS¯I¯k+N¯−βS¯I¯|+βSI−δI−βS¯I¯−δI¯≤r1N−N¯+r1KN2−N¯2+θNSa+bN+cS−N¯S¯a+bN¯+cS¯+r2S2S+q−S¯2S¯+q+pr2SS+q−S¯S¯+q+mS2k+N−S¯2k+N¯+mSIk+N−S¯I¯k+N¯+βSI−S¯I¯+βSI−S¯I¯+δI−I¯≤r1N−N¯+r1KN+N¯N−N¯+θaN+bN¯Na2S−S¯+θaS¯+cS¯Sa2N−N¯+r2q2S¯SS−S¯+r2qS+S¯S−S¯+pr2qS−S¯+mk+N¯S+S¯k2S−S¯+mS¯2k2N−N¯+mk+N¯Ik2S−S¯+mN¯Sk2I−I¯+mS¯Ik2N−N¯+2βSI−I¯+2βIS−S¯+δI−I¯≤M1N−N¯+M2S−S¯+M3I−I¯≤MZ−Z¯,where(26)M1=r1+2r1K+θaΦ+θca2+2mk2Φ2,M2=pr2q+θa+2r2q+3mk+2βΦ+θba2+r2q2+3mk2Φ2,M3=δ+2βΦ+mk2Φ2,M=maxM1,M2,M3.
SinceFZ satisfies the Lipschitz condition with respect to Z, then system (4) with initial condition Z0 has a unique solution Zt∈ϒ for all t≥0. The proof is completed.


### 3.3. Existence and Stability of Equilibrium Points

To study the dynamical behaviors of system (4), we first investigate the equilibrium points of the system (4) which are the constant solutions to the following system:(27)r1N1−NK−θNSa+bN+cS=0,Sr2S−pS+q−mS+Ik+N−βSI=0,βSI−δI=0.Next, the local stability of equilibrium points of system (4) is evaluated by computing the eigenvalues of the Jacobian matrix of system (4) at the equilibrium point N∗,S∗,I∗ which is given as follows:(28)JN∗,S∗,I∗=J1,1J1,20J2,1J2,2J2,30J3,2J3,3,where(29)J1,1=r11−N∗K−r1N∗K−θS∗a+bN∗+cS∗+bθN∗S∗a+bN∗+cS∗2,J1,2=−θN∗a+bN∗+cS∗+cθN∗S∗a+bN∗+cS∗2J2,1=mS∗S∗+I∗k+N∗2J2,2=r2S∗−pS∗+q−mS∗+I∗k+N∗+S∗r2S∗+q−r2S∗−pS∗+q2−mk+N∗−βI∗J2,3=−mS∗k+N∗+βS∗J3,2=βI∗J3,3=βS∗−δ.


### 3.4. Equilibrium Analysis in the Case of the Weak Allee Effect in Predator

In this subsection, the existence and stability of the equilibrium points of the system (4) for the weak Alee effect p<0 are given. By solving system (9) simultaneously, all positive equilibrium points of system (4) are obtained as follows:(a)
Trivial equilibrium point:E0=0,0,0, which always exists.(b)
Axial equilibrium point:EN=K,0,0, which always exists; and ES=0,Sw,0, where Sw=kr2−mq/2m+kr2−mq/2m2−kpr2/m, which always exists.(c)
Planar equilibrium point:ESI=0,S¯w,I¯w where S¯w=δ/β and I¯w=1/m+βkβkr2δ−βp−δmδ+βq/βδ+βq exists if p<min0,δ/β1−mδ+βq/βkr2; and ENS=N˜w,S˜w,0 exists if p>S˜w1−mS˜w+q/kr2, where N˜w=mS˜wS˜w+q/r2S˜w−p−k and S˜w are all positive roots of the following quartic equation:(30)l1S4+l2S3+l3S2+l4S+l5=0,where(31)l1=mb+cr2mr1,l2=−bk−cp+ar2+mbqmr1−mb+cr2k+Kr2−mqr1+θKr22,l3=−p−bk+amr1r2−−bk−cp+ar2+mbqk+Kr2−mqr1+mb+cr2pk+Kr2r1−2θKr22p,l4=pr2k+Kr2−mq−bk+ar1+−bk−cp+ar2+mbqpk+Kr2r1+θKr22p2,l5=−p−bk+apk+Kr2r2r1.(d)
Interior equilibrium point:E∗=Nw∗,Sw∗,Iw∗, where(32)Nw∗=bK−a+cSw∗2b±Δw,Sw∗=δβ,Iw∗=r2Sw∗−pk+Nw∗−mSw∗Sw∗+qSw∗+qm+βk+Nw∗,Δw=bK−a+cSw∗2b2+Kr1a+cSw∗−θSw∗br1.Assume thatk+Nw∗>δmδ+βq/βr2δ−βp, then we have the following results:(a)
IfΔw<0, then system (4) has no interior equilibrium point.(b)
IfΔw=0 and aβ+cδ<bKβ, then system (4) has a unique interior equilibrium point E∗=Nw∗,δ/β,Iw∗ where Nw∗=bKβ−aβ+cδ/2bβ and Iw∗=βr2δ−βpk+Nw∗−δmδ+βq/βδ+βqm+βk+Nw∗.(c)
IfΔw>0 and aβ+cδ<minbKβ,δθ/r1, then system (4) has two positive interior equilibrium points E1,2∗=Nwi∗,δ/β,Iwi∗, i=1,2, where(33)Nw1,2∗=bKβ−aβ+cδ2bβ±Δw,Iw1,2∗=βr2δ−βpk+N^w1,2−δmδ+βqβδ+βqm+βk+N^w1,2,Δw=bKβ−aβ+cδ2bβ2+Kr1aβ+cδ−δθbr1β.Theorem 10.
The local stability of trivial, axial, and planar equilibrium points of system (4) for the weak Alee effect p<0 are summarized as follows.(a)
E0=0,0,0 is always a saddle point.(b)
EN=K,0,0 is always a saddle point.(c)
ES=0,Sw,0 is locally asymptotically stable if r1<θSw/a+cSw, p+q<m/kr2Sw+q2, and Sw<δ/β.(d)
LetΔ^w=βδr2p+q/δ+βq2−mδ/βk2−4δ/βkβkr2δ−βp−mδδ+βq/δ+βq. ESI=0,S¯w,I¯w is locally asymptotically stable if r1<δθ/aβ+cδ and one of the following conditions holds:(i)
βr2p+q/δ+βq2<m/βk<r2δ−βp/δδ+βq, or(ii)
βr2p+q/δ+βq2>m/βk, Δ^w<0, and α<α∗=2/πtan−1Δ^w/βr2p+q/δ+βq2−m/βk.(e)
Suppose that(34)ω1=−r1N˜wK1−br1K−N˜wθKS˜w+r2S˜w−p2−pp+qS˜w+q2,ω2=r1N˜wK1−br1K−N˜wθKS˜wr2S˜w−p2−pp+qS˜w+q2+θN˜wa+bN˜wa+bN˜w+cS˜w2mS˜w2k+N˜w2.ENS=N˜w,S˜w,0 is locally asymptotically stable if S˜w<δ/β and(i)
ω1<0 and ω2>0, or(ii)
ω1>0, ω12−4ω2<0 and α<α∗=2/πtan−1ω12−4ω2/ω1Proof.
(a)
In view of (28), around E0=0,0,0 we have(35)JE0=r1000−r2pq000−δ.The corresponding eigenvalues of (35) are λ1=r1>0,λ2=−r2p/q>0, and λ3=−δ<0. Thus, argλ1,2=0<απ/2 and argλ3=π>απ/2. Therefore, E0 is always a saddle point.(b)
The Jacobian matrix in (28) around EN=K,0,0 is(36)JEN=−r1−θKa+bK00−r2pq000−δ.Then, the corresponding eigenvalues of (36) are λ1=−r1<0,λ2=−r2p/q>0, and λ3=−δ<0. Since argλ2=0<απ/2 and argλ1,3=π>απ/2, therefore, EN is a saddle point.(c)
By substitutingES=0,Sw,0 to the Jacobian matrix in (28), we have(37)JES=a1100a21a22a2300a33,where(38)a11=r1−θSwa+cSw,a21=mSw2k2,a22=kr2Swp+q−mSwSw+q2kSw+q2,a23=−Swmk+β,a33=βSw−δ.The corresponding eigenvalues of (37) are λ1=r1−θSw/a+cSw, λ2=kr2Swp+q−mSwSw+q2/kSw+q2, and λ3=βSw−δ. Therefore, argλ1,2,3=π>απ/2, whenever r1<θSw/a+cSw, p+q<m/kr2Sw+q2, and δ>βSw.(d)
The Jacobian matrix in (28) around ESI=0,S¯w,I¯w is as follows:(39)JESI=b1100b21b22b230b320,where(40)b11=r1−δθaβ+cδ,b21=mδβkδ2+βq+r2δ−βr2pm+βkδ+βq,b22=βδr2p+qδ+βq2−mδβk,b23=−δmβk+1,b32=βkr2δ−βp−mδδ+βqm+βkδ+βq.The corresponding eigenvalues of (39) are λ1=r1−δθ/aβ+cδ and roots of the following equation:(41)λ2−βδr2p+qδ+βq2−mδβkλ+δβkβkr2δ−βp−mδδ+βqδ+βq=0.Equation (41) has the following eigenvalues:(42)λ2,3=12βδr2p+qδ+βq2−mδβk±Δ^w,where(43)Δ^w=βδr2p+qδ+βq2−mδβk2−4δβkβkr2δ−βp−mδδ+βqδ+βq.Notice that ifr1<δθ/aβ+cδ, then argλ1=π>απ/2. Therefore, the stability of ESI depends on λ2,3. If βr2p+q/δ+βq2<m/βk<r2δ−βp/δδ+βq, then argλ2,3>απ/2. Furthermore, if βr2p+q/δ+βq2>m/βk and Δ^w<0, then λ2,3 is a pair of complex conjugate eigenvalues with positive real parts. Thus, argλ2,3>απ/2 is achieved only if α<α∗=2/πtan−1Δ^w/βr2p+q/δ+βq2−m/βk. Using the Matignon’s condition (see Theorem 5), the theorem is completely proven.(e)
The Jacobian matrix (10) calculated atENS=N˜w,S˜w,0 is given by(44)JENS=c11c120c21c22c2300c33,where(45)c11=−r1N˜wK1−br1K−N˜wθKS˜w,c12=−θN˜wa+bN˜wa+bN˜w+cS˜w2,c21=mS˜w2k+N˜w2,c22=−r2S˜w−p2−pp+qS˜w+q2,c23=−S˜wmk+N˜w+β,c33=βS˜w−δ.
The eigenvalues of (44) are λ1=βS˜w−δ and the roots of the quadratic equation λ2−ω1λ+ω2=0 where(46)ω1=−r1N˜wK1−br1K−N˜wθKS˜w+r2S˜w−p2−pp+qS˜w+q2,ω2=r1N˜wK1−br1K−N˜wθKS˜wr2S˜w−p2−pp+qS˜w+q2+θN˜wa+bN˜wa+bN˜w+cS˜w2mS˜w2k+N˜w2.
IfS˜w<δ/β, then argλ1=π>απ/2. Therefore, the stability of ENS depends on λ2,3. If ω1<0 and ω2>0 then argλ2,3>απ/2. Furthermore, if ω1>0 and ω12<4ω2, then λ2,3 is a pair of complex conjugate eigenvalues. Thus, argλ2,3>απ/2 is attained if α<α∗=2/πtan−1ω12−4ω2/ω1. So, by the Matignon’s condition (see Theorem 5) system (4) exhibits locally asymptotically stable behavior around ENS.Theorem 11.
Stability condition of interior equilibrium point for weak Alee effectp<0. Suppose that(47)η1=r1Nw∗K1−br1K−Nw∗2θKSw∗−Sw∗r2p+qSw∗+q2−mk+Nw∗,η2=−Sw∗r1Nw∗K1−br1K−Nw∗2θKSw∗r2p+qSw∗+q2−mk+Nw∗+βIw∗mk+Nw∗+β−mθNw∗a+bNw∗Sw∗+Iw∗a+bNw∗+cSw∗2k+Nw∗2,η3=βr1Nw∗Sw∗Iw∗K1−br1K−Nw∗2θKSw∗mk+Nw∗+β,DwP=18η1η2η3+η1η22−4η3η13−4η23−27η32.E∗=Nw∗,Sw∗,Iw∗ is locally asymptotically stable if(i)
DwP>0,η1>0,η3>0, andη1η2>η3.(ii)
DwP<0,η1≥0,η2≥0,η3>0, and0<α<2/3.(iii)
DwP<0,η1<0,η2<0, and2/3<α<1.(iv)
DwP<0,η1>0,η2>0,η1η2=η3, and0<α<1.Proof.
The Jacobian matrix (10) evaluated at interior equilibrium point E∗=Nw∗,Sw∗,Iw∗ is given by(48)JE∗=d11d120d21d22d230d320,where(49)d11=−r1Nw∗K1−br1K−Nw∗2θKSw∗,d12=−θNw∗a+bNw∗a+bNw∗+cSw∗2,d21=mSw∗Sw∗+Iw∗k+Nw∗2,d22=Sw∗r2p+qSw∗+q2−mk+Nw∗,d23=−Sw∗mk+Nw∗+βd32=βIw∗.
The corresponding eigenvalues of (48) are the roots of the cubic equation Pλ=λ3+η1λ2+η2λ+η3=0 where the discriminant DwP of the cubic equation Pλ is as follows:(50)DwP=18η1η2η3+η1η22−4η3η13−4η23−27η32.
IfDwP>0, all the roots of the cubic equation are real; and there is only one real root and two complex conjugate roots if DwP<0. Therefore, the Routh-Hurwitz criterion for Caputo fractional-order [39] completes the proof of the stability condition for E∗.


### 3.5. Equilibrium Analysis in the Strong Allee Effect in Predator

In this subsection, the existence and stability of the equilibrium points of system (4) for the strong Alee effect p>0 are given. All positive equilibrium points of system (4) are as follows.(a)
Trivial equilibrium point:Ξ0=0,0,0, which always exists,(b)
Axial equilibrium point:ΞN=K,0,0, which always exists; and ΞS=0,Ss,0where Ss1,2=kr2−mq/2m±kr2−mq/2m2−kpr2/m. Assume that kr2>mq, the existence of ΞS is described as follows.(i)
Ifkr2−mq/2m2<kpr2/m, then ΞS does not exist.(ii)
Ifkr2−mq/2m2<kpr2/m, then there exists a unique ΞS=0,kr2−mq/2m,0(iii)
Ifkr2−mq/2m2<kpr2/m, then there exists twoΞS=0,kr2−mq/2m±kr2−mq/2m2−kpr2/m,0(c)
Planar equilibrium point:ΞSI=0,S¯s,I¯s exists if p<S¯s1−mδ+βq/βkr2, where(51)S¯s=δβ,andI¯w=1m+βkβkr2δ−βp−δmδ+βqβδ+βq;andΞNS=N˜s,S˜s,0 exists if S˜s1−mS˜s+q/kr2<p<S˜s, where N˜s=mS˜sS˜s+q/r2S˜s−p−k and S˜s are also all positive roots of quartic (30).(d)
Interior equilibrium point:Ξ∗=Ns∗,Ss∗,Is∗, where(52)Ns1,2∗=bK−a+cSs∗2b±Δs,Ss∗=δβ,Is∗=r2Ss∗−pk+Ns∗−mSs∗Ss∗+qSs∗+qm+βk+Ns∗,Δs=bK−a+cSs∗2b2+Kr1a+cSs∗−θSs∗br1.Assume thatSs∗>p and k+Ns∗>mSs∗Ss∗+q/r2Ss∗−p, then we have the following results.(a)
IfΔs<0, then system (4) has no interior equilibrium point.(b)
IfΔs=0 and aβ+cδ<bKβ, then system (4) has a unique interior equilibrium point Ξ1∗=Ns∗,δ/β,Is∗ where Ns∗=bKβ−aβ+cδ/2bβ and Is∗=βr2δ−βpk+Ns∗−δmδ+βq/βδ+βqm+βk+Ns∗.(c)
IfΔs>0 and aβ+cδ<minbKβ,δθ/r1, then system (4) has two positive interior equilibrium points Ξ1,2∗=Nsi∗,δ/β,Isi∗ for i=1,2, where(53)Ns1,2∗=bKβ−aβ+cδ2bβ±Δs,Is1,2∗=βr2δ−βpk+Ns1,2∗−δmδ+βqβδ+βqm+βkNs1,2∗,Δs=bKβ−aβ+cδ/2bβ2+Kr1aβ+cδ−δθ/br1β.Theorem 12.
The local stability of all equilibrium points of system (4) for the strong Alee effect p>0 is summarized as follows.(a)
Ξ0=0,0,0 is a saddle point.(b)
ΞN=K,0,0 is always locally asymptotically stable.(c)
ΞS=0,Ss,0 is locally asymptotically stable if r1<θSs/a+cSs, p+q<m/kr2Ss+q2, and Ss<δ/β.(d)
LetS¯s>p and Δ^s=βδr2p+q/δ+βq2−mδ/βk2−4δ/βkβkr2δ−βp−mδδ+βq/δ+βq. ΞSI=0,S¯s,I¯s is locally asymptotically stable if r1<δθ/aβ+cδ and(i)
βr2p+q/δ+βq2<m/βk<r2δ−βp/δδ+βq, or(ii)
βr2p+q/δ+βq2>m/βk, Δ^s<0, and α<α∗=2/πtan−1Δ^s/βr2p+q/δ+βq2−m/βk.(e)
Suppose that(54)ωs1=−r1N˜sK1−br1K−N˜sθKS˜s+r2S˜s−p2−pp+qS˜s+q2,ωs2=r1N˜sK1−br1K−N˜sθKS˜sr2S˜s−p2−pp+qS˜s+q2+θN˜sa+bN˜sa+bN˜s+cS˜s2mS˜s2k+N˜s2.ΞNS=N˜s,S˜s,0 is locally asymptotically stable if S˜s<δ/β and(i)
ωs1<0, or(ii)
ωs2>0,ωs12−4ωs2<0, andα<α∗=2/πtan−1ωs12−4ωs2/ωs1.(f)
Suppose that(55)ξ1=r1Ns∗K1−br1K−Ns∗2θKSs∗−Ss∗r2p+qSs∗+q2−mk+Ns∗,ξ2=−Ss∗r1Ns∗K1−br1K−Ns∗2θKSs∗r2p+qSs∗+q2−mk+Ns∗+βIs∗mk+Ns∗+β−mθNs∗a+bNs∗Ss∗+Is∗a+bNs∗+cSs∗2k+Ns∗2ξ3=βr1Ns∗Ss∗Is∗K1−br1K−Ns∗2θKSs∗mk+Ns∗+βDsP=18ξ1ξ2ξ3+ξ1ξ22−4ξ3ξ13−4ξ23−27ξ32.Ξ∗=Ns∗,Ss∗,Is∗ is locally asymptotically stable if(i)
DsP>0,ξ1>0,ξ3>0,andξ1ξ2>ξ3.(ii)
DsP<0,ξ1≥0,ξ2≥0,ξ3>0, and0<α<2/3.(iii)
DsP<0,ξ1<0,ξ2<0, and2/3<α<1.(iv)
DsP<0,ξ1>0,ξ2>0,ξ1ξ2=ξ3, and0<α<1.Theorem12 has similar proof to Theorem 10 and Theorem 11.


### 3.6. Hopf Bifurcation

In this subsection, we study the conditions of a Hopf bifurcation around the equilibrium point of system (4) when a parameter is varied. This bifurcation ensures a stability change when system (4) passes the critical value which coincides with the emergence of the limit cycle. Hopf bifurcation can occur both in the first-order systems and fractional-order systems. The fundamental difference between the two is the convergence of the limit set of solutions, known as the limit cycle, to the solution of the system. In the first-order systems, the limit cycle converges to the periodic solution; in the fractional-order systems, instead of converging to the periodic solution, the limit cycle converges to the periodic signal [40, 41].Let us consider the following three-dimensional fractional order system:(56)DtαCx=gx,0<α≤1,x∈ℝ3.According to Theorem5, the stability of the system dynamics is significantly affected by the order of the derivative α. Suppose that E∗ is an equilibrium point of system (19) and α∗ is the critical value of bifurcation parameter. The conditions under which system (19) undergoes a Hopf bifurcation driven by α near the equilibrium point E∗ are given in [42].(1)
The Jacobian matrix atE∗ has one real negative eigen value λ1<0 and a pair of complex conjugate eigen values λ2,3=ψ±iϕ (where ψ>0),(2)
χα∗=α∗π/2−argλ2,3α=0, and(3)
The transversality condition:∂χα/∂α|α=α∗≠0.Note that the critical value ofα∗ is the solution of χα∗=α∗π/2−argλ2,3α=0, i.e., α∗=2/πtan−1ϕ/ψ. Therefore, we have the following theorem.Theorem 13.
Existence of Hopf bifurcation driven byα around ESI, ENS, or E∗ for the weak Allee effect case.(i)
Suppose thatr1<δθ/aβ+cδ, βr2p+q/δ+βq2>m/βk, and Δ˜w<0. System (4) undergoes a Hopf bifurcation around ESI when α passes through(57)α∗=2πtan−1Δ˜wβr2p+q/δ+βq2−m/βk.(ii)
IfS˜<δ/β, ω1>0, and ω12−4ω2<0. Then, system (4) undergoes a Hopf bifurcation around ENS when parameter bifurcation α passes through(58)α∗=2πtan−1ω12−4ω2ω1.(iii)
Suppose that characteristic equationPλ=λ3+η1λ2+η2λ+η3=0 has one real negative eigen value and a pair of complex conjugate eigen values λ2,3=ζ1±iζ2 where ζ1>0. System (4) undergoes a Hopf bifurcation around the interior equilibrium point E∗=Nw∗,Sw∗,Iw∗ when α passes through(59)α∗=2πtan−1ζ12−4ζ2ζ1.Remark 14.
In addition to the existence of Hopf bifurcation driven by order of the fractional derivativeα, the transmission rate β also can be considered as a parameter bifurcation. Theorem 4 in [42] gives the existence condition for Hopf bifurcation by varying β. Since it is difficult to find out the exact critical values of β=β∗, therefore we numerically investigated the existence of Hopf bifurcation in system (4) driven by parameter β.


## 3.1. Non-Negativity and Boundedness of Solutions

From the biological point of view, we only concern to system (4) in the non-negative and bounded solutions. We consider ℝ+3≔N,S,I∈ℝ3:N≥0,S≥0,I≥0 as the set of all non-negative solutions. The following theorem shows the non-negativity and boundedness of the solutions of system (4).Theorem 8.
All solutions of the system (4) which initiate in ℝ+3 are non-negative and uniformly bounded.Proof.
LetN0>0,S0>0, and I0>0 and assume that Nt≥0 for t∈0,∞. Suppose the assumption is not true, then there exists t∗ which satisfies 0≤t<t∗ such that Nt>0 for t∈0,t∗, Nt∗=0, and Nt<0 for t>t∗. From the first equation of system (4), we have(13)DtαCNt∗|Nt∗=0=0.
According to Lemma3, we get Nt=0 for t>t∗ which contradicts to Nt<0,∀t>t∗. Therefore, we have Nt≥0 for t∈0,∞. Using the similar way, it can be shown that St≥0 and It≥0 for t∈0,∞. Hence, all solutions of the system (4) are non-negative.
Next, we show the boundedness of solutions of the system (4) by using the fractional comparison lemma. From the first equation of system (4), we obtain(14)DtαCNt+Nt=r1N1−NK−θNSa+bN+cS+N≤−r1N2K+r1+1N=−r1KN−Kr1+12r12+Kr1+124r1≤Kr1+124r1.
Based on Lemma4, we have(15)Nt≤N0−Kr1+124r1Eα−tα+Kr1+124r1,where Eα is the Mittag–Leffler function. Lemma 5 and Corollary 6 in [38] implies that Eα−tα⟶0 as t⟶∞, then it follows(16)0≤Nt≤Kr1+124r1=η1,ast⟶∞.
Thus, for anyε>0, Nt with initial condition in ℝ+3 is uniformly bounded to the region(17)Θ1=Nt≤η1+ε,ε>0.
DefineWt=St+It and adding the last two equations of system (4). Then, we get(18)DtαCWt≤r2S−pr2SS+q−mS2k+N−δI.
Forξ>0 and Nt≤η1+ε, we have(19)DtαCWt+ξWt≤−mk+η1+εS−r2+ξk+η1+ε2m2+r2+ξ2k+η1+ε4m+ξ−δI−pr2SS+q.
If we chooseξ<δ, then(20)DtαCWt+ξWt≤−mk+η1+εS−r2+ξk+η1+ε2m2+r2+ξ2k+η1+ε4m−pr2SS+q.≤r2+ξ2k+η1+ε4m−pr2SS+q.
Using the same argument as before and lettingε⟶0, then for the weak Alee effect p<0 we have the following equation:(21)0≤Wt≤r2+ξ2k+η14ξm−pr2ξq=η2,ast⟶∞.
On the other hand, for the strong Alee effectp>0 we have the following equation:(22)0≤Wt≤r2+ξ2k+η14ξm=η3,ast⟶∞.
Thus, for anyε>0, the combination of the solutions St and It with initial conditions started in ℝ+3 are uniformly bounded in the region(23)Θ2=St+It≤η2+ε,ε>0,p<0,St+It≤η3+ε,ε>0,p>0.


## 3.2. Existence and Uniqueness of the Solution

To investigate the existence and uniqueness of solution of system (4), we apply the locally Lipschitz condition in the region 0,T×ϒ, where ϒ=N,S,I∈ℝ3:maxN,S,I≤Φ, for sufficient large Φ and T<∞.Theorem 9.
For eachZ0=N0,S0,I0∈ϒ, there exists a unique solution Zt=Nt,St,It∈ϒ of system (4) with initial condition Z0, which is defined for all t≥0.Proof.
For anyZ,Z¯∈ϒ, we consider a mapping of a continuous function FZ=F1Z,F2Z,F3Z where(24)F1Z=r1N−r1N2K−θNSa+bN+cS,F2Z=r2S2S+q−pr2SS+q−mS2k+N−mSIk+N−βSI,F3Z=βSI−δI.
It follows that(25)FZ−FZ¯=F1Z−F1Z¯+F2Z−F2Z¯+F3Z−F3Z¯=|r1N−r1N2K−θNSa+bN+cS−r1N¯−r1N¯2K−θN¯S¯a+bN¯+cS¯|+|r2S2S+q−pr2SS+q−mS2k+N−mSIk+N−βSI−r2S¯2S¯+q−pr2S¯S¯+q−mS¯2k+N¯−mS¯I¯k+N¯−βS¯I¯|+βSI−δI−βS¯I¯−δI¯≤r1N−N¯+r1KN2−N¯2+θNSa+bN+cS−N¯S¯a+bN¯+cS¯+r2S2S+q−S¯2S¯+q+pr2SS+q−S¯S¯+q+mS2k+N−S¯2k+N¯+mSIk+N−S¯I¯k+N¯+βSI−S¯I¯+βSI−S¯I¯+δI−I¯≤r1N−N¯+r1KN+N¯N−N¯+θaN+bN¯Na2S−S¯+θaS¯+cS¯Sa2N−N¯+r2q2S¯SS−S¯+r2qS+S¯S−S¯+pr2qS−S¯+mk+N¯S+S¯k2S−S¯+mS¯2k2N−N¯+mk+N¯Ik2S−S¯+mN¯Sk2I−I¯+mS¯Ik2N−N¯+2βSI−I¯+2βIS−S¯+δI−I¯≤M1N−N¯+M2S−S¯+M3I−I¯≤MZ−Z¯,where(26)M1=r1+2r1K+θaΦ+θca2+2mk2Φ2,M2=pr2q+θa+2r2q+3mk+2βΦ+θba2+r2q2+3mk2Φ2,M3=δ+2βΦ+mk2Φ2,M=maxM1,M2,M3.
SinceFZ satisfies the Lipschitz condition with respect to Z, then system (4) with initial condition Z0 has a unique solution Zt∈ϒ for all t≥0. The proof is completed.


## 3.3. Existence and Stability of Equilibrium Points

To study the dynamical behaviors of system (4), we first investigate the equilibrium points of the system (4) which are the constant solutions to the following system:(27)r1N1−NK−θNSa+bN+cS=0,Sr2S−pS+q−mS+Ik+N−βSI=0,βSI−δI=0.Next, the local stability of equilibrium points of system (4) is evaluated by computing the eigenvalues of the Jacobian matrix of system (4) at the equilibrium point N∗,S∗,I∗ which is given as follows:(28)JN∗,S∗,I∗=J1,1J1,20J2,1J2,2J2,30J3,2J3,3,where(29)J1,1=r11−N∗K−r1N∗K−θS∗a+bN∗+cS∗+bθN∗S∗a+bN∗+cS∗2,J1,2=−θN∗a+bN∗+cS∗+cθN∗S∗a+bN∗+cS∗2J2,1=mS∗S∗+I∗k+N∗2J2,2=r2S∗−pS∗+q−mS∗+I∗k+N∗+S∗r2S∗+q−r2S∗−pS∗+q2−mk+N∗−βI∗J2,3=−mS∗k+N∗+βS∗J3,2=βI∗J3,3=βS∗−δ.


## 3.4. Equilibrium Analysis in the Case of the Weak Allee Effect in Predator

In this subsection, the existence and stability of the equilibrium points of the system (4) for the weak Alee effect p<0 are given. By solving system (9) simultaneously, all positive equilibrium points of system (4) are obtained as follows:(a)
Trivial equilibrium point:E0=0,0,0, which always exists.(b)
Axial equilibrium point:EN=K,0,0, which always exists; and ES=0,Sw,0, where Sw=kr2−mq/2m+kr2−mq/2m2−kpr2/m, which always exists.(c)
Planar equilibrium point:ESI=0,S¯w,I¯w where S¯w=δ/β and I¯w=1/m+βkβkr2δ−βp−δmδ+βq/βδ+βq exists if p<min0,δ/β1−mδ+βq/βkr2; and ENS=N˜w,S˜w,0 exists if p>S˜w1−mS˜w+q/kr2, where N˜w=mS˜wS˜w+q/r2S˜w−p−k and S˜w are all positive roots of the following quartic equation:(30)l1S4+l2S3+l3S2+l4S+l5=0,where(31)l1=mb+cr2mr1,l2=−bk−cp+ar2+mbqmr1−mb+cr2k+Kr2−mqr1+θKr22,l3=−p−bk+amr1r2−−bk−cp+ar2+mbqk+Kr2−mqr1+mb+cr2pk+Kr2r1−2θKr22p,l4=pr2k+Kr2−mq−bk+ar1+−bk−cp+ar2+mbqpk+Kr2r1+θKr22p2,l5=−p−bk+apk+Kr2r2r1.(d)
Interior equilibrium point:E∗=Nw∗,Sw∗,Iw∗, where(32)Nw∗=bK−a+cSw∗2b±Δw,Sw∗=δβ,Iw∗=r2Sw∗−pk+Nw∗−mSw∗Sw∗+qSw∗+qm+βk+Nw∗,Δw=bK−a+cSw∗2b2+Kr1a+cSw∗−θSw∗br1.Assume thatk+Nw∗>δmδ+βq/βr2δ−βp, then we have the following results:(a)
IfΔw<0, then system (4) has no interior equilibrium point.(b)
IfΔw=0 and aβ+cδ<bKβ, then system (4) has a unique interior equilibrium point E∗=Nw∗,δ/β,Iw∗ where Nw∗=bKβ−aβ+cδ/2bβ and Iw∗=βr2δ−βpk+Nw∗−δmδ+βq/βδ+βqm+βk+Nw∗.(c)
IfΔw>0 and aβ+cδ<minbKβ,δθ/r1, then system (4) has two positive interior equilibrium points E1,2∗=Nwi∗,δ/β,Iwi∗, i=1,2, where(33)Nw1,2∗=bKβ−aβ+cδ2bβ±Δw,Iw1,2∗=βr2δ−βpk+N^w1,2−δmδ+βqβδ+βqm+βk+N^w1,2,Δw=bKβ−aβ+cδ2bβ2+Kr1aβ+cδ−δθbr1β.Theorem 10.
The local stability of trivial, axial, and planar equilibrium points of system (4) for the weak Alee effect p<0 are summarized as follows.(a)
E0=0,0,0 is always a saddle point.(b)
EN=K,0,0 is always a saddle point.(c)
ES=0,Sw,0 is locally asymptotically stable if r1<θSw/a+cSw, p+q<m/kr2Sw+q2, and Sw<δ/β.(d)
LetΔ^w=βδr2p+q/δ+βq2−mδ/βk2−4δ/βkβkr2δ−βp−mδδ+βq/δ+βq. ESI=0,S¯w,I¯w is locally asymptotically stable if r1<δθ/aβ+cδ and one of the following conditions holds:(i)
βr2p+q/δ+βq2<m/βk<r2δ−βp/δδ+βq, or(ii)
βr2p+q/δ+βq2>m/βk, Δ^w<0, and α<α∗=2/πtan−1Δ^w/βr2p+q/δ+βq2−m/βk.(e)
Suppose that(34)ω1=−r1N˜wK1−br1K−N˜wθKS˜w+r2S˜w−p2−pp+qS˜w+q2,ω2=r1N˜wK1−br1K−N˜wθKS˜wr2S˜w−p2−pp+qS˜w+q2+θN˜wa+bN˜wa+bN˜w+cS˜w2mS˜w2k+N˜w2.ENS=N˜w,S˜w,0 is locally asymptotically stable if S˜w<δ/β and(i)
ω1<0 and ω2>0, or(ii)
ω1>0, ω12−4ω2<0 and α<α∗=2/πtan−1ω12−4ω2/ω1Proof.
(a)
In view of (28), around E0=0,0,0 we have(35)JE0=r1000−r2pq000−δ.The corresponding eigenvalues of (35) are λ1=r1>0,λ2=−r2p/q>0, and λ3=−δ<0. Thus, argλ1,2=0<απ/2 and argλ3=π>απ/2. Therefore, E0 is always a saddle point.(b)
The Jacobian matrix in (28) around EN=K,0,0 is(36)JEN=−r1−θKa+bK00−r2pq000−δ.Then, the corresponding eigenvalues of (36) are λ1=−r1<0,λ2=−r2p/q>0, and λ3=−δ<0. Since argλ2=0<απ/2 and argλ1,3=π>απ/2, therefore, EN is a saddle point.(c)
By substitutingES=0,Sw,0 to the Jacobian matrix in (28), we have(37)JES=a1100a21a22a2300a33,where(38)a11=r1−θSwa+cSw,a21=mSw2k2,a22=kr2Swp+q−mSwSw+q2kSw+q2,a23=−Swmk+β,a33=βSw−δ.The corresponding eigenvalues of (37) are λ1=r1−θSw/a+cSw, λ2=kr2Swp+q−mSwSw+q2/kSw+q2, and λ3=βSw−δ. Therefore, argλ1,2,3=π>απ/2, whenever r1<θSw/a+cSw, p+q<m/kr2Sw+q2, and δ>βSw.(d)
The Jacobian matrix in (28) around ESI=0,S¯w,I¯w is as follows:(39)JESI=b1100b21b22b230b320,where(40)b11=r1−δθaβ+cδ,b21=mδβkδ2+βq+r2δ−βr2pm+βkδ+βq,b22=βδr2p+qδ+βq2−mδβk,b23=−δmβk+1,b32=βkr2δ−βp−mδδ+βqm+βkδ+βq.The corresponding eigenvalues of (39) are λ1=r1−δθ/aβ+cδ and roots of the following equation:(41)λ2−βδr2p+qδ+βq2−mδβkλ+δβkβkr2δ−βp−mδδ+βqδ+βq=0.Equation (41) has the following eigenvalues:(42)λ2,3=12βδr2p+qδ+βq2−mδβk±Δ^w,where(43)Δ^w=βδr2p+qδ+βq2−mδβk2−4δβkβkr2δ−βp−mδδ+βqδ+βq.Notice that ifr1<δθ/aβ+cδ, then argλ1=π>απ/2. Therefore, the stability of ESI depends on λ2,3. If βr2p+q/δ+βq2<m/βk<r2δ−βp/δδ+βq, then argλ2,3>απ/2. Furthermore, if βr2p+q/δ+βq2>m/βk and Δ^w<0, then λ2,3 is a pair of complex conjugate eigenvalues with positive real parts. Thus, argλ2,3>απ/2 is achieved only if α<α∗=2/πtan−1Δ^w/βr2p+q/δ+βq2−m/βk. Using the Matignon’s condition (see Theorem 5), the theorem is completely proven.(e)
The Jacobian matrix (10) calculated atENS=N˜w,S˜w,0 is given by(44)JENS=c11c120c21c22c2300c33,where(45)c11=−r1N˜wK1−br1K−N˜wθKS˜w,c12=−θN˜wa+bN˜wa+bN˜w+cS˜w2,c21=mS˜w2k+N˜w2,c22=−r2S˜w−p2−pp+qS˜w+q2,c23=−S˜wmk+N˜w+β,c33=βS˜w−δ.
The eigenvalues of (44) are λ1=βS˜w−δ and the roots of the quadratic equation λ2−ω1λ+ω2=0 where(46)ω1=−r1N˜wK1−br1K−N˜wθKS˜w+r2S˜w−p2−pp+qS˜w+q2,ω2=r1N˜wK1−br1K−N˜wθKS˜wr2S˜w−p2−pp+qS˜w+q2+θN˜wa+bN˜wa+bN˜w+cS˜w2mS˜w2k+N˜w2.
IfS˜w<δ/β, then argλ1=π>απ/2. Therefore, the stability of ENS depends on λ2,3. If ω1<0 and ω2>0 then argλ2,3>απ/2. Furthermore, if ω1>0 and ω12<4ω2, then λ2,3 is a pair of complex conjugate eigenvalues. Thus, argλ2,3>απ/2 is attained if α<α∗=2/πtan−1ω12−4ω2/ω1. So, by the Matignon’s condition (see Theorem 5) system (4) exhibits locally asymptotically stable behavior around ENS.Theorem 11.
Stability condition of interior equilibrium point for weak Alee effectp<0. Suppose that(47)η1=r1Nw∗K1−br1K−Nw∗2θKSw∗−Sw∗r2p+qSw∗+q2−mk+Nw∗,η2=−Sw∗r1Nw∗K1−br1K−Nw∗2θKSw∗r2p+qSw∗+q2−mk+Nw∗+βIw∗mk+Nw∗+β−mθNw∗a+bNw∗Sw∗+Iw∗a+bNw∗+cSw∗2k+Nw∗2,η3=βr1Nw∗Sw∗Iw∗K1−br1K−Nw∗2θKSw∗mk+Nw∗+β,DwP=18η1η2η3+η1η22−4η3η13−4η23−27η32.E∗=Nw∗,Sw∗,Iw∗ is locally asymptotically stable if(i)
DwP>0,η1>0,η3>0, andη1η2>η3.(ii)
DwP<0,η1≥0,η2≥0,η3>0, and0<α<2/3.(iii)
DwP<0,η1<0,η2<0, and2/3<α<1.(iv)
DwP<0,η1>0,η2>0,η1η2=η3, and0<α<1.Proof.
The Jacobian matrix (10) evaluated at interior equilibrium point E∗=Nw∗,Sw∗,Iw∗ is given by(48)JE∗=d11d120d21d22d230d320,where(49)d11=−r1Nw∗K1−br1K−Nw∗2θKSw∗,d12=−θNw∗a+bNw∗a+bNw∗+cSw∗2,d21=mSw∗Sw∗+Iw∗k+Nw∗2,d22=Sw∗r2p+qSw∗+q2−mk+Nw∗,d23=−Sw∗mk+Nw∗+βd32=βIw∗.
The corresponding eigenvalues of (48) are the roots of the cubic equation Pλ=λ3+η1λ2+η2λ+η3=0 where the discriminant DwP of the cubic equation Pλ is as follows:(50)DwP=18η1η2η3+η1η22−4η3η13−4η23−27η32.
IfDwP>0, all the roots of the cubic equation are real; and there is only one real root and two complex conjugate roots if DwP<0. Therefore, the Routh-Hurwitz criterion for Caputo fractional-order [39] completes the proof of the stability condition for E∗.


## 3.5. Equilibrium Analysis in the Strong Allee Effect in Predator

In this subsection, the existence and stability of the equilibrium points of system (4) for the strong Alee effect p>0 are given. All positive equilibrium points of system (4) are as follows.(a)
Trivial equilibrium point:Ξ0=0,0,0, which always exists,(b)
Axial equilibrium point:ΞN=K,0,0, which always exists; and ΞS=0,Ss,0where Ss1,2=kr2−mq/2m±kr2−mq/2m2−kpr2/m. Assume that kr2>mq, the existence of ΞS is described as follows.(i)
Ifkr2−mq/2m2<kpr2/m, then ΞS does not exist.(ii)
Ifkr2−mq/2m2<kpr2/m, then there exists a unique ΞS=0,kr2−mq/2m,0(iii)
Ifkr2−mq/2m2<kpr2/m, then there exists twoΞS=0,kr2−mq/2m±kr2−mq/2m2−kpr2/m,0(c)
Planar equilibrium point:ΞSI=0,S¯s,I¯s exists if p<S¯s1−mδ+βq/βkr2, where(51)S¯s=δβ,andI¯w=1m+βkβkr2δ−βp−δmδ+βqβδ+βq;andΞNS=N˜s,S˜s,0 exists if S˜s1−mS˜s+q/kr2<p<S˜s, where N˜s=mS˜sS˜s+q/r2S˜s−p−k and S˜s are also all positive roots of quartic (30).(d)
Interior equilibrium point:Ξ∗=Ns∗,Ss∗,Is∗, where(52)Ns1,2∗=bK−a+cSs∗2b±Δs,Ss∗=δβ,Is∗=r2Ss∗−pk+Ns∗−mSs∗Ss∗+qSs∗+qm+βk+Ns∗,Δs=bK−a+cSs∗2b2+Kr1a+cSs∗−θSs∗br1.Assume thatSs∗>p and k+Ns∗>mSs∗Ss∗+q/r2Ss∗−p, then we have the following results.(a)
IfΔs<0, then system (4) has no interior equilibrium point.(b)
IfΔs=0 and aβ+cδ<bKβ, then system (4) has a unique interior equilibrium point Ξ1∗=Ns∗,δ/β,Is∗ where Ns∗=bKβ−aβ+cδ/2bβ and Is∗=βr2δ−βpk+Ns∗−δmδ+βq/βδ+βqm+βk+Ns∗.(c)
IfΔs>0 and aβ+cδ<minbKβ,δθ/r1, then system (4) has two positive interior equilibrium points Ξ1,2∗=Nsi∗,δ/β,Isi∗ for i=1,2, where(53)Ns1,2∗=bKβ−aβ+cδ2bβ±Δs,Is1,2∗=βr2δ−βpk+Ns1,2∗−δmδ+βqβδ+βqm+βkNs1,2∗,Δs=bKβ−aβ+cδ/2bβ2+Kr1aβ+cδ−δθ/br1β.Theorem 12.
The local stability of all equilibrium points of system (4) for the strong Alee effect p>0 is summarized as follows.(a)
Ξ0=0,0,0 is a saddle point.(b)
ΞN=K,0,0 is always locally asymptotically stable.(c)
ΞS=0,Ss,0 is locally asymptotically stable if r1<θSs/a+cSs, p+q<m/kr2Ss+q2, and Ss<δ/β.(d)
LetS¯s>p and Δ^s=βδr2p+q/δ+βq2−mδ/βk2−4δ/βkβkr2δ−βp−mδδ+βq/δ+βq. ΞSI=0,S¯s,I¯s is locally asymptotically stable if r1<δθ/aβ+cδ and(i)
βr2p+q/δ+βq2<m/βk<r2δ−βp/δδ+βq, or(ii)
βr2p+q/δ+βq2>m/βk, Δ^s<0, and α<α∗=2/πtan−1Δ^s/βr2p+q/δ+βq2−m/βk.(e)
Suppose that(54)ωs1=−r1N˜sK1−br1K−N˜sθKS˜s+r2S˜s−p2−pp+qS˜s+q2,ωs2=r1N˜sK1−br1K−N˜sθKS˜sr2S˜s−p2−pp+qS˜s+q2+θN˜sa+bN˜sa+bN˜s+cS˜s2mS˜s2k+N˜s2.ΞNS=N˜s,S˜s,0 is locally asymptotically stable if S˜s<δ/β and(i)
ωs1<0, or(ii)
ωs2>0,ωs12−4ωs2<0, andα<α∗=2/πtan−1ωs12−4ωs2/ωs1.(f)
Suppose that(55)ξ1=r1Ns∗K1−br1K−Ns∗2θKSs∗−Ss∗r2p+qSs∗+q2−mk+Ns∗,ξ2=−Ss∗r1Ns∗K1−br1K−Ns∗2θKSs∗r2p+qSs∗+q2−mk+Ns∗+βIs∗mk+Ns∗+β−mθNs∗a+bNs∗Ss∗+Is∗a+bNs∗+cSs∗2k+Ns∗2ξ3=βr1Ns∗Ss∗Is∗K1−br1K−Ns∗2θKSs∗mk+Ns∗+βDsP=18ξ1ξ2ξ3+ξ1ξ22−4ξ3ξ13−4ξ23−27ξ32.Ξ∗=Ns∗,Ss∗,Is∗ is locally asymptotically stable if(i)
DsP>0,ξ1>0,ξ3>0,andξ1ξ2>ξ3.(ii)
DsP<0,ξ1≥0,ξ2≥0,ξ3>0, and0<α<2/3.(iii)
DsP<0,ξ1<0,ξ2<0, and2/3<α<1.(iv)
DsP<0,ξ1>0,ξ2>0,ξ1ξ2=ξ3, and0<α<1.Theorem12 has similar proof to Theorem 10 and Theorem 11.


## 3.6. Hopf Bifurcation

In this subsection, we study the conditions of a Hopf bifurcation around the equilibrium point of system (4) when a parameter is varied. This bifurcation ensures a stability change when system (4) passes the critical value which coincides with the emergence of the limit cycle. Hopf bifurcation can occur both in the first-order systems and fractional-order systems. The fundamental difference between the two is the convergence of the limit set of solutions, known as the limit cycle, to the solution of the system. In the first-order systems, the limit cycle converges to the periodic solution; in the fractional-order systems, instead of converging to the periodic solution, the limit cycle converges to the periodic signal [40, 41].Let us consider the following three-dimensional fractional order system:(56)DtαCx=gx,0<α≤1,x∈ℝ3.According to Theorem5, the stability of the system dynamics is significantly affected by the order of the derivative α. Suppose that E∗ is an equilibrium point of system (19) and α∗ is the critical value of bifurcation parameter. The conditions under which system (19) undergoes a Hopf bifurcation driven by α near the equilibrium point E∗ are given in [42].(1)
The Jacobian matrix atE∗ has one real negative eigen value λ1<0 and a pair of complex conjugate eigen values λ2,3=ψ±iϕ (where ψ>0),(2)
χα∗=α∗π/2−argλ2,3α=0, and(3)
The transversality condition:∂χα/∂α|α=α∗≠0.Note that the critical value ofα∗ is the solution of χα∗=α∗π/2−argλ2,3α=0, i.e., α∗=2/πtan−1ϕ/ψ. Therefore, we have the following theorem.Theorem 13.
Existence of Hopf bifurcation driven byα around ESI, ENS, or E∗ for the weak Allee effect case.(i)
Suppose thatr1<δθ/aβ+cδ, βr2p+q/δ+βq2>m/βk, and Δ˜w<0. System (4) undergoes a Hopf bifurcation around ESI when α passes through(57)α∗=2πtan−1Δ˜wβr2p+q/δ+βq2−m/βk.(ii)
IfS˜<δ/β, ω1>0, and ω12−4ω2<0. Then, system (4) undergoes a Hopf bifurcation around ENS when parameter bifurcation α passes through(58)α∗=2πtan−1ω12−4ω2ω1.(iii)
Suppose that characteristic equationPλ=λ3+η1λ2+η2λ+η3=0 has one real negative eigen value and a pair of complex conjugate eigen values λ2,3=ζ1±iζ2 where ζ1>0. System (4) undergoes a Hopf bifurcation around the interior equilibrium point E∗=Nw∗,Sw∗,Iw∗ when α passes through(59)α∗=2πtan−1ζ12−4ζ2ζ1.Remark 14.
In addition to the existence of Hopf bifurcation driven by order of the fractional derivativeα, the transmission rate β also can be considered as a parameter bifurcation. Theorem 4 in [42] gives the existence condition for Hopf bifurcation by varying β. Since it is difficult to find out the exact critical values of β=β∗, therefore we numerically investigated the existence of Hopf bifurcation in system (4) driven by parameter β.


## 4. Global Stability

By using the appropriate Lyapunov function, we investigate the global stability of the stable equilibrium points of system (4), both for the weak Allee effect and the strong Allee effect cases in the predator.


### 4.1. For the Weak Allee Effect in Predator

Theorem 15.
ES=0,Sw,0 is globally asymptotically stable, if(60)−qSwη2+q−1≤p≤−η2+qr2Swr1K4+mSw2,Sw<δβ,hold.Proof.
LetNt,St,It be any positive solution of system (4). From (20), an ε>0 can be chosen such that(61)pr2Swη2+ε+q≤−r1K4+mSw2,pq≥1−Swη2+ε+q,Sw<δm+β<δβ.
Further, from Theorem8, we have 0≤St≤η2+ε. We now define a positive definite Lyapunov function at ES=0,Sw,0 as follows:(62)L1N,S,I=N+S−Sw−SwlnSSw+I.
By using Lemma6, the fractional time derivative of L1N,S,I along the solutions of system (4) is given by(63)DtαCL1N,S,I≤DtαCN+S−SwSDtαCS+DtαCI=r11−NK−θSa+bN+cSN+S−Swr2S−pS+q−mS+Ik+N−βI+βSI−δI=−r1KN−K22+r1K4−θNSa+bN+cS+r2S2S+q−pr2SS+q−r2SwSS+q+pr2SwS+q−mSwSk+N+mSw2k+N−mk+NS−Sw2−mSIk+N+mSwIk+N+βSwI−δI.≤r1K4+pr2Swη2+ε+q+mSw2+r21−pq−Swη2+ε+qS−mk+NS−Sw2+mSw+βSw−δI.
If condition in (61) is satisfied then we obtain the following:(64)DtαCL1N,S,I≤−mk+NS−Sw2.
In this case,DtαCL1N,S,I≤0,∀N,S,I∈ℝ+3, and DtαCL1N,S,I=0 at ES. Based on that, the only invariant set on which DtαCL1N,S,I=0 is the singleton ES. Using Lemma 7, the sufficient condition for the global asymptotic stability of ES is achieved.Theorem 16.
The planar equilibrium pointESI=0,S¯w,I¯w is globally asymptotically stable, if(65)−βη2+qδpr1K4+δ2mβ2k+δI¯w+δmη2βk≤r2≤pr2q+δr2βη2+q+δmβη1+k+βI¯w,holds.Proof.
Assume thatNt,St,It be the positive solution of system (4), so we can choose ε>0 such that(66)δpr2βη2+ε+q≤−r1K4+δ2mβ2k+δI¯w+δmβkη2+εr2≤pr2q+δr2βη2+ε+q+δmβk+η1+ε+βI¯w.
From Theorem8, we also have that Nt and St are bounded, i.e., 0≤Nt≤η1+ε and 0≤St≤η2+ε. Now, we define a positive definite Lyapunov function at ESI=0,S¯w,I¯w as(67)L2N,S,I=N+S−δβ−δβlnβSδ+I−I¯w−I¯wlnII¯w.
As before, we have the fractional time derivative ofL2N,S,I along the solutions of system (4) as follows:(68)DtαCL2N,S,I≤DtαCN+S−δ/βSDtαCS+I−I¯wIDtαCI=r11−NK−θSa+bN+cSN+S−δβr2S−pS+q−mS+Ik+N−βI≤r1K4+r2S2S+q−pr2SS+q−r2δ/βSS+q+pr2δ/βS+q−mδ/βSk+N+mδ/β2k+N−mS−δ/β2k+N−mSIk+N+mδ/βIk+N−βI¯wS+I¯wδ≤r1K4+δpr2βη2+ε+q+δ2mβ2k+δI¯w+δmβkη2+ε+r2−pr2q−δr2βη2+ε+q−δmβk+η1+ε−βI¯wS−mk+NS−δβ2−mk+NSI.
Observe that when the condition of (66) is fulfilled, then we have the following:(69)DtαCL2N,S,I≤−mk+NS−δβ2−mk+NSI.
One can easily show thatDtαCL2N,S,I≤0,∀N,S,I∈ℝ+3, and DtαCL2N,S,I=0 at ESI. Therefore, the only invariant set on which DtαCL2N,S,I=0 is the singleton ESI. Following Lemma 7, whenever ESI exists and(70)−βη2+qδpr1K4+δ2mβ2k+δI¯w+δmη2βk≤r2≤pr2q+δr2βη2+q+δmβη1+k+βI¯w,then it is globally asymptotically stable.Theorem 17.
The planar equilibrium point of system (4)ENS=N˜w,S˜w,0 is globally asymptotically stable if(71)N˜w≥K,r1≤S˜w−pr2/η2+q−mS˜w/kN˜wN˜w/K−1,r2≤mS˜w/k+η1−θN˜w/a1−S˜w/η2+q−p/q,δ≥mS˜wk+βS˜w.Proof.
LetNt,St,It be any positive solution of system (4). Then, from (24), an ε>0 can be chosen such that(72)1−N˜wK<0,r1N˜w2K−r1N˜w<−pr2S˜wη2+ε+q+mS˜w2k,mS˜wk+βS˜w<δθN˜wa+r2<r2S˜wη2+ε+q+pr2q+mS˜wk+η1+ε.
Also, from Theorem8, we have that Nt and St are bounded such that 0≤Nt≤η1+ε and 0≤St≤η2+ε. Let us consider a positive definite Lyapunov function at ENS=N˜w,S˜w,0 as follows:(73)L3N,S,I=N−N˜w−N˜wlnNN˜w+S−S˜w−S˜wlnNS˜w+I.
Based on Lemma6 and by calculating the fractional time derivative of L3N,S,I along the solutions of system (4),(74)DtαCL3N,S,I≤N−N˜wNDtαCN+S−S˜wSDtαCS+DtαCI=N−N˜wNr11−NK−θSa+bN+cSN+S−S˜wSr2S−pS+q−mS+Ik+N−βIS+βSI−δI=−r1KN−N˜w2+r1N−r1N˜wKN−r1N˜w+r1N˜w2K−θNSa+bN+cS+θN˜wSa+bN+cS+r2S2S+q−r2SS˜wS+q−pr2SS+q+pr2S˜wS+q−mk+NS−S˜w2−mS˜wk+NS+mS˜w2k+N−mISk+N+mIS˜wk+N+βS˜wI−δI≤−r1KN−N˜w2+r11−N˜wKN+r1N˜w2K−r1N˜w+pr2S˜wη2+ε+q+mS˜w2k+θN˜wa+r2−r2S˜wη2+ε+q−pr2q−mS˜wk+η1+εS−mk+NS−S˜w2−mk+NSI+mS˜wk+βS˜w−δI.
If condition in (72) is achieved, then we have the following:(75)DtαCL3N,S,I≤−r1KN−N˜w2−mk+NS−S˜w2−mk+NSI.
It is obvious thatDtαCL3N,S,I≤0,∀N,S,I∈ℝ+3, and DtαCL3N,S,I=0 at ENS. Hence, the singleton ENS is the only invariant set on which DtαCL3N,S,I=0. Thus, we can conclude that ENS is globally asymptotically stable which allowed by Lemma 7. □Theorem 18.
Suppose that(76)Ω1=2ar1Kaβ+bNw∗β−cδaβ+bNw∗β+2bδΩ2=2amk+η1a+bNw∗Ω3=βIw∗qr2−qΩ4=δm+βmη2+β2kIw∗η2+qβkr2.Ifθ≤minΩ1,Ω2 and Ω3≤−p≤Ω4, then E∗=Nw∗,Sw∗,Iw∗ is globally asymptotically stable.Proof.
LetNt,St,It be any positive solution of system (4), where Nt and St are bounded by Theorem 8. Further, one can use the following positive definite Lyapunov function to study the global stability of E∗=Nw∗,Sw∗,Iw∗:(77)L5N,S,I=a+bNw∗+cδβN−Nw∗−Nw∗lnNNw∗+S−δβ−δβlnNδ/β+I−Iw∗−Iw∗lnIIw∗.
Again, using Lemma6, we have the following equation:(78)DtαCL5N,S,I≤a+bNw∗+cδβN−Nw∗NDtαCN+S−δ/βSDtαCS+I−Iw∗IDtαCI=a+bNw∗+cδβN−Nw∗r11−NK−θSa+bN+cS+S−δβr2S−pS+q−mS+Ik+N−βI+I−Iw∗βS−δ≤−r1Ka+bNw∗+cδβN−Nw∗2+θa+bNw∗a+bN+cSN−Nw∗2+S−δ/β22+bθδβa+bN+cSN−Nw∗2−mk+NS−δβ2+r2S2S+q−pr2SS+q−δβr2SS+q+δβpr2S+q−mδβk+NS+δ2β2mk+N−mk+NSI+δβmk+NI−βSIw∗+δIw∗≤−r1Ka+bNw∗−cδβ−θa+bNw∗2a−bθδaβN−Nw∗2+θa+bNw∗2a−mk+η1+εS−δβ2+r2−pr2q−βIw∗S+δpr2βη2+ε+q+δ2mβ2k+δmη2+εβk+δIw∗.
Since we haveθ≤minΩ1,Ω2 and Ω3≤−p≤Ω4, then for an ε>0, we can choose the following condition:(79)r1Ka+bNw∗−cδβ≥θa+bNw∗2a+bθδaβ,θa+bNw∗2a≤mk+η1+ε,r2<pr2q+βIw∗,andδpr2βη2+ε+q≤−δ2mβ2k+δmη2+εβk+δIw∗,which implies DtαCL5N,S,I≤0,∀N,S,I∈ℝ+3, and DtαCL5N,S,I=0 at E∗. Based on that, the only invariant set on which DtαCL5N,S,I=0 is the singleton E∗. Employing the same argument, it follows that E∗ is globally asymptotically stable.


### 4.2. For the Strong Allee Effect in Predator

Theorem 19.
Ifp≥θK+r2η3+q/r2, then ΞN=K,0,0 is globally asymptotically stable.Proof.
LetNt,St,It be any positive solution of system (4). From the hypothesis of Theorem 12, an ε>0 can be chosen such that θK+r2≤r2p/η3+ε+q. In the strong Allee effect in predator, we have 0≤St≤η3+ε. Let us consider a positive definite Lyapunov function at ΞN=K,0,0 as follows:(80)L6N,S,I=N−K−KlnNK+S+I.
By taking the fractional time derivative ofL6N,S,I along the solutions of system (4) and using Lemma 6, one has the following equation:(81)DtαCL6N,S,I≤N−KNDtαCN+DtαCS+DtαCI=N−KNr11−NK−θSa+bN+cSN+Sr2S−pS+q−mS+Ik+N−δI≤−r1KN−K2+θK+r2−r2pη3+qS−δI.
Regarding thatθK+r2≤r2p/η3+ε+q, we have DtαCL6N,S,I≤0 for all N,S,I∈ℝ+3. Further, DtαCL6N,S,I=0 implies that N,S,I=K,0,0. Using Lemma 7, we conclude that E∗ is globally asymptotically stable.


## 4.1. For the Weak Allee Effect in Predator

Theorem 15.
ES=0,Sw,0 is globally asymptotically stable, if(60)−qSwη2+q−1≤p≤−η2+qr2Swr1K4+mSw2,Sw<δβ,hold.Proof.
LetNt,St,It be any positive solution of system (4). From (20), an ε>0 can be chosen such that(61)pr2Swη2+ε+q≤−r1K4+mSw2,pq≥1−Swη2+ε+q,Sw<δm+β<δβ.
Further, from Theorem8, we have 0≤St≤η2+ε. We now define a positive definite Lyapunov function at ES=0,Sw,0 as follows:(62)L1N,S,I=N+S−Sw−SwlnSSw+I.
By using Lemma6, the fractional time derivative of L1N,S,I along the solutions of system (4) is given by(63)DtαCL1N,S,I≤DtαCN+S−SwSDtαCS+DtαCI=r11−NK−θSa+bN+cSN+S−Swr2S−pS+q−mS+Ik+N−βI+βSI−δI=−r1KN−K22+r1K4−θNSa+bN+cS+r2S2S+q−pr2SS+q−r2SwSS+q+pr2SwS+q−mSwSk+N+mSw2k+N−mk+NS−Sw2−mSIk+N+mSwIk+N+βSwI−δI.≤r1K4+pr2Swη2+ε+q+mSw2+r21−pq−Swη2+ε+qS−mk+NS−Sw2+mSw+βSw−δI.
If condition in (61) is satisfied then we obtain the following:(64)DtαCL1N,S,I≤−mk+NS−Sw2.
In this case,DtαCL1N,S,I≤0,∀N,S,I∈ℝ+3, and DtαCL1N,S,I=0 at ES. Based on that, the only invariant set on which DtαCL1N,S,I=0 is the singleton ES. Using Lemma 7, the sufficient condition for the global asymptotic stability of ES is achieved.Theorem 16.
The planar equilibrium pointESI=0,S¯w,I¯w is globally asymptotically stable, if(65)−βη2+qδpr1K4+δ2mβ2k+δI¯w+δmη2βk≤r2≤pr2q+δr2βη2+q+δmβη1+k+βI¯w,holds.Proof.
Assume thatNt,St,It be the positive solution of system (4), so we can choose ε>0 such that(66)δpr2βη2+ε+q≤−r1K4+δ2mβ2k+δI¯w+δmβkη2+εr2≤pr2q+δr2βη2+ε+q+δmβk+η1+ε+βI¯w.
From Theorem8, we also have that Nt and St are bounded, i.e., 0≤Nt≤η1+ε and 0≤St≤η2+ε. Now, we define a positive definite Lyapunov function at ESI=0,S¯w,I¯w as(67)L2N,S,I=N+S−δβ−δβlnβSδ+I−I¯w−I¯wlnII¯w.
As before, we have the fractional time derivative ofL2N,S,I along the solutions of system (4) as follows:(68)DtαCL2N,S,I≤DtαCN+S−δ/βSDtαCS+I−I¯wIDtαCI=r11−NK−θSa+bN+cSN+S−δβr2S−pS+q−mS+Ik+N−βI≤r1K4+r2S2S+q−pr2SS+q−r2δ/βSS+q+pr2δ/βS+q−mδ/βSk+N+mδ/β2k+N−mS−δ/β2k+N−mSIk+N+mδ/βIk+N−βI¯wS+I¯wδ≤r1K4+δpr2βη2+ε+q+δ2mβ2k+δI¯w+δmβkη2+ε+r2−pr2q−δr2βη2+ε+q−δmβk+η1+ε−βI¯wS−mk+NS−δβ2−mk+NSI.
Observe that when the condition of (66) is fulfilled, then we have the following:(69)DtαCL2N,S,I≤−mk+NS−δβ2−mk+NSI.
One can easily show thatDtαCL2N,S,I≤0,∀N,S,I∈ℝ+3, and DtαCL2N,S,I=0 at ESI. Therefore, the only invariant set on which DtαCL2N,S,I=0 is the singleton ESI. Following Lemma 7, whenever ESI exists and(70)−βη2+qδpr1K4+δ2mβ2k+δI¯w+δmη2βk≤r2≤pr2q+δr2βη2+q+δmβη1+k+βI¯w,then it is globally asymptotically stable.Theorem 17.
The planar equilibrium point of system (4)ENS=N˜w,S˜w,0 is globally asymptotically stable if(71)N˜w≥K,r1≤S˜w−pr2/η2+q−mS˜w/kN˜wN˜w/K−1,r2≤mS˜w/k+η1−θN˜w/a1−S˜w/η2+q−p/q,δ≥mS˜wk+βS˜w.Proof.
LetNt,St,It be any positive solution of system (4). Then, from (24), an ε>0 can be chosen such that(72)1−N˜wK<0,r1N˜w2K−r1N˜w<−pr2S˜wη2+ε+q+mS˜w2k,mS˜wk+βS˜w<δθN˜wa+r2<r2S˜wη2+ε+q+pr2q+mS˜wk+η1+ε.
Also, from Theorem8, we have that Nt and St are bounded such that 0≤Nt≤η1+ε and 0≤St≤η2+ε. Let us consider a positive definite Lyapunov function at ENS=N˜w,S˜w,0 as follows:(73)L3N,S,I=N−N˜w−N˜wlnNN˜w+S−S˜w−S˜wlnNS˜w+I.
Based on Lemma6 and by calculating the fractional time derivative of L3N,S,I along the solutions of system (4),(74)DtαCL3N,S,I≤N−N˜wNDtαCN+S−S˜wSDtαCS+DtαCI=N−N˜wNr11−NK−θSa+bN+cSN+S−S˜wSr2S−pS+q−mS+Ik+N−βIS+βSI−δI=−r1KN−N˜w2+r1N−r1N˜wKN−r1N˜w+r1N˜w2K−θNSa+bN+cS+θN˜wSa+bN+cS+r2S2S+q−r2SS˜wS+q−pr2SS+q+pr2S˜wS+q−mk+NS−S˜w2−mS˜wk+NS+mS˜w2k+N−mISk+N+mIS˜wk+N+βS˜wI−δI≤−r1KN−N˜w2+r11−N˜wKN+r1N˜w2K−r1N˜w+pr2S˜wη2+ε+q+mS˜w2k+θN˜wa+r2−r2S˜wη2+ε+q−pr2q−mS˜wk+η1+εS−mk+NS−S˜w2−mk+NSI+mS˜wk+βS˜w−δI.
If condition in (72) is achieved, then we have the following:(75)DtαCL3N,S,I≤−r1KN−N˜w2−mk+NS−S˜w2−mk+NSI.
It is obvious thatDtαCL3N,S,I≤0,∀N,S,I∈ℝ+3, and DtαCL3N,S,I=0 at ENS. Hence, the singleton ENS is the only invariant set on which DtαCL3N,S,I=0. Thus, we can conclude that ENS is globally asymptotically stable which allowed by Lemma 7. □Theorem 18.
Suppose that(76)Ω1=2ar1Kaβ+bNw∗β−cδaβ+bNw∗β+2bδΩ2=2amk+η1a+bNw∗Ω3=βIw∗qr2−qΩ4=δm+βmη2+β2kIw∗η2+qβkr2.Ifθ≤minΩ1,Ω2 and Ω3≤−p≤Ω4, then E∗=Nw∗,Sw∗,Iw∗ is globally asymptotically stable.Proof.
LetNt,St,It be any positive solution of system (4), where Nt and St are bounded by Theorem 8. Further, one can use the following positive definite Lyapunov function to study the global stability of E∗=Nw∗,Sw∗,Iw∗:(77)L5N,S,I=a+bNw∗+cδβN−Nw∗−Nw∗lnNNw∗+S−δβ−δβlnNδ/β+I−Iw∗−Iw∗lnIIw∗.
Again, using Lemma6, we have the following equation:(78)DtαCL5N,S,I≤a+bNw∗+cδβN−Nw∗NDtαCN+S−δ/βSDtαCS+I−Iw∗IDtαCI=a+bNw∗+cδβN−Nw∗r11−NK−θSa+bN+cS+S−δβr2S−pS+q−mS+Ik+N−βI+I−Iw∗βS−δ≤−r1Ka+bNw∗+cδβN−Nw∗2+θa+bNw∗a+bN+cSN−Nw∗2+S−δ/β22+bθδβa+bN+cSN−Nw∗2−mk+NS−δβ2+r2S2S+q−pr2SS+q−δβr2SS+q+δβpr2S+q−mδβk+NS+δ2β2mk+N−mk+NSI+δβmk+NI−βSIw∗+δIw∗≤−r1Ka+bNw∗−cδβ−θa+bNw∗2a−bθδaβN−Nw∗2+θa+bNw∗2a−mk+η1+εS−δβ2+r2−pr2q−βIw∗S+δpr2βη2+ε+q+δ2mβ2k+δmη2+εβk+δIw∗.
Since we haveθ≤minΩ1,Ω2 and Ω3≤−p≤Ω4, then for an ε>0, we can choose the following condition:(79)r1Ka+bNw∗−cδβ≥θa+bNw∗2a+bθδaβ,θa+bNw∗2a≤mk+η1+ε,r2<pr2q+βIw∗,andδpr2βη2+ε+q≤−δ2mβ2k+δmη2+εβk+δIw∗,which implies DtαCL5N,S,I≤0,∀N,S,I∈ℝ+3, and DtαCL5N,S,I=0 at E∗. Based on that, the only invariant set on which DtαCL5N,S,I=0 is the singleton E∗. Employing the same argument, it follows that E∗ is globally asymptotically stable.


## 4.2. For the Strong Allee Effect in Predator

Theorem 19.
Ifp≥θK+r2η3+q/r2, then ΞN=K,0,0 is globally asymptotically stable.Proof.
LetNt,St,It be any positive solution of system (4). From the hypothesis of Theorem 12, an ε>0 can be chosen such that θK+r2≤r2p/η3+ε+q. In the strong Allee effect in predator, we have 0≤St≤η3+ε. Let us consider a positive definite Lyapunov function at ΞN=K,0,0 as follows:(80)L6N,S,I=N−K−KlnNK+S+I.
By taking the fractional time derivative ofL6N,S,I along the solutions of system (4) and using Lemma 6, one has the following equation:(81)DtαCL6N,S,I≤N−KNDtαCN+DtαCS+DtαCI=N−KNr11−NK−θSa+bN+cSN+Sr2S−pS+q−mS+Ik+N−δI≤−r1KN−K2+θK+r2−r2pη3+qS−δI.
Regarding thatθK+r2≤r2p/η3+ε+q, we have DtαCL6N,S,I≤0 for all N,S,I∈ℝ+3. Further, DtαCL6N,S,I=0 implies that N,S,I=K,0,0. Using Lemma 7, we conclude that E∗ is globally asymptotically stable.


## 5. Numerical Simulations

In this section, we demonstrate the numerical simulations based on the Adams–Bashforth–Moulton predictor-corrector method provided by Diethelm et al. [43] to verify the theoretical results established in the previous section. In addition, we also present the complex dynamics of system (4) such as the existence of bistability, forward, backward, saddle-node, and Hopf bifurcations as the effects of the disease and the fractional derivative.The first simulation is given to observe the role of the transmission rateβ to the dynamics of system (4) for the weak Allee effect case. We consider a set of hypothetical parameters in (26), where some of the parameters are taken from [28].(82)r1=0.5,K=5,θ=0.65,a=1,b=1,c=0.1,r2=0.1,k=1,m=0.01,p=−0.46,δ=0.2,q=3.In view of the existence conditions of the equilibrium points, Theorem10 and Theorem 11, we plot a bifurcation diagram for the weak Allee effect case as shown in Figure 1. Here, α is set at α=0.952 and varies the value of β in the range of 0,0.4. With these parameter values, system (4) always has the equilibrium points E0 and EN which is unstable according to Theorem 10 (a) and (b) and does not have a positive planar equilibrium point ENS. Therefore, we only observe the dynamics of system (4) around ES (the solid and dash black curve), ESI (the solid and dash green curve), and E∗ (the solid and dash blue curve).Figure 1
Bifurcation diagram of system (4) driven by β with parameter values (26) and α=0.952.In Figure1, when β<β1∗=0.0263, the axial equilibrium point ES is the only local asymptotically stable equilibrium point of system (4). A phase portrait of this behavior is shown in Figure 2(a), for example, β=0.01. One can easily confirm that 0.5=r1<θSw/a+cSw=2.8078, 2.54=p+q<m/kr2Sw+q2=11.2463, and 7.6049≈Sw<δ/β=20, which satisfy the locally asymptotically stable condition for ES=0,7.60487,0 as stated in Theorem 10(c). If we increase the value of parameter β such that β>β1∗=0.0263, ES loses its stability and the stable planar equilibrium point ESI arises via forward bifurcation. We confirm this behavior by numerical solutions in Figure 2(b). For β=0.05, we have the following conditions: 0.5=r1<δθ/aβ+cδ=1.8571 and 0.10367=βr2p+q/δ+βq2<m/βk=0.2<r2δ−βp/δδ+βq=0.31857 that fulfill the stability conditions of ESI=0,4,0.39524 in Theorem 10di but do not satisfy the stability conditions of ES in Theorem 10c. Since our numerical simulation set α=0.952, the stability of ESI is only attained at interval β∈β1∗,β4∗, where β4∗≈0.24. From the existence conditions of the interior point, we observe that when 0<β<β2∗≈0.1377, the interior equilibrium point E∗ does not exist. Next, based on the existence condition of the interior equilibrium point, two positive interior equilibrium points E1∗ and E2∗ which have different sign appear simultaneously after β across the critical value β2∗≈0.1377 via saddle-node bifurcation. Notice that for the set of parameters and β2∗<β<β3∗, the stability conditions of ESI and E∗ are satisfied. We plot a phase portrait for β=0.14 in Figure 2(c) to present the behavior, where E1∗=2.31322,1.42857,0.26803 and ESI=0,1.428571429,0.1890629800 are stable equilibrium points and E2∗=1.54392,1.42857,0.25727 is the unstable equilibrium point. It means system (4) experiences the bistability phenomenon in that interval.Figure 2
Phase portrait of system (4) with parameter values (26) and α=0.952. (a)β=0.01, (b)β=0.05, (c)β=0.14, (d)β=0.20, (e)β=0.25, and (f)β=0.40.
(a)(b)(c)(d)(e)(f)A further investigation shows that the stable branch interior pointE1∗ experiences a Hopf bifurcation when we increase β passes through β3∗≈0.17475. If we check the transversality condition through Figure 3(a), then we have dχβ/dββ=β3∗≈0.15602π≠0. To confirm Hopf bifurcation, we choose β=0.20 (see Figure 2(d)). One can easily check that the interior point E1∗=3.62407,1,0.16985 has one real negative eigen value, a pair of complex conjugate eigen values with positive real part and α>α∗=0.95044. Therefore, E1∗ fail to keep its stability which is indicated by the existence of limit cycle enclosing the interior point E1∗. Similar behavior around the interior point E1∗ still occur until β reaches β5∗≈0.26407 where the interior equilibrium point E1∗ retrieves again its stability via Hopf bifurcation for β larger than β5∗≈0.26407. Here, we verify the transversality condition using Figure 3(b) which satisfies dχβ/dββ=β5∗≈0.01915π≠0. Therefore, to illustrate this behavior, we depict a phase portrait of system (4) for β=0.4 in Figure 2(f). On the other hand, the unstable branch interior point E2∗ collides with the stable planar equilibrium point ESI when β attains β4∗ via backward bifurcation. Hence, the planar equilibrium point ESI also loses its stability, as shown in Figure 2(e). In this case, all existing equilibrium points of system (4) lose their stability. However, as we can see in Figure 2(e), all solutions of system (4) converge to a limit cycle around the interior point E1∗.Figure 3
Evolution ofχβ versus β with parameter values (26) and α=0.952.
(a)(b)Next, we discuss the contribution of transmission rateβ to the dynamics of system (4) under the strong Alee effect case p=0.4>0 with the same set of parameter values as in [28] subsection 6.2.(83)r1=0.5,K=5,θ=0.4,a=1,b=1,c=0.1,r2=0.1,k=1,m=0.05,p=0.4,δ=0.2,q=3,α=0.98.Notice that, for the strong Alee effect casep>0, ΞN=5,0,0 is always locally asymptotically stable according to Theorem 12b and Ξ0=0,0,0 is always a saddle point as stated in Theorem 12(a). In Figure 4(a), we portray a phase portrait of system (4) for β=0.073. System (4) has the following equilibrium points: Ξ0=0,0,0,ΞN=5,0,0,ΞNS1=4.58193,0.58956,0, and ΞNS2=2.2404,2.40121,0. Since ΞNS1 satisfies the stability conditions in Theorem 12(e) (i), then system (4) has bistability phenomenon. When we increase β to β=0.10, the planar equilibrium point ΞNS2=2.2404,2.40121,0 becomes a saddle point along with the emergence of the locally asymptotically stable interior point Ξ∗=3.16886,2,0.07154. In this case, we still have bistability condition, as we can see in Figure 4(b). Further, when we keep increasing β until β=0.73, the interior point Ξ∗ of system (4) does not exists anymore and the equilibrium point ΞN=5,0,0 becomes a unique stable equilibrium point of system (4), as shown in Figure 4(c).Figure 4
Phase portrait of system (4) for the strong Allee effect case with parameter values (27). (a)β=0.073, (b)β=0.10, and (c)β=0.73.
(a)(b)(c)Remark 20.
According to Figures1, 2, and 4, the transmission rate of disease β has a great impact to the dynamics of system (4). The numerical simulations indicate that system (4) could exhibit the forward, backward, saddle-node, and Hopf bifurcations driven by β. System (4) also performs a bi-stability phenomenon both for weak and strong Alee effects which suggest that solutions of system (4) are sensitive enough to the initial conditions. Figures 2(a), 2(b), and 2(f) confirm the global behavior of system (4) for the weak Allee effect.Remark 21.
From an ecological point of view, as the transmission rate of diseaseβ increases in the predator population, the number of susceptible predator population decreases; on the contrary, the prey population increases. When the Allee effect is weak, there is no possibility of extinction for the predator, and conversely when the Allee effect is strong, there is always a possibility of extinction for the predator. We also found for the same Allee threshold, models with disease would have a higher risk of predator extinction when compared to models without disease. This result suggests that to save all of the species from extinction with disease threatening one of the population, one should control the transmission rate of disease.
In the next simulation, we will show the influence of the order of fractional derivativeα on the dynamics of system (4) using the following parameter values:(84)r1=0.5,K=5,θ=0.4,a=1,b=1,c=0.1,r2=0.1,k=1,m=0.05,p=−0.46,δ=0.2,β=0.073q=3.
For the above parameter values, it is found that system (4) has four equilibrium points, where E0=0,0,0,EN=5,0,0,ES=0,0.58166,0 are saddle points, respectively, and ENS=1.21751,2.31593,0 is a conditionally locally asymptotically stable equilibrium point. It can be observed that system (4) undergoes a Hopf bifurcation at α∗=0.84124, as shown in Figure 5(a). We also can examine for α=α∗ the transversality condition is holds, i.e., ∂χα/∂α|α=α∗=π/2≠0. Then, we can observe that when α<α∗=0.84124, the planar equilibrium point ENS=1.21751,2.31593,0 is locally asymptotically stable as exhibited in Figure 5(b) for α=0.82. Otherwise, when α>α∗, there exists a limit cycle around ENS=1.21751,2.31593,0 which shows that both prey and susceptible predator populations are fluctuating, as shown in Figure 5(c).
We next exhibit the existence of two limit cycles as the solutions of system (4) by varying the order of the derivative α at α=0.85,0.96, and 0.98. Using the set of parameters in (26), except parameter β fixed at β=0.22, system (4) has six equilibrium points, i.e., E0=0,0,0, EN=5,0,0, ES=0,7.60487,0, ESI=0,0.90909,0.11275, E1∗=3.78913,0.90909,0.14915, and E2∗=0.11996,0.90909,0.11753. One can easily check that E0,EN,ES, and E2∗ are unstable, while the stability of ESI and E1∗ depends on the order of the derivative α. When α=0.85, system (4) undergoes the bistability phenomenon which occurs around the interior equilibrium point E1∗ and the planar equilibrium point ESI as we can see in Figure 6(a). This is confirmed by the real negative eigen value of E1∗ and ESI, also the critical value of the equilibrium points α∗ are 0.9505 and 0.9729, respectively, which are greater than α=0.85. The solutions of system (4) with two close enough initial conditions converge to different equilibrium points. When we increase α to 0.96, ESI is still locally asymptotically stable since0.9729=α∗>α=0.96, while E1∗ becomes unstable via a Hopf bifurcation and limit cycle appears enclosing the interior point E1∗, as shown in Figure 6(b). When we continue to increase α at α=0.98, following E1∗, then ESI loses its stability too through a Hopf bifurcation and two stable limit cycles appear instead in system (4), see Figure 6(c). These numerical simulations show that system (4) could exhibit two Hopf bifurcations simultaneously which is controlled by the order of the derivative α.Figure 5
Dynamics of system (4) for the weak Allee effect case with parameter values (28). (a) Bifurcation diagram, (b) α = 0.82, and (c) α = 0.98.
(a)(b)(c)Figure 6
The existence of two limit cycles for the weak Allee effect case with parameter values (26) and β=0.22. (a)α=0.85, (b)α=0.96, and (c)α=0.98.
(a)(b)(c)Remark 22.
From the ecological point of view, a phenomenon which is given in Figure6(c) show that if the initial value of prey population is relatively small then the susceptible and infected predators will oscillate even in the absence of prey population. On the other side, when the initial prey population is relatively large, all populations may oscillate for a long period of time.Remark 23.
The analytical results and numerical simulations of the model in the absence of disease in predator population (system (1)) had been found by Rahmi et al. [28]. Numerically, the recent model with disease in predator population has more rich dynamical behaviors than the former model. One of the reasons comes from varying the transmission rate of disease parameter. When there is no disease, there is only bistability, forward and Hopf bifurcations, meanwhile there is bistability, saddle-node, backward, forward, and Hopf bifurcations for the model with disease. More specifically, the proposed model exhibits a new phenomenon in the fractional order system, which is the existence of two limit cycles driven by the order of derivative as given in Figure 6(c).


## 6. Conclusion

In this paper, we have merged a predator-prey model and epidemiology model into the eco-epidemiological model. In ecological literature, there is much evidence that the double Allee effect may be acting on a single population. Moreover, it is natural to assume that there is a deadly infectious disease spreading within the population. Thus, we have developed an eco-epidemiological model incorporating the double Allee effect and disease spread on the predator population. This eco-epidemiological model has been modelled by a system with fractional-order differential equations in Caputo sense. We have split the predator population into two subpopulations: the susceptible predator and the infected predator. We proved that the solution of system (4) exists uniquely and whenever we started with a positive initial condition then all solutions remain non-negative and bounded. We showed that our system (4) has four types of equilibrium points, i.e., the trivial, two axial, two planar, and the interior equilibrium points. The trivial equilibrium point for both the strong and weak Allee effect cases is always unstable, which means that there is no condition of extinction of all populations in the future. The first axial equilibrium point (the predator extinction point) is always unstable for the weak Allee effect case, while for the strong Alee effect case, it is always stable. It means that the predator population can go to extinction in the strong Allee effect, which is contrary to the weak Allee effect case. The second axial equilibrium point (the prey and infected predator extinction point) is conditionally stable. We have observed that the dynamics of the axial equilibrium point do not depend on the order of the fractional derivative. Then, we have two planar equilibrium points (the prey extinction point and the infected predator extinction point) which are conditionally stable. The interior equilibrium point is also conditionally stable for both weak and strong Allee effects. The stability conditions of the last two types of equilibrium points show that the order of the fractional derivative α could affect its stability. By choosing an order α which is smaller than the critical value of α∗, then the equilibrium point may be stabilized. Numerical simulations show that the dynamics of system (4) have corresponded well with our theoretical results. It has been observed numerically that the transmission rate β takes a vital role to control the dynamical behavior of system (4). The combination of disease and double Allee effects on predators may accelerate the decline in predator populations which ultimately increases the risk of extinction. We also show numerically the existence of the Hopf, forward, backward, and saddle-node bifurcations driven by the transmission rate β. Further, the bistability phenomenon also appears in both the strong and weak Allee effect cases. All numerical results indicate that each population in system (4) could be extinct or survive or oscillate depending on the parameter values and initial conditions. In [28], the authors studied the dynamics of the model with no disease in the predator. They also carried out some numerical simulations to show how the capturing rate, Allee threshold, and the order of the derivative affect the dynamics of the model. The comparative results of these two models, with and without disease in the predator population are given in Remark 23. However, it could be interesting if the eco-epidemiological model with the double Allee effect and disease in predators also considered vaccination strategies as the approach to reduce population extinction risks. For example, there are African wild dogs that are endangered carnivores which get vaccination coverage for rabies virus and canine distemper virus (CDV) [44].



---

*Source: 5030729-2023-02-03.xml*