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

Authors: Victoria Maystruk; Kenzu Abdella
Journal: ISRN Applied Mathematics
 (2011)

---

## 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.


---

## Body

## 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—N
(
t
)
,
B
(
t
)
,
C
N
(
t
)
,
C
B
(
t
), and C
E
(
t
). N
(
t
) represents the animal species density at time t
,
B
(
t
) represents the plant species density at time t
,
C
N
(
t
) represents the concentration of the toxicant in an individual animal at time t
,
C
B
(
t
) represents the concentration of the toxicant in an individual plant at time t, and C
E
(
t
) represents the concentration of the toxicant in the environment at time t.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,b
(
t
), and death rate, d
(
t
), of the populations. In general, we can think of population growth as (2.1)
d
N
d
t
=
(
b
(
t
)
-
d
(
t
)
)
N
,
If we let r
0
=
b
0
-
d
0 where b
0 is the intrinsic birth rate and d
0 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(2.2)
b
(
t
)
=
b
0
-
r
0
N
K
.
Similarly, the general death rate can be thought of as (2.3)
d
(
t
)
=
d
0
.
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. Letb
B
0 denote the intrinsic birth rate of the plant species, r
B
0 the intrinsic growth rate of the plant species, K
B
0 the intrinsic carrying capacity of the plant species, and K
B
C
E 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(2.4)
b
B
(
B
,
C
E
)
=
b
B
0
-
r
B
0
B
K
B
0
-
K
B
C
E
C
E
.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. Ifd
B
0 is the intrinsic death rate of the plant species, let d
B
N and d
B
C
B 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(2.5)
d
B
(
N
,
C
B
)
=
d
B
0
+
d
B
N
N
+
d
B
C
B
C
B
.Usingr
B
0
=
b
B
0
-
d
b
0, the equation governing the plant population is(2.6)
d
B
d
t
=
(
b
B
(
B
,
C
E
)
-
d
B
(
N
,
C
B
)
)
B
,
(2.7)
d
B
d
t
=
(
r
B
0
-
r
B
0
B
K
B
0
-
K
B
C
E
C
E
-
d
B
N
N
-
d
B
C
B
C
B
)
B
.
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. Letb
N
0
,
r
N
0, and K
B
0 represent the intrinsic birth rate, growth rate, and carrying capacity of the animal species, respectively. Additionally, let K
N
B denote the rate the carrying capacity increases relative to the plant biomass and K
N
C
E 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 (2.8)
b
N
(
N
,
B
,
C
E
)
=
b
N
0
-
r
N
0
N
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
.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 letd
N
0 represent the intrinsic death rate of the animal species, d
N
B the rate the animals death rate decreases relative to the biomass of the plant, and d
N
C
N 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(2.9)
d
N
(
B
,
C
N
)
=
d
N
0
-
d
N
B
B
+
d
N
C
N
C
N
.
Clearly, using r
N
0
=
b
N
0
-
d
N
0, we can describe the dynamics of the animal population as(2.10)
d
N
d
t
=
(
b
N
(
N
,
B
,
C
E
)
-
d
N
(
B
,
C
N
)
)
N
,
(2.11)
d
N
d
t
=
(
r
N
0
-
r
N
0
N
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
+
d
N
B
B
-
d
N
C
N
C
N
)
N
.
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 usem
n and m
b to represent the average per capita mass of the animal and plant populations, respectively, and we use m
e to represent the mass of the environment. Then m
b
B
C
B represents the total mass of toxicant in the plant population, m
n
N
C
N represents the total mass of toxicant in the animal population, and m
e
C
E represents the total quantity of toxicant free in the environment. m
e is assumed to be sufficiently large and near constant, and hence variations in m
e do not need to be modeled. The flow of pollution in the system can be visualized in Figure 1 and is described next.Figure 1
Flowchart of the pollution movement through the ecosystem.We assume that pollution first enters the environment at a rateU
(
t
). This pollution naturally dissipates at a rate relative to the amount of pollution present in the environment h
m
e
C
E. 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 g
1
m
b
C
E
B. The loss of environmental pollution due to absorption into the animal population is given as k
2
m
n
C
E
N. Pollution reenters the environment, as part of the pollution cycle as it is egested by the animal population. The egestion occurs at the rate k
3
m
n
C
N
N. 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 m
n
d
N
(
B
,
C
N
)
C
N
N and from dead plants at the rate m
b
d
B
(
N
,
C
B
)
C
B
B. Using these rates pertaining to the amount of pollution in the environment, we get(2.12)
d
m
e
C
E
d
t
=
U
(
t
)
-
h
m
e
C
E
-
g
1
m
b
C
E
B
-
k
2
m
n
C
E
N
+
k
3
m
n
C
N
N
+
m
n
d
N
(
B
,
C
N
)
C
N
N
+
m
b
d
B
(
N
,
C
B
)
C
B
B
;
then, dividing by m
e and using the notation u
(
t
)
=
U
(
t
)
/
m
e, we have
(2.13)
d
C
E
d
t
=
u
(
t
)
-
h
C
E
-
g
1
m
b
m
e
C
E
B
-
k
2
m
n
m
e
C
E
N
+
k
3
m
n
m
e
C
N
N
+
m
n
m
e
(
d
N
0
-
d
N
B
B
+
d
N
C
N
C
N
)
C
N
N
+
m
b
m
e
(
b
N
0
-
r
N
0
N
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
)
C
B
B
.
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 rateg
1
m
b
C
E
B. Pollution inside the plant biomass is naturally metabolized, leaving the system at a rate g
2
m
b
C
B
B. Pollution returns to the environment due to dead biomass at a rate d
B
(
N
,
C
B
)
m
b
C
B
B. 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 k
1
m
b
m
n
C
B
B
N. So we get (2.14)
d
C
B
m
b
B
d
t
=
g
1
m
b
C
E
B
-
g
2
m
b
C
B
B
-
d
B
(
N
,
C
B
)
m
b
C
B
B
-
k
1
m
b
m
n
C
B
B
N
;
dividing by m
b, we have
(2.15)
d
C
B
B
d
t
=
g
1
C
E
B
-
g
2
C
B
B
-
k
1
m
n
C
B
B
N
-
d
B
(
N
,
C
B
)
C
B
B
,
but clearly
(2.16)
d
C
B
B
d
t
=
C
B
d
B
d
t
+
B
d
C
B
d
t
.
Substituting in (2.6) and (2.15), using (2.17)
b
B
(
B
,
C
E
)
=
b
B
0
+
r
B
0
B
K
B
0
-
K
B
C
E
C
E
,
and factoring, we have (2.18)
d
C
B
d
t
=
g
1
C
E
-
k
1
m
n
C
B
N
-
(
g
2
+
b
B
0
-
r
B
0
B
K
B
0
-
K
B
C
E
C
E
)
C
B
.
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 ratek
2
m
n
C
E
N and by the ingestion of plant biomass containing pollution at a rate k
1
m
n
m
b
C
B
B
N. Pollution is released from the animal population as a result of being egested, which occurs at a loss rate of k
3
m
n
C
N
N, metabolized, which occurs at a loss rate of k
4
m
n
C
N
N, and released from dead animals, which occurs at a rate of d
N
(
B
,
C
N
)
m
n
C
N
N. Hence, we obtain (2.19)
d
C
N
m
n
N
d
t
=
k
1
m
n
m
b
C
B
B
N
+
k
2
m
n
C
E
N
-
k
3
m
n
C
N
N
-
k
4
m
n
C
N
N
-
d
N
(
B
,
C
N
)
m
n
C
N
N
.First, we divide bym
n to get(2.20)
d
C
N
N
d
t
=
k
1
m
b
C
B
B
N
+
k
2
C
E
N
-
k
3
C
N
N
-
k
4
C
N
N
-
d
N
(
B
,
C
N
)
C
N
N
.
Then, we use the fact that (2.21)
d
C
N
N
d
t
=
C
N
d
N
d
t
+
N
d
C
N
d
t
.
Substituting in (2.10) and (2.20), using b
N
(
N
,
B
,
C
E
)
=
b
N
0
+
r
N
0
N
/
(
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
), and factoring we determine that (2.22)
d
C
N
d
t
=
k
1
m
b
C
B
B
+
k
2
C
E
-
(
k
3
+
k
4
+
b
N
0
-
r
N
0
N
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
)
C
N
.
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:(2.23)
d
N
d
t
=
(
r
N
0
-
r
N
0
N
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
+
d
N
B
B
-
d
N
C
N
C
N
)
N
,
d
B
d
t
=
(
r
B
0
-
r
B
0
B
K
B
0
-
K
B
C
E
C
E
-
d
B
N
N
-
d
B
C
B
C
B
)
B
,
d
C
N
d
t
=
k
1
m
b
C
B
B
+
k
2
C
E
-
(
k
3
+
k
4
+
b
N
0
-
r
N
0
N
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
)
C
N
,
d
C
B
d
t
=
g
1
C
E
-
k
1
m
n
C
B
N
-
(
g
2
+
b
B
0
-
r
B
0
B
K
B
0
-
K
B
C
E
C
E
)
C
B
,
d
C
E
d
t
=
u
(
t
)
+
m
b
m
e
(
d
B
0
+
d
B
N
N
+
d
B
C
B
C
B
)
C
B
B
+
m
n
m
e
(
k
3
+
d
N
0
-
d
N
B
B
+
d
N
C
N
C
N
)
C
N
N
-
m
b
m
e
g
1
C
E
B
-
m
n
m
e
k
2
C
E
N
-
h
C
E
.


