Abstract and Applied Analysis

Volume 2013 (2013), Article ID 321930, 15 pages

http://dx.doi.org/10.1155/2013/321930

Research Article

Stability and Global Hopf Bifurcation Analysis on a Ratio-Dependent Predator-Prey Model with Two Time Delays

Huitao Zhao,1,2 Yiping Lin,¹ and Yunxian Dai¹

¹Department of Applied Mathematics, Kunming University of Science and Technology, Kunming, Yunnan 650093, China

²Department of Mathematics and Information Science, Zhoukou Normal University, Zhoukou, Henan 466001, China

Received 4 September 2013; Revised 3 November 2013; Accepted 5 November 2013

Academic Editor: István Györi

Copyright © 2013 Huitao Zhao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

A ratio-dependent predator-prey model with two time delays is studied. By means of an iteration technique, sufficient conditions are obtained for the global attractiveness of the positive equilibrium. By comparison arguments, the global stability of the semitrivial equilibrium is addressed. By using the theory of functional equation and Hopf bifurcation, the conditions on which positive equilibrium exists and the quality of Hopf bifurcation are given. Using a global Hopf bifurcation result of Wu (1998) for functional differential equations, the global existence of the periodic solutions is obtained. Finally, an example for numerical simulations is also included.

1. Introduction

The main purpose of this paper is to investigate the bifurcation phenomena from the delays for the following predator-prey system:

where and stand for the population (or density) of the prey and the predator at time , respectively. From the biological sense, we assume that . , and are positive constants, in which denotes the intrinsic growth rate of the prey, is the intraspecific competition rate of the prey, is the capturing rate of the predator, describes the efficiency of the predator in converting consumed prey into predator offspring, is the interference coefficient of the predators, and is the predator mortality rate. The delay denotes the gestation period of the predator; is the hunting delay of the predator to prey.

This model is labeled “ratio-dependent,” which means that the functional and numerical responses depend on the densities of both prey and predators, especially when predator has to search for food. Such a functional response is called a ratio-dependent response function (see [1] for more details). In system (1), the ratio-dependent response function is of the form .

The ratio-dependent predator-prey model has been studied by several researchers recently and very rich dynamics have been observed [2–5]. For example, Xu et al. [4] studied a delayed ratio-dependent predator-prey model with the same ratio-dependent response function of system (1). By means of an iteration technique, they obtained the sufficient conditions for the global attractiveness of the positive equilibrium. By comparison arguments, they proved the global stability of the semitrivial equilibrium. Finally using the theory of functional equation and Hopf bifurcation, they gave the condition on which positive equilibrium exists and the formulae to determine the quality of Hopf bifurcation. But in their work, the global continuation of local Hopf bifurcation was not mentioned.

In general, periodic solutions through the Hopf bifurcation in delay differential equations are local for the values of parameters which are only in a small neighborhood of the critical values (see, e.g., [6, 7]). Therefore we would like to know if these nonconstant periodic solutions obtained through local bifurcation can continue for a large range of parameter values. Recently, a great deal of research has been devoted to the topics [8–12]. One of the methods used in them is the global Hopf bifurcation theorem by Wu [13]. For example, Song et al. [12] studied a predator-prey system with two delays, and using the methods in [13], they get the global existence of periodic solutions.

Motivated by [12], we will study the system (1); special attention is paid to the global continuation of local Hopf bifurcation. We suppose that the initial condition for system (1) takes the form

where , which is the Banach space of continuous functions mapping the interval into , where .

By the fundamental theory of functional differential equations [14], system (1) has a unique solution satisfying initial condition (2).

The rest of the paper is organized as follows. In Section 2, we show the positivity and the boundedness of solutions of system (1) with initial condition (2). In Section 3, we study the existence of Hopf bifurcation for system (1) at the positive equilibrium. In Section 4, using the normal form theory and the center manifold reduction, explicit formulae are derived to determine the direction of bifurcation and the stability and other properties of bifurcating periodic solutions. In Section 5, by means of an iteration technique, sufficient conditions are obtained for the global attractiveness of the positive equilibrium. In Section 6, we consider the global existence of bifurcating periodic solutions and give some numerical simulations. In Section 7, a brief discussion is given.

2. Positivity and Boundedness

In this section, we study the positivity and boundedness of solutions of system (1) with initial conditions (2).

