Bifurcation and chaos analysis for a discrete ecological developmental systems.

This work concentrates on the dynamic analysis including bifurcation and chaos of a discrete ecological developmental systems. Specifically, it is a prey-predator-scavenger (PPS) system, which is derived by Euler discretization method. By choosing the step size h as a bifurcation parameter, we determine the set consists of all system's parameters, in which the system can undergo flip bifurcation (FB) and Neimark-Sacker bifurcation (NSB). The theoretical results are verified by some numerical simulations. It is shown that the discrete systems exhibit more interesting behaviors, including the chaotic sets, quasi-periodic orbits, and the cascade of period-doubling bifurcation in orbits of periods 2, 4, 8, 16. Finally, corresponding to the two bifurcation behaviors discussed, the maximum Lyapunov exponent is numerically calculated, which further verifies the rich dynamic characteristics of the discrete system.

Introduction

In the last decades, the study on revealing the relationship between complex dynamical behavior, such as bifurcation and chaos, and systems’ parameters has never been stopped, especially for ecological developmental system. Ecological developmental system refers to the interaction, mutual restriction, and continuous evolution between creatures and environment in a certain space of nature, so as to achieve a dynamic balance and relatively stable integrated whole. Prey–predator–scavenger systems are typical ecological developmental system. Since the foundational research by Lotka [1] and Volterra [2], prey and predator system as a typical ecological system has always aroused great interest from many researchers. By introducing different biological factors into modeling process, various prey–predator models are discussed. For example, in view of the functional response and selective capture of predator species, the dynamics of a prey–predator fishery model was studied in [3]. In addition, note that strong Allee effect and a protection zone for the prey are important factors that cannot be ignored, the authors in [4] study a reaction–diffusion predator–prey system. In [5], the authors discuss how the population dynamics is influenced by predator dormancy. Other relevant works can also be found in [6, 7].

However, it should be noted that few works have taken a scavenger species into account of the study of prey–predator systems [8, 9]. Scavengers feed on the bodies of other animals, all of which die naturally or are killed by other animals. Therefore, the study of such kinds of population relation model is significant due to the benefit to the cleaning of natural environment. Based on this, by introducing a scavenger species that scavenges the predator, the authors in [10] demonstrate that all solutions of this model are bounded and verify that there exist cascades of period-doubling orbits by numerical simulations. More recently, by considering Michaelis–Menten type of harvesting function, an ecological prey–predator–scavenger model is proposed in [11], and then, the stability and bifurcation are analyzed completely and verified numerically. In [12], by assuming the effect of harvesting and the death from toxicants, the authors propose and investigate a PPS food web model.

On the other hand, it is not difficult to understand that food and energy are of great significance for all species, such as disease spread model including SARS-CoV-2 [13]. Accordingly, it is an interesting research direction by combining the multi-species food web systems with the harvesting [14-16]. As described above, in this work, by considering the effects of harvesting and toxicant on prey, predator, and scavenger, we propose and discuss a PPS food web model. Additionally, many relevant works focus on exploring the interspecific interactions of continuous prey–predator models [17], which is mainly the existence of stable equilibrium or limit cycles. However, it is of great significance to study the discrete PPS model. One reason is that it is easier to simulate and implement by computers from the point of view of experiment or calculation. Another reason is that different from the continuous case, the discrete model can produce more abundant dynamical behaviors [18, 19].

Above all, the main contributions of this work include the following threefold.

(1) The discrete-time prey–predator–scavenger system is obtained by using the forward Euler scheme, and then, some conditions are derived which guarantee the existence and stability of the unique positive equilibrium point.

(2) Based on the flip bifurcation and Neimark–Sacker bifurcation theories of discrete systems, the sets consisting of systems’ parameters are determined, in which the discrete-time prey–predator–scavenger system can undergo such two kinds of bifurcations.

(3) By choosing some values of systems’ parameter, the bifurcation diagram, phase portrait, and the maximum Lyapunov exponents are calculated to further prove the complexity of the system dynamic behavior.

The structure of this work can be listed as follows. Sect.2 presents some conditions which guarantee the stability of the existing unique positive equilibrium. Sect.3 analyzes the flip bifurcation and Neimark–Sacker bifurcation. Then, some simulation examples are given in Sect.4to validate the theoretical results on such bifurcation behaviors. Finally, Sect. 5 draws conclusion.

Existence and stability of positive equilibrium

If x(t), y(t), and z(t) are used to denote the density of prey, predator, and scavenger at time t, respectively, then the PPS model is established aswith , and . The definitions of system parameters are given in Table 1, where all parameters are positive. The forward Euler discretization method is applied to system (1), and we get the discrete PPS system aswhere h is the step length in the process of Euler discretization, and x(n), y(n), and z(n) are approximate values corresponding to x(nh), y(nh), and z(nh) for sufficiently small h. Then, selecting h as the bifurcation parameter, its effect on the complex dynamical behaviors of system (2) will be discussed, which can show that (2) can undergo FB and NSB with the variation of parameters.

Table 1

Parameters and their meanings

Parameter	Description
r	The intrinsic growth rate of the prey species
K	The carrying capacity without predation, harvesting, and toxicant
β	The combined harvesting effort
α_i(i=1,2,3)	The predation and scavenge with positive maximum attack rate
c_i(i=1,2,3)	The positive catch ability coefficients
γ_i(i=1,2,3)	The coefficients of toxicity of prey, predator, and scavenger
d_j(j=1,2)	The predator and scavenger decay with natural death rates
e_j(j=1,2)	The conversion rates of prey to predator and scavenger

