# Dynamics of the Oxygen, Carbon Dioxide, and Water Interaction across the Insect Spiracle

Authors: S. M. Simelane; S. Abelman; F. D. Duncan
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This paper explores the dynamics of respiratory gases interactions which are accompanied by the loss of water through an insect’s spiracle. Here we investigate and analyze this interaction by deriving a system of ordinary differential equations for oxygen, carbon dioxide, and water vapor. The analysis is carried out in continuous time. The purpose of the research is to determine bounds for the gas volumes and to discuss the complexity and stability of the equilibria. Numerical simulations also demonstrate the dynamics of our model utilizing the new conditions for stability and instability.


---

## Body

## 1. Introduction

The Lotka-Volterra equations have been used as a basis for many types of models that involve the interactions of different populations (see [1–3]), mostly the interaction of two populations. Other extensions and modifications such as diffusion, ratio-dependence, and functional response have been added to the Lotka-Volterra equations to gain an even better and deeper understanding of the dynamics of population interactions [2, 4, 5]. The model discussed in this paper is also based on the simple Lotka-Volterra model developed by Månsson and Lundberg [6]. In predator-prey models, both the predator and the prey are being born; the predator feeds on its prey and the predator dies at some point. In a sense O
2 and C
O
2 are being born and die. O
2 is born when it enters the trachea and starts to die as it is used by the cells. O
2 is replaced by C
O
2. In the same way if there is sufficient prey the predator population expands at the expense of the prey. So C
O
2 increases as O
2 is consumed and the point is reached when the populations need to be stable; this is when the spiracle opens and we now have the cycle reversing in that at this point O
2 enters the trachea, is used up by the cells and is then replaced by C
O
2. Just the same way as the classical predator-prey model has an increase in prey then a delayed increase in predators which causes a decrease in prey and then a decrease in the predators and it cycles, so too do the O
2 and C
O
2 in the insect trachea. The motivation behind this is that the model presented here is a competing species model in that the respiratory gases, oxygen, and carbon dioxide are competing for space in the trachea and there is a maximum spiracle area through which the gases can pass. When the spiracle opens for gas exchange, the gases manage to move in opposite directions simultaneously. Thus on the opening of the spiracle, O
2 needs to be replaced and C
O
2 needs to exit through the same area. This analysis prompted us to consider any existing competing species model to try and explain the gases movement. The simple competing predator-prey model provides a starting point for our argument.Studies of gas exchange through the insect tracheal system have been carried out in detail (see [7, 8] for reviews) focusing on partial pressure differences as the driving force for gas exchange. Individuals attempting to quantify gas exchange in insects have been divided into two groups according to Snyder et al. [9]. One group interested in insect physiology focused on gas exchange and the relationship between gas exchange and water balance [8–11]. Several models have been developed for predicting respiratory water loss as a consequence of the opening of the spiracles (see [12, 13]). The other group interested in respiratory physiology focused on the sites of resistance to gas exchange, attempting to determine as to where resistance lies and if gas exchange occurs entirely by diffusion or convection [7, 9, 14].Cockroaches and locusts inhale air through the anterior spiracles and exhale through the more posterior spiracles [15]. Thus the coordinated action of spiracles can direct air flow through an insect. Insects are also known to close their spiracles for periods of time resulting in no gas exchange taking place. The intermittent opening and closing of spiracles is referred to as discontinuous gas exchange (DGC) [16]. The species used for this study is a flightless dung beetle, which uses a single spiracle for gas exchange at rest [17].The purpose of this paper is to analyze and demonstrate the dynamics of various gases in our model system by(a)
determining ultimate bounds for the interacting gas volumes,
(b)
exploring the stability and instability of the equilibrium equations,
(c)
obtaining numerical simulations for the pattern of the dynamics for the interacting gases.In our research, we present a mathematical model that captures how the respiratory gases manage to move in opposite directions simultaneously. The model should be applicable to terrestrial insects under any environmental conditions to trace the DGC, cyclic, and continuous gas flux. It should answer questions such as the following: how do respiratory gases manage to move in opposite directions simultaneously?, How does extended spiracle opening influence the influx of the gases? and What factors are responsible for instability of the system?. To answer all the above questions, we need to study the dynamical behavior of the system and test it with experimental data.In Section2, the discontinuous gas exchange across the insect spiracles is summarized. In Section 3, the model assumptions and the system model are presented. Primary results for the model and boundedness properties for the respiratory gases are given in Section 4. The analysis of equilibrium solutions is also undertaken along with conditions that ensure nonnegativity of the equilibrium solutions in Section 5. Asymptotic stability and instability for the five equilibrium solutions are determined through linearization. In Section 6, graphical display for the dynamics and pattern of the gas volumes is shown by utilizing the conditions for stability and instability. Finally in Section 7, we summarize our findings and explain how the conclusions are made.


## 2. Gas Exchange in Insects

