Bifurcation Analysis of a Singular Bioeconomic Model with Allee Effect and Two Time Delays

Abstract

A singular prey-predator model with time delays is formulated and analyzed. Allee effect is considered on the growth of the prey population. The singular prey-predator model is transformed into its normal form by using differential-algebraic system theory. We study its dynamics in terms of local analysis and Hopf bifurcation. The existence of periodic solutions via Hopf bifurcation with respect to two delays is established. In particular, we study the direction of Hopf bifurcation and the stability of bifurcated periodic solutions by applying the normal form theory and the center manifold argument. Finally, numerical simulations are included supporting the theoretical analysis and displaying the complex dynamical behavior of the model outside the domain of stability.

1. Introduction

The economic theory proposed by Gordon [1] in 1954 was described as follows: net economic revenue (NER) = total revenue (TR) − total cost (TC).

This provides theoretical evidence for the establishment of singular bioeconomic model, which is described by differential-algebraic equations. Recently, there has been extensive literature dealing with such systems, describing the interactions between the different species and harvesting effort regarding activity, stabilities of equilibrium, bifurcations, and other dynamics (see, e.g., [2–5] and the references therein).

In most of ecosystems, since one species does not respond instantaneously to interactions with other species, some delays due to several reasons, such as gestation, hunting, and maturation, are required. To incorporate this idea in a modelling approach, time delays have been introduced into ecosystems. For a long time, it has been recognized that delays can have very complicated impact on the dynamics of a system (see, e.g., monographs by Hale [6], Hale and Verduyn Lunel [7], Kuang [8], Yuan and Song [9], and Wu [10]). Generally speaking, delays can cause the loss of stability and lead to periodic solutions. Some authors like in [11] described the effects of time delay in a prey-predator model incorporating parasite infection for the prey population. The literature [12] described the dynamics of a stage structured population model with time delay in fluctuating environment. The literature [13] investigated dynamics of delayed prey-predator model with harvesting. Chakraborty et al. [14] introduced a single discrete gestation delay in a differential-algebraic biological economic system and showed Hopf bifurcation in the neighborhood of coexisting equilibrium point through considering the delay as a bifurcation parameter. Zhang et al. [15] studied a ratio-dependent prey-predator singular model and analyzed the direction and stability of periodic solutions. Some other authors have discussed dynamical models with multiple delays. While analyzing, the multiply delayed models are mostly simplified by taking equal magnitude for two delays 𝜏₁=𝜏₂=𝜏 [16, 17] or choosing the sum 𝜏 of two delays as a bifurcation parameter [18, 19]. However, the delays appearing in different terms of an ecological system are not always equal. Therefore, it is necessary to discuss dynamical system with different delays.

In addition, numerous examples demonstrate that the growth of natural populations can exhibit Allee effect [20], which describes a positive relation between population density and the per capita growth rate. When populations get larger, the effect will saturate or disappear. Allee effect may result from some causes, such as mate finding, social dysfunction, inbreeding depression, fool exploitation, and predator avoidance of defense. Therefore, in this paper, we consider a differential-algebraic prey-predator model with two time delays and the Allee effect as follows: 𝑑𝑥𝑥𝑑𝑡=𝑥𝐴+𝑥−𝑑₁−𝑥𝑡−𝜏₁−𝑥𝑦1+ℎ₁𝑥,𝑑𝑦𝑑𝑡=𝑦𝑏𝑥𝑡−𝜏₂1+ℎ₁𝑥𝑡−𝜏₂−𝑑₂−𝐸𝑦,0=𝐸𝑝𝑦−𝑐−𝑚, (1) where 𝐴 is the Allee effect constant of the prey species, respectively. 𝑑₁ and 𝑑₂ are the intrinsic mortality rate of the prey and predator population, respectively. ℎ₁ is the handling time. 𝑏 denotes food utilization efficiency. Prey dynamics is delayed by 𝜏₁ due to crowing and the predator takes time 𝜏₂ to convert the food into its growth; 𝐸 is harvesting effort for predator, and 𝑝>0, 𝑐>0, and 𝑚>0 are harvesting reward per unit harvesting effort for unit weight of predator, harvesting cost per unit harvesting effort for predator, and the net economic revenue per unit harvesting effort, respectively. All the parameters are positive constants.

In this paper, we investigate the effects of two delays on the dynamics of singular bioeconomic model with Allee effect. The existence of periodic solution has been explored through Hopf bifurcation. The domain of stability is defined. Then, the formula for determining the properties of Hopf bifurcation is derived by using the normal form method and center manifold theorem. Finally, numerical simulations show the bifurcation plot with respect to time delay and give the effectiveness of the result mentioned above.

2. Stability and Existence of Hopf Bifurcation

There exists a positive equilibrium 𝑃₀=(𝑥₀,𝑦₀,𝐸₀) for model (1), where coordinates 𝑥₀, 𝑦₀, and 𝐸₀ satisfy the following equations: 𝑥𝐴+𝑥−𝑑₁𝑦−𝑥−1+ℎ₁𝑥=0,𝑏𝑥𝑡−𝜏₂1+ℎ₁𝑥𝑡−𝜏₂−𝑑₂𝐸−𝐸=0,𝑝𝑦−𝑐−𝑚=0. (2)

For the sake of the simplicity, denote the two differential equations by 𝑓₁(𝑥,𝑦,𝐸) and 𝑓₂(𝑥,𝑦,𝐸) and the algebraic equation by 𝑔(𝑥,𝑦,𝐸).

