Abstract

A class of Lotka-Volterra mutualistic system with time delays of benefit and feedback delays is introduced. By analyzing the associated characteristic equation, the local stability of the positive equilibrium and existence of Hopf bifurcation are obtained under all possible combinations of two or three delays selecting from multiple delays. Not only explicit formulas to determine the properties of the Hopf bifurcation are shown by using the normal form method and center manifold theorem, but also the global continuation of Hopf bifurcation is investigated by applying a global Hopf bifurcation result due to Wu (1998). Numerical simulations are given to support the theoretical results.

1. Introduction

In recent years, population models have merited a great deal of attention due to their theoretical and practical significance since the pioneering theoretical works by Lotka [1] and Volterra [2]; see [36] and references therein. Generally speaking, there are three kinds of fundamental forms of the interactions between two species such as competition, cooperation, and prey-predation in population biology. Among these interactions, the competition mechanism has been paying extreme attention because it possesses very significant function as a kind of restriction factor in the process of evolvement of biology.

It is well known that time delays of one type or another have been incorporated into mathematical models of population dynamics due to maturation time, capturing time, or other reasons. The effect of the past history on the stability of the system is an important problem in population biology. In general, delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and cause the population to fluctuate. In 2002, Jin and Ma [3] investigated the competition model with four delays: where are maturation times and are called the hunting delays for the prey and the predator species, respectively. The result was that time delay in the competitive Lotka-Volterra system does not harm the boundedness and persistence. Song et al. [4] analyzed the stability of the interior positive equilibrium and the existence of local Hopf bifurcation for system (1) with by taking the delay as the bifurcation parameter. Their results showed that changes of hunting delays for system (1) do not lead to the occurrence of Hopf bifurcation when interspecies competition is weaker than intraspecies competition. Zhang et al. [5] studied the stability and Hopf bifurcation of system (1) with . Zhang [6] also investigated the dynamics of system (1) only when and .

In fact, predator-prey system with time delays has also been investigated by lots of authors [711]. Faria [7] considered the delayed predator-prey system of the form where and can be interpreted as the population densities of the prey and the predator at the time , respectively, and and denote the hunting delay and the time of predator maturation. In [7], taking the single delay, , as a parameter and assuming the ratio to be constant, the author studied the properties of the local Hopf bifurcation. Song et al. [8] obtained the properties of the local Hopf bifurcation as well as the global existence of periodic solutions by choosing the sum as a bifurcation parameter. Yan and Li [9] considered the following delayed predator-prey system: where denotes the feedback time delay of prey species to the growth of species itself, and also for the predator. They found that the unique positive equilibrium of system (3) will no longer be absolutely stable and the switches from stability to instability to stability disappear as the feedback time delay increases monotonously from zero, which had been obtained for system (3) by Song and Wei [12].

However, the research for cooperative systems with time delays is still relatively little because delays in mutualistic systems usually deprive the boundedness and persistence in [13]. But time delays in prey-predator and competitive systems do not harm these properties [3, 14, 15]. Two species cohabit a common habitat and each species enhances the average growth rate of the other; this type of ecological interaction is known as facultative mutualism. Mutual phenomenon can increase viability and make species persistently multiply. As far as we know, the research on mutual system is less than prey-predator system and competitive system. At present, the known results of mutual system mainly focus on stability and persistence [16, 17]. Research on the ecologic system stability of positive equilibrium and existence of periodic solutions is very crucial, which can help us to realize the law for species quantity and predict the trend for species quantity. In 1997, He and Gopalsamy [17] considered the Lotka-Volterra mutualistic system with delay and obtained the stability of the positive equilibrium and the existence of Hopf bifurcation when . Meng and Wei [18] studied the following system with delay: where and are the rates of transmission between two species. They found that there are stability switches and Hopf bifurcation occurring when the delay passes through a sequence of critical values.

