Table of Contents
ISRN Applied Mathematics
Volume 2011, Article ID 643985, 31 pages
http://dx.doi.org/10.5402/2011/643985
Research Article

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

1Department of Mathematics, Trent University, Peterborough, ON, K9J 7B8, Canada
2Department 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 [29]. 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 𝑑𝑁𝑑𝑡=(𝑏(𝑡)𝑑(𝑡))𝑁,(2.1) If we let 𝑟0=𝑏0𝑑0 where 𝑏0 is the intrinsic birth rate and 𝑑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𝑏(𝑡)=𝑏0𝑟0𝑁𝐾.(2.2) Similarly, the general death rate can be thought of as 𝑑(𝑡)=𝑑0.(2.3) 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 𝑏𝐵0 denote the intrinsic birth rate of the plant species, 𝑟𝐵0 the intrinsic growth rate of the plant species, 𝐾𝐵0 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𝑏𝐵𝐵,𝐶𝐸=𝑏𝐵0𝑟𝐵0𝐵𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸.(2.4)

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 𝑑𝐵0 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𝑑𝐵𝑁,𝐶𝐵=𝑑𝐵0+𝑑𝐵𝑁𝑁+𝑑𝐵𝐶𝐵𝐶𝐵.(2.5)

Using 𝑟𝐵0=𝑏𝐵0𝑑𝑏0, the equation governing the plant population is𝑑𝐵=𝑏𝑑𝑡𝐵𝐵,𝐶𝐸𝑑𝐵𝑁,𝐶𝐵𝐵,(2.6)𝑑𝐵=𝑟𝑑𝑡𝐵0𝑟𝐵0𝐵𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸𝑑𝐵𝑁𝑁𝑑𝐵𝐶𝐵𝐶𝐵𝐵.(2.7) 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 𝑏𝑁0,𝑟𝑁0, and 𝐾𝐵0 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 𝑏𝑁𝑁,𝐵,𝐶𝐸=𝑏𝑁0𝑟𝑁0𝑁𝐾𝑁0+𝐾𝑁𝐵𝐵𝐾𝑁𝐶𝐸𝐶𝐸.(2.8)

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 𝑑𝑁0 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𝑑𝑁𝐵,𝐶𝑁=𝑑𝑁0𝑑𝑁𝐵𝐵+𝑑𝑁𝐶𝑁𝐶𝑁.(2.9) Clearly, using 𝑟𝑁0=𝑏𝑁0𝑑𝑁0, we can describe the dynamics of the animal population as𝑑𝑁=𝑏𝑑𝑡𝑁𝑁,𝐵,𝐶𝐸𝑑𝑁𝐵,𝐶𝑁𝑁,(2.10)𝑑𝑁=𝑟𝑑𝑡𝑁0𝑟𝑁0𝑁𝐾𝑁0+𝐾𝑁𝐵𝐵𝐾𝑁𝐶𝐸𝐶𝐸+𝑑𝑁𝐵𝐵𝑑𝑁𝐶𝑁𝐶𝑁𝑁.(2.11) 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 [1113].

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.

643985.fig.001
Figure 1: Flowchart of the pollution movement through the ecosystem.

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 𝑔1𝑚𝑏𝐶𝐸𝐵. The loss of environmental pollution due to absorption into the animal population is given as 𝑘2𝑚𝑛𝐶𝐸𝑁. Pollution reenters the environment, as part of the pollution cycle as it is egested by the animal population. The egestion occurs at the rate 𝑘3𝑚𝑛𝐶𝑁𝑁. 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𝑑𝑚𝑒𝐶𝐸𝑑𝑡=𝑈(𝑡)𝑚𝑒𝐶𝐸𝑔1𝑚𝑏𝐶𝐸𝐵𝑘2𝑚𝑛𝐶𝐸𝑁+𝑘3𝑚𝑛𝐶𝑁𝑁+𝑚𝑛𝑑𝑁𝐵,𝐶𝑁𝐶𝑁𝑁+𝑚𝑏𝑑𝐵𝑁,𝐶𝐵𝐶𝐵𝐵;(2.12) then, dividing by 𝑚𝑒 and using the notation 𝑢(𝑡)=𝑈(𝑡)/𝑚𝑒, we have 𝑑𝐶𝐸𝑑𝑡=𝑢(𝑡)𝐶𝐸𝑔1𝑚𝑏𝑚𝑒𝐶𝐸𝑘𝐵2𝑚𝑛𝑚𝑒𝐶𝐸𝑁+𝑘3𝑚𝑛𝑚𝑒𝐶𝑁𝑚𝑁+𝑛𝑚𝑒𝑑𝑁0𝑑𝑁𝐵𝐵+𝑑𝑁𝐶𝑁𝐶𝑁𝐶𝑁𝑁+𝑚𝑏𝑚𝑒𝑏𝑁0𝑟𝑁0𝑁𝐾𝑁0+𝐾𝑁𝐵𝐵𝐾𝑁𝐶𝐸𝐶𝐸𝐶𝐵𝐵.(2.13) 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 𝑔1𝑚𝑏𝐶𝐸𝐵. Pollution inside the plant biomass is naturally metabolized, leaving the system at a rate 𝑔2𝑚𝑏𝐶𝐵𝐵. 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 𝑘1𝑚𝑏𝑚𝑛𝐶𝐵𝐵𝑁. So we get 𝑑𝐶𝐵𝑚𝑏𝐵𝑑𝑡=𝑔1𝑚𝑏𝐶𝐸𝐵𝑔2𝑚𝑏𝐶𝐵𝐵𝑑𝐵𝑁,𝐶𝐵𝑚𝑏𝐶𝐵𝐵𝑘1𝑚𝑏𝑚𝑛𝐶𝐵𝐵𝑁;(2.14) dividing by 𝑚𝑏, we have 𝑑𝐶𝐵𝐵𝑑𝑡=𝑔1𝐶𝐸𝐵𝑔2𝐶𝐵𝐵𝑘1𝑚𝑛𝐶𝐵𝐵𝑁𝑑𝐵𝑁,𝐶𝐵𝐶𝐵𝐵,(2.15) but clearly 𝑑𝐶𝐵𝐵𝑑𝑡=𝐶𝐵𝑑𝐵𝑑𝑡+𝐵𝑑𝐶𝐵.𝑑𝑡(2.16) Substituting in (2.6) and (2.15), using 𝑏𝐵𝐵,𝐶𝐸=𝑏𝐵0+𝑟𝐵0𝐵𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸,(2.17) and factoring, we have 𝑑𝐶𝐵𝑑𝑡=𝑔1𝐶𝐸𝑘1𝑚𝑛𝐶𝐵𝑔𝑁2+𝑏𝐵0𝑟𝐵0𝐵𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸𝐶𝐵.(2.18) 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 𝑘2𝑚𝑛𝐶𝐸𝑁 and by the ingestion of plant biomass containing pollution at a rate 𝑘1𝑚𝑛𝑚𝑏𝐶𝐵𝐵𝑁. Pollution is released from the animal population as a result of being egested, which occurs at a loss rate of 𝑘3𝑚𝑛𝐶𝑁𝑁, metabolized, which occurs at a loss rate of 𝑘4𝑚𝑛𝐶𝑁𝑁, and released from dead animals, which occurs at a rate of 𝑑𝑁(𝐵,𝐶𝑁)𝑚𝑛𝐶𝑁𝑁. Hence, we obtain 𝑑𝐶𝑁𝑚𝑛𝑁𝑑𝑡=𝑘1𝑚𝑛𝑚𝑏𝐶𝐵𝐵𝑁+𝑘2𝑚𝑛𝐶𝐸𝑁𝑘3𝑚𝑛𝐶𝑁𝑁𝑘4𝑚𝑛𝐶𝑁𝑁𝑑𝑁𝐵,𝐶𝑁𝑚𝑛𝐶𝑁𝑁.(2.19)