Note that if the step size h in the discretization process is sufficiently small, then the equilibrium point of (2) satisfiesSome direct calculations give that system (2) has four equilibrium points, involving and , whereIn addition, the predator-free equilibrium point can be denoted by , whereand .

From a biological perspective, it makes sense to discuss the positive equilibrium of (2). Therefore, the following theorem is given.

Theorem 1

Denote as the positive equilibrium of system (2), ifandhold, and then, is an unique positive equilibrium of (2). Furthermore,whereand is a positive solution of the equation,

Next, the stability of the positive equilibrium point can be investigated. The calculation of Jacobian matrix J of the system (2) corresponding to yieldsThen, the characteristic equation of at can be written aswherewithConsequently, we can obtain Theorem 2.

Theorem 2

Assume that the condition holds, then of system (2) is locally asymptotically stable.

Proof

It is known from the Routh–Hurwitz criterion that equation (6) has three roots with negative real parts provided that and . Then, a direct calculation yields the conclusion of Theorem 2 and thus completes the proof.

Bifurcation behavior analysis

By choosing h as bifurcation parameter, we will discuss FB and NSB. Denote and as the critical values, respectively, when system undergoes FB and NSB. Based on the FB and NSB theory, the following results can be derived.

Theorem 3

System (2) can undergo flip bifurcation when one of the following cases hold, where one is that ifholds, thenand the other is that ifholds, then

Note that if characteristic equation (5) has a solution , then system (2) may undergo flip bifurcation. Therefore, let and substitute it into (6) yielding Moreover, other eigenvalues should satisfy , because and are two roots of the equation Then, by a direct calculation, we can derive the conclusion and Theorem 3 is thus completed.

a Local amplification corresponding to Fig. 1 for ; b the largest Lyapunov exponents associated with a

Fig. 1

(a) Bifurcation diagram of system (2) in (h, x) plane with the initial value (0.4, 0.2, 0.3); b the largest Lyapunov exponents associated with a

Theorem 4

Suppose that andand . Then, system (2) may undergo NSB.

Fig. 6

Phase portrait for different values of bifurcation parameter h

On the basis of the Neimark–Sacker bifurcation theory, it is known that Neimark–Sacker bifurcation occurs when the characteristic equation (5) have two conjugate complex root with module one. Transform (5) into the following formIn addition, it also should satisfy and . By comparing (7) with (5), we can obtain Thus, the critical value of bifurcation parameter should satisfyThat is, the critical value of NSB is the minimum positive real root. It follows from (7) that the three eigenvalues areThen, based on the nondegenerate condition and transversality condition, we can arrive at the conclusion and thus Theorem 4 is completed.

LetandThe positive equilibrium point may undergo FB when parameters are chosen from the sets of or . LetThe positive equilibrium point may undergo NSB when parameters are chosen from .

Illustrative example

This section will verify the theoretical results by using MATLAB simulation analysis tool. Bifurcation diagrams, phase diagram, as well as the largest Lyapunov exponents will be depicted by choosing some parameters from the sets or , and . It can be seen from Theorem 1 that system (2) has only one positive equilibrium point. By computing the equilibrium of the system (2), we have that the FB emerges from (3.4155, 0.9859, 1.2455) at andThus, the correctness of Theorem 4 is further verified.

Varying h in range and fixing ;

Varying h in range and fixing .

It can be seen from Fig. 1a that is stable for and when , it loses the stability, and then, system (2) undergoes flip bifurcation. The largest Lyapunov exponents associated with Fig. 1a are obtained in Fig. 1b. Local amplification corresponding to Fig. 1a is presented when , as well as the corresponding largest Lyapunov exponents. The phase portraits are depicted in Fig. 5.

Fig. 5

a Local amplification corresponding to Fig. 6 for ; b the largest Lyapunov exponents associated with (a)

Based on Theorem 1, it is clear that system (2) has only one positive equilibrium. Through analysis and calculation, NSB generates at the equilibrium point (1.5867, 0.9358, 1.077) when andThus, we verify Theorem 4.

It can be seen from Fig. 4a that the equilibrium point is stable for and loses its stability , in which system (2) undergoes NSB. Moreover, an invariant circle appears when the parameter exceeds 0.338. The largest Lyapunov exponents associated with Fig. 4a are calculated in Fig. 4b. Local amplification corresponding to Fig. 4a is presented when , as well as the corresponding largest Lyapunov exponents. The phase portraits corresponding to Fig. 4a are depicted in Fig. 6.

Fig. 4

a Bifurcation diagram of system (2) in (h, x) plane with the initial value (1.5, 0.2, 0.4); b the largest Lyapunov exponents

Conclusions

This paper has analyzed the dynamics of the discrete PPS model which are discussed. We firstly prove the existence and stability of an unique positive equilibrium point by analyzing the characteristics equation. Based on the bifurcation theory, the sets consisting of system parameters are derived, in which the systems can undergo FB and NSB. Moreover, compared with the continuous case, it shows that the discrete PPS system exhibits many interesting dynamical behaviors, such as period-doubling bifurcation, invariant cycle, quasi-periodic orbit, and chaotic set.

The future works may cover the design and applications of control strategy for ecological system with time delays [20], as well as circuit realization of bifurcations for delay systems [21, 22].