## 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. Letb
B
0 denote the intrinsic birth rate of the plant species, r
B
0 the intrinsic growth rate of the plant species, K
B
0 the intrinsic carrying capacity of the plant species, and K
B
C
E 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(2.4)
b
B
(
B
,
C
E
)
=
b
B
0
-
r
B
0
B
K
B
0
-
K
B
C
E
C
E
.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. Ifd
B
0 is the intrinsic death rate of the plant species, let d
B
N and d
B
C
B 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(2.5)
d
B
(
N
,
C
B
)
=
d
B
0
+
d
B
N
N
+
d
B
C
B
C
B
.Usingr
B
0
=
b
B
0
-
d
b
0, the equation governing the plant population is(2.6)
d
B
d
t
=
(
b
B
(
B
,
C
E
)
-
d
B
(
N
,
C
B
)
)
B
,
(2.7)
d
B
d
t
=
(
r
B
0
-
r
B
0
B
K
B
0
-
K
B
C
E
C
E
-
d
B
N
N
-
d
B
C
B
C
B
)
B
.
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. Letb
N
0
,
r
N
0, and K
B
0 represent the intrinsic birth rate, growth rate, and carrying capacity of the animal species, respectively. Additionally, let K
N
B denote the rate the carrying capacity increases relative to the plant biomass and K
N
C
E 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 (2.8)
b
N
(
N
,
B
,
C
E
)
=
b
N
0
-
r
N
0
N
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
.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 letd
N
0 represent the intrinsic death rate of the animal species, d
N
B the rate the animals death rate decreases relative to the biomass of the plant, and d
N
C
N 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(2.9)
d
N
(
B
,
C
N
)
=
d
N
0
-
d
N
B
B
+
d
N
C
N
C
N
.
Clearly, using r
N
0
=
b
N
0
-
d
N
0, we can describe the dynamics of the animal population as(2.10)
d
N
d
t
=
(
b
N
(
N
,
B
,
C
E
)
-
d
N
(
B
,
C
N
)
)
N
,
(2.11)
d
N
d
t
=
(
r
N
0
-
r
N
0
N
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
+
d
N
B
B
-
d
N
C
N
C
N
)
N
.
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 usem
n and m
b to represent the average per capita mass of the animal and plant populations, respectively, and we use m
e to represent the mass of the environment. Then m
b
B
C
B represents the total mass of toxicant in the plant population, m
n
N
C
N represents the total mass of toxicant in the animal population, and m
e
C
E represents the total quantity of toxicant free in the environment. m
e is assumed to be sufficiently large and near constant, and hence variations in m
e do not need to be modeled. The flow of pollution in the system can be visualized in Figure 1 and is described next.Figure 1
Flowchart of the pollution movement through the ecosystem.We assume that pollution first enters the environment at a rateU
(
t
). This pollution naturally dissipates at a rate relative to the amount of pollution present in the environment h
m
e
C
E. 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 g
1
m
b
C
E
B. The loss of environmental pollution due to absorption into the animal population is given as k
2
m
n
C
E
N. Pollution reenters the environment, as part of the pollution cycle as it is egested by the animal population. The egestion occurs at the rate k
3
m
n
C
N
N. 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 m
n
d
N
(
B
,
C
N
)
C
N
N and from dead plants at the rate m
b
d
B
(
N
,
C
B
)
C
B
B. Using these rates pertaining to the amount of pollution in the environment, we get(2.12)
d
m
e
C
E
d
t
=
U
(
t
)
-
h
m
e
C
E
-
g
1
m
b
C
E
B
-
k
2
m
n
C
E
N
+
k
3
m
n
C
N
N
+
m
n
d
N
(
B
,
C
N
)
C
N
N
+
m
b
d
B
(
N
,
C
B
)
C
B
B
;
then, dividing by m
e and using the notation u
(
t
)
=
U
(
t
)
/
m
e, we have
(2.13)
d
C
E
d
t
=
u
(
t
)
-
h
C
E
-
g
1
m
b
m
e
C
E
B
-
k
2
m
n
m
e
C
E
N
+
k
3
m
n
m
e
C
N
N
+
m
n
m
e
(
d
N
0
-
d
N
B
B
+
d
N
C
N
C
N
)
C
N
N
+
m
b
m
e
(
b
N
0
-
r
N
0
N
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
)
C
B
B
.
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 rateg
1
m
b
C
E
B. Pollution inside the plant biomass is naturally metabolized, leaving the system at a rate g
2
m
b
C
B
B. Pollution returns to the environment due to dead biomass at a rate d
B
(
N
,
C
B
)
m
b
C
B
B. 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 k
1
m
b
m
n
C
B
B
N. So we get (2.14)
d
C
B
m
b
B
d
t
=
g
1
m
b
C
E
B
-
g
2
m
b
C
B
B
-
d
B
(
N
,
C
B
)
m
b
C
B
B
-
k
1
m
b
m
n
C
B
B
N
;
dividing by m
b, we have
(2.15)
d
C
B
B
d
t
=
g
1
C
E
B
-
g
2
C
B
B
-
k
1
m
n
C
B
B
N
-
d
B
(
N
,
C
B
)
C
B
B
,
but clearly
(2.16)
d
C
B
B
d
t
=
C
B
d
B
d
t
+
B
d
C
B
d
t
.
Substituting in (2.6) and (2.15), using (2.17)
b
B
(
B
,
C
E
)
=
b
B
0
+
r
B
0
B
K
B
0
-
K
B
C
E
C
E
,
and factoring, we have (2.18)
d
C
B
d
t
=
g
1
C
E
-
k
1
m
n
C
B
N
-
(
g
2
+
b
B
0
-
r
B
0
B
K
B
0
-
K
B
C
E
C
E
)
C
B
.
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 ratek
2
m
n
C
E
N and by the ingestion of plant biomass containing pollution at a rate k
1
m
n
m
b
C
B
B
N. Pollution is released from the animal population as a result of being egested, which occurs at a loss rate of k
3
m
n
C
N
N, metabolized, which occurs at a loss rate of k
4
m
n
C
N
N, and released from dead animals, which occurs at a rate of d
N
(
B
,
C
N
)
m
n
C
N
N. Hence, we obtain (2.19)
d
C
N
m
n
N
d
t
=
k
1
m
n
m
b
C
B
B
N
+
k
2
m
n
C
E
N
-
k
3
m
n
C
N
N
-
k
4
m
n
C
N
N
-
d
N
(
B
,
C
N
)
m
n
C
N
N
.First, we divide bym
n to get(2.20)
d
C
N
N
d
t
=
k
1
m
b
C
B
B
N
+
k
2
C
E
N
-
k
3
C
N
N
-
k
4
C
N
N
-
d
N
(
B
,
C
N
)
C
N
N
.
Then, we use the fact that (2.21)
d
C
N
N
d
t
=
C
N
d
N
d
t
+
N
d
C
N
d
t
.
Substituting in (2.10) and (2.20), using b
N
(
N
,
B
,
C
E
)
=
b
N
0
+
r
N
0
N
/
(
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
), and factoring we determine that (2.22)
d
C
N
d
t
=
k
1
m
b
C
B
B
+
k
2
C
E
-
(
k
3
+
k
4
+
b
N
0
-
r
N
0
N
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
)
C
N
.
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:(2.23)
d
N
d
t
=
(
r
N
0
-
r
N
0
N
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
+
d
N
B
B
-
d
N
C
N
C
N
)
N
,
d
B
d
t
=
(
r
B
0
-
r
B
0
B
K
B
0
-
K
B
C
E
C
E
-
d
B
N
N
-
d
B
C
B
C
B
)
B
,
d
C
N
d
t
=
k
1
m
b
C
B
B
+
k
2
C
E
-
(
k
3
+
k
4
+
b
N
0
-
r
N
0
N
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
)
C
N
,
d
C
B
d
t
=
g
1
C
E
-
k
1
m
n
C
B
N
-
(
g
2
+
b
B
0
-
r
B
0
B
K
B
0
-
K
B
C
E
C
E
)
C
B
,
d
C
E
d
t
=
u
(
t
)
+
m
b
m
e
(
d
B
0
+
d
B
N
N
+
d
B
C
B
C
B
)
C
B
B
+
m
n
m
e
(
k
3
+
d
N
0
-
d
N
B
B
+
d
N
C
N
C
N
)
C
N
N
-
m
b
m
e
g
1
C
E
B
-
m
n
m
e
k
2
C
E
N
-
h
C
E
.