First, we divide by 𝑚𝑛 to get𝑑𝐶𝑁𝑁𝑑𝑡=𝑘1𝑚𝑏𝐶𝐵𝐵𝑁+𝑘2𝐶𝐸𝑁𝑘3𝐶𝑁𝑁𝑘4𝐶𝑁𝑁𝑑𝑁𝐵,𝐶𝑁𝐶𝑁𝑁.(2.20) Then, we use the fact that 𝑑𝐶𝑁𝑁𝑑𝑡=𝐶𝑁𝑑𝑁𝑑𝑡+𝑁𝑑𝐶𝑁𝑑𝑡.(2.21) Substituting in (2.10) and (2.20), using 𝑏𝑁(𝑁,𝐵,𝐶𝐸)=𝑏𝑁0+𝑟𝑁0𝑁/(𝐾𝑁0+𝐾𝑁𝐵𝐵𝐾𝑁𝐶𝐸𝐶𝐸), and factoring we determine that 𝑑𝐶𝑁𝑑𝑡=𝑘1𝑚𝑏𝐶𝐵𝐵+𝑘2𝐶𝐸𝑘3+𝑘4+𝑏𝑁0𝑟𝑁0𝑁𝐾𝑁0+𝐾𝑁𝐵𝐵𝐾𝑁𝐶𝐸𝐶𝐸𝐶𝑁.(2.22) 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:𝑑𝑁=𝑟𝑑𝑡𝑁0𝑟𝑁0𝑁𝐾𝑁0+𝐾𝑁𝐵𝐵𝐾𝑁𝐶𝐸𝐶𝐸+𝑑𝑁𝐵𝐵𝑑𝑁𝐶𝑁𝐶𝑁𝑁,𝑑𝐵=𝑟𝑑𝑡𝐵0𝑟𝐵0𝐵𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸𝑑𝐵𝑁𝑁𝑑𝐵𝐶𝐵𝐶𝐵𝐵,𝑑𝐶𝑁𝑑𝑡=𝑘1𝑚𝑏𝐶𝐵𝐵+𝑘2𝐶𝐸𝑘3+𝑘4+𝑏𝑁0𝑟𝑁0𝑁𝐾𝑁0+𝐾𝑁𝐵𝐵𝐾𝑁𝐶𝐸𝐶𝐸𝐶𝑁,𝑑𝐶𝐵𝑑𝑡=𝑔1𝐶𝐸𝑘1𝑚𝑛𝐶𝐵𝑔𝑁2+𝑏𝐵0𝑟𝐵0𝐵𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸𝐶𝐵,𝑑𝐶𝐸𝑚𝑑𝑡=𝑢(𝑡)+𝑏𝑚𝑒𝑑𝐵0+𝑑𝐵𝑁𝑁+𝑑𝐵𝐶𝐵𝐶𝐵𝐶𝐵𝐵+𝑚𝑛𝑚𝑒𝑘3+𝑑𝑁0𝑑𝑁𝐵𝐵+𝑑𝑁𝐶𝑁𝐶𝑁𝐶𝑁𝑁𝑚𝑏𝑚𝑒𝑔1𝐶𝐸𝑚𝐵𝑛𝑚𝑒𝑘2𝐶𝐸𝑁𝐶𝐸.(2.23)

3. Equilibrium Analysis

Our model has four nonnegative equilibria: 𝐸0=0,0,𝐶𝑁0,𝐶𝐵0,𝐶𝐸0,𝐸1=0,𝐵1,𝐶𝑁1,𝐶𝐵1,𝐶𝐸1,𝐸2=𝑁2,0,𝐶𝑁2,𝐶𝐵2,𝐶𝐸2,𝐸=𝑁,𝐵,𝐶𝑁,𝐶𝐵,𝐶𝐸.(3.1) 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, 0=𝑟𝑁0𝑟𝑁0𝑁𝐾𝑁0+𝐾𝑁𝐵𝐵𝐾𝑁𝐶𝐸𝐶𝐸+𝑑𝑁𝐵𝐵𝑑𝑁𝐶𝑁𝐶𝑁,(3.2)0=𝑟𝐵0𝑟𝐵0𝐵𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸𝑑𝐵𝑁𝑁𝑑𝐵𝐶𝐵𝐶𝐵,(3.3)0=𝑘1𝑚𝑏𝐶𝐵𝐵+𝑘2𝐶𝐸𝑘3+𝑘4+𝑏𝑁0𝑟𝑁0𝑁𝐾𝑁0+𝐾𝑁𝐵𝐵𝐾𝑁𝐶𝐸𝐶𝐸𝐶𝑁,(3.4)0=𝑔1𝐶𝐸𝑘1𝑚𝑛𝐶𝐵𝑁𝑔2+𝑏𝐵0𝑟𝐵0𝐵𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸𝐶𝐵𝑚𝑏𝑚𝑒𝑔1𝐶𝐸𝐵𝑚,(3.5)0=𝑢(𝑡)+𝑏𝑚𝑒𝑑𝐵0+𝑑𝐵𝑁𝑁+𝑑𝐵𝐶𝐵𝐶𝐵𝐶𝐵𝐵𝑚𝑏𝑚𝑒𝑔1𝐶𝐸𝐵+𝑚𝑛𝑚𝑒𝑘3+𝑑𝑁0𝑑𝑁𝐵𝐵+𝑑𝑁𝐶𝑁𝐶𝑁𝐶𝑁𝑁𝑚𝑛𝑚𝑒𝑘2𝐶𝐸𝑁𝐶𝐸.(3.6)