Let 𝑢=𝑥, 𝑣=𝑦, and 𝐸=(𝑝𝐸₀𝑦/(𝑝𝑦₀−𝑐))+𝐸. Therefore, the system (1) can be rewritten in the following form: 𝑑𝑢𝑢𝑑𝑡=𝑢𝐴+𝑢−𝑑₁−𝑢𝑡−𝜏₁−𝑢𝑣1+ℎ₁𝑢,𝑑𝑣𝑑𝑡=𝑣𝑏𝑢𝑡−𝜏₂1+ℎ₁𝑢𝑡−𝜏₂−𝑑₂−𝐸+𝑝𝐸₀𝑣𝑝𝑦₀,−𝑐0=𝐸−𝑝𝐸₀𝑣𝑝𝑦₀−𝑐𝑝𝑣−𝑐−𝑚. (3) Consider the following local parametric 𝜓 of the third equation of the system (3): 𝑢,𝑣,𝐸^𝑇=𝜓𝜏₁,𝜏₂,𝑌=𝑁^𝑇₀+𝑈₀𝑌+𝑉₀ℎ𝑌(𝑡),𝑔𝜓(𝑌(𝑡))=0, (4) where 𝑈₀=¹⁰⁰¹₀₀, 𝑉₀=⁰⁰₁, 𝑌=(𝑦₁,𝑦₂)^𝑇, 𝑁₀=(𝑢₀,𝑣₀,𝐸₀), and ℎ(𝑦₁,𝑦₂)∶𝑅²→𝑅 is a smooth mapping and has the following expression: ℎ𝑦₁,𝑦₂=𝑝𝐸₀𝑣₀+𝑦₂𝑝𝑣₀+𝑚−𝑐𝑝𝑣₀+𝑦₂−−𝑐𝑝𝐸₀𝑦₀𝑝𝑣₀−𝑐−𝐸₀. (5) Then, the linearization part of the system (3) is 𝐽𝑁₀=𝐷_𝑁𝑓₁𝑁₀𝐷_𝑁𝑓₂𝑁₀𝐷_𝑁𝑔(𝑁₀)𝑈^𝑇₀⁻¹=001001𝐴𝑢₀𝐴+𝑢₀²+ℎ₁𝑢₀𝑣₀1+ℎ₁𝑢₀²−𝑢₀𝑒^−𝜆𝜏1−𝑢₀1+ℎ₁𝑢₀𝑏𝑣₀1+ℎ₁𝑢₀²𝑒^−𝜆𝜏2𝑝𝐸₀𝑣₀𝑝𝑣₀.−𝑐 (6) For the sake of the simplicity, let 𝑢₁=𝑢₁₁+𝑢₁₂,𝑢₁₁=−𝑢₀,𝑢₁₂=𝐴𝑢₀𝐴+𝑢₀²+ℎ𝑢₀𝑣₀1+ℎ𝑢₀²,𝑢₂𝑢=−⁰1+ℎ𝑢₀,𝑢₃=𝑏𝑣₀1+ℎ𝑢₀²,𝑢₄=𝑝𝑣₀𝐸₀𝑝𝑣₀.−𝑚 (7)

According to the Jacobian matrix 𝐽, the characteristic equation of (1) at 𝑃₀ can be written as follows: 𝑅₁(𝜆)=𝜆²−𝑢₄+𝑢₁₂𝜆+𝑢₁₂𝑢₄−𝑢₁₁𝜆−𝑢₄𝑒^−𝜆𝜏1−𝑢₂𝑢₃𝑒^−𝜆𝜏2=0. (8)

When the time delays are not considered, the characteristic polynomial for model (1) can be expressed as the following form: 𝜆²−𝑢₁+𝑢₄𝜆+𝑢₁𝑢₄−𝑢₂𝑢₃=0. (9) It is clear that if the following conditions are satisfied: 𝑢₁+𝑢₄<0,𝑢₁𝑢₄−𝑢₂𝑢₃>0, (10) then 𝑃₀ is locally asymptotically stable for model (1).

Next, we consider the effects of the two time delays on the stability of model (1) in the following cases.

Case 1 (𝜏₁=0 and 𝜏₂>0). Due to 𝜏₁=0, (8) becomes 𝜆²−𝑢₁+𝑢₄𝜆+𝑢₁𝑢₄−𝑢₂𝑢₃=0. (11)

Assume that a purely imaginary solution of the form 𝜆=𝑖𝜔₂ exists in (11). Substituting it into (11) and separating the real and imaginary parts, we have −𝜔²₂+𝑢₁₁+𝑢₁₂𝑢₄=𝑢₂𝑢₃cos𝜔₂𝜏₂,−𝑢₁₁+𝑢₁₂+𝑢₄𝜔₂=−𝑢₂𝑢₃sin𝜔₂𝜏₂. (12) Taking square on both sides of (12) and summing them up, we obtain 𝜔⁴₂+𝑢²₁+𝑢²₄𝜔²₂+𝑢²₁𝑢²₄−𝑢²₂𝑢²₃=0. (13) From condition (10), (13) has unique positive real root 𝜔₂₀ if 𝑢₁𝑢₄+𝑢₂𝑢₃<0. (14) The critical value of the delay corresponding to 𝜔₂₀ is given by 𝜏_2𝑘=1𝜔₂₀arccos−𝜔²₂₀+𝑢₁𝑢₄𝑢₂𝑢₃+2𝑘𝜋𝜔₂₀,𝑘=0,1,2,…. (15) If 𝑢₁𝑢₄+𝑢₂𝑢₃>0, (13) has no real root and model (1) is asymptotically stable for any time delay 𝜏₂>0.

Now, differentiating (8) with respect to 𝜏₂ and substituting the eigenvalue 𝑖𝜔₂₀, it follows that 𝑑sign𝑑𝜏₂Re(𝜆)_𝜏2𝑘=sign𝜔⁴−𝑢²₁𝑢²₄+𝑢²₂𝑢²₃>0. (16)