## 3. Equilibrium Analysis

Our model has four nonnegative equilibria:(3.1)
E
0
=
(
0,0
,
C
N
0
,
C
B
0
,
C
E
0
)
,
E
1
=
(
0
,
B
1
,
C
N
1
,
C
B
1
,
C
E
1
)
,
E
2
=
(
N
2
,
0
,
C
N
2
,
C
B
2
,
C
E
2
)
,
E
*
=
(
N
*
,
B
*
,
C
N
*
,
C
B
*
,
C
E
*
)
.
Since we derived the equations for d
C
B
/
d
t and d
C
N
/
d
t under the assumption that both N and B are nonzero, then the solution we are chiefly concerned with is E
*.


### 3.1. Existence ofE
*
=
(
N
*
,
B
*
,
C
N
*
,
C
B
*
,
C
E
*
)

ForE
*
=
(
N
*
,
B
*
,
C
N
*
,
C
B
*
,
C
E
*
)
,
N
*
,
B
*
,
C
N
*
,
C
B
*, and C
E
* are the positive solutions of, (3.2)
0
=
r
N
0
-
r
N
0
N
*
K
N
0
+
K
N
B
B
*
-
K
N
C
E
C
E
*
+
d
N
B
B
*
-
d
N
C
N
C
N
*
,
(3.3)
0
=
r
B
0
-
r
B
0
B
*
K
B
0
-
K
B
C
E
C
E
*
-
d
B
N
N
*
-
d
B
C
B
C
B
*
,
(3.4)
0
=
k
1
m
b
C
B
*
B
*
+
k
2
C
E
*
-
(
k
3
+
k
4
+
b
N
0
-
r
N
0
N
*
K
N
0
+
K
N
B
B
*
-
K
N
C
E
C
E
*
)
C
N
*
,
(3.5)
0
=
g
1
C
E
*
-
k
1
m
n
C
B
*
N
*
-
(
g
2
+
b
B
0
-
r
B
0
B
*
K
B
0
-
K
B
C
E
C
E
*
)
C
B
*
-
m
b
m
e
g
1
C
E
*
B
*
,
(3.6)
0
=
u
(
t
)
+
m
b
m
e
(
d
B
0
+
d
B
N
N
*
+
d
B
C
B
C
B
*
)
C
B
*
B
*
-
m
b
m
e
g
1
C
E
*
B
*
+
m
n
m
e
(
k
3
+
d
N
0
-
d
N
B
B
*
+
d
N
C
N
C
N
*
)
C
N
N
*
-
m
n
m
e
k
2
C
E
*
N
*
-
h
C
E
*
.Recall thatr
B
0
=
b
B
0
-
d
B
0, and so, by (3.3), (3.7)
b
B
0
-
r
B
0
B
*
K
B
0
-
K
B
C
E
C
E
*
=
d
b
0
+
d
B
N
N
*
+
d
B
C
B
C
B
*
.
Substituting this into (3.5), we have(3.8)
0
=
d
B
C
B
(
C
B
*
)
2
+
(
k
1
m
n
N
*
+
g
2
+
d
b
0
+
d
B
N
N
*
)
C
B
*
-
g
1
C
E
*
.
Hence, (3.9)
C
B
*
=
-
β
1
(
N
*
)
+
β
1
(
N
*
)
2
+
4
g
1
d
B
C
B
C
E
*
2
d
B
C
B
,
where,
(3.10)
β
1
(
N
*
)
=
k
1
m
n
N
*
+
d
B
N
N
*
+
g
2
+
d
B
0
.
Let(3.11)
f
1
(
N
*
,
C
E
*
)
=
-
β
1
(
N
*
)
+
β
1
(
N
*
)
2
+
4
g
1
d
B
C
B
C
E
*
2
d
B
C
B
.
Consider (3.3). Since C
B
=
f
1
(
N
*
,
C
E
*
), (3.12)
B
*
=
(
K
B
0
-
K
B
C
E
C
E
*
)
(
r
B
0
-
d
B
N
N
*
-
d
B
C
B
f
1
(
N
*
,
C
E
*
)
)
r
B
0
.
Let (3.13)
f
2
(
N
*
,
C
E
*
)
=
(
K
B
0
-
K
B
C
E
C
E
*
)
(
r
B
0
-
d
B
N
N
*
-
d
B
C
B
f
1
(
N
*
,
C
E
*
)
)
r
B
0
.
Next, recall that r
N
0
=
b
N
0
-
d
N
0, and so, by (3.2),(3.14)
b
N
0
-
r
N
0
N
*
K
N
0
+
K
N
B
B
*
-
K
N
C
E
C
E
*
=
d
N
0
-
d
N
B
B
*
+
d
N
C
N
C
N
*
.
Substituting this into (3.4), we have(3.15)
0
=
d
N
C
N
(
C
N
*
)
2
+
(
k
3
+
k
4
+
d
N
0
-
d
N
B
B
*
)
C
N
*
-
(
k
1
m
b
C
B
*
B
*
+
k
2
C
E
*
)
.
Using C
B
*
=
f
1
(
N
*
,
C
E
*
) and B
*
=
f
1
(
N
*
,
C
E
*
), we find
(3.16)
C
N
*
=
-
β
2
(
N
*
,
C
E
*
)
+
β
2
(
N
*
,
C
E
*
)
2
+
β
3
(
N
*
,
C
E
*
)
2
d
N
C
N
,
where
(3.17)
β
2
(
N
*
,
C
E
*
)
=
k
3
+
k
4
+
d
N
0
-
d
N
B
f
2
(
N
*
,
C
E
*
)
,
β
3
(
N
*
,
C
E
*
)
=
4
d
N
C
N
(
k
1
m
b
f
1
(
N
*
,
C
E
*
)
f
2
(
N
*
,
C
E
*
)
+
k
2
C
E
*
)
.
Let (3.18)
f
3
(
N
*
,
C
E
*
)
=
-
β
2
(
N
*
,
C
E
*
)
+
β
2
(
N
*
,
C
E
*
)
2
+
β
3
(
N
*
,
C
E
*
)
2
d
N
C
N
.
Next, if we consider (3.2), we have (3.19)
0
=
r
N
0
-
r
N
0
N
*
K
N
0
+
K
N
B
B
*
-
K
N
C
E
C
E
*
+
d
N
B
B
*
-
d
N
C
N
C
N
*
.
Using B
*
=
f
2
(
N
*
,
C
E
*
) and C
N
*
=
f
3
(
N
*
,
C
E
*
), we find (3.20)
N
*
=
β
4
(
N
*
,
C
E
*
)
(
r
N
0
+
d
N
B
f
2
(
N
*
,
C
E
*
)
-
d
N
C
N
f
3
(
N
*
,
C
E
*
)
)
r
N
0
,
where (3.21)
β
4
(
N
*
,
C
E
*
)
=
K
N
0
+
K
N
B
f
2
(
N
*
,
C
E
*
)
-
K
N
C
E
C
E
*
.
Finally, we consider (3.6) and (3.22)
0
=
u
(
t
)
+
m
b
m
e
(
d
B
0
+
d
B
N
N
*
+
d
B
C
B
C
B
*
)
C
B
*
B
*
-
m
b
m
e
g
1
C
E
*
B
*
+
m
n
m
e
(
k
3
+
d
N
0
-
d
N
B
B
*
+
d
N
C
N
C
N
*
)
C
N
N
*
-
m
n
m
e
k
2
C
E
*
N
*
-
h
C
E
*
.
Substituting in C
B
*
=
f
1
(
N
*
,
C
E
*
), B
*
=
f
2
(
N
*
,
C
E
*
), and C
N
*
=
f
3
(
N
*
,
C
E
*
), we find(3.23)
C
E
*
=
m
e
u
+
m
b
(
d
B
0
+
d
B
N
N
*
+
d
B
C
B
f
1
)
f
1
f
2
m
b
g
1
f
2
+
m
n
k
2
N
*
+
m
e
h
+
m
n
(
k
3
+
d
N
0
-
d
N
B
f
2
+
d
N
C
N
f
3
)
f
3
N
*
m
b
g
1
f
2
+
m
n
k
2
N
*
+
m
e
h
.Thus, the nontrivial solution to the system is given by(3.24)
N
*
=
β
4
(
N
*
,
C
E
*
)
(
r
N
0
-
d
N
C
N
f
3
(
N
*
,
C
E
*
)
+
d
N
B
f
2
(
N
*
,
C
E
*
)
)
r
N
0
,
B
*
=
f
2
(
N
*
,
C
E
*
)
,
C
N
*
=
f
3
(
N
*
,
C
E
*
)
,
C
B
*
=
f
1
(
N
*
,
C
E
*
)
,
C
E
*
=
m
e
u
+
m
b
β
6
(
N
*
,
C
E
*
)
f
1
(
N
*
,
C
E
*
)
f
2
(
N
*
,
C
E
*
)
β
5
(
N
*
,
C
E
*
)
+
m
n
β
7
(
N
*
,
C
E
*
)
f
3
(
N
*
,
C
E
*
)
N
*
β
5
(
N
*
,
C
E
*
)
,
where (3.25)
f
1
(
N
*
,
C
E
*
)
=
-
β
1
(
N
*
)
+
β
1
(
N
*
)
2
+
4
g
1
d
B
C
B
C
E
*
2
d
B
C
B
,
f
2
(
N
*
,
C
E
*
)
=
(
K
B
0
-
K
B
C
E
C
E
*
)
(
r
B
0
-
d
B
N
N
*
+
d
B
C
B
f
1
(
N
*
,
C
E
*
)
)
r
B
0
,
f
3
(
N
*
,
C
E
*
)
=
-
β
2
(
N
*
,
C
E
*
)
+
β
2
(
N
*
,
C
E
*
)
2
+
β
3
(
N
*
,
C
E
*
)
2
d
N
C
N
,
where (3.26)
β
1
(
N
*
)
=
k
1
m
n
N
*
+
d
B
N
N
*
+
g
2
+
d
B
0
,
β
2
(
N
*
,
C
E
*
)
=
k
3
+
k
4
+
d
N
0
-
d
N
B
f
2
(
N
*
,
C
E
*
)
,
β
3
(
N
*
,
C
E
*
)
=
4
d
N
C
N
(
k
1
m
b
f
1
(
N
*
,
C
E
*
)
f
2
(
N
*
,
C
E
*
)
+
k
2
C
E
*
)
,
β
4
(
N
*
,
C
E
*
)
=
K
N
0
+
K
N
B
f
2
(
N
*
,
C
E
*
)
-
K
N
C
E
C
E
*
,
β
5
(
N
*
,
C
E
*
)
=
m
b
g
1
f
2
(
N
*
,
C
E
*
)
+
m
n
k
2
N
*
+
m
e
h
,
β
6
(
N
*
,
C
E
*
)
=
d
B
0
+
d
B
N
N
*
+
d
B
C
B
f
1
(
N
*
,
C
E
*
)
,
β
7
(
N
*
,
C
E
*
)
=
k
3
+
d
N
0
-
d
N
B
f
2
(
N
*
,
C
E
*
)
+
d
N
C
N
f
3
(
N
*
,
C
E
*
)
.
Provided that N
* and C
E
* exist, we can find the corresponding C
B
*
,
B
*, and C
N
*.


