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

#### 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 and are being born and die. is born when it enters the trachea and starts to die as it is used by the cells. is replaced by . In the same way if there is sufficient prey the predator population expands at the expense of the prey. So increases as 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 enters the trachea, is used up by the cells and is then replaced by . 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 and 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, needs to be replaced and 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 Section 2, 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 (Figure 1). 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.

**(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 (-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 () within the trachea declines and the partial pressure of carbon dioxide () increases. Owing to the high solubility of , the increase in in the air within the trachea is less than the decline in . As a result, total atmospheric pressure in the tracheal lumina declines slightly during the closed phase [18].

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

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

##### 2.3. Open Spiracle Phase (-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]. and 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, this is stored in the tissues and haemolymph (similar to blood in vertebrates). When this storage is saturated the 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 the is being used by the cells and is stored. When there is sufficient build-up of , this increase in concentration causes the spiracles to open. Upon opening rushes out of the tracheal system and 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 , (ii) Carbon dioxide denoted by . (iii) Water vapor denoted by .

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 -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, , is described by where is the carbon dioxide volume, is the oxygen volume, is the amount of taken in by the respiring cell when is in excess, and 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 where is the respiratory quotient.

Now, we assume that the carbon dioxide elimination rate is a function of oxygen intake rate, . Thus, the production and elimination functions are, respectively, where measures how the cells turn the inspired oxygen into carbon dioxide, measures how carbon dioxide production is physiologically constrained by these cells, measures carbon dioxide maximum per capita elimination rate, and 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 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: with , , , , , , , , and .

, , and denote the gas volumes. is the oxygen growth rate at low oxygen densities, is the tracheal carrying capacity of the oxygen gas, is the maximum spiracular water loss in the thorax and elytral case as a result of the insect respiration, is the water ingress rate into the trachea, and is the water density/volume where water loss is half saturated. Lastly, 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 capacity as where is mean surface area of the spiracle, trachea, and cells and here governs how 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 grow indefinitely with . Instead, the size of decides how 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 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 , , and of the initial value problem where . , , and are continuous in . We have the following.

Lemma 1. *The comparison argument: assume that , , and are continuous in . Comparing any two, if in and , then the respective solutions and of (6) satisfy on .*

Note that the comparison of follows as above.

Theorem 2. *The amount of oxygen, , in our model system satisfies
*

*Proof. *From the nonnegativity of the gas quantities, we have
By the comparison argument in Lemma 1, we know that
which is the solution of the initial value problem
On the other hand, for the upper bound of we have
Once again from the comparison argument, we have
which is the solution of the initial value problem

From the results of Theorem 2, we can then determine the ultimate bounds for the amount of oxygen in the trachea, ; namely,

Next we look at the ultimate bounds for the opportunistic gases; in our model study we can assume both and .

Theorem 3. *When , the gas volumes and satisfy
*

*Proof. *Define , then
where . By the comparison argument in Lemma 1, we have
When then for all . It can also be seen from the model that
Again by the comparison argument, we have

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. *If and (), then
*

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 Table 1.

The entries in Table 1 are defined as follows. For simplicity, we define and , as follows: Finally, In Table 2, we summarize the conditions for a nonnegative equilibrium solution.

We discuss the stability of the equilibrium solutions . 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 , , and in this section only. The Jacobian matrix of our system is where,

The stability analysis here is carried out on the water vapor differential equation parameters: , , , and . 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 solution is unstable.*

*Proof. *By substituting the equilibrium into the Jacobian matrix, , the associated eigenvalue problem is
This is equivalent to
Since , 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 is unstable.

is asymptotically stable when . For ; Corollary 4 indicates that
Note that converges () to which is the global attracting equilibrium in the logistic equation

The next theorems state new results concerning the stability of equilibrium solutions . Due to the complexity of the eigenvalue problems associated with these equilibrium solutions, we only obtained conditions of instability.

Theorem 6. *Assume , where is as defined earlier. The equilibrium solutions and are unstable if
**
where and .*

*Proof. *By substituting into the Jacobian matrix, , and then simplifying the resulting eigenvalue problem, one of the eigenvalues is given by
Now, implies
Hence

Substituting into follows the same analysis.

Theorem 7. *Assume and , the equilibrium solutions and are unstable if either of these conditions hold:
**
where , and .*

*Proof. *By substituting into the Jacobian matrix, , the third eigenvalue is given by
Now, implies
where .

The above inequality gives rise to two cases of instability. The first case which is independent of occurs when . In this case, the inequality (35) will hold.

The second case occurs when . The inequality (35) holds when

By substituting 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 parameters , , , , and 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 . 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: mL, mL, and mL.

Figures 2, 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 shows the stability of when , characterized by a quick convergence of and 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 goes to equilibrium , the tracheal carrying capacity. and go to equilibrium zero. In Figures 2 and 4, for prolonged spiracle opening, the and graphs decay to the long-term steady-state values as outlined in Table 1. Values of and below zero would mean that the insect is not breathing or may be dead from dessication, hence the condition is of interest to biologists.

Cuticular water-loss plays a major role in determining the stability of both and . Figure 3 illustrates the case when equilibria of and are both unstable. The time series is characterized by an outward spiral. For instance, the volume at the first peak when is approximately mL. At the second peak at , volume is approximately mL. For even larger values of , the values reach even higher volumes, while both and oscillate in a stable manner with no convergence. Figure 4 shows the stability of both and for specific parameter values , , , , , and . converges to the equilibrium which is slightly below the tracheal carrying capacity, . The and graphs oscillate in a stable manner initially and finally converge to and , 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 Section 3 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 Section 4, we found that the and 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 volume should be positive. Likewise, prolonged uptake means that respiration is prolonged and thus there is a positive volume in the trachea. The combined tracheal volume of and cannot exceed .

The stability analysis was undertaken on the water vapor differential equation parameters: , , , and 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 , , and were bounded over a large parameter regime compared to the other system parameters such as and . Our system supports stable limit cycles in a reasonably fair parameter range for , , and . Unpredictable behavior was observed outside these studied parameter values in narrow ranges. For instance, the cumulative effect of increased results in an inward shift in the graph as shown in Figure 5. With increased cuticular water loss rate, the dynamics of quickly converge to zero, beyond which we would expect dessication.

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 of and thus unpredictable dynamics will be observed. Figure 6 depicts this observation. For instance, the and graphs oscillate in an unstable manner. As increases, the spiracle remains open for a longer time to increase uptake. The longer the spiracular opening is the more is lost to the atmosphere. Prolonged loss results in unpredictable dynamics.

The dung beetle receives oxygen from the atmosphere when it opens its spiracles. These high oxygen volumes increase the values of and . 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 and 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 , we expect that an increase in results in large quantities of water loss. Dessication results due to large quantities of water being lost; hence the cumulative effect of increased 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, uptake and release follow a cyclical pattern defined by periods of little to no release of . Thus they have a way of keeping , , and low most of the time. The value of (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 .

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 in 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 emission. The equilibria are stable over fairly large parameter regimes and exhibit cycles in which prolonged uptake perpetually rises and falls in a systemic manner.

#### Disclosure

This paper is a collaborative effort.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

Simphiwe Simelane thanks the University of the Witwatersrand for a Postgraduate Merit Award. Supervisors, Professors Shirley Abelman and Frances Duncan thank the NRF South Africa for financial support. The reviewers are thanked for their valuable comments and suggestions which have significantly improved the paper.