Stability Analysis of a Diffusive Three-Species Ecological System with Time Delays

: In this study, the dynamics of a diffusive Lotka–Volterra three-species system with delays were explored. By employing the Galerkin Method, which generates semi-analytical solutions, a partial differential equation system was approximated through mathematical modeling with delay differential equations. Steady-state curves and Hopf bifurcation maps were created and discussed in detail. The effects of the growth rate of prey and the mortality rate of the predator and top predator on the system’s stability were demonstrated. Increase in the growth rate of prey destabilised the system, whilst increase in the mortality rate of predator and top predator stabilised it. The increase in the growth rate of prey likely allowed the occurrence of chaotic solutions in the system. Additionally, the effects of hunting and maturation delays of the species were examined. Small delay responses stabilised the system, whilst great delays destabilised it. Moreover, the effects of the diffusion coefﬁcients of the species were investigated. Alteration of the diffusion coefﬁcients rendered the system permanent or extinct.


Introduction
Recently, the population dynamics of various biological and ecological systems have garnered much attention, promoting the emergence of mathematical models describing population dynamics and species interactions. These mathematical models must incorporate both spatial diffusion and time delays to mirror the dynamic nature of biological systems and the tendency of the species to move to the least densely populated areas. Reaction-diffusion models with time delays, which show oscillatory phenomena, can describe the delayed response to past conduct and the spatial structure of chemical, biological, and ecological systems [1][2][3][4]. In particular, the delay-diffusive logistic equation describes the growth dynamics of a single species, whilst the delay-diffusive Lotka-Volterra predator-prey model describes the dynamics of multiple species [5][6][7].
In recent years, the predator-prey relationship represents the most studied of environmental dynamics. In the 1920s, Alfred Lotka and Vito Volterra [8] proposed a mathematical model for a predator-prey system, which described the fish catch in the Atlantic. This model can also be utilised to study physical systems and chemical reactions, such as the dynamics of resonantly coupled lasers [8,9]. Several studies have examined the stability of the Lotka-Volterra predator-prey model. For instance, Faria [10] investigated the predator-prey system with one or two time delays and a unique positive equilibrium. They examined the dynamics of this system based on the local stability of the equilibrium and the region of the Hopf bifurcation map that has been proven to occur when one of the delays is assumed to be a bifurcation parameter. Yan and Chu [11] explored the stability of the Lotka-Volterra predator-prey model and found conditions for the occurrence of oscillating solutions. The authors also studied the stability of the oscillating solutions. Chen et al. [12] analysed the diffusive Lotka-Volterra predator-prey model with two delays. The authors predator which consumes both U and V. Through the Lotka-Volterra type of interactions and the effects of diffusion and time delays for each species, this system can be modelled by the following delay reaction-diffusion equations: In the equations above , u, v, and w are the corresponding scaled concentrations of the population densities of the prey, predator, and top predator, respectively. At x = 0, zero-flux boundary conditions are applied; thus, it is an impermeable boundary as it indicates that there is no population flux across the boundary. At x = ±1, the Dirichlet boundary conditions are applied, which indicate that the population is fixed [14,42]. So the solution is symmetric about the center of the domain x = 0. Thus, system (1) and (2) is open and allows for the occurrence of steady-state solutions, periodic oscillations, and chaotic behaviours. Here, the intrinsic growth rate of the prey species is denoted by the parameter α, whilst the mortality rates of the predator and top predator species are represented by the parameters β and γ, respectively. The carrying capacities of the populations of the three species, namely the prey u, predator v, and top predator w, are represented by the parameters δ 1 , µ 2 , and κ 3 , respectively. µ 1 and κ 1 indicate the reduction in the prey population as a result of the presence of the predator and top predator, respectively. δ 2 indicates the increase in the predator population as a result of the presence of prey, whilst κ 2 represents the decrease in the predator population as a result of the presence of the top predator. The parameters δ 3 and µ 3 indicate the increase in the top predator population as a result of the presence of the prey and predator, respectively. τ i (i = 1 . . . 3) refers to hunting and maturation delays. The parameters D i (i = 1 . . . 3) indicate the diffusion coefficients of the three species u, v, and w. For physically realistic population models, all parameters are positive. On the other hand, if the top predator population (w) is zero, (where D 3 = γ = µ 3 = κ i = 0; i = 1 . . . 3), then the system (1) reduces to the two-species prey-predator model, which has been investigated in [14].