Recall that 𝑟𝐵0=𝑏𝐵0𝑑𝐵0, and so, by (3.3), 𝑏𝐵0𝑟𝐵0𝐵𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸=𝑑𝑏0+𝑑𝐵𝑁𝑁+𝑑𝐵𝐶𝐵𝐶𝐵.(3.7) Substituting this into (3.5), we have0=𝑑𝐵𝐶𝐵𝐶𝐵2+𝑘1𝑚𝑛𝑁+𝑔2+𝑑𝑏0+𝑑𝐵𝑁𝑁𝐶𝐵𝑔1𝐶𝐸.(3.8) Hence, 𝐶𝐵=𝛽1𝑁+𝛽1(𝑁)2+4𝑔1𝑑𝐵𝐶𝐵𝐶𝐸2𝑑𝐵𝐶𝐵,(3.9) where, 𝛽1𝑁=𝑘1𝑚𝑛𝑁+𝑑𝐵𝑁𝑁+𝑔2+𝑑𝐵0.(3.10) Let𝑓1𝑁,𝐶𝐸=𝛽1𝑁+𝛽1(𝑁)2+4𝑔1𝑑𝐵𝐶𝐵𝐶𝐸2𝑑𝐵𝐶𝐵.(3.11) Consider (3.3). Since 𝐶𝐵=𝑓1(𝑁,𝐶𝐸), 𝐵=𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸𝑟𝐵0𝑑𝐵𝑁𝑁𝑑𝐵𝐶𝐵𝑓1𝑁,𝐶𝐸𝑟𝐵0.(3.12) Let 𝑓2𝑁,𝐶𝐸=𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸𝑟𝐵0𝑑𝐵𝑁𝑁𝑑𝐵𝐶𝐵𝑓1𝑁,𝐶𝐸𝑟𝐵0.(3.13) Next, recall that 𝑟𝑁0=𝑏𝑁0𝑑𝑁0, and so, by (3.2),𝑏𝑁0𝑟𝑁0𝑁𝐾𝑁0+𝐾𝑁𝐵𝐵𝐾𝑁𝐶𝐸𝐶𝐸=𝑑𝑁0𝑑𝑁𝐵𝐵+𝑑𝑁𝐶𝑁𝐶𝑁.(3.14) Substituting this into (3.4), we have0=𝑑𝑁𝐶𝑁𝐶𝑁2+𝑘3+𝑘4+𝑑𝑁0𝑑𝑁𝐵𝐵𝐶𝑁𝑘1𝑚𝑏𝐶𝐵𝐵+𝑘2𝐶𝐸.(3.15) Using 𝐶𝐵=𝑓1(𝑁,𝐶𝐸) and 𝐵=𝑓1(𝑁,𝐶𝐸), we find 𝐶𝑁=𝛽2𝑁,𝐶𝐸+𝛽2𝑁,𝐶𝐸2+𝛽3𝑁,𝐶𝐸2𝑑𝑁𝐶𝑁,(3.16) where 𝛽2𝑁,𝐶𝐸=𝑘3+𝑘4+𝑑𝑁0𝑑𝑁𝐵𝑓2𝑁,𝐶𝐸,𝛽3𝑁,𝐶𝐸=4𝑑𝑁𝐶𝑁𝑘1𝑚𝑏𝑓1𝑁,𝐶𝐸𝑓2𝑁,𝐶𝐸+𝑘2𝐶𝐸.(3.17) Let 𝑓3𝑁,𝐶𝐸=𝛽2𝑁,𝐶𝐸+𝛽2𝑁,𝐶𝐸2+𝛽3𝑁,𝐶𝐸2𝑑𝑁𝐶𝑁.(3.18) Next, if we consider (3.2), we have 0=𝑟𝑁0𝑟𝑁0𝑁𝐾𝑁0+𝐾𝑁𝐵𝐵𝐾𝑁𝐶𝐸𝐶𝐸+𝑑𝑁𝐵𝐵𝑑𝑁𝐶𝑁𝐶𝑁.(3.19) Using 𝐵=𝑓2(𝑁,𝐶𝐸) and 𝐶𝑁=𝑓3(𝑁,𝐶𝐸), we find 𝑁=𝛽4𝑁,𝐶𝐸𝑟𝑁0+𝑑𝑁𝐵𝑓2𝑁,𝐶𝐸𝑑𝑁𝐶𝑁𝑓3𝑁,𝐶𝐸𝑟𝑁0,(3.20) where 𝛽4𝑁,𝐶𝐸=𝐾𝑁0+𝐾𝑁𝐵𝑓2𝑁,𝐶𝐸𝐾𝑁𝐶𝐸𝐶𝐸.(3.21) Finally, we consider (3.6) and 𝑚0=𝑢(𝑡)+𝑏𝑚𝑒𝑑𝐵0+𝑑𝐵𝑁𝑁+𝑑𝐵𝐶𝐵𝐶𝐵𝐶𝐵𝐵𝑚𝑏𝑚𝑒𝑔1𝐶𝐸𝐵+𝑚𝑛𝑚𝑒𝑘3+𝑑𝑁0𝑑𝑁𝐵𝐵+𝑑𝑁𝐶𝑁𝐶𝑁𝐶𝑁𝑁𝑚𝑛𝑚𝑒𝑘2𝐶𝐸𝑁𝐶𝐸.(3.22) Substituting in 𝐶𝐵=𝑓1(𝑁,𝐶𝐸), 𝐵=𝑓2(𝑁,𝐶𝐸), and 𝐶𝑁=𝑓3(𝑁,𝐶𝐸), we find𝐶𝐸=𝑚𝑒𝑢+𝑚𝑏𝑑𝐵0+𝑑𝐵𝑁𝑁+𝑑𝐵𝐶𝐵𝑓1𝑓1𝑓2𝑚𝑏𝑔1𝑓2+𝑚𝑛𝑘2𝑁+𝑚𝑒+𝑚𝑛𝑘3+𝑑𝑁0𝑑𝑁𝐵𝑓2+𝑑𝑁𝐶𝑁𝑓3𝑓3𝑁𝑚𝑏𝑔1𝑓2+𝑚𝑛𝑘2𝑁+𝑚𝑒.(3.23)

