Gas Exchange Models for a Flexible Insect Tracheal System

Abstract

In this paper two models for movement of respiratory gases in the insect trachea are presented. One model considers the tracheal system as a single flexible compartment while the other model considers the trachea as a single flexible compartment with gas exchange. This work represents an extension of Ben-Tal’s work on compartmental gas exchange in human lungs and is applied to the insect tracheal system. The purpose of the work is to study nonlinear phenomena seen in the insect respiratory system. It is assumed that the flow inside the trachea is laminar, and that the air inside the chamber behaves as an ideal gas. Further, with the isothermal assumption, the expressions for the tracheal partial pressures of oxygen and carbon dioxide, rate of volume change, and the rates of change of oxygen concentration and carbon dioxide concentration are derived. The effects of some flow parameters such as diffusion capacities, reaction rates and air concentrations on net flow are studied. Numerical simulations of the tracheal flow characteristics are performed. The models developed provide a mathematical framework to further investigate gas exchange in insects.

Appendices

Appendix 1: Diffusion Equation Modeled by Krogh

The rate of movement of oxygen ($d O_2 /d t$) across muscle tissues under both steady-state and non-steady state conditions conforms to the one-dimensional diffusion equation for gases as modified by Krogh (1919, 1959),

$$\begin{aligned} \frac{d O_2}{d t} = - \frac{\kappa _{O_2} A p_{O_2} }{L}, \end{aligned}$$

(4.1)

where $\kappa _{O_2}$ is the diffusion constant for $\text{O}_2$ which is the product of the diffusion coefficient for $\text{O}_2$ and its solubility. $\kappa _{O_2} = D_{O_2} \times \alpha _{O_2}$: $D_{O_2}$ is the diffusion coefficient for $\text{O}_2$ and $\alpha _{O_2}$ is the solubility coefficient for $\text{O}_2$. A is the surface area through which exchange occurs, $p_{O_2}$ is the $\text{O}_2$ partial pressure gradient across the diffusion path, and L is the length of the diffusion path.

We also note that changing temperature affects both $D_{O_2}$ and $\alpha _{O_2}$. $D_{O_2}$ in water increases by ${\sim}3\,\%$ for each rise in temperature of $1\,^{\circ }{\hbox {C}}$, while solubility decreases with increasing temperature by ${\sim}1.4\,\%\,^{\circ }{\hbox {C}}^{-1}$. Thus, the thermal sensitivity of the diffusion constant ($\kappa _{O_2}$) is $\sim$ $1.6\%\,^{\circ }{\hbox {C}}^{-1}$. $D_{O_2}$ decreases with decreasing temperature due to:

• 1. reduction in the kinematic energy of the system,

• 2. increase in viscosity (cytoplasm viscosity), which is inversely related to cell temperature.

Appendix 2: Absorption Model

The absorption model can be written as a single ordinary differential equation

$$\begin{aligned} \frac{d x}{d t} = g(x , p_o), \end{aligned}$$

(4.2)

where $x = [ \text{MO}_{2} ]$, i.e., x denotes the concentration of $\text{O}_2$-mitochondria combinations. We then define the mitochondria saturation function as

$$\begin{aligned} S = f({x}). \end{aligned}$$

(4.3)

In Eq. 2.33 we thus have that $N_o = 1$ for $\text{O}_2$ uptake by mitochondria. We take the mitochondria saturation function S to be $S(x) = x / T_m$, and furthermore we assume that $d x / dt = g(x,p_o ) = r \sigma p_o (T_m - x) - l x$.

Thus

$$\begin{aligned} \frac{d S}{d t}&= \frac{1}{T_m} \frac{d x}{d t}, \end{aligned}$$

(4.4)

$$\begin{aligned}&= \frac{1}{T_m} (r \sigma p_o (T_m - x) - l x) . \end{aligned}$$

(4.5)

With $N_o = 1$, substitute Eq. 4.4 into Eq. 2.33,

