Abstract

In this paper, we put forward a time-delay ecological competition system with food restriction and diffusion terms under Neumann boundary conditions. For the case without delay, the conditions for local asymptotic stability and Turing instability are constructed. For the case with delay, the existence of Hopf bifurcation is demonstrated by analyzing the root distribution of the corresponding characteristic equations. Furthermore, by using the normal form theory and the center manifold reduction of partial functional differential equations, explicit formulas are obtained to determine the direction of bifurcations and the stability of bifurcating periodic solutions. Finally, some simulation examples are provided to substantiate our analysis.

1. Introduction

Based on the diversity of biological populations and ecosystems, various ecological competition systems have been established and widely studied (see [15]). In ecological competition systems, the interaction between populations is usually reflected by functional response functions. Results have shown that the predation relationship between populations greatly affects the dynamic behavior of predator-prey systems (see [68]). In fact, with the development of the economy, humans will harvest biological populations and develop related biological resources to obtain economic benefits, such as fisheries, forestry, and wildlife management (see [911]). In recent years, many scholars have introduced the harvesting terms into biological systems to study the system modeling and related dynamic characteristics [1215]. May et al. [16] established the following model to analyse the interaction between populations under different harvesting strategies:where and represent the prey and predator densities, respectively. and indicate the intrinsic growth rates. represents the internal growth limit of the population without predators. stands for the maximum predation coefficient of predator. is the quality standard for measuring prey as food. and represent the impact of human harvesting for predator and prey populations. On this basis, scholars have studied the relevant dynamic characteristics of different harvesting terms in system (1). In [17], the authors studied the case of and in system (1), that is, the constant effort harvesting to both prey and predator population. Particularly, they analyzed the maximum sustainable yield of another population under the condition of limiting the harvest of one population. In [18], the authors discussed the case of and for system (1), namely constant yield harvesting for prey and constant effort harvesting for the predator population. At the same time, the relevant dynamic characteristics of the system were analyzed and the effects of different management strategies on the stability of the system were compared. Moreover, in [19, 20], the authors considered the constant yield harvesting of the prey population in system (1), that is, and . They obtained the conditions of Bogdanov-Takens bifurcation and supercritical/subcritical Hopf bifurcation of codimension 1.

Based on the research of the appeal literature, it is found that only the impact of linear capture is taken into account. However, from the perspective of biology and economics, the linear harvesting term cannot accurately reflect the real social activities (harvesting population and developing biological resources) of humans and their impact on the predator-prey system. Thus, the nonlinear harvesting terms have been introduced to model the dynamics of predator-prey systems in recent years (see [2123]). In [24], a realistic harvesting functional form was proposed as follows:where represents the catchability coefficient and stands for the harvesting effort. The conditions for Hopf bifurcation were derived and a nonlinear state feedback controller was designed to control the Hopf bifurcation.

Consider the logistic equation

Here, is the intrinsic growth rate and represents the carrying capacity. Meanwhile, we can see from (3) that the predicted relation of the specific growth rate, , to mass, , is a straight line. However, when considering a measure of the portion of available limiting factors not yet utilized by the population, the linear average growth rate cannot accurately describe the growth trend of the population. Research shows the growth and development of organisms depend on the availability and utilization of food in the living environment, which implies that different degrees of food supply rates will affect the stable age composition of the population and then have an influence on the average growth rate of the population. Thus, when growth limitations are based on the proportion of available resources not utilized, Smith [25] established a food-limited growth functionwhere is positive constant and represents the rate of replacement of mass in the population at saturation. It follows from (4) that the predicted relation of the specific growth rate, , to mass, , is not a straight line but a concave curve.

Considering the growth limit is based on the proportion of unused available resources and setting , , then system (1) becomeswhere represents the functional response function between populations. represents that when the density of prey population is low, the predator population can switch to other prey for capture. Then, we will analyse the cases on and in the discussion that follows, respectively.

