Abstract

We propose a prey predator model with stage structure for prey. A discrete delay and a distributed delay for predator described by an integral with a strong delay kernel are also considered. Existence of two feasible boundary equilibria and a unique interior equilibrium are analytically investigated. By analyzing associated characteristic equation, local stability analysis of boundary equilibrium and interior equilibrium is discussed, respectively. It reveals that interior equilibrium is locally stable when discrete delay is less than a critical value. According to Hopf bifurcation theorem for functional differential equations, it can be found that model undergoes Hopf bifurcation around the interior equilibrium when local stability switch occurs and corresponding stable limit cycle is observed. Furthermore, directions of Hopf bifurcation and stability of the bifurcating periodic solutions are studied based on normal form theory and center manifold theorem. Numerical simulations are carried out to show consistency with theoretical analysis.

1. Introduction

In recent years, much research efforts have been extensively made on interaction and coexistence mechanism of population in prey predator ecosystem by means of Lotka-Volterra dynamical models [13]. Generally, in order to reflect the dynamical behavior of mathematical models depending on the past history, time delay is usually incorporated into model, which can be utilized to mathematically describe hunting delay, maturation delay, and gestation delay for population within prey predator ecosystem. With the introduction of time delay, it may cause the loss of stability and other complicated dynamical behavior of model, such as Hopf bifurcation and saddle node bifurcation. It should be noted that there is a well-developed theory of dynamical models which incorporates time delay into model [4]. In particular, the properties of periodic solutions arising from the Hopf bifurcation are of great interest [4, 5].

Recently, plenty of dynamical models with discrete and distributed delays have been proposed to discuss the population dynamics of prey predator ecosystem due to variation of maturation and hunting factors [613]. Song and Peng [9] proposed a logistic model with discrete and distributed delays where , and are positive constants. Function is called the delayed kernel, which is the weight given to the population at time . Under the assumption that for all and the normalized condition that , which ensures that the steady state of model (1) is unaffected by the delay, the local stability of model (1) around interior equilibrium and existence of Hopf bifurcation are studied. Furthermore, direction of Hopf bifurcations and the stability of bifurcated periodic solutions are investigated by using the normal form theory and center manifold theorem.

By supposing that the predator population at every age stage has the predation ability and the prey population captured by the predator population in the past is all contributing to the predator population at time , the growth dynamics of the two species can be described by the following delayed Lotka-Volterra two species prey predator model with distributed delays: where constants and denote intrinsic growth rate of prey population and death rate of predator population, respectively. , , and represent the carrying capacity of prey population, predator coefficient of predator population, biomass conversion rate of prey population captured by predator population, and intraspecific competition rate of predator population, respectively. It should be noted that for , , , the delay kernel and are bounded nonnegative functions and the following normalized conditions hold:

Model system (2) with various delay kernels and delayed intraspecific competitions have been extensively investigated in [613]. Faria [6] investigated the stability of interior equilibrium of model (2) and Hopf bifurcation of nonconstant periodic solutions around the interior equilibrium. When and , , , model (2) admits two different discrete time delays; Ruan [7] and Yan and Zhang [11] discussed the stability of the interior equilibrium of model (2) and Hopf bifurcation of nonconstant periodic solutions regarding the sum of two delays and as the bifurcation parameter. Furthermore, dynamic effect of intraspecific competition on population dynamics of model (2) is studied in [8, 9]. By assuming that , , the delay kernel function may take the weak generic kernel function and strong generic kernel function (), where the weak generic kernel implies that the importance of events in the past simply decreases exponentially and the further one looks into the past while the strong generic kernel implies that a particular time in the past is more important than any other [1, 2, 5]. Under this assumption, model (2) can be reduced to the following system with a discrete delay and a distributed delay:

When the delay kernel function admits the weak generic kernel, Song and Yuan [10] and Ma et al. [12] discussed the local asymptotical stability of interior equilibrium and Hopf bifurcations of nonconstant periodic solutions based on linearization method and regarding the discrete hunting delay as bifurcation parameter. It is shown that interior equilibrium is asymptotically stable when is less than a certain critical value and Hopf bifurcation occurs at a critical value of the discrete hunting delay. They also studied the direction of the Hopf bifurcations and the stability of bifurcated periodic solutions occurring through Hopf bifurcations. In the case of strong generic kernel, Zhang et al. [13] studied local stability of all boundary and interior equilibria; conditions for existence of Hopf bifurcation are also investigated. Furthermore, an explicit algorithm determining the direction of Hopf bifurcations and stability of bifurcating periodic solutions occurring through Hopf bifurcations is also given.