Thus, the nontrivial solution to the system is given by 𝑁=𝛽4𝑁,𝐶𝐸𝑟𝑁0𝑑𝑁𝐶𝑁𝑓3𝑁,𝐶𝐸+𝑑𝑁𝐵𝑓2𝑁,𝐶𝐸𝑟𝑁0,𝐵=𝑓2𝑁,𝐶𝐸,𝐶𝑁=𝑓3𝑁,𝐶𝐸,𝐶𝐵=𝑓1𝑁,𝐶𝐸,𝐶𝐸=𝑚𝑒𝑢+𝑚𝑏𝛽6𝑁,𝐶𝐸𝑓1𝑁,𝐶𝐸𝑓2𝑁,𝐶𝐸𝛽5𝑁,𝐶𝐸+𝑚𝑛𝛽7𝑁,𝐶𝐸𝑓3𝑁,𝐶𝐸𝑁𝛽5𝑁,𝐶𝐸,(3.24) where 𝑓1𝑁,𝐶𝐸=𝛽1𝑁+𝛽1(𝑁)2+4𝑔1𝑑𝐵𝐶𝐵𝐶𝐸2𝑑𝐵𝐶𝐵,𝑓2𝑁,𝐶𝐸=𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸𝑟𝐵0𝑑𝐵𝑁𝑁+𝑑𝐵𝐶𝐵𝑓1𝑁,𝐶𝐸𝑟𝐵0,𝑓3𝑁,𝐶𝐸=𝛽2𝑁,𝐶𝐸+𝛽2𝑁,𝐶𝐸2+𝛽3𝑁,𝐶𝐸2𝑑𝑁𝐶𝑁,(3.25) where 𝛽1𝑁=𝑘1𝑚𝑛𝑁+𝑑𝐵𝑁𝑁+𝑔2+𝑑𝐵0,𝛽2𝑁,𝐶𝐸=𝑘3+𝑘4+𝑑𝑁0𝑑𝑁𝐵𝑓2𝑁,𝐶𝐸,𝛽3𝑁,𝐶𝐸=4𝑑𝑁𝐶𝑁𝑘1𝑚𝑏𝑓1𝑁,𝐶𝐸𝑓2𝑁,𝐶𝐸+𝑘2𝐶𝐸,𝛽4𝑁,𝐶𝐸=𝐾𝑁0+𝐾𝑁𝐵𝑓2𝑁,𝐶𝐸𝐾𝑁𝐶𝐸𝐶𝐸,𝛽5𝑁,𝐶𝐸=𝑚𝑏𝑔1𝑓2𝑁,𝐶𝐸+𝑚𝑛𝑘2𝑁+𝑚𝑒𝛽,6𝑁,𝐶𝐸=𝑑𝐵0+𝑑𝐵𝑁𝑁+𝑑𝐵𝐶𝐵𝑓1𝑁,𝐶𝐸,𝛽7𝑁,𝐶𝐸=𝑘3+𝑑𝑁0𝑑𝑁𝐵𝑓2𝑁,𝐶𝐸+𝑑𝑁𝐶𝑁𝑓3𝑁,𝐶𝐸.(3.26) Provided that 𝑁 and 𝐶𝐸 exist, we can find the corresponding 𝐶𝐵,𝐵, and 𝐶𝑁.

4. Stability Analysis

Theorem 4.1 (conditions for local stability). Given 𝑐3<14𝑟𝑁0𝑔2+𝑘1𝑚𝑛𝑁+𝑏𝐵𝐵,𝐶𝐸𝐾𝑁𝐵,𝐶𝐸𝑘1𝑚𝑛𝐶𝐵2,𝑐1=𝑐3𝑟𝐵0𝐶𝐵𝑑𝐵𝐶𝐵𝐾𝐵𝐶𝐸,𝑐2<𝑐34𝑘3+𝑘4+𝑏𝑁𝑁,𝐵,𝐶𝐸Ω𝑁,𝐵,𝐶𝑁,𝐶𝐵,𝐶𝐸,(4.1) where Ω𝑁,𝐵,𝐶𝑁,𝐶𝐵,𝐶𝐸𝐶=min𝐵𝑑𝐵𝐶𝐵𝑟𝐵0𝐾𝐵𝐶𝐸𝛿𝑁,𝐵,𝐶𝑁,𝐶𝐵,𝐶𝐸2,𝑘1𝑚𝑛𝑁+𝑔2+𝑏𝐵𝐵,𝐶𝐸𝑘1𝑚𝑏𝐵2,(4.2) where 𝛿𝑁,𝐵,𝐶𝑁,𝐶𝐵,𝐶𝐸=𝑘1𝑚𝑏𝐶𝐵𝐾𝑁𝐵𝑟𝑁0𝑁𝐶𝑁𝐾𝑁𝐵,𝐶𝐸2.(4.3)
if the following inequalities hold:
(1)𝑑𝑁𝐶𝑁𝑐2𝑟𝑁0𝐶𝑁𝐾𝑁𝐵,𝐶𝐸2<𝑐24𝑟𝑁0𝑘3+𝑘4+𝑏𝑁𝑁,𝐵,𝐶𝐸𝐾𝑁𝐵,𝐶𝐸,(4.4)
(2)3𝑐14𝑟𝐵0𝐾𝐵𝐶𝐸>𝑑𝑁𝐵+𝐾𝑁𝐵𝑟𝑁0𝑁/𝐾𝑁𝐵,𝐶𝐸2𝑐1𝑑𝐵𝑁2𝑟𝑁0/𝐾𝑁𝐵,𝐶𝐸+𝑐1𝑚𝑒𝐾𝐵𝐶𝐸𝑟𝐵0𝐵/𝐾𝐵𝐶𝐸2𝑚𝑏𝑑𝐵𝑁,𝐶𝐵𝐶𝐵+𝑚𝑛𝑑𝑁𝐵𝐶𝑁𝑁+𝑚𝑏𝑔1𝐶𝐸2𝑚𝑒𝑚𝑏𝑔1𝐵+𝑚𝑛𝑘2𝑁+𝑚𝑒,(4.5)
(3)
3𝑚𝑒4𝑚𝑏𝑔1𝐵+𝑚𝑛𝑘2𝑁+𝑚𝑒>𝑚𝑒𝐾𝑁𝐶𝐸𝑟𝑁0𝑁/𝐾𝑁𝐵,𝐶𝐸2𝑚𝑏𝑑𝐵𝑁𝐶𝐵𝐵𝑚𝑛𝑘3+𝑑𝑁𝐵,𝐶𝑁𝐶𝑁+𝑚𝑛𝑘2𝐶𝐸2𝑟𝑁0/𝐾𝑁𝐵,𝐶𝐸+𝑐2𝑚𝑒𝑘2+𝐾𝑁𝐶𝐸𝑟𝑁0𝑁𝐶𝑁/𝐾𝑁𝐵,𝐶𝐸2+𝑚𝑛𝑁𝑘3+𝑑𝑁𝐵,𝐶𝑁+𝑑𝑁𝐶𝑁𝐶𝑁2𝑐2𝑘3+𝑘4+𝑏𝑁𝑁,𝐵,𝐶𝐸+𝑐3𝑚𝑒𝑔1+𝐾𝐵𝐶𝐸𝑟𝐵0𝐵𝐶𝐵/𝐾𝐵𝐶𝐸2+𝑚𝑏𝐵𝑑𝐵𝑁,𝐶𝐵+𝑑𝐵𝐶𝐵𝐶𝐵2𝑐3𝑘1𝑚𝑛𝑁+𝑔2+𝑏𝐵𝐵,𝐶𝐸,(4.6) then 𝐸 is a locally stable equilibrium of system (2.23).