In general, the delays appearing in different terms of an ecological system are not equal. Therefore, it is more realistic to consider dynamical system with different delays. Motivated by the references [6, 9, 17], we consider the following two-species Lotka-Volterra mutualistic system with multiple delays: where and are the densities of two species at the time and and are the intrinsic growth rates of the two species. is the feedback delay of one species to its grow and is the feedback delay of the other species; and are the time delays of the benefit. The coefficients are all positive constants and . In addition, we would like to mention the work of F. Q. Zhang and Y. J. Zhang [19], who considered the system (6) only with and . In fact, feedback delay on different species should not be same, and the benefit period should not be equal to the feedback time delay. Thus, how does every delay or the combined two or three delays of the above four delays have an effect on the dynamics of system (6)? The purpose of this paper is to answer this question partially.

This paper is organized as follows. In Section 2, by analyzing the characteristic equation of linearized system of system (6) at the positive equilibrium, some sufficient conditions ensuring the local stability of the positive equilibrium and the existence of Hopf bifurcation are obtained. Some explicit formulas determining the direction and stability of periodic solutions bifurcating from Hopf bifurcations are given by applying the normal form method and center manifold theory due to Hassard et al. [20] in Section 3. In Section 4, we show the global existence of the periodic solutions due to a global Hopf bifurcation result of Wu [21] for functional differential equations. To support our theoretical predictions, some numerical simulations are included in Section 5. A brief discussion is also given in the last section.

2. Local Stability and Hopf Bifurcation

For convenience, letting , and dropping the bars for simplification of notation, system (6) is transformed into where , , , and .

From the point view of biological meaning, we are interested in the positive equilibrium. It is obvious that system (7) has a unique positive equilibrium defined by provided that the following condition(H1)is satisfied.

Next, we will consider the stability of the positive equilibrium and the existence of Hopf bifurcation.

Letting , , system (7) can be transformed into Linearizing system (9) at and rewriting as , system (9) is rewritten as where , , , and .

The associated characteristic equation of system (10) with is It is obvious that is not the root of (11). All roots of (11) have negative real parts since and if (H1) holds. So, the equilibrium point is locally asymptotically stable when (H1) holds.

2.1. Only Considering and

The characteristic equation of system (10) with is In the following, we will discuss the distribution of roots of (12) while and are given different values.

When and , (12) becomes (11). Thus, the positive equilibrium is locally asymptotically stable when (H1) holds.

Case  1A. Consider , . Equation (12) can be written in the form

We first introduce the following result which was proved by Ruan and Wei [22] using Rouch’s theorem.

Lemma 1. Consider the exponential polynomial where , are constants. As vary, the sum of the order of the zeros of on the open right half-plane can change only if a zero appears on or crosses the imaginary axis.

Let be the root of (13). Substituting it into (13) and separating the real and imaginary parts, we have It follows that Equation (16) does not have positive root since and . Thus, Hopf bifurcation does not occur at the positive equilibrium of system (7).

Case  1B. Consider , . Equation (12) is written in the form Equation (17) is similar to (13). Thus, Hopf bifurcation also does not occur at the positive equilibrium of system (7). We omit the corresponding proof.

Case  1C. Consider . Equation (12) is written in the form Equation (18) is also similar to (16). Thus, Hopf bifurcation does not occur at the positive equilibrium of system (7).

Case  1D. Consider and . Let be a root of (12). We can obtain the following form by giving the value and regarding as a parameter: It follows that It is obvious that Hopf bifurcation does not occur at the positive equilibrium of system (7) according to (16).

2.2. Only Considering and

The characteristic equation of system (10) with is In the following, we will discuss the distribution of roots of (21) while and are given different values.

When , , (21) becomes (11) and the equilibrium is locally asymptotically stable when (H1) holds.

Case  2A. Consider , . Equation (21) is written as follows:

Let be the root of (22), and we have that It follows that Since , (24) has one positive root, defined by .