In the natural world, many species have a life history that takes them through two stages, juvenile stage and adult stage. Individuals in each stage are identical in biological characteristics, and some vital rates (rates of survival, development, and reproduction) of individuals in a population almost always depend on stage structure. Furthermore, many complex biological phenomena arising in prey predator ecosystem always depend on the past history of the system, and it has been recognized that time delay may have complicated impact on the dynamics of prey predator ecosystem [2]. In the past several decades, there has been an increasing interest in prey predator model with stage structure and time delay. In the model proposed by Aiello and Freedman [14], stage structure of single population growth with stage structure and time delay representing for maturation of population is considered. Their model predicts a positive steady state as the global attractor, thereby suggesting that stage structure does not generate sustained oscillations frequently observed in single population in the real world. Subsequent work made by other authors [1520] suggest that time delay to adulthood should be state dependent. Al-Omari and Gourley [16] suggested that the time delay to adulthood should be state dependent and careful formulation of such state dependent time delays can lead to models that produce periodic solutions. Xu et al. [21] studied the persistence and stability of a delayed prey predator model with stage structure for predator. Gourley and Kuang [22] and Bandyopadhyay and Banerjee [23] formulated a class of general and robust prey predator models with stage structure and constant maturation time delay and performed a systematic mathematical and computational study. They have shown that there is a window in maturation time delay parameter that generates sustainable oscillatory dynamics. Ration dependent prey predator model is proposed in the work done in [12, 2428] and permanence and stability analysis are also investigated.

In model (4) with strong generic kernel proposed in [13], stage structure of prey population is not considered. It is well known that the immature prey population is usually considered as an easier target for predators compared with mature prey population [14]. Hence, the mature prey population predated by predator population can be ignored [1], and it is necessary to investigate dynamic effect of stage structure on prey predator model with discrete delay and distributed delay. By incorporating stage structure of prey population into model (4) with strong generic kernel, work done in [13] is extended in this paper. The organization of rest sections of this paper is as follows: a stage-structured prey predator model with discrete and distributed delay is proposed in the second section. It should be noted that the delay kernel function takes the strong generic kernel function and stage structure for prey population is considered in this paper. In the third section, qualitative analysis of the proposed model is investigated. Conditions for existence of two feasible boundary equilibria and a unique interior equilibrium are analytically discussed. By analyzing corresponding characteristic equation, local stability of two feasible boundary equilibria and an interior equilibrium are discussed, respectively. By taking discrete time delay as a parameter, Hopf bifurcation around the interior equilibrium is studied due to variation of discrete delay. Furthermore, directions of Hopf bifurcation and stability of the bifurcating periodic solutions are studied based on normal form theory and center manifold theorem. In the fourth section, numerical simulations are carried out to show consistency with theoretical results obtained in this paper. Finally, this paper ends with a conclusion.

2. Model Formulation

In this paper, a stage-structured prey predator model with discrete and distributed delay is investigated based on the following four hypotheses, which are given as follows.(H1)Prey population is divided into two-stage groups, that is, immature prey population and mature prey population . The birth of immature prey population is proportional to the existing mature prey population with proportionality constant at any time , and therefore the term is in the first equation of model (5). and represent the intraspecific competition rate of prey population and predation coefficient of predator population, respectively.(H2)The rate of transformation of the mature prey population is proportional to the existing immature prey population with proportionality constant , which explains the term in the first equation and in the second equation. denotes death rate of mature prey population.(H3) denotes death rate of predator population. It is assumed that the immature prey population is usually considered as an easier target for predators; the mature prey population predated by predator population can be ignored. The predator is considered as a whole group, and only immature prey population is under predation from predator population. and denote biomass conversion rate of immature prey population captured by predator population and intraspecific competition rate of predator population, respectively.(H4)The delay kernel function takes strong generic kernel function (); the strong generic kernel implies that a particular time in the past is more important than any other [1, 2, 5].

Based on hypotheses (H1)–(H4) and model (4), a stage-structured prey predator model with discrete hunting delay and distributed maturation delay is established as follows: where , , and denote the density of immature prey, mature prey, and predator population, respectively. Other parameters share the same interpretations introduced in hypotheses (H1)–(H4).