Many beetles living in arid areas are flightless and have a subelytral cavity. The dung beetle is wingless and has a sealed airspace called the subelytral cavity beneath the wing covers (Figure1). Dung beetles have only one pair of spiracles in the head region and have 7 pairs of spiracles in the body region. The 7 pairs of spiracles open into the subelytral cavity and are all positioned under the elytra. The spiracles are the sites of gas exchange.Summary figure showing two beetles. (a) Shows a dung beetle (photo: F.D. Duncan) and (b) shows the typical morphology of a beetle (see [24]).
(a)
(b)Gas exchange in insects is a continual conflict between the need to support aerobic need and the need to reduce respiratory water loss [7, 9, 11]. The respiratory system of insects consists of a highly branched system of cuticle-lined tubes extending throughout the body. The tubes are filled with air, which greatly facilitates the transport of gases through the body [18]. The trachea opens to the external atmosphere through valve-like spiracles on the surface of the abdomen and thorax. The principal sites for gas exchange are the fine branches (less than a micrometer in diameter) at the tip of the tracheal system, internal to the spiracles [18]. Insects lose water through the cuticle and through the open spiracles, with respiratory water loss contributing a small proportion of the overall water loss in most insect species [16].Many insects exchange respiratory gases with the atmosphere discontinuously, thus permitting gas exchange while minimizing water loss despite their relatively high metabolic rates. While oxygen consumption and carbon dioxide production are continuous, the external elimination of gases is discontinuous [9, 18]. The classical pattern of DGC consists of three phases [16, 19, 20], also observed in the beetle used for this study.


### 2.1. Closed Spiracle Phase (C-Phase)

This phase begins when the spiracles are tightly closed and the insect is essentially hermetically sealed [7, 9]. Respiration occurs here and is driven by the oxygen stores in the cells. The partial pressure of oxygen (P
O
2) within the trachea declines and the partial pressure of carbon dioxide (P
C
O
2) increases. Owing to the high solubility of C
O
2, the increase in P
C
O
2 in the air within the trachea is less than the decline in P
O
2. As a result, total atmospheric pressure in the tracheal lumina declines slightly during the closed phase [18].


### 2.2.  Flutter Phase (F-Phase)

The spiracles begin to flutter, allowing large quantities of air down this pressure gradient. Inward convection during theF-phase is an important means of restricting outward water movement [20, 21]. Here, convective influx of oxygen rich air into the tracheal system allows O
2 uptake with minimal water loss. C
O
2 continues to accumulate throughout the F-phase period.


### 2.3. Open Spiracle Phase (O-Phase)

During this phase, the spiracles remain fully opened and the gases are exchanged with the atmosphere. Oxygen gain and carbon dioxide loss happen during this time [22]. P
O
2 and P
C
O
2 tend towards levels present in the external atmosphere due to the diffusion of these gases down their concentration gradients through the trachea and open spiracles [18, 23]. This phase is important because water loss is high as long as the spiracles remain open [7–9].In the following sections, we develop and analyze a model for the respiratory gases dynamics that considers the discontinuous exchange of these gases. The model includes all of the above features.


## 2.1. Closed Spiracle Phase (C-Phase)

This phase begins when the spiracles are tightly closed and the insect is essentially hermetically sealed [7, 9]. Respiration occurs here and is driven by the oxygen stores in the cells. The partial pressure of oxygen (P
O
2) within the trachea declines and the partial pressure of carbon dioxide (P
C
O
2) increases. Owing to the high solubility of C
O
2, the increase in P
C
O
2 in the air within the trachea is less than the decline in P
O
2. As a result, total atmospheric pressure in the tracheal lumina declines slightly during the closed phase [18].


## 2.2.  Flutter Phase (F-Phase)

The spiracles begin to flutter, allowing large quantities of air down this pressure gradient. Inward convection during theF-phase is an important means of restricting outward water movement [20, 21]. Here, convective influx of oxygen rich air into the tracheal system allows O
2 uptake with minimal water loss. C
O
2 continues to accumulate throughout the F-phase period.


## 2.3. Open Spiracle Phase (O-Phase)

During this phase, the spiracles remain fully opened and the gases are exchanged with the atmosphere. Oxygen gain and carbon dioxide loss happen during this time [22]. P
O
2 and P
C
O
2 tend towards levels present in the external atmosphere due to the diffusion of these gases down their concentration gradients through the trachea and open spiracles [18, 23]. This phase is important because water loss is high as long as the spiracles remain open [7–9].In the following sections, we develop and analyze a model for the respiratory gases dynamics that considers the discontinuous exchange of these gases. The model includes all of the above features.


## 3. The Model System