Theorem 1. Solutions of system (1) with initial condition (2) are positive for all .

Proof. Assume to be a solution of system (1) with initial condition (2). Let us consider for . It follows from the second equation of system (1) that

then, from initial condition (2), we have , for . We derive from the first equation of system (1) that

that is, for . This ends the proof.

For the following discussion of boundedness, we first consider the following ordinary differential equation:

where , and are positive constants. From Lemma 2.1 in , it is easy to verify the following result.

Lemma 2. If , the trivial equilibrium of (5) is globally stable. If , then (5) admits a unique positive equilibrium which is globally asymptotically stable in .

Theorem 3. Positive solutions of system (1) with initial condition (2) are ultimately bounded.

Proof. Let be a positive solution of system (1) with initial condition (2). From the first equation of system (1), we have

which yields

hence, for sufficiently small, there is a such that if , .

We now consider the boundedness of . If , we derive from the second equation of system (1) that

from monotone bounded theorem, it is easy to show that .

Therefore, we assume below that . We derive from the second equation of system (1) that, for ,

noting that , by Lemma 2, a comparison argument shows that

This completes the proof.

3. Local Stability and Hopf Bifurcation

In this section, we discuss the local stability of the positive equilibrium and the semitrivial equilibrium of system (1) and establish the existence of Hopf bifurcation at the positive equilibrium.

It is easy to show that system (1) always has a semitrivial equilibrium . Further, if the following condition holds:

 (H1) ,

then system (1) has a unique positive equilibrium , where

where

For convenience, let us introduce new variables , rewriting as , so that system (1) can be written as the following system with a single delay:

Clearly, system (13) has the same equilibrium as system (1).

The characteristic equation of system (13) at the semitrivial equilibrium is of the form

Clearly, (14) always has a root , and if , the other root of (14) is negative; if , the other root of (14) is positive. Hence the semitrivial equilibrium is locally asymptotically stable (unstable) if ().

The characteristic equation of system (13) at the positive equilibrium is of the form

where

where is defined as (12).

When , (15) becomes

It is easy to show that

Obviously, if (H1) holds, then . Hence, the positive equilibrium of system (13) is locally stable when if

and it is unstable when if

We assume that is a root of (15); this is the case if and only if satisfies the following equation:

Separating the real and imaginary parts, we obtain the following system for :

It follows that

Letting , (42) becomes

By a direct calculation, it follows that

Note that if (H1) holds, then . Hence if (H1) and hold, (24) has no positive roots. Accordingly, if (H1) and hold, the positive equilibrium of system (13) exists and is locally asymptotically stable for all . If (H1) and hold, then (24) has a unique positive root , where

Then, we can get

at which (15) admits a pair of purely imaginary roots of the form .

Let and be defined above. Denote

the root of (15) satisfying

It is not difficult to verify that the following result holds.

Lemma 4. If (H1) and hold, the transversal condition holds.

Proof. Differentiating (15) with respect , we obtain that

it follows that

from (15) and (31), we have

We therefore derive that

Noting that , hence, if (H1) and hold, we have . Accordingly, the transversal condition holds and a Hopf bifurcation occurs at .

By Lemma B in [5], we have the following results.

Theorem 5. Suppose (H1) holds and let be defined in (12), for system (13), one has the following.

(i)If and , then the positive equilibrium is locally asymptotically stable for all .(ii)If and , then there exists a positive number such that the positive equilibrium is locally asymptotically stable if and is unstable if . Further, system (13) undergoes a Hopf bifurcation at when .

4. Direction and Stability of Hopf Bifurcations

In Section 3, we have shown that system (13) admits a periodic solution bifurcated from the positive equilibrium at the critical value . In this section, we derive explicit formulae to determine the direction of Hopf bifurcations and stability of periodic solutions bifurcated from the positive equilibrium at critical value by using the normal form theory and the center manifold reduction (see, e.g., [15, 16]).

Set ; then is a Hopf bifurcation value of system (13). Thus we can consider the problem above in the phase space .

Let  . System (13) is transformed into

where

For the simplicity of notations, we rewrite (34) as

where , is defined by , and , are given, respectively, by

By the Riesz representation theorem, there exists a function of bounded variation for such that

In fact, we can choose

where is the Dirac delta function. For , define

Then when , system (36) is equivalent to

