A fractional-order approach to predator-prey interactions: Modeling fear and disease dynamics with memory effects

TL;DR

This paper introduces a new predator-prey model using fractional calculus to study the effects of fear and disease with memory in ecological systems.

Contribution

The model introduces fear effects on recruitment, reinfection dynamics, selective predation on infected prey, and memory effects via fractional derivatives.

Findings

The model identifies three equilibrium points representing disease and predator extinction, predator-free, and co-existence states.

Local and global dynamics are analyzed using Matignon’s condition, Routh-Hurwitz criterion, and Lyapunov functions.

Numerical simulations confirm forward and Hopf bifurcations showing dynamic changes with parameter variation.

Abstract

The dynamical behaviors of a predator-prey model with fear effect and disease on prey are studied by employing the fractional-order derivative with a power-law kernel as the operator. The proposed model introduces four novel aspects: the impact of fear on a constant recruitment rate, previously unexplored in the literature; reinfection dynamics where infected prey can re-acquire the disease; selective predation by predators on infected prey due to pathogen-induced vulnerability; and the used of Caputo fractional derivative to include the memory effect. A deterministic approach is provided to establish the mathematical model including its validity by showing the existence, uniqueness, non-negativity, and boundedness. Three types of equilibrium points are obtained which represent the extinction of disease and predator, the predator-free, and the co-existence points. The local dynamics are investigated using Matignon’s condition and generalized Routh-Hurwitz criterion for the Caputo fractional-order model. The Volterra, quadratic, and linear functions are utilized to construct the Lyapunov function, as well as the LaSalle invariance principle, to show global dynamics. Some phenomena are provided namely forward and Hopf bifurcations to show the change of the dynamics when a parameter is varied. These results are supported by numerical ways using the predictor-corrector scheme.

Figures

Bifurcation diagram and phase portraits of model (6) using parameter values: Λ = 10, μ = 0.5, β = 1.5, ω = 0.5, δ = 0.5, m = 4, a = 1, n = 1, d = 0.5, α = 0.95 and varying the value of k.

Phase portraits of model (6) using parameter values: Λ = 10, k = 0.4, μ = 0.5, ω = 0.5, δ = 0.5, m = 4, a = 1, n = 1, d = 0.5, α = 0.95, and varying the value of β.

Phase portraits of model (6) using parameter values: λ = 10, k = 0.4, μ = 0.5, β = 1.5, ω = 0.5, δ = 0.5, m = 4, a = 1, n = 1, d = 0.5, and varying the value of α.

Equations

Full text

1 Introduction

In the prey-predator interaction, predation could be seen as direct killing (lethal effects) or indirect killing (nonlethal effects). The study on predation as direct killing in prey-predator systems has been widely published since Lotka and Volterra proposed the first prey-predator model considering various functional responses; see, for example, [1–4]. However, some experimental investigations show that the indirect killing risk, also known as the fear effect, is sometimes more crucial than direct killing [5–7]. Zanette et al. [8] have manipulated song sparrows during an entire breeding season and found that the song sparrows’ reproductive success was reduced to 40% only by broadcasting predator sounds around all protected nests. Fear of predators caused the female song sparrows to lay fewer eggs affecting the birth rate, hatch fewer of those eggs (related to survival rate), and lose more nestlings to death in the nest (related to survival rate). Other evidence given by Elliott, et al. [9] which conducted an observation on the common fruit fly (Drasophila melanogaster) with the scent of its predator, mantid (Tenodera sinensis). The results show that by the end of the non-breeding season, the flies lost 60 % more mass and had a 64 % higher mortality rate while in the breeding season, they produced 20 % fewer offspring with a 9 % lower maturity weight when compared to the control. Moreover, how fear can suppress the reproductive output and offspring quality of species can also be found in wild rodents [10], female mice [11], Brandt’s voles [12], mountain goat [13], and bugs [14]. Based on all of this evidence, it is necessary to include both the direct and indirect killing effects of predators to characterize the dynamics of prey-predator interaction.