Theorem 1. Assume that condition (14) holds; then, there exists a 𝜏₂₀>0 given by 𝜏₂₀=1𝜔₂₀arccos−𝜔²₂₀+𝑢₁𝑢₄𝑢₂𝑢₃+2𝑘𝜋𝜔₂₀, (17) such that the positive equilibrium 𝑃₀ of model (1) is locally asymptotically stable for 𝜏₂<𝜏₂₀ and undergoes a Hopf bifurcation when 𝜏₂=𝜏₂₀. That is, model (1) has a branch of periodic solutions bifurcating from the positive equilibrium 𝑃₀ near 𝜏₂=𝜏₂₀.

Case 2 (𝜏₁>0 and 𝜏₂=0). According to a similar discussion as in Case 1, we can obtain the following result.

Theorem 2. There exists a 𝜏₁₀>0 such that the positive equilibrium 𝑃₀ of model (1) is locally asymptotically stable for 𝜏₁<𝜏₁₀ and undergoes a Hopf bifurcation when 𝜏₁=𝜏₁₀ given by 𝜏₁₀=1𝜔₁₀arccos−𝑢₁₂𝜔²₁₀+𝑢₄−𝑢₁₂𝑢₄+𝑢₂𝑢₃𝑢₁₁𝑢²₄+𝜔²₁₀, (18) where 𝑖𝜔₁₀ is the root of corresponding characteristic equation.

Case 3 (𝜏₁>0 and 𝜏₂∈(0,𝜏₂₀)). This case states that 𝜏₁ is regarded as a parameter and 𝜏₂ is in its stable interval. Assume that (8) has purely imaginary solution of the form 𝜆=𝑖𝜔. Substituting it into (8) and separating the real and imaginary parts, we have −𝜔²+𝑢₁₂𝑢₄−𝑢₂𝑢₃cos𝜔𝜏₂=−𝑢₁₁𝑢₄cos𝜔𝜏₁+𝑢₁₁𝜔sin𝜔𝜏₁,−𝑢₁₂+𝑢₄𝜔+𝑢₂𝑢₃sin𝜔𝜏₂=𝑢₁₁𝜔cos𝜔𝜏₁+𝑢₁₁𝑢₄sin𝜔𝜏₁. (19) Eliminating 𝜏₁ leads to 𝜔⁴+𝑢₄+𝑢₁₂²−2𝑢₁₂𝑢₄−𝑢₂𝑢₃cos𝜔𝜏₂−𝑢²₁₁𝜔²−2𝑢₂𝑢₃𝑢₄+𝑢₁₂⋅sin𝜔𝜏₂𝜔+𝑢²₁₂𝑢²₄+𝑢²₂𝑢²₃−𝑢²₁₁𝑢²₄−2𝑢₁₂𝑢₂𝑢₃𝑢₄cos𝜔𝜏₂=0. (20) It can be seen that there exists at least one real positive root 𝜔^′₁₀ for (20) if the condition (𝑢₁₂𝑢₄−𝑢₂𝑢₃)²<𝑢²₁₁𝑢²₄ holds. Equation (19) can also be written as −𝑢₂𝑢₃cos𝜔𝜏₂=−𝑢₁₁𝑢₄cos𝜔𝜏₁+𝑢₁₁𝜔sin𝜔𝜏₁+𝜔²−𝑢₁₂𝑢₄,𝑢₂𝑢₃sin𝜔𝜏₂=𝑢₁₁𝜔cos𝜔𝜏₁+𝑢₁₁𝑢₄sin𝜔𝜏₁+𝑢₁₂+𝑢₄𝜔. (21) Equation (21) can be solved as 𝜏^′_1𝑘=1𝜔^′₁₀×arccos𝜔²−𝑢₁₂𝑢₄²+𝑢₄+𝑢₁₂²𝜔²+𝑢²₁₁𝑢²₄+𝑢²₁₁𝜔²−𝑢²₂𝑢²₃×−2𝑢₁₁𝑢²₄+𝜔²×𝜔²−𝑢₁₂𝑢₄²+𝜔²𝑢₄+𝑢₁₂²⁻¹+𝜙₁+𝜙₂+2𝑘𝜋𝜔^′₁₀,𝑘=0,1,2,…, (22) where 𝜙₁=arctan(𝑢₄/𝜔^′₁₀) and 𝜙₂=arctan((𝜔^′2₁₀−𝑢₁₂𝑢₄)/𝜔^′₁₀(𝑢₄+𝑢₁₂)).

Now, differentiating (8) with respect to 𝜏₁ and substituting the eigenvalue 𝑖𝜔^′₁₀ and time delay 𝜏₁=𝜏^′₁₀, it follows that 𝑀𝑑(Re𝜆)𝑑𝜏_1𝜏^′₁₀+𝑁𝑑𝜔𝑑𝜏_1𝜏^′₁₀𝑑()=𝑃,−𝑁Re𝜆𝑑𝜏_1𝜏^′₁₀+𝑀𝑑𝜔𝑑𝜏_1𝜏^′₁₀=𝑄, (23) where 𝑀=−𝑢₄+𝑢₁₂−𝑢₁₁cos𝜔^′₁₀𝜏^′₁₀+𝜏^′₁₀𝑢₁₁−𝑢₄cos𝜔^′₁₀𝜏^′₁₀+𝜔^′₁₀sin𝜔^′₁₀𝜏^′₁₀+𝜏₂𝑢₂𝑢₃cos𝜔^′₁₀𝜏₂,𝑁=−2𝜔^′₁₀−𝑢₁₁sin𝜔^′₁₀𝜏^′₁₀−𝜏^′₁₀𝑢₁₁𝜔^′₁₀cos𝜔^′₁₀𝜏^′₁₀+𝑢₄sin𝜔^′₁₀𝜏^′₁₀+𝜏₂𝑢₂𝑢₃sin𝜔^′₁₀𝜏₂,𝑃=𝑢₁₁𝜔^′₁₀𝜔^′₁₀cos𝜔^′₁₀𝜏^′₁₀+𝑢₄sin𝜔^′₁₀𝜏^′₁₀,𝑄=𝑢₁₁𝜔^′₁₀−𝜔^′₁₀sin𝜔^′₁₀𝜏^′₁₀+𝑢₄cos𝜔^′₁₀𝜏^′₁₀. (24)(25)

