Complex Dynamics and PID Control Strategies for a Fractional Three-Population Model

Abstract

In recent decades, there have been many studies on Hopf bifurcation and population stability with time delay. However, the stability and Hopf bifurcation of fractional-order population systems with time delay are lower. In this paper, we discuss the dynamic behavior of a fractional-order three-population model with pregnancy delay using Laplace transform of fractional differential equations, stability and bifurcation theory, and MATLAB software. The specific conditions of local asymptotic stability and Hopf bifurcation for fractional-order time-delay systems are determined. A fractional-order proportional–integral–derivative (PID) controller is applied to the three-population food chain system for the first time. The convergent speed and vibration amplitude of the system can be changed by PID control. For example, after fixing the values of the integral control gain $k_i$ and the differential control gain $k_d$, the amplitude of the system decreases and the convergence speed changes as the proportional control gain $k_p$ decreases. The effectiveness of the PID control strategy in complex ecosystem is proved. The numerical simulation results are in good agreement with the theoretical analysis. The research in this paper has potential application values concerning the management of complex population systems. The bifurcation theory of fractional-order time-delay systems is also enriched.

1. Introducing

The relationship between prey and predators is crucial in the ecosystem, playing a vital role in shaping population structure and ensuring the stability of the ecosystem. Over the recent few years, researchers have dedicated a considerable amount of effort to studying this area, resulting in valuable insights and a plethora of research findings [1,2]. These ecological models have given us a way to gain a deeper understanding of the intricacies of ecosystems and anticipate the results of ecological advancements. They have also been used as a theoretical foundation and reference for the future development of large-scale ecological models. In 1959, a model was proposed [1] that considers the functional response of predators and the growth period of prey. In 1991, Hastings et al. [2] expanded on the Hollin model and found the presence of chaotic behavior in the model of a three-population food chain. Recently, there has been a growing fascination with studying population ecosystems that involve fractional-order dynamics [3,4,5,6]. These studies focus on examining the behavior and stability of ecosystems with three species that interact with each other.

The origins of fractional calculus can be traced back to the early days of traditional calculus, pioneered by Newton and Leibniz. During its early development, fractional calculus faced challenges due to its limited theoretical basis and practical application background. This was not changed until 1983, when fractal theory was introduced in [7], and the application of fractional calculus has been widely developed. It is well known that integer-order calculus primarily focuses on the local characteristics of functions, while fractional-order calculus takes into account long memory and heredity, allowing for a better representation of the global characteristics of functions [8]. Fractional calculus has found extensive application in describing infectious disease models and ecological models [9,10,11]. Ji et al. [12] examined systems with two fractional-order species and presented parameters for a stable system. Das et al. [13] provided a systematic introduction to fractional differential equations and their applications. Tavazoei et al. [14] demonstrated that in systems with a fractional order, the final set of trajectories may not align with the system’s solution, which distinguishes it from integer-order systems. Tavazoei et al. [15] demonstrated that no periodic solutions exist in time-invariant fractional systems. They suggested a system that involves trajectories that are not periodic but eventually converge to periodic signals. Aguila-Camacho et al. [16] demonstrated stability for various fractional systems by utilizing a fractional extension of Lyapunov’s direct method.

Population modeling has a long history of utilizing delay differential equations (DDEs). Time delay is typically seen as the period between the response of a species or process and the interaction of other species or processes. This delay could be caused by physiological processes between species, the rate of transmission of environmental changes, or the rate of population migration, such as gestation, maturity, and feeding times [17,18,19,20,21,22]. Regarding mathematical modeling, delayed systems have been found to be more effective than non-delayed systems [23,24]. Time delays in the predator–prey model lead to a more intricate dynamic behavior. One of the factors to consider is gestation delay, which refers to the time it takes for biomass to be transferred from prey to a predator. As an illustration, Tripathi et al. [25] studied a system that incorporated discrete gestation delays. The study focused on Beddington–Deangelis-type functional response functions to investigate the level of mutual interference among density-dependent predators. Zheng et al. [26] considered a stage-structured internecine model with two delays, one representing the lag effect of negative feedback from prey species and the other representing the pregnancy delay of adult predator populations.

Wang et al. [27] introduced a novel concept for secure telecommunication synchronization leveraging neural networks, with a particular emphasis on memristor-based architectures. This selection stems from the direct correlation between fractional calculus and memristors, particularly in their shared “memory effect”. Notably, in the realm of neural networks, fractional-order memristors can emulate synaptic connections among neurons, thereby enabling the realization of more intricate neural network structures and functionalities. Ouannas et al. [28] employed both linear and non-linear controllers for synchronization purposes in secure telecommunication applications. The enhanced Robinovich chaotic system they utilized is of the fractional type. The numerical findings underscore the efficacy of their proposed approach in safeguarding information confidentiality. In our future work, we will try to use the concept of fractional calculus to make fractional-order memristors used to construct chaotic circuits and cryptographic systems.

Many control methods developed for integer derivative systems cannot be applied to fractional derivative systems due to the complexity of the characteristic equations involving exponential fractional pseudo-polynomial functions. Chen et al. [29] introduced a group of gene regulatory networks with fractional-order dynamics, incorporating time delays and parameter uncertainties. Through the implementation of an adaptive controller, they successfully achieved synchronization in a finite amount of time. Kekha Javan et al. [30] investigated the polymorphic synchronization of chaotic systems with non-identical, unknown, and time-varying delays in the presence of external perturbations and parameter uncertainties. Two Chen hyperchaotic systems were simulated as slave systems and one Chua hyperchaotic system as master systems. The results show that its control effect is better than that of the most advanced technology. A method of fractional-order proportional derivative (PD) control has been thoroughly developed and studied [31,32]. PD control has a tendency to create additional equilibrium points [33]. To address this issue, integer-order proportional–integral–derivative (PID) control was introduced in order to provide a more effective solution [34].

Lu et al. [35] introduced the fractional-order PD control strategy for a delayed fractional-order prey–predator system, and the fractional-order PD controller has been applied to realize Hopf bifurcation control for the fractional-order systems with time delay. The fractional-order PID controller can effectively replace the fractional-order PD controllers by adjusting their control parameters. In this paper, a method using a fractional-order PID controller to solve the optimal control problem of Hopf bifurcation in a three-species food chain model with the fractional-order Caputo derivative is introduced. The fractional-order PID controller is rarely applied to three-population system in previous studies. We use the fractional-order Caputo derivative to describe the dynamic characteristics of the system and introduce the fractional-order PID controller to stabilize the system. This method is significantly useful for solving the Hopf bifurcation problem in food chain models. The work in [35] is enriched, and some interesting conclusions are drawn.

The contributions of this paper are as follows:

(1) In the population system, there are few studies on fractional order. In this paper, a fractional-order three-population model with pregnancy delay is explored, and the stability and Hopf bifurcation conditions of the system are given with the delay as the bifurcation parameter.

(2) In contrast to previous work, we have studied all possible cases of the local stability of the delay-free model in detail.

(3) We demonstrate the feasibility of using PID controllers in ecological systems. The experimental results show that selecting appropriate control parameters in the PID controller can change the system’s convergence rate and vibration amplitude.

In our paper, we present Section 2, where we introduce a three-population ecosystem with Caputo fractional-order sense that involves time delay. Section 3 presents the local stability analysis of the system and provides the conditions for local stability. In Section 4, it is shown that the occurrence of Hopf bifurcation is dependent on the delay parameter exceeding a critical value. Section 5 introduces the PID controller to the system, analyzing the eigenvalues of the controlled system to obtain the dynamic stability criterion and Hopf bifurcation condition. Finally, Section 6 includes numerical simulations that confirm the theoretical analysis.

2. Introduction and Model Description

2.1. Caputo Fractional-Order Derivative

In mathematical research, there are different definitions of fractional calculus, but three common forms exist: Fractional derivatives such as Riemann–Liouville, Grunwald–Letnikov, and Caputo are all of the same size. The Laplace transform of the Caputo derivative is solely determined by the initial value of the integer derivative, with the calculation of the initial value of the Riemann–Liouville derivative posing challenges. It is straightforward to determine the initial value of the integer derivative [8]. Therefore, the Caputo derivative is extensively utilized in practice, and the system in this paper employs the Caputo fractional derivative for discussion.

Definition 1.

The definition of the Caputo fractional derivative of a continuous function $y\left(t\right)$ with respect to ε is defined as follows:

(1)$\begin{array}{cc}\hfill {}_a{}^{C}D_t^{\epsilon }f\left(t\right)& =I^{n-\epsilon }\left[{\displaystyle \frac{{\mathrm{d}}^{\mathrm{n}}}{\mathrm{d}{t}^n}}f\left(t\right)\right]={\displaystyle \frac{1}{\Gamma (n-\epsilon )}}{\int }_a^t{f}^n\left(\tau \right){(t-\tau )}^{n-\epsilon -1}\mathrm{d}\tau ,\hfill \end{array}$

where $n-1<\epsilon \le n$, $n\in {\mathbb{Z}}^{+}$, Γ is the $Gamma$ function in the form $\Gamma \left(s\right)={\int }_0^{\infty }{t}^{s-1}{e}^{-t}\mathrm{d}s$.

Definition 2.

The definition of the Riemann–Liouville derivative is

(2) $\begin{array}{cc}\hfill {}_0{}^{RL}D_t^{\epsilon }f\left(t\right)& ={\displaystyle \frac{{\mathrm{d}}^{\mathrm{k}}}{\mathrm{d}{t}^k}}\left({}_{0}J_t^{\epsilon -k}\right)={\displaystyle \frac{1}{\Gamma (k-\epsilon )}}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{k}}}{\mathrm{d}{t}^k}}\left[{\int }_0^{t}f\left(\tau \right){(t-\tau )}^{n-\epsilon -1}\mathrm{d}\tau \right],\hfill \end{array}$

where $n-1<\epsilon \le n$, $n\in {\mathbb{Z}}^{+}$, $t>0$.

Lemma 1.

The Laplace transform of the Riemann–Liouville derivative is represented by

(3) $L\left[{}_0{}^{RL}D_t^{\epsilon }f\left(t\right)\right]=s^{\epsilon F\left(s\right)}-\sum _0^{k-1}{\left[{}_0{}^{RL}D_t^{\epsilon -j-1}{f}^{\left(j\right)}\left(t\right)\right]}_{t=0}{s}^j.$

From (1) and (2), we can derive the following:

(4) $\begin{array}{cc}\hfill {}_0{}^{C}D_x^{\epsilon }{[}_0{J}_x^{-\epsilon }f\left(x\right)]& ={}_{0}J_x^{\epsilon -k}\left\{{\displaystyle \frac{{\mathrm{d}}^{\mathrm{k}}}{\mathrm{d}{t}^k}}{[}_0{J}_x^{-\epsilon }f\left(x\right)]\right\}\hfill \\ & ={}_{0}J_x^{\epsilon -k}\left[{}_0{}^{RL}D_t^{k-\epsilon }f\left(x\right)\right]\hfill \\ & =f\left(x\right)-{\left[{}_{0}J_x^{k-\epsilon -1}f\left(x\right)\right]}_{x=0}{\displaystyle \frac{x^{k-\epsilon -1}}{\Gamma (k-\epsilon )}}\hfill \\ & =f\left(x\right).\hfill \end{array}$

Theorem 1

([14]). Consider the following fractional-order system:

(5)${}_a{}^{C}D_t^{\epsilon }y=f\left(y\right),$