From (23), if we denote then are a pair of purely imaginary root of (22) with .

Define . Let be the root of (22) near satisfying , . We first check whether the transversality condition is satisfied.

Lemma 2. the following transversality condition is satisfied:

Proof. This will show that there exists at least one eigenvalue with positive real part for . Moreover, the conditions for the existence of a Hopf bifurcation [15, 22] are satisfied to yield a periodic solution. Differentiating (22) with respect to , it follows that which implies that
For simplifying, define as and as , and we can obtain Thus, the transversality condition holds and Hopf bifurcation occurs at . We have the following theorem.

Theorem 3. For system (7), if (H1) holds, then there exists a positive number such that the coexistence equilibrium is locally asymptotically stable for and unstable for . Further, system (7) undergoes a Hopf bifurcation at the equilibrium for .

Case  2B. Consider , . Equation (21) can be written in the form Since (30) has the similar form of (22), the corresponding results are omitted.

Case  2C. Consider . Equation (21) can be written in the form Multiplying the into the both sides of (31) and letting to be the root of (31), we have that

Suppose that . We have when (H1) holds. From the second equation of (32), we can get . Thus, we can get , . Define , . Let be the root of (31) near , satisfying , . Further, (31) has a pair of purely imaginary roots .

If , then . If , then the second equation of (32) becomes If , (33) has one positive root . Thus, we can get , . Define , . Let be the root of (31) near , satisfying , . Further, (31) has a pair of purely imaginary roots . If , and (33) has two positive roots . It is obvious that . From (32), we denote , . When (resp., ), then (resp., ) are a pair of purely imaginary roots of (31). Further, we can check that there exists a positive integral such that and . Let be the root of (31) near , satisfying , . If , then the second equation of (32) becomes , which does not have positive root.

Lemma 4. Suppose that (H1) holds.(i)If , then all roots of (31) have strictly negative real parts; if , then all roots of (31), except for ±, have strictly negative real parts.(ii)If and , then all roots of (31) have strictly negative real parts; if , then all roots of (31), except for , have strictly negative real parts.(iii)If , there exist two sequences and . Further, there exists , when , and all the roots of (31) have negative real part; when , (31) has at least one root with positive real part. When and , (31) has a pair of pure imaginary roots.

Lemma 5. The following transversality condition is satisfied:

Proof. Multiplying the into the both sides of (31) and differentiating (31) with respect to , we can obtain that For simplifying, define as and as , and we can obtain Thus, the transversality condition is satisfied. This ends the proof.

Remark 6. When , the transversality condition is satisfied as follows:

Theorem 7. For system (7), consider the following.(i)If (H1) (resp., (H1) and ) holds, then there exists a positive number (resp., ) such that the positive equilibrium is locally asymptotically stable for (resp., ) and unstable for (resp., ). Further, system (7) undergoes a Hopf bifurcation at the positive equilibrium for   .(ii)If (H1) and hold, then there is a positive integer , such that the positive equilibrium switches times from stability to instability to stability; that is, the positive equilibrium of system (7) is locally asymptotically stable when and unstable when . Hopf bifurcation occurs at the positive equilibrium when and .

Case  2D. Consider and . Regarding as a parameter, we consider (21) with in its stable interval. Without loss of generality, we consider system (7) under Case  2A. Let be a root of (21), and then we obtain, by selecting delay , Eliminating the parameter , (38) leads to where and .

Denote . It is obvious that and when (H1) holds. Without loss of generality, we assume that (39) has at least finite positive roots, which are defined by . From (38), for every fixed , we have

Define , . Equation (39) has a pair of purely imaginary roots for .

In the following, differentiating Equation (21) with respect to and substituting , we can get where

If the following assumption(H2)holds, then the following results on stability and bifurcation of system (7) are obtained by the general Hopf bifurcation theorem for FDEs in Hale [14].

