Open AccessArticle

Modeling and Dynamical Analysis of a Fractional-Order Predator–Prey System with Anti-Predator Behavior and a Holling Type IV Functional Response

Baiming Wang

and

Xianyi Li

Department of Big Data Science, School of Science, Zhejiang University of Science and Technology, Hangzhou 310023, China

Author to whom correspondence should be addressed.

Fractal Fract. 2023, 7(10), 722; https://doi.org/10.3390/fractalfract7100722

Submission received: 13 August 2023 / Revised: 25 September 2023 / Accepted: 26 September 2023 / Published: 30 September 2023

(This article belongs to the Special Issue Fractional Differential Equations: Computation and Modelling with Applications)

Download

Browse Figures

Versions Notes

Abstract

We here investigate the dynamic behavior of continuous and discrete versions of a fractional-order predator–prey system with anti-predator behavior and a Holling type IV functional response. First, we establish the non-negativity, existence, uniqueness and boundedness of solutions to the system from a mathematical analysis perspective. Then, we analyze the stability of its equilibrium points and the possibility of bifurcations using stability analysis methods and bifurcation theory, demonstrating that, under specific parameter conditions, the continuous system exhibits a Hopf bifurcation, while the discrete version exhibits a Neimark–Sacker bifurcation and a period-doubling bifurcation. After providing numerical simulations to illustrate the theoretically derived conclusions and by summarizing the various analytical results obtained, we finally present four interesting conclusions that can contribute to better management and preservation of ecological systems.

Keywords:

Caputo fractional derivative; memory effect; supercritical Hopf bifurcation; Neimark–Sacker bifurcation; period-doubling bifurcation

MSC:

39A28; 39A30

1. Introduction

The study of predator–prey systems can be traced back to the 18th century. However, the establishment of predator–prey systems in the modern sense is primarily attributed to the work of Alfred J. Lotka and Vito Volterra in the early 20th century. They independently proposed models for predator–prey systems and conducted in-depth research on their dynamical behavior. Their research laid an important theoretical foundation for the dynamics of predator–prey systems. Lotka–Volterra systems describe the interactions between predator and prey, and quickly became a hot topic in dynamical research. Even today, studying the dynamical relationship between predator and prey remains an important subject. To comprehend the intricate dynamical properties presented in predator–prey systems, numerous researchers have dedicated themselves to studying predator–prey models in depth. During the course of their research, they have found multitudes of fascinating dynamical properties among various systems. Berryman [1] pointed out that the original Lotka–Volterra predator–prey equations were built around the principle of mass action, and equations derived in this way lead to the paradoxes of enrichment and biological control. Wang and Chen [2] established the condition for the permanence of populations and sufficient conditions under which positive equilibrium of the model is globally stable. References [3,4] investigated the complex dynamical behavior of discrete predator–prey systems. Qi and Meng [5] found that in a predator–prey system with prey refuge and the fear effect, the survival rate of prey can be improved by increasing the strength of the refuge, decreasing the cost of fear or reducing the intensity of white noise. Blasius, Rudolf, et al. [6] found through experiments that a long-term cyclic persistence exists in a simple predator–prey model. Mukherjee [7] and Qiu and Guo [8] investigated the complex dynamics of a predator–prey system with the fear effect and a predator–prey system with prey taxis, respectively. Although we may subconsciously assume that prey are inherently disadvantaged in a predator–prey system, there are many instances where prey can resist predation and cause harm to predators, even leading to the death of the predator. The existence of such scenarios underscores the significance of determining which entities hold the advantage of prey in a predator–prey system [9].

The anti-predator behavior of prey is widely observed in the natural world. Many scholars have conducted research on the anti-predator behavior of prey and have identified two main ways in which prey exhibit such behaviors: (1) through morphological or behavioral changes [10,11], or (2) by actively attacking the predator [12,13,14].

In 1987, Ives and Dobson [15] proposed the following system to simulate anti-predator behavior (1):

