ISRN Applied Mathematics

Volume 2011, Article ID 643985, 31 pages

http://dx.doi.org/10.5402/2011/643985

## Modelling the Effects of Pollution on a Population and a Resource in a Polluted Environment

^{1}Department of Mathematics, Trent University, Peterborough, ON, K9J 7B8, Canada^{2}Department of Mathematics, Statistics and Physics, Qatar University, P.O. Box 2713, Doha, Qatar

Received 18 July 2011; Accepted 23 August 2011

Academic Editors: X. Meng and E. Skubalska-Rafajlowicz

Copyright © 2011 Victoria Maystruk and Kenzu Abdella. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

A model for the effect of pollution on an animal population partially dependent on a plant resource is examined. Using a system of ordinary differential equations, the model tracks and relates changes in an animal population and its internal pollution levels, a plant population and its internal pollution levels, and the overall environmental pollution level. The model system is analysed using standard mathematical techniques, including the direct Lyapunov method and numerical simulations. Criteria for the stability of the system are found and numerically tested. Three inequalities are sufficient to establish global stability, and a parameter range exists in which these criteria are satisfied. The stability criteria dictate that the system will be globally stable provided that the removal rate of the pollution from the environment, the intrinsic growth rate of the plant population, and the rate the animal population relieves itself of its pollution are all sufficiently large.

#### 1. Introduction

Industrial production does not currently occur without emitting waste and byproducts. The industrial emissions enter the environment surrounding the industrial site and affect the local ecosystem. It has been shown that industrial wastes can decrease the birth rates, increase the death rates, and reduce the carrying capacity of flora and fauna proximal to the source of the industrialization [1]. Moreover, the local species can also affect the concentration of pollution in the environment; the species may absorb, digest, and eventually metabolize the pollution.

Modeling how pollution emitted into an ecosystem travels through that environment as well as modeling the subsequent effects of the pollution on the surrounding populations is an important step toward the conservation and preservation of populations living near industrial sites. Accurate simulation of the pollution-population system is crucial to mitigating the potential loss of biodiversity. The majority of previous research endeavours that have mathematically modelled the effects of industrial pollution on species surrounding the industrial site have used the rate of uptake of pollution by the species as one of the variables in the model [2–9]. In these instances, the rate of uptake has been modeled as a linear process. While this approach simplifies some of the mathematical analysis, it does not track the movement of pollution through the ecosystem.

Recently, research modelling the movement of pollution through an ecosystem using a mass-balance approach for the pollution level was published [10]. This paper only deals with a single species living in a polluted environment; however, a similar approach can be used for a multi species model.

This paper draws on the multi species models previously developed by Dubey, Shukla, and their collaborators, while implementing a mass-balance approach for the pollution levels in the system presented by He and Wang to produce a new population-pollution model which is analysed using classical stability techniques.

We begin by proposing a new model to study the effects of a pollutant on local animal and plant biomass. We conduct local and global stability analysis on the nontrivial equilibrium solution of the model using the direct Lyapunov method. And finally, we conduct numerical simulations in order to compare the analytical results with numerical computations, in order to confirm that there is a parameter range for which the results are relevant.

#### 2. Mathematical Model

Consider an ecosystem in which there are two principal biological species—a plant species and an animal species. The animal species consumes and is partially dependent on the plant species. Pollution is also being inputted into the ecosystem. We assume that the ecosystem is a spatially homogeneous environment and also that there is no migration to or emigration from the ecosystem. Additionally, we assume that all of the individuals within either species population are identical.

To model this system we will use five state variables—, and . represents the animal species density at time represents the plant species density at time represents the concentration of the toxicant in an individual animal at time represents the concentration of the toxicant in an individual plant at time , and represents the concentration of the toxicant in the environment at time .

We assume that the growth dynamics of the animal species and the plant species can be described using a modified version of the logistic growth equations.

When discussing the dynamics of population growth it can be useful to discuss the dynamics in terms of the birth rate, , and death rate, , of the populations. In general, we can think of population growth as If we let where is the intrinsic birth rate and is the intrinsic death rate of the population, then we can think of the general birth rate of a population growing according to a logistic growth model as Similarly, the general death rate can be thought of as In the development of our model, we will first consider the population dynamics.

##### 2.1. Population Dynamics

When developing an equation for the carrying capacity of the plant population, we assume that the carrying capacity of the plant species decreases due to the presence of pollution in the environment. Additionally, we assume that this change is proportional to the concentration of pollution in the environment. Let denote the intrinsic birth rate of the plant species, the intrinsic growth rate of the plant species, the intrinsic carrying capacity of the plant species, and the rate the carrying capacity decreases relative to the concentration of pollution in the environment. The birth rate of the plant species can be written as