3. Qualitative Analysis of Model System

In this section, conditions for existence of two feasible boundary equilibria and a unique interior equilibrium are analytically investigated. By analyzing corresponding characteristic equation, local stability of model around boundary equilibrium and interior equilibrium is discussed, respectively. By taking discrete time delay as parameter, conditions for existence of Hopf bifurcation are discussed based on Hopf bifurcation theorem for functional differential equations. Furthermore, directions of Hopf bifurcation and stability of the bifurcating periodic solutions are studied by using normal form theory and center manifold theorem.

3.1. Stability Analysis of Equilibria

For model (5), it follows from normalized condition (3) that there are two feasible boundary equilibria and a unique interior equilibrium, which are given as follows:(i)two feasible boundary equilibria: , ;(ii)a unique interior equilibrium: exists provided that where = , = , and = .

Theorem 1. The local stability of model (5) around two-boundary equilibria is as follows:(i) is a saddle node of model (5);(ii) is a saddle node of model (5) provided that exists;(iii) is locally stable if and .

Proof. By simple computation, the characteristic equation of model (5) around takes the following form:
By solving (7), three eigenvalues of model (5) around can be obtained as follows: , , and . According to introduced in (H2), it derives that , , and . Hence, is a saddle node of model (5).
By simple computation, the characteristic equation of model (5) around takes the following form:
By solving (8), three eigenvalues of model (5) around are as follows: , which derives that if exists; that is, holds; if .
Furthermore, the other two eigenvalues of (8) can be determined by the following equation:
It follows from Routh-Hurwitz criterion [29] that , provided that .
Consequently, is a saddle node of model (5) provided that exists; is locally stable if and .

Define new variables and as follows:

According to the law of solving the derivative for an integral with parameterized variables, it can be obtained that

By virtue of (11), model (5) can be transformed into the following five-dimensional system of FDEs with a discrete delay:

It follows from model (12) that for the interior equilibrium of model (12), which will be denoted as in the following part of this paper. Furthermore, the interior equilibrium of model (5) is transformed into the interior equilibrium of model (12). Hence, the qualitative local stability analysis around of model (12) is equivalent to the qualitative local stability analysis around interior equilibrium of model (5).

According to the leading matrix of model (12) around the interior equilibrium , the characteristic equation takes the following form: which derives that where

When , (14) takes the following form:

Based on the values of , , , , , , , , defined in (14), it is easy to show that and .

Define quantities , as follows:

Further computations show that for . According to the above analysis and Routh-Hurwitz criteria [2], it follows from , , and for that model (12) is locally stable around in the case of , which implies that model (5) is locally stable around .

Subsequently, the dynamic effect of discrete delay on model dynamics is investigated. In the case of positive discrete delay , let () be a root of (14). Separating the real and imaginary parts of (14) with yields the following equations:

By adding up the squares of the corresponding sides of (18), it yields the following algebraic equation with respect to : where , , , , and .

Let ; (19) can be rewritten as follows:

By virtue of derivative of with respect to , it gives that

By simple computation, it is easy to show that , , , , , which derives that on and is strictly monotonically increasing on . Consequently, has a unique positive root if and does not have a positive root if .

Suppose , which derives that and then there exists a unique root of , which is denoted by , and the unique positive root of (14) is .

Suppose , which derives that and then (14) does not have a positive root.

By eliminating , it follows from (18) that and then (14) has a pair of purely imaginary roots when for .

By using Butler’s lemma [29], the interior equilibrium remains stable for . Based on the above analysis, local stability analysis of model (5) around the interior equilibrium can be concluded in the following Theorem 2.

Theorem 2. If condition (6) holds, then the following results hold.(i)If , then model (5) is locally stable around for any values .(ii)If , then model (5) is locally stable around for .
Subsequently, the conditions for Hopf bifurcation in [29] are utilized to investigate whether there is a phenomenon of Hopf bifurcation as increases through .

Theorem 3. If , then model (5) undergoes Hopf bifurcation at the interior equilibrium when , . Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium when and .

Proof. As mentioned above, let represent the purely imaginary root of (14). It follows from (14) that , which determines a set of possible values of .
In the following part, we determine the direction of motion of as is varied; namely, we determine
By differentiating (14) with respect to , it can be obtained that
If , then which signifies that
Hence, the transversality condition [29] holds and model (5) undergoes Hopf bifurcation at the interior equilibrium when , . Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium when and .