Internally within an insect, oxygen is taken into the cells and through the process of cellular respiration carbon dioxide is produced as a waste product. In insects, thisC
O
2 is stored in the tissues and haemolymph (similar to blood in vertebrates). When this storage is saturated the C
O
2 enters the trachea (hollow tubes) in order to be expelled to the exterior of the insect. Air enters and leaves the insect through the spiracles, which are the opening of the trachea tube to the exterior.In insects using DGC the spiracles are closed while theO
2 is being used by the cells and C
O
2 is stored. When there is sufficient build-up of C
O
2, this increase in C
O
2 concentration causes the spiracles to open. Upon opening C
O
2 rushes out of the tracheal system and O
2 rushes in. Both these gases and water vapor have to leave or enter through the small spiracle opening. Thus these gases are competing for space [7].The model presented looks at the movement of the respiratory gases across the spiracle opening, across the tracheal system. We consider the standard predator-prey model which consists of two nonlinear ordinary differential equations (see [1–6] for reviews). We extend this model to include a third nonlinear differential equation, hence a system of three nonlinear ordinary differential equations. We assume there are three components. (i) Oxygen denoted by O
2, (ii) Carbon dioxide denoted by C
O
2. (iii) Water vapor denoted by H
2
O.Predator-prey models are an essential tool in mathematical modeling and for our understanding of interacting species in the natural environment. A better description of this model is a competing species model. Thus throughout the paper we will refer to a competing species model. Traditional predator-prey models are based on the principle of mass action. The mass action principle predicts that the rate at which an individual predator consumes prey should depend only on prey density. In applying the mass action principle there is also the implicit assumption that the system is well mixed so that inhomogeneities do not develop. In developing our model of the gas interaction, the principle of mass action was the main driving principle. The rate at which respiring cells take in oxygen depends on the oxygen density in the trachea. This parallel analysis prompted the use of basic predator-prey models in an attempt to trace the movement of the respiratory gases. Assumptions in building the competing species model are as follows:(1)
if one gas is moving, the rate of change of the gas follows a logistic model (model with tracheal size threshold);
(2)
the rate of change of each gas decreases at a rate proportional to the amount of gas flow (competition).These assumptions are reasonable because when the spiracle opens, the rate of change of oxygen growth is an accelerating function of oxygen until atmospheric pressure is reached. Again, the rate of change of oxygen once the spiracle closes decreases at a rate proportional to the amount of oxygen used up for respiration and thus carbon dioxide production.Many models that trace the movement of respiratory gases exist in the literature. These models are based on both Fick’s laws of diffusion, which use the fact that most insects use both convective and diffusive exchanges. For instance, Welch and Tracy [12] developed a model of convective gas flux as a function of the total gas, the inspired, and expired fractions of the gases. Even though the model is widely accepted, it ignores known complications to insect physiologists. This includes the effects of end diffusion, interference between adjacent respiratory tubes, and interactions between water vapor leaving and gas entering through the tubes. It also fails to address the fact that carbon dioxide elimination takes longer than the ingress of oxygen.The failure of the existing models to provide generality and unification because of their taxon-specific design prompted us to develop a novel general model. The assumptions in our model derivation are analogous to those in predator-prey models. With these assumptions, we develop a ratio-dependent model and derive a model for oxygen-carbon dioxide interactions with ratio-dependent interaction rate. Ratio-dependent functional response is a better starting point for modeling population dynamics. There is a legitimate predictive use of ratio-dependent models in the literature (see [25–28] for reviews). Evidence suggests that even when functional responses are measured in the short term, they are closer to ratio-dependent than prey-dependent. This is supported by the argument that on a periodic time scale, the surrounding cells limit each other’s oxygen supply through sharing (i.e., when oxygen is in short supply during the C-phase). Berryman [29] developed a ratio-dependent model by extending the logistic growth model. The new model is developed by extending Berryman’s [29] model (see [6, 29–31]). Here, the oxygen intake by the respiring cells, i, is described by
(1)
i
=
a
1
+
b
(
C
O
2
/
O
2
)
,
where C
O
2 is the carbon dioxide volume, O
2 is the oxygen volume, a is the amount of O
2 taken in by the respiring cell when O
2 is in excess, and b is the saturating factor deciding how quickly the respiring cells intake decelerates as oxygen density decreases in the trachea.It follows from (1) that the risk of being used by a respiring cell for each unit of oxygen is
(2)
a
(
C
O
2
/
O
2
)
1
+
b
(
C
O
2
/
O
2
)
,
where C
O
2
/
O
2 is the respiratory quotient.Now, we assume that the carbon dioxide elimination rate is a function of oxygen intake rate,i. Thus, the C
O
2 production and elimination functions are, respectively,
(3)
c
i
1
+
g
i
,
d
maxC
O
2
-
δ
C
O
2
i
,
where c measures how the cells turn the inspired oxygen into carbon dioxide, g measures how carbon dioxide production is physiologically constrained by these cells, d
maxC
O
2 measures carbon dioxide maximum per capita elimination rate, and δ
C
O
2 adjusts carbon dioxide elimination rate proportionality to the respiring cell intake rate. We model the cells carbon dioxide production rate as a nonlinear increasing function of oxygen intake rate while the elimination rate is a linearly decreasing function. The nonlinearity of the carbon dioxide production function reflects a physiological constraint on the rate of new carbon dioxide. It involves trade-offs with, for instance, cell growth and is likely to be limited since the quantities of carbon dioxide are finite. Note that the C
O
2 elimination could be modeled as any kind of decreasing function of intake, but for simplicity we used a linear function.The following differential equations model the dynamics of the gas volumes:(4)
d
O
2
d
t
=
u
O
2
(
1
-
O
2
+
C
O
2
K
)
-
a
O
2
C
O
2
O
2
+
b
C
O
2
,
d
H
2
O
d
t
=
f
H
2
O
-
e
C
O
2
H
2
O
H
2
O
+
h
-
s
H
2
O
,
d
C
O
2
d
t
=
a
c
O
2
C
O
2
O
2
(
1
+
g
a
)
+
b
C
O
2
-
d
maxC
O
2
C
O
2
+
δ
C
O
2
a
O
2
C
O
2
O
2
+
b
C
O
2
.
with u, a, b, c, e, f, g, s, and K
>
0.O
2, C
O
2, and H
2
O denote the gas volumes. u is the oxygen growth rate at low oxygen densities, K is the tracheal carrying capacity of the oxygen gas, e is the maximum spiracular water loss in the thorax and elytral case as a result of the insect respiration, f is the water ingress rate into the trachea, and h is the water density/volume where water loss is half saturated. Lastly, s is the rate of water loss due to other forms of water loss besides through the spiracular openings.In the system modeled here, the mean surface area of the spiracle, trachea and respiratory cells affect both oxygen intake, carbon dioxide production, and water loss. Worth noting is the fact that these affect carbon dioxide production only indirectly. Hence, the carrying capacity may be defined as a function of the spiracle size, tracheal length, and the respiring cells size. This would mean that the trachea can sustain more oxygen when airflow (air-supply) is higher. In addition, for a given density, oxygen intake rate is higher when the size of both the cells and spiracle is larger. We then define the carrying capacityK as
(5)
K
=
i
Ψ
1
+
m
Ψ
,
where Ψ is mean surface area of the spiracle, trachea, and cells and i here governs how K increases with Ψ during the different phases of DGC. To allow for other limiting factors such as the ambient temperature, humidity, spiracle location, and body size, we do not let K grow indefinitely with Ψ. Instead, the size of m decides how K is limited due to the constraints of these factors.Ratio-dependence, as opposed to dependence on total size ofΨ, is a mechanistic way of modeling oxygen intake in that a unit of oxygen is affected by the amount of oxygen molecules surrounding it through sharing the navigation space. This model for oxygen is captured in our system. To simplify the model, we allow the respiring cells to take up only new oxygen as only a very small volume of O
2 remains in the trachea from the previous inhalation.