For the death rate of the plant population, we assume that the death rate increases proportional to the animal species population and proportional to the concentration of toxicant present within the plant. These dynamics are akin to the animal population eating and the pollution poisoning the plant population. If is the intrinsic death rate of the plant species, let and denote the rates the death rate of the plant species increases relative to biomass of the animal species and the concentration of pollution in the plant species, respectively. Then, we describe the death rate of the plant population as

Using , the equation governing the plant population is Now, consider the animal population dynamics.

Similar to the plant population, we assume that the animal species carrying capacity decreases in proportion to the quantity of pollution in the environment. Additionally, for the animal population, we assume that the carrying capacity increases in proportion to the population of the plant species. This increase is due to our assumption that the animals are partially dependent on the plants. Let , and represent the intrinsic birth rate, growth rate, and carrying capacity of the animal species, respectively. Additionally, let denote the rate the carrying capacity increases relative to the plant biomass and the rate the carrying capacity decreases relative to the concentration of pollution in the environment. Then, we can write the birth rate of the animal species as

To formulate the death rate of the animal species, we assume that the death rate decreases in proportion to the quantity of plant biomass and increases in proportion to the concentration of pollution present within the animal itself. We let represent the intrinsic death rate of the animal species, the rate the animals death rate decreases relative to the biomass of the plant, and the rate the animals death rate increases relative to the concentration of pollution present within the animal. We can write the animal death rate as Clearly, using , we can describe the dynamics of the animal population as Continuing the derivation of our model, we next develop the equations for the concentration of pollution dynamics via balance arguments.

##### 2.2. Pollution Dynamics

The use of balance equations to model movements within a system is based upon the first law of thermodynamics [11–13].

We use and to represent the average per capita mass of the animal and plant populations, respectively, and we use to represent the mass of the environment. Then represents the total mass of toxicant in the plant population, represents the total mass of toxicant in the animal population, and represents the total quantity of toxicant free in the environment. is assumed to be sufficiently large and near constant, and hence variations in do not need to be modeled. The flow of pollution in the system can be visualized in Figure 1 and is described next.

We assume that pollution first enters the environment at a rate . This pollution naturally dissipates at a rate relative to the amount of pollution present in the environment . The environmental pollution level decreases due to absorption into both the plant and animal populations. The loss of environmental pollution due to absorption to the plant population is given as . The loss of environmental pollution due to absorption into the animal population is given as . Pollution reenters the environment, as part of the pollution cycle as it is egested by the animal population. The egestion occurs at the rate . Pollution also reenters the general environment as it is released from dead animals and plants. The pollution is released from dead animals at the rate and from dead plants at the rate . Using these rates pertaining to the amount of pollution in the environment, we get then, dividing by and using the notation , we have Next, we consider the quantity of toxicant in the plant biomass.

We have already established that the plant biomass absorbs pollutant from the environment, increasing the plant toxicant at a rate . Pollution inside the plant biomass is naturally metabolized, leaving the system at a rate . Pollution returns to the environment due to dead biomass at a rate . Finally, pollution from the plant population is transferred to the animal population due to ingestion of plant matter by the animal population at a rate . So we get dividing by , we have but clearly Substituting in (2.6) and (2.15), using and factoring, we have We now consider the quantity of toxicant in the animal biomass.

Pollution enters the animal population via two routes—by the absorption of pollution from the environment at a rate and by the ingestion of plant biomass containing pollution at a rate . Pollution is released from the animal population as a result of being egested, which occurs at a loss rate of , metabolized, which occurs at a loss rate of , and released from dead animals, which occurs at a rate of . Hence, we obtain

First, we divide by to get Then, we use the fact that Substituting in (2.10) and (2.20), using , and factoring we determine that After the previous assumptions and derivations we can combine (2.7), (2.11), (2.13), (2.18), and (2.22) and arrive at our completed model:

#### 3. Equilibrium Analysis

Our model has four nonnegative equilibria: Since we derived the equations for and under the assumption that both and are nonzero, then the solution we are chiefly concerned with is .

##### 3.1. Existence of

For , and are the positive solutions of,

Recall that , and so, by (3.3), Substituting this into (3.5), we have Hence, where, Let Consider (3.3). Since , Let Next, recall that , and so, by (3.2), Substituting this into (3.4), we have Using and , we find where Let Next, if we consider (3.2), we have Using and , we find where Finally, we consider (3.6) and Substituting in , , and , we find

Thus, the nontrivial solution to the system is given by where where Provided that and exist, we can find the corresponding , and .

#### 4. Stability Analysis

Theorem 4.1 (conditions for local stability). *Given
**
where
**
where
**if the following inequalities hold: **(1)*

(2)

(3)

then is a locally stable equilibrium of system (2.23).

*Theorem 4.2 (conditions for global stability). Using the notation
let , and be such that
Then, given
where
if the following inequalities hold: (1)(2)(3)
then is a globally stable equilibrium of system (2.23).*

*5. Numerical Simulation*

*5. Numerical Simulation**To study the applicability of the model we need to test a range of parameter values for which exists and is stable. *