In 2016, Wang et al. [15] formulated the first predator-prey model that investigated the cost of fear in the prey population. They showed that high levels of fear could stabilize the predator-prey system by ruling out the existence of periodic solutions, whereas low levels of fear could create multiple limit cycles. The same result was also performed by Sarkar and Khajanchi [1]. Pal et al. [16] developed and analyzed a predator-prey model with hunting cooperation among predators that induces fear in the prey population. They showed how an increase in hunting cooperation might cause the system to destabilize and lead to periodic solutions through a Hopf-bifurcation. Moreover, Das and Samanta [17] work on the dynamical behavior of the prey-predator model with fear effect and Holling type IV functional response incorporating the memory effect. Analytical results show the local stability with the occurrence of saddle-node and Hopf bifurcation. Wu et al. [18] investigated the Leslie-Gower model incorporating the fear effect and nonlinear harvesting on prey species, and dynamical analysis shows the complex results, such as the existence of saddle-node, Hopf and Bogdanov-Taken bifurcation. Generally, the prey-predator models under consideration can be stabilized or destabilized by the fear effect.

Recently, the evolution of mathematical modeling in the context of prey-predator interactions and how pathogen transmission affects one [19,20] or both species [21,22] has become the interesting topic in mathematical biology. In ecological communities, the predators can either positively or negatively influence disease dynamics through their effects on prey density and traits [23]. The positive impact occurs when the pathogen in infected prey could infect the predator as well. Conversely, if predation selectively removes the infected prey without affecting the predator then there is the possibility to reduce the pathogen in the prey population [24,25]. In this case, disease has a negative effect on the dynamics of the predator-prey model. From the view of dynamical behavior, the spread of pathogens in a Lotka-Volterra predator-prey model causes the existing prey-predator cycle to be disrupted, and the interactions will become more complex [26].

The further exploration study of the intermingling between a prey-predator model incorporating fear effect and the spread of disease either in prey or predator population has gained widespread concern in order to protect endangered species. Sarkar and Khajanchi [27] construct the model by assuming that the prey population grows logistically in the absence of predator species and that predators feed both susceptible and infected prey. They revealed that increasing levels of fear could stabilize the system toward the equilibrium point from a limit cycle oscillation condition. In this case, the disease still exists in the system but the amplitude of the infected prey decreases. Further, as the rate of infection in prey increases to a certain level the predator becomes extinct though predators are not infected by the disease. Hossain, et al. [28] also investigated the impact of fear in an eco-epidemiological model which assumes the prey grows logistically. The transmission rate of disease in the prey population follows a general disease transmission term. Kashyap et al. [29] work on the scenario of how the fear response in prey with disease could affect the mass mortality of predators (case in the Salton sea between tilapia fish and pelican birds) using the Caputo fractional derivative. On the other hand, Barman et al. [30] presented the eco-epidemiological model with infected disease in predator populations whether the fear induced by susceptible and infected predators at different levels. Zhang [26] also established the prey-predator model with standard incidence term as the transmission mechanism of infectious disease in the predator population.

Motivated by the above literature, in this paper, we propose a model addressing these key contributions:

The impact of fear on the constant recruitment rate which never considered previous work given by Moustafa et al. [31].Infected prey can be able to reinfection the disease.Since infected prey typically have pathogen-induced morbidity that makes them more vulnerable to predation, predators exclusively feed the infected ones [32–34].The Caputo fractional differential equation is employed to involve memory and hereditary properties which are nonlocal [35–40].

The following is the structure for the following content in this paper: Some preliminaries about the fractional calculus are recalled in Sect 2. In Sect 3, we formulate the eco-epidemiological model under the influence of an SIS infectious disease in the prey population by using the Caputo FDE. In Sect 4, we explore some qualitative analysis of the model by showing their biological validity, the existence of equilibrium points, local and global dynamics, and the occurrence of bifurcation phenomena. The numerical simulation of the model is given in Sect 5. We end the paper by summarizing the results in Sect 6.

2 Preliminaries