## 4. Ultimate Bounds for the Gas Volumes

The following comparison argument [32] is employed in the proofs associated with the boundedness of the gas volumes. Consider the respective solutions ν
1, ν
2, and ν
3 of the initial value problem
(6)
ν
j
′
=
Φ
j
(
t
,
ν
j
)
in
(
0
,
T
)
,
ν
j
(
0
)
=
ν
j
,
0
,
where j
=
1,2
,
3. Φ
1, Φ
2, and Φ
3 are continuous in [
0
,
T
]
×
R. We have the following.Lemma 1.
The comparison argument: assume that∂
Φ
1
/
∂
t, ∂
Φ
2
/
∂
t, and ∂
Φ
3
/
∂
t are continuous in [
0
,
T
]
×
R. Comparing any two, if Φ
1
(
t
,
ν
)
≤
Φ
2
(
t
,
ν
) in (
0
,
T
]
×
R and ν
1,0
≤
ν
2,0, then the respective solutions ν
1 and ν
2 of (6) satisfy ν
1
(
t
)
≤
ν
2
(
t
) on [
0
,
T
].Note that the comparison ofΦ
2
(
t
,
ν
)
≤
Φ
3
(
t
,
ν
) follows as above.Theorem 2.
The amount of oxygen,O
2, in our model system satisfies
(7)
[
(
1
O
2
(
0
)
+
u
a
K
)
e
a
t

-

u
a
K
]
-
1

≤
O
2
(
t
)

≤
[
(
1
O
2
(
0
)

-

1
K
)
e
-
u
t

+

1
K
]
-
1
.Proof.
From the nonnegativity of the gas quantities, we have(8)
d
O
2
d
t
≥
-
u
O
2
2
K
-
a
O
2
.
By the comparison argument in Lemma 1, we know that
(9)
O
2
(
t
)
≥
[
(
1
O
2
(
0
)
+
u
a
K
)
e
a
t
-
u
a
K
]
-
1
,
which is the solution of the initial value problem
(10)
ν
′
=
-
u
ν
2
K
-
b
ν
in
(
0
,
∞
)
,
ν
(
0
)
=
O
2
(
0
)
.
On the other hand, for the upper bound of O
2
(
t
) we have
(11)
d
O
2
d
t
≤
u
O
2
-
u
O
2
2
K
.
Once again from the comparison argument, we have
(12)
O
2
(
t
)
≤
[
(
1
O
2
(
0
)
-
1
K
)
e
-
u
t
+
1
K
]
-
1
,
which is the solution of the initial value problem
(13)
ν
′
=
u
ν
-
u
ν
2
K
in
(
0
,
∞
)
,
ν
(
0
)
=
O
2
(
0
)
.From the results of Theorem2, we can then determine the ultimate bounds for the amount of oxygen in the trachea, O
2
(
t
); namely,
(14)
0
≤
lim
⁡
t
→
∞
inf
⁡
(
O
2
(
t
)
)
≤
lim
⁡
t
→
∞
sup
⁡
(
O
2
(
t
)
)
≤
K
.Next we look at the ultimate bounds for the opportunistic gases; in our model study we can assume bothC
O
2 and H
2
O.Theorem 3.
WhenN
j
≤
K  (
j
=
1,2
), the gas volumes C
O
2 and H
2
O satisfy
(15)
(
C
O
2
(
0
)
+
H
2
O
(
0
)
)
e
-
μ
t


≤
C
O
2
(
t
)
+
H
2
O
(
t
)


≤

(
C
O
2
(
0
)

+

H
2
O
(
0
)
)