Solving (23), we obtain 𝑑(Re𝜆)𝑑𝜏_1𝜏^′₁₀=𝑀𝑃−𝑁𝑄𝑀²+𝑁²≠0. (26)

Theorem 3. Assume (𝑢₁₂𝑢₄−𝑢₂𝑢₃)²<𝑢²₁₁𝑢²₄ and 𝜏₂∈[0,𝜏₂₀). The positive equilibrium 𝑃₀ of model (1) is asymptotically stable for 𝜏₁∈[0,𝜏^′₁₀) and undergoes Hopf bifurcation at 𝜏₁=𝜏^′₁₀, where 𝜏^′₁₀=1𝜔^′₁₀×arccos𝜔^′2₁₀−𝑢₁₂𝑢₄²+𝑢₄+𝑢₁₂²𝜔^′2₁₀+𝑢²₁₁𝑢²₄+𝑢²₁₁𝜔^′2₁₀−𝑢²₂𝑢²₃×−2𝑢₁₁𝑢²₄+𝜔^′2₁₀×𝜔^′2₁₀−𝑢₁₂𝑢₄²+𝜔^′2₁₀𝑢₄+𝑢₁₂²⁻¹+𝜙₁+𝜙₂, (27) where 𝜙₁=arctan(𝑢₄/𝜔^′₁₀) and 𝜙₂=arctan((𝜔^′2₁₀−𝑢₁₂𝑢₄)/𝜔^′₁₀(𝑢₄+𝑢₁₂)).

Case 4 (𝜏₁∈(0,𝜏₁₀) and 𝜏₂>0). According to a similar discussion as in Case 3, we can obtain the following result.

Theorem 4. The positive equilibrium 𝑃₀ of model (1) is asymptotically stable for 𝜏₂∈[0,𝜏^′₂₀) and undergoes Hopf bifurcation at 𝜏₂=𝜏^′₂₀, where 𝜏^′₂₀=1𝜔^′₂₀×arccos𝜔^′2₂₀−𝑢₁₂𝑢₄²+𝑢₄+𝑢₁₂²𝜔^′2₂₀+𝑢²₁₂𝑢²₃−𝑢²₁₁𝑢²₄−𝑢²₁₁𝜔^′2₂₀×−2𝑢₂𝑢₃×𝑢₄+𝑢₁₂²𝜔^′2₂₀+𝜔^′2₂₀−𝑢₁₂𝑢₄²⁻¹−𝜓, (28) where 𝜓=arctan((𝑢₄+𝑢₁₂)𝜔^′₂₀/(𝜔^′2₂₀−𝑢₁₂𝑢₄)).

3. Direction and the Stability of Hopf Bifurcation

In the previous section, we obtain the conditions under which a family of periodic solutions bifurcate from the positive equilibrium at the critical values of time delays 𝜏₁ and 𝜏₂. As pointed out by Hassard et al., by employing the normal form and center manifold theory, the formulae for determining the directions, stability, and period of Hopf bifurcation can be presented. Following the ideas of Hassard et al., this section discusses the directions, stability, and period of Hopf bifurcation. These properties are studied with respect to 𝜏₁ for fixed 𝜏₂∈(0,𝜏₂₀). Throughout this section, it is considered that the system (1) undergoes Hopf bifurcation at 𝜏₁=𝜏^′₁₀, 𝜏^∗₂∈(0,𝜏₂₀) at the equilibrium 𝑃₀. Without loss of generality, this section assumes that 𝜏^∗₂<𝜏^′₁₀.