Some advantageous tools are provided below which consist of definitions, lemmas, and theorems related to fractional calculus and their application to mathematical modeling. We also propose some references such as [41–44] to obtain more in-depth information about the Caputo FDE.

Definition 1. [43] The Caputo FDE of order $[eqn]$ is defined by

$$[eqn]$$

where the function $[eqn]$ with $[eqn]$ is a space of n times continuously differentiable function and $[eqn]$ is the Gamma function.

Lemma 1. [45] Consider the following Caputo FDE.

$$[eqn]$$

with initial condition $[eqn]$ and $[eqn]$ . The Caputo FDE (1) has a unique solution in Φ when w(z,t) satisfies the local Lipschitz condition i.e

$$[eqn]$$

for any $[eqn]$ where $[eqn]$ is a Lipschitz constant.

Theorem 1. [46] The Caputo FDE (1) has non-negative solution for any non-negative initial condition as $[eqn]$ if (i) the function $[eqn]$ satisfies locally Lipschitz condition; and (ii) for $[eqn]$ , the Caputo FDE (1) has the non-negative solution.

Lemma 2. [47] Let $[eqn]$ satisfies

$$[eqn]$$

Then, the following inequality holds.

$$[eqn]$$

Theorem 2 (Matignon Condition [43,48]). A point $[eqn]$ is called the equilibrium point of the Caputo FDE (1) if $[eqn]$ . Moreover, the equilibrium point $[eqn]$ is said locally asymptotically stable (LAS) if all eigenvalues $[eqn]$ of the Jacobian matrix $[eqn]$ satisfy $[eqn]$ . It is said a saddle-point if there exist $[eqn]$ , $[eqn]$ such that $[eqn]$ and $[eqn]$ .

Theorem 3. [49] Consider the Caputo FDE (1) with n = 3 and the Jacobian matrix evaluated at the equilibrium point $[eqn]$ has polynomial characteristic

$$[eqn]$$

where has the discriminant $[eqn]$ . All eigenvalues of the Jacobian matrix given by the polynomial characteristic (2) satisfies $[eqn]$ if

(i) $[eqn]$ , c1>0*, c3>0, and* $[eqn]$ , or(ii) $[eqn]$ , $[eqn]$ , $[eqn]$ , c3>0 and $[eqn]$ , or(iii) $[eqn]$ , c1>0*, c2>0, and* $[eqn]$ .

If $[eqn]$ , c1<0*, c2<0, and* $[eqn]$ then all eigen values satisfies $[eqn]$ .

Lemma 3 (Quadratic-type Lyapunov functions [50]). For $[eqn]$ , the following inequality holds.

$$[eqn]$$

for all $[eqn]$ and $[eqn]$ .

Lemma 4 (Volterra-type Lyapunov functions [51]). Let $[eqn]$ and $[eqn]$ . For $[eqn]$ , the following inequality holds.

$$[eqn]$$

for all $[eqn]$ and $[eqn]$ .

Lemma 5. [52] Suppose that Φ is a bounded closed set and every solution of the Caputo FDE (1) which has the initial condition in Φ for $[eqn]$ . If there exists $[eqn]$ with continuous first order partial derivatives satisfies

$$[eqn]$$

then every solution z(t) originating in Ψ tends to $[eqn]$ as $[eqn]$ , where $[eqn]$ is the largerst invariant set of $[eqn]$ . If $[eqn]$ is the equilibrium point of eq. 1 which satisfies this condition, then $[eqn]$ is said globally asymptotically stable (GAS).

3 Model formulation

In this paper, we construct the prey-predator model where the prey population follows the SIS (Susceptible-Infected-Susceptible) type and the predator-induced fear effect on the prey population. We modify the eco-epidemic model proposed by Moustafa et al. [31]. The model assumes that there are two populations namely prey and predator where prey is divided into two compartments namely the susceptible and infected prey which respectively denoted by S and I and the predator denoted by P which only targeting the infected prey for food with Holling type-2 for the predatio functional response. As a result, we have the following model.

$$[eqn]$$