e
(
a
K
/
(
a
+
K
)
-
μ
)
t
.Proof.
DefineP
=
C
O
2
+
H
2
O, then
(16)
d
P
d
t
=
a
O
2
C
O
2
a
+
O
2
-
μ
C
O
2
+
a
C
O
2
H
2
O
a
+
C
O
2
-
μ
H
2
O
≥
-
μ
P
,
where μ
=
s
×
d
ma
x
C
O
2. By the comparison argument in Lemma 1, we have
(17)
C
O
2
(
t
)
+
H
2
O
(
t
)
=
P
(
t
)
≥
(
C
O
2
(
0
)
+
H
2
O
(
0
)
)
e
-
μ
t
.
When N
j
(
0
)
≤
K then N
j
(
t
)
≤
K for all t. It can also be seen from the model that
(18)
d
P
d
t
=
a
O
2
C
O
2
a
+
O
2
-
μ
C
O
2
+
a
C
O
2
H
2
O
a
+
C
O
2
-
μ
H
2
O
≤
a
K
[
C
O
2
a
+
K
+
H
2
O
a
+
K
]

-

μ
(
C
O
2
+
H
2
O
)

=

a
K
P
a
+
K
-
μ
P
≤
(
a
K
a
+
K
-
μ
)
P
.
Again by the comparison argument, we have
(19)
C
O
2
(
t
)
+
H
2
O
(
t
)
=
P
(
t
)
≤
P
(
0
)
e
(
a
K
/
(
a
+
K
)
-
μ
)
t
=
(
C
O
2
(
0
)
+
H
2
O
(
0
)
)
e
(
a
K
/
(
a
+
K
)
-
μ
)
t
.The following corollary gives a condition for the extinction of the opportunistic gases and will be used later for the stability of a semitrivial equilibrium solution.Corollary 4.
Ifa
K
/
(
a
+
K
)
<
μ and N
j
(
0
)
≤
K (j
=
1,2), then
(20)
lim
⁡
t
→
∞
(
C
O
2
(
t
)
+
H
2
O
(
t
)
)
=
0
.This result shows that the persistence or extinction of the dependent gases is directly affected by the carrying capacity and initial volume of the tracheal gas volume, as well as other ecological parameters such as the ingress rate, respiration rate, and saturation value of the functional response.


## 5. Equilibrium Solutions and Stability