Theorem 4.2 (conditions for global stability). Using the notation 𝐾𝑁𝐵,𝐶𝐸=𝐾𝑁0+𝐾𝑁𝐵𝐵𝐾𝑁𝐶𝐸𝐶𝐸,𝐾𝐵𝐶𝐸=𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸,𝐷𝑁𝐵,𝐶𝑁=𝑑𝑁0𝑑𝑁𝐵𝐵+𝑑𝑁𝐶𝑁𝐶𝑁,𝐷𝐵𝑁,𝐶𝐵=𝑑𝐵0+𝑑𝐵𝑁𝑁+𝑑𝐵𝐶𝐵𝐶𝐵,(4.7) let 𝐾𝑁(𝐵,𝐶𝐸),𝐾𝐵(𝐶𝐸),𝐷𝑁(𝐵,𝐶𝑁), and 𝐷𝐵(𝑁,𝐶𝐵) be such that 𝐾𝑁=𝐾𝑁0𝐾𝑁𝐶𝐸𝐾𝑁𝐵,𝐶𝐸𝐾𝑁0+𝐾𝑁𝐵𝐾𝐵0=𝐾𝑁,𝐾𝐵=𝐾𝐵0𝐾𝐵𝐶𝐸𝐾𝐵𝐶𝐸𝐾𝐵0=𝐾𝐵,||𝐷𝑁𝐵,𝐶𝑁||𝑑𝑁0+𝑑𝑁𝐶𝑁=𝐷𝑁,||𝐷𝐵𝑁,𝐶𝐵||𝑑𝐵0+𝑑𝐵𝑁𝑁𝑀+𝑑𝐵𝐶𝐵=𝐷𝐵.(4.8) Then, given 0<𝑐2<𝑟𝑁0𝑔2+𝑘1𝑚𝑛𝑁4𝑘1𝑚𝑛2𝐾𝑁,0<𝑐1<𝑘3+𝑘44𝑟max𝐵0𝐾𝐵𝑘1𝑚𝑏+𝑟𝑁0𝑁𝐶𝑁𝜂12𝑐2𝑔2+𝑘1𝑚𝑛𝑁𝑘1𝑚𝑏𝐵2,(4.9) where 𝜂1=𝐾𝑁𝐵𝐾𝑁0+𝐵𝐾𝑁𝐵𝐾𝑁𝐶𝐸2,𝜂2=𝐾𝑁𝐶𝐸𝐾𝑁0+𝐵𝐾𝑁𝐵𝐾𝑁𝐶𝐸𝐶𝐸2,𝐾𝜇=𝐵𝐶𝐸𝐾𝐵0𝐾𝐵𝐶𝐸𝐶𝐸2,(4.10) if the following inequalities hold: (1)𝑐14𝑟𝑁0(𝑘3+𝑘4)𝐾𝑁>𝑑𝑁𝐶𝑁+𝑐1𝑟𝑁0𝐾𝑁2,(4.11)(2)34𝑟𝑏0𝐾𝐵>𝑑𝑁𝐵+𝑁𝑟𝑁0𝜂1+𝑑𝐵𝑁2𝐾𝑁𝑟𝑁0+𝑑𝐵𝐶𝐵+𝑐2𝑟𝐵0/𝐾𝐵2𝑐2𝑔2+𝑘1𝑚𝑛𝑁+𝐵𝑟𝐵0𝜇𝑚𝑒+𝑚𝑏𝑔1+𝑚𝑏𝐷𝐵+𝑚𝑛𝑁𝐶𝑁𝑑𝑁𝐵2𝑚𝑒𝑚𝑒+𝑚𝑏𝑔1𝐵+𝑚𝑛𝑘2𝑁,(4.12)(3)3𝑚𝑒4𝑚𝑒+𝑚𝑏𝑔1𝐵+𝑚𝑛𝑘2𝑁>𝐾𝑁𝑚𝑒𝑁𝑟𝑁0𝜂2+𝑚𝑛𝑘3+𝑚𝑛𝑘2+𝑚𝑛𝐷𝑁+𝑚𝑏𝐵𝐶𝐵𝑑𝐵𝑁2𝑟𝑁0+𝑐1𝑚𝑒𝑘2+𝑟𝑁0𝑁𝐶𝑁𝜂2+𝑁𝑚𝑛𝑘3+𝐷𝑁+𝐶𝑁𝑑𝑁𝐶𝑁2𝑐1𝑘3+𝑘4+𝑐2𝑚𝑒(𝑔1+𝑟𝐵0𝐵𝐶𝐵𝜇)+𝑚𝑏𝐵(𝐷𝐵+𝐶𝐵𝑑𝐵𝐶𝐵)2𝑐2𝑔2+𝑘1𝑚𝑛𝑁,(4.13) then 𝐸 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 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 𝑁.9919,𝐵0.4934,𝐶𝑁0.0007,𝐶𝐵0.0003,𝐶𝐸0.0095.(5.1)
Using 𝑐1=0.0751, 𝑐2=1, and 𝑐3=1, 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 𝑐1=1 and 𝑐2=1, 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 (𝑁0,𝐵0,𝐶𝑁0,𝐶𝐵0,𝐶𝐸0)=(0.5,0.25,0,0,0) was used for the simulation. Both time series were achieved using the MATLAB solver ode15s.

