Mathematical Analysis of a Single-Species Population Model in a Polluted Environment with Discrete Time Delays

We have discussed the dynamical behaviour of a single-species population model in a polluted environment which describes the effect of toxicants on ecological system. Boundedness, positivity, and stability analysis of the model at various equilibrium points is discussed thoroughly. We have also studied the effect of single discrete delay as well as double discrete delays on the population model. Existence conditions of theHopf bifurcation for single time delay are investigated.The length of delay preserving the stability is also estimated. The direction and the stability criteria of the bifurcating periodic solutions are determined by using the normal form theory and the center manifold theorem. The stability of the model with double time delays is investigated by using the Nyquist criteria. By choosing one of the delays as a bifurcation parameter, the model is found to undergo a Hopf bifurcation. Some numerical simulations for justifying the theoretical results are also illustrated by usingMATLAB, which shows the reliability of our model from the practical point of view.


Introduction
In the world today, the pollution of the environment is a threatening problem due to the rapid development of industrialization.The presence of toxicants in the environment decreases the growth rate of species and its carrying capacity.In recent years, the most serious problem to society is the change in the environment caused by pollution, affecting the long-term survival of species, human life style, and biodiversity of the habitat.A large quantity of the toxicants and contaminants enter into the ecosystem one after another which seriously threaten the survival of the exposed population including humans.Therefore, the study of the effects of toxicants on the population and the assessment of the risk to populations are becoming quite important.
The problem of estimating qualitatively the effect of a toxicant on a population by mathematical models is a very effective way.The study of deterministic dynamic population models with toxicant effect was proposed by Hallam and his colleagues in 1980s [1][2][3].This model was revisited by many researchers [4,5].He and Ma in [6,7] discussed the survival of a single-species in a polluted environment, considering the organism's uptake of toxicant from the environment and egestion of toxicant into the environment.In [4] Ma studied a Leslie resource-consumer model and obtained the threshold between the persistence and extinction of the consumer.In [8], Buonomo et al. studied the effect of variation of the population on the toxicant concentration in the organism and environment.Samanta and Maiti [9] studied a dynamical model of a single-species system in a polluted environment under two cases: constant exogenous input of toxicant and rapidly fluctuating random exogenous input of toxicant into the environment by means of ordinary and stochastic differential equations.There are many other works on the effects of a single toxicant on various ecosystems by using mathematical models [10][11][12][13][14][15].
Time delays have been incorporated in ecological models by many authors [16][17][18][19][20][21][22][23][24][25][26][27].Introducing such delays in ecological models makes them more realistic and reasonable.Delay may have very complicated impact on the dynamics of a system since delay can cause the change of stability and can bifurcate periodic solutions.
In this paper, we have developed a single-species population model in polluted environment.In Section 2, we present a brief sketch of the construction of the model.Boundedness and positivity analysis of the solutions is shown which implies that the system is ecologically well behaved.In the next section, we have discussed the existence and stability analysis of various equilibrium points under zero-exogenous input and nonzero constant exogenous input.Next, we obtain necessary conditions for the existence of interior equilibrium  * and local and global stability of the system at  * .The analysis of single delayed model is described in Section 4. This analysis shows that there is a critical value  * of delay  below which the system is stable and above which the system becomes unstable at the interior equilibrium  * .The system undergoes a Hopf bifurcation around  * at  =  * .In Section 5, we have estimated the length of delay to preserve the stability around  * .The direction and the stability criteria of the bifurcation periodic solutions by using the normal form theory and the center manifold theorem are discussed in Section 6.In the next section, the analysis of double delayed model is discussed.The stability of the model with double delays is investigated by the Nyquist criteria.By choosing one of the delays as a bifurcation parameter, the model is found to undergo a Hopf bifurcation.Some of the important analytic results are numerically verified by using MATLAB in Section 8. Finally, Section 9 contains the general discussions of the paper and ecological implications of our mathematical findings.