Our model system has five equilibrium solutions which are summarized in Table1.Table 1
Five equilibrium solutions to our system.
Equilibrium
O2
CO2
H2O
E
1
N
1
(
N
d
+
N
1
W
1
)
/
N
f
0
E
2
N
1
(
N
d
+
N
1
W
1
)
/
N
f
(
N
e
+
N
1
W
2
)
/
N
g
E
3
N
2
(
N
d
+
N
2
W
1
)
/
N
f
0
E
4
N
2
(
N
d
+
N
2
W
1
)
/
N
f
(
N
e
+
N
2
W
2
)
/
N
g
E
5
K
0
0The entries in Table1 are defined as follows. For simplicity, we define p
=
d
maxC
O
2 and q
=
δ
C
O
2, as follows:
(21)
N
1
=
(
N
-
K
N
a
N
b
)
2
N
c
,
N
2
=
(
N
+
K
N
a
N
b
)
2
N
c
,
N
=
-
a
2
c
g
K
q
u
+
a
2
g
2
K
p
q
u
+
a
2
(
-
g
)
K
q
2
u
+
a
b
c
g
K
q
u
2
-
a
b
g
2
K
p
q
u
2
-
2
a
b
g
K
p
q
u
+
a
b
g
K
q
2
u
2
+
2
a
c
g
K
p
u
-
2
a
c
K
q
u
-
K
u
(
a
c
)
2
-
a
g
2
K
p
2
u
+
2
a
g
K
p
q
u
-
a
K
q
2
u
+
2
b
2
g
K
p
q
u
2
-
b
c
K
p
u
+
b
g
K
p
2
u
-
2
b
g
K
p
q
u
2
-
b
K
p
q
u
+
c
K
p
u
-
g
K
p
2
u
+
K
p
q
u
,
N
a
=
c
2
-
2
c
g
p
+
2
c
q
+
g
2
p
2
+
2
g
p
q
+
q
2
,
N
b
=
u
(
a
2
g
q
-
a
b
g
q
u
+
a
c
-
a
g
p
+
a
q
+
b
p
-
p
)
,
N
c
=
-
a
2
g
2
q
2
u
2
+
a
b
c
g
q
u
2
-
a
b
g
2
p
q
u
2
+
a
b
g
q
2
u
2
-
a
c
g
q
u
2
+
a
g
2
p
q
u
2
-
a
g
q
2
u
2
+
b
2
g
p
q
u
2
-
2
b
g
p
q
u
2
+
g
p
q
u
2
,
N
d
=
a
3
(
-
g
)
K
q
+
2
a
2
b
g
K
q
u
-
a
2
c
K
+
a
2
g
K
p
-
a
2
K
q
-
a
b
2
g
K
q
u
2
+
a
b
c
K
u
-
a
b
g
K
p
u
+
a
b
K
q
u
+
a
K
p
,
N
e
=
a
3
(
-
e
)
g
K
q
+
2
a
2
b
e
g
K
q
u
-
a
2
b
f
g
q
h
u
+
a
2
b
g
q
s
h
u
-
a
2
c
e
K
+
a
2
e
g
K
p
-
a
2
e
K
q
-
a
b
2
e
g
K
q
u
2
+
a
b
2
f
g
q
h
u
2
-
a
b
2
g
q
s
h
u
2
+
a
b
c
e
K
u
-
a
b
c
f
h
u
+
a
b
c
s
h
u
-
a
b
e
g
K
p
u
+
a
b
e
K
q
u
+
a
b
f
g
p
h
u
-
a
b
f
q
h
u
-
a
b
g
p
s
h
u
+
a
b
q
s
h
u
+
a
e
K
p
-
b
2
f
p
h
u
+
b
2
p
s
h
u
+
b
f
p
h
u
-
b
p
s
h
u
,
N
f
=
a
2
b
g
q
u
-
a
b
2
g
q
u
2
+
a
b
c
u
-
a
b
g
p
u
+
a
b
q
u
+
b
2
p
u
-
b
p
u
,
N
g
=
a
2
b
f
g
q
u
-
a
2
b
g
q
s
u
-
a
b
2
f
g
q
u
2
+
a
b
2
g
q
s
u
2
+
a
b
c
f
u
-
a
b
c
s
u
-
a
b
f
g
p
u
+
a
b
f
q
u
+
a
b
g
p
s
u
-
a
b
q
s
u
+
b
2
f
p
u
-
b
2
p
s
u
-
b
f
p
u
+
b
p
s
u
.
Finally,
(22)
W
1
=
a
2
(
-
b
)
g
q
u
-
a
2
g
p
q
u
+
a
b
2
g
q
u
2
-
a
c
u
+
a
g
p
u
-
a
q
u
-
b
p
u
+
p
u
,
W
2
=
a
2
(
-
b
)
e
g
q
u
-
a
2
e
g
q
u
+
a
b
2
e
g
q
u
2
+
a
c
e
u
+
a
e
g
p
u
-
a
e
q
u
-
b
e
p
u
+
e
p
u
.
In Table 2, we summarize the conditions for a nonnegative equilibrium solution.Table 2
Conditions for a nonnegative equilibrium.
Equilibrium
Condition
E
1 & E
2
(1)N
>
K
N
a
N
b & N
c
>
0
(2)N
d
>
-
N
1
W
1 & N
f
>
0
(3)N
e
>
-
N
1
W
2 & N
g
>
0
E
3 & E
4
(1)N
>
-
K
N
a
N
b & N
c
>
0
(2)N
d
>
-
N
2
W
1 & N
f
>
0
(3)N
e
>
-
N
2
W
2 & N
g
>
0
E
5
(1)K
>
0We discuss the stability of the equilibrium solutionsE
1
–
E
5. Through the method of linearization, it is known that an equilibrium solution is asymptotically stable if all the eigenvalues of the Jacobian matrix have negative real parts. For ease of computation, define x
=
O
2, y
=
C
O
2, and z
=
H
2
O in this section only. The Jacobian matrix of our system is
(23)
J
=
(
η
1
η
2
η
3
φ
1
φ
2
φ
3
ξ
1
ξ
2
ξ
3
)
,
where,
(24)
η
1
=
-
a
y
b
y
+
x
+
a
(
x
y
)
(
b
y
+
x
)
2
-
a
x
K
+
u
(
1
-
x
+
y
K
)
,
η
2
=
-
a
b
x
K
-
a
x
b
y
+
x
+
a
(
x
y
)
(
b
y
+
x
)
2
,
η
3
=
0
,
φ
1
=
a
c
y
x
(
a
g
+
1
)
+
b
y
-
a
c
(
(
a
g
+
1
)
(
x
y
)
)
(
x
(
a
g
+
1
)
+
b
y
)
2
+
a
q
y
b
y
+
x
-
a
q
(
x
y
)
(
b
y
+
x
)
2
,
φ
2
=
a
c
x
x
(
a
g
+
1
)
+
b
y
-
a
b
c
(
x
y
)
(
x
(
a
g
+
1
)
+
b
y
)
2
+
a
q
x
b
y
+
x
-
a
b
q
(
x
y
)
(
b
y
+
x
)
2
-
p
,
φ
3
=
0
,
ξ
1
=
0
,
ξ
2
=
-
e
z
h
+
z
,
ξ
3
=
-
e
y
h
+
z
+
e
(
y
z
)
(
h
+
z
)
2
+
(
f
-
s
)
.The stability analysis here is carried out on the water vapor differential equation parameters:f, s, e, and h. This is motivated by the fact that insect survival solely depends on its ability to minimize water loss to resist dessication. Hence, the model system stability depends heavily on these parameters.Theorem 5.
The equilibrium solutionE
5 is unstable.Proof.
By substituting the equilibriumE
5 into the Jacobian matrix, J, the associated eigenvalue problem is
(25)
det
⁡
|
J
-
λ
I
|
=
0
.
This is equivalent to
(26)
(
-
a
-
λ
)
(
f
-
λ
-
s
)
(
a
c
a
g
+
1
+
a
q
-
λ
-
p
)
=
0
.
Since f
>
s, at least one eigenvalue is positive. As long as the insect is alive, the insect water ingress rate should be greater than the cuticular water loss. This indicates that E
5 is unstable.
E
5 is asymptotically stable when f
<
s. For O
2
(
0
)
≤
K; Corollary 4 indicates that
(27)
lim
⁡
t
→
∞
(
C
O
2
(
t
)
+
H
2
O
(
t
)
)
=
0
.
Note that O
2
(
t
) converges (t
→
∞) to K which is the global attracting equilibrium in the logistic equation
(28)
d
O
2
d
t
=
u
O
2
(
1
-
O
2
K
)
.The next theorems state new results concerning the stability of equilibrium solutionsE
1
–
E
4. Due to the complexity of the eigenvalue problems associated with these equilibrium solutions, we only obtained conditions of instability.Theorem 6.
AssumeK
>
N
>
0, where N is as defined earlier. The equilibrium solutions E
1 and E
3 are unstable if
(29)
f
-
s
e
>
N
d
+
N
j
W
1
h
N
f
,
where j
=
1,2 and h
>
0.Proof.
By substitutingE
1 into the Jacobian matrix, J, and then simplifying the resulting eigenvalue problem, one of the eigenvalues is given by
(30)
λ
=
f
-
s
-
e
N
d
+
N
1
W
1
h
N
f
.
Now, λ
>
0 implies
(31)
f
-
s
-
e
N
d
+
N
1
W
1
h
N
f
>
0
.
Hence
(32)
f
-
s
e
>
N
d
+
N
1
W
1
h
N
f
.SubstitutingE
3 into J follows the same analysis.Theorem 7.
AssumeK
>
0 and e
>
0, the equilibrium solutions E
2 and E
4 are unstable if either of these conditions hold:
(33)
(
i
)
N
e
+
N
1
W
1
>
N
g
(
h
+
G
G
j
)
,
(
i
i
)
N
d
+
N
j
W
1
N
f
(
h
+
G
G
j
)
[
N
e
+
N
j
W
2
N
g
(
h
+
G
G
j
)
-
1
]
>
s
-
f
e
,
where G
G
j
=
(
N
e
+
N
j
W
2
)
/
N
g, j
=
1,2 and h
>
0.Proof.
By substitutingE
2 into the Jacobian matrix, J, the third eigenvalue is given by
(34)
λ
=
e
(
N
1
W
1
+
N
d
)
(
N
1
W
2
+
N
e
)
N
f
N
g
(
(
N
1
W
2
+
N
e
)
/
N
g
+
h
)
2
-
e
(
N
1
W
1
+
N
d
)
N
f
(
(
N
1
W
2
+
N
e
)
/
N
g
+
h
)
+
f
-
s
.
Now, λ
>
0 implies
(35)
f
-
s
+
e
N
1
W
1
+
N
d
N
f
(
h
+
G
G
1
)
[
N
e
+
N
1
W
2
N
g
(
h
+
G
G
1
)
-
1
]
>
0
,
where G
G
1
=
(
N
e
+
N
1
W
2
)
/
N
g.
The above inequality gives rise to two cases of instability. The first case(
i
) which is independent of e occurs when N
e
+
N
1
W
2
>
N
g
(
h
+
G
G
1
). In this case, the inequality (35) will hold.
The second case(
ii
) occurs when N
g
(
h
+
G
G
1
)
>
N
e
+
N
1
W
2. The inequality (35) holds when
(36)
N
1
W
1
+
N
d
N
f
(
h
+
G
G
1
)
[
N
e
+
N
1
W
2
N
g
(
h
+
G
G
1
)
-
1
]
>
s
-
f
e
.By substitutingE
4 into the Jacobian matrix and then simplifying the resulting matrix as above, similar results are found for the stability.