tab1
Table 1
fig2
Figure 2: Numerical approximation to the time series solutions for the system of differential equations (2.23). The parameter values used were 𝑏𝑛0=𝑑𝑏0=3, 𝑑𝑛0=2, 𝑘𝑛0=𝑘𝑏0=.5, 𝑏𝑏0=80, 𝑑𝑛𝑐𝑛=0.01, 𝑘1=𝑔2==100, 𝑘𝑛𝑐𝑒=0.1667, 𝑘𝑏𝑐𝑒=0.0125, 𝑘2=3.01, 𝑑𝑛𝑏=0, 𝑘4=52.01, 𝑔1=5.5, and 𝑘𝑛𝑏=𝑑𝑏𝑛=𝑑𝑏𝑐𝑏=𝑚𝑛=𝑚𝑏=𝑚𝑒=𝑘3=𝑈=1. The initial condition used was (𝑁0,𝐵0,𝐶𝑁0,𝐶𝐵0,𝐶𝐸0)=(0.5,0.25,0,0,0).

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 𝑁.7634,𝐵0.4236,𝐶𝑁0.0749,𝐶𝐵0.02925,𝐶𝐸0.9795.(5.2)
Using 𝑐1=0.6104,𝑐2=1, and 𝑐3=1, 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 𝑐1=1 and 𝑐2=1, the conditions given by (4.11)–(4.13) for Theorem 4.2 are satisfied and 𝐸 is globally asymptotically stable.

