1 Introduction

The dynamic relationship between predators and prey has long been, and will continue to be, one of the dominant themes in both ecology and mathematical ecology due to its universal existence and importance (Berryman, [1]). In most of ecosystems, the population of one species does not respond instantaneously to interactions with other species. To incorporate this idea in a modeling approach, time delay models have been developed. In most cases, time delays have a destabilizing effect towards dynamical behavior and often time delays are responsible for oscillations of various species. The question of global stability and uniform persistence of individual species involved with the model under consideration is important in a delay differential equation model. There are several publications which explain from mathematical and ecological points of view the necessity of delay differential equation models (Gopalsamy, [2]; Kuang, [3]).

Time delays of one type or another have been incorporated into biological models by many researchers. Freedman and Rao [4] obtained criteria for local stability of predator–prey model with delays. Freedman and Waltman [5] consider a general model of two predators competing for a single prey. They derived criteria for strong persistence in terms of conditions on system parameters. Kuang [6] studied global stability results obtained from comparison analysis, Bendixson–Dulac criterion or limit cycle stability analysis for the general, Gauss-type, predator–prey system without delay. The obtained criteria involve restrictions on the functions (such as prey species growth rate in the absence of predation and predator functional response). Delay models have also been investigated by Hale and Waltman, [7]; Waltman, [8] and Wang and Ma, [9]; for Lotka–Volterra systems. Lu and Takeuchi [10] have proved that a two species Lotka–Volterra delayed competition system is permanent under any delay effect provided that the corresponding undelayed system has a globally stable positive equilibrium. They have also obtained conditions for global stability of positive equilibrium.

Modeling of population ecological interactions involving time delay is being dealt by Kuang [3]. Aziz-Alaoui and Daher Okiye [11], Cao and Freedman [12], Upadhyay and Rai [13], and Upadhyay and Iyenger [14] consider prey predator models and find some significant results. Xiao and Chen [15] consider a system of retarded functional differential equations as a predator prey model with disease in the prey. Permanence and global stability are analyzed. They show that positive equilibrium is locally asymptotically stable when the time delay is suitably small, while a loss of stability by Hopf-bifurcation can occur as the delay increases. Mukherjee and Roy [16] proposed a generalized prey–predator system with time delay and find the conditions for uniform persistence and global stability. Recently, Kar [17] studied a Gaussian-type prey–predator model with selective harvesting and introduced a time delay in the harvesting term. In general, delay differential equations exhibit much more complicated dynamics than ordinary differential equations since the time delay could cause a stable equilibrium to become unstable and cause the population to fluctuate.

In this article, we have considered a two predator-one prey system with time delay due to gestation. The response function is of the Holling type II.

Before we introduce the model and its analysis we would like to present a brief sketch of the construction of the model which may indicate the biological relevance of it:

(i) There are three populations namely, two predators whose population densities are y and z, and one prey whose population density is denoted by x;

(ii) In absence of predation, the prey population grows according to a logistic law of growth with intrinsic growth rate r and carrying capacity K;

(iii) One predator species consumes the prey with the functional response ${\alpha }_{1}xy/(a_1+x)$, known as the Holling-type II functional response and contributes to its growth rate ${\alpha }_1{\beta }_{1}xy/(a_1+x)$, another predator consumes the prey with the functional response ${\alpha }_{2}xy/(a_2+x)$, and contributes to its growth rate ${\alpha }_2{\beta }_{2}xy/(a_2+x)$. Here ${\beta }_1$ and ${\beta }_2$ are conversion of biomass constants, ${\alpha }_1$ is the maximum value of the per capita reduction rate of x due to y and ${\alpha }_2$ is the maximum value of per capita reduction rate of x due to z;

(iv) Mortality rates of predators are assumed to be proportional to their populations. We have also considered density dependent mortality rate of both the consumer species as $\gamma y^2$ and $\delta z^2$. These terms describes either a self limitation of consumers or the influence of predation. Self limitation can occur if there is some other factor (other than food) which becomes limiting at high population densities. Predation on consumers can increase as $\gamma y^2$ and $\delta z^2$ if higher consumer densities attract more attention from predators or if consumers become more vulnerable at higher densities (see Ruan et al. [18] and references there in).

Several researchers at the present time claimed that the effect of time delay must be taken into account to form a biologically meaningful mathematical model (MacDonald, [19]). Form this view point we have introduced the delay in our model and this delay is referred to as the gestation period.

So our proposed model is as follows:

In our system, all the parameters are positive constants. There is a time delay of time τ in the prey species; γ and δ denote the intraspecific competition coefficients of the predators; ${\beta }_1$, ${\beta }_2$ are the conversion of biomass constant; $d_1$, $d_2$ are the death rate of first and second predator species respectively; ${\alpha }_1$ is the maximum values of per capita reduction rate of x due to y and ${\alpha }_2$ is the maximum values of per capita reduction rate of x due to z; $a_1$, $a_2$ are half saturation constants.

2 Equilibria, stability and Hopf bifurcation

Theorem 2.2

Assume that the positive equilibrium$F_{xyz}(x_5^{\ast },y_5^{\ast },z_5^{\ast })$of system(1)–(3)is locally stable. If${\eta }_1(0)>0$, then it is globally asymptotically stable (see Fig. 1).

3 Simulation and discussion

First we consider the case with $\tau =0$. To illustrate the analytical results numerically, let us take $r=3.5$, $K=150$, ${\alpha }_1=2$, ${\alpha }_2=2$, $a_1=35$, $a_2=45$, $d_1=1.3$, $d_2=0.05$, ${\beta }_1=2$, ${\beta }_2=3$, $\delta =0.2$. For these values of parameters a super critical Hopf-bifurcation occurs when ${\gamma }^{\ast }=0.00608308$. Now if we gradually increase the value of γ, keeping other parameters fixed, then $F_{xyz}$ achieves stability from instability as γ crosses its critical value ${\gamma }^{\ast }=0.00608308$ (see Figs. 2, 3 and 4).

Next consider ${\alpha }_1$ (the maximum value of the per capita reduction rate of x due to y) as the bifurcation parameter. With parameter values $r=3.5$, $K=150$, $a_1=35$, ${\alpha }_2=2$, $a_2=45$, $d_1=1.3$, $b_1=2$, $\gamma =0.006$, $d_2=0.05$, $b_2=4$, $\delta =0.2$, a Hopf-bifurcation occurs when ${\alpha }_1^{\ast }=2.071026$ (see Figs. 5, 6, and 7).

Now we would like to mention that the stability criteria of the system without delay do not necessarily guarantee the stability of the system with delay. It has been shown that the positive equilibrium which is stable without delay, remains stable under certain conditions when the time delay is less than the threshold value, otherwise the stable equilibrium become unstable. To illustrate the results numerically, choose $r=2.5$, $K=100$, ${\alpha }_1=2$, ${\alpha }_2=2$, $a_1=30$, $a_2=35$, $d_1=0.03$, $d_2=0.02$, $b_1=5$, $b_2=6$, $\gamma =0.1$, $\delta =0.2$ (Figs. 8 and 9).

For the above choices of parameters $F_{xyz}(x_5^{\ast },y_5^{\ast },z_5^{\ast })$ is locally asymptotically stable in the absence of delay. Now for the same values of parameters, it is seen from the Theorem 2.3, that there exists a critical value of $\tau ={\tau }_0=0.711633041$ and $F_{xyz}$ losses its stability as τ crosses the critical value ${\tau }_0$ (Fig. 10).