In nature, the survival and development of a population often depend on the amount of food and the living space available. Importantly, the greater the population density, the higher the requirements for the living environment. At the same time, the acquisition of food mainly depends on the living environment of the population, which shows that the change of the population’s living space affects the survival and development of the population to a great extent. Therefore, the population will instinctively migrate and diffuse in space to seek a more suitable environment for survival and development. It is necessary to consider the influence of the diffusion effect on population dynamics in predator-prey systems. Mathematically, the nonlinear system with diffusion will show complex dynamic properties [2629]. In the reaction-diffusion system proposed by Turing [30, 31], the spatial heterogeneity caused by the internal reaction-diffusion characteristics of the system results in the loss of system symmetry and makes the system self-organize to produce some spatial patterns. The process of pattern formation is called Turing instability (Turing bifurcation). The symmetry of the system is broken, leading to the formation of Turing patterns. Therefore, we call this phenomenon “Turing instability caused by diffusive reaction” [32].

On the other hand, time delay has become a factor that cannot be ignored in many biological dynamic systems. A large number of studies have revealed that time delay has an important impact on the dynamic characteristics of biological systems and it is common in predator-prey systems, mainly including mature time delay, capture time delay, and pregnancy time delay [3335]. Local stability of the system means that if the initial state is adjacent to the equilibrium state, the system will not vibrate, and its state trajectory will eventually fall to the equilibrium state. In particular, Hopf bifurcation is a dynamic bifurcation phenomenon, which shows that when the parameters change near the critical value, the stability of the equilibrium point will change and periodic solutions will be generated in its small neighborhood. Meanwhile, it is found that time delay as a bifurcation parameter can induce Hopf bifurcation [29, 36, 37]. Therefore, in this paper, we consider the pregnancy delay of the predator population and analyse the dynamic characteristics of ecological competition systems.

Based on the discussion above, we introduce the diffusion effect of the population [2629] and the pregnancy delay of the predator population [3335] into system (5) to explore its impacts on the dynamic characteristics of the ecological competition system, which can be described bywith Neumann boundary conditions and initial conditions

Here, stand for the population density of prey and predator at a spatial location and time , respectively. represent the diffusion coefficients associated to and , respectively. denotes the Laplacian operator in . Suppose is a bounded domain with a smooth boundary .

Lemma 1. For any solution of system (6) without delay,

Proof. Let be a nonnegative solution of system (6) without delay. Note that the functional response . Then,We can estimate the upper limits of and due to the standard comparison principle:In other words, for arbitrary , there exists positive constants such thatfor , , andfor , . So, the conclusion follows immediately.
Compared with the models proposed in [15], model (6) is a new measure of the population with a nonlinear average growth rate based on food restriction. At the same time, we consider the reaction-diffusion factors for system (6) to study their influence on the dynamic behaviors of systems. In addition, on the basis of references [1620], we introduce the nonlinear harvesting term of the harvesting-effort coefficient into system (6) and study the stability and related dynamic characteristics of the system under Holling I and Holling II functional response functions. Importantly, the addition of pregnancy delay can more accurately reflect the evolution of the population and make the system show more complex dynamic characteristics than the model without delay, which is also a widely concerned direction in the research of biological systems.
The main contributions of this paper can be stated as follows. In Section 2, the effects of diffusion on the dynamic behavior of the systems without time delay are investigated and some conditions for system stability and Turing instability are determined. It is found that the appropriate diffusion coefficients will lead to Turing instability. In Section 3, we analyse the stability of equilibrium and Hopf bifurcation in the predator-prey system with time delay as the bifurcation parameter. The condition for Hopf bifurcation is constructed and the expression of the bifurcation threshold is given. In Section 4, we calculated the direction of the bifurcations to get more information about the bifurcations. Section 5 uses some numerical examples to verify the correctness of the previous derivation. Section 6 gives the conclusion of this paper.

2. Equilibrium Stability and Turing Instability Analysis

Assume that the predator-prey relationship in system (6) satisfies Holling type I functional response, that is , then we make the following nondimensional scaling transformation [22]:

Dropping the tildes, system (6) can be rewritten by

Consider the case of no time delay in system (9), namely, , then system (9) becomes

Considering the practical significance of the ecosystem, we are interested in the coexisting equilibrium. In order to obtain the positive equilibrium of system (10), let

