The Gestation Delay : A Factor Causing Complex Dynamics in Gause-Type Competition Models

In this paper, we consider a Gause-type model system consisting of two prey and one predator. Gestation period is considered as the time delay for the conversion of both the prey and predator. Bobcats and their primary prey rabbits and squirrels, found in North America and southern Canada, are taken as an example of an ecological system. It has been observed that there are stability switches and the system becomes unstable due to the effect of time delay. Positive invariance, boundedness, and local stability analysis are studied for the model system. Conditions under which both delayed and nondelayed model systems remain globally stable are found. Criteria which guarantee the persistence of the delayed model system are derived. Conditions for the existence of Hopf bifurcation at the nonzero equilibrium point of the delayed model system are also obtained. Formulae for the direction, stability, and period of the bifurcating solution are conducted using the normal form theory and center manifold theorem. Numerical simulations have been shown to analyze the effect of each of the parameters considered in the formation of the model system on the dynamic behavior of the system. The findings are interesting from the application point of view.


Introduction
Persistence, structure, and dynamics of biological species have received great attention from ecologists as well as mathematical biologists.The interspecific interference such as competition and predation dominates the dynamics of the system.Competitive exclusion principle by Gause states that when two species compete for the same resources within the environment, one of them will eventually outcompete and displace the other.The displaced population may become locally extinct by either migration or death, or it may adapt to a sufficiently distinct niche within the environment [1].Upadhyay and Iyengar [2] studied the extinction and coexistence of two competing prey species of a Gause-type ecological model in which the predator species is influenced by the damage effect caused by crowding from the members of its own population.
In ecosystem modeling, the type of predation is very important, and consecutively, the dynamics become richer with more realistic assumptions [3].The functional response (rate of prey consumption by an average predator) is one of the important features that govern the dynamics of the prey-predator system.Holling type II saturated functional response has been widely used to investigate the behavior of the prey-predator system.In real situation, conversion of prey biomass into predator biomass does not take place instantaneously.There is always a time lag of gestation.To make the system considered in [2] more realistic, we introduced the gestation period as a time delay and studied the effect of delay on the dynamics of the model system.Kar and Batabyal [4] have studied a two-prey one-predator system in the presence of time delay and derived the conditions for persistence and global stability of the model system.Recently, Deka et al. [5] considered a general Gause-type two-prey one-predator model system and studied the effect of predation on two competing prey species.The authors observed Hopf bifurcation by taking the death rate of the predator population as a bifurcation parameter.Jana et al. [6] also studied a two-prey onepredator model system in which the prey species is divided into weaker and stronger classes due to predator's catching efficiency.The evolutionary effects of selective disturbances on an evolving trait (e.g., body size and maturation age) of the predator individuals in a one-predator two-prey community are investigated by Meng et al. [7].Using methods of adaptive dynamics, they constructed an invasion fitness function and obtained the conditions for evolutionary branching and stability under selective disturbances in both monomorphic and dimorphic populations.Lu et al. [8] developed a rigorous analysis to investigate how the dual phenomena of prey infection and hunting strategies of the predators influence the disease dynamics in a predatorprey model with infected prey and stage-structured predator.A three-species model system is also studied by Gakkhar and Gupta [9] in which two competing prey species grow logistically but the third species acts as a predator as well as host, and they observed that the possibility of survival of all three species is due to commensalism of a prey and predator population.The dynamic behavior of a twoprey one-predator system with impulsive effects concerning periodic biological and chemical control is investigated by Pei et al. [10], and it was observed that there exist asymptotically stable pest-eradication periodic solutions if the impulsive period is less than some threshold.Wang et al. [11] studied a two-prey one-predator system with Watttype functional response and impulsive control strategy and established conditions for extinction of prey and predator population.The problem of adaptive control and chaotic behavior of a continuous-time two-prey and onepredator system is studied by El-Gohary and Al-Ruzaiza [12] with nonlinear feedback.
Many different food chain models including three species are constructed.Das et al. [13] studied chaotic dynamics of a three-species competition model with noise.Dubey and Upadhyay [14] studied the interaction of two predators competing for a single prey species with ratiodependent-type functional response and established the conditions for persistence and extinction of the species.Song et al. [15] explored the combined influence of multiple delays with two different scenarios, i.e., the intrinsic growth rate of the resource is less/more than those of the consumer and predator in a three-species model with consumer-eatresource and predator-eat-consumer.The significance of top predator interference and gestation delay is studied by Jana et al. [16] on a three-species food chain model involving intermediate and top predator populations, and they observed the subcritical Hopf bifurcation phenomena.Wang and his group explored the dynamics of different stochastic epidemic models like SIRS and FIV with different infectious forces under intervention strategies and environmental variability [17][18][19][20].
Motivated by the above works, we explain the dynamical behavior of the undertaken delayed model system and study the effect of control parameters with time delay on the dynamics of the model system.Since the model represents the very common food chain model of realworld systems, we are taking one of them as an example-a species found most throughout North America and from southern Canada to central Mexico, known as bobcats.Its main prey are rabbits and squirrels, which depend on grass and leaves.Kittens may be taken as prey by adult bobcats when prey levels are low, but this is very rare and does not influence much the dynamics of the population [21].Some species like golden eagle and coyotes used to predate on bobcats.Diseases, accidents, hunters, and starvation are the leading causes of death of the species.A bobcat is able to survive for long periods without food but will eat heavily when prey is abundant and used to store their hunted prey sometimes to eat later.The gestation period of the species is generally from 60 to 70 days.We study this time lag by introducing the time delay.In our work, interspecific competition coefficient of both prey species along with the intraspecific interference coefficient of the predator is considered.The effect of all of these parameters is studied.We have considered Holling type II functional response as predation of both the prey and studied the effect of predation rates b 13 , b 23 with the conversion rate of prey into predator β 31 , β 32 on the model system.The death rate of predators is also affecting the dynamics of the model system significantly.All these parameters are considered in numerical simulations to show the behavior and changes in the dynamics of the model system.
The rest of the organization of this paper is as follows: In Section 2, the model system is formulated.Also, positive invariance and boundedness have been shown in this section.Stability analysis for both delayed and nondelayed model systems has been made in Section 3. The conditions for the existence of Hopf bifurcation are also performed here.Next, Section 4 is devoted to establish the formulae to determine the direction, period, and stability of Hopf bifurcation of the bifurcating periodic solution and some numerical simulations are carried out to illustrate the analytical results in Section 5. Finally, we end the paper with a brief conclusion and discussion in Section 6.

