Research Article | Open Access

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

**Academic Editor:**Liwei Zhang

#### Abstract

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 the Hopf 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 using MATLAB, which shows the reliability of our model from the practical point of view.

#### 1. 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–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–15].

Time delays have been incorporated in ecological models by many authors [16–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.

#### 2. 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 and the death rate is , where , , , and are assumed to be positive constants.

We consider the model with initial data , , and . 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 , 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 where , , , , 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.*

* Proof. *From the first equation of the system (3) we get

From the second equation of (3), we get

Similarly, from the third equation of (3), we get

Therefore, , , and , for all .

From the first equation of (3), we have

By a standard comparison theorem and by (5), we have,

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 . This completes the proof.

#### 3. Stability Behaviour of the Model

##### 3.1. Case I: Zero Exogenous Input

Theorem 2. * If , then the model (3) has nonnegative equilibria which is unstable and which is locally asymptotically stable. The interior equilibrium is not feasible. *

* Proof. *The variational matrix of system (3) at is

The eigen values are , , and . So, it is obvious that, is unstable (hyperbolic saddle).

Now the varitiational matrix of system (3) at is

The characteristic equation of is , where and , since .

The eigen values are and .

Since , and , therefore the signs of the real parts of , are negative. Hence is locally asymptotically stable.

It is noted here that the other equilibrium point is not feasible, since
are opposite in sign since , by (2).

This completes the proof.

##### 3.2. Case II: Nonzero Exogenous Input

When , the model (3) has two nonnegative equilibria, and . The variational matrix of system (3) at is given by The characteristic equation of is .

So, is locally asymptotically stable if and only if and if , and then becomes unstable.

The interior equilibrium point of system (3) is given by 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.

Theorem 3. * If , then the equilibrium point of the system (3) is locally asymptotically stable if and only if and the unique interior equilibrium point of system (3) exists if and only if the following two conditions are satisfied:
*

##### 3.3. Local Stability of

The variational matrix at is given by where The characteristic equation is where By the Routh-Hurwitz criterion [28–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.*

Furthermore, the conditions arising through Hopf bifurcation cannot be derived explicitly in terms of system parameters. Let us assume that stands for any parameter involved with the model system (3). If there exists a critical magnitude such that , , , and , then undergoes a Hopf bifurcation at (Liu’s criterion) [28–30].

##### 3.4. Global Stability of

is not always globally asymptotically stable. The conditions which guarantee the global stability of are stated in the following theorem.

Theorem 5. * Since , , and are bounded, there exist positive constants () such that , , and . If the following inequalities hold:
**
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–30]. This competes the proof.

#### 4. 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 . For a biological meaning, we further assume that .

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 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 in the presence of discrete delay . 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 , , . 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 . 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 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 . 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 is the first positive root of (41) and hence is the first positive root of (40), and then Thus, the transversality condition is satisfied and the steady state becomes unstable when . Moreover, a Hopf bifurcation [28–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 and as . Therefore, the cubic equation (41) always has at least one positive root. Consequently the stability criteria of the system for will not necessarily ensure the stability of the system for . In the following theorem, we have given a criterion for switching the stability behaviour of .

Theorem 6 (see [28–30]). *If exists with the conditions (24) and be a positive root of (41), then there exists a such that*(i)* is locally asymptotically stable for ,*(ii)* is unstable for ,*(iii)*the system (31) undergoes a Hopf-bifurcation around at ,
**and the minimum is taken over all positive such that is a solution of (41).*

#### 5. Estimation of the Length of Delay to Preserve Stability

We linearize the system (31) about its interior equilibrium and get where Taking Laplace transform of the system (47), we get where and , , and are the Laplace transform of , and , 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] where and 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 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 , independent of so that (53) holds for all values of , and hence in particular at .

We rewrite (54) as

Maximizing subject to , , we obtain

Hence, if then clearly from (56) and (57), we have .

From the inequality (53), we get

As is locally asymptotically stable for (provided that inequalities given in (24) are satisfied), for sufficiently small , (58) will continue to hold. Substituting (55) into (58) and rearranging, we get Using the bounds we obtain

Now, from (59) and (60), we get where then stability is preserved for .

Summarizing the previous discussions we come to the following theorem.

Theorem 7. * The delayed model (31) and (32) will be locally asymptotically stable at with the conditions (24) if the delay lies within the interval where is given by (63).*

#### 6. 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 for . , and we can rewrite the system (31) and (32) as follows: where , are given, respectively, by where By the Riesz representation theorem, there exists a matrix function whose components are bounded variation function in such that In fact we can choose where is the Dirac delta function.

For , we define Then the system (64) is equivalent to For , where is the 3-dimensional space of row vectors, the adjoint operator of A is defined as For and , we define a bilinear inner product taking the first argument to be conjugate linear and the second to be linear of the form where . It is easy to verify that and are a pair of adjoint operators.

From the discussions in the previous section, we know that are the eigenvalues of . Hence, the eigenvalues of are .

Suppose that is eigenvector of corresponding to , and then we have It follows from (67), (68), and the definition of that where is identity matrix of order , that is, Solving (76), we have , where Now, suppose that is the eigenvector of corresponding to , and then we have It follows from (67), (68), and the definition of that, where is identity matrix of order , that is, Solving (80), we have , where Now, Thus in order to assume , we can choose Also since , we have Therefore, .

Now, we first compute the coordinates to describe the center manifold at . Let be the solution of (64) when . Define On the center manifold , we have where and are local coordinates of center manifold in the direction of and . Note that is real if is real. We only consider real solutions.

From (85), we get For the solution of (64) and since , we have where Since , we have where

Then,

It follows that