Then system (10) has a positive equilibrium , where and satisfies the following quadratic equation:where

Lemma 2. For equation (12), we come to the following results:(1)If , then equation (12) has no positive roots.(2)If , then equation (12) has a unique positive root .

Hence, when , system (10) has a unique positive equilibrium , where .

Consider system (10) without diffusion

The Jacobian matrix of (19) at is , where

Lemma 3. When , then we can deduce .

Proof. Notice thatThus, we can obtain when . As , then we can deduce .
Then from the Jacobian matrix of (13), we get the characteristic equationwhere

Theorem 1. If , then of system (13) is locally asymptotically stable.

Proof. By Lemma 3, we know when . Together with and , we can get and . Thus, by the Routh-Hurwitz criterion, of system (13) is locally asymptotically stable.
It follows from [38] that the Laplacian operator has the eigenvalue under the homogeneous Neumann boundary condition. And the corresponding eigenfunctions are , where and construct a basis of the phase space and is defined bywith the inner product . Thus, for system (10), the characteristic equation at iswhereHence,

Theorem 2. If , then of system (10) is asymptotically stable.

Proof. It follows from Lemma 3 that when . Notice that and . For equation (14), we obtainThus, all roots of (25) have negative real parts for , which implies that of system (10) is asymptotically stable.

Remark 1. Assume that there are no time delays and for system (6). It can be seen from Theorems 1 and 2 that when , the introduction of diffusion terms does not change the stability of , which means that Turing instability does not occur.
Next, we consider the predator-prey relationship among populations in system (6) satisfying Holling type II functional response, that is , where represents the maximum per capita reduction rate of the prey population and stands for the average saturation rate [3]. Then we make the following nondimensional scaling transformation [22]: Dropping the tildes, then system (6) without time delay turns intoSystem (15) has a positive equilibrium , where and satisfies the following quadratic equation:where

Lemma 4. For equation (16), we have the following results:(1)If and , then equation (16) has positive roots(2)If and , then equation (16) has a unique positive root .(3)If , then equation (16) has a unique positive root .

Theorem 3. For system (15), we come to the following results:(1)If and , then system (15) has two positive equilibrium points and , where .(2)If and , then system (15) has a unique positive equilibrium point , where .(3)If , then system (15) has a unique positive equilibrium point , where .

The linearization of system (15) at can be represented bywhere , and

Thus, the characteristic equation iswhere

Obviously, for (36), we have

We make the following assumptions:

Theorem 4. If hold and , then system (15) is locally asymptotically stable at .

Proof. If and hold, for equation (17), we haveThus, all roots of (36) have negative real parts. By the Routh-Hurwitz criterion, system (15) is locally asymptotically stable at .
Then we can describe system (15) without diffusion term by the following equations:The characteristic equation of system (18) at iswhere

Lemma 5. If , and are satisfied, then of system (18) is locally asymptotically stable.

Let

We make the following assumptions:

Lemma 6. If , and hold, then has two positive roots .

Theorem 5. Assume that and , and hold. Then, diffusion-driven instability occurs for system (15) at if there exists a such that for .

Proof. It follows from Lemmas 5 and 6 that system (18) is stable at and has two positive roots when and , and hold. Then for the spatial system (15), the corresponding characteristic equation (17) has an eigenvalue with positive real part if there exists a such that for . This implies the spatial system (15) is unstable at , that is, the diffusion-driven instability occurs.

Remark 2. Assume that there are no time delays and for system (6). It can be seen from Theorem 5 that the introduction of diffusion terms may change the stability of , resulting in Turing instability.

3. Hopf Bifurcation Analysis

In this section, we consider the effect of time delay on the dynamics of the system and get the conditions for Hopf bifurcation of system (9).

Linearizing system (9) at [3436, 39], we havewhere , , . Thus, the characteristic equation iswhere stands for identity matrix and . (47) is equivalent towhere