Let 𝑥₁=𝑦₁−𝑦^∗₁, 𝑥₂=𝑦₂−𝑦^∗₂, 𝑡=𝑡/𝜏₁, 𝑥_𝑖(𝑡)=𝑥_𝑖(𝜏₁𝑡), and 𝜏₁=𝜏^′₁₀+𝜇 and dropping the bars for simplification of notations, the system (1) can be written as a functional differential equation in 𝐶=𝐶([−1,0],𝑅²)𝑥̇(𝑡)=𝐿_𝜇𝑥_𝑡+𝐹𝜇,𝑥_𝑡, (29) where 𝑥(𝑡)=(𝑥₁(𝑡),𝑥₂(𝑡))^𝑇∈𝑅², and 𝐿_𝜇∶𝐶→𝑅, 𝐹∶𝑅×𝐶→𝑅 are given, respectively, by 𝐿_𝜇𝜙=𝜏^′₁₀×+𝜇𝐴𝑢₀𝐴+𝑢₀²+ℎ₁𝑢₀𝑣₀1+ℎ₁𝑢₀²−𝑢₀1+ℎ₁𝑢₀0𝑝𝐸₀𝑣₀𝑝𝑣₀(0)+−𝑐×𝜙00𝑏𝑣₀1+ℎ₁𝑢₀²0𝜙−𝜏^∗₂𝜏^′₁₀+−𝑢₀0𝜙(),𝐹=𝜏00−1𝜇,𝜙^′₁₀×𝐴+𝜇₂𝐴+𝑢₀³+ℎ₁𝑣₀1+ℎ₁𝑢₀³𝜙²₁(0)−𝜙₁(0)𝜙₁()−1−11+ℎ₁𝑢₀²𝜙₁(0)𝜙₂(0)−𝐴²𝐴+𝑢₀⁴+ℎ²₁𝑣₀1+ℎ₁𝑢₀⁴⋅𝜙³₁(0)+ℎ₁1+ℎ₁𝑢₀³𝜙²₁(0)𝜙₂(0)−+⋯𝑏ℎ₁𝑣₀1+ℎ₁𝑢₀³𝜙²₁−𝜏^∗₂𝜏^′₁₀+𝑏1+ℎ₁𝑢₀²𝜙₁−𝜏^∗₂𝜏^′₁₀𝜙₂(0)−𝑝𝐸₀𝑐𝑝𝑣₀−𝑐²𝜙²₂(0)+𝑏𝑣₀ℎ²₁1+ℎ₁𝑢₀⁴𝜙³₁−𝜏^∗₂𝜏^′₁₀−𝑏ℎ₁1+ℎ₁𝑢₀³𝜙²₁−𝜏^∗₂𝜏^′₁₀𝜙₂(0)+𝑚𝑝²𝑐𝑝𝑣₀𝜙−𝑐³²(0)+⋯, (30) where 𝜙(𝜃)=(𝜙₁(𝜃),𝜙₂(𝜃))^𝑇∈𝐶. By the Riesz representation theorem, there exists a 2×2 matrix 𝜂(𝜃,𝜇) whose components are bounded variation for 𝜃∈[−1,0], such that 𝐿_𝜇𝜙=⁰₋₁𝑑𝜂𝜃,𝜇𝜙(𝜃),for𝜙∈𝐶, (31) where we choose 𝜂=𝜏𝜃,𝜇_′10×+𝜇𝐴𝑢₀𝐴+𝑢₀²+ℎ₁𝑢₀𝑣₀1+ℎ₁𝑢₀²−𝑢₀−𝑢₀1+ℎ₁𝑢₀𝑏𝑣₀1+ℎ₁𝑢₀²𝑝𝐸₀𝑣₀𝑝𝑣₀,𝜏−𝑐𝜃=0,_′10+𝜇−𝑢₀0𝑏𝑣₀1+ℎ₁𝑢₀²0−𝜏,𝜃∈^∗2𝜏^′₁₀,𝜏,0_′10+𝜇−𝑢₀0𝜏00,𝜃∈−1,−^∗2𝜏^′₁₀,0_2×2,𝜃=−1. (32) For 𝜙∈𝐶¹([−1,0],𝑅²), define 𝐴𝜇(𝜃)𝜙=𝑑𝜙[),∫𝑑𝜃,𝜃∈−1,0₀₋₁𝜙(𝑠)𝑅𝜇[),𝐹𝑑𝜂𝑠,𝜇,𝜃=0,𝜙=0,𝜃∈−1,0𝜇,𝜙,𝜃=0. (33) Then, system (29) is equivalent to 𝑥̇_𝑡=𝐴𝜇𝑥_𝑡+𝑅𝜇𝑥_𝑡, (34) where 𝑥_𝑡(𝜃)=𝑥(𝑡+𝜃) for 𝜃∈[−1,0].

For 𝜓∈𝐶¹([0,1],(𝑅²)^∗), define 𝐴^∗𝜓(𝑠)=−𝑑𝜓(𝑠)(],∫𝑑𝑠,𝑠∈0,1₀₋₁𝑑𝜂^𝑇(𝑡,0)𝜓(−𝑡),𝑠=0, (35) and bilinear form 𝜓(𝑠),𝜙(𝜃)=𝜓(0)𝜙(0)−⁰₋₁^𝜃_𝜉=0𝜓(𝜉−𝜃)𝑑𝜂(𝜃)𝜙(𝜉)𝑑𝜉, (36) where 𝜂(𝜃)=𝜂(𝜃,0), 𝐴=𝐴(0), and 𝐴^∗ are adjoint operators. Since ±𝑖𝜔^′₁₀𝜏^′₁₀ are eigenvalues of 𝐴(0), they are also eigenvalues of 𝐴^∗. Next, we need to compute the eigenvector of 𝐴(0) and 𝐴^∗ corresponding to 𝑖𝜔^′₁₀𝜏^′₁₀ and −𝑖𝜔^′₁₀𝜏^′₁₀, respectively.

In terms of the discussion mentioned above, we see that ±𝑖𝜔^′₁₀𝜏^′₁₀ are eigenvalues of 𝐴(0) and 𝐴^∗. Next, we calculate the eigenvector 𝑞(𝜃) of 𝐴(0) corresponding to 𝑖𝜔^′₁₀𝜏^′₁₀ and the eigenvector 𝑞^∗(𝑠) of 𝐴^∗ belonging to −𝑖𝜔^′₁₀𝜏^′₁₀, respectively. By the definition of 𝐴(0) and 𝐴^∗, it is not difficult to obtain that 𝛼=𝑖𝜔^′₁₀+𝑢₀𝑒^{−𝑖𝜔′10𝜏′10}−𝐴𝑢₀𝐴+𝑢₀²−ℎ₁𝑢₀𝑣₀1+ℎ₁𝑢₀²×−𝑢₀1+ℎ₁𝑢₀⁻¹. (37) Similarly, let 𝑞^∗(𝑠)=𝐷(1,𝛼^∗)𝑒^{𝑖𝜔′10𝜏′𝑠10} be the eigenvector of 𝐴^∗ corresponding to −𝑖𝜔^′₁₀𝜏^′₁₀, and we can obtain 𝛼^∗=−𝑖𝜔^′₁₀+𝑢₀𝑒^{𝑖𝜔′10𝜏′10}−𝐴𝑢₀𝐴+𝑢₀²−ℎ₁𝑢₀𝑣₀1+ℎ₁𝑢₀²×𝑏𝑣₀1+ℎ₁𝑢₀²𝑒^{𝑖𝜔′10𝜏∗2−1}. (38) It follows from (36) that 𝑞^∗(𝑠)(𝜃=,𝑞)𝐿1,𝛼^∗(1,𝛼)^𝑇−⁰₋₁^𝜃_𝜉=0𝐿1,𝛼^∗𝑒^{−𝑖𝜔′10𝜏′10(𝜉−𝜃)}𝑑𝜂(𝜃)(1,𝛼)^𝑇𝑒^{𝑖𝜔′10𝜏′10𝜉}=𝑑𝜉𝐿1+𝛼^∗𝛼−𝑢₀𝑒^{−𝑖𝜔′10𝜏′10}+𝑏𝑣₀𝜏^∗₂𝛼^∗𝜏^′₁₀1+ℎ₁𝑢₀²𝑒^{−𝑖𝜔′10𝜏∗2}. (39) Then, we obtain 𝐿=1+𝛼^∗𝛼−𝑢₀𝑒^{−𝑖𝜔′10𝜏′10}+𝑏𝑣₀𝜏^∗₂𝛼^∗𝜏^′₁₀1+ℎ₁𝑢₀²𝑒^{−𝑖𝜔′10𝜏∗2−1}, (40) such that ⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=1, ⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=0.