Galerkin Method
By applying the Galerkin method, the semi-analytical model (1) was obtained and developed in a one dimensional (1 − D) spatial domain. This method is based on the approximation of the spatial structure of the population density profile using a set of basis functions, converting the governing PDE (1) and boundary conditions (2) for formulation with the simplest possible laws of the fundamental equations of the DDE system. Here, we use the expansion function to represent the two-term semi-analytical method as follows: This expansion function (3) fulfills the given boundary conditions (2) and also has the property that the concentrations of the population densities at the impermeable boundary x = 0 are u = Σu i , v = Σv i and w = Σw i , i = 1, 2. The averaged versions of the governing PDEs, weighted using the basis functions, must be evaluated to find the free parameters that exist in (3). This technique yields the DDEs, see (A1) in the Appendix A.
The DDEs (A1) are obtained by truncating series (3) using a two-term combination. The use of the two-term method is preferred due to its accuracy without the need to lengthen the expression [14,28]. Moreover, the one-term solution (when u 2 = v 2 = w 2 = 0) was calculated to ensure accuracy in comparison.
The numerical simulation results of the DDEs (A1) is obtained using a fourth-order Runge-Kutta method, while the numerical solutions of the PDEs model (1) and (2) are found using a Crank-Nicholson finite-difference scheme. The percentage error, which is the absolute value of the difference between the approximate (one-or two-terms) and exact (numerical of PDEs) values divided by the exact value multiplied by 100, is used to calculate the difference between the one-and two-term solutions and the numerical solution of (1) and (2).

Positive Steady-State Analysis and Profiles
This section discusses the steady-state solutions of the system. All time derivative terms in (A1) are zero at the steady-state, yielding a set of six transcendental equations.
and inserted these terms back into the Equation (A1). As a result, we obtained six transcendental equations, f i = 0, i = 1, . . . , 6, which can be solved numerically and also make use of (3). For the one-term semi-analytical model case (u is = v is = w is = 0), we obtained three transcendental equations. The system yielded five non-negative steadystate solutions (ss i ). For the one-term semi-analytical model, the steady-state solutions at the origin steady-state point ss 1 = (u 1 , v 1 , w 1 ) = (0, 0, 0) always exist and the axial steady-state point ) and , 0) and one interior steadystate point ss 5 = (u * 1 , v * 1 , w * 1 ) are given by: In Figure 1a-c, the steady-state population density profiles of the prey u, the predator v, and the top predator w versus x are shown, respectively. These profiles illustrate the numerical solution of systems (1) and (2) and the solutions of the one and two terms method with the parameters α = 0.5,  (1), the two-term method yielded an excellent approximation. For the prey, predator, and top predator population density, the errors were <0.6%. The errors were marginally higher for the one-term method, but did not exceed 7%. The two-term method was superior to the one-term profile, because it more effectively modelled flat u, v, and w population density profiles.  In Figure 2, the effects of different growth rates (α) of the prey on the population density profiles of the species are presented. The figure shows the results of the two-term method at different values of α (1, 1.5, and 2); the other parameters are the same as those in Figure 1. These results suggest that an increase in the prey growth rate α is consistent with increase in the population density of the prey u and top predator w and decrease in the population density of the predator v. In other words, the top predator consumes the predator. At x = 0, the peaks of the population density of the species reach the following values: (u, v, w) = (4.517, 0.287, 3.757) when α = 1, (u, v, w) = (6.933, 0.275, 6.320) when α = 1.5, and (u, v, w) = (9.280, 0.265, 8.761) when α = 2. Furthermore, as the growth rate of prey increases, the population density of prey migrates from the center and moves to the boundary of the domain.