The Mathematical Model
Consider a Gause-type model to draw the scenario of an ecological system where a predator (bobcats) species X t predates two competing prey species N 1 (squirrels) and N 2 (rabbits) (see Figure 1).It is well known that gestation delay, as a time lag, plays a complicated role on the 2 Complexity dynamics of the prey-predator system.Thus, the considered system can be governed by the following delay differential equations [2,5]: with the initial conditions where ϕ −τ, 0 → ℜ 3 with the norm Here, all the parameters are positive.r 1 , K 1 , α 12 , r 2 , K 2 , α 21 are the intrinsic growth rate, carrying capacity, and interspecific competition coefficient of the prey species N 1 , N 2 , respectively, b 13 , β 31 , b 23 , β 32 are the feeding rates per predator per unit prey consumed, assimilation, or conversion efficiency of the predator of the prey species N 1 , N 2 , respectively, and δ, γ are the death rate and intraspecific interference coefficient of the predator species X in the absence of the prey.The physical interpretation of a i is that a i = 1/A i , where A i is the half saturation constant of the prey species N i .The values of all the parameters are positive.

Positive Invariance.
It is necessary to insure that the species always survive.For this, we need to show that the system is positively invariant.Theorem 2.1.All the solutions of (1) are positive with positive initial conditions.

Boundedness
Theorem 2.2.For any positive solution ϕ t = N 1 t , N 2 t , X t of system (1), there exists a T > 0 such that 0 < N 1 t ≤ M 1 , 0 < N 2 t ≤ M 2 and 0 < X t ≤ M 3 for t ≥ T, where Proof.Obviously, from Theorem 2.1, N 1 t , N 2 t , X t remains nonnegative.So, we need to show that From the first equation of model system (1), we get Therefore, by using the standard comparison rule, we have Similarly, from the second equation of model system (1), we get Therefore, by using the standard comparison rule, we have Now, from the third equation of model system (1), we have 3 Complexity Therefore, by using the standard comparison rule, we have