Basic Mathematical Model
In this paper, we analyze a model which describes the effect of toxicant on a single-species population.Here we assume the following: (): concentration of the population biomass, (): concentration of the toxicant in the population, (): concentration of the toxicant in the environment.
The model satisfies the following assumptions.
(A1) There is a given toxicant in the environment and the living organisms absorb into their bodies part of this toxicant so that the dynamics of the population is affected by the toxicant.
(A2) For the growth rate of population, we assume that the birth rate is  0 − () and the death rate is  0 + (), where  0 , ,  0 , and  are assumed to be positive constants.
We consider the model with initial data (0) > 0, (0) > 0, and (0) > 0. Here, : depletion rate of toxicant in the environment due to its intake made by the population, : depletion rate of toxicant in the population due to egestion, : depletion rate of toxicant in the population due to metabolization process, ℎ: depletion rate of toxicant in the environment, (): exogenous toxicant input rate which is assumed to be a smooth bounded nonnegative function of .
We can see that if  0 −  0 − () ≤ 0, then () will be going to extinct.So we suppose that ( The model we have just specified has nine parameters, which make the analysis difficult.To reduce the number of parameters, we make the following transformations:  = ,  = , and  = .So, system (1) becomes with initial data: where  =  0 −  0 ,  =  +  +  0 ,  = /,  =  +  0 , and () = ().Therefore, from (2), we have Theorem 1.Each component of the solution of system (3) subject to initial conditions (4) is positive and bounded.
From the first equation of (3), we have By a standard comparison theorem and by (5), we have, lim sup From the third equation of (3), we have So, we conclude that each component of the solution of system (3) subject to ( 4) is positive and bounded for all  ≥ 0. This completes the proof.
Proof.The variational matrix of system (3) at  0 is The eigen values are , −, and −ℎ.So, it is obvious that,  0 is unstable (hyperbolic saddle).Now the varitiational matrix of system (3) at  1 is The characteristic equation of ( 1 ) is The eigen values are  1 = − < 0 and Since  > 0, and  > 0, therefore the signs of the real parts of  2 ,  3 are negative.Hence  1 is locally asymptotically stable.
It is noted here that the other equilibrium point ( * ,  * , and  * ) is not feasible, since are opposite in sign since  > 0, by (2).This completes the proof.

3.2.
Case II: Nonzero Exogenous Input (() =  > 0).When (() =  > 0, a constant), the model (3) has two nonnegative equilibria,  2 (0, /ℎ, /ℎ) and  * ( * ,  * ,  * ).The variational matrix of system (3) at  2 is given by The characteristic equation of ( where It can be provided that the unique interior equilibrium point  * ( * ,  * ,  * ) of system (3) exists if and only if the following two conditions are satisfied: Summarizing the above analysis we come to the following theorem. where The characteristic equation is where By the Routh-Hurwitz criterion [28][29][30], it follows that all eigenvalues of characteristic equation have negative real part if and only if Therefore, we come to the following theorem.
Theorem 4.  * is locally asymptotically stable if and only if the inequalities (24) are satisfied.

Global
then  * is globally asymptotically stable.
Proof.We consider the following positive definite function which has a strict minimum at  * : Differentiating both sides with respect to  along the solution of (3), we get (after some simple calculations) where Now, sufficient conditions for / to be negative definite are that is, Hence  is a Lyapunov function and so  * is globally asymptotically stable [28][29][30].This competes the proof.

Model with Discrete Delay
It is already mentioned that time delay is an important factor in biological system.It is also reasonable to assume that the effect of the environmental toxicant on the population growth will not be instantaneous, but mediated by some discrete time lag  required for incubation.So, we modify system (3) and (4) as follows: with initial conditions where () is a nonnegative continuous functions on  ∈ [−, 0].For a biological meaning, we further assume that (0) > 0.
All parameters are the same as in system (3) except that the positive constant  represents the activation period or reaction time of the toxicant in the population and () = (> 0) is a constant.
The system (31) has the same equilibria as in art.3.2.The main purpose of this section is to study the stability behaviour of  * ( * ,  * , and  * ) in the presence of discrete delay ( > 0).We linearize system (31) by using the following transformations: Then transformed linear system is given by where Let us take the solution of the system (34) as () =  1   , V() =  2   , () =  3   .This leads to the following characteristic equation: where A necessary condition for a stability change of  * is that the characteristic equation should have purely imaginary solutions.In order to see the delay induced instability, let us consider  as bifurcation parameter and assume a purely imaginary solution of (3) in the form  =  ( ̸ = 0).Therefore, substituting  =  in (37) and separating real and imaginary parts, we get Eliminating  by squaring and adding (38), we get the equation for determining  as where Substituting  2 =  in (40), we get a cubic equation given by Now, from (38), we get Thus when  =   , the characteristic equation (37) has a pair of purely imaginary roots ± 0 .Let () = () + () be a root of (37) such that the following two conditions hold: Differentiating both sides of (37) with respect to , we get Now, if  =  2 0 is the first positive root of (41) and hence  =  0 (> 0) is the first positive root of (40), and then Re (   ) Thus, the transversality condition is satisfied and the steady state becomes unstable when  >  * .Moreover, a Hopf bifurcation [28][29][30] occurs when  passes through the critical value  * , where  * is the smallest positive value of   given in (42).We notice that, in (41), Φ is continuous everywhere with Φ(0) < 0 and Φ() → ∞ as  → ∞.Therefore, the cubic equation (41) always has at least one positive root.Consequently the stability criteria of the system for  = 0 will not necessarily ensure the stability of the system for  ̸ = 0.In the following theorem, we have given a criterion for switching the stability behaviour of  * .Theorem 6 (see [28][29][30]).If  * exists with the conditions (24) and  0 =  2 0 be a positive root of (41), then there exists a  =  * such that (i)  * is locally asymptotically stable for 0 ≤  <  * , (ii)  * is unstable for  >  * , (iii) the system (31) undergoes a Hopf-bifurcation around and the minimum is taken over all positive  0 such that  2 0 is a solution of (41).

Estimation of the Length of Delay to Preserve Stability
We linearize the system (31) about its interior equilibrium  * ( * ,  * , and  * ) and get where Taking Laplace transform of the system (47), we get where and  1 (),  1 (), and  1 () are the Laplace transform of  1 (),  1 () and  1 (), respectively.Now, we will use the "Nyquist theorem" [21] which states that if  is the arc length along a curve encircling the right half of the plane, then a curve Γ() will encircle the origin a number of times equal to the difference between the number of poles and the number of zeros of Γ() in the right half of the plane.Using "Nyquist theorem" [21,31], it can be shown that the conditions for local asymptotic stability of  * ( * ,  * ,  * ) are given by [32] Im  ( 0 ) > 0, Re  ( 0 ) = 0, where () =  3 +  1  2 +  2  +  3 + ( 4  +  5 ) − and  0 is the smallest positive root of (52).
We have already shown that  * ( * ,  * ,  * ) is stable in absence of delay (provided that inequalities given in (24) are satisfied).Hence, by continuity, all eigenvalues will continue to have negative real parts for sufficiently small  > 0 provided that one can guarantee that no eigenvalue with positive real part bifurcates from infinity as  increases from zero.This can be proved by using Butler's lemma [32].
In this case, conditions (51) and (52) give Now, (53) and (54), if satisfied simultaneously, are sufficient conditions to guarantee stability.We will utilize them to get an estimate on the length of delay.Our aim is to find an upper bound  + on  0 , independent of  so that (53) holds for all values of , 0 ≤  ≤  + and hence in particular at  =  0 .
Summarizing the previous discussions we come to the following theorem.

Direction of Hopf Bifurcation
In this section, we will derive the direction of the Hopf bifurcation [17,21,33] and sufficient conditions of the stability of bifurcating periodic solution from the interior equilibrium  * of the system (31) at the critical value  =  * .We will employ the approach of the normal form method and center manifold theorem introduced by Hassard et al. [33]. Let ), and we can rewrite the system (31) and (32) as follows: where   : R × C → R 3 ,  : R × C → R 3 are given, respectively, by where In fact we can choose where () is the Dirac delta function.
From the discussions in the previous section, we know that ± 0 are the eigenvalues of (0).Hence, the eigenvalues of  * are ∓ 0 .
From the definition of  and (101) for  = 0, we get Notice that Substituting ( 106) and ( 111) into (109), we obtain which leads to that is, ) . ( On simplification, we get where Similarly substituting (107) and ( 112) into (110), we get which leads to that is, ) .
On simplification, we get Thus,  20 () and  11 () can be determined from (90) and (93).Furthermore,  21 can also be determined.Therefore, based on the previous analysis,   can be determined in terms of the parameters and delays involved in the system.Thus, we can compute the following values: where sign of  2 determines the direction of the Hopf bifurcation, sign of  2 determines the stability of the bifurcating periodic solution, and sign of  2 determines the period of the bifurcating periodic solution in the center manifold at the critical value  * .Therefore, by the results of Hassard et al. [33], we can summarize the properties of the Hopf bifurcation at the critical value  =  * in the following theorem.
Nyquist criterion [21,31] implies that if  0 is a solution of (131), then ] . (135) Then, clearly,  + ≥  0 .Equation (130) indicates that the following inequality should hold: where If we can find φ( 1 ,  2 ) > ψ( 1 ,  2 ) such that where 0 <  0 < Next, we will study the existence of Hopf bifurcation of system (125a) and (125b) by choosing one of the delays as a bifurcation parameter.Let us take  2 as a bifurcation parameter.The characteristic equation corresponding to system (125a) and (125b) is Differentiating equation ( 146) implicitly with respect to  2 , we obtain   1. also satisfied and consequently  * is globally asymptotically stable.The phase diagram is shown in Figure 8.