3.2. Direction and Stability of Hopf Bifurcation

By using normal theory and center manifold theorem [29, 30], directions of Hopf bifurcation and stability of the bifurcating periodic solutions are discussed in this section. As analyzed in Section 3.1, some transformations associated with the interior equilibrium of model (12) are given as follows: and then is the Hopf bifurcation value of model (12). Bars of variables are dropped for simplicity of notations; model (12) is transformed to a functional differential equation in , where is the Banach space of continuous functions mapping into , , for , and , are defined as follows, respectively:

It is easy to show that is a continuous linear function mapping into . According to Riesz representation theorem [4], there exists a matrix function of bounded variation for such that where .

In fact, we can choose where denotes the Dirac delta function.

If is any given function in and is the unique solution of the linearized equation of (29) with initial function at zero, then the solution operator is defined by

It follows from Lemma  7.1.1 in [29] that , is a strongly continuous semigroup of linear transformation on and the infinitesimal generator of , is as follows: for , the space of functions mapping the interval into which have continuous first derivative, and also define and then model (29) is equivalent to

For , the space of functions mapping interval into the five-dimensional row vectors which have continuous first derivative, define and a bilinear inner product where . It follows from the above analysis that and are adjoint operators.

By virtue of discussion in Section 3.1, are eigenvalues of . Hence, they are also eigenvalues of . In the following, eigenvectors of and correspond to and , respectively.

Suppose is the eigenvectors of corresponding to , which derives that . By using the definition of , (30), and (31), it gives that

For , then it can be obtained that

Similarly, it follows from simple computation that eigenvector of corresponds to ; that is,

In order to assume , we need to determine the value of in the following part.

By virtue of , it derives that where

Hence, we can choose as follows:

Nextly, we will compute the coordinate to describe the centre manifold at . Let be the solution of (37) when .

Define

On the center manifold , it derives that where and are local coordinates for center manifold in the direction of and .

It is noted that is real if is real, and we only consider real solutions. For solution of (37), since , it derives that

The above equation can be rewritten as follows: where

It follows from (46) and (48) that

By virtue of (31), (32), and (33), it derives that

By comparing the coefficients with (51), it gives that

Since is associated with and , further attempts should be carried out to compute and , which can be found in the appendix.

Furthermore, we can compute based on (54). Hence, the following values can be computed as follows:

Some properties of bifurcating periodic solution in the center manifold at the critical values follow from [29]. Based on the analysis in Section 3.1 of this paper, the following theorem can be concluded.

Theorem 4. The properties of Hopf bifurcation are determined by values computed in (55).(i) determines directions of Hopf bifurcation: if   , then Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for .(ii) determines the stability of bifurcating periodic solutions: bifurcating periodic solutions are stable (unstable) if .(iii) determines the period of bifurcating periodic solutions: period increases (decreases) if .

4. Numerical Simulation

In this section, values of parameters are partially taken from numerical simulation of [13] and set in appropriate units, which can be found as follows: , , , , , , , , and . Numerical simulations are provided to support the theoretical findings obtained in Section 3 of this paper.

By virtue of given values of parameters, it can be computed that which implies that (6) holds and there exists a unique interior equilibrium of model (5), that is, . Furthermore, it follows from given values of parameters that which shows that (22) holds and (19) has a unique positive root. Based on (24), the corresponding can be computed. According to Theorem 2, model is locally stable around the interior equilibrium when . Dynamical responses of model (5) with are plotted in Figure 1 and phase portrait of model (5) with is plotted in Figure 2. It should be noted that in Figures 1 and 2 is randomly selected in the interval , which is enough to merit the above mathematical study. It follows from that , , and . As increases through , a periodic solution caused by the phenomenon of Hopf bifurcation occurs, that is, a family of periodic solutions bifurcate from the interior equilibrium based on Theorem 3. Since and , the Hopf bifurcation is supercritical, the directions of the Hopf bifurcation are , and these bifurcating periodic solutions from the interior equilibrium at are stable based on Theorem 4. Dynamical responses of model (5) with are plotted in Figure 3, since , which shows that model (5) is unstable around the interior equilibrium and stable periodic solutions bifurcate from the interior equilibrium . Furthermore, Figure 4 indicates a limit cycle corresponding to the stable periodic solutions plotted in Figure 3.