Stability Analysis and Existence of Hopf Bifurcation
Model system (1) has seven possible nonnegative equilibrium points, namely, E 0 0,0,0 , Here, we are omitting the existence criteria of all these equilibrium points as it has been shown in Upadhyay and Iyenger [2].
3.1.Linear Stability Analysis.In this section, we will discuss the stability criteria of all abovementioned equilibrium points.
(i) The eigenvalues of the variational matrix around E 0 are r 1 , r 2 , −δ.Hence, E 0 is a saddle point (ii) The eigenvalues of the variational matrix around (iv) The predator-free equilibrium point is given as One eigenvalue of the variational matrix around and the other two eigenvalues are roots of the following equation: where We note that B > 0 under condition (12).Thus, the equilibrium point E 3 N 1 , N 2 , 0 is locally asymptotically stable if λ 1 < 0 and E 3 exists under condition (12).
Again, E 3 is unstable in the positive direction orthogonal to the N 1 − N 2 plane, i.e., the X direction if λ 1 > 0 (v) The equilibrium point is given as E 4 N1 , 0, X where X = 1/γ β 31 N1 /1 + a 1 N1 − δ and N1 is a positive solution of the following equation: where By Descarte's rule of sign, (17) has a unique positive solution One eigenvalue of the variational matrix around E 4 is and the other two eigenvalues are the roots of the following characteristic equation where Thus, the other two eigenvalues have negative real parts if It follows that the equilibrium E 4 N1 , 0, X is locally asymptotically stable if λ 1 (given by ( 20)) is negative and condition (23) holds.Also E 4 is unstable in the positive direction orthogonal to the N 1 − X plane, i.e., in the a positive solution of the following equation: where Similar to E 4 , it can be checked that the equilibrium Again, E 5 is locally asymptotically stable if X * , we discuss both the cases separately Case 1.When τ = 0.
The eigenvalues of the "J" matrix around E * are the roots of the following equation: where From the Routh-Hurwitz criterion, we can state the following theorem.
Separating real and imaginary parts, we have The following equation is obtained by some mathematical calculations As sin 2 ωτ + cos 2 ωτ = 1, from (38), we have where Without loss of generality, assume that it has at least one positive real root ω 0 , for which from (40), we have where k = 0, 1, 2, … , and then ±iω 0 are a pair of purely imaginary roots of (35) with τ k .Let us define Also, differentiating (34) with respect to τ, we get

Re dλ dτ
where Hence, the transversality condition is satisfied if U sin ωτ + V cos ωτ − D 2 ω 2 > 0. Thus, we can now state the following theorem.
Theorem 3.2.The equilibrium point E * of model system (1) is locally asymptotically stable when τ ∈ 0, τ 0 if the conditions of Theorem 3.1 hold and are unstable for τ > τ 0 .Furthermore, the system undergoes Hopf bifurcation at E * when the value of τ crosses through τ 0 provided that U sin ωτ + V cos ωτ − D 2 ω 2 > 0.

Global Stability Analysis
Theorem 3.3.The positive equilibrium point E * of system (1) is globally asymptotically stable provided that 6 Complexity where Proof.To prove that the equilibrium point E * x * , y * , z * of model system ( 1) is globally asymptotically stable, the method of Lyapunov function is used.We derive a sufficient condition which guarantees that E * is globally asymptotically stable.For mathematical convenience, we make the following transformations of the following variables: These coordinate change transforms the positive equilibrium E * into the trivial equilibrium x t = y t = z t = 0 for all t > 0. Due to the above transformations, model system (1) is reduced as follows: Computing the upper right derivative of V 1 t along the solution of (1), we get Again, let V 2 = y t .Computing the upper right derivative of V 2 t along the solution of (1), we obtain Now, let V 3 = z t .Computing the upper right derivative of V 3 t along the solution of (1), we obtain Again, due to the structure of (51), we consider the following functional: whose upper right derivative along the solution of system (1) is given by According to the assumptions of the theorem, we know that η = η ij 3×3 is an M matrix; hence, there exist positive constants h i i = 1,2,3 such that, Let us consider the Lyapunov function V defined by Then from (50), (51), and (54), we obtain It follows from (55) that, for all t ≥ T + r, we have Since system (1) is positively invariant, one can see that there exist positive constants m 1 , m 2 , m 3 and , and X * e z t = X t ≥ m 3 for t ≥ T * .Also, system (1) is bounded, and hence, there exist positive constants Using the mean-valued theorem, one obtains, where N * 1 e θ 1 t lies between N 1 t and N * 1 , N * 2 e θ 2 t lies between N 2 t and N * 2 , and X * e θ 3 t lies between X t and X * .
Letting α = min m 1 l 1 , m 2 l 2 , m 3 l 3 , then it follows that, for t ≥ T * , Noting that the Lyapunov function is such that Hence, by applying the method of Lyapunov function and (60), we can conclude that the zero solution of system (49) is globally asymptotically stable; therefore, the positive equilibrium of system (1) is globally asymptotically stable.The proof is complete.
Note.The positive equilibrium point E * will also be globally asymptotically stable for nondelayed model system (1).The conditions under which E * remains globally asymptotically stable are stated in the following lemma.Lemma 3.1.If the positive equilibrium point E * is locally asymptotically stable, then it is also globally asymptotically stable in the interior of the positive octant under the following conditions: where Proof.The proof of the lemma is referred to in Appendix.
Proof.We prove this theorem by the method of suitable Lyapunov function [22].Let the Lyapunov function for system (1) be where p 1 , p 2 , and p 3 are positive constants.Clearly, σ is a nonnegative C 1 function defined in R 3 + .Then we have By computing ψ at the boundary equilibria, we have If we choose p 1 and p 2 large enough so that p 1 r 1 + p 2 r 2 − p 3 δ > 0, hence, ψ is positive at E 0 .Also, if conditions (12), (19), and (23) satisfy, then p 2 L 1 + p 3 L 2 > 0 and p 1 L 3 + p 3 L 4 > 0, i.e., ψ is positive at E 1 and E 2 .And, if inequalities in (65) hold, then ψ is positive at all the remaining boundary equilibria.Therefore, system (1) is uniformly persistent, which follows from Theorem 3.12 of [22].The proof is complete.