Theoretical Framework
This section explains the methods for determining the stability and Hopf bifurcation points of the DDEs system (A1). A common phenomenon in biological, chemical, and physical models is Hopf bifurcation, which arises when periodic solutions occur through a local change in the stability of a steady state as a result of crossing a conjugated pair of eigenvalues over an imaginary axis. The standard literature on bifurcation theory and dynamical systems [43][44][45][46] has clarified the Hopf bifurcation theory. Moreover, many researchers discussed the techniques to find the Jacobian matrix, the conditions for stability, and stability switching for prey-predator models [47][48][49][50][51]. The stability of the semi-analytical model was studied here, and the findings were used to investigate the effects of the three time delays of the species on the stability of the system (1) and (2).
The Hopf bifurcation points can be obtained by expanding a Taylor series around the steady-state solution, as follows: u jτ = u js + m j e −λt e λτ 3 , v jτ = v js + n j e −λt e λτ 1 , w jτ = w js + p j e −λt e λτ 2 ; j = 1, 2.
In addition to being linearised around the steady state, Equation (6) can be inserted into the DDE system (A1). Here, Equation (6) was inserted in the DDE system (A1) and linearised around the steady state. The Jacobian matrix eigenvalues distinguish the perturbations in the system. Therefore, a characteristic equation can be obtained for λ. In this characteristic equation, when λ = iω, the real ( e ) and imaginary (I m ) parts can be separated. At points where λ is strictly imaginary, the Hopf bifurcation points occur. Hence, the following conditions must be fulfilled to obtain the Hopf bifurcation points: f j = e = I m = 0, j = 1, . . . , 6. (6)