$$\begin{aligned} \frac{d p_o}{d t}&= \frac{D_o}{v_{ct} \sigma } ( p_{to} - p_o ) - \frac{T_m}{\sigma } \left( \frac{1}{T_m} (r \sigma p_o (T_m - x) - l x) \right) , \end{aligned}$$

(4.6)

$$\begin{aligned}&= \frac{D_o}{v_{ct}\sigma } ( p_{to} - p_o ) - \frac{1}{\sigma } (r \sigma p_o (T_m - x) - l x) . \end{aligned}$$

(4.7)

The above equations can then be written in their non-dimensional form. Recall

$$\begin{aligned} \frac{d p_o}{d t}&= \frac{D_o }{v_{ct} \sigma } ( p_{to} - p_{o} ) - \frac{N_o T_m}{\sigma } \frac{d S}{d t}, \end{aligned}$$

(4.8)

$$\begin{aligned} \frac{d x}{d t}&= g(x,p_o), \end{aligned}$$

(4.9)

$$\begin{aligned} S&= S(x). \end{aligned}$$

(4.10)

Assume that the tracheal partial pressure of oxygen, $p_{to}$, is constant. Divide Eq. 4.8 through by $p_{to} \ne 0$

$$\begin{aligned} \frac{1}{p_{to}} \frac{d p_o}{d t} = \frac{D_o }{v_{ct} \sigma } \left( 1 - \frac{p_o}{p_{to}} \right) - \frac{N_o T_m}{\sigma p_{to}} \frac{d S}{d t}. \end{aligned}$$

(4.11)

Let $p_r = p_o / p_{to} \rightarrow \dfrac{d p_r }{d t} = \dfrac{1}{p_{to}} \dfrac{d p_o}{d t}$. Therefore

$$\begin{aligned} \frac{d p_r}{d t}&= \frac{D_o }{v_{ct} \sigma } (1 - p_r ) - \frac{N_o T_m}{\sigma p_{to}} \frac{d S}{d t}, \end{aligned}$$

(4.12)

$$\begin{aligned} \frac{v_{ct} \sigma }{D_o } \frac{d p_r}{d t}&= 1 - p_r - \frac{N_o T_m}{p_{to}} \frac{v_{ct}}{D_o } \frac{d S}{d t} \end{aligned}$$

(4.13)

Let $\tau = (D_o / v_{ct} \sigma ) t$ such that $d\tau = (D_o / v_{ct} \sigma ) dt$. Thus

$$\begin{aligned} \frac{d p_r}{d \tau } = 1 - p_r - \frac{N_o T_m}{\sigma p_{to}} \frac{d S}{d \tau }. \end{aligned}$$

(4.14)

Let $\alpha = N_o T_m / \sigma p_{to}$. Hence

$$\begin{aligned} \frac{d p_r}{d \tau } = 1 - p_r - \alpha \frac{d S}{d \tau }. \end{aligned}$$

(4.15)

Recall $d x /dt = g(x,p_o)$, $p_o = p_{to} p_r$ and $d\tau = (D_o / v_{ct} \sigma ) dt,$ thus

$$\begin{aligned} \frac{d x}{d t}&= (r \sigma p_o (T_m - x) - l x) , \end{aligned}$$

(4.16)

$$\begin{aligned} \frac{D_o}{v_{ct} \sigma }\frac{d x}{d \tau }&= (r \sigma p_{to} p_r (T_m - x) - l x) , \end{aligned}$$

(4.17)

$$\begin{aligned} \frac{d x}{d \tau }&= \frac{v_{ct} \sigma }{D_o} (r \sigma p_{to} p_r (T_m - x) - l x), \end{aligned}$$

(4.18)

$$\begin{aligned} \frac{1}{T_m} \frac{d x}{d \tau }&= \frac{v_{ct} \sigma }{D_o} \left( r \sigma p_{to} p_r \left( 1 - \frac{x}{T_m}\right) - l \frac{x}{T_m}\right) . \end{aligned}$$

(4.19)

Let $\tilde{x} = x / T_m$. Thus

$$\begin{aligned} \frac{d \tilde{x}}{d \tau }&= \frac{v_{ct} \sigma }{D_o} \left( r \sigma p_{to} p_r \left( 1 - \tilde{x}\right) - l \tilde{x}\right) , \end{aligned}$$

