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 𝑑𝑁𝑑𝑑=(𝑏(𝑑)βˆ’π‘‘(𝑑))𝑁,(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 [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 𝑔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π‘βˆ—πΆβˆ—π‘πœ‚1ξ€Έ2𝑐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.

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.

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.

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𝑉𝐢𝑁,𝐡,𝑁,𝐢𝐡,𝐢𝐸=12ξ‚Έ1π‘βˆ—ξ‚€ξπ‘βˆ’π‘βˆ—ξ‚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ξ‚€ξπ΅βˆ’π΅βˆ—ξ‚ξ‚Ά2βˆ’12ξƒ©π‘Ž22βˆ’π‘Ž212π‘Ž11ξƒͺξ‚€ξπ΅βˆ’π΅βˆ—ξ‚2βˆ’12π‘Ž11ξ‚΅(ξπ‘βˆ’π‘βˆ—π‘Ž)+13π‘Ž11ξ‚€ξπΆπ‘βˆ’πΆβˆ—π‘ξ‚ξ‚Ά2βˆ’12ξƒ©π‘Ž33βˆ’π‘Ž213π‘Ž11ξƒͺξ‚€ξπΆπ‘βˆ’πΆβˆ—π‘ξ‚2βˆ’12π‘Ž11ξ‚΅ξ‚€ξπ‘βˆ’π‘βˆ—ξ‚+π‘Ž14π‘Ž11ξ‚€ξπΆπ΅βˆ’πΆβˆ—π΅ξ‚ξ‚Ά2βˆ’12ξƒ©π‘Ž44βˆ’π‘Ž214π‘Ž11ξƒͺξ‚€ξπΆπ΅βˆ’πΆβˆ—π΅ξ‚2βˆ’12π‘Ž11ξ‚΅ξ‚€ξπ‘βˆ’π‘βˆ—ξ‚+π‘Ž15π‘Ž11ξ‚€ξπΆπΈβˆ’πΆβˆ—πΈξ‚ξ‚Ά2βˆ’12ξƒ©π‘Ž55βˆ’π‘Ž215π‘Ž11ξƒͺξ‚€ξπΆπΈβˆ’πΆβˆ—πΈξ‚2βˆ’12π‘Ž22ξ‚΅ξ‚€ξπ΅βˆ’π΅βˆ—ξ‚βˆ’π‘Ž23π‘Ž22ξ‚€ξπΆπ‘βˆ’πΆβˆ—π‘ξ‚ξ‚Ά2βˆ’12ξƒ©π‘Ž33βˆ’π‘Ž223π‘Ž22ξƒͺξ‚€ξπΆπ‘βˆ’πΆβˆ—π‘ξ‚2βˆ’12π‘Ž44ξ‚΅ξ‚€ξπΆπ΅βˆ’πΆβˆ—π΅ξ‚+π‘Ž24π‘Ž44ξ‚€ξπ΅βˆ’π΅βˆ—ξ‚ξ‚Ά2βˆ’12ξƒ©π‘Ž22βˆ’π‘Ž224π‘Ž44ξƒͺξ‚€ξπ΅βˆ’π΅βˆ—ξ‚2βˆ’12π‘Ž55ξ‚΅ξ‚€ξπΆπΈβˆ’πΆβˆ—πΈξ‚+π‘Ž25π‘Ž55ξ‚€ξπ΅βˆ’π΅βˆ—ξ‚ξ‚Ά2βˆ’12ξƒ©π‘Ž22βˆ’π‘Ž225π‘Ž55ξƒͺξ‚€ξπ΅βˆ’π΅βˆ—ξ‚2βˆ’12π‘Ž33ξ‚΅ξ‚€ξπΆπ‘βˆ’πΆβˆ—π‘ξ‚βˆ’π‘Ž34π‘Ž33ξ‚€ξπΆπ΅βˆ’πΆβˆ—π΅ξ‚ξ‚Ά2βˆ’12ξƒ©π‘Ž44βˆ’π‘Ž234π‘Ž33ξƒͺξ‚€ξπΆπ΅βˆ’πΆβˆ—π΅ξ‚2βˆ’12π‘Ž33ξ‚΅ξ‚€ξπΆπ‘βˆ’πΆβˆ—π‘ξ‚βˆ’π‘Ž35π‘Ž33ξ‚€ξπΆπΈβˆ’πΆβˆ—πΈξ‚ξ‚Ά2βˆ’12ξƒ©π‘Ž55βˆ’π‘Ž235π‘Ž33ξƒͺξ‚€ξπΆπΈβˆ’πΆβˆ—πΈξ‚2βˆ’12π‘Ž44ξ‚΅ξ‚€ξπΆπ΅βˆ’πΆβˆ—π΅ξ‚βˆ’π‘Ž45π‘Ž44ξ‚€ξπΆπΈβˆ’πΆβˆ—πΈξ‚ξ‚Ά2βˆ’12ξƒ©π‘Ž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π‘βˆ—πΆβˆ—π‘/πΎπ‘ξ€·π΅βˆ—,πΆβˆ—πΈξ€Έ22.(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π‘βˆ—πΆβˆ—π‘/πΎπ‘ξ€·π΅βˆ—,πΆβˆ—πΈξ€Έ22⎞⎟⎟⎟⎠.(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π‘šπ‘›π‘βˆ—+𝑔2+π‘π΅ξ€·π΅βˆ—,πΆβˆ—πΈ.ξ€Έξ€Έ(A.20)

These conditions are equivalent to those given in (4.1)–(4.6).

B. Proof of Theorem 4.2

Consider the positive definite function 𝑉1𝑁,𝐡,𝐢𝑁,𝐢𝐡,𝐢𝐸=π‘βˆ’π‘βˆ—βˆ’π‘βˆ—ξ‚€π‘lnπ‘βˆ—ξ‚+π΅βˆ’π΅βˆ—βˆ’π΅βˆ—ξ‚€π΅lnπ΅βˆ—ξ‚+𝑐12ξ€·πΆπ‘βˆ’πΆβˆ—π‘ξ€Έ2+𝑐22ξ€·πΆπ΅βˆ’πΆβˆ—π΅ξ€Έ2+12ξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έ2.(B.1)

Then, differentiating 𝑉1 with respect to 𝑑, we get 𝑑𝑉1=ξ€·π‘‘π‘‘π‘βˆ’π‘βˆ—ξ€Έπ‘π‘‘π‘+ξ€·π‘‘π‘‘π΅βˆ’π΅βˆ—ξ€Έπ΅π‘‘π΅π‘‘π‘‘+𝑐1ξ€·πΆπ‘βˆ’πΆβˆ—π‘ξ€Έπ‘‘πΆπ‘π‘‘π‘‘+𝑐2ξ€·πΆπ΅βˆ’πΆβˆ—π΅ξ€Έπ‘‘πΆπ΅+ξ€·πΆπ‘‘π‘‘πΈβˆ’πΆβˆ—πΈξ€Έπ‘‘πΆπΈ.𝑑𝑑(B.2) Using (2.23), we find𝑑𝑉1=ξ€·π‘‘π‘‘π‘βˆ’π‘βˆ—ξ€Έξƒ©π‘π‘0βˆ’π‘Ÿπ‘0𝑁𝐾𝑁𝐡,πΆπΈξ€Έβˆ’π·π‘ξ€·π΅,𝐢𝑁ξƒͺ+ξ€·π΅βˆ’π΅βˆ—ξ€Έξƒ©π‘π΅0βˆ’π‘Ÿπ΅0π΅πΎπ΅ξ€·πΆπΈξ€Έβˆ’π·π΅ξ€·π‘,𝐢𝐡ξƒͺ+𝑐1ξ€·πΆπ‘βˆ’πΆβˆ—π‘ξ€Έξƒ©π‘˜1π‘šπ‘πΆπ΅π΅+π‘˜2πΆπΈβˆ’ξ€·π‘˜3+π‘˜4+𝑏𝑁0𝐢𝑁+π‘Ÿπ‘0𝑁𝐢𝑁𝐾𝑁𝐡,𝐢𝐸ξƒͺ+𝑐2ξ€·πΆπ΅βˆ’πΆβˆ—π΅ξ€Έξƒ©π‘”1πΆπΈβˆ’π‘˜1π‘šπ‘πΆπ΅ξ€·π‘”π‘βˆ’2+𝑏𝐡0𝐢𝐡+π‘Ÿπ΅0𝐡𝐢𝐡𝐾𝐡𝐢𝐸ξƒͺ+ξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έξ‚΅π‘šπ‘’+π‘π‘šπ‘’π·π΅ξ€·π‘,𝐢𝐡𝐢𝐡𝐡+π‘šπ‘›π‘šπ‘’ξ€·π‘˜3+𝐷𝑁𝐡,πΆπ‘πΆξ€Έξ€Έπ‘π‘šπ‘βˆ’π‘π‘šπ‘’π‘”1πΆπΈπ‘šπ΅βˆ’π‘›π‘šπ‘’π‘˜2πΆπΈπ‘βˆ’β„ŽπΆπΈξ‚Ά.(B.3)

Using the notationπœ‚1𝐡,𝐢𝐸=⎧βŽͺβŽͺ⎨βŽͺβŽͺβŽ©ξƒ©1𝐾𝑁𝐡,πΆπΈξ€Έβˆ’1πΎπ‘ξ€·π΅βˆ—,𝐢𝐸ξƒͺ/ξ€·π΅βˆ’π΅βˆ—ξ€Έ,π΅β‰ π΅βˆ—,βˆ’1𝐾2π‘ξ€·π΅βˆ—,πΆπΈξ€Έπœ•πΎπ‘ξ€·π΅βˆ—,πΆπΈξ€Έπœ•π΅,𝐡=π΅βˆ—,πœ‚2ξ€·π΅βˆ—,𝐢𝐸=⎧βŽͺβŽͺ⎨βŽͺβŽͺβŽ©ξƒ©1πΎπ‘ξ€·π΅βˆ—,πΆπΈξ€Έβˆ’1πΎπ‘ξ€·π΅βˆ—,πΆβˆ—πΈξ€Έξƒͺ/ξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έ,πΆπΈβ‰ πΆβˆ—πΈ,βˆ’1𝐾2π‘ξ€·π΅βˆ—,πΆβˆ—πΈξ€Έπœ•πΎπ‘ξ€·π΅βˆ—,πΆβˆ—πΈξ€Έπœ•πΆπΈ,𝐢𝐸=πΆβˆ—πΈ,πœ‡ξ€·πΆπΈξ€Έ=⎧βŽͺβŽͺ⎨βŽͺβŽͺβŽ©ξƒ©1πΎπ΅ξ€·πΆπΈξ€Έβˆ’1πΎπ΅ξ€·πΆβˆ—πΈξ€Έξƒͺ/ξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έ,πΆπΈβ‰ πΆβˆ—πΈ,βˆ’1𝐾2π΅ξ€·πΆβˆ—πΈξ€Έπœ•πΎπ΅ξ€·πΆβˆ—πΈξ€Έπœ•πΆπΈ,𝐢𝐸=πΆβˆ—πΈ,πœ“1𝐡,𝐢𝑁=⎧βŽͺ⎨βŽͺβŽ©π·π‘ξ€·π΅,πΆπ‘ξ€Έβˆ’π·π‘ξ€·π΅βˆ—,πΆπ‘ξ€Έπ΅βˆ’π΅βˆ—,π΅β‰ π΅βˆ—,πœ•π·π‘ξ€·π΅,πΆπ‘ξ€Έπœ•π΅,𝐡=π΅βˆ—,πœ“2ξ€·π΅βˆ—,𝐢𝑁=⎧βŽͺβŽͺ⎨βŽͺβŽͺβŽ©π·π‘ξ€·π΅βˆ—,πΆπ‘ξ€Έβˆ’π·π‘ξ€·π΅βˆ—,πΆβˆ—π‘ξ€ΈπΆπ‘βˆ’πΆβˆ—π‘,πΆπ‘β‰ πΆβˆ—π‘,πœ•π·π‘ξ€·π΅βˆ—,πΆπ‘ξ€Έπœ•πΆπ‘,𝐢𝑁=πΆβˆ—π‘,𝜁1𝑁,𝐢𝐡=⎧βŽͺ⎨βŽͺβŽ©π·π΅ξ€·π‘,πΆπ΅ξ€Έβˆ’π·π΅ξ€·π‘βˆ—,πΆπ΅ξ€Έπ‘βˆ’π‘βˆ—,π‘β‰ π‘βˆ—,πœ•π·π΅ξ€·π‘,πΆπ΅ξ€Έπœ•π‘,𝑁=π‘βˆ—,𝜁2ξ€·π‘βˆ—,𝐢𝐡=⎧βŽͺβŽͺ⎨βŽͺβŽͺβŽ©π·π΅ξ€·π‘βˆ—,πΆπ΅ξ€Έβˆ’π·π΅ξ€·π‘βˆ—,πΆβˆ—π΅ξ€ΈπΆπ΅βˆ’πΆβˆ—π΅,πΆπ΅β‰ πΆβˆ—π΅,πœ•π·π΅ξ€·π‘βˆ—,πΆπ΅ξ€Έπœ•πΆπ΅,𝐢𝐡=πΆβˆ—π΅(B.4)

and some algebraic manipulation, we can rewite the systems as𝑑𝑉11𝑑𝑑=βˆ’2π‘Ž11ξ€·π‘βˆ’π‘βˆ—ξ€Έ2+π‘Ž12ξ€·π‘βˆ’π‘βˆ—ξ€Έξ€·π΅βˆ’π΅βˆ—ξ€Έβˆ’12π‘Ž22ξ€·π΅βˆ’π΅βˆ—ξ€Έ2βˆ’12π‘Ž11ξ€·π‘βˆ’π‘βˆ—ξ€Έ2+π‘Ž13ξ€·π‘βˆ’π‘βˆ—πΆξ€Έξ€·π‘βˆ’πΆβˆ—π‘ξ€Έβˆ’12π‘Ž33ξ€·πΆπ‘βˆ’πΆβˆ—π‘ξ€Έ2βˆ’12π‘Ž11ξ€·π‘βˆ’π‘βˆ—ξ€Έ2+π‘Ž14ξ€·π‘βˆ’π‘βˆ—πΆξ€Έξ€·π΅βˆ’πΆβˆ—π΅ξ€Έβˆ’12π‘Ž44ξ€·πΆπ΅βˆ’πΆβˆ—π΅ξ€Έ2βˆ’12π‘Ž11ξ€·π‘βˆ’π‘βˆ—ξ€Έ2+π‘Ž15ξ€·π‘βˆ’π‘βˆ—πΆξ€Έξ€·πΈβˆ’πΆβˆ—πΈξ€Έβˆ’12π‘Ž55ξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έ2βˆ’12π‘Ž22ξ€·π΅βˆ’π΅βˆ—ξ€Έ2+π‘Ž23ξ€·π΅βˆ’π΅βˆ—πΆξ€Έξ€·π‘βˆ’πΆβˆ—π‘ξ€Έβˆ’12π‘Ž33ξ€·πΆπ‘βˆ’πΆβˆ—π‘ξ€Έ2βˆ’12π‘Ž22ξ€·π΅βˆ’π΅βˆ—ξ€Έ2+π‘Ž24ξ€·π΅βˆ’π΅βˆ—πΆξ€Έξ€·π΅βˆ’πΆβˆ—π΅ξ€Έβˆ’12π‘Ž44ξ€·πΆπ΅βˆ’πΆβˆ—π΅ξ€Έ2βˆ’12π‘Ž22ξ€·π΅βˆ’π΅βˆ—ξ€Έ2+π‘Ž25ξ€·π΅βˆ’π΅βˆ—πΆξ€Έξ€·πΈβˆ’πΆβˆ—πΈξ€Έβˆ’12π‘Ž55ξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έ2βˆ’12π‘Ž33ξ€·πΆπ‘βˆ’πΆβˆ—π‘ξ€Έ2+π‘Ž34ξ€·πΆπ‘βˆ’πΆβˆ—π‘πΆξ€Έξ€·π΅βˆ’πΆβˆ—π΅ξ€Έβˆ’12π‘Ž44ξ€·πΆπ΅βˆ’πΆβˆ—π΅ξ€Έ2βˆ’12π‘Ž33ξ€·πΆπ‘βˆ’πΆβˆ—π‘ξ€Έ2+π‘Ž35ξ€·πΆπ‘βˆ’πΆβˆ—π‘πΆξ€Έξ€·πΈβˆ’πΆβˆ—πΈξ€Έβˆ’12π‘Ž55ξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έ2βˆ’12π‘Ž44ξ€·πΆπ΅βˆ’πΆβˆ—π΅ξ€Έ2+π‘Ž45ξ€·πΆπ΅βˆ’πΆβˆ—π΅πΆξ€Έξ€·πΈβˆ’πΆβˆ—πΈξ€Έβˆ’12π‘Ž55ξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έ2,(B.5) where π‘Ž11=π‘Ÿπ‘02𝐾𝑁𝐡,𝐢𝐸,π‘Ž22=π‘Ÿπ΅02𝐾𝐡𝐢𝐸,π‘Ž33=𝑐12ξƒ©π‘˜3+π‘˜4+𝑏𝑁0βˆ’π‘Ÿπ‘0π‘βˆ—πΎπ‘ξ€·π΅,𝐢𝐸ξƒͺ,π‘Ž44=𝑐22𝑔2+π‘˜1π‘šπ‘›π‘βˆ—+𝑏𝐡0βˆ’π‘Ÿπ΅0π΅βˆ—πΎπ΅ξ€·πΆπΈξ€Έξƒͺ,π‘Ž55=12ξ‚΅π‘šβ„Ž+π‘π‘šπ‘’π‘”1π΅βˆ—+π‘šπ‘›π‘šπ‘’π‘˜2π‘βˆ—ξ‚Ά,π‘Ž12=βˆ’πœ“1𝐡,πΆπ‘ξ€Έβˆ’π‘βˆ—π‘Ÿπ‘0πœ‚1𝐡,πΆπΈξ€Έβˆ’πœ1𝑁,𝐢𝐡,π‘Ž13=βˆ’πœ“2ξ€·π΅βˆ—,𝐢𝑁+𝑐1π‘Ÿπ‘0𝐢𝑁𝐾𝑁𝐡,𝐢𝐸,π‘Ž14=βˆ’π‘2π‘˜1π‘šπ‘›πΆπ΅,π‘Ž15=βˆ’π‘βˆ—π‘Ÿπ‘0πœ‚2ξ€·π΅βˆ—,𝐢𝐸+ξ‚€π‘šπ‘›ξ€·π‘˜π‘šπ‘’3+𝐷𝑁𝐡,πΆπ‘πΆξ€Έξ€Έπ‘βˆ’π‘šπ‘›π‘˜π‘šπ‘’2𝐢𝐸+π‘šπ‘π΅π‘šπ‘’βˆ—πΆβˆ—π΅πœ1𝑁,𝐢𝐡,π‘Ž23=𝑐1ξ€·π‘˜1π‘šπ‘πΆπ΅+π‘Ÿπ‘0π‘βˆ—πΆβˆ—π‘πœ‚1𝐡,𝐢𝐸,π‘Žξ€Έξ€Έ24=βˆ’πœ2ξ€·π‘βˆ—,𝐢𝐡+𝑐2π‘Ÿπ΅0𝐢𝐡𝐾𝐡𝐢𝐸,π‘Ž25=βˆ’π΅βˆ—π‘Ÿπ΅0πœ‡ξ€·πΆπΈξ€Έ+ξ‚΅π‘šπ‘π‘šπ‘’π·π΅ξ€·π‘,πΆπ΅ξ€ΈπΆπ΅βˆ’π‘šπ‘π‘šπ‘’π‘”1𝐢𝐸+π‘šπ‘›π‘šπ‘’π‘βˆ—πΆβˆ—π‘πœ“1𝐡,𝐢𝑁,π‘Ž34=𝑐1π‘˜1π‘šπ‘π΅βˆ—,π‘Ž35=𝑐1ξ€·π‘˜2+π‘Ÿπ‘0π‘βˆ—πΆβˆ—π‘πœ‚2ξ€·π΅βˆ—,𝐢𝐸+π‘šξ€Έξ€Έπ‘›π‘šπ‘’π‘βˆ—ξ€·π‘˜3+𝐷𝑁𝐡,𝐢𝑁+πΆβˆ—π‘πœ“2ξ€·π΅βˆ—,𝐢𝑁,π‘Žξ€Έξ€Έ45=𝑐2𝑔1+π‘Ÿπ‘0π΅βˆ—πΆβˆ—π΅πœ‡ξ€·πΆπΈ+π‘šξ€Έξ€Έπ‘π‘šπ‘’π΅βˆ—ξ€·π·π΅ξ€·π‘,𝐢𝑁+πΆβˆ—π΅πœ2ξ€·π‘βˆ—,𝐢𝐡.ξ€Έξ€Έ(B.6)

Note that||||πœ•πΎπ‘ξ€·π΅,𝐢𝐸||||πœ•π΅=𝐾𝑁𝐡,||||πœ•πΎπ‘ξ€·π΅,πΆπΈξ€Έπœ•πΆπΈ||||=𝐾𝑁𝐢𝐸,||||πœ•πΎπ΅ξ€·πΆπΈξ€Έπœ•πΆπΈ||||=𝐾𝐡𝐢𝐸,||||πœ•π·π‘ξ€·π΅,𝐢𝑁||||πœ•π΅=𝑑𝑁𝐡,||||πœ•π·π‘ξ€·π΅,πΆπ‘ξ€Έπœ•πΆπ‘||||=𝑑𝑁𝐢𝑁,||||πœ•π·π΅ξ€·π‘,𝐢𝐡||||πœ•π‘=𝑑𝐡𝑁,||||πœ•π·π΅ξ€·π‘,πΆπ΅ξ€Έπœ•πΆπ΅||||=𝑑𝐡𝐢𝐡.(B.7)

Additionally, note that, by mean value theorem and the previous definitions, ||πœ‚1||≀𝐾𝑁𝐡𝐾𝑁0+π΅βˆ—πΎπ‘π΅βˆ’πΎπ‘πΆπΈξ€Έ2=πœ‚1,||πœ‚2||≀𝐾𝑁𝐢𝐸𝐾𝑁0+π΅βˆ—πΎπ‘π΅βˆ’πΎπ‘πΆπΈπΆβˆ—πΈξ€Έ2=πœ‚2,||πœ‡||≀𝐾𝐡𝐢𝐸𝐾𝐡0βˆ’πΎπ΅πΆπΈπΆβˆ—πΈξ€Έ2=||πœ“πœ‡,1||≀𝑑𝑁𝐡,||πœ“2||≀𝑑𝑁𝐢𝑁,||𝜁1||≀𝑑𝐡𝑁,||𝜁2||≀𝑑𝐡𝐢𝐡.(B.8) Factoring 𝑑𝑉1/𝑑𝑑 as 𝑑𝑉11𝑑𝑑=βˆ’2π‘Ž11ξ‚΅ξ€·π‘βˆ’π‘βˆ—ξ€Έβˆ’π‘Ž12π‘Ž11ξ€·π΅βˆ’π΅βˆ—ξ€Έξ‚Ά2βˆ’12ξƒ©π‘Ž22βˆ’π‘Ž212π‘Ž11ξƒͺξ€·π΅βˆ’π΅βˆ—ξ€Έ2βˆ’12π‘Ž11ξ‚΅ξ€·π‘βˆ’π‘βˆ—ξ€Έβˆ’π‘Ž13π‘Ž11ξ€·πΆπ‘βˆ’πΆβˆ—π‘ξ€Έξ‚Ά2βˆ’12ξƒ©π‘Ž33βˆ’π‘Ž213π‘Ž11ξƒͺξ€·πΆπ‘βˆ’πΆβˆ—π‘ξ€Έ2βˆ’12π‘Ž11ξ‚΅ξ€·π‘βˆ’π‘βˆ—ξ€Έβˆ’π‘Ž14π‘Ž11ξ€·πΆπ΅βˆ’πΆβˆ—π΅ξ€Έξ‚Ά2βˆ’12ξƒ©π‘Ž44βˆ’π‘Ž214π‘Ž11ξƒͺξ€·πΆπ΅βˆ’πΆβˆ—π΅ξ€Έ2βˆ’12π‘Ž11ξ‚΅ξ€·π‘βˆ’π‘βˆ—ξ€Έβˆ’π‘Ž15π‘Ž11ξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έξ‚Ά2βˆ’12ξƒ©π‘Ž55βˆ’π‘Ž215π‘Ž11ξƒͺξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έ2βˆ’12π‘Ž22ξ‚΅ξ€·π΅βˆ’π΅βˆ—ξ€Έβˆ’π‘Ž23π‘Ž22ξ€·πΆπ‘βˆ’πΆβˆ—π‘ξ€Έξ‚Ά2βˆ’12ξƒ©π‘Ž33βˆ’π‘Ž223π‘Ž22ξƒͺξ€·πΆπ‘βˆ’πΆβˆ—π‘ξ€Έ2βˆ’12π‘Ž44ξ‚΅ξ€·πΆπ΅βˆ’πΆβˆ—π΅ξ€Έβˆ’π‘Ž24π‘Ž44ξ€·π΅βˆ’π΅βˆ—ξ€Έξ‚Ά2βˆ’12ξƒ©π‘Ž22βˆ’π‘Ž224π‘Ž44ξƒͺξ€·π΅βˆ’π΅βˆ—ξ€Έ2βˆ’12π‘Ž55ξ‚΅ξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έβˆ’π‘Ž25π‘Ž55ξ€·π΅βˆ’π΅βˆ—ξ€Έξ‚Ά2βˆ’12ξƒ©π‘Ž22βˆ’π‘Ž225π‘Ž55ξƒͺξ€·π΅βˆ’π΅βˆ—ξ€Έ2βˆ’12π‘Ž33ξ‚΅(πΆπ‘βˆ’πΆβˆ—π‘π‘Ž)βˆ’34π‘Ž33(πΆπ΅βˆ’πΆβˆ—π΅)ξ‚Ά2βˆ’12ξƒ©π‘Ž44βˆ’π‘Ž234π‘Ž33ξƒͺξ€·πΆπ΅βˆ’πΆβˆ—π΅ξ€Έ2βˆ’12π‘Ž33ξ‚΅ξ€·πΆπ‘βˆ’πΆβˆ—π‘ξ€Έβˆ’π‘Ž35π‘Ž33ξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έξ‚Ά2βˆ’12ξƒ©π‘Ž55βˆ’π‘Ž235π‘Ž33ξƒͺξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έ2βˆ’12π‘Ž44ξ‚΅ξ€·πΆπ΅βˆ’πΆβˆ—π΅ξ€Έβˆ’π‘Ž45π‘Ž44ξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έξ‚Ά2βˆ’12ξƒ©π‘Ž55βˆ’π‘Ž245π‘Ž44ξƒͺξ€·πΆπΈβˆ’πΆβˆ—πΈξ€Έ2,(B.9) we see that sufficient conditions for 𝑑𝑉1/𝑑𝑑 to be negative definite are: π‘Ž214<π‘Ž11π‘Ž44,π‘Ž(B.10)223<π‘Ž22π‘Ž33π‘Ž,(B.11)234<π‘Ž33π‘Ž44π‘Ž,(B.12)213<π‘Ž11π‘Ž33π‘Ž,(B.13)212π‘Ž11+π‘Ž224π‘Ž44+π‘Ž225π‘Ž55<3π‘Ž22π‘Ž,(B.14)215π‘Ž11+π‘Ž235π‘Ž33+π‘Ž245π‘Ž44<3π‘Ž55.(B.15) To simplify the conditions, we can find an upper bound for the left-hand sides of the inequalities and lower bounds for the right-hand sides of the inequalities.

First consider (B.10). This can be used to place a restriction on 𝑐2: π‘Ž214<π‘Ž11π‘Ž44,𝑐2<π‘Ÿπ‘0𝑔2+π‘˜1π‘šπ‘›π‘βˆ—ξ€Έ4πΎπ‘ξ€·π‘˜1π‘šπ‘›ξ€Έ2.(B.16)

Equation (B.11) can be used to place a restriction on 𝑐1: π‘Ž223<π‘Ž22π‘Ž33,𝑐1<π‘Ÿπ΅0ξ€·π‘˜3+π‘˜4ξ€Έ4πΎπ΅ξ€·π‘˜1π‘šπ‘+π‘Ÿπ‘0π‘βˆ—πΆβˆ—π‘πœ‚1ξ€Έ2.(B.17)

Considering (B.12), an additional restriction can be placed on 𝑐1: π‘Ž234<π‘Ž33π‘Ž44,𝑐1<𝑐2ξ€·π‘˜3+π‘˜4𝑔2+π‘˜1π‘šπ‘›π‘βˆ—ξ€Έ4ξ€·π‘˜1π‘šπ‘π΅βˆ—ξ€Έ2.(B.18)

Since we have established 𝑐1<π‘Ÿπ΅0ξ€·π‘˜3+π‘˜4ξ€Έ4πΎπ΅ξ€·π‘˜1π‘šπ‘+π‘Ÿπ‘0π‘βˆ—πΆβˆ—π‘πœ‚1ξ€Έ2,𝑐1<𝑐2ξ€·π‘˜3+π‘˜4𝑔2+π‘˜1π‘šπ‘›π‘βˆ—ξ€Έ4ξ€·π‘˜1π‘šπ‘π΅βˆ—ξ€Έ2,(B.19) both conditions will be satisfied if 𝑐1ξƒ©π‘Ÿ<min𝐡0ξ€·π‘˜3+π‘˜4ξ€Έ4πΎπ΅ξ€·π‘˜1π‘šπ‘+π‘Ÿπ‘0π‘βˆ—πΆβˆ—π‘πœ‚1ξ€Έ2,𝑐2ξ€·π‘˜3+π‘˜4𝑔2+π‘˜1π‘šπ‘›π‘βˆ—ξ€Έ4ξ€·π‘˜1π‘šπ‘π΅βˆ—ξ€Έ2ξƒͺ.(B.20)

Next consider (B.13) to establish the first condition for global stability: π‘Ž213<π‘Ž11π‘Ž33,𝑑𝑁𝐢𝑁+𝑐1π‘Ÿπ‘0𝐾𝑁ξƒͺ2<𝑐1π‘Ÿπ‘0ξ€·π‘˜3+π‘˜4ξ€Έ4𝐾𝑁.(B.21)

Considering (B.14), we develop the second sufficient condition for global stability: 3π‘Ž22>π‘Ž212π‘Ž11+π‘Ž224π‘Ž44+π‘Ž225π‘Ž5534π‘Ÿπ΅0𝐾𝐡>𝐾𝑁𝑑𝑁𝐡+π‘βˆ—π‘Ÿπ‘0πœ‚1+𝑑𝐡𝑁2π‘Ÿπ‘0+𝑑𝐡𝐢𝐡+𝑐2π‘Ÿπ΅0/𝐾𝐡2𝑐2𝑔2+π‘˜1π‘šπ‘›π‘βˆ—ξ€Έ+ξ‚€π‘šπ‘’π΅βˆ—π‘Ÿπ΅0πœ‡+π‘šπ‘π·π΅+π‘šπ‘π‘”1+π‘šπ‘›π‘βˆ—πΆβˆ—π‘π‘‘π‘π΅ξ‚2π‘šπ‘’ξ€·π‘šπ‘’β„Ž+π‘šπ‘π‘”1π΅βˆ—+π‘šπ‘›π‘˜2π‘βˆ—ξ€Έ.(B.22)

Finally, considering (B.15), 3π‘Ž55>π‘Ž215π‘Ž11+π‘Ž235π‘Ž33+π‘Ž245π‘Ž44,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π‘šπ‘›π‘βˆ—ξ€Έ.(B.23) These conditions are equivalent to the conditions for global stability given in Theorem 4.2.

Acknowledgment

This work was supported by the National Science and Engineering Research Council of Canada.