## 3.1. Existence ofE
*
=
(
N
*
,
B
*
,
C
N
*
,
C
B
*
,
C
E
*
)

ForE
*
=
(
N
*
,
B
*
,
C
N
*
,
C
B
*
,
C
E
*
)
,
N
*
,
B
*
,
C
N
*
,
C
B
*, and C
E
* are the positive solutions of, (3.2)
0
=
r
N
0
-
r
N
0
N
*
K
N
0
+
K
N
B
B
*
-
K
N
C
E
C
E
*
+
d
N
B
B
*
-
d
N
C
N
C
N
*
,
(3.3)
0
=
r
B
0
-
r
B
0
B
*
K
B
0
-
K
B
C
E
C
E
*
-
d
B
N
N
*
-
d
B
C
B
C
B
*
,
(3.4)
0
=
k
1
m
b
C
B
*
B
*
+
k
2
C
E
*
-
(
k
3
+
k
4
+
b
N
0
-
r
N
0
N
*
K
N
0
+
K
N
B
B
*
-
K
N
C
E
C
E
*
)
C
N
*
,
(3.5)
0
=
g
1
C
E
*
-
k
1
m
n
C
B
*
N
*
-
(
g
2
+
b
B
0
-
r
B
0
B
*
K
B
0
-
K
B
C
E
C
E
*
)
C
B
*
-
m
b
m
e
g
1
C
E
*
B
*
,
(3.6)
0
=
u
(
t
)
+
m
b
m
e
(
d
B
0
+
d
B
N
N
*
+
d
B
C
B
C
B
*
)
C
B
*
B
*
-
m
b
m
e
g
1
C
E
*
B
*
+
m
n
m
e
(
k
3
+
d
N
0
-
d
N
B
B
*
+
d
N
C
N
C
N
*
)
C
N
N
*
-
m
n
m
e
k
2
C
E
*
N
*
-
h
C
E
*
.Recall thatr
B
0
=
b
B
0
-
d
B
0, and so, by (3.3), (3.7)
b
B
0
-
r
B
0
B
*
K
B
0
-
K
B
C
E
C
E
*
=
d
b
0
+
d
B
N
N
*
+
d
B
C
B
C
B
*
.
Substituting this into (3.5), we have(3.8)
0
=
d
B
C
B
(
C
B
*
)
2
+
(
k
1
m
n
N
*
+
g
2
+
d
b
0
+
d
B
N
N
*
)
C
B
*
-
g
1
C
E
*
.
Hence, (3.9)
C
B
*
=
-
β
1
(
N
*
)
+
β
1
(
N
*
)
2
+
4
g
1
d
B
C
B
C
E
*
2
d
B
C
B
,
where,
(3.10)
β
1
(
N
*
)
=
k
1
m
n
N
*
+
d
B
N
N
*
+
g
2
+
d
B
0
.
Let(3.11)
f
1
(
N
*
,
C
E
*
)
=
-
β
1
(
N
*
)
+
β
1
(
N
*
)
2
+
4
g
1
d
B
C
B
C
E
*
2
d
B
C
B
.
Consider (3.3). Since C
B
=
f
1
(
N
*
,
C
E
*
), (3.12)
B
*
=
(
K
B
0
-
K
B
C
E
C
E
*
)
(
r
B
0
-
d
B
N
N
*
-
d
B
C
B
f
1
(
N
*
,
C
E
*
)
)
r
B
0
.
Let (3.13)
f
2
(
N
*
,
C
E
*
)
=
(
K
B
0
-
K
B
C
E
C
E
*
)
(
r
B
0
-
d
B
N
N
*
-
d
B
C
B
f
1
(
N
*
,
C
E
*
)
)
r
B
0
.
Next, recall that r
N
0
=
b
N
0
-
d
N
0, and so, by (3.2),(3.14)
b
N
0
-
r
N
0
N
*
K
N
0
+
K
N
B
B
*
-
K
N
C
E
C
E
*
=
d
N
0
-
d
N
B
B
*
+
d
N
C
N
C
N
*
.
Substituting this into (3.4), we have(3.15)
0
=
d
N
C
N
(
C
N
*
)
2
+
(
k
3
+
k
4
+
d
N
0
-
d
N
B
B
*
)
C
N
*
-
(
k
1
m
b
C
B
*
B
*
+
k
2
C
E
*
)
.
Using C
B
*
=
f
1
(
N
*
,
C
E
*
) and B
*
=
f
1
(
N
*
,
C
E
*
), we find
(3.16)
C
N
*
=
-
β
2
(
N
*
,
C
E
*
)
+
β
2
(
N
*
,
C
E
*
)
2
+
β
3
(
N
*
,
C
E
*
)
2
d
N
C
N
,
where
(3.17)
β
2
(
N
*
,
C
E
*
)
=
k
3
+
k
4
+
d
N
0
-
d
N
B
f
2
(
N
*
,
C
E
*
)
,
β
3
(
N
*
,
C
E
*
)
=
4
d
N
C
N
(
k
1
m
b
f
1
(
N
*
,
C
E
*
)
f
2
(
N
*
,
C
E
*
)
+
k
2
C
E
*
)
.
Let (3.18)
f
3
(
N
*
,
C
E
*
)
=
-
β
2
(
N
*
,
C
E
*
)
+
β
2
(
N
*
,
C
E
*
)
2
+
β
3
(
N
*
,
C
E
*
)
2
d
N
C
N
.
Next, if we consider (3.2), we have (3.19)
0
=
r
N
0
-
r
N
0
N
*
K
N
0
+
K
N
B
B
*
-
K
N
C
E
C
E
*
+
d
N
B
B
*
-
d
N
C
N
C
N
*
.
Using B
*
=
f
2
(
N
*
,
C
E
*
) and C
N
*
=
f
3
(
N
*
,
C
E
*
), we find (3.20)
N
*
=
β
4
(
N
*
,
C
E
*
)
(
r
N
0
+
d
N
B
f
2
(
N
*
,
C
E
*
)
-
d
N
C
N
f
3
(
N
*
,
C
E
*
)
)
r
N
0
,
where (3.21)
β
4
(
N
*
,
C
E
*
)
=
K
N
0
+
K
N
B
f
2
(
N
*
,
C
E
*
)
-
K
N
C
E
C
E
*
.
Finally, we consider (3.6) and (3.22)
0
=
u
(
t
)
+
m
b
m
e
(
d
B
0
+
d
B
N
N
*
+
d
B
C
B
C
B
*
)
C
B
*
B
*
-
m
b
m
e
g
1
C
E
*
B
*
+
m
n
m
e
(
k
3
+
d
N
0
-
d
N
B
B
*
+
d
N
C
N
C
N
*
)
C
N
N
*
-
m
n
m
e
k
2
C
E
*
N
*
-
h
C
E
*
.
Substituting in C
B
*
=
f
1
(
N
*
,
C
E
*
), B
*
=
f
2
(
N
*
,
C
E
*
), and C
N
*
=
f
3
(
N
*
,
C
E
*
), we find(3.23)
C
E
*
=
m
e
u
+
m
b
(
d
B
0
+
d
B
N
N
*
+
d
B
C
B
f
1
)
f
1
f
2
m
b
g
1
f
2
+
m
n
k
2
N
*
+
m
e
h
+
m
n
(
k
3
+
d
N
0
-
d
N
B
f
2
+
d
N
C
N
f
3
)
f
3
N
*
m
b
g
1
f
2
+
m
n
k
2
N
*
+
m
e
h
.Thus, the nontrivial solution to the system is given by(3.24)
N
*
=
β
4
(
N
*
,
C
E
*
)
(
r
N
0
-
d
N
C
N
f
3
(
N
*
,
C
E
*
)
+
d
N
B
f
2
(
N
*
,
C
E
*
)
)
r
N
0
,
B
*
=
f
2
(
N
*
,
C
E
*
)
,
C
N
*
=
f
3
(
N
*
,
C
E
*
)
,
C
B
*
=
f
1
(
N
*
,
C
E
*
)
,
C
E
*
=
m
e
u
+
m
b
β
6
(
N
*
,
C
E
*
)
f
1
(
N
*
,
C
E
*
)
f
2
(
N
*
,
C
E
*
)
β
5
(
N
*
,
C
E
*
)
+
m
n
β
7
(
N
*
,
C
E
*
)
f
3
(
N
*
,
C
E
*
)
N
*
β
5
(
N
*
,
C
E
*
)
,
where (3.25)
f
1
(
N
*
,
C
E
*
)
=
-
β
1
(
N
*
)
+
β
1
(
N
*
)
2
+
4
g
1
d
B
C
B
C
E
*
2
d
B
C
B
,
f
2
(
N
*
,
C
E
*
)
=
(
K
B
0
-
K
B
C
E
C
E
*
)
(
r
B
0
-
d
B
N
N
*
+
d
B
C
B
f
1
(
N
*
,
C
E
*
)
)
r
B
0
,
f
3
(
N
*
,
C
E
*
)
=
-
β
2
(
N
*
,
C
E
*
)
+
β
2
(
N
*
,
C
E
*
)
2
+
β
3
(
N
*
,
C
E
*
)
2
d
N
C
N
,
where (3.26)
β
1
(
N
*
)
=
k
1
m
n
N
*
+
d
B
N
N
*
+
g
2
+
d
B
0
,
β
2
(
N
*
,
C
E
*
)
=
k
3
+
k
4
+
d
N
0
-
d
N
B
f
2
(
N
*
,
C
E
*
)
,
β
3
(
N
*
,
C
E
*
)
=
4
d
N
C
N
(
k
1
m
b
f
1
(
N
*
,
C
E
*
)
f
2
(
N
*
,
C
E
*
)
+
k
2
C
E
*
)
,
β
4
(
N
*
,
C
E
*
)
=
K
N
0
+
K
N
B
f
2
(
N
*
,
C
E
*
)
-
K
N
C
E
C
E
*
,
β
5
(
N
*
,
C
E
*
)
=
m
b
g
1
f
2
(
N
*
,
C
E
*
)
+
m
n
k
2
N
*
+
m
e
h
,
β
6
(
N
*
,
C
E
*
)
=
d
B
0
+
d
B
N
N
*
+
d
B
C
B
f
1
(
N
*
,
C
E
*
)
,
β
7
(
N
*
,
C
E
*
)
=
k
3
+
d
N
0
-
d
N
B
f
2
(
N
*
,
C
E
*
)
+
d
N
C
N
f
3
(
N
*
,
C
E
*
)
.
Provided that N
* and C
E
* exist, we can find the corresponding C
B
*
,
B
*, and C
N
*.