where the biological interpretation of each variable and parameter is given in Table 1. To involve the fear effect, we multiply the recruitment rate into susceptible prey population by the fear factor term $[eqn]$ where the parameter k represents the level of fear that reduces the birth rate of prey [53,54]. As a result, we have

Table 1: Biological interpretation of variables and parameters of model (6).

$$[eqn]$$

To include the memory effect to the model, we replace the first-order derivative at the left-hand side of model (4) with the fractional-order derivative $[eqn]$ along by rescale the dimension of time of the parameters to maintain their consistency. Therefore, we have

$$[eqn]$$

Using Caputo fractional derivative given by Definition 1 and rewriting the parameter values such that $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ , we have

$$[eqn]$$

In the next section, we will study the validity of the model, the dynamical behaviors, and some exploration in numerical simulations.

4 Qualitative analysis

4.1 Biological validity

We present this subsection to give the biological validity of the model by showing that each solution given by model (6) with non-negative initial condition will always exists, unique, non-negative, and bounded in the biological region

$$[eqn]$$

We provide the following theorems to show the existence, uniqueness, non-negativity, and boundedness of the solution of model (6).

Theorem 4. Let

$$[eqn]$$

Model (6) with initial condition in Φ has unique solution in $[eqn]$ .

Proof: Let $[eqn]$ and $[eqn]$ are solutions of model (6). For $[eqn]$ , we have

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ . According to Lemma 1, the function $[eqn]$ satisfies the locally Lipschitz condition with respect to the variables $[eqn]$ with $[eqn]$ is the Lipschitz constant. Thus, the model (6) has the unique solution in $[eqn]$ for each initial condition in Φ. □

Theorem 5. The solution S(t), I(t), and P(t) of model (6) always stays on $[eqn]$ for any $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ .

Proof: For $[eqn]$ , model (6) can be rewritten as follows.

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , and $[eqn]$ . From the second and the third equations on (7), we get solutions

$$[eqn]$$ $$[eqn]$$

Since the initial conditions $[eqn]$ and $[eqn]$ , we confirm that $[eqn]$ and $[eqn]$ as $[eqn]$ . From the non-negativity of I(t) and P(t), we derive the first equation on (7) as

$$[eqn]$$ $$[eqn]$$

Thus, we have $[eqn]$ which always non-negative since $[eqn]$ . Furthermore, based on the proof of Theorem 4, we also have $[eqn]$ satisfies the locally Lipschitz condition with respect to the variables $[eqn]$ . We assure from Theorem 1 that the model (6) always stays on $[eqn]$ for any non-negative initial condition. □

Theorem 6. Let $[eqn]$ . The solution of model (6) will always bounded.

Proof: We define the positive function as follows.

$$[eqn]$$

by computing the Caputo fractional derivative of eq. 8, we obtain

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Let $[eqn]$ . we achieve

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

and hence we have $[eqn]$ . According to Lemma 2, the following inequality holds.

$$[eqn]$$

where E is the Mittag-Leffler function. According to the convergence of the Mittag-Leffler function [43], we have $[eqn]$ as $[eqn]$ . Thus, for any non-negative initial condition, the solution is confined to $[eqn]$ where

$$[eqn]$$

This completes the proof. □

4.2 Feasible equilibrium points

In this subsection, all feasible equilibrium points are investigated including their existence condition. Following Theorem 2, the equilibrium points are obtained by solving $[eqn]$ as well as finding the positive solutions of the following equations.

$$[eqn]$$

As a result, we have three equilibrium points as follows.

(i) The prey-disease-predator-free point (PDPF) given by

$$[eqn]$$

which describes the condition where only the susceptible prey exists in the ecosystem while the disease disappears and the predator is extinct. To include the impact of the spread of the disease on prey, we compute the basic reproduction number $[eqn]$ using the next generation matrix as in [55,56]. We first consider the new infection term $[eqn]$ and the remaining transfer terms $[eqn]$ given by the second equation on (6) as follows.

$$[eqn]$$

Therefore, we have the basic reproduction number by solving $[eqn]$ at $[eqn]$ . The following term is acquired.

$$[eqn]$$