Hopf Bifurcation Maps and Limit Cycle
Here, semi-analytical maps including Hopf bifurcation points were created. In addition, the effects of both time delays and diffusion coefficients were investigated. Figure 4 depicts the Hopf bifurcation curve map of the prey growth rate α versus the predator mortality rate β. The results of the one-and two-term models, as well as the numerical solution are shown. The parameter values utilised were γ = 0.4, µ 1 = δ 3 = κ 3 = 0.5, µ 2 = δ 2 = κ 2 = 0.3, µ 3 = δ 1 = κ 1 = 0.1, D i = 0.05, and τ i = 1; i = 1, 2, 3. In this figure, there are two portions, as shown below. The upper portion of the Hopf bifurcation curve represents a stable region, whereas the lower portion of the curve represents an unstable region. In general, the increase in the growth rate of the prey α destabilised the system, whereas the increase in the mortality rate of the predator β stabilised it, and this is in agreement with what was mentioned in [49]. The Hopf bifurcation points can only occur for α ≥ 1.21 and 1.29 for the one-term and two-term methods, respectively. Comparison of the two-term results with the numerical solution of system (1) and (2) revealed close concordance. As such, at α = 2.3, the numerical estimate at which Hopf bifurcations can occur is β = 0.088, whereas the one-and two-term values are β = 0.102 and 0.091, respectively, with errors of 16 and 3.5%. Figure 5 further shows the Hopf bifurcation curve in the α-γ plane. The figure illustrates the one-and two-term model results with the numerical solution of the PDE system. The other parameters are the same as those in Figure 4, and β = 0.04. Once again, the Hopf bifurcation curve divides the plane into upper and lower parts. Limit cycles and unstable solutions occur at the top, whilst limit cycles disappear and stable solutions occur below the curve. The results in the figure show that the increase in the growth rate α of the prey is offset by the decrease in the mortality rate γ of the top predator. In other words, the increase in the growth rate of the prey and mortality rate of the top predator destabilise the system. There was a close concordance between the results of the two-term model and the numerical solution, with an error of <2%.   In Figure 6, the effects of hunting and maturation time delay of the three species on the stability of the system are examined. Figure 6a,b present the Hopf bifurcation curves in the α-β and α-γ planes, respectively, at different values of τ i = 1, 1.2, and 1.5, (i = 1, . . . , 3). Results of the two-term method are shown. The other variables are the same as those in Figures 4 and 5. The results show that the increase in time delay shifts the Hopf bifurcation curve upward in the α-β and downward in the in the α-γ domain. As a result of this increase, the zones of the limit cycle broadened and the stable solutions decline, leading to system destabilisation. For instance, the points (α, β) = (2, 1.5) and (α, γ) = (5, 0.1) are stable at τ i = 1 but not at τ i = 1.2 and 1.5, i = 1, 2, 3.   Figure 4. The curves separate the plane into two parts: upper and lower. Limit cycles and unstable solutions occur above the curves, whereas only stable solutions occur below them. The figure also exhibits that the increase in the delay in the maturation time of the top predator extends the zone of instability of the system. Generally, small delays result in stability whereas large delays, which suggest feedback from the distant past, result in instability [1,14]. The diffusion coefficients play an important role in determining species persistence and extinction [52]. Figure 8 exhibits the Hopf bifurcation regions for the two-term method in the D 1 -D 2 phase plane at different values of diffusion coefficient (D 3 = 0.02, 0.05, and 0.1).
The parameter values are α = 2 and β = 0.02, and the other parameters are the same as those in Figure 4. Here, the inside region shows the oscillating solutions and the limit cycles, whilst the outside region shows the stable solutions. In addition, as the value of the diffusion coefficient for the top predator (D 3 ) increases, the system dynamics vary from stable solution to limit cycle. Briefly, the alteration of the diffusion coefficients rendered the system either stable or unstable [50,51].    Figure 4. The parameters applied in this case correspond to the zone below the Hopf bifurcation curve in Figure 4, resulting in a limit cycle. The figures present the solutions of both one-and two-term methods, as well as the numerical solution of system (1) and (2). The limit cycle period in the numerical solution was estimated to be 6.84, whilst those in the one-and two-term methods were 6.86 and 6.85, respectively. Moreover, the limit cycle amplitudes of the species were (u, v, w) = (9.155, 1.517, 8. (1) and (2). The parameter values are α = 1.5 and β = 0.15, and the other parameters are the same as those in Figure 4. Since all parameters are located above the Hopf bifurcation curve in Figure 4, a stable solution is possible. When the time is sufficient (t > 60), the solution of the system converges to a coexistence steady state. This coexistence steady-state solution is given by (u, v, t) = (8.37, 0.133, 7.173) (one-term), (u, v, w) = (7.320, 0.121, 6.678) (twoterm), and (u, v, w) = (7.552, 0.122, 6.719) (numerical). These results indicate that the two-term solution is superior to the numerical solution of the PDE system, with an error of <3.1% for all species.

Bifurcation Diagrams
Utilisation of a bifurcation diagram is considered a standard way to exhibit behaviour in non-linear dynamic systems including complex periodic and chaotic behaviours [53,54]. This section compares and investigates the dynamics of the PDE system and the semianalytical DDE model. The bifurcation diagrams indicated long-term solutions for the amplitude of the steady-state branch, as well as the maximum and minimum amplitudes of the oscillatory solutions. As chaotic solutions were obtained for values of the population densities of the three species against the growth rate of the prey α, the bifurcation diagrams were drawn at the centre of the domain (x = 0) for large values of time. Figure 11 shows the bifurcation diagrams of the population densities of the prey u, predator v, and top predator w versus α. β = 0.05, and the other parameters are the same as those in Figure 4. The diagram presents the results of the two-term method for u (Figure 11a), v (Figure 11c), and w (Figure 11e), as well as the PDE numerical solution of the species u (Figure 11b), v (Figure 11d), and w (Figure 11f). In this figure, the system is initially stable at α < 1.63 and then the Hopf bifurcation occurs at α = 1.63 (twoterm) and 1.64 (numerical), providing periodic solutions in the range 1.63 ≤ α < 4.1 (two-term) and 1.64 ≤ α < 4 (numerical). At α = 4.1 (two-term) and 4 (numerical), attracting period-doubling bifurcation occurs. At α > 4.1, for each branch of the periodic solutions, the system induces a route to chaos, replacing the periodic solution cascades. The estimates of the two-term method were consistent with the numerical results of the PDE system, with an error of < 2.5%. In conclusion, the increase in the growth rate of prey in a three-species system may lead to the occurrence of chaotic solutions, as reported previously [55].

Conclusions
In this study, the dynamics of a diffusive Lotka-Volterra two predator-prey system with delays were explored. By employing the Galerkin Method, which generates semianalytical solutions, the PDE system was approximated using a DDE model. Steady-state curves and Hopf bifurcation maps were created and discussed in detail.
The effects of the growth rate of the prey and the mortality rate of the predator and top predator on system stability were demonstrated. As such, an increase in the growth rate of prey destabilised the system, whilst an increase in the mortality of both predator and top predator stabilised the system. In addition, the impact of hunting and maturation delays of the species were examined. Small delays stabilised the system, whilst great delays destabilised it. Furthermore, the effects of the diffusion coefficients of the species were investigated. The alteration of the diffusion coefficients rendered the system either permanent or extinct.
These results can serve as the foundation for future research. The present study demonstrated the high accuracy and utility of the semi-analytical model for studying the dynamics of specific biological systems, such as delayed diffusive multispecies Lotka-Volterra food chains or two-prey/one-predator systems [22]. Furthermore, in the light of harvesting in non-diffusive models [47,49], one more interesting investigation would be to vary the mortality rate as a control parameter, and examine the impacts of harvesting upon the diffusive system (1) and (2).
Funding: This research received no external funding.