*Example 5.1. *Consider the set of parameter values for our numerical calculations in Table 1.

With these parameter values, the nontrivial equilibrium, , of the model exists and is

Using , , and , it can be verified that the conditions given by (4.4)–(4.6) for Theorem 4.1 are satisfied. Hence, is locally asymptotically stable.

Moreover, using and , it can also be verified that the conditions given by (4.11)–(4.13) for Theorem 4.2 are satisfied. Hence, is globally asymptotically stable.

Figure 2(a) depicts the numerical approximation to the time series for the populations, and Figure 2(b) depicts the numerical approximation to the time series for the pollution levels. An initial condition of was used for the simulation. Both time series were achieved using the MATLAB solver ode15s.

*Example 5.2. *If we increase the amount of pollution being input to the system, reduce the birthrate of the plant population, and rerun the simulations using the set of parameter values in Table 2.The nontrivial equilibrium, , of the model exists and is

Using , and , the conditions given by (4.4)–(4.6) for Theorem 4.1 are satisfied, and therefore is proven to be locally asymptotically stable.

Similarly, using and , the conditions given by (4.11)–(4.13) for Theorem 4.2 are satisfied and is globally asymptotically stable.

*Figure 3(a) depicts the numerical approximation to the time series for the populations, and Figure 3(b) depicts the numerical approximation to the time series for the pollution levels. An initial condition of was used for the simulation. Both time series were achieved using the MATLAB solver ode15s. Notice that the plant population approaches the equilibrium much slower in Example 5.2 than in Example 5.1. Also, it appears as though with the parameter settings used in Example 5.2 that the pollution has a greater impact on the animal population than the plant population; the magnitude of the decrease in size of the animal population at the equilibrium between the two examples is greater than the decrease seen in the plant population. It is also interesting to note that the pollution values are two degrees of magnitude larger in the second example than in the first example. The increase in the concentration of the pollution levels is of the same level of magnitude as the increase in the rate of pollution input into the system. It was found that if the simulation is rerun with , then the second global stability criterion given in (4.12) fails; however, numerical solution remains similar to that computed in Example 5.2. This implies that the stability criteria are sufficient but not necessary.*

*6. Conclusions*

*6. Conclusions**In this paper a new system of differential equations designed to study the effects of pollution on two species—one plant and one animal—in a closed environment was developed. The model was analysed using standard methods for nonlinear systems of differential equations. The problem was chosen as much work in the field has been done in the past; however, simplifications have always been introduced in order to help make the problem mathematically tractable. Our work strives to incorporate another level of biological realism to the model by fully accounting for the movement of pollution in the system using a mass-balance approach. In addition, the work was designed to include practical parameters that would be measurable in a real-world context.*

*It is clear from Examples 5.1 and 5.2 that the nontrivial solution for our system is stable for at least some range of parameter values and there is also a parameter range in which the stability criteria are satisfied. By analysing the model using the Lyapunov direct method we were able to derive criteria for the stability of the system. For both the local and global stability, three sufficient conditions for stability were determined. For local stability, these are given in (4.4), (4.5), and (4.6). The three criteria sufficient for establishing global stability are given in (4.11), (4.12), and (4.13). We can interpret the biological meaning of these criteria to gain further insight into our model.*

*Upon analysing the parameters included in stability criteria it was found that the system is locally stable provided that(1)pollution degrades or is removed from the environment at a fast enough rate, (2)the plant population experiences net growth at a fast enough rate, (3)(a) the animal population relieves itself of pollution at a fast enough rate(b)the animal carrying capacity is large enough,(c)the animal population births at a fast enough rate. *

*The system is globally stable provided that (1)pollution degrades or is removed from the environment at a fast enough rate,(2)the plant population experiences net growth at a fast enough rate, (3)the animal population is able to relieve itself of its pollution burden at a great enough.*

*Further work that could be done with the model includes parameterizing the model using parameters for a known situation and comparing the predictions of the model to the real-world phenomena.*

*Further extension that could be made based upon the model includes incorporating features such as age structure, delay, or diffusion.*

*Appendices*

*Appendices**A. Proof of Theorem 4.1*

*A. Proof of Theorem 4.1**We will work with the linearized form of the model, which can be written as
where represents the absolute value from the row and column of the Jacobian matrix evaluated at .*

*Consider the positive definite function*

*Differentiating with respect to , we get*

*Using the notation
and some algebraic manipulation, we can rewrite the system as *

*Then, factoring the equations, we obtain
*

*Sufficient conditions for to be negative definite are*

*Examining (A.7), we find
Considering (A.8), we find:
Considering (A.9), we find an additional condition for choosing our coefficient :*

*Note that, if both of the conditions relating to are satisfied, then*

*Looking at (A.10), we find the first general condition for local stability: *

*Equation (A.11) gives the second condition for stability:
If we let , this reduces to*

*Finally, (A.12) gives the third and final condition for local stability:*