Theorem 8. For system (7), suppose that (H1) and (H2) hold and . Then, system (7) undergoes a Hopf bifurcation at the positive equilibrium when . That is, system (7) has a branch of periodic solution bifurcation from the zero solution near .

2.3. Only Considering and

The characteristic equation of system (7) with is In the following, we will discuss the distribution of roots of (43) while and are given different values.

When , , (43) becomes (11) and the equilibrium is locally asymptotically stable when (H1) holds.

Case  3A. Consider . Equation (43) can be written in the form Equation (44) is the same as (22), so we do not prove it.

Case  3B. Consider , . Equation (43) can be written in the form Equation (45) is the same as (17), so we also do not prove it.

Case  3C. Consider . Equation (43) is written in the form Suppose that is the root of (46) and we have that It obtains that It is obvious that (48) has only one positive root . From (47), we obtain that Equation (46) has a pair of purely imaginary roots . Define , . Let be the root of (46) near satisfying , .

Lemma 9. If the condition (H1) holds, then the following transversality condition is satisfied:

Proof. Differentiating (46) with respect to and noticing that is a function of , we can obtain that For simplifying, define as and as , and we can obtain Thus, if the condition (H1) holds, the transversality condition is satisfied. This ends the proof.

Theorem 10. For system (7), if (H1) holds, then there exists a positive number such that the interior equilibrium is locally asymptotically stable for and unstable for . Further, system (7) undergoes a Hopf bifurcation at the equilibrium for .

Case  3D. Consider and . Regarding as a parameter, we consider (43) with in its stable interval. Without loss of generality, we consider system (7) under Case  3A. Let be a root of (43). We can obtain the following results by selecting the value of delay : Eliminating the parameter , (53) leads to where , and .

Suppose that(H3)Equation (53) has at least finite positive roots.

If (H3) holds, the roots of (53) are defined by . From (53), for every fixed , there exists a sequence , where and then are a pair of purely imaginary roots of (43). Let , . When , (43) has a pair of purely imaginary roots for . Let be the root of (43) near satisfying , .

In the following, differentiating Equation (43) with respect to and substituting , we get where

Assume that(H4).

Therefore, by the general Hopf bifurcation theorem for FDEs in Hale [14], we have the following results on stability and bifurcation of system (7).

Theorem 11. For system (7), suppose that (H1), (H3), and (H4) hold and . Then, system (7) undergoes a Hopf bifurcation at the positive equilibrium when . That is, system (7) has a branch of periodic solutions bifurcating from the zero solution near .

Remark 12. If we only consider system (7) with delays , we have the characteristic equation which is similar to (43). Thus, we do not discuss the roots of (58). Similarly, if we only consider system (7) with delays or , we can obtain the corresponding characteristic equation which is also similar to (43). The detailed discussions are also omitted.

2.4. Considering and

We can obtain that the characteristic equation of system (7) with and is

When , (59) becomes (11) and the equilibrium is locally asymptotically stable when (H1) holds.

Multiplying into the both sides of (59) and letting be the root of (59), we have Since , we obtain that where , , , , , and . Without losses of generality, suppose that (61) has six different positive roots, defined by .

From (60), if we denote then are a pair of purely imaginary roots of (59) with .

Define , . Let be the root of (59) near satisfying , . Next, differentiating (59) with respect to , it follows that For simplifying, define as and as , and we can obtain where , , and .

Since , we suppose that(H5). If (H5) holds, then the transversality condition holds and a Hopf bifurcation occurs at . We have the following theorem.

Theorem 13. For system (7), if (H1) and (H5) hold, then there exists a positive number such that the positive equilibrium is locally asymptotically stable for and unstable for . Further, system (7) undergoes a Hopf bifurcation at the equilibrium for .

Further, the corresponding characteristic equation with and is

When , (66) becomes (11). The positive equilibrium is locally asymptotically stable when holds.