where for .

For , define

and a bilinear inner product,

where and denotes the conjugate complex of . Then and are adjoint operators. By the discussion in Section 3, we know that are eigenvalues of . Thus, they are also eigenvalues of . We first need to compute the eigenvector of and corresponding to and , respectively.

Suppose that is the eigenvector of corresponding to . Then . From the definition of , it is easy to get .

Similarly, let be the eigenvector of corresponding to . By the definition of , we can compute .

In order to assure , we need to determine the value of . From (44) and the definitions of and , we have such that and .

In the following, we first compute the coordinates to describe the center manifold at . Define

On the center manifold , we have

where and are local coordinates for center manifold in the directions of and . Note that is real if is real. We consider only real solutions. For the solution , since , we have

where

By (45), we have

It follows from (38) and (48) that

In order to assure the value of , we need to compute and . By (42) and (45), we have

where

Notice that near the origin on the center manifold , we have

thus, we have

By (51), for , we have

Comparing the coefficients with (51) gives that

From (56), (54), and the definition of , we can get

Notice that ; we have

where is a constant vector. In the same way, we can also obtain

where is also a constant vector. In what follows, we will compute and . From the definition of and (54), we have

where .

From (51), (58), and (60) and noting that

we have

where

From (52), (59), and (61) and noting that

we have

where

Thus, we can determine and from (58) and (59). Furthermore, we can determine each . Therefore, each is determined by the parameters and delay in (13). Thus, we can compute the following values [15]:

which determine the quantities of bifurcating periodic solutions in the center manifold at the critical value ; that is, determines the directions of the Hopf bifurcation: if , then the Hopf bifurcation is supercritical (subcritical) and the bifurcation exists for ; determines the stability of the bifurcation periodic solutions: the bifurcating periodic solutions are stable (unstable) if ; and determines the period of the bifurcating periodic solutions: the period increases (decreases) if .

5. Global Attractiveness

In this section, following Chaplygin [17], taking into account the upper and lower solution technique and using monotone iterative methods [18, 19], we discuss the global attractiveness of the positive equilibrium and the global stability of the semitrivial equilibrium of system (1), respectively.

Theorem 6. Suppose (H1) holds and let be defined above, then the positive equilibrium of system (1) is globally attractive provided that the following holds:

 (H2) ,

Proof. Let be any positive solution of system (1) with initial conditions (2).

Let

Using iteration method, we will proof that .

From the first equation of system (1), we have

by comparison, it follows that

hence, for sufficiently small, there exists a such that if .

From the second equation of system (1), we have, for ,

Consider the following auxiliary equation:

Since (H1) holds, by Lemma 2, it follows from (73) that

where is defined in (12). By comparison, we obtain that

since this inequality holds true for arbitrary sufficiently small, it follows that , where

Hence, for sufficiently small, there is a such that if .

For sufficiently small, noting that , we derive from the first equation of system (1) that, for ,

by comparison, it follows that

hence, for sufficiently small, there is a , such that if , .

For sufficiently small, we derive from the second equation of system (1) that, for ,

Consider the following auxiliary equation:

Since (H1) holds, by Lemma (5), it follows from (80) that

by comparison we derive that

Since this inequality holds true for arbitrary sufficiently small, we conclude that , where

Therefore, for sufficiently small, there is a such that if , .

Again, for sufficiently small, it follows from the first equation of system (1) that, for ,

by comparison we derive that

Since the above inequality holds true for arbitrary sufficiently small, it follows that , where

hence, for sufficiently small, there is a such that if , .

It follows from the second equation of system (1) that, for ,

By Lemma 2 and a comparison argument we derive from (87) that

since this inequality holds true for sufficiently small, we get , where

hence, for sufficiently small, there is a such that if , .

For sufficiently small, it follows from the first equation of system (1) that, for ,

by comparison, we can obtain that

Since the above inequality holds true for arbitrary sufficiently small, it follows that , where

therefore, for sufficiently small, there is a such that if , .

For sufficiently small, we derive from the second equation of system (1) that, for ,

Since (H1) holds, by Lemma 2 and a comparison argument, it follows (93) that

since, for arbitrary sufficiently small, this inequality holds true, we conclude that , where

Continuing this process, we obtain four sequences such that, for ,

where is defined in (12). It is readily seen that