## 6. Numerical Simulation

Numerical simulations for the dynamics and pattern of the gas volumes are discussed here. The dynamics of the model system change as we vary the parameter values. This behavior may eventually settle on a steady-state or may exhibit some unstable behavior. In each of the figures, the biological parametersK, f, s, e, and h are chosen according to the conditions of stability or instability of the equilibrium given in the theorems in Sections 4 and 5. A stable steady-state in the plots of the dynamic system suggests that the system trajectories tend to an equilibrium state as t
→
∞. That is, the system is characterized by a fixed loop or an inward spiral. Instability is characterized by unstable behavior, where the trajectories drift in an unbounded fashion and the system never forms an orbit around that point [33]. The initial gas volumes were taken from the sample data received earlier: O
2
(
0
)
=
0.25 mL, C
O
2
(
0
)
=
0.75 mL, and H
2
O
(
0
)
=
1.75 mL.Figures2, 3, and 4 show the time series diagrams of the system of differential equations for different values of the system parameters. The simulations were carried out to determine the ranges in the parameter spaces which support the dynamic behavior of this system. These ranges were then compared with the actual parameter values from the data received. For the model to be applicable, the ranges of the parameter values must be biologically feasible and must contain or enclose the parameter values from the actual data. The system stability was studied by varying one of the critical parameter values while fixing the others.Figure 2
This figure shows the stability of equilibriumE
5 when f
<
s, which is given in Theorem 5. We choose f
=
0.5 and s
=
1. As can be seen from the graph, the oxygen volume goes to equilibrium K
=
1.9 and the carbon dioxide and water-vapor go to equilibrium 0.Figure 3
This figure illustrates the case when equilibriaE
1 and E
3 are both unstable. The condition for instability is (
f
-
s
)
/
e
>
(
N
d
+
N
j
W
1
)
/
h
N
f. For this case we choose K
=
5, u
=
0.7, f
=
0.5, s
=
0.03974, e
=
0.8, and h
=
2.7.Figure 4
This figure shows the stability of bothE
2 and E
4. We choose the parameters a
=
0.6, K
=
1.9, u
=
1, f
=
0.5, s
=
1, e
=
0.5, and h
=
1.Figure2 shows the stability of E
5 when f
<
s, characterized by a quick convergence of C
O
2 and H
2
O towards their equilibrium values. However, this condition is not feasible biologically while the insect is alive. The model stability is achieved when the water ingress rate into the trachea is less than the water loss due to the cyclic openings. This would mean that the insect is losing water and thus the insect eventually dies from dessication. The time series shows that O
2 goes to equilibrium K
=
1.9, the tracheal carrying capacity. C
O
2 and H
2
O go to equilibrium zero. In Figures 2 and 4, for prolonged spiracle opening, the C
O
2 and H
2
O graphs decay to the long-term steady-state values as outlined in Table 1. Values of C
O
2 and H
2
O below zero would mean that the insect is not breathing or may be dead from dessication, hence the condition f
>
s is of interest to biologists.Cuticular water-loss plays a major role in determining the stability of bothE
1 and E
3. Figure 3 illustrates the case when equilibria of E
1 and E
3 are both unstable. The time series is characterized by an outward spiral. For instance, the O
2 volume at the first peak when t
=
9.3 is approximately 3.96 mL. At the second peak at t
=
22.5, O
2 volume is approximately 4.46 mL. For even larger values of t, the O
2 values reach even higher volumes, while both C
O
2 and H
2
O oscillate in a stable manner with no convergence. Figure 4 shows the stability of both E
2 and E
4 for specific parameter values K, u, f, s, e, and h. O
2 converges to the equilibrium N
1 which is slightly below the tracheal carrying capacity, K. The C
O
2 and H
2
O graphs oscillate in a stable manner initially and finally converge to (
N
d
+
N
1
W
1
)
/
N
f and (
N
e
+
N
1
W
2
)
/
N
g, respectively.