Regarding as a parameter, we consider (66) with in its stable interval. Without loss of generality, we consider system (7) under the above case. Multiplying into the both sides of (66) and letting be the root of (66), we have that by selecting the value of delay since , we obtain that where , , , and .

Suppose that(H6)Equation (68) has at least finite positive roots.If (H6) holds, then (68) has finite positive roots, given by . From (67), for every fixed , if there exists a sequence such that (68) holds, where then are a pair of purely imaginary root of (66).

Define , . Let be the root of (66) near satisfying , . Next, we will check whether the following transversality condition is true.

Differentiating (66) with respect to , it follows that For simplifying, define as and as , and we can obtain where , , , and .

Suppose that(H7). If (H7) holds, then the transversality condition holds and a Hopf bifurcation occurs at , . We have the following theorem.

Theorem 14. For system (7), if (H1), (H6), and (H7) hold, then there exists a positive number such that the positive equilibrium is locally asymptotically stable for and unstable for . Further, system (7) undergoes a Hopf bifurcation at the equilibrium for .

Remark 15. In addition, if we consider system (7) with and , we can obtain the same characteristic equation as (66) and the same result for the time delay . So, we omit it.

3. Direction and Stability of the Hopf Bifurcation

In this section, we will study the direction of the Hopf bifurcation and the stability of bifurcating periodic solution of system (7). The approach employed here is the normal form method and center manifold theorem introduced by Hassard et al. [20]. More precisely, we will compute the reduced system on the center manifold with the pair of conjugate complex, purely imaginary solution of the characteristic equation of system (10). By this reduction, we can determine the Hopf bifurcation direction, that is, to answer the question of whether the bifurcation branch of periodic solution exists locally for supercritical bifurcation or subcritical bifurcation.

Without loss of generality, we only consider the case with and . For convince, we denote any one of these critical values by , at which the characteristic equation has a pair of purely imaginary roots . Let ; then, is the Hopf bifurcation value of system (7). We first let , and normalize the delay with the scaling , and then (7) is transformed into an FDE in as where and are given, respectively, by where here , , , .

Turning to the linear problem by the Reisz representation theorem, there exists matrix-valued function such that

In fact, we can choose where denotes the Dirac delta function. Then, (78) is satisfied.

Next, for , we define the operator as Since , then system (73) is equivalent to the following operator equation: where , for , which is an equation of the form we desired.

For , we further define the adjoint of as and a bilinear form where . Then, and are adjoint operators. From the above analysis, we obtain that are the eigenvalues of and therefore they are also the eigenvalues of . Let be eigenvector of corresponding to and let be eigenvector of corresponding to , and then we have

By some simple computation, we can obtain where Then, , .

In the remainder of this section, by using the same notations as in Hassard et al. [20] and a computation process similar to that in [11], we obtain the coefficients determining the important quantities of the periodic solution: where

However, where and are both two-dimensional vectors and can be determined by

Furthermore, we can see that each in (87) is determined by parameters and delays in system (7). Thus, we can compute the following quantities: which determine the properties of bifurcation periodic solutions in the center manifold at the critical value . By the results of Hassard et al. [20], we have the following results.

Theorem 16. In (91), the following results hold.(i)The sign of determines the direction of the Hopf bifurcation: if (), then the Hopf bifurcation is supercritical (subcritical) and the bifurcation periodic solutions exist for ().(ii)The sign of determines the stability of the bifurcating periodic solutions: the bifurcation periodic solutions are stable (unstable) if ().(iii)The sign of determines the period of the bifurcating periodic solutions: the period increases (decreases) if ().

4. Global Continuation of Local Hopf Bifurcation

In this part, we will study the global continuation of periodic solutions bifurcating from the positive equilibrium , by applying a global Hopf bifurcation result due to Wu [21].

Let , ,, and , be defined as , for and . Since and denote the density of one species and the other species, respectively. The positive solution of system (7) is interior and its periodic solutions only arise in the first quadrant. Thus, throughout this section, we consider system (7) only in the domain .