\{\begin{matrix} \frac{d x}{d t} = α x (1 - \frac{x}{k}) - β - \frac{e^{- γ β} q x y}{1 + a x}, \\ \frac{d y}{d t} = \frac{b e^{- γ β} q x y}{1 + a x} - d y, \end{matrix}

(1.1)

where the meanings of all parameters are presented in Table 1.

The prey requirements in anti-predator behavior (2) are higher, as they not only demand that adult prey can resist predation by predators but also require adult prey to have the ability to kill the juveniles of the predator. However, there have been few studies on anti-predator behavior (2). In 2015, Tang and Xiao [16] proposed a system to simulate anti-predator behavior (2), and considered the following Holling type IV functional response function system:

\{\begin{matrix} \frac{d u}{d t} = a u (1 - \frac{u}{k}) - \frac{b u v}{h + u^{2}}, \\ \frac{d v}{d t} = \frac{c b u v}{h + u^{2}} - d v - g u v, \end{matrix}

(1.2)

where the meanings of all parameters are given in Table 2. It is worth noting that the Holling type IV functional response function used in this paper was originally proposed by W. Sokol [17] and has been widely applied in predator–prey systems for several decades. It primarily describes a nonlinear interaction between a predator and its prey. The predation rate of the predator adjusts to changes in the population density of prey, exhibiting a saturation tendency.

The concept of fractional derivatives can be traced back to the 18th century, and the mathematician who first proposed fractional derivatives was Liouville [18]. In the 20th century, the mathematician Riesz made the initial reference to the concept of fractional derivatives and conducted research on their properties in reference [19], combining the studies of Liouville and Riesz to establish the Riesz–Liouville definition of fractional derivatives that is used today. Subsequently, the mathematician Caputo introduced the Caputo definition of fractional derivatives in reference [20].

Definition 1.

Denote

_{0}^{C} D_{t}^{α} f (t) = J^{l - α} f^{(l)} (t), α > 0,

where

f^{(l)}

denotes the derivative of f with order l, l is the nearest integer value of α, and

J^{q}

is the operator of the Riemann–Liouville integral of q order:

J^{q} h (t) = \frac{\int_{0}^{t} {(t - τ_{e})}^{q - 1} h (τ_{e}) d τ_{e}}{Γ (q)},

where

Γ (q)

is Euler’s Gamma function. The

a l p h a

-order Caputo differential operator is the term used to describe the operator

_{0}^{C} D_{t}^{α}

From a biological perspective, considering a fractional-order predator–prey system makes logical sense; for most organisms in the natural world, their current behaviors are influenced by historical context. In fractional calculus, the rate of change at any given moment, i.e., the fractional-order derivative, depends on the population density over a certain period of time. Therefore, fractional-order predator–prey systems have unique advantages in describing memory effects within populations. Over the past two decades, owing to the advantages of fractional derivatives in studying various ecological systems’ memory effects, numerous mathematicians have turned their attention to investigating fractional-order ecological systems, finding many interesting dynamical properties presented in them [21,22,23,24,25,26,27,28]. At present, a relatively comprehensive research framework has been established for mathematical models of integer-order ecosystems, while the study of fractional-order ecosystems is still in its early stages. Hence, the authors of this paper intend to introduce the Caputo fractional derivation to system (1.2) and extend it to a fractional-order ecosystem. We intend to employ the Caputo definition of fractional derivatives to analyze how anti-predator behavior and the Holling type IV functional response function in a fractional-order ecosystem will impact the dynamics of the system. As a result, we introduce the following fractional-order predator–prey system with a Holling type IV functional response and anti-predator behaviors:

\{\begin{matrix} _{0}^{C} D_{t}^{α} u (t) = a u (1 - \frac{u}{k}) - \frac{b u v}{h + u^{2}}, \\ _{0}^{C} D_{t}^{α} v (t) = \frac{c b u v}{h + u^{2}} - d v - g u v, \end{matrix}

(1.3)

where the meanings of all parameters are presented in Table 2. For the method of introducing the Caputo fractional differential equation into an ecosystem model, reference [29] may be consulted.

There is a current lack of comprehensive dynamic analysis methods for continuous fractional-order predator–prey systems. For example, in the literature [30], the analysis of fractional-order systems has mainly focused on Hopf bifurcations, while more extensive research has been dedicated to analyzing integer-order cases.

In references [31,32], the authors, respectively, conducted dynamical analyses of different discretized predator–prey models and found that discrete predator–prey models, in comparison to their continuous counterparts, exhibit a greater variety of dynamical behaviors and hold advantages in numerical simulations. In reference [33], the authors employed the piecewise constant approximation (PCA) method to discretize a continuous fractional-order predator–prey system and analyzed the dynamical properties of and discussed the types of bifurcations present in this system. Their work motivates us to consider the discrete counterpart of system (1.3). In recent years, many researchers have studied the dynamical behavior of discrete fractional-order predator–prey systems and have discovered numerous intriguing dynamical properties within these systems [34,35,36].

Hence, in order to better understand the properties of system (1.3), we here also consider discretizing system (1.3) for further dynamical analyses and comparing its properties with the continuous model (1.3), finding that there are many differences in dynamical properties between system (1.3) and its discrete version (1.6). This sufficiently shows that it is very helpful to consider the same problems from different angles.

We use the PCA method to discretize model (1.3), and the specific steps are as follows:

Assume that the initial conditions of system (1.3) are

u (0) = u_{0}

and

v (0) = v_{0}

. The discretized version of model (1.3) is given as

\{\begin{matrix} _{0}^{C} D_{t}^{α} u (t) = a u ([\frac{t}{ρ}]) (1 - \frac{u ([\frac{t}{ρ}])}{k}) - \frac{b u ([\frac{t}{ρ}]) v ([\frac{t}{ρ}])}{h + u {([\frac{t}{ρ}])}^{2}}, \\ _{0}^{C} D_{t}^{α} v (t) = \frac{c b u ([\frac{t}{ρ}]) v ([\frac{t}{ρ}])}{h + u {([\frac{t}{ρ}])}^{2}} - d v ([\frac{t}{ρ}]) - g u ([\frac{t}{ρ}]) v ([\frac{t}{ρ}]) . \end{matrix}

First, let

t \in [0, ρ)

, then

\frac{t}{ρ} \in [0, 1)

. Thus, we obtain

\{\begin{matrix} _{0}^{C} D_{t}^{α} u (t) = a u_{0} (1 - \frac{u_{0}}{k}) - \frac{b u_{0} v_{0}}{h + u_{0}^{2}}, \\ _{0}^{C} D_{t}^{α} v (t) = \frac{c b u_{0} v_{0}}{h + u_{0}^{2}} - d v_{0} - g u_{0} v_{0} . \end{matrix}

(1.4)

The answer to (1.4) is simplified to

\begin{matrix} u_{1} (t) & = u_{0} + J^{α} (a u_{0} (1 - \frac{u_{0}}{k}) - \frac{b u_{0} v_{0}}{h + u_{0}^{2}}) \\ = u_{0} + \frac{t^{α}}{α Γ (α)} (a u_{0} (1 - \frac{u_{0}}{k}) - \frac{b u_{0} v_{0}}{h + u_{0}^{2}}), \\ v_{1} (t) & = v_{0} + J^{α} (\frac{c b u_{0} v_{0}}{h + u_{0}^{2}} - d v_{0} - g u_{0} v_{0}) \\ = v_{0} + \frac{t^{α}}{α Γ (α)} (\frac{c b u_{0} v_{0}}{h + u_{0}^{2}} - d v_{0} - g u_{0} v_{0}) . \end{matrix}

Second, let

t \in [ρ, 2 ρ)

, so

\frac{t}{ρ} \in [1, 2)

. Then,

\{\begin{matrix} _{0}^{C} D_{t}^{α} u (t) = a u_{1} (1 - \frac{u_{1}}{k}) - \frac{b u_{1} v_{1}}{h + u_{1}^{2}}, \\ _{0}^{C} D_{t}^{α} v (t) = \frac{c b u_{1} v_{1}}{h + u_{1}^{2}} - d v_{1} - g u_{1} v_{1} . \end{matrix}

(1.5)

After simplifying (1.5), we can obtain the following solution

\begin{matrix} u_{1} (t) & = u_{1} (ρ) + J_{ρ}^{α} (a u_{1} (1 - \frac{u_{1}}{k}) - \frac{b u_{1} v_{1}}{h + u_{1}^{2}}) \\ = u_{1} (ρ) + \frac{{(t - ρ)}^{α}}{α Γ (α)} (a u_{1} (1 - \frac{u_{1}}{k}) - \frac{b u_{1} v_{1}}{h + u_{1}^{2}}), \\ v_{1} (t) & = v_{1} (ρ) + J_{ρ}^{α} (\frac{c b u_{1} v_{1}}{h + u_{1}^{2}} - d v_{1} - g u_{1} v_{1}) \\ = v_{1} (ρ) + \frac{{(t - ρ)}^{α}}{α Γ (α)} (\frac{c b u_{1} v_{1}}{h + u_{1}^{2}} - d v_{1} - g u_{1} v_{1}), \end{matrix}

where

J_{ρ}^{α} = \frac{1}{Γ (α)} \int_{ρ}^{t} {(t - τ_{e})}^{α - 1} d τ_{e}

0 < α < 1

. After n repetitions, we obtain

\begin{matrix} u_{n + 1} (t) & = u_{n} (n ρ) + \frac{{(t - n ρ)}^{α}}{α Γ (α)} (a u_{n} (n ρ) (1 - \frac{u_{n} (n ρ)}{k}) - \frac{b u_{n} (n ρ) v_{n} (n ρ)}{h + u_{n} {(n ρ)}^{2}}), \\ v_{n + 1} (t) & = v_{n} (n ρ) + \frac{{(t - n ρ)}^{α}}{α Γ (α)} (\frac{c b u_{n} (n ρ) v_{n} (n ρ)}{h + u_{n} {(n ρ)}^{2}} - d v_{n} (n ρ) - g u_{n} (n ρ) v_{n} (n ρ)), \end{matrix}

where

t \in [n ρ, (n + 1) ρ)

. For

t \to (n + 1) ρ,

the system above becomes

\begin{matrix} \{\begin{matrix} u_{n + 1} = u_{n} + \frac{{(ρ)}^{α}}{Γ (α + 1)} (a u_{n} (1 - \frac{u_{n}}{k}) - \frac{b u_{n} v_{n}}{h + u_{n}^{2}}), \\ v_{n + 1} = v_{n} + \frac{{(ρ)}^{α}}{Γ (α + 1)} (\frac{c b u_{n} v_{n}}{h + u_{n}^{2}} - d v_{n} - g u_{n} v_{n}) . \end{matrix} \end{matrix}

(1.6)

The overall structure of this paper is described as follows: In Section 2, some preliminaries are provided for some definitions, lemmas and theorems that will be used to analyze the dynamical properties of systems (1.3) and (1.6). In Section 3, the well-posedness of system (1.3) is analyzed. In Section 4, the existence and stability of the equilibrium points of systems (1.3) and (1.6) are investigated, respectively. In Section 5, we demonstrate that, under certain parameter conditions, system (1.3) exhibits a Hopf bifurcation, while system (1.6) exhibits a Neimark–Sacker bifurcation and a period-doubling bifurcation. In Section 6, numerical simulations are performed to validate the results of our theoretical analysis. In Section 7, interesting conclusions are drawn based on some findings in the previous sections.

2. Preliminaries

In this section, we primarily introduce the definition and some conclusions of Caputo fractional derivatives that are necessary for our subsequent research.

Definition 2

([37]). Under the definition of Caputo fractional derivatives, the fractional derivative of function

f (ξ) \in A C^{n} ([0, + \infty], R)

is given as

_{0}^{C} D_{ξ}^{α} f (ξ) = \int_{0}^{ξ} \frac{f^{(n)} (ϑ)}{Γ (n - α) {(ξ - ϑ)}^{α - n + 1}} d ϑ,

where α represents the order of the fractional derivative.

When

n = 1

, the fractional derivative

_{0}^{C} D_{ξ}^{α} f (ξ)

takes the form of

_{0}^{C} D_{ξ}^{α} f (ξ) = \int_{0}^{ξ} \frac{f (ϑ)}{Γ (1 - α) {(ξ - ϑ)}^{α}} d ϑ .

Definition 3

([37]). The Mittag–Leffler function

M_{i}

, when the order i of

M_{i}

is positive, is defined as

M_{i} (ζ) = \sum_{j = 0}^{\infty} \frac{ζ_{j}}{Γ (j i + 1)}, ζ_{j} \in C,

as the sequence converges.

Definition 4

([38]). Let

Q (u, v)

be a fixed piont of system (1.6) with multipliers

λ_{1}

and

λ_{2}

(i)

| λ_{1} | < 1

and

| λ_{2} | < 1

, the fixed point

Q (u, v)

is called a sink, and the sink is locally asymptotically stable.

(i i)

| λ_{1} | > 1

and

| λ_{2} | > 1

, the fixed point

Q (u, v)

is called a source, and the source is locally asymptotically unstable.

(i i i)

| λ_{1} | < 1

and

| λ_{2} | > 1

(or

| λ_{1} | > 1

and

| λ_{2} | < 1

), the fixed point

Q (u, v)

is called a saddle.

(i v)

If either

| λ_{1} | = 1

| λ_{2} | = 1

, the fixed point

Q (u, v)

is called non-hyperbolic.

Lemma 1

([39]). For

_{0}^{C} D_{ξ}^{α} f (ξ) \in A C ([0, Ξ], R)

, if

f (ξ_{1}) = 0

and

f (ξ_{0}) > 0

(all

ξ_{0} \in [0, ξ_{1}))

, then

_{0}^{C} D_{ξ}^{α} f (ξ_{1}) < 0

Lemma 2

([40]). For the fractional-order system

_{0}^{C} D_{ξ}^{α} Y (ξ) = Z (ξ, Y), ξ \geq 0,

with initial condition

Y (0) = (u (0), v (0))

, where

0 < α \leq 1

Z : [0, + \infty) \times τ \to R^{n}

τ \subseteq R

, if

Z (ξ, Y)

fulfills the local Lipschitz condition for

Z \in R_{n}

∥ Z (ξ, Y) - Z (ξ, \tilde{Y}) ∥ \leq Δ \cdot ∥ Y - \tilde{Y} ∥,

then the system has a unique solution on

[0, + \infty) \times τ

, and

∥ Y (y_{1}, y_{2}, y_{3}, \dots, y_{n}) - Y (\tilde{y_{1}}, \tilde{y_{2}}, \tilde{y_{3}}, \dots, \tilde{y_{n}}) ∥ \leq \sum_{i = 1}^{n} | y_{i} - \tilde{y_{i}} |,

for

i = 1, 2, 3, \dots, n, y_{i}, \tilde{y_{i}} \in R .

Lemma 3

([41]). Let

F (λ) = λ^{2} + B λ + C

, where B and C are two real constants. Suppose

λ_{1}

and

λ_{2}

are two roots of

F (λ) = 0

. Then, the following statements hold.

(i)

F (1) > 0,

then

(i . 1)

| λ_{1} | < 1

and

| λ_{2} | < 1

if and only if

F (- 1) > 0

and

C < 1

;

(i . 2)

λ_{1} = - 1

and

λ_{2} \neq - 1

if and only if

F (- 1) = 0

and

B \neq 2

;

(i . 3)

| λ_{1} | < 1

and

| λ_{2} | > 1

if and only if

F (- 1) < 0

;

(i . 4)

| λ_{1} | > 1

and

| λ_{2} | > 1

if and only if

F (- 1) > 0

and

C > 1

;

(i . 5)

λ_{1}

and

λ_{2}

are a pair of conjugate complex roots and

| λ_{1} | = | λ_{2} | = 1

if and only if

- 2 < B < 2

and

C = 1

;

(i . 6)

λ_{1} = λ_{2} = - 1

if and only if

F (- 1) = 0

and

B = 2

(i i)

F (1) = 0,

namely, 1 is a root of

F (λ) = 0

, then the another root λ satisfies

| λ | = (<, >) 1

if and only if

| C | = (<, >) 1 .

(i i i)

F (1) < 0,

then

F (λ) = 0

has one root lying in

(1, \infty)

. Moreover,

(i i i . 1)

The other root λ satisfies

λ < (=) - 1

if and only if

F (- 1) < (=) 0

;

(i i i . 2)

The other root

- 1 < λ < 1

if and only if

F (- 1) > 0

Theorem 1

([42]). The Laplace transform of

_{0}^{C} D_{ξ}^{α} f (ξ)

L [_{0}^{C} D_{ξ}^{α} f (ξ)] = ϑ^{α} F (ϑ) - \sum_{j = 0}^{n - 1} ϑ^{α - j - 1} f^{j} (0),

where

F (ϑ) = L [f (ξ)]

α \in (n - 1, n), n \in Z_{+}

Theorem 2

([43]). Assume

α > 0

β > 0

and

K \in C^{n \times n}

, then

L [ξ^{β - 1} E_{α, β} (K ξ^{α})] = \frac{ϑ^{α - β}}{ϑ^{α} - K},

for

R e (ϑ) > {∥ K ∥}^{\frac{1}{α}}

, where

R e (ϑ)

is the real part of complex number ϑ and

E_{α, β}

is the Mittag–Leffler function.

Theorem 3

([44]). For the following fractional-order system

_{0}^{C} D_{ξ}^{α} f (ξ) = g (f (ξ)), f (0) = f_{0} \in R^{N}, α \in (0, 1),

where

f (ξ) = {(f_{1} (ξ), f_{2} (ξ), f_{3} (ξ), \dots, f_{n} (ξ))}^{T} \in R^{n}

and

g = {(g_{1}, g_{2}, g_{3}, \dots, g_{n})}^{T} : R^{n} \to R^{n}

. If

g (f^{*}) = 0

, then

f^{*}

is an equilibrium point. Set

J (f^{*})

as the Jacobian matrix

J = \frac{\partial g}{\partial f} = \frac{\partial (g_{1}, g_{2}, g_{3}, \dots, g_{n})}{\partial (f_{1}, f_{2}, f_{3}, \dots, f_{n})}

for

f = f^{*}

. If the characteristic values

λ_{i} (i = 1, 2, 3, \dots, n)

J (f^{*})

meet

| a r g (λ_{i}) | > \frac{α π}{2} (i = 1, 2, 3, \dots, n)

, then

f^{*}

is locally asymptotically stable.

Theorem 4

([45]). We say that a fractional-order system undergoes a fractional Hopf bifurcation if there exists a critical value

β = β_{c}

such that the following conditions are satisfied:

1.: $λ_{1} (β_{c})$ and $λ_{2} (β_{c})$ satisfy $| a r g (λ_{i} (β_{c})) | = \frac{π α}{2}, (i = 1, 2)$ ;
2.: $| a r g (λ_{i} (β_{c})) | \neq \frac{π α}{2}, (i = 3, 4, 5, \dots, n)$ ;
3.: $\frac{d}{d β} | a r g (λ_{i} (β)) {| |}_{β = β_{c}} \neq 0, (i = 1, 2)$ ,

where λ represents the eigenvalues of the Jacobian matrix of the system.

3. Analysis of the Well-Posedness of System (1.3)

In this section, we consider the uniqueness, non-negativity and boundedness of the solutions of system (1.3).

Theorem 5.

For the initial condition

(u (0), v (0)) \in A,

system (1.3) has a unique solution

U (t) = (u (t), v (t)) \in A

for all

t \geq 0

, where

A = {(u, v) \in R^{2} : m a x {| u |, | v s . |} < γ_{1},

m i n {| u |, | v s . |} > γ_{2}}

Proof.

Consider the time interval

[0, t_{1}], t_{1} < + \infty

. Construct a mapping

G (U) = (G_{1} (U),

G_{2} {(U))}^{T}

, where

U = {(u, v)}^{T}

and

\{\begin{matrix} G_{1} (U) = a u (1 - \frac{u}{k}) - \frac{b u v}{h + u^{2}}, \\ G_{2} (U) = \frac{c b u v}{h + u^{2}} - d v - g u v . \end{matrix}

For

U, \tilde{U} \in A

, we have

\begin{matrix} ∥ G (U) - G (\tilde{U}) ∥ \\ = | G_{1} (U) - G_{1} (\tilde{U}) | + | G_{2} (U) - G_{2} (\tilde{U}) | \\ = | a u (1 - \frac{u}{k}) - \frac{b u v}{h + u^{2}} - a \tilde{u} (1 - \frac{\tilde{u}}{k}) + \frac{b \tilde{u} \tilde{v}}{h + {\tilde{u}}^{2}} | \end{matrix}

\begin{matrix} + | \frac{c b u v}{h + u^{2}} - d v - g u v - \frac{c b \tilde{u} \tilde{v}}{h + {\tilde{u}}^{2}} - d \tilde{v} - g \tilde{u} \tilde{v} | \\ \leq a | u - \tilde{u} | + \frac{a}{k} | u + \tilde{u} | | u - \tilde{u} | + | \frac{b u v (h + {\tilde{u}}^{2}) - b \tilde{u} \tilde{v} (h + u^{2})}{(h + u^{2}) (h + {\tilde{u}}^{2})} | \\ + | \frac{c b u v (h + {\tilde{u}}^{2}) - c b \tilde{u} \tilde{v} (h + u^{2})}{(h + u^{2}) (h + {\tilde{u}}^{2})} | + d | v - \tilde{v} | + g | u v - \tilde{u} \tilde{v} | \end{matrix}

\begin{matrix} = a | u - \tilde{u} | + \frac{a}{k} | u + \tilde{u} | | u - \tilde{u} | \\ + (1 + c) | \frac{b h u (v - \tilde{v}) + b h \tilde{v} (u - \tilde{u})}{(h + u^{2}) (h + {\tilde{u}}^{2})} + \frac{b u^{2} \tilde{u} (v - \tilde{v}) - b u \tilde{u} v (u - \tilde{u})}{(h + u^{2}) (h + {\tilde{u}}^{2})} | \\ + d | v - \tilde{v} | + g | u (v - \tilde{v}) + \tilde{v} (u - \tilde{u}) | \\ \leq (a + g γ_{1} + \frac{2 a γ_{1}}{k} + \frac{b (1 + c) (h γ_{2} + γ_{2}^{3})}{{(h + γ_{2}^{2})}^{2}}) | u - \tilde{u} | \\ + (d + g γ_{1} + \frac{b (1 + c) (h γ_{2} + γ_{2}^{3})}{{(h + γ_{2}^{2})}^{2}}) | v - \tilde{v} | \\ = L_{1} | u - \tilde{u} | + L_{2} | v - \tilde{v} | \leq L ∥ U - \tilde{U} ∥, \end{matrix}

where

L = m a x {L_{1}, L_{2}}

with

L_{1} = a + g γ_{1} + \frac{2 a γ_{1}}{k} + \frac{b (1 + c) (h γ_{2} + γ_{2}^{3})}{{(h + γ_{2}^{2})}^{2}}

and

L_{2} = d + g γ_{1} + \frac{b (1 + c) (h γ_{2} + γ_{2}^{3})}{{(h + γ_{2}^{2})}^{2}}

. □

Hence,

G (U)

conforms to the local Lipschitz condition, and system (1.3) has a unique solution by Lemma 2.

Theorem 6.

All solutions of system (1.3) initiating from

(u (0), v (0)) \in R^{+}

are non-negative and bounded in the region W, where

W = \{(u (t), v (t) \in R_{+}^{2}) ∥ 0 \leq u (t) + \frac{v (t)}{c} \leq \frac{k {(a + d)}^{2}}{4 a d}\} .

Proof.

First, let us prove the non-negativity of the solution. Assume that there exists a constant

μ (> 0)

that satisfies

\{\begin{matrix} u (t) > 0, t \in [0, μ), \\ u (μ) = 0, \\ u (μ^{+}) < 0 . \end{matrix}

We can easily find that

_{0}^{C} D_{t}^{α} {u (t) |}_{μ} = 0

, and derive

u (μ^{+}) = 0

from Lemma 1, which obviously contradicts

u (μ^{+}) < 0

. Thus,

u (t) > 0

for any

t \in [0, + \infty)

. Similarly, we can prove

v (t) > 0

for

\forall t \in [0, + \infty)

Now, construct a function

X (t) = u (t) + \frac{v (t)}{c}

, which will help us prove the boundedness of the solution. The Caputo fractional derivative of

X (t)

with order

α

\begin{matrix} _{0}^{C} D_{t}^{α} X (t) & = a u (t) - \frac{a u {(t)}^{2}}{k} - \frac{d v (t)}{c} - \frac{g u v}{c} \\ = a u (t) (1 - \frac{u (t)}{k}) - \frac{g u (t) v (t)}{c} - d (X (t) - u (t)) . \end{matrix}

Then,

\begin{matrix} _{0}^{C} D_{t}^{α} X (t) + d X (t) & = a u (t) (1 - \frac{u (t)}{k}) - \frac{g u (t) v (t)}{c} + d u (t) \\ \leq \frac{k {(a + d)}^{2}}{4 a} - {(u (t) \sqrt{\frac{a}{k}} - \frac{(a + d)}{2} \sqrt{\frac{k}{a}})}^{2} \\ \leq \frac{k {(a + d)}^{2}}{4 a}, \end{matrix}

i.e.,

\begin{matrix} _{0}^{C} D_{t}^{α} X (t) + d X (t) \leq \frac{k {(a + d)}^{2}}{4 a} . \end{matrix}

Applying Theorem 1 and taking the Laplace transform on both sides of the above inequality at the same time, one has

L [_{0}^{C} D_{t}^{α} X (t) + d X (t)] \leq L [\frac{k {(a + d)}^{2}}{4 a}] .

This leads to

ϑ^{α} F (ϑ) - ϑ^{α - 1} X (0) + d F (ϑ) \leq \frac{1}{ϑ} \frac{k {(a + d)}^{2}}{4 a},

where

F (ϑ) = L [X (t)]

. Hence,

\begin{matrix} F (ϑ) \leq \frac{k {(a + d)}^{2}}{4 ϑ a (ϑ^{α} + d)} + \frac{ϑ^{α - 1}}{ϑ^{α} + d} X (0) . \end{matrix}

By using the inverse Laplace transform on both sides of the above inequality, we may derive

\begin{matrix} L^{- 1} [F (ϑ)] \leq T L^{- 1} [\frac{1}{ϑ (ϑ^{α} + d)}] + X (0) L^{- 1} [\frac{ϑ^{α - 1}}{ϑ^{α} + d}] \\ \Rightarrow & X (t) \leq T L^{- 1} [\frac{ϑ^{- 1}}{ϑ^{α} + d}] + X (0) L^{- 1} [\frac{ϑ^{α - 1}}{ϑ^{α} + d}], \end{matrix}

where

T = \frac{k {(a + d)}^{2}}{4 a}

. From Theorem 2, one obtains

\begin{matrix} X (t) \leq T t^{α} E_{α, α + 1} (- d t^{α}) + X (0) E_{α, 1} (- d t^{α}) . \end{matrix}

According to the properties of the Mittag–Leffer function, we get

\begin{matrix} E_{α, 1} (- d t^{α}) = - d t^{α} E_{α, α + 1} (- d t^{α}) + \frac{1}{Γ (1)}, \end{matrix}

i.e.,

\begin{matrix} - \frac{1}{d} [E_{α, 1} (- d t^{α}) - 1] = t^{α} E_{α, α + 1} (- d t^{α}), \end{matrix}

which displays

\begin{matrix} X (t) \leq (X (0) - \frac{T}{d}) E_{α, 1} (- d t^{α}) + \frac{T}{d} . \end{matrix}

Note that

E_{α, 1} \to 0

when

t \to \infty

. Thus, we have

X (t) \leq \frac{T}{d}

for large t, i.e.,

X (t) \leq \frac{k {(a + d)}^{2}}{4 a d}

for large t. Accordingly, all solutions of system (1.3) are bounded in the region

\begin{matrix} W = \{(u (t), v (t)) \in R_{+}^{2} ∥ 0 \leq u (t) + \frac{v (t)}{c} \leq \frac{k {(a + d)}^{2}}{4 a d}\} . \end{matrix}

The proof is over. □

4. Local Stability of Systems (1.3) and (1.6)

In this section, we first identify the equilibrium points of systems (1.3) and (1.6), which are the same, then analyze their local stability.

4.1. Existence of an Equilibrium Point

We first can easily observe that the two points

Q_{0} (0, 0)

and

Q_{k} (k, 0)

always are equilibrium points of systems (1.3) and (1.6).

Next, we consider the positive equilibrium points of systems (1.3) and (1.6). It is evident that the positive equilibrium solutions of systems (1.3) and (1.6) satisfy the following equations:

\{\begin{matrix} a (1 - \frac{u}{k}) - \frac{b v}{h + u^{2}} = 0, \\ \frac{c b u}{h + u^{2}} - d - g u = 0 . \end{matrix}

(4.1)

By performing a transformation on the second equation, we find that the component u of positive equilibrium point (u,v) meets the following equation

p (u) = g u^{3} + d u^{2} + (g h - c b) u + d h = 0,

(4.2)

while the positive component

v = \frac{a}{b} (1 - \frac{u}{k}) (h + u^{2})

. Therefore, the problem of finding positive equilibrium points of systems (1.3) and (1.6) is transferred to that of solving the positive solutions of Equation (4.2). It is easy to derive

p^{'} (u) = 3 g u^{2} + 2 d u + g h - c b, p^{″} (u) = 6 g u + 2 d .

(4.3)

Obviously,

p^{″} (u) > 0

always holds for

u > 0

. This implies

p^{'} (u)

is monotonically increasing for

u > 0

. Now, consider the solutions of

p (u) = 0

according to the following two cases.

Case 1:

g h - c b \geq 0

. Then,

p^{'} (u) > 0

, indicating that

p (u)

is monotonically increasing. Again,

p (0) = d h > 0

. Therefore, there are no positive solutions of

p (u) = 0

for

g h - c b \geq 0

, which then implies that system (1.3) has no positive equilibrium points.

Case 2:

g h - c b < 0

. Then,

p^{'} (u) = 0

has a unique positive solution, denoted by

u_{*}

, where

u_{*} = \frac{- d + \sqrt{d^{2} - 3 g (g h - c b)}}{3 g}

. Furthermore, since

p^{'} (u)

is monotonically increasing, we can conclude that

p (u)

is monotonically decreasing in the interval

(0, u_{*})

, whereas it is monotonically increasing in the interval

(u_{*}, + \infty)

. Thus, the function

p (u)

has a minimum at

u = u_{*}

for

u \in (0, \infty)

. Substituting

u_{*}

into (4.2), we obtain

\begin{matrix} p (u_{*}) = & \frac{1}{27 g^{2}} [R^{3} + 3 (3 h g^{2} - 3 b c g - d^{2}) R + 3 (d^{2} + 3 b c d g + 6 d h g^{2})] \\ = & \frac{1}{27 g^{2}} [- 2 R^{3} + 3 (d^{2} + 3 b c d g + 6 d h g^{2})] = \frac{2}{27 g^{2}} (R_{0}^{3} - R^{3}), \end{matrix}

where

R = \sqrt{d^{2} - 3 g (g h - c b)}

and

R_{0} = \sqrt[3]{\frac{3}{2} (d^{2} + 3 b c d g + 6 d h g^{2})}

. Then, we can discuss the positive solution of

p (u) = 0

in view of the following three subcases:

Subcase 1.

p (u_{*}) > 0 ⟺ R < R_{0}

. This means that the equation

p (u) = 0

has no positive roots.

Subcase 2.

p (u_{*}) = 0 ⟺ R = R_{0}

. This indicates that there is only one positive solution

u_{*}

of the equation

p (u) = 0

Subcase 3.

p (u_{*}) < 0 ⟺ R > R_{0}

. This means that there are two positive roots of the equation

p (u) = 0

, denoted by

u_{1}

and

u_{2}

. Namely,

\begin{matrix} u_{1} = \frac{- d + \sqrt{R} (c o s (\frac{o}{3}) - \sqrt{3} s i n (\frac{o}{3}))}{3 g}, \\ u_{2} = \frac{- d + \sqrt{R} (c o s (\frac{o}{3}) + \sqrt{3} s i n (\frac{o}{3}))}{3 g}, \end{matrix}

where

o = a r c c o t (J)

J = \frac{2 d R^{2} - 3 g T}{2 R^{3}} (J \in (- 1, 1))

and

T = d (g h - c b) - 9 g h d

. Evidently,

0 < u_{1} < u_{*} < u_{2}

Denote the two positive equilibria as

Q_{i} (u_{i}, v_{i})

u_{i} < k, i = 1, 2

. Summarizing the above analysis, we can obtain the following result.

Theorem 7.

Let

R, R_{0}, u_{*}, u_{1}, u_{2}

be, respectively, defined in Case 2 and Case 3. For the existence of an equilibrium point of systems (1.3) and (1.6), the following statements hold.

1.

Regardless of the value of the parameters, systems (1.3) and (1.6) always have a trivial equilibrium point

Q_{0} (0, 0)

and a boundary equilibrium point

Q_{k} (k, 0)

2.

When

g h - c b \geq 0

, systems (1.3) and (1.6) do not have positive equilibrium points.

3.

When

g h - c b < 0

, we further have the following conclusions.

(a): If $R_{0} > R$ , then systems (1.3) and (1.6) do not have positive equilibrium points.
(b): If $R_{0} = R$ , then, for $0 < k \leq u_{*}$ , systems (1.3) and (1.6) do not have positive equilibrium points; for $u_{*} < k$ , systems (1.3) and (1.6) have one positive equilibrium point $Q_{*} (u_{*}, v_{*})$ .
(c): If $R_{0} < R$ , then, for $0 < k \leq u_{1}$ , systems (1.3) and (1.6) do not have positive equilibrium points; for $u_{1} < k \leq u_{2}$ , systems (1.3) and (1.6) have only one positive equilibrium point $Q_{1} (u_{1}, v_{1})$ ; for $u_{2} < k$ , systems (1.3) and (1.6) have two positive equilibrium points $Q_{1} (u_{1}, v_{1})$ and $Q_{2} (u_{2}, v_{2})$ .

Next, we analyze the stability of these equilibrium points of systems (1.3) and (1.6).

4.2. Stability Analysis of Equilibrium Points of System (1.3)

The Jacobian matrix of system (1.3) at any equilibrium

Q (u, v)

is as follows

J_{c} (u, v) = (\begin{matrix} a (1 - \frac{2 u}{k}) - \frac{b v}{h + u^{2}} + \frac{2 b u^{2} v}{{(h + u^{2})}^{2}} & - \frac{b u}{h + u^{2}} \\ \frac{b c v (h - u^{2})}{{(h + u^{2})}^{2}} - g v & \frac{b c u}{h + u^{2}} - d - g u \end{matrix}) .

4.2.1. The Stability of the Trivial Equilibrium Point $Q_{0} (0, 0)$

Theorem 8.

The trivial equilibrium point

Q_{0} (0, 0)

is a saddle.

Proof.

Substituting the trivial equilibrium point

Q_{0} (0, 0)

into the Jacobian matrix

J_{c} (u, v)

, we obtain

J_{c} (Q_{0}) = (\begin{matrix} a & 0 \\ 0 & - d \end{matrix}),

and it is easy to see that the Jacobian matrix

J_{c} (Q_{0})

has two eigenvalues:

λ_{1} = a > 0

and

λ_{2} = - d < 0

. Since

| a r g (λ_{1}) | = 0 < \frac{α π}{2}

and

| a r g (λ_{2}) | = π > \frac{α π}{2}

, the trivial equilibrium point

Q_{0}

is a saddle. □

4.2.2. The Stability of Boundary Equilibrium Point $Q_{k} (k, 0)$

Theorem 9.

The boundary equilibrium point

Q_{k} (k, 0)

is a stable node for

d > \frac{c b k}{h + k^{2}} - g k

, while it is a saddle for

d < \frac{c b k}{h + k^{2}} - g k

Proof.

Substituting boundary equilibrium point

Q_{k} (k, 0)

into the Jacobian matrix

J_{c} (u, v)

, we have

J_{c} (Q_{k}) = (\begin{matrix} - a & \frac{- b k}{h + k^{2}} \\ 0 & \frac{c b k}{h + k^{2}} - d - g k \end{matrix})

Now, consider the following two cases:

Case 1:

d < \frac{c b k}{h + k^{2}} - g k

. Then, we obtain the two eigenvalues of the Jacobian matrix

J_{c} (Q_{k})

λ_{1} = - a < 0

and

λ_{2} = \frac{c b k}{h + k^{2}} - g k - d > 0

. Therefore,

| a r g (λ_{1}) | = π > \frac{α π}{2}

and

| a r g (λ_{2}) | = 0 < \frac{α π}{2}

. Thus, the boundary equilibrium point

Q_{k}

is a saddle.

Case 2:

d > \frac{c b k}{h + k^{2}} - g k

. Then, the two eigenvalues of the Jacobian matrix

J_{c} (Q_{k})

are

λ_{1} = - a < 0

and

λ_{2} = \frac{c b k}{h + k^{2}} - d - g k < 0

. As a result of

| a r g (λ_{1}) | = π > \frac{α π}{2}

and

| a r g (λ_{2}) | = π > \frac{α π}{2}

, the boundary equilibrium point

Q_{k}

is a stable node.

□

4.2.3. The Stability of Positive Equilibrium Points $Q_{i} (u_{i}, v_{i}) (i = 1, 2)$

Theorem 10.

The positive equilibrium point

Q_{1} (u_{1}, v_{1})

is stable for

k > \frac{h + 3 u_{1}^{2}}{2 u_{1}}

and unstable for

k < \frac{h + 3 u_{1}^{2}}{2 u_{1}}

; the postive equilibrium point

Q_{2} (u_{2}, v_{2})

is always a saddle point.

Proof.

For a better comprehension, let us begin to analyze the stability of the positive equilibrium point

Q_{2} (u_{2}, v_{2})

Substituting the equilibrium point

Q_{2} (u_{2}, v_{2})

into the Jacobian matrix

J_{c} (u, v)

, one obtains

J_{c} (Q_{2}) = (\begin{matrix} - a \frac{u_{2}}{k} + \frac{2 b u_{2}^{2} v_{2}}{{(h + u_{2}^{2})}^{2}} & - \frac{b u_{2}}{h + u_{2}^{2}} \\ \frac{b c v_{2} (h - u_{2}^{2})}{{(h + u_{2}^{2})}^{2}} - g v_{2} & 0 \end{matrix}),

from which we can easily derive the following result:

\begin{matrix} | J_{c} (Q_{2}) | = \frac{b u_{2} v_{2}}{h + u_{2}^{2}} (\frac{b c (h - u_{2}^{2})}{{(h + u_{2}^{2})}^{2}} - g) . \end{matrix}

(4.4)

From (4.2), we can deduce

h = \frac{b c u_{2}}{g u_{2} + d} - u_{2}^{2} .

(4.5)

Substituting (4.5) into (4.4) obtains

| J_{c} (Q_{2}) | = \frac{(g u_{2} + d) v_{2}}{b c^{2} u_{2}} (b c d - 2 g^{2} u_{2}^{3} - 4 g d u_{2}^{2} - 2 d^{2} u_{2}) .

(4.6)

Let

q (u_{2}) = b c d - 2 g^{2} u_{2}^{3} - 4 g d u_{2}^{2} - 2 d^{2} u_{2}

. Since

p^{'} (u_{*}) = 0

p (u_{*}) < 0

and

q (u_{2})

is monotonically decreasing for

u_{2} > 0

, we can obtain

\begin{matrix} q (u_{2}) < q (u_{*}) & = b c d - 2 g^{2} {u_{*}}^{3} - 4 g d {u_{*}}^{2} - 2 d^{2} u_{*} \\ = g (d h - 2 g {u_{*}}^{3} - d {u_{*}}^{2}) \\ < g (- 2 g {u_{*}}^{3} - d {u_{*}}^{2} - g {u_{*}}^{3} - d {u_{*}}^{2} - (g h - c b) u_{*}) \\ = g u_{*} (- 3 g {u_{*}}^{2} - 2 d u_{*} - g h + c b) = 0 . \end{matrix}

(4.7)

This verifies that

| J_{c} (Q_{2}) | < 0

holds if

Q_{2}

exists, which reads

λ_{1} λ_{2} < 0

. Accordingly,

λ_{1} > (<) 0 \Rightarrow | a r g (λ_{1}) | < (>) \frac{α π}{2}

and

λ_{2} < (>) 0 \Rightarrow | a r g (λ_{1}) | > (<) \frac{α π}{2}

. Thus, the postive equilibrium point

Q_{2} (u_{2}, v_{2})

is always a saddle.

Similarly, for the postive equilibrium point

Q_{1} (u_{1}, v_{1})

, we have

\begin{matrix} q (u_{1}) & = b c d - 2 g^{2} u_{1}^{3} - 4 g d u_{1}^{2} - 2 d^{2} u_{1} \\ > g (d h - 2 g {u_{*}}^{3} - d {u_{*}}^{2}) \\ > g u_{*} (- 3 g {u_{*}}^{2} - 2 d u_{*} - g h + c b) = 0 . \end{matrix}

(4.8)

Thus,

| J_{c} (Q_{1}) | > 0

, which reads

λ_{1} λ_{2} > 0

. In order to determine the signs of

λ_{1}

and

λ_{2}

, we need to further consider the sign of the trace of matrix

J_{c} (Q_{1})

. Note that the trace of

J_{c} (Q_{1})

t r (J_{c} (u_{1}, v_{1})) = - a \frac{u_{1}}{k} + \frac{2 b u_{1}^{2} v_{1}}{{(h + u_{1})}^{2}}

. Note that

\frac{c b u_{1} v_{1}}{h + u_{1}^{2}} - d v_{1} - g u_{1} v_{1} = 0

and

v_{1} = \frac{a}{b} (1 - \frac{u_{1}}{k}) (h + u_{1}^{2})

. So,

t r (J_{c} (u_{1}, v_{1})) = \frac{2 a u_{1}^{2}}{k (h + u_{1}^{2})} (k - \frac{h + 3 u_{1}^{2}}{2 u_{1}}) > (=, <) 0 \Leftrightarrow k > (=, <) \frac{h + 3 u_{1}^{2}}{2 u_{1}} .

Therefore, we can conclude that if

k < \frac{h + 3 u_{1}^{2}}{2 u_{1}}

, then

λ_{1} < 0 (\Rightarrow | a r g (λ_{1}) | > \frac{α π}{2}

) and

λ_{2} < 0 (\Rightarrow | a r g (λ_{2}) | > \frac{α π}{2}

); thus, system (1.3) is stable at

Q_{1} (u_{1}, v_{1})

. If

k > \frac{h + 3 u_{1}^{2}}{2 u_{1}}

, then

λ_{1} > 0 (\Rightarrow | a r g (λ_{1}) | = 0 < \frac{α π}{2}

) and

λ_{2} > 0 (\Rightarrow | a r g (λ_{2}) | = 0 < \frac{α π}{2}

); hence, system (1.3) is unstable at

Q_{1} (u_{1}, v_{1})

. □

For readers’ convenience, we summarize the stability of the equilibrium points of system (1.3) in Table 3.

4.3. Stability Analysis of the Equilibrium Points of System (1.6)

The Jacobian matrix of system (1.6) at any equilibrium point

Q (u, v)

is as follows:

J_{d} (u, v) = (\begin{matrix} 1 + \frac{ρ^{α}}{Γ (α + 1)} (a (1 - \frac{2 u}{k}) - \frac{b v}{h + u^{2}} + \frac{2 b u^{2} v}{{(h + u^{2})}^{2}}) & \frac{ρ^{α}}{Γ (α + 1)} (- \frac{b u}{h + u^{2}}) \\ \frac{ρ^{α}}{Γ (α + 1)} (\frac{b c v (h - u^{2})}{{(h + u^{2})}^{2}} - g v) & 1 + \frac{ρ^{α}}{Γ (α + 1)} (\frac{b c u}{h + u^{2}} - d - g u) \end{matrix}) .

The characteristic polynomial of the Jacobian matrix

J (Q)

reads

F (λ) = λ^{2} - \hat{P} λ + \hat{Q},

where

\begin{matrix} \hat{P} & = 2 + \frac{ρ^{α}}{Γ (α + 1)} (- \frac{a u}{k} + \frac{2 b u^{2} v^{2}}{{(h + u^{2})}^{2}}), \\ \hat{Q} & = 1 + \frac{ρ^{α}}{Γ (α + 1)} (- \frac{a u}{k} + \frac{2 b u^{2} v^{2}}{{(h + u^{2})}^{2}}) + \frac{ρ^{2 α} b u v}{{(Γ (α + 1))}^{2} h + u^{2}} (\frac{b c (h - u^{2})}{{(h + u^{2})}^{2}} - g) . \end{matrix}

Denote

Δ = \frac{ρ^{α}}{Γ (α + 1)}

Ω_{1} = \frac{b u v}{h + u^{2}} (\frac{b c (h - u^{2})}{{(h + u^{2})}^{2}} - g)

, and

Ω_{2} = - \frac{a u}{k} + \frac{2 b u^{2} v^{2}}{{(h + u^{2})}^{2}}

4.3.1. The Stability of Trivial Equilibrium Point $Q_{0} (0, 0)$

Theorem 11.

The following statements about the trivial equilibrium point

Q_{0} (0, 0)

of system (1.6) are true.

If $d < \frac{2}{Δ}$ , then $Q_{0} (0, 0)$ is a saddle.
If $d = \frac{2}{Δ}$ , then $Q_{0} (0, 0)$ is non-hyperbolic.
If $d > \frac{2}{Δ}$ , then $Q_{0} (0, 0)$ is a stable node.

Proof.

Substituting the trivial equilibrium point

Q_{0} (0, 0)

into the Jacobian matrix

J_{d} (u, v)

, we obtain

J_{d} (Q_{0}) = (\begin{matrix} 1 + Δ a & 0 \\ 0 & 1 - Δ d \end{matrix}),

and it is easy to see that the Jacobian matrix

J_{d} (Q_{0})

has two eigenvalues:

λ_{1} = 1 + Δ a > 1

and

| λ_{2} | = | 1 - Δ d | < (=, >) 1

for

d < (=, >) \frac{2}{Δ}

. By using Definition 4, we can derive Theorem 14. □

4.3.2. The Stability of Boundary Equilibrium Point $Q_{k} (k, 0)$

Theorem 12.

The following conclusions for the fixed point

Q_{k} (k, 0)

of system (1.6) are true.

1.

a < \frac{2}{Δ}

, then,

(a): For $d < \frac{c b k}{h + k^{2}} - g k$ or $d > \frac{2}{Δ} + \frac{c b k}{h + k^{2}} - g k$ , $Q_{k} (k, 0)$ is is a saddle;
(b): For $d = \frac{c b k}{h + k^{2}} - g k$ or $d = \frac{2}{Δ} + \frac{c b k}{h + k^{2}} - g k$ , $Q_{k} (k, 0)$ is non-hyperbolic;
(c): For $\frac{c b k}{h + k^{2}} - g k < d < \frac{2}{Δ} + \frac{c b k}{h + k^{2}} - g k$ , $Q_{k} (k, 0)$ is a stable node, i.e., a sink.

2.

a = \frac{2}{Δ}

, then

Q_{k} (k, 0)

is non-hyperbolic.

3.

a > \frac{2}{Δ}

, then,

(a): For $d < \frac{c b k}{h + k^{2}} - g k$ or $d > \frac{2}{Δ} + \frac{c b k}{h + k^{2}} - g k$ , $Q_{k} (k, 0)$ is an unstable node, i.e., a source;
(b): For $d = \frac{c b k}{h + k^{2}} - g k$ or $d = \frac{2}{Δ} + \frac{c b k}{h + k^{2}} - g k$ , $Q_{k} (k, 0)$ is non-hyperbolic;
(c): For $\frac{c b k}{h + k^{2}} - g k < d < \frac{2}{Δ} + \frac{c b k}{h + k^{2}} - g k$ , $Q_{k} (k, 0)$ is a saddle.

Proof.

Substituting boundary equilibrium point

Q_{k} (k, 0)

into the Jacobian matrix

J_{d} (u, v)

, we have

J_{d} (Q_{k}) = (\begin{matrix} 1 - Δ a & - Δ \frac{b k}{h + k^{2}} \\ 0 & 1 + Δ (\frac{c b k}{h + k^{2}} - d - g k) \end{matrix}) .

Now, consider the following three cases.

Case 1:

a < \frac{2}{Δ}

. Then, the two eigenvalues of the Jacobian matrix

J_{c} (Q_{k})

satisfy

| λ_{1} | = | 1 - Δ a | < 1

and

| λ_{2} | = | 1 + Δ (\frac{c b k}{h + k^{2}} - d - g k) | > 1

for

d < \frac{c b k}{h + k^{2}} - g k

d > \frac{2}{Δ} + \frac{c b k}{h + k^{2}} - g k

, so

Q_{k} (k, 0)

is a saddle;

| λ_{2} | = 1

for

d = \frac{c b k}{h + k^{2}} - g k

d = \frac{2}{Δ} + \frac{c b k}{h + k^{2}} - g k

, then

Q_{k} (k, 0)

is non-hyperbolic; and

| λ_{2} | < 1

for

\frac{c b k}{h + k^{2}} - g k < d < \frac{2}{Δ} + \frac{c b k}{h + k^{2}} - g k

, thus

Q_{k} (k, 0)

is a stable node, i.e., a sink.

Case 2:

a = \frac{2}{Δ}

. Then, one eigenvalue of the Jacobian matrix

J_{c} (Q_{k})

meets

| λ_{1} | = | 1 - Δ a | = 1

. Thus,

Q_{k} (k, 0)

is non-hyperbolic.

Case 3:

a > \frac{2}{Δ}

. Then, the two eigenvalues of the Jacobian matrix

J_{c} (Q_{k})

meet

| λ_{1} | = | 1 - Δ a | > 1

and

| λ_{2} | = | 1 + Δ (\frac{c b k}{h + k^{2}} - d - g k) | > 1

for

d < \frac{c b k}{h + k^{2}} - g k

d > \frac{2}{Δ} + \frac{c b k}{h + k^{2}} - g k

, then

Q_{k} (k, 0)

is an unstable node, i.e., a source;

| λ_{2} | = 1

for

d = \frac{c b k}{h + k^{2}} - g k

d = \frac{2}{Δ} + \frac{c b k}{h + k^{2}} - g k

, then

Q_{k} (k, 0)

is non-hyperbolic; and

| λ_{2} | < 1

for

\frac{c b k}{h + k^{2}} - g k < d < \frac{2}{Δ} + \frac{c b k}{h + k^{2}} - g k

, then at this time,

Q_{k} (k, 0)

is a saddle.

The proof is over. □

4.3.3. The Stability of Positive Equilibrium Points $Q_{i} (u_{i}, v_{i}) (i = 1, 2)$

Theorem 13.

When the positive equilibrium

Q_{i} (u_{i}, v_{i}) (i = 1, 2)

exists, let

Ω_{1} (i) = \frac{b u_{i} v_{i}}{h + u_{i}^{2}} (\frac{b c (h - u_{i}^{2})}{{(h + u_{i}^{2})}^{2}} - g)

Ω_{2} (i) = - \frac{a u_{i}}{k} + \frac{2 b u_{i}^{2} v_{i}^{2}}{{(h + u_{i}^{2})}^{2}} (i = 1, 2)

. Then, the results for the positive fixed points

Q_{i} (u_{i}, v_{i}) (i = 1, 2)

of system (1.6) summarized in Table 4 and Table 5 are true.

Proof.

Substituting the positive equilibrium points

Q_{i} (u_{i}, v_{i}) (i = 1, 2)

into the Jacobian matrix

J_{d} (u, v)

, one can see

J_{d} (Q_{i}) = (\begin{matrix} 1 + Δ (a (1 - \frac{2 u_{i}}{k}) - \frac{b v_{i}}{h + u_{i}^{2}} + \frac{2 b u_{i}^{2} v_{i}}{{(h + u_{i}^{2})}^{2}}) & Δ (- \frac{b u_{i}}{h + u_{i}^{2}}) \\ Δ (\frac{b c v_{i} (h - u_{i}^{2})}{{(h + u_{i}^{2})}^{2}} - g v_{i}) & 1 + Δ (\frac{b c u_{i}}{h + u_{i}^{2}} - d - g u_{i}) \end{matrix}) .

We first analyze the stability of the positive equilibrium point

Q_{2} (u_{2}, v_{2})

. By computing

J_{d} (Q_{2})

, we can obtain its characteristic polynomial

F_{2} (λ) = λ^{2} - {\hat{P}}_{2} λ + {\hat{Q}}_{2},

where

\begin{matrix} {\hat{P}}_{2} & = 2 + Δ (- \frac{a u_{2}}{k} + \frac{2 b u_{2}^{2} v_{2}^{2}}{{(h + u_{2}^{2})}^{2}}), \\ {\hat{Q}}_{2} & = 1 + Δ (- \frac{a u_{2}}{k} + \frac{2 b u_{2}^{2} v_{2}^{2}}{{(h + u_{2}^{2})}^{2}}) + Δ^{2} \frac{b u_{2} v_{2}}{h + u_{2}^{2}} (\frac{b c (h - u_{2}^{2})}{{(h + u_{2}^{2})}^{2}} - g) . \end{matrix}

Denote

Ω_{1} (2) = \frac{b u_{2} v_{2}}{h + u_{2}^{2}} (\frac{b c (h - u_{2}^{2})}{{(h + u_{2}^{2})}^{2}} - g)

and

Ω_{2} (2) = - \frac{a u_{2}}{k} + \frac{2 b u_{2}^{2} v_{2}^{2}}{{(h + u_{2}^{2})}^{2}}

. Then,

F_{2} (1) = Δ^{2} Ω_{1} (2), F_{2} (- 1) = Ω_{1} (2) Δ^{2} + 2 Ω_{2} (2) Δ + 4 .

From (4.7), we know that

Ω_{1} (2) < 0

; hence,

F_{2} (1) < 0

, which shows that

F_{2} (λ) = 0

has one root lying in

(1, \infty)

. For the value of

F_{2} (- 1)

, we consider the following three cases for discussion:

Case 1:

0 < Δ < \frac{- Ω_{2} (2) - \sqrt{Ω_{2}^{2} (2) - 4 Ω_{1} (2)}}{Ω_{1} (2)}

. Then,

F_{2} (- 1) > 0

. Lemma 3 reads that

| λ_{1} | > 1

and

| λ_{2} | < 1

, which shows that

Q_{2}

is a saddle.

Case 2:

Δ = \frac{- Ω_{2} (2) - \sqrt{Ω_{2}^{2} (2) - 4 Ω_{1} (2)}}{Ω_{1} (2)}

. Then,

F_{2} (- 1) = 0

. This shows that

Q_{2}

is non-hyperbolic.

Case 3:

Δ > \frac{- Ω_{2} (2) - \sqrt{Ω_{2}^{2} (2) - 4 Ω_{1} (2)}}{Ω_{1} (2)}

. Then,

F_{2} (- 1) < 0

. Lemma 3 tells us that

| λ_{1} | > 1

and

| λ_{2} | > 1

, so

Q_{2}

is an unstable node.

Next, we discuss the stability of the positive equilibrium point

Q_{1} (u_{1}, v_{1})

. Similarly, after calculating

J_{d} (Q_{1})

, we obtain its characteristic polynomial

F_{1} (λ) = λ^{2} - {\hat{P}}_{1} λ + {\hat{Q}}_{1},

where

\begin{matrix} {\hat{P}}_{1} & = 2 + Δ (- \frac{a u_{1}}{k} + \frac{2 b u_{1}^{2} v_{1}^{2}}{{(h + u_{1}^{2})}^{2}}), \\ {\hat{Q}}_{1} & = 1 + Δ (- \frac{a u_{1}}{k} + \frac{2 b u_{1}^{2} v_{1}^{2}}{{(h + u_{1}^{2})}^{2}}) + Δ^{2} \frac{b u_{1} v_{1}}{h + u_{1}^{2}} (\frac{b c (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g) . \end{matrix}

Let

Ω_{1} (1) = \frac{b u_{1} v_{1}}{h + u_{1}^{2}} (\frac{b c (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g)

and

Ω_{2} (1) = - \frac{a u_{1}}{k} + \frac{2 b u_{1}^{2} v_{1}^{2}}{{(h + u_{1}^{2})}^{2}}

. Then,

\begin{matrix} F_{1} (1) = Δ^{2} Ω_{1} (1), F_{1} (- 1) = Ω_{1} (1) Δ^{2} + 2 Ω_{2} (1) Δ + 4 . \end{matrix}

From (4.8), we know that

Ω_{1} (1) > 0

, so

F_{1} (1) > 0

. Note also that

{\hat{Q}}_{1} > (=, <) 1 ⟺ Δ > (=, <) - \frac{Ω_{2} (1)}{Ω_{1} (1)}

. Now, consider the following three cases for discussion.

Case 1:

Ω_{2}^{2} (1) - 4 Ω_{1} (1) < 0

. This implies that

F_{1} (- 1) > 0

. Thus, we consider the following three subcases for discussion:

Subcase 1.

Δ < - \frac{Ω_{2} (1)}{Ω_{1} (1)}

. Then,

{\hat{Q}}_{1} < 1 .

Lemma 3 reads

| λ_{1} | < 1

and

| λ_{2} | < 1

. Thus,

Q_{1}

is a stable node, i.e., a sink.

Subcase 2.

Δ = - \frac{Ω_{2} (1)}{Ω_{1} (1)}

. Then,

{\hat{Q}}_{1} = 1

Subcase 3.

Δ > - \frac{Ω_{2} (1)}{Ω_{1} (1)} .

Then,

{\hat{Q}}_{1} > 1 .

Thus,

| λ_{1} | > 1

| λ_{2} | > 1

, and hence

Q_{1}

is an unstable node, i.e., a source.

Case 2:

Ω_{2}^{2} (1) - 4 Ω_{1} (1) = 0

. Consider the following three subcases:

Subcase 1.

Δ < - \frac{Ω_{2} (1)}{Ω_{1} (1)} .

Then,

F_{1} (- 1) > 0, {\hat{Q}}_{1} < 1

. Hence,

| λ_{1} | < 1

| λ_{2} | < 1

, and so

Q_{1}

is a stable node, i.e., a sink.

Subcase 2.

Δ = - \frac{Ω_{2} (1)}{Ω_{1} (1)}

. Then,

F_{1} (- 1) = 0

. Thus,

Q_{1}

is non-hyperbolic.

Subcase 3.

Δ > - \frac{Ω_{2} (1)}{Ω_{1} (1)} .

Then,

F_{1} (- 1) > 0

and

{\hat{Q}}_{1} > 1 .

It follows from Lemma 3 that

| λ_{1} | > 1

and

| λ_{2} | > 1

. Thus,

Q_{1}

is an unstable node, i.e., a source.

Case 3:

Ω_{2}^{2} (1) - 4 Ω_{1} (1) > 0

. Then,

Ω_{2} (1) > 2 \sqrt{Ω_{1} (1)}

Ω_{2} (1) < - 2 \sqrt{Ω_{1} (1)}

. Consider the following two subcases:

Subcase 1.

Ω_{2} (1) > 2 \sqrt{Ω_{1} (1)}

. Then,

Δ > 0 > \frac{- Ω_{2} (1) + \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)} .

So,

F_{1} (- 1) > 0, {\hat{Q}}_{1} > 1

, which, in view of Lemma 3, implies

| λ_{1} | > 1

and

| λ_{2} | > 1

. Thus,

Q_{1}

is an unstable node, i.e., a source.

Subcase 2.

Ω_{2} (1) < - 2 \sqrt{Ω_{1} (1)}

. We further consider the following five subsubcases.

Subcase 1.

0 < Δ < \frac{- Ω_{2} (1) - \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)} < - \frac{Ω_{2} (1)}{Ω_{1} (1)} .

Then,

F_{1} (- 1) > 0, {\hat{Q}}_{1} < 1

, indicating

| λ_{1} | < 1

| λ_{2} | < 1

. Hence,

Q_{1}

is a stable node, i.e., a sink.

Subcase 2.

Δ = \frac{- Ω_{2} (1) - \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)} .

Then,

F_{1} (- 1) = 0

and so

Q_{1}

is non-hyperbolic.

Subcase 3.

\frac{- Ω_{2} (1) - \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)} < Δ < \frac{- Ω_{2} (1) + \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)}

. Then,

F_{1} (- 1) < 0

. In light of Lemma 3,

| λ_{1} | > (<) 1

| λ_{2} | < (>) 1

. Thus,

Q_{1}

is a saddle.

Subcase 4.

Δ = \frac{- Ω_{2} (1) + \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)}

. Then,

F_{1} (- 1) = 0

and so

Q_{1}

is non-hyperbolic.

Subcase 5.

Δ > \frac{- Ω_{2} (1) + \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)} > - \frac{Ω_{2} (1)}{Ω_{1} (1)} .

Then,

F_{1} (- 1) > 0, {\hat{Q}}_{1} > 1

and so

| λ_{1} | > 1

| λ_{2} | > 1

. Therefore,

Q_{1}

is an unstable node, i.e., a source.

The proof is finished. □

5. Bifurcation Analysis

In this section, we, respectively, analyze the existence of bifurcations in the positive equilibrium point

Q_{1} (u_{1}, v_{1})

of systems (1.3) and (1.6).

5.1. Bifurcation Analysis of the Positive Equilibrium Point $Q_{1} (u_{1}, v_{1})$ in System (1.3)

In Section 3, we see that the Jacobian matrix of system (1.3) at the positive equilibrium point

Q_{1} (u_{1}, v_{1})

is as follows:

J (Q_{1}) = (\begin{matrix} - a \frac{u_{1}}{k} + \frac{2 b u_{1}^{2} v_{1}}{{(h + u_{1}^{2})}^{2}} & - \frac{b u_{1}}{h + u_{1}^{2}} \\ \frac{c b v_{1} (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g v_{1} & 0 \end{matrix})

The characteristic equation of the Jacobian matrix

J (Q_{1})

is given by

λ^{2} - (- a \frac{u_{1}}{k} + \frac{2 b u_{1}^{2} v_{1}}{{(h + u_{1}^{2})}^{2}}) λ + \frac{b u_{1}}{h + u_{1}^{2}} (\frac{c b v_{1} (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g v_{1}) = 0 .

(5.1)

Substituting

v_{1} = \frac{a}{b} (1 - \frac{u_{1}}{k}) (h + u_{1}^{2})

into Equation (5.1), we have

\begin{matrix} λ^{2} - a u_{1} (\frac{2 k u_{1} - 3 u_{1}^{2} - h}{k (h + u_{1}^{2})}) λ + a u_{1} (1 - \frac{u_{1}}{k}) (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g) = 0 . \end{matrix}

(5.2)

Take k as the bifurcation parameter of system (1.3). If k takes a critical value,

k_{0} > 0

, such that the corresponding eigenvalues are

λ_{k_{0}} = r e^{i γ}

, where

γ = \pm \frac{α π}{2}

, then a bifurcation occurs. Now, we look for

k_{0}

such that

λ_{k_{0}}

satisfies Equation (5.2).

Substituting

λ_{k_{0}}

into (5.2), we can obtain the following equation:

\begin{matrix} r^{2} e^{2 i γ} - a u_{1} (\frac{2 k_{0} u_{1} - 3 u_{1}^{2} - h}{k_{0} (h + u_{1}^{2})}) r e^{i γ} + a u_{1} (1 - \frac{u_{1}}{k_{0}}) (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g) = 0 . \end{matrix}

Namely,

\begin{matrix} r^{2} (c o s (2 γ) + i s i n (2 γ)) - a u_{1} (\frac{2 k_{0} u_{1} - 3 u_{1}^{2} - h}{k_{0} (h + u_{1}^{2})}) r (c o s γ + i s i n γ) \\ + a u_{1} (1 - \frac{u_{1}}{k_{0}}) (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g) = 0 . \end{matrix}

Hence,

\begin{matrix} \begin{matrix} \{\begin{matrix} r^{2} c o s (2 γ) - a u_{1} (\frac{2 k_{0} u_{1} - 3 u_{1}^{2} - h}{k_{0} (h + u_{1}^{2})}) r c o s γ + a u_{1} (1 - \frac{u_{1}}{k_{0}}) (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g) = 0, \\ r^{2} s i n (2 γ) - a u_{1} (\frac{2 k_{0} u_{1} - 3 u_{1}^{2} - h}{k_{0} (h + u_{1}^{2})}) r s i n γ = 0 . \end{matrix} \end{matrix} . \end{matrix}

(5.3)

Since we are interested in non-zero solutions for r in (5.3), from the second equation of (5.3) we can derive

r = a u_{1} (\frac{2 k_{0} u_{1} - 3 u_{1} - h}{2 c o s γ k_{0} (h + u_{1}^{2})})

. After substitution into the first equation of (5.3), one has

{(\frac{2 k_{0} u_{1} - 3 u_{1}^{2} - h}{k_{0} (h + u_{1}^{2})})}^{2} \frac{- a^{2} u_{1}^{2}}{4 c o s^{2} γ} + a u_{1} (1 - \frac{u_{1}}{k_{0}}) (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g) = 0 .

So,

a u_{1} {(\frac{2 k_{0} u_{1} - 3 u_{1}^{2} - h}{k_{0} (h + u_{1}^{2})})}^{2} = 4 c o s^{2} γ (1 - \frac{u_{1}}{k_{0}}) (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g),

namely,

\begin{matrix} ω_{1} k_{0}^{2} - ω_{2} k_{0} + ω_{3} = 0, \end{matrix}

(5.4)

where

\begin{matrix} ω_{1} = 4 a u_{1}^{3} - 4 c o s^{2} γ (c b (h - u_{1}^{2}) - g {(h + u_{1}^{2})}^{2}), \\ ω_{2} = 4 a u_{1}^{2} (3 u_{1}^{2} + h) - 4 c o s^{2} γ u_{1} (c b (h - u_{1}^{2}) - g {(h + u_{1}^{2})}^{2}), \\ ω_{3} = a u_{1} {(3 u_{1}^{2} + h)}^{2} > 0 . \end{matrix}

| J (Q_{1}) | > 0

implies

c b (h - u_{1}^{2}) - g {(h + u_{1}^{2})}^{2} > 0

. Let

γ_{0} = \frac{a {u_{1}}^{3}}{c b (h - u_{1}^{2}) - g {(h + u_{1}^{2})}^{2}}

. Then,

γ_{0} > 0

. After a lengthy and tedious calculation, we can classify the following three cases for further discussion:

Case 1:

c o s^{2} γ > γ_{0}

. Then, we have

ω_{1} < 0

. Due to

k_{0} > 0

, we can obtain

k_{0} = \frac{ω_{2} - \sqrt{ω_{2}^{2} - 4 ω_{1} ω_{3}}}{2 ω_{1}}

Case 2:

c o s^{2} γ = γ_{0}

. Then,

ω_{1} = 0

. Noticing

k_{0} > 0

, we can obtain

k_{0} = \frac{ω_{3}}{ω_{2}}

Case 3:

c o s^{2} γ < γ_{0}

. Then, we have

ω_{1} > 0

and

ω_{2} > 0

. Calculate

ω_{2}^{2} - 4 ω_{1} ω_{3}

to obtain

ω_{2}^{2} - 4 ω_{1} ω_{3} = 16 u_{1} c o s^{2} γ [u_{1} c o s γ^{2} (c b (h - u_{1}^{2}) - g (h + u_{1}^{2})) + a (3 u_{1}^{2} + h) (u_{1}^{2} + h)] > 0 .

Then, we can derive that

k_{0}

has two values:

\frac{ω_{2} \pm \sqrt{ω_{2}^{2} - 4 ω_{1} ω_{3}}}{2 ω_{1}}

In any case, the critical value

k_{0}

always exists. Next, we prove that system (1.3) satisfies the conditions of Theorem 4 at the positive equilibrium point

Q_{1} (u_{1}, v_{1})

From the existence of

k_{0}

, we see that

| a r g (λ_{i} (k_{0})) | = \frac{α π}{2} (i = 1, 2)

; hence, the first condition in Theorem 4 holds true. The Jacobian matrix of system (1.3) has only two eigenvalues; thus, we do not need to consider the second condition in Theorem 4. Next, we focus on proving that system (1.3) satisfies the third condition of Theorem 4. Take the derivative of Equation (5.2) with respect to k to obtain

\begin{matrix} 2 λ \frac{d λ}{d k} - a u_{1} (\frac{2 u_{1}}{h + u_{1}^{2}} - \frac{3 u_{1}^{2} + h}{k (h + u_{1}^{2})}) \frac{d λ}{d k} - a u_{1} \frac{3 u_{1}^{2} + h}{k^{2} (h + u_{1}^{2})} λ + \\ \frac{a u_{1}^{2}}{k^{2}} (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g) = 0 . \end{matrix}

(5.5)

Thus,

\begin{matrix} \frac{d λ}{d k} = \frac{\frac{a u_{1} (3 u_{1}^{2} + h)}{k^{2} (h + u_{1}^{2})} λ - \frac{a u_{1}^{2}}{k^{2}} (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g)}{2 λ + \frac{3 a u_{1}^{3} + a u_{1} h - 2 a k u_{1}^{2}}{k (h + u_{1}^{2})}} . \end{matrix}

(5.6)

It suffices for us to verify

λ_{1, k_{0}} = r e^{i \frac{α π}{2}}

. Substituting

k = k_{0}

and

λ_{1, k_{0}} = r e^{i \frac{α π}{2}}

into the right-hand side of (5.6) obtains

\begin{matrix} \frac{d λ}{d k} |_{k = k_{0}} & = \frac{\frac{a u_{1} (3 u_{1}^{2} + h)}{k_{0}^{2} (h + u_{1}^{2})} r e^{i γ} - \frac{a u_{1}^{2}}{k_{0}^{2}} (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g)}{2 r e^{i γ} + \frac{3 a u_{1}^{3} + a u_{1} h - 2 a k_{0} u_{1}^{2}}{k_{0} (h + u_{1}^{2})}} \\ = \frac{\frac{a u_{1} (3 u_{1}^{2} + h)}{k_{0}^{2} (h + u_{1}^{2})} r (c o s γ + i s i n γ) - \frac{a u_{1}^{2}}{k_{0}^{2}} (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g)}{2 r (c o s γ + i s i n γ) + \frac{3 a u_{1}^{3} + a u_{1} h - 2 a k_{0} u_{1}^{2}}{k_{0} (h + u_{1}^{2})}} \\ = \frac{ψ_{1} + ψ_{2} i}{ψ_{3}^{2} + ψ_{4}^{2}}, \end{matrix}

(5.7)

where

\begin{matrix} ψ_{1} = & (\frac{a u_{1} (3 u_{1}^{2} + h)}{k_{0}^{2} (h + u_{1}^{2})} r c o s γ - \frac{a u_{1}^{2}}{k_{0}} (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g)) \\ \times (2 r c o s γ + \frac{3 a u_{1}^{3} + a u_{1} h - 2 a k_{0} u_{1}^{2}}{k_{0} (h + u_{1}^{2})}) + \frac{a u_{1} (3 u_{1}^{2} + h)}{k_{0}^{2} (h + u_{1}^{2})} 2 r^{2} s i n^{2} γ, \\ ψ_{2} = & (\frac{a (3 u_{1}^{2} + h) (3 u_{1}^{2} + h - 2 k_{0} u_{1})}{k_{0}^{3} {(h + u_{1}^{2})}^{2}} + \frac{2}{k_{0}} (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g)) a u_{1}^{2} r s i n γ, \\ ψ_{3} = & 2 r c o s γ + \frac{3 a u_{1}^{3} + a u_{1} h - 2 a k_{0} u_{1}^{2}}{k_{0} (h + u_{1}^{2})}, ψ_{4} = 2 r s i n γ . \end{matrix}

Denote

λ (k) = m (k) + i n (k)

, then

a r g (λ (k)) = a r c t a n (\frac{n}{m})

. By differentiating

a r g (λ (k))

with respect to k, we get

\begin{matrix} \frac{d}{d k} a r g (λ (k)) = \frac{m n^{'} - n m^{'}}{m^{2} + n^{2}} = \frac{W (m, n)}{{| λ (k) |}^{2}} \end{matrix}

(5.8)

where

W (m, n) = |\begin{matrix} m & n \\ m^{'} & n^{'} \end{matrix}|

Therefore,

\begin{matrix} \frac{d}{d k} a r g (λ (k)) |_{k = k_{0}} & = \frac{W (m (k_{0}), n (k_{0}))}{| λ (k_{0}) |^{2}} \\ = \frac{1}{| λ (k_{0}) |^{2}} |\begin{matrix} m (k_{0}) & n (k_{0}) \\ m^{'} (k_{0}) & n^{'} (k_{0}) \end{matrix}| \\ = \frac{1}{| λ (k_{0}) |^{2}} |\begin{matrix} r c o s γ & r s i n γ \\ \frac{ψ_{1}}{ψ_{3}^{2} + ψ_{4}^{2}} & \frac{ψ_{2}}{ψ_{3}^{2} + ψ_{4}^{2}} \end{matrix}| . \end{matrix}

(5.9)

We can easily deduce that when

k_{0}

exists,

| λ (k_{0}) |^{2} \neq \infty

holds true.

Next, we prove the conditions under which

W (m, n) (k_{0}) \neq 0

holds true.

In fact,

\begin{matrix} W (m, n) (k_{0}) = |\begin{matrix} r c o s γ & r s i n γ \\ \frac{ψ_{1}}{ψ_{3}^{2} + ψ_{4}^{2}} & \frac{ψ_{2}}{ψ_{3}^{2} + ψ_{4}^{2}} \end{matrix}| \neq 0 \\ ⟺ \frac{1}{ψ_{3}^{2} + ψ_{4}^{2}} (r c o s γ ψ_{2} - r s i n γ ψ_{1}) \neq 0 \\ ⟺ c o s γ ψ_{2} - s i n γ ψ_{1} \neq 0 \\ ⟺ 4 u_{1} r c o s γ (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g) + \\ u_{1} (\frac{c b (h - u_{1}^{2})}{{(h + u_{1}^{2})}^{2}} - g) (\frac{3 a u_{1}^{3} + a u_{1} h - 2 a k_{0} u_{1}^{2}}{k_{0} (h + u_{1}^{2})}) - 2 r^{2} \frac{3 u_{1}^{2} + h}{k_{0} (h + u_{1}^{2})} \neq 0 \\ ⟺ k_{0} \neq k_{h}, \end{matrix}

(5.10)

where

k_{h} = \frac{2 r^{2} (3 u_{1}^{2} + h) {(h + u_{1}^{2})}^{2} - u_{1} (c b (h - u_{1}^{2}) - g (h + u_{1}^{2})) (3 a u_{1}^{3} + a u_{1} h)}{u_{1} (c b (h - u_{1}^{2}) - g (h + u_{1}^{2})) (4 r c o s γ (h + u_{1}^{2}) - 2 a u_{1}^{2})} .

(5.11)

This is true by adding the assumption

k_{0} \neq k_{h}

. So, summarizing the above analysis, one has the following results.

Theorem 14.

Suppose that all parameters in system (1.3) are positive. Let

R_{0}, R, k_{0}, k_{h}

be defined as above. If

g h - c b < 0

R_{0} < R

u_{1} < k

k_{0} \neq k_{h}

, then system (1.3) undergoes a fractional Hopf bifurcation at the positive equilibrium point

Q_{1} (u_{1}, v_{1})

5.2. Bifurcation Analysis of the Positive Equilibrium Point $Q_{1} (u_{1}, v_{1})$ in System (1.6)

In this subsection, we study the bifurcation problems of system (1.6) at the positive equilibrium point

Q_{1} (u_{1}, v_{1})

by using the center manifold theorem and local bifurcation theory.

5.2.1. Neimark–Sacker Bifurcation at the Fixed Point $Q_{1} (u_{1}, v_{1})$

From Case 1 in the proof of Theorem 13 for the stability of the positive equilibrium point

Q_{1} (u_{1}, v_{1})

, we see that the dimension numbers for the stable manifold and unstable manifold of system (1.6) at the positive equilibrium point

Q_{1} (u_{1}, v_{1})

change when

Δ

varies in the vicinity of

Δ_{0}

(correspondingly,

ρ

varies in the vicinity of

ρ_{0}

) for

Δ < - \frac{4}{Ω_{2} (1)}

and

Ω_{2}^{2} (1) - 4 Ω_{1} (1) < 0

, where

Δ_{0} = = - \frac{Ω_{2} (1)}{Ω_{1} (1)}, ρ_{0} = {(Γ (1 + α) Δ_{0})}^{\frac{1}{α}} .

(5.12)

Thus, a bifurcation, to be shown to be Neimark–Sacker, may occur. Let

{\hat{N S}}_{Q_{1}} = \{(a, b, c, d, g, h, k, α, ρ) \in R_{+}^{7} | Δ < - \frac{4}{Ω_{2} (1)}, Ω_{2}^{2} (1) - 4 Ω_{1} (1) < 0\} .

To analyze the Neimark–Sacker bifurcation, we perform the following.

Let

x_{n} = u_{n} - u_{1}

y_{n} = v_{n} - v_{1}

, which transforms the fixed point

Q_{1} (u_{1}, v_{1})

to the origin

O (0, 0)

. Assume that

ρ^{*}

is a small perturbation of

ρ

with

| ρ^{*} | ≪ 1

. After shifting and perturbation, system (1.6) takes the following form:

\begin{matrix} \{\begin{matrix} x_{n + 1} & = x_{n} + \frac{{(ρ + ρ^{*})}^{α}}{Γ (α + 1)} (a (x_{n} + u_{1}) (1 - \frac{(x_{n} + u_{1})}{k}) - \frac{b (x_{n} + u_{1}) (y_{n} + v_{1})}{h + {(x_{n} + u_{1})}^{2}}), \\ y_{n + 1} & = y_{n} + \frac{{(ρ + ρ^{*})}^{α}}{Γ (α + 1)} (\frac{c b (x_{n} + u_{1}) (y_{n} + v_{1})}{h + {(x_{n} + u_{1})}^{2}} - d (y_{n} + v_{1}) - g (x_{n} + u_{1}) (y_{n} + v_{1})) . \end{matrix} \end{matrix}

(5.13)

Using the Taylor series expansion of system (5.13) at

O (0, 0)

to the third order results in the following system:

\{\begin{matrix} x_{n + 1} = & ϵ_{10} x_{n} + ϵ_{01} y_{n} + ϵ_{20} x_{n}^{2} + ϵ_{11} x_{n} y_{n} + ϵ_{02} y_{n}^{2} \\ + ϵ_{30} x_{n}^{3} + ϵ_{21} x_{n}^{2} y_{n} + ϵ_{12} x_{n} y_{n}^{2} + ϵ_{03} y_{n}^{3} + o (ρ_{1}^{3}), \\ y_{n + 1} = & ζ_{10} x_{n} + ζ_{01} y_{n} + ζ_{20} x_{n}^{2} + ζ_{11} x_{n} y_{n} + ζ_{02} y_{n}^{2} \\ + ζ_{30} x_{n}^{3} + ζ_{21} x_{n}^{2} y_{n} + ζ_{12} x_{n} y_{n}^{2} + ζ_{03} y_{n}^{3} + o (ρ_{1}^{3}), \end{matrix}

(5.14)

where

ρ_{1} = \sqrt{x_{n}^{2} + y_{n}^{2}}

\begin{matrix} ϵ_{10} = & 1 + \frac{ρ^{α}}{Γ (α + 1)} (a (1 - \frac{2 u_{1}}{k}) - b v_{1} (\frac{h - u_{1}^{2}}{{(u_{1}^{2} + h)}^{2}})), ϵ_{01} = - \frac{ρ^{α} b u_{1}}{Γ (α + 1) (u_{1}^{2} + h)}, \\ ϵ_{20} = & \frac{ρ^{α}}{Γ (α + 1)} ((3 h - u_{1}^{2}) \frac{b u_{1} v_{1}}{{(u_{1}^{2} + h)}^{3}} - \frac{a}{k}), ϵ_{02} = 0, ϵ_{11} = \frac{ρ^{α} b}{Γ (α + 1)} (\frac{u_{1}^{2} - h}{{(u_{1}^{2} + h)}^{2}}), \\ ϵ_{30} = & \frac{ρ^{α} b v_{1}}{Γ (α + 1) {(u_{1}^{2} + h)}^{3}} (h + \frac{2 u_{1}^{2} (3 u_{1}^{2} - h)}{u_{1}^{2} + h} - 5 u_{1}^{2}), ϵ_{03} = 0, \\ ϵ_{21} = & \frac{ρ^{α} b u_{1}}{Γ (α + 1) {(u_{1}^{2} + h)}^{2}} (2 - \frac{3 (u_{1}^{2} - h)}{u_{1}^{2} + h}), ϵ_{12} = 0, \\ ζ_{10} = & \frac{ρ^{α} v_{1}}{Γ (α + 1)} (\frac{b c}{u_{1}^{2} + h} (1 - \frac{2 u_{1}^{2}}{u_{1}^{2} + h}) - g), ζ_{01} = 1 - \frac{ρ^{α}}{Γ (α + 1)} (d + g u_{1} - \frac{b c u_{1}}{u_{1}^{2} + h}), \\ ζ_{20} = & \frac{ρ^{α} b c u_{1} v_{1}}{Γ (α + 1) {(u_{1}^{2} + h)}^{2}}, ζ_{02} = 0, ζ_{11} = \frac{ρ^{α}}{Γ (α + 1)} (\frac{b c}{u_{1}^{2} + h} (1 - \frac{2 u_{1}^{2}}{u_{1}^{2} + h}) - g), \\ ζ_{30} = & \frac{ρ^{α} b c v_{1}}{Γ (α + 1) {(u_{1}^{2} + h)}^{2}} (\frac{4 u_{1}^{2}}{u_{1}^{2} + h} + 4 u_{1}^{2} (1 - 2 u_{1}^{2}) - 1), ζ_{03} = 0 \\ ζ_{21} = & \frac{ρ^{α} b c u_{1}}{Γ (α + 1) {(u_{1}^{2} + h)}^{2}} (\frac{4 u_{1}^{2}}{u_{1}^{2} + h} - 3), ζ_{12} = 0 . \end{matrix}

The characteristic equation of the linearized equation of system (5.14) is

λ^{2} + {\hat{P}}_{1} (ρ^{*}) λ + {\hat{Q}}_{1} (ρ^{*}) = 0,

where

{\hat{P}}_{1} (ρ^{*}) = 2 + Δ (ρ^{*}) Ω_{2} (1)

{\hat{Q}}_{1} (ρ^{*}) = 1 + Δ (ρ^{*}) Ω_{2} (1) + Δ^{2} (ρ^{*}) Ω_{1} (1)

Δ (ρ^{*}) = \frac{{(ρ + ρ^{*})}^{α}}{Γ (α + 1)}

. Noting that

Δ (0) = Δ_{0}

and

{\hat{P}}_{1}^{2} (0) - 4 {\hat{Q}}_{1} (0) = Δ_{0}^{2} (Ω_{2}^{2} (1) - 4 Ω_{1} (1)) < 0

, the two roots of the characteristic equation are

λ_{1, 2} (ρ^{*}) = \frac{- {\hat{P}}_{1} (ρ^{*}) \pm i \sqrt{4 {\hat{Q}}_{1} (ρ^{*}) - {\hat{P}}_{1}^{2} (ρ^{*})}}{2} .

Moreover,

(| λ_{1, 2} (ρ^{*}) |) |_{ρ^{*} = 0} = \sqrt{{\hat{Q}}_{1} (ρ^{*})} |_{ρ^{*} = 0} = 1,

(\frac{d | λ_{1, 2} (ρ^{*}) |}{d ρ^{*}}) |_{ρ^{*} = 0} = - \frac{ρ^{α - 1} Ω_{2} (1)}{2 Γ (α)} \neq 0 .

It is obvious that

λ_{1, 2}^{i} (0) \neq 1,

for

i = 1, 2, 3, 4

. Thus, the transversal and nondegenerate conditions hold for a Neimark–Sacker bifurcation to occur.

In order to derive the normal form of system (5.14), let

T = (\begin{matrix} 0 & ϵ_{01} \\ μ & 1 - ω \end{matrix}),

where

ω = - \frac{{\hat{P}}_{1} (ρ^{*})}{2}

μ = \frac{\sqrt{4 {\hat{Q}}_{1} (ρ^{*}) - {\hat{P}}_{1}^{2} (ρ^{*})}}{2}

. Then, we have

T^{- 1} = (\begin{matrix} \frac{ω - 1}{μ ϵ_{01}} & \frac{1}{μ} \\ \frac{1}{ϵ_{01}} & 0 \end{matrix}) .

Change the variables to

{(x, y)}^{T} = T {(U, V)}^{T},

then, system (5.14) changes to the following form:

(\begin{matrix} U \\ V \end{matrix}) \to (\begin{matrix} ω & - μ \\ μ & ω \end{matrix}) (\begin{matrix} U \\ V \end{matrix}) + (\begin{matrix} F (U, V) + o (ρ_{2}^{3}) \\ G (U, V) + o (ρ_{2}^{3}) \end{matrix}),

where

ρ_{1} = \sqrt{x_{n}^{2} + y_{n}^{2}}

\begin{matrix} F (U, V) = η_{20} x^{2} + η_{11} x y + η_{02} y^{2} + η_{30} x^{3} + η_{21} x^{2} y + η_{12} x y^{2} + η_{03} y^{3}, \\ G (U, V) = θ_{20} x^{2} + θ_{11} x y + θ_{02} y^{2} + θ_{30} x^{3} + θ_{21} x^{2} y + θ_{12} x y^{2} + θ_{03} y^{3}, \\ x = ϵ_{01} V, y = μ U + (1 - ω) V, \\ η_{20} = \frac{ϵ_{20} (ω - 1)}{μ ϵ_{01}} + \frac{ζ_{20}}{μ}, η_{02} = \frac{ϵ_{02} (ω - 1)}{μ ϵ_{01}} + \frac{ζ_{02}}{μ}, η_{11} = \frac{ϵ_{11} (ω - 1)}{μ ϵ_{01}} + \frac{ϵ_{11}}{μ}, \\ η_{30} = \frac{ϵ_{30} (ω - 1)}{μ ϵ_{01}} + \frac{ζ_{30}}{μ}, η_{03} = \frac{ϵ_{03} (ω - 1)}{μ ϵ_{01}} + \frac{ζ_{03}}{μ}, η_{21} = \frac{ϵ_{21} (ω - 1)}{μ ϵ_{01}} + \frac{ζ_{21}}{μ}, \\ η_{12} = \frac{ϵ_{12} (ω - 1)}{μ ϵ_{01}} + \frac{ζ_{12}}{μ}, \\ θ_{20} = \frac{ϵ_{20}}{ϵ_{01}}, θ_{02} = \frac{ϵ_{11}}{ϵ_{01}}, θ_{11} = \frac{ϵ_{02}}{ϵ_{01}}, θ_{30} = \frac{ϵ_{30}}{ϵ_{01}}, θ_{03} = \frac{ϵ_{03}}{ϵ_{01}}, θ_{21} = \frac{ϵ_{21}}{ϵ_{01}}, θ_{12} = \frac{ϵ_{12}}{ϵ_{01}} . \end{matrix}

Furthermore,

\begin{matrix} F_{U U} {|_{(0, 0)} = 2 η_{02} μ^{3}, F_{U V} |}_{(0, 0)} = η_{11} ϵ_{01} μ + 2 η_{02} μ (1 - ω), \\ F_{V V} {|_{(0, 0)} = 2 η_{02} ϵ_{01}^{2} + 2 η_{11} ϵ_{01} (1 - ω), F_{U U U} |}_{(0, 0)} = 6 η_{03} μ^{3}, \\ F_{U U V} |_{(0, 0)} = 2 η_{21} ϵ_{01} μ^{2} + 6 η_{03} μ^{2} (1 - ω), \\ F_{U V V} |_{(0, 0)} = 2 η_{21} ϵ_{01}^{2} μ + 4 η_{12} ϵ_{01} μ (1 - ω) + 6 η_{03} μ {(1 - ω)}^{2}, \\ F_{V V V} |_{(0, 0)} = 4 {(1 - ω)}^{3} + 6 η_{30} ϵ_{01}^{3} + 4 η_{21} ϵ_{01}^{2} (1 - ω) + 6 η_{12} ϵ_{01} {(1 - ω)}^{2}, \\ G_{U U} {|_{(0, 0)} = 2 θ_{02} μ^{3}, G_{U V} |}_{(0, 0)} = θ_{11} ϵ_{01} μ + 2 θ_{02} μ (1 - ω), \\ G_{V V} {|_{(0, 0)} = 2 θ_{02} ϵ_{01}^{2} + 2 θ_{11} ϵ_{01} (1 - ω), G_{U U U} |}_{(0, 0)} = 6 θ_{03} μ^{3}, \end{matrix}

\begin{matrix} G_{U U V} |_{(0, 0)} = 2 θ_{21} ϵ_{01} μ^{2} + 6 θ_{03} μ^{2} (1 - ω), \\ G_{U V V} |_{(0, 0)} = 2 θ_{21} ϵ_{01}^{2} μ + 4 θ_{12} ϵ_{01} μ (1 - ω) + 6 θ_{03} μ {(1 - ω)}^{2}, \\ G_{V V V} |_{(0, 0)} = 4 {(1 - ω)}^{3} + 6 θ_{30} ϵ_{01}^{3} + 4 θ_{21} ϵ_{01}^{2} (1 - ω) + 6 θ_{12} ϵ_{01} {(1 - ω)}^{2} . \end{matrix}

To determine the stability and direction of the bifurcated closed orbit of system (1.4), the following discriminating quantity L should be calculated and not to be zero, where

L = - R e (\frac{(1 - 2 λ_{1}) λ_{2}^{2}}{1 - λ_{1}} τ_{20} τ_{11}) - \frac{1}{2} | τ_{11} |^{2} - {| τ_{02} |}^{2} + R e (λ_{2} τ_{21}),

(5.15)

\begin{matrix} τ_{20} & = \frac{1}{8} [F_{X X} - F_{Y Y} + 2 G_{X Y} + i (G_{X X} - G_{Y Y} - 2 F_{X Y})] |_{(0, 0)}, \\ τ_{11} & = \frac{1}{4} [F_{X X} + F_{Y Y} + i (G_{X X} + G_{Y Y})] |_{(0, 0)}, \\ τ_{02} & = \frac{1}{8} [F_{X X} - F_{Y Y} - 2 G_{X Y} + i (G_{X X} - G_{Y Y} + 2 F_{X Y})] |_{(0, 0)}, \\ τ_{21} & = \frac{1}{16} [F_{X X X} + F_{X Y Y} + G_{X X Y} + G_{Y Y Y} + i (G_{X X X} + G_{X Y Y} - F_{X X Y} \\ - F_{Y Y Y} {)] |}_{(0, 0)} . \end{matrix}

We now come to the following conclusion as a result of the analysis derived above.

Theorem 15.

Suppose that the positive equilibrium point

Q_{1} (u_{1}, v_{1})

of system (1.6) exists. Let the parameters

(a, b, c, d, g, h, k, α, ρ) \in {\hat{N S}}_{Q_{1}}

and

Δ_{0}

and

ρ_{0}

be defined as in (5.12). If the parameter ρ varies in a vicinity of

ρ_{0}

(correspondingly, Δ varies around

Δ_{0}

) and

L \neq 0

, then system (1.6) undergoes a Neimark–Sacker bifurcation at the equilibrium point

Q_{1} (u_{1}, v_{1})

. Moreover, if

L < (>) 0

, a stable (an unstable) smooth closed invariant curve can be bifurcated out and the bifurcation is supercritical (subcritical).

5.2.2. Period-Doubling Bifurcation at the Fixed Point $Q_{1} (u_{1}, v_{1})$

From Case 3 in the proof of Theorem 13 for the stability of the positive equilibrium point

Q_{1} (u_{1}, v_{1})

, one can see that the dimension numbers for the stable manifold and unstable manifold of system (1.6) at the equilibrium point

Q_{1} (u_{1}, v_{1})

change when

Δ

varies in the vicinity of

Δ_{0}

(correspondingly,

ρ

varies in the vicinity of

ρ_{0}

) for

Ω_{2} (1) < - 2 \sqrt{Ω_{1} (1)}

, where

Δ_{0} = \frac{- Ω_{2} (1) \pm \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)}, ρ_{0} = {(Γ (1 + α) Δ_{0})}^{\frac{1}{α}} .

(5.16)

Hence, a bifurcation may occur. Noting that

λ_{1} = - 1

and

| λ_{2} | \neq 1

for

Δ = Δ_{0}

, we show that this bifurcation is a period-doubling one. Let

{\hat{P D}}_{Q_{1}} = \{(a, b, c, d, g, h, k, α, ρ) \in R_{+}^{7} | Ω_{2} (1) < - 2 \sqrt{Ω_{1} (1)}\} .

To analyze the period-doubling bifurcation of system (1.6) at the fixed point

Q_{1} (u_{1}, v_{1})

, it suffices for us to consider

Δ_{0} = \frac{- Ω_{2} (1) + \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)}

. The proof for the case

Δ_{0} = \frac{- Ω_{2} (1) - \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)}

is completely similar and will be omitted here. Now, proceed in the following way.

Let

x_{n} = u_{n} - u_{1}

y_{n} = v_{n} - v_{1}

, which transforms the fixed point

Q_{1} (u_{1}, v_{1})

to the origin

O (0, 0)

. Consider

ρ^{*}

as a small perturbation of

ρ

, i.e.,

ρ^{*} = ρ - ρ_{0},

with

| ρ^{*} | ≪ 1

. After the perturbation, system (1.6) takes the following form:

\begin{matrix} \{\begin{matrix} x_{n + 1} = x_{n} + \frac{{(ρ + ρ^{*})}^{α}}{Γ (α + 1)} (a (x_{n} + u_{1}) (1 - \frac{(x_{n} + u_{1})}{k}) - \frac{b (x_{n} + u_{1}) (y_{n} + v_{1})}{h + {(x_{n} + u_{1})}^{2}}), \\ y_{n + 1} = y_{n} + \frac{{(ρ + ρ^{*})}^{α}}{Γ (α + 1)} (\frac{c b (x_{n} + u_{1}) (y_{n} + v_{1})}{h + {(x_{n} + u_{1})}^{2}} - d (y_{n} + v_{1}) - g (x_{n} + u_{1}) (y_{n} + v_{1})) . \end{matrix} \end{matrix}

(5.17)

Set

ρ_{n + 1}^{*} = ρ_{n}^{*} = ρ^{*}

, then (5.17) can be seen as

\begin{matrix} \{\begin{matrix} x_{n + 1} & = x_{n} + \frac{{(ρ + ρ^{*})}^{α}}{Γ (α + 1)} (a (x_{n} + u_{1}) (1 - \frac{(x_{n} + u_{1})}{k}) - \frac{b (x_{n} + u_{1}) (y_{n} + v_{1})}{h + {(x_{n} + u_{1})}^{2}}), \\ y_{n + 1} & = y_{n} + \frac{{(ρ + ρ^{*})}^{α}}{Γ (α + 1)} (\frac{c b (x_{n} + u_{1}) (y_{n} + v_{1})}{h + {(x_{n} + u_{1})}^{2}} - d (y_{n} + v_{1}) - g (x_{n} + u_{1}) (y_{n} + v_{1})) \\ ρ_{n + 1}^{*} & = ρ_{n}^{*}, \end{matrix} \end{matrix}

(5.18)

Taylor expanding system (5.18) at

(x_{n}, v_{n}, ρ_{n}^{*}) = (0, 0, 0)

results in

\begin{matrix} (\begin{matrix} x_{n} \\ y_{n} \\ ρ_{n}^{*} \end{matrix}) \to (\begin{matrix} ϵ_{100} & ϵ_{010} & 0 \\ ζ_{100} & ζ_{010} & 0 \\ 0 & 0 & 1 \end{matrix}) (\begin{matrix} x_{n} \\ y_{n} \\ ρ^{*} \end{matrix}) + (\begin{matrix} M (x_{n}, y_{n}, ρ_{n}^{*}) + o (ρ_{3}^{3}) \\ N (x_{n}, y_{n}, ρ_{n}^{*}) + o (ρ_{3}^{3}) \\ 0 \end{matrix}), \end{matrix}

(5.19)

where

ρ_{3} = \sqrt{x_{n}^{2} + y_{n}^{2} + {ρ_{n}^{*}}^{2}}

\begin{matrix} M (x_{n}, y_{n}, ρ_{n}^{*}) = & ϵ_{200} x_{n}^{2} + ϵ_{020} y_{n}^{2} + ϵ_{002} {ρ_{n}^{*}}^{2} + ϵ_{110} x_{n} y_{n} + ϵ_{101} x_{n} ρ^{*} + ϵ_{011} y_{n} ρ_{n}^{*} \\ + ϵ_{300} x_{n}^{3} + ϵ_{030} y_{n}^{3} + ϵ_{003} {ρ_{n}^{*}}^{3} + ϵ_{210} x_{n}^{2} y_{n} + ϵ_{120} x_{n} y_{n}^{2} \\ + ϵ_{021} y_{n}^{2} ρ_{n}^{*} + ϵ_{012} y_{n} {ρ_{n}^{*}}^{2} + ϵ_{201} x_{n}^{2} ρ_{n}^{*} + ϵ_{102} x_{n} {ρ_{n}^{*}}^{2} + ϵ_{111} x_{n} y_{n} ρ_{n}^{*}, \\ N (x_{n}, y_{n}, ρ_{n}^{*}) = & ζ_{200} x_{n}^{2} + ζ_{020} y_{n}^{2} + ζ_{002} {ρ_{n}^{*}}^{2} + ζ_{110} x_{n} y_{n} + ζ_{101} x_{n} ρ^{*} + ζ_{011} y_{n} ρ_{n}^{*} \\ + ζ_{300} x_{n}^{3} + ζ_{030} y_{n}^{3} + ζ_{003} {ρ_{n}^{*}}^{3} + ζ_{210} x_{n}^{2} y_{n} + ζ_{120} x_{n} y_{n}^{2} \\ + ζ_{021} y_{n}^{2} ρ_{n}^{*} + ζ_{012} y_{n} {ρ_{n}^{*}}^{2} + ζ_{201} x_{n}^{2} ρ_{n}^{*} + ζ_{102} x_{n} {ρ_{n}^{*}}^{2} + ζ_{111} x_{n} y_{n} ρ_{n}^{*}, \end{matrix}

\begin{matrix} ϵ_{100} & = 1 + \frac{ρ^{α}}{Γ (α + 1)} (a (1 - \frac{2 u_{1}}{k}) - b v_{1} (\frac{h - u_{1}^{2}}{{(u_{1}^{2} + h)}^{2}})), ϵ_{010} = - \frac{ρ^{α} b u_{1}}{Γ (α + 1) (u_{1}^{2} + h)}, \\ ϵ_{200} & = \frac{ρ^{α}}{Γ (α + 1)} ((3 h - u_{1}^{2}) \frac{b u_{1} v_{1}}{{(u_{1}^{2} + h)}^{3}} - \frac{a}{k}), ϵ_{020} = 0, \\ ϵ_{002} & = \frac{α (1 - α ρ^{α - 2} u_{1})}{2 Γ (α + 1)} (a (\frac{u_{1}}{k} - 1) + \frac{b v_{1}}{u_{1}^{2} + h}), ϵ_{110} = \frac{ρ^{α} b}{Γ (α + 1)} (\frac{u_{1}^{2} - h}{{(u_{1}^{2} + h)}^{2}}), \end{matrix}

\begin{matrix} ϵ_{101} & = \frac{α ρ^{α - 1}}{Γ (α + 1)} (\frac{b v_{1}}{u_{1}^{2} + h} (\frac{2 u_{1}^{2} v_{1}}{u_{1}^{2}} - 1) + a (1 - \frac{2 u_{1}}{k})), ϵ_{011} = - \frac{α ρ^{α - 1} b u_{1}}{Γ (α + 1) (u_{1}^{2} + h)}, \\ ϵ_{300} & = \frac{ρ^{α} b v_{1}}{Γ (α + 1) {(u_{1}^{2} + h)}^{3}} (h + \frac{2 u_{1}^{2} (3 u_{1}^{2} - h)}{u_{1}^{2} + h} - 5 u_{1}^{2}), ϵ_{030} = 0, \\ ϵ_{003} & = \frac{α (3 α - α^{2} - 2) ρ^{α - 3}}{6 Γ (α + 1)} (a u_{1} (\frac{u_{1}}{k} - 1) + \frac{b u_{1} v_{1}}{u_{1}^{2} + h}), \\ ϵ_{210} & = \frac{ρ^{α} b u_{1}}{Γ (α + 1) {(u_{1}^{2} + h)}^{2}} (2 - \frac{3 (u_{1}^{2} - h)}{u_{1}^{2} + h}), ϵ_{120} = 0, ϵ_{021} = 0, \\ ϵ_{012} & = \frac{α (1 - α) ρ^{α - 2} b u_{1}}{2 Γ (α + 1) (u_{1}^{2} + h)}, ϵ_{201} = \frac{α ρ^{α - 1}}{Γ (α + 1)} ((3 h - u_{1}^{2}) \frac{b u_{1} v_{1}}{{(u_{1}^{2} + h)}^{3}} - \frac{a}{k}), \\ ϵ_{102} & = \frac{α (1 - α) ρ^{α - 2} b u_{1}}{2 Γ (α + 1) (u_{1}^{2} + h)}, ϵ_{111} = \frac{α ρ^{α - 1} b (u_{1}^{2} - h)}{Γ (α + 1) {(u_{1}^{2} + h)}^{2}}, \\ ζ_{100} = & \frac{ρ^{α} v_{1}}{Γ (α + 1)} (\frac{b c}{u_{1}^{2} + h} (1 - \frac{2 u_{1}^{2}}{u_{1}^{2} + h}) - g), ζ_{010} = 1 - \frac{ρ^{α}}{Γ (α + 1)} (d + g u_{1} - \frac{b c u_{1}}{u_{1}^{2} + h}), \\ ζ_{200} = & \frac{ρ^{α} b c u_{1} v_{1}}{Γ (α + 1) {(u_{1}^{2} + h)}^{2}}, ζ_{020} = 0, ζ_{002} = \frac{α (1 - α) ρ^{α - 2} v 1}{2 Γ (α + 1)} (d - g u_{1} - \frac{b c u_{1}}{u_{1}^{2} + h}), \\ ζ_{110} = & \frac{ρ^{α}}{Γ (α + 1)} (\frac{b c}{u_{1}^{2} + h} (1 - \frac{2 u_{1}^{2}}{u_{1}^{2} + h}) - g), ζ_{101} = \frac{α ρ^{α - 1} v_{1}}{Γ (α + 1)} (\frac{b c}{u_{1}^{2} + h} (1 - \frac{2 u_{1}}{u_{1}^{2} + h}) - g), \\ ζ_{011} = & \frac{α ρ^{α - 1}}{Γ (α + 1)} (\frac{b c u_{1}}{u_{1}^{2} + h} - d - g u_{1}), ζ_{300} = \frac{ρ^{α} b c v_{1}}{Γ (α + 1) {(u_{1}^{2} + h)}^{2}} (\frac{4 u_{1}^{2}}{u_{1}^{2} + h} + 4 u_{1}^{2} (1 - 2 u_{1}^{2}) - 1), \\ ζ_{030} = & 0, ζ_{003} = \frac{α (3 α - α^{2} - 2) ρ^{α - 3}}{6 Γ (α + 1)} (d + g u_{1} - \frac{b c u_{1}}{u_{1}^{2} + h}), \\ ζ_{210} = & \frac{ρ^{α} b c u_{1}}{Γ (α + 1) {(u_{1}^{2} + h)}^{2}} (\frac{4 u_{1}^{2}}{u_{1}^{2} + h} - 3), ζ_{120} = 0, ζ_{021} = 0, \\ ζ_{012} = & \frac{α (1 - α ρ^{α - 2})}{2 Γ (α + 1)} (d + g u_{1} - \frac{b c u_{1}}{u_{1}^{2} + h}), ζ_{201} = \frac{α ρ^{α - 1} b c u_{1} v_{1}}{Γ (α + 1)} (\frac{4 u_{1}^{2}}{u_{1}^{2} + h} - 1), \\ ζ_{102} = & \frac{α (1 - α) ρ^{α - 2} v_{1}}{2 Γ (α + 1)} (g - \frac{b c}{u_{1}^{2} + h} (1 - \frac{2 u_{1}}{u_{1}^{2} + h})), \\ ζ_{111} = & \frac{α ρ^{α - 1}}{Γ (α + 1)} (\frac{b c}{u_{1}^{2} + h} (1 - \frac{2 u_{1}^{2}}{u_{1}^{2} + h} - g)) . \end{matrix}

Take

T = (\begin{matrix} ϵ_{010} & ϵ_{010} \\ - 1 - ϵ_{100} & λ_{2} - ϵ_{100} \end{matrix}),

which is invertible. Now, using the transformation

(\begin{matrix} x_{n} \\ y_{n} \end{matrix}) = T (\begin{matrix} \bar{u} \\ \bar{v} \end{matrix}),

system (5.19) becomes

\begin{matrix} \{\begin{matrix} {\bar{u}}_{n + 1} = - {\bar{u}}_{n} + M (x_{n}, y_{n}, ρ_{n}^{*}), \\ {\bar{v}}_{n + 1} = λ_{2} {\bar{v}}_{n} + N (x_{n}, y_{n}, ρ_{n}^{*}) . \end{matrix} \end{matrix}

(5.20)

System (5.20) has a center manifold

W^{c} (0, 0, 0)

(0, 0)

in the neighborhood of

ρ^{*} = 0

, which can be deduced using the center manifold theorem and is essentially expressed as follows:

W^{c} (0, 0, 0) = \{({\bar{u}}_{n}, {\bar{v}}_{n}, ρ^{*}) \in R^{3} : {\bar{v}}_{n} = η_{1} {\bar{u}}_{n}^{2} + η_{2} {\bar{u}}_{n} ρ^{*} + o ({(| {\bar{u}}_{n} | + | ρ^{*} |)}^{2})\}

where

\begin{matrix} η_{1} = & \frac{ϵ_{010} ((1 + ϵ_{100}) ϵ_{200} + ϵ_{010} ζ_{200}) + ζ_{020} {(1 + ϵ_{100})}^{2} - (1 + ϵ_{100}) (ϵ_{110} (1 + ϵ_{100}) + ϵ_{010} ζ_{110})}{1 - λ_{2}^{2}}, \\ η_{2} = & \frac{(1 + ϵ_{100}) (ϵ_{011} (1 + ϵ_{100}) + ϵ_{010} ζ_{011})}{ϵ_{010} {(1 + λ_{2})}^{2}} - \frac{(1 + ϵ_{100}) (ϵ_{101} + ϵ_{010} ζ_{101})}{{(1 + λ_{2})}^{2}} . \end{matrix}

So, system (5.20) restrained on the center manifold

W^{c} (0, 0, 0)

has the following form:

{\bar{u}}_{n + 1} = - {\bar{u}}_{n} + θ_{1} {\bar{u}}_{n}^{2} + θ_{2} {\bar{u}}_{n} ρ^{*} + θ_{3} {\bar{u}}_{n}^{2} ρ^{*} + θ_{4} {\bar{u}}_{n} {ρ^{*}}^{2} + θ_{5} {\bar{u}}_{n}^{3} + o ({(| {\bar{u}}_{n} | + | ρ^{*} |)}^{3}) \equiv Z ({\bar{u}}_{n}, ρ^{*}),

where

\begin{matrix} θ_{1} = & \frac{η_{2} ((λ_{2} - η_{1}) ϵ_{200} - η_{2} ζ_{200})}{1 + λ_{2}} - \frac{ζ_{020} {(1 + η_{1})}^{2}}{1 + λ_{2}} - \frac{(1 + η_{1}) ((λ_{2} - η_{1}) ϵ_{110} - η_{2} ζ_{110})}{1 + λ_{2}}, \\ θ_{2} = & \frac{(λ_{2} - η_{1}) ϵ_{101} - η_{2} ζ_{101}}{1 + λ_{2}} - \frac{(1 + η_{1}) (λ_{2} - η_{1}) ϵ_{011} - η_{2} ζ_{011}}{η_{2} (1 + λ_{2})}, \\ θ_{3} = & \frac{(λ_{2} - ϵ_{100}) η_{1} ϵ_{101} - ϵ_{010} ζ_{101}}{1 + λ_{2}} + \frac{((λ_{2} - ϵ_{100}) ϵ_{011} - η_{2} ζ_{011}) (λ_{2} - ϵ_{100}) η_{1}}{ϵ_{010} (1 + λ_{2})} \\ - \frac{(1 - ϵ_{100}) ((λ_{2} - ϵ_{100}) ϵ_{111} - ϵ_{010} ζ_{111})}{1 + λ_{2}} + \frac{ϵ_{010} ((λ_{2} - ϵ_{100}) ϵ_{201} - ϵ_{010} ζ_{201})}{1 + λ_{2}} \\ - \frac{ζ_{021} {(1 - ϵ_{100})}^{2}}{1 + λ_{2}} + \frac{2 ϵ_{100} η_{2} ((λ_{2} - ϵ_{100}) ϵ_{200} - ϵ_{010} ζ_{200})}{1 + λ_{2}} \\ - \frac{2 ζ_{020} η_{2} (1 + ϵ_{100}) (λ_{2} - ϵ_{100})}{1 + λ_{2}} + \frac{η_{2} ((λ_{2} - ϵ_{100}) ϵ_{110} - ϵ_{010} ζ_{110}) (λ_{2} - 1 - 2 ϵ_{100})}{1 + λ_{2}}, \\ θ_{4} = & \frac{η_{2} ((λ_{2} - ϵ_{100}) ϵ_{101} - ϵ_{010} ζ_{101})}{1 + λ_{2}} + \frac{((λ_{2} - ϵ_{100}) ϵ_{011} - ϵ_{010} ζ_{011}) (λ_{2} - ϵ_{100}) η_{2}}{ϵ_{010} (1 + λ_{2})} \\ + \frac{2 ϵ_{010} η_{2} ((λ_{2} - ϵ_{100}) ϵ_{200} - ϵ_{010} ζ_{200})}{1 + λ_{2}} + \frac{2 ζ_{020} η_{2} (1 + ϵ_{100}) (λ_{2} - ϵ_{100})}{1 + λ_{2}} \\ + \frac{η_{2} (λ_{2} - 1 - 2 ϵ_{100}) ((λ_{2} - ϵ_{100}) ϵ_{110} - ϵ_{010} ζ_{110})}{1 + λ_{2}}, \\ θ_{5} = & \frac{2 ϵ_{010} η_{1} ((λ_{2} - ϵ_{100}) ϵ_{200} - ϵ_{010} ζ_{200})}{1 + λ_{2}} + \frac{η_{1} (λ_{2} - 1 - 2 ϵ_{100}) ((λ_{2} - ϵ_{100}) ϵ_{200} - ϵ_{010} ζ_{200})}{2 + λ_{2}} \\ + \frac{2 ζ_{020} η_{1} (λ_{2} - ϵ_{100}) (1 + ϵ_{100})}{1 + λ_{2}} + \frac{η_{2}^{2} ((λ_{2} - ϵ_{100}) ϵ_{300} - ϵ_{010} ζ_{300})}{1 + λ_{2}} \\ - \frac{η_{2} (1 + ϵ_{100}) ((λ_{2} - ϵ_{100}) ϵ_{210} - ϵ_{010} ζ_{210})}{1 + λ_{2}} . \end{matrix}

In order for the period-doubling bifurcation to occur, the two determinating quantities

ξ_{1}

and

ξ_{2}

must both be nonzero, where

ξ_{1} = (\frac{\partial^{2} Z}{\partial \bar{u} \partial ρ^{*}} + \frac{1}{2} \frac{\partial Z}{\partial ρ^{*}} \frac{\partial^{2} Z}{\partial {\bar{u}}^{2}}) {|_{(0, 0)}, ξ_{2} = (\frac{1}{6} \frac{\partial^{3} Z}{\partial {\bar{u}}^{3}} + {(\frac{1}{2} \frac{\partial^{2} Z}{\partial {\bar{u}}^{2}})}^{2}) |}_{(0, 0)} .

Finally, the outcome of the analysis above is as follows.

Theorem 16.

Suppose that the positive equilibrium point

Q_{1} (u_{1}, v_{1})

of system (1.6) exists. Let the parameters

(a, b, c, d, g, h, k, α, ρ) \in {\hat{P D}}_{Q_{1}}

and

Δ_{0}

and

ρ_{0}

be defined as in (5.16). If the parameter ρ varies in a neighbourhood of

ρ_{0}

(correspondingly, Δ varies around

Δ_{0}

) and

ξ_{1} ξ_{2} \neq 0

, then system (1.6) undergoes a period-doubling bifurcation at the equilibrium point

Q_{1} (u_{1}, v_{1})

. Furthermore, for

ξ_{2} > (<) 0

, the period-two orbit that bifurcates from

Q_{1} (u_{1}, v_{1})

is stable (unstable).

6. Numerical Simulation

In this section, we perform numerical simulations of the dynamical behavior of systems (1.3) and (1.6) using Matlab, aiming to provide readers with a more intuitive understanding to the dynamics of systems (1.3) and (1.6).

In Figure 1, the parameter values in system (1.3) are

c = 0.015

k = 0.3

a = 0.1

b = 0.9

d = 0.3

h = 0.2

g = 0.1

and

α = 0.98

. Figure 1a displays the trajectories of system (1.3) starting from different points. Although it can be observed that system (1.3) exhibits a saddle at the origin, it is not entirely clear. To provide a more clear representation of the behavior of system (1.3) at the origin, we constructed streamline plots, depicted in Figure 1b. From Figure 1b, it is evident that system (1.3) possesses a saddle at the origin.

In Figure 2a,b, the parameter values of system (1.3) are

c = 0.9

k = 10

b = 0.2

a = 3

h = 0.9

d = 0.2

g = 0.2

and

α = 0.98

, which satisfy

d > \frac{c b k}{h + k^{2}} - g k

. Figure 2a shows that the behavior of system (1.3), regardless of whether it starts from the point

(30, 10)

(30, 12)

(30, 14)

, will eventually converge to the point (10, 0). Figure 2b demonstrates how the populations of prey and predator change over time when starting from the point

(30, 10)

. We can observe that as time increases, the population of prey tends to 10, while the predator becomes extinct. In Figure 2c, the parameter values of system (1.3) are

c = 0.6

k = 1

b = 0.5

a = 0.1

h = 0.6

d = 0.1

g = 0.1

and

α = 0.98

, which satisfy

d < \frac{c b k}{h + k^{2}} - g k

. We can clearly see that system (1.3) exhibits a saddle at the boundary equilibrium point

(1, 0)

For the positive equilibrium point of system (1.3), we are interested in its bifurcation behavior. In Figure 3, the parameter values of system (1.3) are

c = 0.6

b = 0.5

a = 0.1

h = 0.6

d = 0.1

g = 0.1

and

α = 0.98

. Figure 3a,b shows that the positive equilibrium point

Q_{1} (u_{1}, v_{1})

is stable and unstable when

k = 1

and

k = 10

, respectively. Furthermore, we can see from Figure 3b that when k crosses the critical value, a stable limit cycle emerges, indicating the occurrence of a supercritical Hopf bifurcation in system (1.3).

With values of

α = 0.75, b = 0.9, a = 0.2, k = 3, h = 5, c = 0.7, d = 0.2

and

g = 0.1

, Figure 4 is the bifurcation diagram of system (1.6) starting from the point

(u_{0}, v_{0}) = (0.6, 0.1)

, and we can clearly observe that system (1.6) undergoes a period-doubling bifurcation at the critical value. With values of

α = 0.75, b = 0.9, a = 0.2, k = 3, h = 5, c = 1, d = 0.2, g = 0.1

, Figure 5 is the bifurcation diagram of system (1.6) starting from the point

(u_{0}, v_{0}) = (0.5, 0.1)

, and it is clear that system (1.6) undergoes a Neimark–Sacker bifurcation at the critical value.

Figure 6 depicts the phase diagram of system (1.6) starting from the point

(u_{0}, v_{0}) = (0.5888, 4.324)

with parameters

α = 0.75, b = 0.9, a = 0.2, k = 3, h = 5,

c = 1,

d = 0.2, g = 0.1

. We can observe that as

ρ

increases, the equilibrium point gradually transitions from a stable focus to an unstable focus, and a stable limit cycle emerges.

7. Conclusions

In this paper, we propose a fractional-order predator–prey model with a Holling type IV functional response and anti-predator behavior. According to the discrete and continuous versions, from two different perspectives we analyzed their dynamical behavior in detail, including the feasibility, existence and stability of equilibrium points and the possibility of local bifurcations. Our main aim is to provide readers with a better understanding of the dynamics of the system. As there is currently a lack of effective ways and methods to study the dynamics of fractional-order differential systems, in this paper, we propose an effective way to consider this problem from different angles—both continuous and discrete. This is the novelty of this paper. Indeed, we find that there exist some differences in the dynamics of the system between the continuous version and the discrete version. Numerical simulations also illustrate corresponding theoretical results. By analyzing the dynamical behavior of systems (1.3) and (1.6), respectively, we can deduce the following conclusions:

(1) By analyzing the stability of the equilibrium point

Q_{0} (0, 0)

and conducting numerical simulations, we can determine that the equilibrium point

Q_{0} (0, 0)

is a saddle point. This implies that under any conditions existing in nature, the simultaneous extinction of predator and prey does not occur.

(2) Through the study of the dynamical behavior of the boundary equilibrium point

Q_{k} (k, 0)

and numerical simulations, we have found that when d is large, it leads to the extinction of predator. In this case, the prey population tends towards a stable density. On the other hand, when d is small, the extinction of the predator does not occur, and the prey population tends to a stable state. This indicates that when a detrimental condition for the survival of predator and prey arises in nature, the predator may tend towards extinction, while the prey population, although it may decrease, does not tend towards extinction. Instead, it stabilizes at a certain level.

(3) Based on the analysis and numerical simulations of the positive equilibrium point

Q_{1} (u_{1}, v_{1})

, we can draw the following conclusions: When the parameter k exceeds a critical value, the system exhibits a stable limit cycle. This implies that the interaction between predator and prey leads to periodic oscillations. The presence of this limit cycle indicates that the system exhibits rich dynamic behavior, and under specific conditions, the populations of predator and prey undergo periodic fluctuations. Therefore, we can achieve a steady coexistence state and eliminate the limit cycle by reducing the environmental carrying capacity to the prey.

(4) Through the bifurcation analysis used in this paper, we find that the analysis methods for bifurcation problems are applicable to other types of fractional differential systems. It is well known that the current analysis methods for bifurcation problems in fractional-order dynamic systems are not well developed. Thus, in order to better understand the dynamics of this system, we discretize the fractional-order system to study its dynamics from a different angle. A richer set of dynamical properties is obtained, indicating that investigations after discretizing this system are indeed more valuable and helpful to understanding the properties of this system.

Author Contributions

Methodology, X.L.; Software, B.W.; Formal analysis, B.W.; Resources, B.W.; Data curation, B.W.; Writing—original draft, B.W.; Writing—review & editing, X.L.; Supervision, X.L.; Project administration, X.L.; Funding acquisition, X.L. All authors contributed equally and significantly to writing this paper. All authors have read and agreed to the final manuscript.

Funding

This work was partly supported by the National Natural Science Foundation of China (61473340), the Distinguished Professor Foundation of Qianjiang Scholar in Zhejiang Province (F703108L02), and the Natural Science Foundation of Zhejiang University of Science and Technology (F701108G14).

Data Availability Statement

There are no applicable data associated with this manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

Berryman, A.A. The orgins and evolution of predator-prey theory. Ecology 1992, 73, 1530–1535. [Google Scholar] [CrossRef]
Wang, W.; Chen, L. A Predator-Prey System with Stage-Structure for Predator. Comput. Math. Appl. 1997, 33, 83–91. [Google Scholar] [CrossRef]
Ba, Z.; Li, X. Period-doubling bifurcation and Neimark-Sacker bifurcation of a discrete predator-prey model with Allee effect and cannibalism. Electron. Res. Arch. 2023, 31, 1405–1438. [Google Scholar] [CrossRef]
Liu, Y.; Li, X. Dynamics of a discrete predator-prey model with Holling-II functional response. Int. J. Biomath. 2021, 14, 2150068. [Google Scholar] [CrossRef]
Qi, H.; Meng, X. Threshold behavior of a stochastic predator–prey system with prey refuge and fear effect. Appl. Math. Lett. 2021, 113, 106846. [Google Scholar] [CrossRef]
Blasius, B.; Rudolf, L.; Weithoff, G.; Gaedke, U.; Fussmann, G.F. Long-term cyclic persistence in an experimental predator–prey system. Nature 2020, 577, 226–230. [Google Scholar] [CrossRef]
Mukherjee, D. Role of fear in predator–prey system with intraspecific competition. Math. Comput. Simul. 2020, 177, 263–275. [Google Scholar] [CrossRef]
Qiu, H.; Guo, S.; Li, S. Stability and Bifurcation in a Predator–Prey System with Prey-Taxis. Int. J. Bifurc. Chaos 2020, 30, 2050022. [Google Scholar] [CrossRef]
Perc, M.; Grigolini, P. Collective behavior and evolutionary games–An introduction. Chaos Solitons Fractals 2013, 56, 1–5. [Google Scholar] [CrossRef]
Lima, S.L. Stress and Decision Making under the Risk of Predation: Recent Developments from Behavioral, Reproductive, and Ecological Perspectives. Adv. Study Behav. 1998, 27, 215–290. [Google Scholar]
Tollrian, R. Predator-Induced Morphological Defenses: Costs, Life History Shifts, and Maternal Effects in Daphnia Pulex. Ecology 1995, 76, 1691–1705. [Google Scholar] [CrossRef]
Choh, Y.; Lgnacio, M.; Sabelis, M.W.; Janssen, A. Predator-prey role reversals, juvenile experience and adult antipredator behavior. Lab. Appl. Entomol. 2012, 2, 728. [Google Scholar]
Janssen, A.; Faraji, F.; Hammen, T.V.D.; Magalhaes, S.; Sabelis, M.W. Interspecific infanticide deters predators. Ecol. Lett. 2002, 5, 490–494. [Google Scholar] [CrossRef]
Saito, Y. Prey kills predator: Counter-attack success of a spider mite against its specific phytoseiid predator. Exp. Appl. Acarol. 1986, 2, 47–62. [Google Scholar] [CrossRef]
Ives, A.; Dobson, A. Antipredator behavior and the population dynamics of simple predator-prey systems. Am. Nat. 1987, 130, 431–447. [Google Scholar] [CrossRef]
Tang, B.; Xiao, Y. Bifurcation analysis of a predator–prey model with anti-predator behavior. Chaos Solitons Fractals 2015, 70, 58–68. [Google Scholar] [CrossRef]
Sokol, W. Oxidation of an Inhibitory Substrate by Washed Celts (Oxidation of Phenol by Pseudornonas putida). Biotechnol. Bioeng. 1987, 30, 921–927. [Google Scholar] [CrossRef]
Liouville, J. Troisième mémoire sur le développement des fonctions ou parties de fonctions en séries dont les divers termes sont assujettis à satisfaire à une même équation différentielle du second ordre, contenant un paramètre variable. J. Math. Pures Appl. 1837, 2, 418–436. [Google Scholar]
Riesz, M. L’intégrale de Riemann–Liouville et le problème de Cauchy pour l’équation des ondes. Bull. Soc. Math. Fr. 1939, 67, 153–170. [Google Scholar] [CrossRef]
Caputo, M. Linear Models of Dissipation whose Q is almost Frequency Independent-II. Geophys. J. R. Astron. Soc. 1976, 13, 529–539. [Google Scholar] [CrossRef]
Debnath, S.; Majumdar, P.; Sarkar, S.; Ghosh, U. Memory effect on prey–predator dynamics: Exploring the role of fear effect, additional food and anti-predator behaviour of prey. J. Comput. Sci. 2023, 66, 101929. [Google Scholar] [CrossRef]
Singh, A.; Sharma, V.S. Bifurcations and chaos control in a discrete-time prey–predator model with Holling type-II functional response and prey refuge. J. Comput. Appl. Math. 2023, 418, 114666. [Google Scholar] [CrossRef]
Kumar, S.; Kumar, R.; Cattani, C.; Samet, B. Chaotic behaviour of fractional predator-prey dynamical system. Chaos Solitons Fractals 2020, 135, 109811. [Google Scholar] [CrossRef]
Das, M.; Samanta, G. A delayed fractional order food chain model with fear effect and prey refuge. Math. Comput. Simul. 2020, 178, 218–245. [Google Scholar] [CrossRef]
Barman, D.; Roy, J.; Alrabaiah, H.; Panja, P.; Mondal, S.P.; Alam, S. Impact of predator incited fear and prey refuge in a fractional order prey predator model. Chaos Solitons Fractals 2021, 142, 110420. [Google Scholar] [CrossRef]
Das, M.; Samanta, G. A prey-predator fractional order model with fear effect and group defense. Int. J. Dyn. Control 2021, 9, 334–349. [Google Scholar] [CrossRef]
Yousef, F.; Yousef, A.; Maji, C. Effects of fear in a fractional-order predator-prey system with predator density-dependent prey mortality. Chaos Solitons Fractals 2021, 145, 110711. [Google Scholar] [CrossRef]
Zhang, N.; Kao, Y.; Xie, B. Impact of fear effect and prey refuge on a fractional order prey–predator system with Beddington–DeAngelis functional response. Chaos 2022, 32, 043125. [Google Scholar] [CrossRef]
Majumdar, P.; Mondal, B.; Debnath, S.; Ghosh, U. Controlling of periodicity and chaos in a three dimensional prey predator model introducing the memory effect. Chaos Solitons Fractals 2022, 164, 112585. [Google Scholar] [CrossRef]
Xie, B.; Zhang, Z. Impact of Allee and fear effects in a fractional order prey–predator system incorporating prey refuge. Chaos 2023, 33, 013131. [Google Scholar] [CrossRef]
Salmana, S.M.; Yousef, A.M.; Elsadanyc, A.A. Stability, bifurcation analysis and chaos control of a discrete predator-prey system with square root functional response. Chaos Solitons Fractals 2016, 93, 20–31. [Google Scholar] [CrossRef]
Yousef, A.M.; Rida, S.Z.; Ali, H.M.; Zaki, A.S. Stability, co-dimension two bifurcations and chaos control of a host-parasitoid model with mutual interference. Chaos Solitons Fractals 2023, 166, 112923. [Google Scholar] [CrossRef]
Yousef, A.M.; Rida, S.Z.; Gouda, Y.G.; Zaki, A.S. Dynamical Behaviors of a Fractional-Order Predator–Prey Model with Holling Type IV Functional Response and Its Discretization. Int. J. Nonlinear Sci. Numer. Simul. 2019, 20, 125–136. [Google Scholar] [CrossRef]
Singh, A.; Elsadany, A.A.; Elsonbaty, A. Complex dynamics of a discrete fractional-order Leslie-Gower predator-prey model. Math. Methods Appl. Sci. 2019, 42, 3992–4007. [Google Scholar] [CrossRef]
Elsonbaty, A.; Elsadany, A.A. On discrete fractional order Lotka Volterra model based on the Caputo diference discrete operator. Math. Sci. 2023, 17, 67–79. [Google Scholar] [CrossRef]
Uddin, M.J.; Rana, S.M.S.; Işık, S.; Kangalgil, F. On the qualitative study of a discrete fractional order prey–predator model with the effects of harvesting on predator population. Chaos Solitons Fractals 2023, 175, 113932. [Google Scholar] [CrossRef]
Kai, D. The Analysis of Fractional Differential Equations (An Application-Oriented Exposition Using Differential Operators of Caputo Type); Springer Scienc & Business Media: Berlin/Heidelberg, Germany, 2010; Volume 2004, pp. 237–244. [Google Scholar]
Li, W.; Li, X.Y. Neimark–Sacker bifurcation of a semi-discrete hematopoiesis model. J. Appl. Anal. Comput. 2018, 8, 1679–1693. [Google Scholar]
Cieck, M.; Yakar, C.; Oğur, B. Stability, Boundedness, and Lagrange Stability of Fractional Differential Equations with Initial Time Difference. Sci. World J. 2014, 2014, 939027. [Google Scholar]
Odibat, Z.; Shawagfeh, N. Generalized Taylor’s formula. Appl. Math. Comput. 2007, 186, 286–293. [Google Scholar] [CrossRef]
Wang, C.; Li, X.Y. Stability and Neimark–Sacker bifurcation of a semi-discrete population model. J. Appl. Anal. Comput. 2014, 4, 419–435. [Google Scholar]
Liang, S.; Wu, R.; Chen, L. Laplace transform of fractional order differential equations. Electron. J. Differ. Equ. 2015, 2015, 1–15. [Google Scholar]
Kexue, L.; Jigen, P. Laplace transform and fractional differential equations. Appl. Math. Lett. 2011, 24, 2019–2023. [Google Scholar] [CrossRef]
Petráš, I. Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation; Springer: Berlin/Heidelberg, Germany, 2011. [Google Scholar]
Deshpande, A.; Daftardar-Gejji, V.; Sukale, Y.V. On Hopf bifurcation in fractional dynamical systems. Chaos Solitons Fractals 2017, 98, 189–198. [Google Scholar] [CrossRef]

$Fractalfract 07 00722 g001a$ $Fractalfract 07 00722 g001b$

Figure 1. (a) shows the trajectories of system (1.3) starting from points

(0.8, 0.4)

(0.7, 0.3)

and

(0.6, 0.35)

, respectively; (b) represents the streamline plots of system (1.3) at the origin.

Figure 1. (a) shows the trajectories of system (1.3) starting from points

(0.8, 0.4)

(0.7, 0.3)

and

(0.6, 0.35)

, respectively; (b) represents the streamline plots of system (1.3) at the origin.

$Fractalfract 07 00722 g001a$ $Fractalfract 07 00722 g001b$

$Fractalfract 07 00722 g002$

Figure 2. (a) The properties of system (1.3) at the boundary equilibrium point (10,0) for

d > \frac{c b k}{h + k^{2}} - g k

; (b) the quantities of prey and predators starting from the point

(30, 10)

over time; (c) the properties of system (1.3) at the boundary equilibrium point (1,0) for

d < \frac{c b k}{h + k^{2}} - g k

Figure 2. (a) The properties of system (1.3) at the boundary equilibrium point (10,0) for

d > \frac{c b k}{h + k^{2}} - g k

; (b) the quantities of prey and predators starting from the point

(30, 10)

over time; (c) the properties of system (1.3) at the boundary equilibrium point (1,0) for

d < \frac{c b k}{h + k^{2}} - g k

$Fractalfract 07 00722 g002$

$Fractalfract 07 00722 g003$

Figure 3. The existence of a supercritical Hopf bifurcation of system (1.3) with the parameter values

c = 0.6

b = 0.5

a = 0.1

h = 0.6

d = 0.1

g = 0.1

and

α = 0.98

for

k = 1

(a) and

k = 10

(b).

Figure 3. The existence of a supercritical Hopf bifurcation of system (1.3) with the parameter values

c = 0.6

b = 0.5

a = 0.1

h = 0.6

d = 0.1

g = 0.1

and

α = 0.98

for

k = 1

(a) and

k = 10

(b).

$Fractalfract 07 00722 g003$

$Fractalfract 07 00722 g004$

Figure 4. The existence of a period-doubling bifurcation of system (1.6) with the k taking values from 8 to

11.6

Figure 4. The existence of a period-doubling bifurcation of system (1.6) with the k taking values from 8 to

11.6

$Fractalfract 07 00722 g004$

$Fractalfract 07 00722 g005$

Figure 5. The existence of a Neimark–Sacker bifurcation of system (1.6) with the k taking values from

10.6

11.4

Figure 5. The existence of a Neimark–Sacker bifurcation of system (1.6) with the k taking values from

10.6

11.4

$Fractalfract 07 00722 g005$

$Fractalfract 07 00722 g006$

Figure 6. Phase portraits of system (1.6) with

α = 0.75, b = 0.9, a = 0.2, k = 3, h = 5, c = 1,

d = 0.2,

g = 0.1

and different

ρ

when the initial value

(x_{0}, y_{0}) = (0.5888, 0.4324)

Figure 6. Phase portraits of system (1.6) with

α = 0.75, b = 0.9, a = 0.2, k = 3, h = 5, c = 1,

d = 0.2,

g = 0.1

and different

ρ

when the initial value

(x_{0}, y_{0}) = (0.5888, 0.4324)

$Fractalfract 07 00722 g006$

Table 1. Biological meanings of all parameters in system (1.1).

Parameter	Meaning
x	Prey population density
y	Predator population density
$α > 0$	Natality of prey population
$k > 0$	Carrying capacity of the environment to prey
$β > 0$	Cost incurred by the prey as a result of anti-predator behavior
$γ > 0$	Effects resulting from the anti-predator behavior of prey
$\frac{q x}{1 + a x}, q > 0, a > 0$	Holling II functional response function
$b > 0$	Conversion rate of prey consumed by predator
$d > 0$	Death rate of predator population

Table 2. Biological meanings of parameters in systems (1.2) and (1.3).

Parameter	Meaning
u	Prey population density
v	Predator population density
$a > 0$	Natality of prey population
$k > 0$	Carrying capacity of the environment to prey
$b > 0$	Predator’s capture rate
$c > 0$	Conversion rate of prey into predator
$\frac{b u}{h + u^{2}}, h > 0$	Holling type IV functional response function
$d > 0$	Death rate of predator
$g > 0$	Mortality rate of predator due to the anti-predator effects of prey
$0 < α \leq 1$	Order of fractional-order derivative

Table 3. The stability of the equilibrium points of system (1.3).

Point	Conditions	Properties
$Q_{0} (0, 0)$		saddle
$Q_{k} (k, 0)$	$d > \frac{b c k}{h + k^{2}} - g k$	stable
	$d < \frac{b c k}{h + k^{2}} - g k$	saddle
$Q_{1} (u_{1}, v_{1})$	$k < \frac{h + 3 u_{1}^{2}}{2 u_{1}}$	unstable
	$k > \frac{h + 3 u_{1}^{2}}{2 u_{1}}$	stable
$Q_{2} (u_{2}, v_{2})$		saddle

Table 4. Properties of the fixed point

Q_{1} (u_{1}, v_{1})

of system (1.6).

Table 4. Properties of the fixed point

Q_{1} (u_{1}, v_{1})

of system (1.6).

Conditions			Eigenvalues	Properties
$Ω_{2}^{2} (1) - 4 Ω_{1} (1) < 0$	$Δ < - \frac{Ω_{2} (1)}{Ω_{1} (1)}$		$\| λ_{1} \| < 1,$ $\| λ_{2} \| < 1$	sink
	$Δ = - \frac{Ω_{2} (1)}{Ω_{1} (1)}$		$\| λ_{1} \| = 1$ or $\| λ_{2} \| = 1$	non-hyperbolic
	$Δ > - \frac{Ω_{2} (1)}{Ω_{1} (1)}$		$\| λ_{1} \| > 1,$ $\| λ_{2} \| > 1$	source
$Ω_{2}^{2} (1) - 4 Ω_{1} (1) = 0$	$Δ < - \frac{Ω_{2} (1)}{Ω_{1} (1)}$		$\| λ_{1} \| < 1,$ $\| λ_{2} \| < 1$	sink
	$Δ = - \frac{Ω_{2} (1)}{Ω_{1} (1)}$		$\| λ_{1} \| = 1$ or $\| λ_{2} \| = 1$	non-hyperbolic
	$Δ > - \frac{Ω_{2} (1)}{Ω_{1} (1)}$		$\| λ_{1} \| > 1,$ $\| λ_{2} \| > 1$	source
$Ω_{2}^{2} (1) - 4 Ω_{1} (1)$ >0	$Ω_{2} (1) > 2 \sqrt{Ω_{1} (1)}$		$\| λ_{1} \| > 1,$ $\| λ_{2} \| > 1$	source
	$Ω_{2} (1) < - 2 \sqrt{Ω_{1} (1)}$	$Δ < \frac{- Ω_{2} (1) - \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)}$	$\| λ_{1} \| < 1$ , $\| λ_{2} \| < 1$	sink
		$Δ = \frac{- Ω_{2} (1) - \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)}$	$\| λ_{1} \| = (\neq) 1$ , $\| λ_{2} \| \neq (=) 1$	non-hyperbolic
		$\frac{- Ω_{2} (1) - \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)} < Δ < \frac{- Ω_{2} (1) + \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)},$	$\| λ_{1} \| < (>) 1$ , $\| λ_{2} \| > (<) 1$	saddle
		$Δ = \frac{- Ω_{2} (1) + \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)}$	$\| λ_{1} \| \neq (=) 1$ , $\| λ_{2} \| = (\neq) 1$	non-hyperbolic
		$Δ > \frac{- Ω_{2} (1) - \sqrt{Ω_{2}^{2} (1) - 4 Ω_{1} (1)}}{Ω_{1} (1)}$	$\| λ_{1} \| > 1$ , $\| λ_{2} \| > 1$	source

Table 5. Properties of the fixed point

Q_{2} (u_{2}, v_{2})

of system (1.6).

Table 5. Properties of the fixed point

Q_{2} (u_{2}, v_{2})

of system (1.6).

Conditions	Eigenvalues	Properties
$Δ < \frac{- Ω_{2} (2) - \sqrt{Ω_{2}^{2} (2) - 4 Ω_{1} (2)}}{Ω_{1} (2)}$	$\| λ_{1} \| > 1$ , $\| λ_{2} \| < 1$	saddle
$Δ = \frac{- Ω_{2} (2) - \sqrt{Ω_{2}^{2} (2) - 4 Ω_{1} (2)}}{Ω_{1} (2)}$	$\| λ_{1} \| > 1$ , $\| λ_{2} \| = 1$	non-hyperbolic
$Δ > \frac{- Ω_{2} (2) - \sqrt{Ω_{2}^{2} (2) - 4 Ω_{1} (2)}}{Ω_{1} (2)}$	$\| λ_{1} \| > 1$ , $\| λ_{2} \| > 1$	source

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Wang, B.; Li, X. Modeling and Dynamical Analysis of a Fractional-Order Predator–Prey System with Anti-Predator Behavior and a Holling Type IV Functional Response. Fractal Fract. 2023, 7, 722. https://doi.org/10.3390/fractalfract7100722

AMA Style

Wang B, Li X. Modeling and Dynamical Analysis of a Fractional-Order Predator–Prey System with Anti-Predator Behavior and a Holling Type IV Functional Response. Fractal and Fractional. 2023; 7(10):722. https://doi.org/10.3390/fractalfract7100722

Chicago/Turabian Style

Wang, Baiming, and Xianyi Li. 2023. "Modeling and Dynamical Analysis of a Fractional-Order Predator–Prey System with Anti-Predator Behavior and a Holling Type IV Functional Response" Fractal and Fractional 7, no. 10: 722. https://doi.org/10.3390/fractalfract7100722

APA Style

Wang, B., & Li, X. (2023). Modeling and Dynamical Analysis of a Fractional-Order Predator–Prey System with Anti-Predator Behavior and a Holling Type IV Functional Response. Fractal and Fractional, 7(10), 722. https://doi.org/10.3390/fractalfract7100722

Article Menu

Modeling and Dynamical Analysis of a Fractional-Order Predator–Prey System with Anti-Predator Behavior and a Holling Type IV Functional Response

Abstract

1. Introduction

2. Preliminaries

3. Analysis of the Well-Posedness of System (1.3)

4. Local Stability of Systems (1.3) and (1.6)

4.1. Existence of an Equilibrium Point

4.2. Stability Analysis of Equilibrium Points of System (1.3)

4.2.1. The Stability of the Trivial Equilibrium Point $Q_{0} (0, 0)$

4.2.2. The Stability of Boundary Equilibrium Point $Q_{k} (k, 0)$

4.2.3. The Stability of Positive Equilibrium Points $Q_{i} (u_{i}, v_{i}) (i = 1, 2)$

4.3. Stability Analysis of the Equilibrium Points of System (1.6)

4.3.1. The Stability of Trivial Equilibrium Point $Q_{0} (0, 0)$

4.3.2. The Stability of Boundary Equilibrium Point $Q_{k} (k, 0)$

4.3.3. The Stability of Positive Equilibrium Points $Q_{i} (u_{i}, v_{i}) (i = 1, 2)$

5. Bifurcation Analysis

5.1. Bifurcation Analysis of the Positive Equilibrium Point $Q_{1} (u_{1}, v_{1})$ in System (1.3)

5.2. Bifurcation Analysis of the Positive Equilibrium Point $Q_{1} (u_{1}, v_{1})$ in System (1.6)

5.2.1. Neimark–Sacker Bifurcation at the Fixed Point $Q_{1} (u_{1}, v_{1})$

5.2.2. Period-Doubling Bifurcation at the Fixed Point $Q_{1} (u_{1}, v_{1})$

6. Numerical Simulation

7. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI

Article Menu

Modeling and Dynamical Analysis of a Fractional-Order Predator–Prey System with Anti-Predator Behavior and a Holling Type IV Functional Response

Abstract

1. Introduction

2. Preliminaries

3. Analysis of the Well-Posedness of System (1.3)

4. Local Stability of Systems (1.3) and (1.6)

4.1. Existence of an Equilibrium Point

4.2. Stability Analysis of Equilibrium Points of System (1.3)

4.2.1. The Stability of the Trivial Equilibrium Point Q 0 ( 0 , 0 )

4.2.2. The Stability of Boundary Equilibrium Point Q k ( k , 0 )

4.2.3. The Stability of Positive Equilibrium Points Q i ( u i , v i ) ( i = 1 , 2 )

4.3. Stability Analysis of the Equilibrium Points of System (1.6)

4.3.1. The Stability of Trivial Equilibrium Point Q 0 ( 0 , 0 )

4.3.2. The Stability of Boundary Equilibrium Point Q k ( k , 0 )

4.3.3. The Stability of Positive Equilibrium Points Q i ( u i , v i ) ( i = 1 , 2 )

5. Bifurcation Analysis

5.1. Bifurcation Analysis of the Positive Equilibrium Point Q 1 ( u 1 , v 1 ) in System (1.3)

5.2. Bifurcation Analysis of the Positive Equilibrium Point Q 1 ( u 1 , v 1 ) in System (1.6)

5.2.1. Neimark–Sacker Bifurcation at the Fixed Point Q 1 ( u 1 , v 1 )

5.2.2. Period-Doubling Bifurcation at the Fixed Point Q 1 ( u 1 , v 1 )

6. Numerical Simulation

7. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI

4.2.1. The Stability of the Trivial Equilibrium Point $Q_{0} (0, 0)$

4.2.2. The Stability of Boundary Equilibrium Point $Q_{k} (k, 0)$

4.2.3. The Stability of Positive Equilibrium Points $Q_{i} (u_{i}, v_{i}) (i = 1, 2)$

4.3.1. The Stability of Trivial Equilibrium Point $Q_{0} (0, 0)$

4.3.2. The Stability of Boundary Equilibrium Point $Q_{k} (k, 0)$

4.3.3. The Stability of Positive Equilibrium Points $Q_{i} (u_{i}, v_{i}) (i = 1, 2)$

5.1. Bifurcation Analysis of the Positive Equilibrium Point $Q_{1} (u_{1}, v_{1})$ in System (1.3)

5.2. Bifurcation Analysis of the Positive Equilibrium Point $Q_{1} (u_{1}, v_{1})$ in System (1.6)

5.2.1. Neimark–Sacker Bifurcation at the Fixed Point $Q_{1} (u_{1}, v_{1})$

5.2.2. Period-Doubling Bifurcation at the Fixed Point $Q_{1} (u_{1}, v_{1})$