Discussions and Conclusions
In this paper, we have developed a single-species population model under the influence of toxicant in the population and in the environment by means of ordinary differential equations in terms of their concentrations.It is shown that the model system (3) with ( 4) is bounded and the solutions are positive, which in turn implies that the system is ecologically well behaved.The analysis of the existence of various parameter values are given in Table 1.parameter values are given in Table 1.
equilibrium points under zero exogenous toxicant input and nonzero exogenous toxicant input leads us to Theorems 2 and 3, respectively.Theorem 2 states that under zero exogenous input the trivial equilibrium  0 is unstable, the interior equilibrium  * is not feasible, and only the axial equilibrium  1 is locally asymptotically stable.On the other hand, when exogenous input is nonzero, the existence of the interior equilibrium  * depends on some conditions and the equilibrium  2 is locally asymptotically stable under some condition.The stability criteria for  * given in Theorems 4 and 5 are the conditions for locally and globally stable coexistence of the population biomass, population toxicant, and environmental toxicant.It is proved by several researchers that the effect of time delay must be taken into account in the population models to    (4,3,5), (6,4,8), (5,3,8), (7,8,4), parameter values are given in Table 1.
make them ecologically more meaningful and useful.It seems also reasonable that the effect of the environmental toxicant on the organismal activities and the population growth will not be instantaneous.It must be mediated with some discrete time lag required for incubation.From this point of view, we have formulated single and double time delayed model given in (31) and (125a) and (125b), respectively, where the delays may be looked upon as the reaction time of the environment toxicant on the population biomass and the population toxicant.
The analysis of the delayed models (31) and (125a) and (125b) shows more complicated behaviour than their nondelayed counterpart.The rigorous analysis of the single delayed model (31) leads us to Theorem 6 which mentions that the stability criteria in absence of delay are no longer enough to guarantee the stability in the presence of delay, but rather there is a critical value  * of the gestation delay  such that the system is stable for  <  * and becomes unstable for  >  * at the interior equilibrium  * ( * ,  * ,  * ).From the next part of our analysis, we obtain the estimated length of delay to preserve stability at  * and the direction of Hopf bifurcation by using normal form theory and center manifold theorem.The stability of the double delayed model (125a) and (125b) is investigated by the Nyquist criteria which leads us to Theorem 9.By choosing one of the delays as a bifurcation parameter, the model (125a) and (125b) is found to undergo a Hopf bifurcation under some conditions stated in Theorem 10.
Analytical studies can never be completed without numerical verification of the results.So, some numerical simulations for justifying the important theoretical results are also illustrated by using MATLAB, which shows the mathematical and ecological reliability of the proposed model.
Stability of  * . * is not always globally asymptotically stable.The conditions which guarantee the global stability of  * are stated in the following theorem.