Based on the ideas in Hassard et al. [21], we can obtain 𝑔(𝑧,𝑧)=𝑞^∗(0)𝐹₀(𝑧,𝑧) and the following expressions: 𝑥_1𝑡(0)=𝑊⁽¹⁾(𝑡,0)+𝑧+𝑥𝑧,_2𝑡(0)=𝑊⁽²⁾(𝑡,0)+𝛼𝑧+𝛼𝑥𝑧,_1𝑡−𝜏^∗₂𝜏^′₁₀=𝑊⁽¹⁾𝜏𝑡,−^∗2𝜏^′₁₀+𝑧𝑒^{−𝑖𝜔′10𝜏∗2}+𝑧𝑒^{𝑖𝜔′10𝜏∗2},𝑥_2𝑡−𝜏^∗₂𝜏^′₁₀=𝑊⁽²⁾𝜏𝑡,−^∗2𝜏^′₁₀+𝛼𝑧𝑒^{−𝑖𝜔′10𝜏∗2}+𝛼𝑧𝑒^{𝑖𝜔′10𝜏∗2},𝑥_1𝑡(−1)=𝑊⁽¹⁾(𝑡,−1)+𝑧𝑒^{−𝑖𝜔′10𝜏′10}+𝑧𝑒^{𝑖𝜔′10𝜏′10},𝑥_2𝑡(−1)=𝑊⁽²⁾(𝑡,−1)+𝛼𝑧𝑒^{−𝑖𝜔′10𝜏′10}+𝛼𝑧𝑒^{𝑖𝜔′10𝜏′10}, (41) where those important parameters can be computed as follows: 𝑔₂₀=2𝐷𝜏^′₁₀𝐴²𝐴+𝑢₀³+ℎ₁𝑣₀1+ℎ₁𝑢₀³−𝑒^{−𝑖𝜔′10𝜏′10}−𝛼1+ℎ₁𝑢₀²+𝛼^∗𝑏𝛼1+ℎ₁𝑢₀²𝑒^{−𝑖𝜔′10𝜏∗2}−𝑏ℎ₁𝑣₀1+ℎ₁𝑢₀³𝑒^{−2𝑖𝜔′10𝜏∗2}−𝑝𝐸₀𝑐𝛼²𝑝𝑣₀−𝑐²,𝑔₁₁=2𝐷𝜏^′₁₀𝐴²𝐴+𝑢₀³+ℎ₁𝑣₀1+ℎ₁𝑢₀³−Re𝑒^{𝑖𝜔′10𝜏′10}−Re{𝛼}1+ℎ₁𝑢₀²+𝛼^∗𝑏1+ℎ₁𝑢₀²Re𝛼𝑒^{𝑖𝜔′10𝜏∗2}−𝑏ℎ₁𝑣₀1+ℎ₁𝑢₀³−𝑝𝐸₀𝑐𝛼𝛼𝑝𝑣₀−𝑐²,𝑔₀₂=2𝐷𝜏^′₁₀𝐴²𝐴+𝑢₀³+ℎ₁𝑣₀1+ℎ₁𝑢₀³−𝑒^{𝑖𝜔′10𝜏′10}−𝛼1+ℎ₁𝑢₀²+𝛼^∗𝑏𝛼1+ℎ₁𝑢₀²𝑒^{𝑖𝜔′10𝜏∗2}−𝑏ℎ₁𝑣₀1+ℎ₁𝑢₀³𝑒^{2𝑖𝜔′10𝜏∗2}−𝑝𝐸₀𝑐𝛼²𝑝𝑣₀−𝑐²,𝑔₂₁=2𝐷𝜏^′₁₀𝐴²𝐴+𝑢₀³+ℎ₁𝑣₀1+ℎ₁𝑢₀³×2𝑊⁽¹⁾₁₁(0)+𝑊⁽¹⁾₂₀(0)−𝑊⁽¹⁾₁₁()+𝑊−1⁽¹⁾²⁰(−1)2+𝑊¹₂₀(0)2𝑒^{𝑖𝜔′10𝜏′10}+𝑊⁽¹⁾₁₁(0)𝑒^{−𝑖𝜔′10𝜏′10}−11+ℎ₁𝑢₀²×𝑊²₁₁(0)+𝑊⁽²⁾₂₀(0)2+𝑊⁽¹⁾₂₀(0)2𝛼+𝑊⁽¹⁾₁₁(0)𝛼𝐴−3²𝐴+𝑢₀⁴+𝑣₀ℎ²₁1+ℎ₁𝑢₀⁴+ℎ₁1+ℎ₁𝑢₀³()+𝛼+2𝛼𝛼^∗𝑏1+ℎ₁𝑢₀²×𝑊⁽²⁾₁₁(0)𝑒^{−𝑖𝜔′10𝜏∗2}+𝑊⁽²⁾₂₀(0)2𝑒^{𝑖𝜔′10𝜏∗2}+𝑊⁽¹⁾₂₀−𝜏^∗₂/𝜏^′₁₀2𝛼+𝑊⁽¹⁾₁₁−𝜏^∗₂𝜏^′₁₀𝛼−𝑏ℎ₁𝑣₀1+ℎ₁𝑢₀³×2𝑊⁽¹⁾₁₁−𝜏^∗₂𝜏^′₁₀𝑒^{−𝑖𝜔′10𝜏∗2}+𝑊⁽¹⁾₂₀−𝜏^∗₂𝜏^′₁₀𝑒^{𝑖𝜔′10𝜏∗2}−𝑝𝐸₀𝑐𝑝𝑣₀−𝑐²2𝑊⁽²⁾₁₁(0)𝛼+𝑊⁽²⁾₂₀(0)𝛼+3𝑏𝑣₀ℎ²₁1+ℎ₁𝑢₀⁴𝑒^{−𝑖𝜔′10𝜏∗2}−𝑏ℎ₁1+ℎ₁𝑢₀³×𝛼𝑒^{−2𝑖𝜔′10𝜏∗2}+2𝛼+3𝑚𝑝²𝑐𝛼²𝛼𝑝𝑣₀−𝑐⁴. (42) According to the computation process similar to Hassard’s, parameters of 𝑔₂₁ can be written as follows: 𝑆₁=𝐴²𝐴+𝑢₀³+ℎ₁𝑣₀1+ℎ₁𝑢₀³−𝑒^{−𝑖𝜔′10𝜏′10}−𝛼1+ℎ₁𝑢₀²𝑏𝛼1+ℎ₁𝑢₀²𝑒^{−𝑖𝜔′10𝜏∗2}−𝑏ℎ₁𝑣₀1+ℎ₁𝑢₀³𝑒^{−2𝑖𝜔′10𝜏∗2}−𝑝𝐸₀𝑐𝛼²𝑝𝑣₀−𝑐²,𝑆₂=𝐴²𝐴+𝑢₀³+ℎ₁𝑣₀1+ℎ₁𝑢₀³−Re𝑒^{𝑖𝜔′10𝜏′10}−Re(𝛼)1+ℎ₁𝑢₀²𝑏1+ℎ₁𝑢₀²Re𝛼𝑒^{𝑖𝜔10𝜏∗2}−𝑏ℎ₁𝑣₀1+ℎ₁𝑢₀³−𝑝𝐸₀𝑐𝛼𝛼𝑝𝑣₀−𝑐²,𝑅₁=22𝑖𝜔^′₁₀−𝐴𝑢₀𝐴+𝑢₀²−ℎ₁𝑢₀𝑣₀1+ℎ₁𝑢₀²+𝑢₀𝑒^{−2𝑖𝜔′10𝜏′10}𝑢₀1+ℎ₁𝑢₀−𝑏𝑣₀1+ℎ₁𝑢₀𝑒^{−2𝑖𝜔′10𝜏∗2}2𝑖𝜔^′₁₀−𝑝𝐸₀𝑣₀𝑝𝑣₀−𝑐⁻¹𝑆₁,𝑅₂𝑢=2⁰−𝐴𝑢₀𝐴+𝑢₀²−ℎ₁𝑢₀𝑣₀1+ℎ₁𝑢₀²𝑢₀1+ℎ₁𝑢₀−𝑏𝑣₀1+ℎ₁𝑢₀²−𝑝𝐸₀𝑣₀𝑝𝑣₀−𝑐⁻¹𝑆₂,𝑊₂₀(𝜃)=𝑖𝑔₂₀𝜔^′₁₀𝜏^′₁₀𝑞(0)𝑒^{𝑖𝜔′10𝜏′10𝜃}+𝑖𝑔₀₂3𝜔^′₁₀𝜏^′₁₀𝑞(0)𝑒^{−𝑖𝜔′10𝜏′10𝜃}+𝑅₁𝑒^{2𝑖𝜔′10𝜏′10𝜃},𝑊₁₁(𝜃)=−𝑖𝑔₁₁𝜔^′₁₀𝜏^′₁₀𝑞(0)𝑒^{𝑖𝜔′10𝜏′10𝜃}+𝑖𝑔₁₁𝜔^′₁₀𝜏^′₁₀𝑞(0)𝑒^{−𝑖𝜔′10𝜏′10𝜃}+𝑅₂. (43)