(4.20)

$$\begin{aligned}&= \frac{r \sigma ^{2} v_{ct} p_{to}}{D_o} \left( p_r \left( 1 - \tilde{x}\right) - \frac{l}{r \sigma p_{to}} \tilde{x}\right) . \end{aligned}$$

(4.21)

Let $\beta = r \sigma ^{2} v_{ct} p_{to} / D_o$ and $\gamma = l / r \sigma p_{to}.$

$$\begin{aligned} \therefore ~~~ \frac{d \tilde{x}}{d \tau } = \beta \left( p_r \left( 1 - \tilde{x}\right) - \gamma \tilde{x}\right) . \end{aligned}$$

(4.22)

Appendix 3: Large “Rate Factor ”($\beta \gg 1$)

To understand what happens when $\beta \gg 1$, we consider two cases: when $\alpha = 0$ and $\alpha \ne 0$. The system of Eqs. 2.40−2.42 can be analyzed in the plane ($p_r , \tilde{x}$). Solving the system at equilibrium (i.e. $d p_r / d\tau = 0$ and $d \tilde{x} / d\tau = 0$) we have $\left( p_r , \tilde{x} \right)$ $=$ $\left( 1, \dfrac{1}{ 1 + \gamma }\right)$. At the equilibrium, the Jacobian matrix is given by

$$\begin{aligned} J = \left( \begin{array}{cc} -1 - \alpha \beta \left( \frac{\gamma }{1 + \gamma } \right) &{} \alpha \beta (1 + \gamma ) \\ \beta \left( \frac{\gamma }{1 + \gamma } \right) &{} - \beta (1 + \gamma ) \end{array} \right) . \end{aligned}$$

The eigenvalues are

$$\begin{aligned} \lambda _1&= -\frac{1}{2} \left( \mathcal {P} + 1 + \frac{\alpha \beta \gamma }{(1+\gamma )} \right) + \frac{1}{2} \left( (\mathcal {P} - 1)^{2} + \frac{\alpha \beta \gamma \mathcal {Q}}{(1 + \gamma )} \right) ^{1/2} , \end{aligned}$$

(4.23)

$$\begin{aligned} \lambda _2&= -\frac{1}{2} \left( \mathcal {P} + 1 + \frac{\alpha \beta \gamma }{(1+\gamma )} \right) - \frac{1}{2} \left( (\mathcal {P} - 1)^{2} + \frac{\alpha \beta \gamma \mathcal {Q}}{(1 + \gamma )} \right) ^{1/2} , \end{aligned}$$

(4.24)

where $\mathcal {P} =\beta (1 + \gamma )$ and $\mathcal {Q} = 2\beta (1 + \gamma ) + 2 + (\alpha \beta \gamma ) / (1 + \gamma )$.

Case 1: When $\alpha = 0$, $\lambda _1 = -1$ and $\lambda _2 = -\beta (1 + \gamma )$. The corresponding eigenvectors are

$$\begin{aligned} \mathbf {v_1}&= \left[ 1 , \frac{-\beta \gamma }{(1 + \gamma ) - \beta (1 + \gamma )^{2}} \right] ^{T} , ~~\text{ and } \end{aligned}$$

(4.25)

$$\begin{aligned} \mathbf {v_2}&= [0 , 1 ]^{T}. \end{aligned}$$

(4.26)

For $\beta \gg 1$, $|\lambda _2 | \gg |\lambda _1 |$ and $\mathbf{{v_1}} \approx [1 , \gamma / (1 + \gamma )^{2}]$. That is, at equilibrium, $\mathbf{{v_1}}$ is tangent to the saturation curve $\tilde{S} = p_r / (p_r + \gamma )$ which is one of the nullclines in the system.

Case 2: When $\alpha \ne 0$, trajectories in time get trapped between the nullclines

$$\begin{aligned} \tilde{x}&= \frac{p_r}{p_r + \gamma } ~~ (\text{ at } \text{ which } ~~ d \tilde{x} / d \tau = 0 ), ~~ \text{ and } \end{aligned}$$