## 4. Stability Analysis

Theorem 4.1 (conditions for local stability).
Given(4.1)
c
3
<
1
4
r
N
0
(
g
2
+
k
1
m
n
N
*
+
b
B
(
B
*
,
C
E
*
)
)
K
N
(
B
*
,
C
E
*
)
(
k
1
m
n
C
B
*
)
2
,
c
1
=
c
3
r
B
0
C
B
*
d
B
C
B
K
B
(
C
E
*
)
,
c
2
<
c
3
4
(
k
3
+
k
4
+
b
N
(
N
*
,
B
*
,
C
E
*
)
)
Ω
(
N
*
,
B
*
,
C
N
*
,
C
B
*
,
C
E
*
)
,
where
(4.2)
Ω
(
N
*
,
B
*
,
C
N
*
,
C
B
*
,
C
E
*
)
=
min
⁡
⁡
(
C
B
*
d
B
C
B
(
r
B
0
K
B
(
C
E
*
)
δ
(
N
*
,
B
*
,
C
N
*
,
C
B
*
,
C
E
*
)
)
2
,
k
1
m
n
N
*
+
g
2
+
b
B
(
B
*
,
C
E
*
)
(
k
1
m
b
B
*
)
2
)
,
where
(4.3)
δ
(
N
*
,
B
*
,
C
N
*
,
C
B
*
,
C
E
*
)
=
k
1
m
b
C
B
*
-
K
N
B
r
N
0
N
*
C
N
*
K
N
(
B
*
,
C
E
*
)
2
.
if the following inequalities hold:
(1)
(4.4)
(
d
N
C
N
-
c
2
r
N
0
C
N
*
K
N
(
B
*
,
C
E
*
)
)
2
<
c
2
4
r
N
0
(
k
3
+
k
4
+
b
N
(
N
*
,
B
*
,
C
E
*
)
)
K
N
(
B
*
,
C
E
*
)
,
(2)
(4.5)
3
c
1
4
r
B
0
K
B
(
C
E
*
)
>
(
d
N
B
+
K
N
B
r
N
0
N
*
/
K
N
(
B
*
,
C
E
*
)
2
-
c
1
d
B
N
)
2
r
N
0
/
K
N
(
B
*
,
C
E
*
)
+
(
c
1
m
e
K
B
C
E
r
B
0
B
*
/
K
B
(
C
E
*
)
2
-
m
b
d
B
(
N
*
,
C
B
*
)
C
B
*
+
m
n
d
N
B
C
N
*
N
*
+
m
b
g
1
C
E
*
)
2
m
e
(
m
b
g
1
B
*
+
m
n
k
2
N
*
+
h
m
e
)
,
(3)
(4.6)
3
m
e
4
(
m
b
g
1
B
*
+
m
n
k
2
N
*
+
m
e
h
)
>