Therefore, we can compute 𝑔₂₁ by the parameters and delays. Then, the following results can be expressed as 𝑐₁(0)=𝑖2𝜔^′₁₀𝜏^′₁₀𝑔₂₀𝑔₁₁−2𝑔₁₁²−13𝑔₀₂²+𝑔₂₁2,𝜇₂𝑐=−Re₁(0)Re𝜆^′𝜏^′₁₀,𝛽₂=2Re𝑐₁(0,𝑇₂)𝑐=−Im₁(0)+𝜇₂Im𝜆^′𝜏^′₁₀𝜔^′₁₀𝜏^′₁₀. (44) These parameters give a description of the bifurcating periodic solutions in the center manifold of system (1) at critical values 𝜏₁=𝜏^′₁₀, which can be expressed as follows.(i)𝜇₂ determines the direction of the Hopf bifurcation: if 𝜇₂>0(𝜇₂<0), the Hopf bifurcation is supercritical (subcritical).(ii)𝛽₂ gives the stability of periodic solution: if 𝛽₂<0(𝛽₂>0), the periodic solution is stable (unstable).(iii)𝑇₂ expresses the period of the bifurcating periodic solutions: if 𝑇₂>0(𝑇₂<0), the period of the bifurcating periodic solution increases (decreases).

4. Numerical Simulation