where $0<\epsilon \le 1$, $t>0$ and $y\in {\mathcal{R}}^n$. The equilibrium points of System (5) are locally asymptotically stable if all eigenvalues λ of the Jacobian matrix $J\left(y\right)={\displaystyle \frac{\partial f}{\partial y}}$ satisfy the condition $|arg\left(\lambda \right)|>{\displaystyle \frac{\epsilon \pi }{2}}$. Figure 1 shows the stability bounds $|arg\left(\lambda \right)|={\displaystyle \frac{\epsilon \pi }{2}}$.

Meanwhile, the Caputo fractional derivative exhibits the following properties:

Property 1.

Let $y\left(t\right)\in {\mathcal{R}}^n$ be a continuous and derivable function; then, for any time instant $t\ge a$, we have

(6) ${\displaystyle \frac{1}{2}}{}_a{}^{C}D_t^{\epsilon }{y}^2\left(t\right)\le y\left(t\right){}_a{}^{C}D_t^{\epsilon }y\left(t\right).$

Property 2.

When $x\left(t\right)\in C[0,T]$, $\gamma ,\epsilon \in {\mathcal{R}}^{+}$ and $\gamma +\epsilon \le 1$, for some $T>0$, then

(7) $\begin{array}{c}\hfill D^{\gamma }{D}^{\epsilon }x\left(t\right)=D^{\epsilon }{D}^{\gamma }x\left(t\right)=D^{\gamma +\epsilon }x\left(t\right),\gamma +\epsilon <1,t\in [0,T].\end{array}$

Property 3.

For $\forall \epsilon >0$ and for $A_1,A_2\in \mathcal{R}$, we have

(8) ${}_a{}^{C}D_t^{\epsilon }(A_{1}x\left(t\right)+A_{2}y\left(t\right))=A_1{}_a{}^{C}D_t^{\epsilon }x\left(t\right)+A_2{}_a{}^{C}D_t^{\epsilon }y\left(t\right).$

Property 4.

For any constant q, we have ${}_a{}^{C}D_t^{\epsilon }q=0$.

2.2. Model Descriptions

In [36], the following model is proposed:

(9)$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}U}{\mathrm{d}T}}=& R_{0}U{\displaystyle \frac{1+E_{1}U}{1+E_{1}U+F_{1}V}}(1-{\displaystyle \frac{U}{K}})-C_1{\displaystyle \frac{A_{1}UV}{B_1+U}}\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}V}{\mathrm{d}T}}=& {\displaystyle \frac{A_{1}UV}{B_1+U}}{\displaystyle \frac{1+E_{2}V}{1+E_{2}V+F_{2}W}}-{\displaystyle \frac{A_{2}VW}{B_2+V}}-D_{1}V\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}W}{\mathrm{d}T}}=& {\displaystyle \frac{C_2{A}_{2}VW}{B_2+V}}-D_{2}W.\hfill \end{array}$

We reduce the parameters of the aforementioned system, System (9), using the method described in [36]. System (9) is then transformed by $U=Kx$, $V={\displaystyle \frac{K}{C_1}}y$, $W={\displaystyle \frac{C_{2}K}{C_1}}z$, $T={\displaystyle \frac{t}{R_0}}$, ${\alpha }_1={\displaystyle \frac{A_{1}K}{R_0{B}_1}}$, ${\alpha }_2={\displaystyle \frac{A_2{C}_{2}K}{R_0{B}_2{C}_1}}$, ${\beta }_1={\displaystyle \frac{K}{B_1}}$, ${\beta }_2={\displaystyle \frac{K}{C_1{B}_2}}$, ${\delta }_1={\displaystyle \frac{D_1}{R_0}}$, ${\delta }_2={\displaystyle \frac{D_3}{R_0}}$, $f_1={\displaystyle \frac{F_{1}K}{C_1}}$, $f_2={\displaystyle \frac{F_2{C}_{2}K}{C_1}}$, $e_1=E_{1}K$ and $e_2={\displaystyle \frac{E_{2}K}{C_1}}$. As a result, the following system is obtained:

(10)$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}}& =x\left((1-x)·{\displaystyle \frac{1+e_{1}x}{1+e_{1}x+f_{1}y}}-{\displaystyle \frac{{\alpha }_{1}y}{1+{\beta }_{1}x}}\right)\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}}& =y\left({\displaystyle \frac{{\alpha }_{1}x}{1+{\beta }_{1}x}}·{\displaystyle \frac{1+e_{2}y}{1+e_{2}y+f_{2}z}}-{\delta }_1-{\displaystyle \frac{{\alpha }_{2}z}{1+{\beta }_{2}y}}\right)\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}z\left(t\right)}{\mathrm{d}t}}& =z\left({\displaystyle \frac{{\alpha }_{2}y}{1+{\beta }_{2}y}}-{\delta }_2\right).\hfill \end{array}$

Time delays are included in the model to create a more accurate representation, as shown in [19,20,21,22,23,24]. We will discover that time delays have a significant impact on the dynamics of ecological models, often resulting in oscillations, bifurcations, and even chaos. However, fractional derivatives possess memory and heritability, making them extensively utilized in the examination of population models [4,5,6]. This paper explores the impact of time delay on the system and incorporates the Caputo fractional-order function into the ordinary differential equation of System (10). The resulting fractional differential equation is described as follows:

(11)$\begin{array}{cc}\hfill {}_0{}^{C}D_t^{\epsilon }x\left(t\right)& =x\left(t\right)\left((1-x\left(t\right))·{\displaystyle \frac{1+e_{1}x\left(t\right)}{1+e_{1}x\left(t\right)+f_{1}y\left(t\right)}}-{\displaystyle \frac{{\alpha }_{1}y\left(t\right)}{1+{\beta }_{1}x\left(t\right)}}\right)\hfill \\ \hfill {}_0{}^{C}D_t^{\epsilon }y\left(t\right)& ={\displaystyle \frac{{\alpha }_{1}x\left(t\right)y\left(t\right)}{1+{\beta }_{1}x\left(t\right)}}·{\displaystyle \frac{1+e_{2}y\left(t\right)}{1+e_{2}y\left(t\right)+f_{2}z\left(t\right)}}-{\delta }_{1}y\left(t\right)-{\displaystyle \frac{{\alpha }_{2}y\left(t\right)z\left(t\right)}{1+{\beta }_{2}y\left(t\right)}}\hfill \\ \hfill {}_0{}^{C}D_t^{\epsilon }z\left(t\right)& ={\displaystyle \frac{{\alpha }_{2}y(t-\tau )z(t-\tau )}{1+{\beta }_{2}y(t-\tau )}}-{\delta }_{2}z\left(t\right),\hfill \end{array}$

where $\epsilon \in (0,1]$. Let $C([-\tau ,0],{\mathcal{R}}_{+0}^3)$, which belongs to the Banach space of smooth functions mapping interval $[-\tau ,0]$ into ${\mathcal{R}}_{+0}^3$. For any $\vartheta \in C\left([-\tau ,0]\right),{\mathcal{R}}_{+0}^3$, we define the norm $\Vert \vartheta \Vert =\underset{-\tau \le t\le 0}{sup}\mid \left(\theta \right)\mid$ with the initial condition of System (11) being

$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill {\vartheta }_x\left(\theta \right)& \ge 0,{\vartheta }_y\left(\theta \right)\ge 0,{\vartheta }_z\left(\theta \right)\ge 0,\vartheta \in [-\tau ,0],\hfill \\ \hfill x\left(0\right)& \ge 0,y\left(0\right)\ge 0,z\left(0\right)\ge 0.\hfill \end{array}\right.\end{array}$

$\tau$ represents the gestation delay of the senior predator. In the prey, primary predator, and senior predator models, it is assumed that prey and primary predators have shorter reproductive cycles compared to senior predators. This implies that the population number does not immediately reflect the number of predators, even after the birth of a predator. Therefore, this delay effect can be simulated by incorporating a time delay in the predator population growth term.

In [37], Cui et al. studied the dynamical behavior of System (11) at $\tau =0$.

(12)$\begin{array}{cc}\hfill {}_0{}^{C}D_t^{\epsilon }x\left(t\right)& =x\left(t\right)\left((1-x\left(t\right))·{\displaystyle \frac{1+e_{1}x\left(t\right)}{1+e_{1}x\left(t\right)+f_{1}y\left(t\right)}}-{\displaystyle \frac{{\alpha }_{1}y\left(t\right)}{1+{\beta }_{1}x\left(t\right)}}\right)\hfill \\ \hfill {}_0{}^{C}D_t^{\epsilon }y\left(t\right)& =y\left(t\right)\left({\displaystyle \frac{{\alpha }_{1}x\left(t\right)}{1+{\beta }_{1}x\left(t\right)}}·{\displaystyle \frac{1+e_{2}y\left(t\right)}{1+e_{2}y\left(t\right)+f_{2}z\left(t\right)}}-{\delta }_1-{\displaystyle \frac{{\alpha }_{2}z\left(t\right)}{1+{\beta }_{2}y\left(t\right)}}\right)\hfill \\ \hfill {}_0{}^{C}D_t^{\epsilon }z\left(t\right)& =z\left(t\right)\left({\displaystyle \frac{{\alpha }_{2}y\left(t\right)}{1+{\beta }_{2}y\left(t\right)}}-{\delta }_2\right).\hfill \end{array}$

3. Analysis of Non-Delayed System

3.1. Equilibrium Point

By applying an identical approach used for determining the equilibrium points of the integer-order system, we can derive the four equilibrium points of System (12) as follows: $E_0(0,0,0)$, $E_1(1,0,0)$, $E_2(x_1^{*},y_1^{*},0)$, and $E_3(x_2^{*},y_2^{*},z_2^{*})$.

Remark 1.

The new variables are explained in Appendix A.

Remark 2.

We find the possible equilibrium points of the system, including the no-higher predator equilibrium, the predator-free equilibrium, and the population-free equilibrium. The equilibrium point of the system is determined by a combination of factors such as the interaction of biological populations, environmental capacity, growth and mortality rates, the characteristics of fractional derivatives, external disturbances, and initial conditions. These factors reflect the interactions between biological populations in ecosystems and their relationship with the environment. It is not difficult to find that any extreme conditions can easily lead to the extinction of the population. Therefore, in the process of protecting biodiversity and maintaining ecological balance, it is necessary to comprehensively consider the impact of these factors and take corresponding conservation measures to reduce the risk of population extinction.

3.2. Local Stability Analysis of Non-Time-Delay System

We analyze the local stability of System (12) by evaluating the Jacobian matrix at the equilibrium point. This section discusses the local stability and the parameters that affect stability.

(13)$J\left(E(x,y,z)\right)=\left[\begin{array}{ccc}{A}_{11}& A_{12}& A_{13}\\ A_{21}& A_{22}& A_{23}\\ A_{31}& A_{32}& A_{33}\end{array}\right].$

At $E_0$, the Jacobian matrix is given by

(14)$J\left(E_0(0,0,0)\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& -{\delta }_1& 0\\ 0& 0& -{\delta }_2\end{array}\right].$

The eigenvalues of the system at the equilibrium point $E_0$ are ${\lambda }_1=1$, ${\lambda }_2=-{\delta }_1$ and ${\lambda }_3=-{\delta }_2$. When $0<\epsilon \le 1$, $arg\left({\lambda }_1\right)=0<{\displaystyle \frac{\epsilon \pi }{2}}$ and $arg\left({\lambda }_{2,3}\right)=\pi >{\displaystyle \frac{\epsilon \pi }{2}}$. According to Theorem 1, $E_0$ is unstable.

The Jacobian matrix at $E_1$ is as follows:

(15)$J\left(E_1(1,0,0)\right)=\left[\begin{array}{ccc}-1& 0& 0\\ 0& -{\delta }_2& 0\\ 0& 0& -{\displaystyle \frac{{\delta }_1-{\alpha }_1+{\beta }_1{\delta }_1}{{\beta }_1+1}}\end{array}\right].$

The eigenvalues of the matrix are ${\lambda }_1=-1$, ${\lambda }_2=-{\delta }_2$ and ${\lambda }_3=-{\displaystyle \frac{{\delta }_1-{\alpha }_1+{\beta }_1{\delta }_1}{{\beta }_1+1}}$. If ${\delta }_1>{\displaystyle \frac{{\alpha }_1}{{\beta }_1+1}}$, $arg\left({\lambda }_{1,2}\right)=\pi >{\displaystyle \frac{\epsilon \pi }{2}}$ and $arg\left({\lambda }_3\right)=\pi >{\displaystyle \frac{\epsilon \pi }{2}}$. According to Theorem 1, $E_1$ is locally asymptotically stable.

The Jacobian matrix at $E_2(x_1^{*},y_1^{*},0)$ is given by

(16)$J\left(E_2(x_1^{*},y_1^{*},0)\right)=\left[\begin{array}{ccc}{B}_1& -B_2& 0\\ B_3& 0& -B_4\\ 0& 0& B_5\end{array}\right].$

The characteristic equation is $({\lambda }^2-B_1\lambda +B_2{B}_3)(B_5-\lambda )=0$. Then, we ordered that $b_1={\displaystyle \frac{B_1}{2}}$ and $b_2={\displaystyle \frac{\sqrt{B_1^2-4B_2{B}_3}}{2}}$. Now, the roots of the characteristic equation are ${\lambda }_1=B_5$ and ${\lambda }_{2,3}=b_1\pm b_2$. Therefore, $arg\left({\lambda }_1\right)=\pi >{\displaystyle \frac{\epsilon \pi }{2}}$ and $arg\left({\lambda }_{2,3}\right)=\pi >{\displaystyle \frac{\epsilon \pi }{2}}$, if $B_1<0$, $B_2>0$ and $B_5<0$. According to Theorem 1, the senior predator-free equilibrium $E_2$ is asymptotically stable.

The Jacobian matrix at $E_3$ is as follows:

(17)$J\left(E_3(x_2^{*},y_2^{*},z_2^{*})\right)=\left[\begin{array}{ccc}{C}_1& -C_2& 0\\ C_3& C_4& -C_5\\ 0& C_6& 0\end{array}\right].$

The characteristic equation of Matrix (17) is

(18)$f\left(\lambda \right)={\lambda }^3+R_1{\lambda }^2+R_2\lambda +R_3=0,$

where $R_1=-(C_1+C_4)$, $R_2=C_1{C}_4+C_2{C}_3+C_5{C}_6$, and $R_3=-C_1{C}_5{C}_6$.

Using the method in [38], we introduce the following discriminant:

(19)$\Delta =18R_1{R}_2{R}_3-4R_1^3{R}_3+R_1^2{R}_2^2-4R_2^3-27R_3^2.$

According to the content of [38], we can obtain the following:

(1). If $\Delta >0$, $R_1>0$, $R_2>0$ and $R_3>0$, then $f\left(\lambda \right)$ has three distinct negative real roots.

(2). If $\Delta >0$, $R_1>0$ and $R_3<0$, then $f\left(\lambda \right)$ has one positive real root and two negative real roots.

(3). If $\Delta >0$, $R_1>0$, $R_2<0$ and $R_3>0$, then $f\left(\lambda \right)$ has one negative real root and two positive real roots.

(4). If $\Delta >0$, $R_1<0$, $R_2>0$ and $R_3<0$, then $f\left(\lambda \right)$ has three positive real roots.

(5). If $\Delta <0$, $R_1>0$, $R_2>0$ and $R_3>0$, then $f\left(\lambda \right)$ has one negative real root and a pair of complex conjugate roots with negative real parts.

(6). If $\Delta <0$, $R_1>0$ and $R_2>0$, then $f\left(\lambda \right)$ has one negative real root and a pair of complex conjugate roots with positive real parts.

(7). If $\Delta <0$, $R_1>0$ and $R_3<0$, then $f\left(\lambda \right)$ has one positive real root and a pair of complex conjugate roots with negative real parts.

Based on the above analysis, we obtain the following seven regions:

${\Omega }_1$ = {$(R_1,R_2,R_3)$, $\Delta >0$, $R_1>0$, $R_2>0$, $R_3>0$}

${\Omega }_2$ = {$(R_1,R_2,R_3)$, $\Delta >0$, $R_1>0$, $R_3<0$}

${\Omega }_3$ = {$(R_1,R_2,R_3)$, $\Delta >0$, $R_1>0$, $R_2<0$, $R_3>0$}

${\Omega }_4$ = {$(R_1,R_2,R_3)$, $\Delta >0$, $R_1<0$, $R_2>0$, $R_3<0$}

${\Omega }_5$ = {$(R_1,R_2,R_3)$, $\Delta <0$, $R_1>0$, $R_2>0$, $R_3>0$}

${\Omega }_6$ = {$(R_1,R_2,R_3)$, $\Delta <0$, $R_1>0$, $R_2>0$}

${\Omega }_7$ = {$(R_1,R_2,R_3)$, $\Delta <0$, $R_1>0$, $R_3<0$}.

Table 1 presents the correlation between the roots of the characteristic equation and the coefficients.

Theorem 2.

Set up the following parameters:

$\begin{array}{cc}\hfill p& =-{\displaystyle \frac{R_1-3R_2}{3}}, q={\displaystyle \frac{2R_1-9R_1{R}_2+27R_3}{27}},\hfill \\ \hfill {\Delta }_1& =-{\displaystyle \frac{q}{2}}+\sqrt{{\left({\displaystyle \frac{q}{2}}\right)}^2+{\left({\displaystyle \frac{p}{3}}\right)}^3}, {\Delta }_2=-{\displaystyle \frac{q}{2}}-\sqrt{{\left({\displaystyle \frac{q}{2}}\right)}^2+{\left({\displaystyle \frac{p}{3}}\right)}^3}.\hfill \end{array}$

For System (12), the following results hold:

(1) If $(R_1,R_2,R_3)\in {\Omega }_1\bigcup {\Omega }_5$, then the equilibrium point $E_3$ of System (12) is locally asymptotically stable for $\epsilon \in (0,1]$.

(2) If $(R_1,R_2,R_3)\in {\Omega }_2\bigcup {\Omega }_3\bigcup {\Omega }_4\bigcup {\Omega }_7$, then the equilibrium point $E_3$ of System (12) is unstable for $\epsilon \in (0,1]$.

(3) If $(R_1,R_2,R_3)\in {\Omega }_6$ and

(20) $arctan{\displaystyle \frac{\sqrt{3}(\sqrt[3]{{\Delta }_2}-\sqrt[3]{{\Delta }_1})}{\sqrt[3]{{\Delta }_1}+\sqrt[3]{{\Delta }_2}+\frac{2}{3}{R}_1}}>{\displaystyle \frac{\epsilon \pi }{2}},$

then the equilibrium point $E_3$ of System (12) is locally asymptotically stable.

(4) If $(R_1,R_2,R_3)\in {\Omega }_6$ and

(21) $arctan{\displaystyle \frac{\sqrt{3}(\sqrt[3]{{\Delta }_2}-\sqrt[3]{{\Delta }_1})}{\sqrt[3]{{\Delta }_1}+\sqrt[3]{{\Delta }_2}+\frac{2}{3}{R}_1}}<{\displaystyle \frac{\epsilon \pi }{2}},$

then the equilibrium point $E_3$ of System (12) is unstable.

Proof. 

The condition of characteristic equation’s roots is evident according to Table 1.

(1) If $(R_1,R_2,R_3)\in {\Omega }_1\bigcup {\Omega }_5$, then all eigenvalues have negative real parts. This signifies that all eigenvalues of (18) are distributed in Region 1. According to Theorem 1, the equilibrium point is locally asymptotically stable at $E_3$.

(2) If $(R_1,R_2,R_3)\in {\Omega }_2\bigcup {\Omega }_3\bigcup {\Omega }_4\bigcup {\Omega }_7$, then at least one eigenvalue of Equation (18) has a positive real part. Consequently, the equilibrium point $E_3$ is unstable.

(3) If $(R_1,R_2,R_3)\in {\Omega }_6$, (18) has a negative real root ${\lambda }_1$ and a pair of conjugate complex roots ${\lambda }_{2,3}$ with positive real parts, as follows:

(22)$\begin{array}{ccc}\hfill {\lambda }_1& =& \sqrt[3]{{\Delta }_1}+\sqrt[3]{{\Delta }_2}-{\displaystyle \frac{R_1}{3}},\hfill \end{array}$

(23)$\begin{array}{ccc}\hfill {\lambda }_{2,3}& =& -{\displaystyle \frac{1}{2}}(\sqrt[3]{{\Delta }_1}+\sqrt[3]{{\Delta }_2})-{\displaystyle \frac{R_1}{3}}\pm {\displaystyle \frac{\sqrt{3}}{2}}(\sqrt[3]{{\Delta }_1}-\sqrt[3]{{\Delta }_2})\mathrm{i}.\hfill \end{array}$

Based on fractional stability theory, the following conclusions can be drawn:

(24)$\begin{array}{ccc}\hfill |arg{\lambda }_1|& =& \pi >{\displaystyle \frac{\pi \epsilon }{2}},\hfill \end{array}$

(25)$\begin{array}{ccc}\hfill |arg{\lambda }_{2,3}|& =& arctan{\displaystyle \frac{\sqrt{3}(\sqrt[3]{{\Delta }_2}-\sqrt[3]{{\Delta }_1})}{\sqrt[3]{{\Delta }_1}+\sqrt[3]{{\Delta }_2}+\frac{2}{3}{R}_1}},\hfill \end{array}$

and if Condition (20) is satisfied, according to Theorem 1, the equilibrium point $E_3$ of System (12) is locally asymptotically stable.

(4) This conclusion follows directly from (21). □

4. Analysis of Time-Delayed System ($\mathbf{\tau }>\mathbf{0}$)

System (11) and System (12) have an identical equilibrium point. The existence, uniqueness, boundedness, and non-negativity conditions are the same for both systems [39]. However, the stability criterion of the delay system is altered because of the presence of the delay parameter $\tau$. In the upcoming study, we will discover that Hopf bifurcation occurs within a specific range of the delay $\tau$, while the solution remains asymptotically stable outside this interval.

To analyze the stability of System (11), the Jacobian matrix of System (11) at $E^{*}(x^{*},y^{*},z^{*})$ is given by

(26)$J{\mid }_{E=E^{*}}=\left[\begin{array}{ccc}{s}^{\epsilon }-a_{11}& -a_{12}& 0\\ -a_{21}& s^{\epsilon }-a_{22}& -a_{23}\\ 0& -b_{32}{e}^{-s\tau }& s^{\epsilon }-a_{33}-b_{33}{e}^{-s\tau }\end{array}\right].$

The corresponding characteristic equation is

(27)${\varphi }_1\left(s\right)+{\varphi }_2\left(s\right)e^{-s\tau }=0,$

where