Properties of Hopf Bifurcation
In this section, the formulae for determining the direction, stability, and period of the bifurcating periodic solutions of system (1) are presented here.The method that we used is based on the normal form theory and center manifold theorem as described in Hassard et al. [23] to determine the abovementioned properties of Hopf bifurcation.In the previous section, we noted that for the critical value of time delay, τ = τ 0 , the family of periodic solutions bifurcates from the steady state under certain condition.Without loss of generality, we can say that (35) has a pair of purely imaginary roots ±iω 0 at these critical values of τ and the system 9 Complexity undergoes Hopf bifurcation.Hence, for any root of (35) of the form λ τ = ν τ − iω τ , ν τ k = 0 and ω τ k = ω 0 and ∂ν/∂τ | τ=τ k ≠ 0. Let τ = τ 0 + μ, μ ∈ R, so that μ = 0 is a Hopf bifurcation value for the system.Define the space of continuous real-valued functions as , and x i t = x i τt for i = 1,2,3; delay system (1) then transforms to functional differential equation in C as

87
which gives where Proceeding in the same manner as Hassard et al. [23], we compute the coordinates to describe the center manifold C 0 at μ = 0. Let x t be the solution of (78) when μ = 0. Define z t = q * , x t , W t, θ = x t θ − 2Re z t q θ 90 On the center manifold C 0 , we have

91
where z and z are local coordinates for center manifold C 0 in the direction of q * and q * .We now consider only the real solution x t ∈ C 0 of (78), which gives where It follows from ( 90) and ( 91) that

94
so that 11 Complexity Now, from (92), we have Comparing the coefficient with (93), we have In order to compute g 21 , we need to compute W 20 θ and W 11 θ .

104
which in comparing the coefficients with (100) gives From (102), (105), and the definition of A, we have Note that q θ = q 0 e iω 0 τ 0 θ ; hence, Similarly, from (103), (106), and the definition of A, we have where and where η θ = η 0, θ .From ( 98) and (100), we get and where Using (108) and ( 113) in (111) and noticing that iω 0 τ 0 I − 0 −1 e iω 0 τ 0 θ dη θ q 0 = 0, 116 and 13 Complexity Solving this system for K 1 , we obtain  110).Furthermore, g 21 can be expressed by the parameters and delay term.Thus, we can compute the following values: which determine the qualities of bifurcating periodic solution in the center manifold at the critical value τ 0 .Now, we state the results in the following theorem.
Theorem 4. Suppose that μ 2 determines the direction of the Hopf bifurcation.If μ 2 > 0, then the Hopf bifurcation is supercritical and the bifurcating periodic solutions exit for τ > τ 0 .If μ 2 < 0, then the Hopf bifurcation is subcritical and the bifurcating periodic solutions exit for τ < τ 0 .β 2 determines the stability of the bifurcating periodic solutions.
The bifurcating periodic solutions are stable if β 2 < 0 and unstable if β 2 > 0 , and T 2 determines the period of the bifurcating periodic solutions.The period increases if T 2 > 0 and decreases if T 2 < 0.