It is easy to know that the sequences are not increasing and the sequences are not decreasing; from accumulation point theorem, the limit of each sequence in , and exists, Denote

We therefore obtain from (96) and (98) that

To complete the proof, it is sufficient to prove that . It follows from (99) that

Letting (100) minus (101), we have

If , we derive from (102) that

Letting , we derive from (103) that

It follows from (104) that

noting that , we derive from (105) that

Substituting (78) into (106), it follows that

Hence, if (H2) holds, we have ; this is a contradiction. Accordingly, we have . Therefore, from (99), we have . Hence, the positive equilibrium is globally attractive. The proof is complete.

Using the same methods in [4, 20], we can also get a similar result.

Theorem 7. If and , the semitrivial equilibrium of system (1) is globally asymptotically stable.

6. Global Continuation of Local Hopf Bifurcations

In this section, we study the global continuation of periodic solutions bifurcating from the positive equilibrium of system (13). Throughout this section, we follow closely the notations in [13]. For simplification of notations, setting , we may rewrite system (13) as the following functional differential equation:

where . It is obvious that if (H1) holds, then system (13) has a semitrivial equilibrium and a positive equilibrium . Following the work of [13], we need to define

Let denote the connected component passing through in , where is defined by (26). From Theorem 5, we know that is nonempty.

We first state the global Hopf bifurcation theory due to Wu [13] for functional differential equations.

Lemma 8. Assume that is an isolated center satisfying the hypotheses ()–() in [13]. Denote by the connected component of in . Then either

(i) is unbounded or(ii) is bounded; is finite and

for all , where is the crossing number of if or it is zero if otherwise.

Clearly, if (ii) in Lemma 8 is not true, then is unbounded. Thus, if the projections of onto -space and onto -space are bounded, then the projection onto -space is unbounded. Further, if we can show that the projection of onto -space is away from zero, then the projection of onto -space must include interval . Following this ideal, we can prove our results on the global continuation of local Hopf bifurcation.

Lemma 9. If condition (H1) holds, then all nonconstant periodic solutions of (13) with initial conditions,

are uniformly bounded.

Proof. Suppose that are nonconstant periodic solutions of system (13) and define

It follows from system (13) that

which implies that the solutions of system (13) cannot cross the -axis and -axis. Thus the nonconstant periodic orbits must be located in the interior of each quadrant. It follows from initial conditions of system (13) that . From system (13), we can get

Since , it follows from the first equation of (114) that

on the other hand, by the second equation of (114) and (115), we have

where is defined in (12). From the discussion above, the lemma follows immediately.

Lemma 10. If conditions (H1) and (H2) hold, then system (13) has no nonconstant periodic solution with period .

Proof. Suppose for a contradiction that system (13) has nonconstant periodic solution with period . Then the following system (117) of ordinary differential equations has nonconstant periodic solution:

which has the same equilibria as system (13), that is, and a positive equilibrium . Note that -axis and -axis are the invariable manifold of system (13) and the orbits of system (13) do not intersect each other. Thus, there is no solution crossing the coordinate axis. On the other hand, note the fact that if system (117) has a periodic solution, then there must be the equilibrium in its interior and are located on the coordinate axis. Thus, we conclude that the periodic orbit of system (117) must lie in the first quadrant. From the proof of Theorem 6, we known that if (H1) and (H2) hold, the positive equilibrium is asymptotically stable and globally attractive; thus, there is no periodic orbit in the first quadrant. This ends the proof.

Theorem 11. Suppose the conditions (H1) and (H2) hold; let and be defined in (26). If , then system (13) has at least periodic solutions for every .

Proof. It is sufficient to prove that the projection of onto -space is for each , where .

The characteristic matrix of (108) at an equilibrium takes the following form:

is called a center if and . A center is said to be isolated if it is the only center in some neighborhood of . It follows from (118) that

where , , and are defined as in Section 3. From the discussion in Section 3, each of (119) and (120) has no purely imaginary root provided that . Thus, we conclude that (108) has no the center of the form as and . On the other hand, from the discussion in Section 3 about the local Hopf bifurcation, it is easy to verify that is an isolated center, and there exist  , , and a smooth curve such that for all and

Let

It is easy to verify that, on ,

Therefore, the hypotheses ()–() in [13] are satisfied. Moreover, if we define

then we have the crossing number of isolated center as follows:

Thus, we have

where has all or parts of the form . It follows from Lemma 8 that the connected component through in is unbounded. From (26), we can know that if (H1) holds, for ,

Now we prove that the projection of onto -space is , where . Clearly, it follows from the proof of Lemma 10 that system (13) with has no nontrivial periodic solution. Hence, the projection of onto -space is away from zero.

For a contradiction, we suppose that the projection of onto -space is bounded; this means that the projection of onto -space is included in an interval . Noticing and applying Lemma 10 we have for belonging to . Applying Lemma 9, we know that the projection of onto -space is bounded. So the component of is bounded. This contradicts our conclusion that is unbounded. The contradiction implies that the projection of onto -space is unbounded above.

Hence, system (13) has at least periodic solution for every . This completes the proof.

Example 12. In system (1), we first choose , and . As depicted in Figure 1, a bifurcation diagram is given for system (1) with respect to the parameters and . By the discussion in Section 3, system (1) always has a semitrivial equilibrium , and if , is asymptotically stable; otherwise, is unstable. So if we choose , as depicted in Figure 1, is always unstable. In domains II, V, and VI, the positive equilibrium is not feasible. In domains I, III, and IV, system (1) has a unique positive equilibrium; it is locally asymptotically stable in domain I and is unstable in domain IV. In domain III, system (1) undergoes a Hopf bifurcation at the positive equilibrium at some . Further, we choose , ,. In this case, system (1) has a positive equilibrium . By computation, we have , and . From Theorem 5, is stable when as illustrated by numerical simulations (see Figure 2). When passes through the critical value , the equilibrium loses its stability and a Hopf bifurcation occurs; that is, a family of periodic solution bifurcates from . By the algorithm derived in Section 3 and Section 4, we have , which implies that , and . Thus, by the discussion in Section 4, the Hopf bifurcation is supercritical for , the bifurcating periodic solutions from at are asymptotically stable, and the period of these periodic solutions is increasing with the increasing of , which are depicted in Figures 3, 4, and 5. Furthermore, Figure 5 shows that the local Hopf bifurcation implies the global Hopf bifurcation after the second critical value of .

Figure 1: The bifurcation diagram of system (1) with , and , where , and .

(a)

(b)

(c)

Figure 2: The trajectories and phase graphs of system (1) with .

(a)

(b)

(c)

Figure 3: The trajectories and phase graphs of system (1) with .

(a)

(b)

(c)

Figure 4: The trajectories and phase graphs of system (1) with .

(a)

(b)

(c)

Figure 5: The trajectories and phase graphs of system (1) with .

7. Discussion

In this paper, we have studied a ratio-dependent predator-prey model with two time delays. By analyzing the corresponding characteristic equation, the local stability of the positive equilibrium and the semitrivial equilibrium of system (1) was discussed. We have obtained the estimated length of gestation delay which would not affect the stable coexistence of both prey and predator species at their equilibrium values. The existence of Hopf bifurcation for system (1) at the positive equilibrium was also established. From theoretical analysis it was shown that the larger values of gestation time delay cause fluctuation in individual population density and hence the system becomes unstable. As the estimated length of delay to preserve stability and the critical length of time delay for Hopf bifurcation are dependent upon the parameters of system, it is possible to impose some control, which will prevent the possible abnormal oscillation in population density. The global attractiveness result in Theorem 6 implied that system (1) is permanent if the intrinsic growth rate of the prey and the conversion rate and the interference rate of the predator are high, and the death rate of the predator is low. From Theorem 7 we see that if the death rate of the predator is greater than the conversion rate of the predator, the predator population become extinct for any gestation delay. In particular, the results about boundedness and attractiveness are similar to the results of [4]. From the discussion in Sections 3 and 4, we see that if the values of , and are given, we can get the Hopf bifurcation value of , and further we may determine the direction of Hopf bifurcation and the stability of periodic solutions bifurcating from the positive equilibrium at the critical point . Furthermore, we show that the local Hopf bifurcation implies the global Hopf bifurcation after the second critical value of delay.

Acknowledgment

The authors would like to thank the anonymous referee for very helpful suggestions and comments which led to improvements of their original paper. This work is supported by the National Natural Science Foundation of China (no. 11061016), Science and Technology Department of Henan Province (no. 122300410417) and Education Department of Henan Province (no. 13A110108).