$\begin{array}{cc}\hfill {\varphi }_1\left(s\right)& =s^{3\epsilon }+G_1{s}^{2\epsilon }+G_2{s}^{\epsilon }+G_3, {\varphi }_2\left(s\right)=H_1{s}^{2\epsilon }+H_2{s}^{\epsilon }+H_3,\hfill \\ \hfill G_1& =-(a_{11}+a_{22}+a_{33}), G_2=a_{11}{a}_{22}+a_{11}{a}_{33}+a_{22}{a}_{33}-a_{12}{a}_{21},\hfill \\ \hfill G_3& =a_{12}{a}_{21}{a}_{33}-a_{11}{a}_{22}{a}_{33}, H_1=-b_{33}, H_2= a_{11}{b}_{33}+b_{33}{a}_{22}-a_{23}{b}_{32},\hfill \\ \hfill H_3& =a_{12}{b}_{32}{a}_{23}-a_{11}{a}_{22}{b}_{33}+a_{12}{a}_{21}{b}_{33}.\hfill \end{array}$

Theorem 3.

When $\tau >0$, if $E_{23}^2<F_{23}^2$ and $O_1{O}_3+O_2{O}_4>0$ hold, System (11) is locally asymptotically stable at the equilibrium point $E(x^{*},y^{*},z^{*})$ for all $\tau \in (0,{\tau }_0)$. When $\tau >{\tau }_0$, System (11) becomes unstable. Furthermore, Hopf bifurcation takes place at the equilibrium point $E(x^{*},y^{*},z^{*})$, when $\tau ={\tau }_0$, where

(28) ${\tau }_{1k}^{\left(j\right)}={\displaystyle \frac{1}{{\omega }_k}}\left(arccos{\displaystyle \frac{{\gamma }_{11}{\eta }_{11}+{\gamma }_{12}{\eta }_{12}}{{{\eta }_{11}}^2+{{\eta }_{12}}^2}}{|}_{\omega ={\omega }_k}+2j\pi \right),$

$k=1,2,3,4,5,6$, $j=0,1,2,\dots$, and ${\omega }_k$ are defined in the proof.

Proof. 

When $\tau >0$, the characteristic Equation (27) of System (11) becomes

(29)$s^{3\epsilon }+G_1{s}^{2\epsilon }+G_2{s}^{\epsilon }+G_3+(H_1{s}^{2\epsilon }+H_2{s}^{\epsilon }+H_3)e^{-s\tau }=0.$

We substitute $s=i\omega =\omega (cos{\displaystyle \frac{\pi }{2}}+isin{\displaystyle \frac{\pi }{2}})$ into (29); then, we have $s^{3\epsilon }={\omega }^{3\epsilon }(cos{\displaystyle \frac{3\epsilon \pi }{2}}+isin{\displaystyle \frac{3\epsilon \pi }{2}})$, $s^{2\epsilon }={\omega }^{2\epsilon }(cos\left(\epsilon \pi \right)+isin(\epsilon \pi ))$, $s^{\epsilon }={\omega }^{\epsilon }(cos{\displaystyle \frac{\epsilon \pi }{2}}+isin{\displaystyle \frac{\epsilon \pi }{2}})$ and $e^{-s\tau }=cos\left(\omega \tau \right)-isin\left(\omega \tau \right)$. Therefore, (29) is transformed into the following:

(30)$\begin{array}{cc}& {\omega }^{3\epsilon }(cos{\displaystyle \frac{3\epsilon \pi }{2}}+isin{\displaystyle \frac{3\epsilon \pi }{2}})+G_1{\omega }^{2\epsilon }(cos\left(\epsilon \pi \right)+isin(\epsilon \pi ))+G_2{\omega }^{\epsilon }(cos{\displaystyle \frac{\epsilon \pi }{2}}+isin{\displaystyle \frac{\epsilon \pi }{2}})\hfill \\ & +\left[H_1{\omega }^{2\epsilon }(cos\left(\epsilon \pi \right)+isin\left(\epsilon \pi \right))+H_2{\omega }^{\epsilon }(cos{\displaystyle \frac{\epsilon \pi }{2}}+isin{\displaystyle \frac{\epsilon \pi }{2}})+H_3\right][cos\left(\omega \tau \right)-isin(\omega \tau )]\hfill \\ & +G_3=0.\hfill \end{array}$

Detaching the imaginary and real components of the aforementioned equations yields the following:

(31)$\begin{array}{cc}\hfill {\eta }_{11}cos\left(\omega \tau \right)+{\eta }_{12}sin\left(\omega \tau \right)& ={\gamma }_{11},\hfill \end{array}$

(32)$\begin{array}{cc}\hfill {\eta }_{12}cos\left(\omega \tau \right)-{\eta }_{11}sin\left(\omega \tau \right)& ={\gamma }_{12},\hfill \end{array}$

and

$\begin{array}{cc}\hfill {\eta }_{11}& =H_1{\omega }^{2\epsilon }cos\left(\epsilon \pi \right)+H_2{\omega }^{\epsilon }cos{\displaystyle \frac{\epsilon \pi }{2}}+H_3, {\eta }_{12}=H_1{\omega }^{2\epsilon }sin\left(\epsilon \pi \right)+H_2{\omega }^{\epsilon }sin{\displaystyle \frac{\epsilon \pi }{2}},\hfill \\ \hfill {\gamma }_{11}& =-{\omega }^{3\epsilon }cos{\displaystyle \frac{3\epsilon \pi }{2}}-G_1{\omega }^{2\epsilon }cos\left(\epsilon \pi \right)-G_2{\omega }^{\epsilon }cos{\displaystyle \frac{\epsilon \pi }{2}}-G_3,\hfill \\ \hfill {\gamma }_{12}& =-{\omega }^{3\epsilon }sin{\displaystyle \frac{3\epsilon \pi }{2}}-G_1{\omega }^{2\epsilon }sin\left(\epsilon \pi \right)-G_2{\omega }^{\epsilon }sin{\displaystyle \frac{\epsilon \pi }{2}},\hfill \end{array}$

and by rearranging (31) and (32), we obtain

(33)$\begin{array}{cc}\hfill cos\left(\omega \tau \right)=& {\displaystyle \frac{{\gamma }_{11}{\eta }_{11}+{\gamma }_{12}{\eta }_{12}}{{{\eta }_{11}}^2+{{\eta }_{12}}^2}},\hfill \end{array}$

(34)$\begin{array}{cc}\hfill sin\left(\omega \tau \right)=& {\displaystyle \frac{{\gamma }_{11}{\eta }_{12}-{\gamma }_{12}{\eta }_{11}}{{{\eta }_{11}}^2+{{\eta }_{12}}^2}}.\hfill \end{array}$

As a consequence of ${sin}^2\left(\omega \tau \right)+{cos}^2\left(\omega \tau \right)=1$, we have

(35)${\omega }^{6\epsilon }+V_{11}{\omega }^{5\epsilon }+V_{12}{\omega }^{4\epsilon }+V_{13}{\omega }^{3\epsilon }+V_{14}{\omega }^{2\epsilon }+V_{15}{\omega }^{\epsilon }+V_{16}=0.$

Let us assume that

(36)$H_1\left(\omega \right)={\omega }^{6\epsilon }+V_{11}{\omega }^{5\epsilon }+V_{12}{\omega }^{4\epsilon }+V_{13}{\omega }^{3\epsilon }+V_{14}{\omega }^{2\epsilon }+V_{15}{\omega }^{\epsilon }+V_{16},$

when $G_3^2<H_3^2$ and $\underset{\omega \to +\infty }{lim}{H}_1\left(\omega \right)=+\infty$, there exists at least one positive ${\omega }_0$ satisfying $H_1\left({\omega }_0\right)=0$. Suppose that Equation (36) has six positive roots $\omega ={\omega }_k(i=1,2,3,4,5,6)$, with the largest positive root being defined as ${\omega }_0$.

According to (33), we can deduce that the bifurcation point is

(37)${\tau }_{1k}^{\left(j\right)}={\displaystyle \frac{1}{{\omega }_k}}\left(arccos{\displaystyle \frac{{\gamma }_{11}{\eta }_{11}+{\gamma }_{12}{\eta }_{12}}{{{\eta }_{11}}^2+{{\eta }_{12}}^2}}{|}_{\omega ={\omega }_k}+2j\pi \right),$

where $k=1,2,3,4,5,6$ and $j=0,1,2,\dots$.

We define

(38)${\tau }_0=min\{{\tau }_{11}^0,{\tau }_{12}^0,{\tau }_{13}^0,{\tau }_{14}^0,{\tau }_{15}^0,{\tau }_{16}^0\}, {\omega }_0{=\omega |}_{\tau ={\tau }_0}.$

Following this, to validate the transversality condition, we differentiate both sides of Equation (29) with respect to $\tau$

(39)$\begin{array}{cc}\hfill {\left[{\displaystyle \frac{\mathrm{d}s\left(\tau \right)}{\mathrm{d}\tau }}\right]}_{\tau ={\tau }_{1k}^{\left(j\right)}}^{-1}& ={\displaystyle \frac{3\epsilon s^{3\epsilon -1}+2\epsilon G_1{s}^{2\epsilon -1}+\epsilon G_2{s}^{\epsilon -1}}{s(H_1{e}^{2\epsilon }+H_2{e}^{\epsilon }+H_3)e^{-s\tau }}}+{\displaystyle \frac{2\epsilon H_1{s}^{2\epsilon -1}+\epsilon H_2{s}^{\epsilon -1}}{s(H_1{e}^{2\epsilon }+H_2{e}^{\epsilon }+H_3)}}-{\displaystyle \frac{\tau }{s}}.\hfill \end{array}$

By replacing the expression with $s=i\omega$, the following result can be obtained:

(40)$\begin{array}{c}\hfill Re{\left[{\displaystyle \frac{\mathrm{d}s\left(\tau \right)}{\mathrm{d}\tau }}\right]}_{{\tau }_1={\tau }_{1k}^{\left(j\right)}}^{-1}={\displaystyle \frac{O_1{O}_3+O_2{O}_4}{O_1^2+O_2^2}}.\end{array}$

Clearly, $sign\left(Re{\left[{\displaystyle \frac{\mathrm{d}s\left(\tau \right)}{\mathrm{d}\tau }}\right]}_{\tau ={\tau }_{1k}^{\left(j\right)}}^{-1}\right)=sign(O_1{O}_3+O_2{O}_4)$. If $O_1{O}_3+O_2{O}_4>0$ holds, then $Re{\left[{\displaystyle \frac{\mathrm{d}s\left(\tau \right)}{\mathrm{d}\tau }}\right]}_{\tau ={\tau }_{1k}^{\left(j\right)}}^{-1}>0$. Therefore, the transversality condition for generating Hopf bifurcation is met, resulting in a Hopf bifurcation at the equilibrium point $E(x^{*},y^{*},z^{*})$ of System (11) when $\tau ={\tau }_0$. □

5. Chaos Control

In this section, we utilize a fractional-order PID controller to regulate the Hopf bifurcation of a three-population food chain model with Caputo fractional derivatives.

(41)$\begin{array}{c}\hfill u_1\left(t\right)=k_p{e}_1\left(t\right)+k_i {}_{0}J_t^{-\epsilon }{e}_1\left(t\right)+k_d {}_0{}^{C}D_t^{\epsilon }{e}_1\left(t\right)\\ \hfill u_2\left(t\right)=k_p{e}_2\left(t\right)+k_i {}_{0}J_t^{-\epsilon }{e}_2\left(t\right)+k_d {}_0{}^{C}D_t^{\epsilon }{e}_2\left(t\right),\end{array}$

where $e_1\left(t\right)=y-y^{*}$ and $e_2\left(t\right)=z-z^{*}$, and $k_p$, $k_i$ and $k_d$ are the control gains.

Remark 3.