Numerical Simulation
In order to study the dynamic behavior of model system (1), we consider a set of parameter values as With respect to the considered parameter set, the nonzero equilibrium point is calculated as E * 77 3512,38 32 01, 50 7297 .Also, corresponding with ω 0 and the critical value of time delay, τ 0 are given as 0 799559 and 0 20757, respectively.In Figure 2, we can see that the nondelay system is stable for this set of parameter values.Also, the system remains stable for τ = 0 2 < τ 0 = 0 20757 and becomes unstable for τ = 0 25 > τ 0 as τ 0 increases.The system bifurcates at the critical value of delay, i.e., at τ 0 .From Figure 2, we can also see the effect of time delay that destabilizes the system as it becomes larger than the critical value.Now, we are interested to study the effect of the other parameters as well on the dynamics of the model system.We found that the carrying capacity K 1 of the first prey N 1 does not affect the dynamics of the system (see Figure 3(a)) whereas the carrying capacity K 2 of the second prey N 2 changes the nature of the model system significantly.The system shows the limit cycle behavior from the stable state due to increase in K 2 (see Figure 3(b)).Now, we have plotted the graph to show the relation between these parameters (K 1 and K 2 ) and the critical value of time delay.In both the cases, we observed that the value of τ 0 decreases as K 1 or K 2 increases (see Figures 3(c) and 3(d)).It is interesting to note that the values of τ 0 vary in the same ranges of K 1 and K 2 .It is easy to check in Figure 4 that the system changes its dynamics as the value of time delay changes.In Figures 4(a) and 4(c), we can see that for some values of K 1 or K 2 , the system shows stable focus and for some other values, it shows the limit cycle behavior accordingly as the value of τ = 0 2 is less or greater than the respective value of τ 0 .Since τ = 0 3 and τ = 0 5 are more than the values of all τ 0 for the respective K 1 and K 2 taken in Figures 4(b) and 4(d), the system shows the limit cycle behavior for each cases.In the same manner, we can see that the growth rate of N 2 only affects the dynamics of the system but not of N 1 .System dynamics change from the stable state to the limit cycle state as we increase the value of r 2 but system remains unchanged when r 1 is changed (see Figure 5).Now, we are checking the effect with the time delay and it is easy to see by comparing Figures 5(c) and 5(d) with Figure 5(a).The system remains stable at τ = 0 15 while it turns into stable limit cycle when τ = 0 3 is taken for all the values of r 1 in both cases.While at τ = 0 2 and τ = 0 5 in Figures 5(e) and 5(f), resp., we observed that the system shows cyclic behavior for some r 2 and remains stable for other values of r 2 as in Figure 5(b).

Effect of Predation Rate of N 1 and Conversion Efficiency
of Predators b 13 , β 31 .In Figures 6(a) and 6(b), we can see that the system changes its behavior from stable focus to limit cycle as we increased the value of β 31 from 1 3 to 1 55 whereas b 13 does not have effect on the dynamics of the system; it remains at stable focus.Moreover, it is interesting to see the reverse effect of increasing β 31 and b 13 in a particular range on the critical values of time delay (see Figures 6(c) and 6(d)).The respective critical value increases with an increase of b 13 whereas it decreases with an increase of β 31 .Also, we observed that the decreasing value of β 31 maintains the system to the stable state for larger time lag.
But, the predation rate of N 2 and conversion efficiency of predator b 23 , β 32 do not follow any particular pattern for the respective critical values in any particular range.So, we study the effect of these parameters on the nondelayed model system only.In Figure 7(a), we can see that the system changes its behavior from stable focus to limit cycle as we are increasing the value of β 32 from 1 05 to 1 35.Also, due to the increase in b 23 from 0 15 to 0 45, the system loses its stability and shows limit cycle behavior (Figure 7(b)).We observed that the system becomes stable as we increase the death rate of the predator species (Figure 8(a)).As we increase δ from 5 1 to 5 4, the system becomes stable from limit cycle behavior.We also study the effect of intraspecific interference coefficient of the predator on the model system and found that the system becomes stable from limit cycle behavior as we increase the value (see Figure 8(b)).For a particular set of parameter values, the system shows limit cycle behavior for δ < 5 3 and becomes stable for δ ≥ 5 3. Similarly, to keep the system stable, we also need γ ≥ 0 005; otherwise, the system shows limit cycle behavior for γ < 0 005.