5. Conclusion

In this paper, a dynamical prey predator model with discrete hunting delay and distributed maturation delay is proposed to investigate the dynamic effect of time delay and stage structure on population dynamics. Many species in the natural world have a life history that takes them through two stages, juvenile stage and adult stage. Individuals in each stage are identical in biological characteristics, and some vital rates (rates of survival, development, and reproduction) of individuals in a population almost always depend on stage structure. In model (4) with strong generic kernel [13], stage structure of prey population is not considered. It is well known that the immature prey population is usually considered as an easier target for predators compared with mature prey population. Hence, the mature prey population predated by predator population can be ignored [1], and it is necessary to investigate the dynamic effect of stage structure on population dynamics of model (4).

By incorporating stage structure for prey population into the model, the work done in [13] is extended in this paper, where discrete delay and distributed delay for predator population described by an integral with a strong delay kernel are also studied in the proposed model. In Section 3 of this paper, conditions for existence of two boundary equilibria and a unique interior equilibrium are discussed. Local stability around two feasible boundary equilibria is analyzed in Theorem 1. By taking the discrete delay as a parameter, local stability analysis around the unique interior equilibrium is investigated due to variation of discrete delay in Theorem 2. It reveals that discrete delay has dynamic effect on population dynamics, which shows that model (5) is locally stable around interior equilibrium for any nonnegative discrete delay provided that and model (5) remains locally stable around interior equilibrium within under the assumption . According to Hopf bifurcation theorem for functional differential equations, it can be found that Hopf bifurcation occurs when the discrete delay crosses through a sequence of critical values and corresponding stable limit cycle is also observed, which can be found in Theorem 3. Further attempts are made to study the directions of Hopf bifurcation and stability of the bifurcating periodic solutions in Theorem 4.

Compared with the previous related work [13], theoretical results obtained in this paper concentrate on stage-structured prey predator model with discrete hunting delay and distributed maturation delay, which are potentially beneficial for administrative agency to have a deep insight on the coexistence and interaction mechanism of prey predator ecosystem in the real world. It makes work done in this paper has some positive and new feature.

Appendix

The detailed mathematical arguments about computation of and in Theorem 4 can be found in the following part.

By virtue of (37) and (46), we have where

By substituting the corresponding series into (A.1) and comparing the coefficients, we have

It follows from (A.1) that for ,

By comparing coefficients in (A.2) with those in (54), it derives that

Based on the definition of and (A.3) and (A.5), it can be obtained that

For , where is a constant vector.

Similarly, it follows from (A.3) and (A.6) that where is a constant vector.

Subsequently, values of and should be computed. By using the definition of and (A.1), we have where . Based on (A.1), it derives that in the case of , which follows that

By virtue of (31), it gives that

By virtue of (46), it can be obtained that and then we have

According to (A.13) and (A.16), we have

Since is the eigenvalue of and is the corresponding eigenvector, we obtain that where is the identity matrix.

By substituting (A.8) and (A.17) into (A.10), it can be obtained that which can be rewritten as follows:

Based on Grammar Law [2], , , , , and can be obtained as follows:where

Similarly, substituting (A.9) and (A.18) into (A.11), it can be obtained that

Based on Grammar Law [2], , , and can be obtained as follows:where

It follows from the above computation and analysis that and can be determined based on (A.8) and (A.9).

Conflict of Interests

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

Acknowledgments

Authors gratefully acknowledge anonymous reviewers and editors valuable comments. Authors also gratefully acknowledge editors consideration and patient work. This work is supported by National Natural Science Foundation of China, Grant no. 61104003, Grant no. 61273008, and Grant no. 61104093, Research Foundation for Doctoral Program of Higher Education of Education Ministry, Grant no. 20110042120016, Hebei Province Natural Science Foundation, Grant no. F2011501023, Fundamental Research Funds for the Central Universities, Grant no. N120423009, and Research Foundation for Science and Technology Pillar Program of Northeastern University at Qinhuangdao, Grant no. XNK201301. This work is supported by State Key Laboratory of Integrated Automation of Process Industry, Northeastern University, and Hong Kong Admission Scheme for Mainland Talents and Professionals, Hong Kong Special Administrative Region.