When choosing $\epsilon =1$, $k_p=k_i=k_d=0$, Controller (41) degrades to integer-order controller. When $k_i=k_d=0$, Controller (41) transforms into a linear feedback controller. When $k_p=k_i=0$, Controller (41) turns into a neutral delay feedback controller. When $k_i=0$, Controller (41) degrades to a fractional-order PD controller. When $k_d=0$, Controller (41) reduces to a fractional-order PI controller. When $\epsilon =1$, Controller (41) becomes an integer-order PID controller.

Remark 4.

In practice, we have the flexibility to use single or multiple controllers depending on the specific characteristics and control requirements of the system.

Applying Controller (41) to System (12) results in the following controlled system:

(42)$\begin{array}{cc}\hfill {}_0{}^{C}D_t^{\epsilon }x\left(t\right)& =x\left(t\right)\left[(1-x\left(t\right))·{\displaystyle \frac{1+e_{1}x\left(t\right)}{1+e_{1}x\left(t\right)+f_{1}y\left(t\right)}}-{\displaystyle \frac{{\alpha }_{1}y\left(t\right)}{1+{\beta }_{1}x\left(t\right)}}\right]\hfill \\ \hfill {}_0{}^{C}D_t^{\epsilon }y\left(t\right)& ={\displaystyle \frac{{\alpha }_{1}x\left(t\right)y\left(t\right)}{1+{\beta }_{1}x\left(t\right)}}·{\displaystyle \frac{1+e_{2}y\left(t\right)}{1+e_{2}y\left(t\right)+f_{2}z\left(t\right)}}-{\delta }_{1}y\left(t\right)-{\displaystyle \frac{{\alpha }_{2}y\left(t\right)z\left(t\right)}{1+{\beta }_{2}y\left(t\right)}}+u_1\left(t\right)\hfill \\ \hfill {}_0{}^{C}D_t^{\epsilon }z\left(t\right)& ={\displaystyle \frac{{\alpha }_{2}y\left(t\right)z\left(t\right)}{1+{\beta }_{2}y\left(t\right)}}-{\delta }_{2}z\left(t\right)+u_2\left(t\right).\hfill \end{array}$

Introduce the new variables as follows:

$\begin{array}{c}\hfill v\left(t\right)={}_{0}J_t^{-\epsilon }{e}_1\left(t\right)\\ \hfill w\left(t\right)={}_{0}J_t^{-\epsilon }{e}_2\left(t\right).\end{array}$

From (4), we can derive the following:

$\begin{array}{c}\hfill {}_0{}^{C}D_t^{\epsilon }v\left(t\right)=e_1\left(t\right)\\ \hfill {}_0{}^{C}D_t^{\epsilon }w\left(t\right)=e_2\left(t\right).\end{array}$

Therefore, the following fractional controlled system is obtained:

(43)$\begin{array}{cc}\hfill {}_0{}^{C}D_t^{\epsilon }x\left(t\right)& =x\left(t\right)\left((1-x\left(t\right))·{\displaystyle \frac{1+e_{1}x\left(t\right)}{1+e_{1}x\left(t\right)+f_{1}y\left(t\right)}}-{\displaystyle \frac{{\alpha }_{1}y\left(t\right)}{1+{\beta }_{1}x\left(t\right)}}\right)\hfill \\ \hfill {}_0{}^{C}D_t^{\epsilon }y\left(t\right)& =[{\displaystyle \frac{{\alpha }_{1}x\left(t\right)y\left(t\right)}{1+{\beta }_{1}x\left(t\right)}}·{\displaystyle \frac{1+e_{2}y\left(t\right)}{1+e_{2}y\left(t\right)+f_{2}z\left(t\right)}}-{\delta }_{1}y\left(t\right)-{\displaystyle \frac{{\alpha }_{2}y\left(t\right)z\left(t\right)}{1+{\beta }_{2}y\left(t\right)}}+k_p(y\left(t\right)-y^{*})\hfill \\ & +k_{i}v\left(t\right)]{\displaystyle \frac{1}{1-k_d}}\hfill \\ \hfill {}_0{}^{C}D_t^{\epsilon }z\left(t\right)& =\left[{\displaystyle \frac{{\alpha }_{2}y\left(t\right)z\left(t\right)}{1+{\beta }_{2}y\left(t\right)}}-{\delta }_{2}z\left(t\right)+k_p(z\left(t\right)-z^{*})+k_{i}w\left(t\right)\right]{\displaystyle \frac{1}{1-k_d}}\hfill \\ \hfill {}_0{}^{C}D_t^{\epsilon }v\left(t\right)& =y\left(t\right)-y^{*}\hfill \\ \hfill {}_0{}^{C}D_t^{\epsilon }w\left(t\right)& =z\left(t\right)-z^{*}\hfill \end{array}$

Remark 5.

It is evident that $(x^{*},y^{*},z^{*},0,0)$ represents a favorable equilibrium point of Controlled System (43). It is worth mentioning that the first three components of the vector $(x^{*},y^{*},z^{*},0,0)$ consist of $x^{*}$, $y^{*}$, and $z^{*}$, which serve as positive equilibrium points for Uncontrolled System (12). Therefore, our proposed fractional PID controller, Controller (41), guarantees that the equilibrium point of System (43) remains unchanged.

Stability and Hopf Bifurcation of a Controlled System

The characteristic equation of Equation (43) is derived by

(44)$\left|\begin{array}{ccccc}{s}^{\epsilon }-a_{11}& -a_{12}& 0& 0& 0\\ -{\displaystyle \frac{a_{21}}{1-k_d}}& s^{\epsilon }-{\displaystyle \frac{k_p+a_{22}}{1-k_d}}& -{\displaystyle \frac{a_{23}}{1-k_d}}& {\displaystyle \frac{-k_i}{1-k_d}}& 0\\ 0& -{\displaystyle \frac{a_{32}}{1-k_d}}& s^{\epsilon }-{\displaystyle \frac{k_p+a_{33}}{1-k_d}}& 0& {\displaystyle \frac{-k_i}{1-k_d}}\\ 0& -1& 0& s^{\epsilon }& 0\\ 0& 0& -1& 0& s^{\epsilon }\end{array}\right|,$

where we have

(45)$s^{5\epsilon }+m_1{s}^{4\epsilon }+m_2{s}^{3\epsilon }+m_3{s}^{2\epsilon }+m_4{s}^{\epsilon }+m_5=0.$

Lemma 2.

The control parameters for the uncontrolled system parameters’ selection formula, Formula (46), are $k_p=-0.2$, $k_i=-0.2$ and $k_d=0.2$. Currently, the non-time-delay system exhibits asymptotic stability at the equilibrium point $(0.7628,0.125,5.2434)$.

Proof. 

By substituting Formula (46) and control parameters $k_p=-0.2$, $k_i=-0.2$ and $k_d=0.2$ into (45), we get $m_1>0$, $m_2>0$, $m_3>0$, $m_4>0$ and $m_5>0$. Therefore, all the roots of Formula (45) are situated in the left half-plane of the complex plane. In this case, the system is asymptotically stable at the equilibrium point $(0.7628,0.125,5.2434)$ without any delay. □

Lemma 3.

In the case of the uncontrolled system parameter selection formula, Formula (46), we choose the control parameters $k_p=k_i=k_d=-0.2<0$. Substituting the above parameters into (45), we can obtain $m_1>0$, $m_2>0$, $m_3>0$, $m_4>0$ and $m_5>0$. The fractional-order controlled system, System (43), has a locally asymptotically stable equilibrium, as displayed in Figure 11.

Proof. 

The proof follows the same steps as Lemma 2. □

6. Numerical Simulation

In this paper, we perform numerical simulations using MATLAB R2021a. The predictor–corrector scheme described in [8] is used to numerically simulate a fractional-order system with delay. According to [37], we choose the parameters’ values as follows:

(46)$\begin{array}{cc}\hfill {\alpha }_1& =5, {\alpha }_2=0.1, {\beta }_1=3, {\beta }_2=2, {\delta }_1=0.4,\hfill \\ \hfill {\delta }_2& =0.01, f_1=3.5, f_2=0.1, e_1=1, e_2=2.1.\hfill \end{array}$

Figure 2 shows that $\epsilon$ can affect the stability of the system, being characterized by long memory and heritability.

Using the above parameters, we can calculate $E_3(0.7628,0.125,5.2434)$. According to Theorem 1, the eigenvalues are ${\lambda }_1=-0.3484$, ${\lambda }_2=0.0099+0.0831i$ and ${\lambda }_3=0.0099-0.0831i$, when $\epsilon <0.9245$. Therefore, System (12) is locally asymptotically stable at $\epsilon =0.9$, as shown in Figure 3. For System (12) with the initial value of (0.7, 0.1, 6), the equilibrium point $E_3(0.7628,0.125,5.2434)$ is unstable when $\epsilon =0.95$, as displayed in Figure 4.

Remark 6.

Fractional derivatives can be used to characterize the dynamics of population size. In ecosystems, the dynamics of population size are often influenced by a variety of factors, including climate, food availability, and the number of natural enemies. These factors are often characterized by cyclical or seasonal variations. Therefore, when these cyclical factors are taken into account, the dynamic changes in population size may take on the characteristics of periodicity, i.e., there is a periodic solution.

We illustrate the impact of time-delay $\tau$ on system stability through two examples. The first example is that we choose the stable parameter $e_2=1$ and $\epsilon =0.99$, and the remaining parameters are the same as in (46). Thus, we can obtain $E_3(0.7628,0.125,5.0228)$. Based on Theorem 3, the value of ${\omega }_0$ is estimated to be approximately 0.0799. Consequently, the value of ${\tau }_0$ is found to be nearly 2.15. The system is locally asymptotic stability when $\tau =1.5<{\tau }_0$, as show in Figure 5. At the same time, Figure 6 shows that the system is unstable when $\tau =2.5>{\tau }_0$. The second example is that we choose the parameters in (46) and $\epsilon =0.9$. The value of ${\omega }_0$ is determined to be around 0.0634, while the critical delay $\tau$ is estimated to be around 1.21. As illustrated in Figure 7, the system is locally asymptotic stability when $\tau =1.1<{\tau }_0$. Figure 8 shows that the system is unstable when $\tau =1.5>{\tau }_0$. As depicted in Figure 9, the decrease in the fractional order is accompanied by a gradual increase in the delay threshold ${\tau }_0$.

We choose the control parameters $k_d=0.2$, $k_i=-0.2$ and $k_p=-0.2$, and we choose the remaining parameters in accordance with Equation (46). The control system satisfies $m_1>0$, $m_2>0$, $m_3>0$, $m_4>0$ and $m_5>0$, resulting in local asymptotic stability at the internal equilibrium point $E_3(0.7628,0.125,5.2434)$, as shown in Figure 10. The control system attains local asymptotically stable conditions at the internal equilibrium point $E_3(0.7628,0.125,5.2434)$ through the selection of the control parameters $k_d=-0.2$, $k_i=-0.2$ and $k_p=-0.2$, as illustrated in Figure 11.

According to the previous theoretical analysis, it is evident that the fractional order significantly influences the bifurcation point ${\tau }_0$. Reasonably adjusting the controller parameters $k_p$, $k_i$ and $k_d$ can effectively advance or delay the bifurcation point when the order of $\epsilon$ is consistent. By choosing the controller parameters mentioned above, we are able to stabilize the unstable trajectory at the equilibrium point $E_3(0.7628,0.125,5.2434)$.

We choose the initial values $\left(x\right(0),y(0),z(0),v(0),w(0\left)\right)=(0.7,0.1,6,0.1,6)$. As shown in Figure 12, when $k_i$ and $k_d$ are fixed, the convergence speed of Controlled System (43) increases as $k_p$ decreases, while the amplitude decreases. Figure 13 shows that when $k_p=-2$ and $k_d=-2$, by reducing the integral gain $k_i$, the amplitude of oscillation is decreased. When $k_p=-2,k_i=0$, the oscillations of Controlled System (43) disappear, and with the decrease in the differential gain $k_d$, the convergence speed may also decrease. This is illustrated in Figure 14.