Remark 3. It should be pointed out that this paper adopts the linearization method [29, 36, 37, 39,40] to deal with the dynamics analysis of system (9), including the local stability, Turing instability, and Hopf bifurcation. It is common knowledge that Lyapunov’s second method is important to the stability theory of dynamical systems and control theory. However, this method is not suitable for investigating the dynamics of the ecological competitive system with delay and diffusion proposed in this paper. The Lyapunov stability criterion can only give a sufficient condition for the stability of a system. In this paper, not only the condition of the local stability is established, but also the boundary of stability (the onset of Hopf bifurcation) is determined.

Lemma 7. If , then all roots of equation (21) have negative real parts when .

When , assume is a pair of pure imaginary roots of (48), satisfies the following equation:

By separating real and imaginary parts, we obtain

Add the squares of both sides of equation (23) to getwhere

Let , (52) can be converted to

Then we make the following assumptions:

Theorem 6. If and hold, then system (9) is locally asymptotically stable at for all .

Proof. By Theorem 2, we know when . If is satisfied, we obtainThus, (54) has no positive roots, which implies (48) has no purely imaginary roots. By Lemma 7, we deduce that all roots of (48) have negative real parts for all , that is, system (9) is locally asymptotically stable at for all .
LetIt follows that when and hold, has a unique positive root having the following form:Denote

Lemma 8. If and are satisfied, equation (21) has a pair of purely imaginary roots at , and

Theorem 7. If and hold, we have

Proof. From equation (28), we obtainBy Lemma 3, we know that when . Together with and , we can get that , is increasing and is decreasing in when hold. Thus,Notice that and is increasing in with . We can deduce that is decreasing in . Hence, is increasing in .
As , we obtainLet be the root of (48) satisfying and .

Theorem 8. For , we have

Proof. By derivation of equation (21) with respect to the delay , we havethenIt follows from (50) and (51) thatNoticing , we obtainTherefore, the transversality condition holds.

Theorem 9. Suppose that and hold. The following statements are valid:(1)System (9) is locally asymptotically stable at when and unstable for .(2)System (9) undergoes a Hopf bifurcations at when .

4. Direction and Stability of Hopf Bifurcation

The previous analysis has shown that system (9) admits a series of periodic solutions bifurcating from the trivial uniform steady state at some critical values. Then, in this section, we are concerned with the direction of Hopf bifurcations and the stability of bifurcating periodic solution.

Denote by and introduce the new parameter . Normalizing the delay by the timing-scaling , then system (9) can be transformed in the phase space where and are given, respectively, bywith , andfor .

From the previous discussions, we know that are a pair of simple purely imaginary eigenvalues of the linear equation

By the Riesz representation theorem, there exists a matrix function such that

Actually, we can take

Let denote the infinitesimal generator of the semigroup and be the formal adjoint of under the bilinear pairingfor . Then, we can give the definition of is

From the above analysis, we get that has a pair of simple purely imaginary eigenvalues . Supposing that , then we can calculate the eigenfunction of corresponding to and the eigenfunction of corresponding to . From (78), we derive

Then, we obtain

Let , wherefor , andfor . Let , then by (77), we can calculatewhere

Next, we define and construct a new basis

Then, we obtain . Denote and , where . Thus, the center space of the (74) is given by having the following form:wherefor .

Let be the infinitesimal generator induced by the solution of (74), then system (9) can be rewritten bywhere

From the decomposition and (38), where denotes the complement subspace of in , the solutions of system (39) can be written aswhere

The bifurcation direction of this part is based on , so the solution of system (31) on the center manifold is given by

Let and notice , then (93) becomeswhere

By [40], it is easy to know satisfieswith

By Taylor expansion of and , we get

From (94) and (98), we obtainwhere

Thus, we can getwithwhere

Denoteand notice that

Then we deduce

Thus, by equations (43), (45), and (46), we have for . If , we get

Also, for , we can get

So far, we have given the expressions of .

It follows from (98) that

By [40], we know satisfieswhere

Substituting (96) into the derivative of (98) and comparing the coefficients with (114) and (115), we have

From (115), we have that for ,

Hence, for ,where

Thus, we can get from (78) and (116) thatfor . That is,where and are constant vectors. From the above discussion and (116), we have

Then, we obtain

Therefore, each can be determined and we can evaluate the following values determining the direction and stability of bifurcating periodic orbits:

Theorem 10. For any critical value , we have(1)When , the Hopf bifurcation is supercritical (subcritical), that is, the bifurcating periodic solutions exist for .(2)When , the bifurcating periodic solutions on the center manifold are orbitally asymptotically unstable (stable).(3)When , the period of bifurcating periodic solutions increases (decreases).

5. Numerical Simulations

In the following, we will carry out numerical simulations of three impacts of diffusion, time delay, and harvesting effort to illustrate our theoretical findings in the previous sections.

5.1. The Impact of Diffusion

We confirm that system (13) has a unique equilibrium with the parameters . By a simple verification, we can see that is satisfied. By Theorem 1, we deduce that system (13) is locally asymptotically stable at . Meanwhile, it is revealed from Theorem 2 that system (10) is still stable at , which implies Turing instability will not occur in system (10) (see Figures 13).

Next, we consider system (15) and set . When , system (15) can be transformed into the nondiffusion system (18). By calculation, we get , and , and are satisfied. Thus, from Lemma 5, system (18) is asymptotically stable at . Then we set and other variables are the same as above. Clearly, and are reached. It follows from Theorem 4 that system (15) is still stable at (see Figure 4). When and other variables are the same as above. By a simple calculation, we can see that , , , , and are satisfied. By Theorem 5, the diffusion-driven instability occurs for system (15) at , which is shown in Figures 5 and 6.

5.2. The Impact of Time Delay

For system (9), we take . By calculation, we get . In the meantime, and are satisfied. By Theorem 6, system (9) is asymptotically stable at for all , as displayed in Figures 7 and 8.

Next, we set for system (9), which has a unique positive equilibrium . It is verified that and are satisfied. Meanwhile, we get , , implying that the cross-section condition holds. By Theorem 9, system (9) is locally asymptotically stable at when (see Figures 9 and 10) and the Hopf bifurcation occurs when crosses the critical value (see Figures 11 and 12).

Furthermore, from formula (125), we obtain . It follows from Theorem 10 that the direction of the bifurcation in system (9) at is forward, the bifurcating period solutions are stable, and the period of bifurcating periodic solutions is increasing.

5.3. The Impact of Harvesting Effort

Considering the impact of harvesting efforts, we fix and let vary in . The stability and instability regions for system (9) are depicted by mapping the nonlinear harvesting to the critical value in Figure 13. Obviously, the critical value of bifurcation increases with the increase of , that is, the stable region of system (9) is expanded with the increase of .

6. Conclusion

We investigate the spatiotemporal dynamic evolution of a time-delay ecological competition system with food restriction and diffusion terms under Neumann boundary conditions. The conditions of asymptotic stability and Turing instability at the positive equilibrium of delay-free systems are obtained under different functional response functions. Compared with the classical population growth model described by the logistic equation, the model with food restriction studied in this paper is more in line with actual biological competition systems. The results show that the diffusion phenomenon caused by the change of population position in space will seriously affect the stability of biological competition systems, eventually resulting in the appearance of Turing instability. Then, by selecting the time delay as the bifurcation parameter, we reveal that the delay can cause very complex dynamic phenomena. When the delay is less than the bifurcation critical value, the system maintains asymptotic stability at the positive equilibrium point, while the system becomes unstable and produces a Hopf bifurcation when the delay is greater than the bifurcation critical value. Meanwhile, by using the central manifold method, we derive the conditions for determining the bifurcation direction and the stability of the bifurcation periodic solution. The harvesting effort has a major influence on the stability and Hopf bifurcation. As the harvesting effort increases in the appropriate range, the bifurcation critical value increases; that is, the stable region of the system is expanded. Therefore, we can conclude that the appropriate harvesting of biological populations and the rational development of biological resources can not only meet human economic needs but also play a beneficial role in biological systems.

Data Availability

The data in relation to the findings of this study are available upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Nos. 62073172 and 61573194), the Open Research Project of the State Key Laboratory of Industrial Control Technology of Zhejiang University (No. ICT2022B43), and the Postgraduate Research and Practice Innovation Program of Jiangsu Province (No. KYCX21_0754)