(ii) The predator-free point (PFP) represents the condition when only susceptible and infected prey exists in the ecosystem while the predator becomes extinct. The PFP is defined by

$$[eqn]$$

(iii) The co-existence point (CEP) is given by $[eqn]$ where $[eqn]$ , $[eqn]$ , and $[eqn]$ is the positive root respect to P of the quadratic equation $[eqn]$ with

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

which describes the condition when susceptible prey, infected prey, and predator exist in the ecosystem.

The existence conditions for each equilibrium point are given by the following theorem.

Theorem 7. For each equilibrium point $[eqn]$ , we have

(i) The PDPF $[eqn]$ is always exists;(ii) The PFP $[eqn]$ exists only if $[eqn]$ ; and(i) The CEP $[eqn]$ exists and unique if n>ad and: (i) $[eqn]$ , or (ii) $[eqn]$ and $[eqn]$ .

Proof: It is clear that the PDPF always exists since $[eqn]$ . For PFP, we emphasize that $[eqn]$ only if $[eqn]$ . Next, the CEP $[eqn]$ only when $[eqn]$ , $[eqn]$ , and $[eqn]$ . We identify that $[eqn]$ and $[eqn]$ are always positive numbers. It is clear that $[eqn]$ only when n>ad. Moreover, $[eqn]$ when $[eqn]$ . Thus, the existence of CEP depends on the value of $[eqn]$ . By identifying the positive root of P, we obtain

$$[eqn]$$ $$[eqn]$$

Therefore, three conditions occur: (i) when $[eqn]$ , $[eqn]$ and CEP does not exist, (ii) when $[eqn]$ , both P1 and P2 become complex numbers and CEP does not exist, (iii) when $[eqn]$ , P1 always negative and P2<0 only when $[eqn]$ . When $[eqn]$ , the value of $[eqn]$ always negative. When $[eqn]$ , $[eqn]$ only when $[eqn]$ . This shows that only one equilibrium point may exist in the interior of the model (6). □

4.3 Local stability analysis

We study the dynamical behavior around the equilibrium point of model (6). The local stability theorem given by Matignon’s (see Theorem 2) is employed. Thus the following theorems are successfully established.

Theorem 8. The PDPF is LAS if $[eqn]$ and a saddle point if $[eqn]$ .

Proof: Applying linearization, we have the following Jacobian matrix for $[eqn]$ .

$$[eqn]$$

Therefore, the following eigenvalues hold: $[eqn]$ , $[eqn]$ , and $[eqn]$ . Since $[eqn]$ , we confirm $[eqn]$ and hence the stability of $[eqn]$ depends on the sign of $[eqn]$ . It is verified that $[eqn]$ when $[eqn]$ and $[eqn]$ when $[eqn]$ . Following Matignon’s condition (see Theorem 2), the validity of the theorem has been proven. □

Theorem 9. Let $[eqn]$ , $[eqn]$ , and $[eqn]$ . The PFP is LAS if $[eqn]$ and a saddle point if $[eqn]$ .

Proof: The Jacobian matrix of model (6) at $[eqn]$ is

$$[eqn]$$

which gives eigenvalues: $[eqn]$ and $[eqn]$ . Since $[eqn]$ for any $[eqn]$ , the local stability properties depend on the argument of $[eqn]$ . We compute that $[eqn]$ when $[eqn]$ and $[eqn]$ when $[eqn]$ . Following Theorem 2, the LAS and saddle point properties given by Thoerem 9 are proven. □

Theorem 10. Let

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The CEP $[eqn]$ is LAS if the following statements hold.

(i) $[eqn]$ , c1>0*, c3>0, and* $[eqn]$ , or(ii) $[eqn]$ , $[eqn]$ , $[eqn]$ , c3>0 and $[eqn]$ , or(iii) $[eqn]$ , c1>0*, c2>0, and* $[eqn]$ .

If $[eqn]$ , c _ 1 _ <0, c _ 2 _ <0, and $[eqn]$ then CEP is unstable.

Proof: We obtain the Jacobian matrix at the CEP as follows.