## 7. Discussion

Our main objective was to obtain a mathematical model that depicts how the respiratory gases manage to move in opposite directions simultaneously. We developed the model described by equations in Section3 and studied its existence for different ranges of system parameters. Additionally in Section 6 we investigated the variable parameter ranges to control unstable dynamics, if any. The simulation results suggest that the model system may show different interesting dynamical behavior starting from stable limit cycle to unstable cycles.In Section4, we found that the O
2 and C
O
2 gases in the trachea were bounded below by zero and above by twice the tracheal carrying capacity. Hence, as long as the insect is alive, tracheal O
2 volume should be positive. Likewise, prolonged O
2 uptake means that respiration is prolonged and thus there is a positive C
O
2 volume in the trachea. The combined tracheal volume of O
2 and C
O
2 cannot exceed 2
K.The stability analysis was undertaken on the water vapor differential equation parameters:f, s, e, and h in Section 6. From the simulation results, these parameters are responsible for determining the stability of our system and thus the insect survival. We found that f, s, and h were bounded over a large parameter regime compared to the other system parameters such as a and e. Our system supports stable limit cycles in a reasonably fair parameter range for f, s, and h. Unpredictable behavior was observed outside these studied parameter values in narrow ranges. For instance, the cumulative effect of increased e results in an inward shift in the H
2
O graph as shown in Figure 5. With increased cuticular water loss rate, the dynamics of H
2
O quickly converge to zero, beyond which we would expect dessication.Figure 5
H
2
O profile in the beetle for varying values of e. Here, all other parameter values are fixed.When the oxygen ingress rate in the trachea is high, the volume of oxygen in the insect respiratory system rises and so does the elimination of carbon dioxide. An increase in the elimination of carbon dioxide results in large volumes of water loss into the atmosphere. A further increase in this oxygen ingress results in further elimination ofC
O
2 and thus unpredictable dynamics will be observed. Figure 6 depicts this observation. For instance, the O
2 and H
2
O graphs oscillate in an unstable manner. As K increases, the spiracle remains open for a longer time to increase O
2 uptake. The longer the spiracular opening is the more H
2
O is lost to the atmosphere. Prolonged H
2
O loss results in unpredictable dynamics.Figure 6
Phase diagram with unstable trajectories. HereK
=
9 and u
=
0.7.The dung beetle receives oxygen from the atmosphere when it opens its spiracles. These high oxygen volumes increase the values ofK and u. The carrying capacity (the insect tracheal system maximum load) limits the amount of oxygen that the insect can take in at any given time. Thus when the values of u and K are high, and keeping other parameters constant, the model system is expected to exhibit some unstable dynamics. Since water is lost as a consequence during the elimination of C
O
2, we expect that an increase in a results in large quantities of water loss. Dessication results due to large quantities of water being lost; hence the cumulative effect of increased a is to force the dynamics of the system to settle on an unstable state.Insects adapt and control for this unstable behavior by performing DGC. During DGC,O
2 uptake and C
O
2 release follow a cyclical pattern defined by periods of little to no release of C
O
2. Thus they have a way of keeping K, u, and a low most of the time. The value of a (respiration rate) is bounded on a narrow range where the system exhibits stable behavior. Search for food and flight from danger demand an increase in the value for a.The competing species model is a useful starting point to illustrate that the respiratory gases compete for space in both the spiracle opening and the tracheal system. The model developed here reveals that the change inC
O
2 elimination is directly related to the level of oxygen uptake into the trachea, while the change in water levels in the trachea is directly proportional to C
O
2 emission. The equilibria are stable over fairly large parameter regimes and exhibit cycles in which prolonged O
2 uptake perpetually rises and falls in a systemic manner.



---

*Source: 157573-2014-09-03.xml*