7. Research Limitations and Future Work

The model used in this study assumes more ideal conditions, such as specific changes in control gains, which may ignore more complex factors in real-world ecosystems, such as environmental fluctuations and non-linear relationships between species. Although the model proves the effectiveness of PID control strategies in complex ecosystems, its scope of application may be limited by specific ecological conditions and parameter ranges. Some of the analyses rely on theoretical derivation and lack detailed data support in actual ecosystems, such as dynamic change data of specific populations and monitoring data of environmental factors. In practice, determining the parameters of a PID controller, such as the proportional control gain $k_p$, can be difficult, as the parameters in the ecosystem are often difficult to accurately measure and predict. Although the stability conditions of the system are analyzed, the stability of the system may be disturbed by a variety of factors in practical applications due to the complexity and uncertainty of the ecosystem, resulting in a poor control effect.

Future research can aim to further improve the model and consider more practical factors, such as environmental fluctuations, species adaptation, etc., to improve the accuracy and practicability of the model; strengthen the collection and analysis of data in the actual ecosystem to provide more reliable data support for models; and strengthen the collection and analysis of data in the actual ecosystem to provide more reliable data support for models. This research is innovative at the theoretical level, but it may face a variety of challenges and limitations in practical applications. Thus, it needs to be further improved and expanded in future research.

8. Conclusions

In the case of non-delay system, System (12), we investigated its local and global stability conditions, providing the local stability condition and the global stability of the system using the Lyapunov function method. In the case of System (11) with time delay, we tested the effect of pregnancy delay on its stability by varying the delay parameter $\tau$. Theorem 3 gives us the conditions under which the system is stable. By selecting appropriate parameter values, a threshold ${\tau }_0$ can be obtained if $\tau$ increases continuously. If $\tau <{\tau }_0$, the system is stable. Otherwise, it becomes unstable, resulting in Hopf bifurcation at $\tau ={\tau }_0$. These observations confirm the importance of delay in this system. We also find that the critical value of delay gradually decreases as the fractional-order value increases.

The experimental results show that selecting appropriate control parameters in the PID controller can change the system’s convergence rate and vibration amplitude. It also demonstrates the feasibility of using PID controllers in ecological systems.

Author Contributions

Y.Z.: conceptualization, data curation and funding acquisition; Z.C.: formal analysis and writing—original draft preparation; R.L.: writing—review and editing. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Footnotes

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.

Figures and Table

Figure 1. Stability region of Fractional System (5) with [Forumla omitted. See PDF.].

Figure 2. The time series (a–c) and phase diagram (d) for System (12) when [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.], where [Forumla omitted. See PDF.]. Other parameters are the same as in (46).