(4.27)

$$\begin{aligned} \tilde{x}&= \frac{p_r (1 + \alpha \beta ) - 1}{\alpha \beta (p_r + \gamma )} ~~ ( \text{ at } \text{ which } ~~ d p_r / d \tau = 0 ). \end{aligned}$$

(4.28)

For $\beta \gg 1$, the two nullclines are close to each other (i.e. $\forall$ $\alpha \beta \gg 1$, we have $(p_r (1 + \alpha \beta ) - 1) / (\alpha \beta (p_r + \gamma )) \approx p_r / (p_r + \gamma )$). This means that for large $\beta,$ trajectories will be near the saturation curve. When $\alpha \ne 0$, $|\lambda _1 | < 1$ and $|\lambda _2 | > \beta (1 + \gamma );$ therefore the transient towards the saturation curve is faster than for the case when $\alpha = 0$. As $\alpha \rightarrow \infty$, in the limit, $\lambda _1 \rightarrow 0$ since J becomes singular. Hence, $\lambda _1$ increases monotonically and therefore is never positive (Ben-Tal 2004).

Appendix 4: Model for the Chemical Reactions

We consider the reaction (2.56) given by

$$\begin{aligned} \text{CO}_2 + \text{H}_2 \text{O} \rightleftharpoons _{\delta l_2}^{\delta r_2} \text{H}^{+} + \text{HCO}_3^{-}. \end{aligned}$$

(4.29)

The chemical reaction considered above is not coupled with either diffusion and ventilation, thereby mimicking an experiment in-vitro.

The rate equations for each of the reactions is given by

$$\begin{aligned} \bar{r}_1&= \delta r_2 [\text{CO}_2 (t)] [\text{H}_2 \text{O} (t)] , \end{aligned}$$

(4.30)

$$\begin{aligned} \bar{r}_2&= \delta l_2 [\text{HCO}_3^{-} (t)] [\text{H}^{+}] , \end{aligned}$$

(4.31)

where $\text{H}^{+}$ is assumed to be constant.

The governing equations for the rate of change of the concentrations of each chemical species depend on the rate equations and the reactions

$$\begin{aligned}{}[\text{CO}_2 (t)]^{\prime }&= \bar{r}_2 - \bar{r}_1 , \end{aligned}$$

(4.32)

$$\begin{aligned}{} [{\text{H}}_2 {\text{O}}(t)]^{\prime }&= \bar{r}_2 - \bar{r}_1 , \end{aligned}$$

(4.33)

$$\begin{aligned}{} [{\text{HCO}}_3^{-}(t)]^{\prime }&= \bar{r}_1 - \bar{r}_2 . \end{aligned}$$

(4.34)

Appendix 5: Parameters and Variables

This section lists the parameters and variables that appear frequently in this manuscript (unless stated otherwise in the text). Table 1 lists the parameters and their values used in the model simulations. Note that the parameter values are adapted from literature data and in cases where they could not be found, they were estimated. Even though some were assumed, there are also some parameters that were calculated. The diffusion capacities, $D_o$ and $D_c$, depend on the surface area and the thickness of the membrane across which the diffusion takes place. It should be noted that no attempt has been carried out to carefully fit parameters to a specific experiment since the developed models are simple and cannot be expected to match experimental data exactly. Variables are presented in Table 2. Note that t is an independent variable here. For the duration of the DGC phenomena in insects, $\nu$ would be lower than used in the simulations but the results would not change qualitatively. During periods of flutter, the pressure within the tracheal system is very close to atmospheric, although slightly negative. The negative pressure causes an inward mass flow of air. When the spiracles completely close, the pressure in the tracheal system drops precipitously. Thus for $P_s$ and the initial value of $P_T$ we use 760 mmHg which is the atmospheric pressure at sea level.

Table 1 Parameter list

Note that $\delta$ is interpreted as a (dimensionless) quantity which affects the buffering capacity of the hemolymph. $\mathcal {T}$ is given as mean ± standard deviation (SD).

Table 2 Variables initial values list