It is easy to see that system (7) can be rewritten as the following functional differential equation: parameterized by two nonnegative real parameters , where . It is obvious that (92) has only an equilibria in under the assumption (H1), and it is easy to see that the mapping is completely continuous. If we restrict to the subspace of constant function , then is identified with and thus we obtain a mapping .

Let be the constant mapping with value . The point is called a stationary solution of (92) if . From system (7), we know easily that the following condition regarding the mapping holds:(A1).

It follows from system (7) that Then, under the assumption (H1), we have Therefore, we have the following condition on the linear operator :(A2) is an isomorphism at the equilibrium ; here, .

In addition, we can easily observe that the following result is true:(A3) is differential with respect to .

The characteristic matrix of (92) at a stationary solution , where , takes the following form: That is, The zeros of are called the characteristic roots. Note that (A2) is equivalent to which is not a characteristic root of any equilibrium of (92). Clearly, the characteristic matrix is continuous in .

A stationary solution of (92) is called a center if for some positive integer or has purely imaginary characteristic roots of the form for some positive integer . A center is said to be isolated if it is the only center in some neighborhood of .

From (96), we can see

Note that the above equation is the same as (31); therefore, it is easy to verify from the discussion regarding the local Hopf bifurcation in Section 2 that ,   is an isolated center and there exists a smooth curve : such that , for all and

For the above , we define the set It is easy to verify on that the following condition holds:(A4) if and only if , , and , .

Thus, the conditions (A1)–(A4) in [21] are satisfied. Moreover, if we define then we know on . Thus, the first crossing number of isolated center can be defined as follows:

In what follows, we define And let denote the connected component of in .

From the above discussion, we have where , in fact, takes the form , . Thus, the connected component through in is nonempty. Since the first crossing number of each center is always −1, by Theorem 3.3 [21], we conclude that is unbounded. Thus, we have proved the following lemma.

Lemma 17. is unbounded for each center .

Lemma 18. All nonconstant periodic solutions of system (7) with initial data , , , and , , are uniformly bounded when is bounded.

Proof. Suppose that are nonconstant periodic solutions of system (7) and define Then, it follows from system (7) that which implies that the solutions of system (7) cannot cross the -axis and -axis. Thus, the nonconstant periodic orbits must be located in the interior of each quadrant. The following results are deduced according to Theorem 4.9.1 in Kuang [15] and the standard comparison principle [14].
If is a solution of system (7) with , , then it is easy to say that there exists a such that, for all , . In addition, from the second equation of (7), we can get which leads to From the second equation of system (7), we also have which implies It follows from (104) and (107) that
On the other hand, we get, from the first equation of (7), It follows that
This shows that the nonconstant periodic solution of system (7) is uniformly bounded when is bounded. This completes the proof.

Lemma 19. System (7) has no nonconstant periodic solutions with period when the condition (H1) is satisfied.

Proof. Let’s suppose that there exist a contradiction in system (7); that is, system (7) has nonconstant periodic solution with period. Then, the following system (113) of ordinary differential equations has nonconstant periodic solution: which has the same equilibria as system (7), that is, a unique positive equilibrium . Note that -axis, and -axis are the invariable manifold of system (113) and the orbits of system (113) do not intersect with each other. Thus, there is no solution crosses the coordinate axes. On the other hand, note the fact that if system (113) has a periodic solution, then there must be equilibrium in its interior. Thus, we conclude that the periodic orbit of system (113) must lie in the first quadrant. It is well known that the positive equilibrium is global asymptotically stable to the first quadrant (see, e.g., [23, 24]). Thus, there is not periodic orbit in the first quadrant, too. This ends the proof.

In the following, we state and prove our main result in this section.

Theorem 20. Suppose that the condition (H1) holds. Then, for each , system (7) has at least periodic solutions.