(
m
e
K
N
C
E
r
N
0
N
*
/
K
N
(
B
*
,
C
E
*
)
2

-

m
b
d
B
N
C
B
*
B
*

-

m
n
(
k
3
+
d
N
(
B
*
,
C
N
*
)
)
C
N
*

+

m
n
k
2
C
E
*
)
2
r
N
0
/
K
N
(
B
*
,
C
E
*
)
+

(
c
2
m
e
(
k
2

+

K
N
C
E
r
N
0
N
*
C
N
*
/
K
N
(
B
*
,
C
E
*
)
2
)

+
m
n
N
*
(
k
3

+

d
N
(
B
*
,
C
N
*
)

+

d
N
C
N
C
N
*
)
)
2
c
2
(
k
3
+
k
4
+
b
N
(
N
*
,
B
*
,
C
E
*
)
)
+
(
c
3
m
e
(
g
1
+
K
B
C
E
r
B
0
B
*
C
B
*
/
K
B
(
C
E
*
)
2
)
+
m
b
B
*
(
d
B
(
N
*
,
C
B
*
)
+
d
B
C
B
C
B
*
)
)
2
c
3
(
k
1
m
n
N
*
+
g
2
+
b
B
(
B
*
,
C
E
*
)
)
,
then E
* is a locally stable equilibrium of system (2.23).Theorem 4.2 (conditions for global stability).
Using the notation(4.7)
K
N
(
B
,
C
E
)
=
K
N
0
+
K
N
B
B
-
K
N
C
E
C
E
,
K
B
(
C
E
)
=
K
B
0
-
K
B
C
E
C
E
,
D
N
(
B
,
C
N
)
=
d
N
0
-
d
N
B
B
+
d
N
C
N
C
N
,
D
B
(
N
,
C
B
)
=
d
B
0
+
d
B
N
N
+
d
B
C
B
C
B
,
let K
N
(
B
,
C
E
)
,
K
B
(
C
E
)
,
D
N
(
B
,
C
N
), and D
B
(
N
,
C
B
) be such that
(4.8)
K
N
̲
=
K
N
0
-
K
N
C
E
≤
K
N
(
B
,
C
E
)
≤
K
N
0
+
K
N
B
K
B
0
=
K
N
¯
,
K
B
̲
=
K
B
0
-
K
B
C
E
≤
K
B
(
C
E
)
≤
K
B
0
=
K
B
¯
,
|
D
N
(
B
,
C
N
)
|
≤
d
N
0
+
d
N
C
N
=
D
N
¯
,
|
D
B
(
N
,
C
B
)
|
≤
d
B
0
+
d
B
N
N
M
+
d
B
C
B
=
D
B
¯
.
Then, given
(4.9)
0
<
c
2
<
r
N
0
(
g
2
+
k
1
m
n
N
*
)
4
(
k
1
m
n
)
2
K
N
¯
,
0
<
c
1
<
k
3
+
k
4
4
max
⁡
(
r
B
0
K
B
¯
(
k
1
m
b
+
r
N
0
N
*
C
N
*
η
1
¯
)
2
c
2
(
g
2
+
k
1
m
n
N
*
)
(
k
1
m
b
B
*
)
2
)
,
where
(4.10)
η
1
¯
=
K
N
B
(
K
N
0
+
B
*
K
N
B
-
K
N
C
E
)
2
,
η
2
¯
=
K
N
C
E
(
K
N
0
+
B
*
K
N
B
-
K
N
C
E
C
E
*
)
2
,
μ
¯
=
K
B
C
E
(
K
B
0
-
K
B
C
E
C
E
*
)
2
,
if the following inequalities hold: (1)
(4.11)
c
1
4
r
N
0
(
k
3
+
k
4
)
K
N
¯
>
(
d
N
C
N
+
c
1
r
N
0
K
N
̲
)
2
,
(2)
(4.12)
3
4
r
b
0
K
B
¯
>
(
d
N
B
+
N
*
r
N
0
η
1
¯
+
d
B
N
)
2
K
N
¯
r
N
0
+
(
d
B
C
B
+
c
2
r
B
0
/
K
B
̲
)
2
c
2
(
g
2
+
k
1
m
n
N
*
)
+
(
B
*
r
B
0
μ
¯
m
e
+
m
b
g
1
+
m
b
D
B
¯
+
m
n
N
*
C
N
*
d
N
B
)
2
m
e
(
m
e
h
+
m
b
g
1
B
*
+
m
n
k
2
N
*
)
,
(3)
(4.13)
3
m
e
4
(
m
e
h
+
m
b
g
1
B
*
+
m
n
k
2
N
*
)
>
K
N
¯
(
m
e
N
*
r
N
0
η
2
¯
+
m
n
k
3
+
m
n
k
2
+
m
n
D
N
¯
+
m
b
B
*
C
B
*
d
B
N
)
2
r
N
0
+
(
c
1
m
e
(
k
2
+
r
N
0
N
*
C
N
*
η
2
¯
)
+
N
*
m
n
(
k
3
+
D
N
¯
+
C
N
*
d
N
C
N
)
)
2
c
1
(
k
3
+
k
4
)
+
(
c
2
m
e
(
g
1
+
r
B
0
B
*
C
B
*
μ
¯
)
+
m
b
B
*
(
D
B
¯
+
C
B
*
d
B
C
B
)
)
2
c
2
(
g
2
+
k
1
m
n
N
*
)
,
then E
* is a globally stable equilibrium of system (2.23).