$$[eqn]$$

Thus, the polynomial characteristic of the Jacobian matrix (16) is $[eqn]$ . Using Generalized Routh-Hurwitz for Caputo fractional-order derivative as in Theorem 3, all statements given by Theorem 10 are well proven. □

4.4 Global stability analysis

We provide the global stability analysis of model (6) by employing the linear, quadratic, and Volterra function to construct the Lyapunov function as well as applying the generalized LaSalle invariance principle given by Lemmas 3 to 5. As a result, the following theorem is presented.

Theorem 11. Let $[eqn]$ and hence $[eqn]$ . The PDPF is GAS if $[eqn]$ and $[eqn]$ .

Proof: The model (6) can be written as

$$[eqn]$$

Now, we define a positive definite quadratic-linear Lyapunov function as follows.

$$[eqn]$$

Using Lemma 3, the following inequality can be computed.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Using the fact $[eqn]$ , we confirm that $[eqn]$ . It also can be ensured that $[eqn]$ only when $[eqn]$ . Obeying Lemma 5, the GAS properties of PDPF are justified. □

Theorem 12. Let $[eqn]$ , $[eqn]$ , and hence $[eqn]$ . The PFP is GAS if $[eqn]$

Proof: We rewrite model (6) into

$$[eqn]$$

We establish the positive definite Volterra-quadratic-linear Lyapunov function as follows.

$$[eqn]$$

Using Lemmas 3 and 4, we obtain the following inequality.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Since $[eqn]$ , we have $[eqn]$ . We also have $[eqn]$ only when $[eqn]$ . Following Lemma 5, we confirm the PFP is GAS. □

Theorem 13. Let

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The CEP $[eqn]$ is GAS if $[eqn]$ for $[eqn]$ .

Proof: The model (6) can be rewritten as

$$[eqn]$$

Moreover, we construct a positive definite Volterra-quadratic Lyapunov function as follows.

$$[eqn]$$

which according to Lemmas 3 and 4, and the the upper bound γ given by Theorem 6, it can be derived as follows.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Since $[eqn]$ affect that $[eqn]$ and $[eqn]$ only when $[eqn]$ , we confirm from Lemma 5 that the CEP $[eqn]$ satisfies GAS behavior. □

4.5 The existence of bifurcations

Bifurcations are a remarkable occurrence in mathematical modeling using fractional derivatives. In this subsection, we will show the occurrence of bifurcations which indicate the change of the dynamical behaviors when a parameter is varied. Several dynamics can occur such as the change of the number of equilibrium points, the change of the stability of the equilibrium points, and the occurrence of the limit cycle. The first bifurcation that occurs in this model is given by forward bifurcation. Forward bifurcation is a condition when a single stable (or unstable) equilibrium point changes into an unstable (or stable) equilibrium point accompanied by the occurrence of the new branch of stable (or unstable) equilibrium point [57,58]. The existence of forward bifurcation is given by the following theorem.

Theorem 14. The PDPF undergoes a forward bifurcation when the basic reproduction number exceeds $[eqn]$ .

Proof: It is clear that PDPF always exists. Furthermore, according to Theorems 7 and 8, when $[eqn]$ , the PDPF is LAS and the PFP does not exist. When $[eqn]$ passes through $[eqn]$ , the PDPF becomes a saddle point while PFP is LAS. This ends the proof. □

Now, we identify the existence of the next bifurcation phenomenon called Hopf bifurcation. The occurrence of Hopf bifurcation was marked by the emergence of a limit cycle and the change of stability of an equilibrium point simultaneously driven by a parameter [59,60]. For the model with fractional derivative as the operator, the occurrence of Hopf bifurcation is not only driven by the parameters on the right-hand side of model (6) but also by the order of the derivative $[eqn]$ . Since the local stability of the equilibrium points lies on the argument of the eigenvalues $[eqn]$ of the Jacobian matrix evaluated at the equilibrium point, the Hopf bifurcation driven by α occurs by the following conditions [61,62]:

$[eqn]$ and $[eqn]$ where $[eqn]$ ; $[eqn]$ ; $[eqn]$ .

If α passes through $[eqn]$ , the stability of the equilibrium point change sign along by the occurrence of a stable limit cycle. In this condition, we call α the bifurcation parameter and $[eqn]$ is the bifurcation point.

5 Numerical simulation

We provide some numerical simulations to support our analytical results, including giving new dynamical behaviors given by model (6) that cannot be described analytically. The Caputo Adam-Bashforth-Moulton (also called the generalized predictor-corrector) numerical scheme provided by Diethelm et al. [63] is employed to obtain the numerical solutions of the model. In this section, we identify the dynamical behaviors of model (6) when the fear level, infection rate, and memory index are varied. Since we do not study a specific ecological system that involves a specific species, all parameter values are provided hypothetically which is adjusted to analytical results. We provide them as follows.

Now, we first investigate the dynamics of model (6) by studying the impact of fear. By setting the parameter values as in Table 2, and varying the fear level (k) in the interval [3,3.6], we have the bifurcation diagram and the phase portraits of model (6) around the CEP. In Fig 1(a), we find that when 1<k < k  for $[eqn]$ , the solution in the interior of model (6) converges to the limit cycle while the CEP is unstable. Furthermore, if we increase the value of k passes through the critical point k , the CEP becomes asymptotically stable while the limit-cycle disappears. This circumstance indicates the occurrence of Hopf bifurcation. We confirm that the parameter k as the fear level becomes the bifurcation parameter while k  is the bifurcation point. From Fig 1(a) we also have a set of stable limit cycles for 1<k < k . The diameter of the limit cycle decreases when the value of k rises. To explore more, we portray the phase portrait of this condition at k = 1 and k = 2.5. When k = 1, we have Fig 1(b) which describes the condition when the solution near the CEP moves away from it and converges to the limit-cycle while the solution far from both the limit-cycle and CEP directly converges to the limit-cycle. When k = 2.5, these solutions converge to the CEP, and the limit cycle disappears. From the biological point of view, the given numerical simulations show the possibility of population density for susceptible prey, infected prey, and the predator when the fear level changes. If the fear level reduces, the density of all populations will change periodically with its period increasing. Conversely, if the effect of fear increases, the density of all populations will converge to an equilibrium point namely the co-existence point. All given simulations show that the density of all populations is always maintained in two different ways, one periodically and the other constantly converging to a single point.

Bifurcation diagram and phase portraits of model (6) using parameter values: Λ = 10, μ = 0.5, β = 1.5, ω = 0.5, δ = 0.5, m = 4, a = 1, n = 1, d = 0.5, α = 0.95 and varying the value of k.

Table 2: Parameter set of numerical simulations.

For the next simulation, we investigate the impact of the disease on the existence of the populations. We set the parameter values as in Table 2 as well as varying the value of the infection rate $[eqn]$ . As a result, we have Fig 2(a) as the bifurcation diagram. For $[eqn]$ $[eqn]$ , only an asymptotically stable PDPF $[eqn]$ exists while other does not exist. We set $[eqn]$ and hence $[eqn]$ which is by the analytical result that the PDFP is locally asymptotically stable, see the phase portrait given by Fig 2(b). It shows that the solutions converge to $[eqn]$ . Furthermore, when β crosses $[eqn]$ , the PDPF loses its stability via forward bifurcation and a new branch of an asymptotically stable equilibrium point occurs namely the PFP. The stability of PFP is maintained for interval $[eqn]$ $[eqn]$ . If we set $[eqn]$ , we obtain Fig 2(c) as the phase portrait which shows that the PDPF $[eqn]$ becomes a saddle point while the PFP $[eqn]$ is asymptotically stable. Similarly, the PFP also loses its stability via forward bifurcation and an asymptotically stable CEP occurs simultaneously when β crosses $[eqn]$ . This condition holds for $[eqn]$ $[eqn]$ , For example, we choose $[eqn]$ and hence we have Fig 2(d). The PDPF and PFP become saddle points and the solution converges to the CEP $[eqn]$ . When β passes through $[eqn]$ , the PFP loses its stability via Hopf bifurcation and the solution will converge to a limit-cycle. We set $[eqn]$ and we can see from Fig 2(e) that a limit-cycle occurs as the impact of Hopf bifurcation while the PFP $[eqn]$ unstable. The instability of CEP is maintained for $[eqn]$ $[eqn]$ . Again, a Hopf bifurcation occurs as a result of the increase in infection rate. When β crosses $[eqn]$ , the limit cycle disappears and the CEP becomes an asymptotically stable equilibrium point, see Fig 2(f) for the example of the phase portrait which shows that the solution converges to $[eqn]$ . From those simulations given by Fig 2, we conclude that the infection rate has a global impact on all populations. For the $[eqn]$ , the infected prey and the predator will become extinct, where only the susceptible prey exists in the ecosystem. If the infection rate increases, the infected prey occurs and the predator extinct. When the infection rate increases again, the susceptible prey, infected prey, and predator will always exist in the ecosystem in two ways, directly tend to a co-existence point or change periodically around the co-existence point. It means the disease impacts not only the prey’s existence but also the predator’s existence.