With the help of MATLAB, we present some numerical simulations to verify and extend our theoretical analysis proved in Sections 2 and 3. The parameters of model (1) are chosen as 𝐴=0.02, 𝑑₁=0.01, ℎ₁=1.25, 𝑏=7.5, 𝑑₂=2, 𝑝=35, 𝑐=1, and 𝑚=0.5.

According to the parameters given above, we can obtain the positive equilibrium point 𝑃₀(0.4055,0.8099,0.0183) of model (1). It is easy to obtain that 𝑃₀ is asymptotically stable in the absence of delay (see Figure 1). Moreover, since 𝑢₁𝑢₄+𝑢₂𝑢₃=−0.7233<0, there is a positive root 𝜔₂₀=0.8389 for (13). Hence, when 𝜏₁=0, the critical value of the delay corresponding to 𝜔₂₀ is 𝜏₂₀=0.2247, which is shown in Figure 2. When 𝜏₂=0, by a similar computation, we can obtain that the critical delay for 𝜏₁ is 𝜏₁₀=1.272.

(a)

(b)

(c)

(d)

Figure 1 

Without any time delay, the positive equilibrium point 𝑃0 is asymptotically stable.

Figure 2 

When 𝜏2=0.288>𝜏20, bifurcating periodic solutions from the equilibrium point 𝑃0 occur.

Based on the results of Cases 3 and 4 in Section 2, we draw a critical curve 𝜏₂=𝑓(𝜏₁) with respect to two parameters 𝜏₁ and 𝜏₂ (see Figure 3). The main algorithm for drawing the critical curve is searching the values 𝜏^′₁₀(𝜏^′₂₀) for fixed positive values of 𝜏₂∈(0,𝜏₂₀)  (𝜏₁∈(0,𝜏₁₀)). From Figure 3, it is clear that the domain surrounded by the critical curve 𝑓(𝜏₁) and the two axes, that is, the gray region, is stable for model (1).

Figure 3 

The critical curve of stable domain with respect to 𝜏1 and 𝜏2.

In addition, according to the algorithm derived in Section 3, we can obtain the following values: 𝑐₁(0)=−1.9243−3.0887𝑖, 𝜔^′₁₀=0.7538, 𝜏^′₁₀=0.6546, 𝜆^′(𝜏^′₁₀)=0.0879+0.1826𝑖, 𝜇₂=21.8851, 𝛽₂=−3.8487, and 𝑇₂=−1.8402. By the discussion in Section 3, we know that model (1) can undergo a supercritical Hopf bifurcation at the positive equilibrium 𝑃₀ and the bifurcating periodic solution occurs when 𝜏₁ crosses 𝜏^′₁₀ to the right with fixed 𝜏^∗₂∈(0,𝜏₂₀) and the bifurcating periodic solution is stable. Figure 4 plots the bifurcating periodic solution of the model (1) at 𝜏₁=𝜏^′₁₀=0.6546 and 𝜏^∗₂=0.15.

(a)

(b)

(c)

(d)

Figure 4 

When at 𝜏1=𝜏′10=0.6546 and 𝜏∗2=0.15, bifurcating periodic solutions of model (1) from 𝑃0 occur.

To explore the possibility of occurrence of chaos, bifurcation diagrams are plotted for the key parameter 𝜏₁. Bifurcation diagram in Figure 5 of model (1) in (𝜏₁−𝑥) plane is given for the fixed 𝜏₂=0.15. We can see that model (1) experiences the processes of periodic, periodic doubling cascade and chaos in the region 0.98<𝜏₁<1.018. When 0.98<𝜏₁<0.996, model (1) exhibits a stable limit cycle. When 𝜏₁ is increasing, the limit cycle is unstable and there is a cascade of periodic doubling bifurcation (see Figure 6) leading to chaos.

Figure 5 

Bifurcation diagram of model (1) showing the effect of 𝜏1 on the dynamical behavior.

(a)

(b)

(c)

Figure 6 

The projection of the attracting set of model (1) in the 𝑥-𝑦 plane for (a) 𝜏1=1.01, (b) 𝜏1=1.0125, and (c) 𝜏1=1.017, respectively.

5. Conclusions

In this paper, dynamical complexity of a multiple-delayed predator-prey model with Allee effect is analyzed by using differential-algebraic system theory. The prey-predator model with single delay has been investigated by many researchers (see, e.g., [22–25]). However, to our best knowledge, there are few papers on the bifurcations of differential-algebraic population dynamics with two or multiple delays. On the other hand, it should be noted that almost the existing bioeconomic models (see [2–5, 14, 15]) only investigate the simplest case of a logistic prey growth function. Compared with these works, the introduction of Allee effect makes the work studied in this paper novel.

This paper provides a new and efficient method for the qualitative analysis of the Hopf bifurcation of differential-algebraic biological economic system with multiple delays. From the analysis of the proposed model, we have obtained some interesting and useful results. The main contribution of this paper is that the stability domain with respect to two delays is defined and plotted. By choosing the time delay 𝜏₁ as a bifurcation parameter, the direction of Hopf bifurcation and stability of the bifurcating periodic orbit are discussed by using normal form and center manifold. Moreover, numerical simulations are carried out which substantiate the theoretical results. At the same time, bifurcation diagrams and attracting sets are plotted, which validate the existence of chaos due to periodic doubling bifurcation. However, those theoretical analyses, such as the existence of chaos and conditions of periodic doubling bifurcation, will be the topics of the future research.

Conflict of Interests

The authors do not have conflict of interests with any commercial identities.

Acknowledgments

This work was supported by the National Science Foundation of China (61273008), Doctor Startup Fund of Liaoning Province (20131026), Fundamental Research Funds for the Central University (N120405009), and Soft Science Research Project of Hubei Province (2012GDA01309). The authors also gratefully thank those anonymous authors whose work largely constitutes this sample paper.