## 5. Numerical Simulation

To study the applicability of the model we need to test a range of parameter values for whichE
* exists and is stable.Example 5.1.
Consider the set of parameter values for our numerical calculations in Table1.
With these parameter values, the nontrivial equilibrium,E
*, of the model exists and is
(5.1)
N
*
≈
.
9919
,
B
*
≈
0.4934
,
C
N
*
≈
0.0007
,
C
B
*
≈
0.0003
,
C
E
*
≈
0.0095
.
Usingc
1
=
0.0751, c
2
=
1, and c
3
=
1, it can be verified that the conditions given by (4.4)–(4.6) for Theorem 4.1 are satisfied. Hence, E
* is locally asymptotically stable.
Moreover, usingc
1
=
1 and c
2
=
1, it can also be verified that the conditions given by (4.11)–(4.13) for Theorem 4.2 are satisfied. Hence, E
* is globally asymptotically stable.
Figure2(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 (
N
0
,
B
0
,
C
N
0
,
C
B
0
,
C
E
0
)
=
(
0.5,0.25,0
,
0,0
) was used for the simulation. Both time series were achieved using the MATLAB solver ode15s.Table 1
Population parameters
b
n
0
=
3
d
n
0
=
2
k
n
0
=
.
5
b
b
0
=
80
d
b
0
=
3
k
b
0
=
.
5
k
n
b
=
1
k
n
c
e
=
0.1667
d
n
b
=
0
d
n
c
n
=
0.01
k
b
c
e
=
0.0125
d
b
n
=
1
d
b
c
b
=
1
Mass parameters
m
n
=
1
m
b
=
1
m
e
=
1
Pollution parameters
k
1
=
3.01
k
2
=
1
k
3
=
1
k
4
=
52.01
g
1
=
5.5
g
2
=
100
U
=
1
h
=
100Numerical approximation to the time series solutions for the system of differential equations (2.23). The parameter values used were b
n
0
=
d
b
0
=
3, d
n
0
=
2, k
n
0
=
k
b
0
=
.
5, b
b
0
=
80, d
n
c
n
=
0.01, k
1
=
g
2
=
h
=
100, k
n
c
e
=
0.1667, k
b
c
e
=
0.0125, k
2
=
3.01, d
n
b
=
0, k
4
=
52.01, g
1
=
5.5, and k
n
b
=
d
b
n
=
d
b
c
b
=
m
n
=
m
b
=
m
e
=
k
3
=
U
=
1. The initial condition used was (
N
0
,
B
0
,
C
N
0
,
C
B
0
,
C
E
0
)
=
(
0.5,0.25,0
,
0,0
).
(a)
Numerical approximation ofN and B
(b)
Numerical approximation ofC
N
,
C
B and C
EExample 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 Table2.The nontrivial equilibrium, E
*, of the model exists and is
(5.2)
N
*
≈
.
7634
,
B
*
≈
0.4236
,
C
N
*
≈
0.0749
,
C
B
*
≈
0.02925
,
C
E
*
≈
0.9795
.
Usingc
1
=
0.6104
,
c
2
=
1, and c
3
=
1, the conditions given by (4.4)–(4.6) for Theorem 4.1 are satisfied, and therefore E
* is proven to be locally asymptotically stable.
Similarly, usingc
1
=
1 and c
2
=
1, the conditions given by (4.11)–(4.13) for Theorem 4.2 are satisfied and E
* is globally asymptotically stable.Table 2
Population parameters
b
n
0
=
3
d
n
0
=
2
k
n
0
=
.
5
b
b
0
=
9
d
b
0
=
3
k
b
0
=
.
5
k
n
b
=
1
k
n
c
e
=
0.1667
d
n
b
=
0
d
n
c
n
=
0.01
k
b
c
e
=
0.0125
d
b
n
=
1
d
b
c
b
=
1
Mass parameters
m
n
=
1
m
b
=
1
m
e
=
1
Pollution parameters
k
1
=
3.01
k
2
=
1
k
3
=
1
k
4
=
52.01
g
1
=
5.5
g
2
=
100
U
=
100
h
=
100Figure3(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 (
N
0
,
B
0
,
C
N
0
,
C
B
0
,
C
E
0
)
=
(
0.5,0.25,0
,
0,0
) 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 b
B
0
=
8, 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.Numerical approximation to the time series solutions for the system of differential equations (2.23). The parameter values used were b
n
0
=
d
b
0
=
3, d
n
0
=
2, k
n
0
=
k
b
0
=
.
5, b
b
0
=
9, d
n
c
n
=
0.01, k
1
=
g
2
=
h
=
U
=
100, k
n
c
e
=
0.1667, k
b
c
e
=
0.0125, k
2
=
3.01, d
n
b
=
0, k
4
=
52.01, g
1
=
5.5, and k
n
b
=
d
b
n
=
d
b
c
b
=
m
n
=
m
b
=
m
e
=
k
3
=
1. The initial condition used was (
N
0
,
B
0
,
C
N
0
,
C
B
0
,
C
E
0
)
=
(
0.5,0.25,0
,
0,0
).
(a)
Numerical approximation ofN and B
(b)
Numerical approximation ofC
N
,
C
B and C
E


## 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 Examples5.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.



---

*Source: 643985-2011-11-03.xml*