Proof. It is sufficient to prove that the connected component onto is for each , where .
From the discussion in Section 2, we have Thus, one can get , for . From Lemma 19, we know that system (7) with has no nontrivial periodic solution. Consequently, the projection of onto is away from zero.
Suppose that the projection of onto is bounded, and then there exits such that the projection of onto in interval . Since and applying Lemma 18, one can obtain for . Therefore, the projection of onto is also bounded. Thus, we get, together with Lemma 18, that the connected component is bounded. This contradicts Lemma 17. The proof is completed.

5. Numerical Example

We will give some numerical results of system (7) to support the analytic results in this section. We consider the following system: and get the positive equilibrium . System (115) without any time delay is locally asymptotically stable (see Figure 1). Firstly, when only consider delays and , and we find that Hopf bifurcation do not occur at the positive equilibrium of system (115) with any value of them (see Figure 2).

Secondly, we have that the critical value of delay is when we only consider system (115) with delays and . From Theorem 3, we obtain the corresponding waveform shown in Figure 3. Similarly, we have that in Case  2B and in Case  2C, respectively. The corresponding waveforms are shown in Figure 4. Regarding as a parameter with , we can obtain . From Theorem 8, the positive equilibrium is locally asymptotically stable for and unstable for (see Figure 5).

Thirdly, we can obtain when only consider . That is, when increases from zero to the critical value , the positive equilibrium is locally asymptotically stable and then it loses its stability and a Hopf bifurcation occurs once (see Figure 6). Similarly, when we regard as a parameter and let , we have that . From Theorem 11, system (115) is locally asymptotically stable for and unstable for (see Figure 7). Though, we only give the associated characteristic equation (58) when only we consider delays and of system (115) with , we give the phase plots when the delay increases from zero to its critical value and the value of belongs to its stable interval (see Figure 8).

Finally, by algorithm (91) derived in Section 3, we have , , and . Thus, the Hopf bifurcation is supercritical, and the bifurcating periodic solutions are asymptotically stable, which is illustrated in Figure 9.

6. Discussion

He and Gopalsamy [17] considered a Lotka-Volterra mutualistic system (4) with delay, in which all delays take the same value; Meng and Wei only considered a mutualistic system (5) with the same value of feedback delay for each population in [18] while we not only introduce feedback delays of two populations to grow of the species in system (6) but also take time delay of each species benefit into account. By theory analysis and simulations, we find that time delays of the benefit and can not affect the stability of the system (6) when we only choose different values of with . That is, two mutualistic species could coexist with delay. However, we find that the feedback delays for two species play almost the same role. The critical value of feedback delay is 0.332 in the system (115) when we only consider it with while the critical value of feedback delay is 0.293 when we only consider it. We also find that feedback delay could destroy the stability of solutions of system (115). That is, the solutions of system (115) turn to oscillate as increasing monotonously from zero. In addition, the solutions of system (115) are stable at first, then oscillate, at last turn to be periodic when time delay of benefit is greater than its critical value and feedback delay belongs to its stable interval (see the right figure of Figure 7), which is different from the result of the case of considering feedback time delay greater than its critical value and feedback delay belonging to its stable interval (see the right figure of Figure 5). This is maybe because time delay of benefit could make the two species stable. This is very valuable in the view of ecology.

It is definitely an interesting future work to consider the following mutual system with different harvesting: where (>0) are the catchability coefficient of two species. The catch rate functions and are based on CPUE (catch per unit effort) hypothesis [25]. Analyzing the complicated bifurcations, we leave this work in the future.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was partially supported by the NNSF of China (10961018), the NSF of Gansu Province of China (1212RJYA011, 2011GS04318), the NSF for Distinguished Young Scholars of Gansu Province of China (1111RJDA003), the Development Program for HongLiu Distinguished Young Scholars in Lanzhou University of Technology, and the Development Program for Outstanding Young Teachers in Lanzhou University of Technology (Q201308).