Phase portraits of model (6) using parameter values: Λ = 10, k = 0.4, μ = 0.5, ω = 0.5, δ = 0.5, m = 4, a = 1, n = 1, d = 0.5, α = 0.95, and varying the value of β.

The final numerical simulations are done to show the dynamical behaviors of model (6) by considering the index of the memory which is stated by the order of the derivative $[eqn]$ . The Table 2 are used as the parameter values. By varying the order in the interval $[eqn]$ , we obtain the bifurcation diagram given by Fig 3(a). It shows that when $[eqn]$ , the CEP $[eqn]$ is asymptotically stable. When α crosses $[eqn]$ , a set of limit cycles occurs until $[eqn]$ and the CEP becomes unstable. The diameter of the limit cycle also increases when α rises. This phenomenon is called Hopf bifurcation where α is the bifurcation parameter and $[eqn]$ is the bifurcation point. Biologically, this circumstance indicates that the memory index affects the dynamics of the population. For a strong memory, the susceptible prey, infected prey, and the predator will directly converge to a constant population point and successfully maintain their existence. For a weak memory, the all population also gets their existence but in different ways. Each of them will eventually change periodically without losing their existence.

Phase portraits of model (6) using parameter values: λ = 10, k = 0.4, μ = 0.5, β = 1.5, ω = 0.5, δ = 0.5, m = 4, a = 1, n = 1, d = 0.5, and varying the value of α.

6 Conclusions

The dynamical behaviors of model (6) which represent an interaction between prey and its predator have been studied. Some assumptions are provided to include more realistic phenomena in the model such as the fear factor, a disease in prey, and the memory effect. The mathematical validity has been given by showing the existence, uniqueness, non-negativity, and boundedness of the model. The dynamics of the model have been studied both analytically and numerically. We have found three equilibrium points along with their dynamical behaviours both locally and globally by applying the Matignon condition, Lyapunov function, and LaSalle invariance principle. The occurrence of Hopf bifurcation driven by the memory index has been proven analytically. The dynamical behaviors have been finally explored numerically by showing the occurrence of forward and Hopf bifurcation driven by the fear level, infection rate, and memory effect using the bifurcation diagrams under the generalized Adam-Bashforth-Moulton numerical scheme. To support all results, some phase portraits are given for each interval in bifurcation diagrams which describe the change of the dynamical behaviors when the bifurcation parameter passes through the bifurcation point. From the biological point of view, some scenarios may occur such as the extinction of both infected prey and predator, the extinction of predator only, or the existence of all populations by maintaining their population densities in two ways, directly converging to a single point or will eventually change periodically. Although the model is completelly analyzed both analytically and numerically, we can develop the model by considering the real phenomena given by a specific species by considering their unique biological components, for example the way the predator hunting the prey, the prey defence, the climate change, and so forth. These conditions will give a significant impact to the mathematical model as well as the results.