Figure 3. The time series (a) and phase diagram (b) for System (12) when the initial value [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. System (12) is locally asymptotically stable at the equilibrium point [Forumla omitted. See PDF.]. Other parameters are the same as in (46).

Figure 4. The time series (a) and phase diagram (b) of the interior equilibrium point [Forumla omitted. See PDF.] result in periodic oscillations for System (12) with the initial value [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. Other parameters are same as in (46).

Figure 5. The time series (a) and phase diagram (b) of the interior equilibrium point [Forumla omitted. See PDF.] are locally asymptotically stable for System (11) with the initial value [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 6. The time series (a) and phase diagram (b) of the interior equilibrium point [Forumla omitted. See PDF.] are locally asymptotically stable for System (11) with the initial value [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.].

Figure 7. The time series (a) and phase diagram (b) of the interior equilibrium point [Forumla omitted. See PDF.] are locally asymptotically stable for System (11) with the initial value [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 8. The time series (a) and phase diagram (b) of the interior equilibrium point [Forumla omitted. See PDF.] are locally asymptotically stable for System (11) with the initial value [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.].

Figure 9. Relation of fractional-order [Forumla omitted. See PDF.] to [Forumla omitted. See PDF.].

Figure 10. The time series (a) and phase diagram (b) of the interior equilibrium point [Forumla omitted. See PDF.] are locally asymptotically stable for Control System (43) with the initial value [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 11. The time series (a) and phase diagram (b) of the interior equilibrium point [Forumla omitted. See PDF.] are locally asymptotically stable for the Control System (43) with the initial value [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 12. The time series (a–c) of controlled System (43) with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and different values of [Forumla omitted. See PDF.].

Figure 13. The time series (a–c) of controlled System (43) with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and different values of [Forumla omitted. See PDF.].

Figure 14. The time series (a–c) of controlled System (43) with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and different values of [Forumla omitted. See PDF.].

Table 1

Distribution of (18) roots.

${\Omega }_1$	${\lambda }_1<0$, ${\lambda }_2<0$, ${\lambda }_3<0$
${\Omega }_2$	${\lambda }_1>0$, ${\lambda }_2<0$, ${\lambda }_3<0$
${\Omega }_3$	${\lambda }_1<0$, ${\lambda }_2>0$, ${\lambda }_3>0$
${\Omega }_4$	${\lambda }_1>0$, ${\lambda }_2>0$, ${\lambda }_3>0$
${\Omega }_5$	${\lambda }_1<0$, ${\lambda }_{2,3}=a\pm b\mathrm{i}(a<0)$
${\Omega }_6$	${\lambda }_1<0$, ${\lambda }_{2,3}=a\pm b\mathrm{i}(a>0)$
${\Omega }_7$	${\lambda }_1>0$, ${\lambda }_{2,3}=a\pm b\mathrm{i}(a<0)$

Appendix A

$$\begin{array}{cc}\hfill x_1^{*}& ={\displaystyle \frac{{\delta }_1}{{\alpha }_1-{\delta }_1{\beta }_1}}\hfill \\ \hfill y_1^{*}& ={\displaystyle \frac{\sqrt{{(1+e_1{x}_1^{*})}^2+4{\alpha }_1^{-1}{f}_1(1-x_1^{*})(1+e_1{x}_1^{*})(1+{\beta }_1{x}_1^{*})}}{2f_1}}+{\displaystyle \frac{-(1+e_1{x}_1^{*})}{2f_1}}\hfill \end{array}$$$x_2^{*}$ satisfies $e_1{\beta }_1{x_2^{*}}^3+k_1{x_2^{*}}^2+k_2{x}_2^{*}+k_3=0$, $k_1=e_1+{\beta }_1-e_1{\beta }_1$, $k_2={\alpha }_1{e}_1{y}_2^{*}+1-e_1-{\beta }_1$ and $k_3={\alpha }_1{y}_2^{*}(1+f_1{y}_2^{*})-1$, $y_2^{*}={\displaystyle \frac{{\delta }_2}{{\alpha }_2-{\beta }_2{\delta }_2}}$, while $z_2^{*}$ satisfies $$\begin{array}{cc}& -(1+e_2{y}_2^{*})\left({\displaystyle \frac{{\alpha }_1{x}_2^{*}}{1+{\beta }_1{x}_2^{*}}}-{\delta }_1\right)+\left({\displaystyle \frac{{\alpha }_2(1+e_2{y}_2^{*})}{1+{\beta }_2{y}_2^{*}}}+{\delta }_1{f}_2\right)z_2^{*}+{\displaystyle \frac{{\alpha }_2{f}_2}{1+{\beta }_2{y}_2^{*}}}{z_2^{*}}^2=0.\hfill \\ \hfill A_{11}& ={\displaystyle \frac{1-2x+2e_{1}x-3e_1{x}^2}{1+e_{1}x+f_{1}y}}-{\displaystyle \frac{e_{1}x(1-x)(1+e_{1}x)}{{(1+e_{1}x+f_{1}y)}^2}}-{\displaystyle \frac{{\alpha }_{1}y}{{(1+{\beta }_{1}x)}^2}}\hfill \\ \hfill A_{12}& =-{\displaystyle \frac{f_{1}x(1-x)(1+e_{1}x)}{{(1+e_{1}x+f_{1}y)}^2}}-{\displaystyle \frac{{\alpha }_{1}x}{1+{\beta }_{1}x}}, A_{13}=0\hfill \\ \hfill A_{21}& ={\displaystyle \frac{{\alpha }_1(1+e_{2}y)y}{(1+e_{2}y+f_{2}z){(1+{\beta }_{1}x)}^2}}\hfill \\ \hfill A_{22}& ={\displaystyle \frac{{\alpha }_{1}x({(1+e_{2}y)}^2+f_{2}z(1+2e_{2}y))}{(1+{\beta }_{1}x)(1+e_{2}y+f_{2}z)}}-{\displaystyle \frac{{\alpha }_{2}z}{{(1+{\beta }_{2}y)}^2}}-{\delta }_1\hfill \\ \hfill A_{23}& =-{\displaystyle \frac{{\alpha }_1{f}_2(1+e_{2}y)xy}{(1+{\beta }_{1}x){(1+e_{2}y+f_{2}z)}^2}}-{\displaystyle \frac{{\alpha }_{2}y}{1+{\beta }_{2}y}}, A_{31}=0, A_{32}={\displaystyle \frac{{\alpha }_{2}z}{{(1+{\beta }_{2}y)}^2}}\hfill \\ \hfill A_{33}& ={\displaystyle \frac{{\alpha }_{2}y}{1+{\alpha }_{2}y}}-{\delta }_2, B_1={\displaystyle \frac{1-2x_1^{*}+2e_1{x}_1^{*}-3e_1{x_1^{*}}^2}{1+e_1{x}_1^{*}+f_1{y}_1^{*}}}-{\displaystyle \frac{e_1{x}_1^{*}(1-x_1^{*})(1+e_1{x}_1^{*})}{{(1+e_1{x}_1^{*}+f_1{y}_1^{*})}^2}}\hfill \\ & -{\displaystyle \frac{{\alpha }_1{y}_1^{*}}{{(1+{\beta }_1{y}_1^{*})}^2}}\hfill \\ \hfill B_2& ={\displaystyle \frac{f_1{x}_1^{*}(1-x_1^{*})(1+e_1{x}_1^{*})}{{(1+e_1{x}_1^{*}+f_1{y}_1^{*})}^2}}+{\displaystyle \frac{{\alpha }_1{x}_1^{*}}{1+{\beta }_1{y}_1^{*}}}, B_3={\displaystyle \frac{{\alpha }_1{y}_1^{*}}{{(1+{\beta }_1{x}_1^{*})}^2}},\hfill \\ \hfill B_4& ={\displaystyle \frac{{\alpha }_1{f}_2{x}_1^{*}{y}_1^{*}}{(1+e_2{y}_1^{*})(1+{\beta }_1{x}_1^{*})}}+{\displaystyle \frac{{\alpha }_2{y}_1^{*}}{(1+{\beta }_2{y}_1^{*})}}, B_5={\displaystyle \frac{{\alpha }_2{y}_1^{*}}{1+{\beta }_2{y}_1^{*}}}-{\delta }_2\hfill \\ \hfill C_1& ={\displaystyle \frac{1-2x_2^{*}+2e_1{x}_2^{*}-3e_1{x_2^{*}}^2}{1+e_1{x}_2^{*}+f_1{y}_2^{*}}}-{\displaystyle \frac{e_1{x}_2^{*}(1-x_2^{*})(1+e_1{x}_2^{*})}{{(1+e_1{x}_2^{*}+f_1{y}_2^{*})}^2}}-{\displaystyle \frac{{\alpha }_1{y}_2^{*}}{{(1+{\beta }_1{x}_2^{*})}^2}}\hfill \\ \hfill C_2& ={\displaystyle \frac{f_1{x}_2^{*}(1-x_2^{*})(1+e_1{x}_2^{*})}{{(1+e_1{x}_2^{*}+f_1{y}_2^{*})}^2}}+{\displaystyle \frac{{\alpha }_1{x}_2^{*}}{1+{\beta }_1{x}_2^{*}}}, C_3={\displaystyle \frac{{\alpha }_1(1+e_2{y}_2^{*})y_2^{*}}{(1+e_2{y}_2^{*}+f_2{z}_2^{*}){(1+{\beta }_1{x}_2^{*})}^2}}\hfill \\ \hfill C_4& ={\displaystyle \frac{{\alpha }_1{e}_2{f}_2{x}_2^{*}{y}_2^{*}{z}_2^{*}}{(1+{\beta }_1{x}_2^{*}){(1+e_2{y}_2^{*}+f_2{z}_2^{*})}^2}}+{\displaystyle \frac{{\alpha }_2{\beta }_2{y}_2^{*}{z}_2^{*}}{{(1+{\beta }_2{y}_2^{*})}^2}}\hfill \\ \hfill C_5& ={\displaystyle \frac{{\alpha }_1{f}_2(1+e_2{y}_2^{*})x_2^{*}{y}_2^{*}}{(1+{\beta }_1{x}_2^{*}){(1+e_2{y}_2^{*}+f_2{z}_2^{*})}^2}}+{\displaystyle \frac{{\alpha }_2{y}_2^{*}}{1+{\beta }_2{y}_2^{*}}}, C_6={\displaystyle \frac{{\alpha }_2{z}_2^{*}}{{(1+{\beta }_2{y}_2^{*})}^2}}\hfill \\ \hfill a_{11}& ={\displaystyle \frac{1-2x^{*}+2e_1{x}^{*}-3e_1{x^{*}}^2}{1+e_1{x}^{*}+f_1{y}^{*}}}-{\displaystyle \frac{e_1{x}^{*}(1-x^{*})(1+e_1{x}^{*})}{{(1+e_1{x}^{*}+f_1{y}^{*})}^2}}-{\displaystyle \frac{{\alpha }_1{y}^{*}}{{(1+{\beta }_1{x}^{*})}^2}},\hfill \\ \hfill a_{12}& ={\displaystyle \frac{f_1[-e_1{x^{*}}^3+(e_1-1){x^{*}}^2+x^{*}]}{{(1+e_1{x}^{*}+f_1{y}^{*})}^2}}-{\displaystyle \frac{{\alpha }_1{x}^{*}}{1+{\beta }_1{x}^{*}}},\hfill \\ \hfill a_{21}& ={\displaystyle \frac{{\alpha }_1{y}^{*}}{{(1+{\beta }_1{x}^{*})}^2}}·{\displaystyle \frac{1+e_2{y}^{*}}{1+e_2{y}^{*}+f_2{z}^{*}}},\hfill \\ \hfill a_{22}& ={\displaystyle \frac{{\alpha }_1{x}^{*}(1+e_2{y}^{*})}{(1+{\beta }_1{x}^{*})(1+e_2{y}^{*}+f_2{z}^{*})}}-{\delta }_1-{\displaystyle \frac{{\alpha }_2{z}^{*}}{1+{\beta }_2{y}^{*}}}+{\displaystyle \frac{{\alpha }_1{f}_2{e}_2{x}^{*}{y}^{*}{z}^{*}}{(1+{\beta }_1{x}^{*}){(1+e_2{y}^{*}+f_2{z}^{*})}^2}},\hfill \\ \hfill a_{23}& =-{\displaystyle \frac{{\alpha }_1{x}^{*}{y}^{*}}{1+{\beta }_1{x}^{*}}}{\displaystyle \frac{f_2(1+e_2{y}^{*})}{{(1+e_2{y}^{*}+f_2{z}^{*})}^2}}-{\displaystyle \frac{{\alpha }_2{y}^{*}}{{(1+{\beta }_2{y}^{*})}^2}},\hfill \\ \hfill a_{33}& =-{\delta }_2, b_{32}={\displaystyle \frac{{\alpha }_2{z}^{*}}{{(1+{\beta }_2{y}^{*})}^2}} b_{33}={\displaystyle \frac{{\alpha }_2{y}^{*}}{1+{\beta }_2{y}^{*}}}.\hfill \\ \hfill V_{11}& =2G_{1}cos{\displaystyle \frac{\pi \epsilon }{2}}, V_{12}=G_1^2+2G_{2}cos\pi \epsilon -H_1^2\hfill \\ \hfill V_{13}& =2G_{3}cos{\displaystyle \frac{3\epsilon \pi }{2}}-2H_1{H}_{2}cos{\displaystyle \frac{\pi \epsilon }{2}}+2G_1{G}_{2}cos{\displaystyle \frac{\epsilon \pi }{2}}\hfill \\ \hfill V_{14}& =G_2^2+2G_1{G}_{3}cos\epsilon \pi -H_2^2-2H_1{H}_{3}cos\epsilon \pi \hfill \\ \hfill V_{15}& =2G_2{G}_{3}cos{\displaystyle \frac{\epsilon \pi }{2}}-2H_2{H}_{3}cos{\displaystyle \frac{\epsilon \pi }{2}}, V_{16}=G_3^2-H_3^2\hfill \\ \hfill O_1& =(H_1{\omega }_0^{2\epsilon }cos\epsilon \pi +H_2{\omega }_0^{\epsilon }cos{\displaystyle \frac{\epsilon \pi }{2}}+H_3)cos{\omega }_0\tau +sin{\omega }_0\tau (H_1{\omega }_0^{2\epsilon }sin\epsilon \pi +H_2{\omega }_0^{\epsilon }sin{\displaystyle \frac{\epsilon \pi }{2}})\hfill \\ \hfill O_2& =-sin{\omega }_0\tau (H_1{\omega }_0^{2\epsilon }cos\epsilon \pi +H_2{\omega }_0^{\epsilon }cos{\displaystyle \frac{\epsilon \pi }{2}}+H_3)+cos{\omega }_0\tau (H_1{\omega }_0^{2\epsilon }sin\epsilon \pi +H_2{\omega }_0^{\epsilon }sin{\displaystyle \frac{\epsilon \pi }{2}})\hfill \\ \hfill O_3& =3\epsilon {\omega }_0^{3\epsilon -2}cos{\displaystyle \frac{(3\epsilon -1)\pi }{2}}+2\epsilon G_1{\omega }_0^{2\epsilon -2}cos{\displaystyle \frac{(2\epsilon -1)\pi }{2}}+\epsilon G_2{\omega }_0^{\epsilon -2}cos{\displaystyle \frac{(\epsilon -1)\pi }{2}}\hfill \\ & +2\epsilon H_1{\omega }_0^{2\epsilon -2}cos({\omega }_0\tau -{\displaystyle \frac{(2\epsilon -1)\pi }{2}})+\epsilon H_2{\omega }_0^{\epsilon -2}cos({\omega }_0\tau -{\displaystyle \frac{(\epsilon -1)\pi }{2}})\hfill \\ \hfill O_4& =3\epsilon {\omega }_0^{3\epsilon -2}sin{\displaystyle \frac{(3\epsilon -1)\pi }{2}}+2\epsilon G_1{\omega }_0^{2\epsilon -2}sin{\displaystyle \frac{(2\epsilon -2)\pi }{2}}+\epsilon G_2{\omega }_0^{\epsilon -2}sin{\displaystyle \frac{(\epsilon -1)\pi }{2}}\hfill \\ & +2\epsilon H_1{\omega }_0^{2\epsilon -2}sin({\omega }_0\tau -{\displaystyle \frac{(2\epsilon -1)\pi }{2}})+\epsilon H_2{\omega }_0^{\epsilon -2}sin({\omega }_0\tau -{\displaystyle \frac{(\epsilon -1)\pi }{2}})\hfill \\ \hfill m_1& =[-a_{11}-a_{22}-a_{33}-2k_p+2a_{11}{k}_d+a_{22}{k}_d+2a_{33}{k}_d\hfill \\ & +3k_p{k}_d-a_{11}{k}_d^2-a_{33}{k}_d^2-k_d^2{k}_p-b_{33}{(1-k_d)}^2]{\displaystyle \frac{1}{{(1-k_d)}^2}}\hfill \\ \hfill m_2& =[k_p^2-2k_i+a_{11}{a}_{22}-a_{12}{a}_{21}+a_{11}{a}_{33}+a_{22}{a}_{33}+2a_{11}{k}_p+a_{22}{k}_p+a_{33}{k}_p+2k_d{k}_i\hfill \\ & -k_d{k}_p^2+a_{11}{k}_d^2{k}_p-a_{11}{a}_{22}{k}_d+a_{12}{a}_{21}{k}_d-2a_{11}{a}_{33}{k}_d-a_{22}{a}_{33}{k}_d-3a_{11}{k}_d{k}_p\hfill \\ & -a_{22}{k}_d{k}_p-a_{33}{k}_d{k}_p+a_{11}{a}_{33}{k}_d^2+a_{11}{b}_{33}+a_{22}{b}_{33}-a_{23}{b}_{32}+b_{33}{k}_p-2a_{11}{b}_{33}{k}_d\hfill \\ & -a_{22}{b}_{33}{k}_d-b_{33}{k}_d{k}_p+a_{11}{b}_{33}{k}_d^2]{\displaystyle \frac{1}{{(1-k_d)}^2}}\hfill \\ \hfill m_3& =[2a_{11}{k}_i+a_{22}{k}_i+a_{33}{k}_i+2k_i{k}_p-a_{11}{k}_p^2+a_{11}{k}_d{k}_p^2-a_{11}{a}_{22}{a}_{33}+a_{12}{a}_{21}{a}_{33}\hfill \\ & -a_{11}{a}_{22}{k}_p+a_{12}{a}_{21}{k}_p-a_{11}{a}_{33}{k}_p-2a_{11}{k}_d{k}_i-a_{33}{k}_d{k}_i-k_d{k}_i{k}_p+a_{11}{a}_{22}{k}_d{k}_p\hfill \\ & -a_{12}{a}_{21}{k}_d{k}_p+a_{11}{a}_{33}{k}_d{k}_p+a_{11}{a}_{22}{a}_{33}{k}_d-a_{12}{a}_{21}{a}_{33}{k}_d+b_{33}{k}_i-2a_{11}{a}_{22}{b}_{33}+a_{11}{a}_{23}{b}_{32}\hfill \\ & +a_{12}{a}_{21}{b}_{33}-a_{11}{b}_{33}{k}_p-b_{33}{k}_d{k}_i+a_{11}{a}_{22}{b}_{33}{k}_d-a_{12}{a}_{21}{b}_{33}{k}_d+a_{11}{b}_{33}{k}_d{k}_p\left)\right]{\displaystyle \frac{1}{{(1-k_d)}^2}}\hfill \\ \hfill m_4& =[k_i^2-a_{11}{a}_{22}{k}_i+a_{12}{a}_{21}{k}_i-a_{11}{a}_{33}{k}_i-2a_{11}{k}_i{k}_p+a_{11}{b}_{33}{k}_d{k}_i+a_{11}{k}_d{k}_i{k}_p\hfill \\ & +(a_{11}{b}_{33}{k}_d{k}_i-a_{11}{b}_{33}{k}_i)]{\displaystyle \frac{1}{{(1-k_d)}^2}}\hfill \\ \hfill m_5& ={\displaystyle \frac{-a_{11}{k}_i^2}{{(1-k_d)}^2}}\hfill \end{array}$$

References

1. Holling, C.S. The components of predation as revealed by a study of small-mammal predation of the European Pine Sawfly. Can. Entomol.; 1959; 91, pp. 293-320. [DOI: https://dx.doi.org/10.4039/Ent91293-5]

2. Hastings, A.; Powell, T. Chaos in a three-species food chain. Ecology; 1991; 72, pp. 896-903. [DOI: https://dx.doi.org/10.2307/1940591]

3. Zheng, Z.; Wang, K.; Yang, D.; Yin, F.; Sun, D.; Yu, W.; Lin, J.; Liu, Y.; Chen, C.; Yang, Z. et al. Factors affecting the transmission of dengue fever in Haikou city in 2019. Math. Biosci. Eng.; 2023; 20, pp. 16045-16059. [DOI: https://dx.doi.org/10.3934/mbe.2023716]

4. Mishra, A.K.; Kumar, R.; Yadav, V.K. Stability Analysis and Controlling Chaos of Fractional-Order Three-Species Food Chain Model with Fear. Int. J. Appl. Comput. Math.; 2021; 7, 214. [DOI: https://dx.doi.org/10.1007/s40819-021-01150-y]

5. Alidousti, J.; Mostafavi Ghahfarokhi, M. Dynamical behavior of a fractional three-species food chain model. Nonlinear Dyn.; 2019; 95, pp. 1841-1858. [DOI: https://dx.doi.org/10.1007/s11071-018-4663-6]

6. Wang, B.; Sajjadi, S.S.; Jahanshahi, H.; Karaca, Y.; Hou, D.; Pi, L.; Xia, W.F.; Aly, A.A. Predictive control of the variable-order fractional chaotic ecological system. Fractals; 2022; 30, 2240178. [DOI: https://dx.doi.org/10.1142/S0218348X22401788]

7. Mandelbrot, B.B. The Fractal Geometry of Nature. Am. J. Phys.; 1983; 51, pp. 286-287.

8. Bhalekar, S.; Daftardar-Gejji, V. A predictor-corrector scheme for solving nonlinear delay differential equations of fractional order. J. Fract. Calc. Appl; 2011; 1, pp. 1-9.

9. Xu, C.; Yu, Y.; Ren, G.; Sun, Y.; Si, X. Stability analysis and optimal control of a fractional-order generalized SEIR model for the COVID-19 pandemic. Appl. Math. Comput.; 2023; 457, 128210. [DOI: https://dx.doi.org/10.1016/j.amc.2023.128210]

10. Liu, C.; Chen, Y.; Yu, Y.; Wang, Z. Bifurcation and Stability Analysis of a New Fractional-Order Prey-Predator Model with Fear Effects in Toxic Injections. Mathematics; 2023; 11, 4367. [DOI: https://dx.doi.org/10.3390/math11204367]

11. Xiao, M.; Zheng, W.; Jiang, G.; Cao, J. Undamped oscillations generated by Hopf bifurcations in fractional-order recurrent neural networks with Caputo derivative. IEEE Trans. Neural Netw. Learn. Syst.; 2015; 26, pp. 3201-3214. [DOI: https://dx.doi.org/10.1109/TNNLS.2015.2425734]

12. Ji, G.; Ge, Q.; Xu, J. Dynamic behaviors of a fractional order two-species cooperative systems with harvesting. Chaos Solitons Fractals; 2016; 92, pp. 51-55. [DOI: https://dx.doi.org/10.1016/j.chaos.2016.09.014]

13. Das, S. Functional Fractional Calculus; Springer: Berlin/Heidelberg, Germany, 2011; Volume 1.

14. Tavazoei, M.S.; Haeri, M.; Attari, M.; Bolouki, S.; Siami, M. More details on analysis of fractional-order Van der Pol oscillator. J. Vib. Control; 2009; 15, pp. 803-819. [DOI: https://dx.doi.org/10.1177/1077546308096101]

15. Tavazoei, M.S.; Haeri, M. A proof for non existence of periodic solutions in time invariant fractional order systems. Automatica; 2009; 45, pp. 1886-1890. [DOI: https://dx.doi.org/10.1016/j.automatica.2009.04.001]

16. Aguila-Camacho, N.; Duarte-Mermoud, M.A.; Gallegos, J.A. Lyapunov functions for fractional order systems. Commun. Nonlinear Sci. Numer. Simul.; 2014; 19, pp. 2951-2957. [DOI: https://dx.doi.org/10.1016/j.cnsns.2014.01.022]

17. Banerjee, M.; Takeuchi, Y. Maturation delay for the predators can enhance stable coexistence for a class of prey–predator models. J. Theor. Biol.; 2017; 412, pp. 154-171. [DOI: https://dx.doi.org/10.1016/j.jtbi.2016.10.016]

18. Dubey, B.; Kumar, A.; Maiti, A.P. Global stability and Hopf-bifurcation of prey-predator system with two discrete delays including habitat complexity and prey refuge. Commun. Nonlinear Sci. Numer. Simul.; 2019; 67, pp. 528-554. [DOI: https://dx.doi.org/10.1016/j.cnsns.2018.07.019]

19. Zhang, J.F. Bifurcation analysis of a modified Holling–Tanner predator–prey model with time delay. Appl. Math. Model.; 2012; 36, pp. 1219-1231. [DOI: https://dx.doi.org/10.1016/j.apm.2011.07.071]

20. Hu, D.; Li, Y.; Liu, M.; Bai, Y. Stability and Hopf bifurcation for a delayed predator–prey model with stage structure for prey and Ivlev-type functional response. Nonlinear Dyn.; 2020; 99, pp. 3323-3350. [DOI: https://dx.doi.org/10.1007/s11071-020-05467-z]

21. Bai, Y.; Li, Y. Stability and Hopf bifurcation for a stage-structured predator–prey model incorporating refuge for prey and additional food for predator. Adv. Differ. Equ.; 2019; 2019, pp. 1-20. [DOI: https://dx.doi.org/10.1186/s13662-019-1979-6]

22. Chen, F.; You, M. Permanence for an integrodifferential model of mutualism. Appl. Math. Comput.; 2007; 186, pp. 30-34. [DOI: https://dx.doi.org/10.1016/j.amc.2006.07.085]

23. Beretta, E.; Kuang, Y. Geometric stability switch criteria in delay differential systems with delay dependent parameters. SIAM J. Math. Anal.; 2002; 33, pp. 1144-1165. [DOI: https://dx.doi.org/10.1137/S0036141000376086]

24. Deng, L.; Wang, X.; Peng, M. Hopf bifurcation analysis for a ratio-dependent predator–prey system with two delays and stage structure for the predator. Appl. Math. Comput.; 2014; 231, pp. 214-230. [DOI: https://dx.doi.org/10.1016/j.amc.2014.01.025]

25. Tripathi, J.P.; Abbas, S.; Thakur, M. A density dependent delayed predator–prey model with Beddington–DeAngelis type function response incorporating a prey refuge. Commun. Nonlinear Sci. Numer. Simul.; 2015; 22, pp. 427-450. [DOI: https://dx.doi.org/10.1016/j.cnsns.2014.08.018]

26. Zheng, T.; Zhang, L.; Luo, Y.; Zhou, X.; Li, H.L.; Teng, Z. Stability and Hopf bifurcation of a stage-structured cannibalism model with two delays. Int. J. Bifurc. Chaos; 2021; 31, 2150242. [DOI: https://dx.doi.org/10.1142/S0218127421502424]

27. Wang, W.; Jia, X.; Luo, X.; Kurths, J.; Yuan, M. Fixed-time synchronization control of memristive MAM neural networks with mixed delays and application in chaotic secure communication. Chaos Solitons Fractals; 2019; 126, pp. 85-96. [DOI: https://dx.doi.org/10.1016/j.chaos.2019.05.041]

28. Ouannas, A.; Bendoukha, S.; Volos, C.; Boumaza, N.; Karouma, A. Synchronization of fractional hyperchaotic Rabinovich systems via linear and nonlinear control with an application to secure communications. Int. J. Control. Autom. Syst.; 2019; 17, pp. 2211-2219. [DOI: https://dx.doi.org/10.1007/s12555-018-0216-5]

29. Chen, S.; Yang, J.; Li, Z.; Li, H.L.; Hu, C. New results for dynamical analysis of fractional-order gene regulatory networks with time delay and uncertain parameters. Chaos Soliton Fractals; 2023; 175, 114041. [DOI: https://dx.doi.org/10.1016/j.chaos.2023.114041]

30. Kekha Javan, A.A.; Shoeibi, A.; Zare, A.; Hosseini Izadi, N.; Jafari, M.; Alizadehsani, R.; Moridian, P.; Mosavi, A.; Acharya, U.R.; Nahavandi, S. Design of Adaptive-Robust Controller for Multi-State Synchronization of Chaotic Systems with Unknown and Time-Varying Delays and Its Application in Secure Communication. Sensors; 2021; 21, 254. [DOI: https://dx.doi.org/10.3390/s21010254] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33401741]

31. Xiao, M.; Zheng, W.X.; Lin, J.; Jiang, G.; Zhao, L.; Cao, J. Fractional-order PD control at Hopf bifurcations in delayed fractional-order small-world networks. J. Frankl. Inst.; 2017; 354, pp. 7643-7667. [DOI: https://dx.doi.org/10.1016/j.jfranklin.2017.09.009]

32. Tang, Y.; Xiao, M.; Jiang, G.; Lin, J.; Cao, J.; Zheng, W.X. Fractional-order PD control at Hopf bifurcations in a fractional-order congestion control system. Nonlinear Dyn.; 2017; 90, pp. 2185-2198. [DOI: https://dx.doi.org/10.1007/s11071-017-3794-5]

33. Tao, B.; Xiao, M.; Jiang, G.; Sun, Q. Bifurcation control of small-world networks with delays VIA PID controller. Asian J. Control; 2020; 22, pp. 818-830. [DOI: https://dx.doi.org/10.1002/asjc.1958]

34. Xiao, M.; Tao, B.; Zheng, W.X.; Jiang, G. Fractional-Order PID Controller Synthesis for Bifurcation of Fractional-Order Small-World Networks. IEEE Trans. Syst. Man Cybern. Syst.; 2021; 51, pp. 4334-4346. [DOI: https://dx.doi.org/10.1109/TSMC.2019.2933570]

35. Lu, L.; Huang, C.; Song, X. Bifurcation control of a fractional-order PD control strategy for a delayed fractional-order prey–predator system. Eur. Phys. J. Plus; 2023; 138, pp. 1-11. [DOI: https://dx.doi.org/10.1140/epjp/s13360-023-03708-9]

36. Ramasamy, S.; Banjerdpongchai, D.; Park, P. Chaos Control of a Delayed Tri-Trophic Food Chain Model with Fear and Its Carry Over Effects. Symmetry; 2023; 15, 484. [DOI: https://dx.doi.org/10.3390/sym15020484]

37. Cui, Z.; Zhou, Y.; Li, R. Complex Dynamics Analysis and Chaos Control of a Fractional-Order Three-Population Food Chain Model. Fractal Fract.; 2023; 7, 548. [DOI: https://dx.doi.org/10.3390/fractalfract7070548]

38. Irving, R.S. Integers, Polynomials, and Rings: A Course in Algebra; Springer: Berlin/Heidelberg, Germany, 2004.

39. Das, M.; Samanta, G. A delayed fractional order food chain model with fear effect and prey refuge. Math. Comput. Simul.; 2020; 178, pp. 218-245. [DOI: https://dx.doi.org/10.1016/j.matcom.2020.06.015]