References

R. Arditi and L. R. Ginzburg, “Coupling in predator-prey dynamics: ratio-dependence,” Journal of Theoretical Biology, vol. 139, no. 3, pp. 311–326, 1989.

D. Xiao, W. Li, and M. Han, “Dynamics in a ratio-dependent predator-prey model with predator harvesting,” Journal of Mathematical Analysis and Applications, vol. 324, no. 1, pp. 14–29, 2006.

M. Xiao and J. Cao, “Hopf bifurcation and non-hyperbolic equilibrium in a ratio-dependent predator-prey model with linear harvesting rate: analysis and computation,” Mathematical and Computer Modelling, vol. 50, no. 3-4, pp. 360–379, 2009.

R. Xu, Q. Gan, and Z. Ma, “Stability and bifurcation analysis on a ratio-dependent predator-prey model with time delay,” Journal of Computational and Applied Mathematics, vol. 230, no. 1, pp. 187–203, 2009.

Y. Kuang and E. Beretta, “Global qualitative analysis of a ratio-dependent predator-prey system,” Journal of Mathematical Biology, vol. 36, no. 4, pp. 389–406, 1998.

S. Ruan and J. Wei, “Periodic solutions of planar systems with two delays,” Proceedings of the Royal Society of Edinburgh A, vol. 129, no. 5, pp. 1017–1032, 1999.

S. Ruan and J. Wei, “On the zeros of transcendental functions with applications to stability of delay differential equations with two delays,” Dynamics of Continuous, Discrete & Impulsive Systems A, vol. 10, no. 6, pp. 863–874, 2003.

B. Zheng, Y. Zhang, and C. Zhang, “Global existence of periodic solutions on a simplified BAM neural network model with delays,” Chaos, Solitons and Fractals, vol. 37, no. 5, pp. 1397–1408, 2008.

C. Sun and M. Han, “Global Hopf bifurcation analysis on a BAM neural network with delays,” Mathematical and Computer Modelling, vol. 45, no. 1-2, pp. 61–67, 2007.

F. Botelho and V. A. Gaiko, “Global analysis of planar neural networks,” Nonlinear Analysis. Theory, Methods & Applications, vol. 64, no. 5, pp. 1002–1011, 2006.

X.-Y. Meng, H.-F. Huo, and X.-B. Zhang, “Stability and global Hopf bifurcation in a delayed food web consisting of a prey and two predators,” Communications in Nonlinear Science and Numerical Simulation, vol. 16, no. 11, pp. 4335–4348, 2011.

Y. Song, Y. Peng, and J. Wei, “Bifurcations for a predator-prey system with two delays,” Journal of Mathematical Analysis and Applications, vol. 337, no. 1, pp. 466–479, 2008.

J. Wu, “Symmetric functional-differential equations and neural networks with memory,” Transactions of the American Mathematical Society, vol. 350, no. 12, pp. 4799–4838, 1998.

J. K. Hale and S. M. Verduyn Lunel, Introduction to Functional-Differential Equations, vol. 99 of Applied Mathematical Sciences, Springer, New York, NY, USA, 1993.

B. D. Hassard, N. D. Kazarinoff, and Y. H. Wan, Theory and Applications of Hopf Bifurcation, vol. 41 of London Mathematical Society Lecture Note Series, Cambridge University Press, Cambridge, Mass, USA, 1981.

T. Faria and L. T. Magalhães, “Normal forms for retarded functional-differential equations and applications to Bogdanov-Takens singularity,” Journal of Differential Equations, vol. 122, no. 2, pp. 201–224, 1995.

S. A. Chaplygin, Complete Works III, Akademia Nauk SSSR, Leningrad, Russia, 1935.

N. S. Kurpel and V. I. Grechko, “On some modifications of Chaplygin's method for equations in partially orderedspaces,” Ukrainian Mathematical Journal, vol. 25, no. 1, pp. 30–36, 1973.

A. Dzieliński, “Stability of a class of adaptive nonlinear systems,” International Journal of Applied Mathematics and Computer Science, vol. 15, no. 4, pp. 455–462, 2005.

X. Song and L. Chen, “Optimal harvesting and stability for a two-species competitive system with stage structure,” Mathematical Biosciences, vol. 170, no. 2, pp. 173–186, 2001.