tab2
Table 2

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 (𝑁0,𝐵0,𝐶𝑁0,𝐶𝐵0,𝐶𝐸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 𝑏𝐵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.

fig3
Figure 3: Numerical approximation to the time series solutions for the system of differential equations (2.23). The parameter values used were 𝑏𝑛0=𝑑𝑏0=3, 𝑑𝑛0=2, 𝑘𝑛0=𝑘𝑏0=.5, 𝑏𝑏0=9, 𝑑𝑛𝑐𝑛=0.01, 𝑘1=𝑔2==𝑈=100, 𝑘𝑛𝑐𝑒=0.1667, 𝑘𝑏𝑐𝑒=0.0125, 𝑘2=3.01, 𝑑𝑛𝑏=0, 𝑘4=52.01, 𝑔1=5.5, and 𝑘𝑛𝑏=𝑑𝑏𝑛=𝑑𝑏𝑐𝑏=𝑚𝑛=𝑚𝑏=𝑚𝑒=𝑘3=1. The initial condition used was (𝑁0,𝐵0,𝐶𝑁0,𝐶𝐵0,𝐶𝐸0)=(0.5,0.25,0,0,0).

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

A. Proof of Theorem 4.1

We will work with the linearized form of the model, which can be written as𝑑𝑁𝑑𝑡=𝑗11𝑁+𝑗12𝐵𝑗13𝐶𝑁𝑗15𝐶𝐸,𝑑𝐵𝑑𝑡=𝑗21𝑁𝑗22𝐵𝑗24𝐵𝐶𝐵𝑗25𝐶𝐸,𝑑𝐶𝑁𝑑𝑡=𝑗31𝑁+𝑗32𝐵𝑗33𝐶𝑁+𝑗34𝐶𝐵+𝑗35𝐶𝐸,𝑑𝐶𝐵𝑑𝑡=𝑗41𝑁+𝑗42𝐵𝑗44𝐶𝐵+𝑗45𝐶𝐸,𝑑𝐶𝐸𝑑𝑡=𝑗51𝑁+𝑗52𝐵+𝑗53𝐶𝑁+𝑗54𝐶𝐵𝑗55𝐶𝐸,(A.1) where 𝑗𝑖𝑗 represents the absolute value from the 𝑖th row and 𝑗th column of the Jacobian matrix evaluated at 𝐸=(𝑁,𝐵,𝐶𝑁,𝐶𝐵,𝐶𝐸).

Consider the positive definite function𝑉𝐶𝑁,𝐵,𝑁,𝐶𝐵,𝐶𝐸=121𝑁𝑁𝑁2+𝑐1𝐵𝐵𝐵2+𝑐2𝐶𝑁𝐶𝑁2+𝑐3𝐶𝐵𝐶𝐵2+𝐶𝐸𝐶𝐸2.(A.2)

Differentiating 𝑉 with respect to 𝑡, we get𝑑𝑉=1𝑑𝑡𝑁𝑁𝑁𝑑𝑁+𝑐𝑑𝑡1𝐵𝐵𝐵𝑑𝐵𝑑𝑡+𝑐2𝐶𝑁𝐶𝑁𝑑𝐶𝑁𝑑𝑡+𝑐3𝐶𝐵𝐶𝐵𝑑𝐶𝐵+𝐶𝑑𝑡𝐸𝐶𝐸𝑑𝐶𝐸.𝑑𝑡(A.3)

Using the notation 2𝑎11=𝑗11=𝑟𝑁0𝐾𝑁𝐵,𝐶𝐸,2𝑎22=𝑐1𝑗22=𝑐1𝑟𝐵0𝐾𝐵𝐶𝐸,2𝑎33=𝑐2𝑗33=𝑐2𝑘3+𝑘4+𝑏𝑁𝑁,𝐵,𝐶𝐸,2𝑎44=𝑐3𝑗44=𝑐3𝑘1𝑚𝑛𝑁+𝑔2+𝑏𝐵𝐵,𝐶𝐸,2𝑎55=𝑗55=𝑚𝑏𝑚𝑒𝑔1𝐵+𝑚𝑛𝑚𝑒𝑘2𝑁𝑎+,12=𝑗12𝑁𝑐1𝑗21𝐵=𝑑𝑁𝐵+𝐾𝑁𝐵𝑟𝑁0𝑁𝐾𝑁𝐵,𝐶𝐸2𝑐1𝑑𝐵𝑁,𝑎13=𝑗13𝑁𝑐2𝑗31=𝑑𝑁𝐶𝑁𝑐2𝑟𝑁0𝐶𝑁𝐾𝑁𝐵,𝐶𝐸,𝑎14=𝑐3𝑗14=𝑐3𝑘1𝐶𝐵,𝑎15=𝑗15𝑁𝑗51=𝐾𝑁𝐶𝐸𝑟𝑁0𝑁𝐾𝑁𝐵,𝐶𝐸2𝑚𝑏𝑚𝑒𝑑𝐵𝑁𝐶𝐵𝐵𝑚𝑛𝑚𝑒𝑘3+𝑑𝑁𝐵,𝐶𝑁𝐶𝑁+𝑚𝑛𝑚𝑒𝑘2𝐶𝐸,𝑎23=𝑐2𝑗32=𝑐2𝑘1𝑚𝑏𝐶𝐵𝐾𝑁𝐵𝑟𝑁0𝑁𝐶𝑁𝐾𝑁𝐵,𝐶𝐸2,𝑎24=𝑐1𝑗24𝐵𝑐3𝑗42=𝑐1𝑑𝐵𝐶𝐵𝑐3𝑟𝐵0𝐶𝐵𝐾𝐵𝐶𝐸,𝑎25=𝑐1𝑗25𝐵𝑗52=𝑐1𝐾𝐵𝐶𝐸𝑟𝐵0𝐵𝐾𝐵𝐶𝐸2𝑚𝑏𝑚𝑒𝑑𝐵𝑁,𝐶𝐵𝐶𝐵+𝑚𝑛𝑚𝑒𝑑𝑁𝐵𝐶𝑁𝑁+𝑚𝑏𝑚𝑒𝑔1𝐶𝐸,𝑎34=𝑐2𝑗34=𝑐2𝑘1𝑚𝑏𝐵,𝑎35=𝑐2𝑗35+𝑗53=𝑐2𝑘2+𝐾𝑁𝐶𝐸𝑟𝑁0𝑁𝐶𝑁𝐾𝑁𝐵,𝐶𝐸2,+𝑚𝑛𝑚𝑒𝑁𝑘3+𝑑𝑁𝐵,𝐶𝑁+𝑑𝑁𝐶𝑁𝐶𝑁,𝑎45=𝑐3𝑗45+𝑗54=𝑐3𝑔1+𝐾𝐵𝐶𝐸𝑟𝐵0𝐵𝐶𝐵𝐾𝐵𝐶𝐸2,+𝑚𝑏𝑚𝑒𝐵𝑑𝐵𝑁,𝐶𝐵+𝑑𝐵𝐶𝐵𝐶𝐵.(A.4) and some algebraic manipulation, we can rewrite the system as 𝑑𝑉𝑎𝑑𝑡=112𝑁𝑁2+𝑎12𝑁𝑁𝐵𝐵𝑎222𝐵𝐵2𝑎112𝑁𝑁2𝑎13𝑁𝑁𝐶𝑁𝐶𝑁𝑎332𝐶𝑁𝐶𝑁2𝑎112𝑁𝑁2𝑎14𝑁𝑁𝐶𝐵𝐶𝐵𝑎442𝐶𝐵𝐶𝐵2𝑎112𝑁𝑁2𝑎15𝑁𝑁𝐶𝐸𝐶𝐸𝑎552𝐶𝐸𝐶𝐸2𝑎222𝐵𝐵2+𝑎23𝐵𝐵𝐶𝑁𝐶𝑁𝑎332𝐶𝑁𝐶𝑁2𝑎222𝐵𝐵2𝑎24𝐵𝐵𝐶𝐵𝐶𝐵𝑎442𝐶𝐵𝐶𝐵2𝑎222𝐵𝐵2𝑎25𝐵𝐵𝐶𝐸𝐶𝐸𝑎552𝐶𝐸𝐶𝐸2𝑎332𝐶𝑁𝐶𝑁2+𝑎34𝐶𝑁𝐶𝑁𝐶𝐵𝐶𝐵𝑎442𝐶𝐵𝐶𝐵2𝑎332𝐶𝑁𝐶𝑁2+𝑎35𝐶𝑁𝐶𝑁𝐶𝐸𝐶𝐸𝑎552𝐶𝐸𝐶𝐸2𝑎442𝐶𝐵𝐶𝐵2+𝑎45𝐶𝐵𝐶𝐵𝐶𝐸𝐶𝐸𝑎552𝐶𝐸𝐶𝐸2.(A.5)

Then, factoring the equations, we obtain 𝑑𝑉1𝑑𝑡=2𝑎11𝑁𝑁𝑎12𝑎11𝐵𝐵212𝑎22𝑎212𝑎11𝐵𝐵212𝑎11(𝑁𝑁𝑎)+13𝑎11𝐶𝑁𝐶𝑁212𝑎33𝑎213𝑎11𝐶𝑁𝐶𝑁212𝑎11𝑁𝑁+𝑎14𝑎11𝐶𝐵𝐶𝐵212𝑎44𝑎214𝑎11𝐶𝐵𝐶𝐵212𝑎11𝑁𝑁+𝑎15𝑎11𝐶𝐸𝐶𝐸212𝑎55𝑎215𝑎11𝐶𝐸𝐶𝐸212𝑎22𝐵𝐵𝑎23𝑎22𝐶𝑁𝐶𝑁212𝑎33𝑎223𝑎22𝐶𝑁𝐶𝑁212𝑎44𝐶𝐵𝐶𝐵+𝑎24𝑎44𝐵𝐵212𝑎22𝑎224𝑎44𝐵𝐵212𝑎55𝐶𝐸𝐶𝐸+𝑎25𝑎55𝐵𝐵212𝑎22𝑎225𝑎55𝐵𝐵212𝑎33𝐶𝑁𝐶𝑁𝑎34𝑎33𝐶𝐵𝐶𝐵212𝑎44𝑎234𝑎33𝐶𝐵𝐶𝐵212𝑎33𝐶𝑁𝐶𝑁𝑎35𝑎33𝐶𝐸𝐶𝐸212𝑎55𝑎235𝑎33𝐶𝐸𝐶𝐸212𝑎44𝐶𝐵𝐶𝐵𝑎45𝑎44𝐶𝐸𝐶𝐸212𝑎55𝑎245𝑎44𝐶𝐸𝐶𝐸2.(A.6)

Sufficient conditions for 𝑑𝑉/𝑑𝑡 to be negative definite are𝑎214<𝑎11𝑎44,𝑎(A.7)223<𝑎22𝑎33𝑎,(A.8)234<𝑎33𝑎44𝑎,(A.9)213<𝑎11𝑎33𝑎,(A.10)212𝑎11+𝑎224𝑎44+𝑎225𝑎55<3𝑎22𝑎,(A.11)215𝑎11+𝑎235𝑎33+𝑎245𝑎44<3𝑎55.(A.12)

Examining (A.7), we find𝑎214<𝑎11𝑎44,𝑐3<14𝑟𝑁0𝑔2+𝑘1𝑁+𝑏𝐵𝐵,𝐶𝐸𝐾𝑁𝐵,𝐶𝐸𝑘1𝑚𝑛𝐶𝐵2.(A.13) Considering (A.8), we find: 𝑎223<𝑎22𝑎33,𝑐2<𝑐14𝑟𝐵0𝑘3+𝑘4+𝑏𝑁𝑁,𝐵,𝐶𝐸𝐾𝐵𝐶𝐸𝑘1𝑚𝑏𝐶𝐵𝐾𝑁𝐵𝑟𝑁0𝑁𝐶𝑁/𝐾𝑁𝐵,𝐶𝐸22.(A.14) Considering (A.9), we find an additional condition for choosing our coefficient 𝑐2:𝑎234<𝑎33𝑎44,𝑐2<𝑐34𝑘3+𝑘4+𝑏𝑁𝑁,𝐵,𝐶𝐸𝑘1𝑚𝑛𝑁+𝑔2+𝑏𝐵𝐵,𝐶𝐸𝑘1𝑚𝑏𝐵2.(A.15)

Note that, if both of the conditions relating to 𝑐2 are satisfied, then𝑐2𝑐<min34𝑘3+𝑘4+𝑏𝑁𝑁,𝐵,𝐶𝐸𝑘1𝑚𝑛𝑁+𝑔2+𝑏𝐵𝐵,𝐶𝐸𝑘1𝑚𝑏𝐵2,𝑐14𝑟𝐵0𝑘3+𝑘4+𝑏𝑁𝑁,𝐵,𝐶𝐸𝐾𝐵𝐶𝐸𝑘1𝑚𝑏𝐶𝐵𝐾𝑁𝐵𝑟𝑁0𝑁𝐶𝑁/𝐾𝑁𝐵,𝐶𝐸22.(A.16)

Looking at (A.10), we find the first general condition for local stability: 𝑎213<𝑎11𝑎33,𝑑𝑁𝐶𝑁𝑐2𝑟𝑁0𝐶𝑁𝐾𝑁𝐵,𝐶𝐸2<𝑐24𝑟𝑁0𝑘3+𝑘4+𝑏𝑁𝑁,𝐵,𝐶𝐸𝐾𝑁𝐵,𝐶𝐸.(A.17)

Equation (A.11) gives the second condition for stability: 3𝑎22>𝑎212𝑎11+𝑎224𝑎44+𝑎225𝑎55,3𝑐14𝑟𝐵0𝐾𝐵𝐶𝐸>𝑑𝑁𝐵+𝐾𝑁𝐵𝑟𝑁0𝑁/𝐾𝑁𝐵,𝐶𝐸2𝑐1𝑑𝐵𝑁2𝑟𝑁0/𝐾𝑁𝐵,𝐶𝐸+𝑐1𝑑𝐵𝐶𝐵𝑐3𝑟𝐵0𝐶𝐵/𝐾𝐵𝐶𝐸2𝑐3𝑘1𝑚𝑛𝑁+𝑔2+𝑏𝐵𝐵,𝐶𝐸+𝑐1𝑚𝑒𝐾𝐵𝐶𝐸𝑟𝐵0𝐵/𝐾𝐵𝐶𝐸2𝑚𝑏𝑑𝐵𝑁,𝐶𝐵𝐶𝐵+𝑚𝑛𝑑𝑁𝐵𝐶𝑁𝑁+𝑚𝑏𝑔1𝐶𝐸2𝑚𝑒𝑚𝑏𝑔1𝐵+𝑚𝑛𝑘2𝑁+𝑚𝑒.(A.18) If we let 𝑐1=𝑐3(𝑟𝐵0𝐶𝐵/𝑑𝐵𝐶𝐵𝐾𝐵(𝐶𝐸)), this reduces to3𝑐14𝑟𝐵0𝐾𝐵𝐶𝐸>𝑑𝑁𝐵+𝐾𝑁𝐵𝑟𝑁0𝑁/𝐾𝑁𝐵,𝐶𝐸2𝑐1𝑑𝐵𝑁2𝑟𝑁0/𝐾𝑁𝐵,𝐶𝐸+𝑐1𝑚𝑒𝐾𝐵𝐶𝐸𝑟𝐵0𝐵/𝐾𝐵𝐶𝐸2𝑚𝑏𝑑𝐵𝑁,𝐶𝐵𝐶𝐵+𝑚𝑛𝑑𝑁𝐵𝐶𝑁𝑁+𝑚𝑏𝑔1𝐶𝐸2𝑚𝑒𝑚𝑏𝑔1𝐵+𝑚𝑛𝑘2𝑁+𝑚𝑒.(A.19)

Finally, (A.12) gives the third and final condition for local stability:3𝑎55>𝑎215𝑎11+𝑎235𝑎33+𝑎245𝑎44,3𝑚𝑒4𝑚𝑏𝑔1𝐵+𝑚𝑛𝑘2𝑁+𝑚𝑒>𝑚𝑒𝐾𝑁𝐶𝐸𝑟𝑁0𝑁/𝐾𝑁𝐵,𝐶𝐸2𝑚𝑏𝑑𝐵𝑁𝐶𝐵𝐵𝑚𝑛𝑘3+𝑑𝑁𝐵,𝐶𝑁𝐶𝑁+𝑚𝑛𝑘2𝐶𝐸2𝑟𝑁0/𝐾𝑁𝐵,𝐶𝐸+𝑐2𝑚𝑒𝑘2+𝐾𝑁𝐶𝐸𝑟𝑁0𝑁𝐶𝑁/𝐾𝑁𝐵,𝐶𝐸2+𝑚𝑛𝑁𝑘3+𝑑𝑁(𝐵,𝐶𝑁)+𝑑𝑁𝐶𝑁𝐶𝑁2𝑐2𝑘3+𝑘4+𝑏𝑁𝑁,𝐵,𝐶𝐸+𝑐3𝑚𝑒𝑔1+𝐾𝐵𝐶𝐸𝑟𝐵0𝐵𝐶𝐵/𝐾𝐵𝐶𝐸2+𝑚𝑏𝐵𝑑𝐵(𝑁,𝐶𝐵)+𝑑𝐵𝐶𝐵𝐶𝐵2𝑐3𝑘1𝑚𝑛𝑁