Conclusions and Discussions
In this work, we have studied a Gause-type two-prey onepredator model system.We have discussed the stability properties of equilibrium points for the model system including positive invariance and boundedness.Nonlinear stability analysis for both delayed and nondelayed model systems has been established.Conditions for the existence of Hopf bifurcation at the nonzero equilibrium point of the delayed model system have been obtained.Criteria for the system to be in uniform persistence have been derived, and it is observed that delay due to gestation does not affect uniform persistence of the system.The formulae for the direction, stability, and period of bifurcating solution have been calculated using the normal form theory and center manifold theorem.Numerical simulations are presented to verify the analytical predictions.It illustrates that the time delay destabilizes the system depending on its strength.We have studied the effect of all the control parameters on the dynamics of the model system.It has been observed in our simulation experiments that some of the parameters change the dynamics of the system either into a stable state or to an unstable state as explained in the following: (i) The conversion rate β 31 of prey N 1 into predator is sufficient to break the stability of the nondelayed model system as it crosses through β 31 = 1 49 (ii) To maintain the stability of the nondelayed model system, we need to take (for a particular set of parameter) intraspecific interference coefficient, γ ≥ 0 005, death rate of predator species, δ ≥ 5 3, and carrying capacity of the second prey population, K 2 ≤ 114 (iii) Time lag increases (up to which the system remains in a stable state) with increase in the predation rate of N 1 , 0 01 ≤ b 13 ≤ 0 5, whereas it decreases with an increase in the conversion rate of prey N 1 into predator, 1 2 ≤ β 31 ≤ 1 48 (iv) Time lag decreases with an increase in the carrying capacities of both prey populations (K 1 , K 2 ) in the range of 80-110, but time delay up to which the system remains stable is larger for decreasing K 2

Complexity
The rest of the parameters does not change the nature of the dynamical system much.We have also studied the effect of the growth rates of the prey population on the critical value of time delay and observed that the value of τ 0 decreases first up to a certain level for a particular value of r 1 and r 2 , and then it starts to increase due to the increase in the growth rate of the prey population.The population of the bobcats depends primarily on the population of its prey, and it has been found generally stable and healthy [24].Since the population is leading the wildcat species in the skin trade, the bobcat is legally hunted in small numbers.In some areas, they are quite rare, while in others, they have stable and sometimes dense populations.Hence, some states allow regulated hunting, while in others, they are protected.a regional level, the bobcat is totally protected in ten USA states, in Canada, hunting and trade is regulated, and in Mexico, hunting is regulated in five states [25].Using the above results, we can reach to the goal of making a healthy stable ecosystem and plan for regulating the bobcat population in North America and southern Canada.It is obvious to satisfy conditions (A.6) for the values of the constants c 1 and c 2 as given in the theorem.Further, conditions (A.5) are satisfied from conditions given in (62).Therefore, V is negative and definite under condition (62), and then V is a Lyapunov function with respect to E * .Hence, E * is globally asymptotically stable.

Figure 3 :
Figure 3: Phase plot diagram of nondelayed model system (1) for different (a) K 1 and (b) K 2 .Relation between the critical value of time delay and (c) K 1 and (d) K 2 .Rest of the parameters is taken as given in the text.

Figure 5 :
Figure 5: Phase plot diagram of nondelayed model system (1) for different (a) r 1 and (b) r 2 .Phase plot diagram of system (1) for different r 1 at (c) τ = 0 15 and (d) τ = 0 3. Phase plot diagram of system (1) for different r 2 at (e) τ = 0 2 and (f) τ = 0 5.The rest of the parameters is taken as given in the text.

Figure 7 :
Figure 7: Phase plot diagram of nondelayed model system (1) for different (a) β 32 and (b) b 23 .The rest of the parameters is taken as given in the text.

Figure 6 :
Figure 6: Phase plot diagram of nondelayed model system (1) for different (a) β 31 and (b) b 13 .Relation between the critical value of time delay and (c) β 31 and (d) b 13 .The rest of the parameters is taken as given in the text.

1 훾Figure 8 :
Figure 8: Phase plot diagram of nondelayed model system (1) for different (a) δ and (b) γ.The rest of the parameters are taken as given in the text.