A joint computational respiratory neural network-biomechanical model for breathing and airway defensive behaviors

Data-driven computational neural network models have been used to study mechanisms for generating the motor patterns for breathing and breathing related behaviors such as coughing. These models have commonly been evaluated in open loop conditions or with feedback of lung volume simply represented as a filtered version of phrenic motor output. Limitations of these approaches preclude assessment of the influence of mechanical properties of the musculoskeletal system and motivated development of a biomechanical model of the respiratory muscles, airway, and lungs using published measures from human subjects. Here we describe the model and some aspects of its behavior when linked to a computational brainstem respiratory network model for breathing and airway defensive behavior composed of discrete “integrate and fire” populations. The network incorporated multiple circuit paths and operations for tuning inspiratory drive suggested by prior work. Results from neuromechanical system simulations included generation of a eupneic-like breathing pattern and the observation that increased respiratory drive and operating volume result in higher peak flow rates during cough, even when the expiratory drive is unchanged, or when the expiratory abdominal pressure is unchanged. Sequential elimination of the model’s sources of inspiratory drive during cough also suggested a role for disinhibitory regulation via tonic expiratory neurons, a result that was subsequently supported by an analysis of in vivo data. Comparisons with antecedent models, discrepancies with experimental results, and some model limitations are noted.

Introduction

The neural mechanisms that regulate and coordinate breathing and respiratory-related behaviors such as coughing are not well understood. This lack of knowledge hampers elucidation of pathophysiological deficits in airway protection and impedes development of new therapeutic approaches for dystussia that occur with neurological disorders (Suárez et al., 2002; McCool, 2006). Computational neural network models for breathing and the neurogenesis of cough inferred from in vivo experiments have iteratively aided prediction and refinement of hypotheses for further in vivo testing (Shannon et al., 1998, 2000; Baekey et al., 2001; Rybak et al., 2008; Poliaček et al., 2011). Such data-driven models, based in part on elements and connectivity inferred from simultaneous extracellular recordings of many brainstem neurons (e.g., Segers et al., 2008; Ott et al., 2012), have largely been evaluated in either open loop conditions or with feedback of lung volume simply represented as a filtered version of the motor output (Lindsey et al., 2012). While useful, these approaches have precluded model-based assessment of the potential influence of mechanical properties of the musculoskeletal system on respiratory motor pattern generation and cough effectiveness (Smith et al., 2012). More generally, neuromechanical models can provide a framework for estimating and predicting the extent to which motor patterns are constrained and influenced by mechanical properties and muscle synergies (Chiel et al., 2009).

A model that relates a respiratory neural output to mechanical outputs has been available for some time (Riddle and Younes, 1981; Younes and Riddle, 1981; Younes et al., 1981) and remains an important element in contemporary models of the respiratory system (Cheng et al., 2010; Cheng and Khoo, 2012). However, the Younes–Riddle model with its single inspiratory neural output and single state variable (lung volume) lacks features essential for model-based assessment of the respective contributions and interactions of neural and biomechanical mechanisms during cough. We have developed a respiratory neural network model with inspiratory (phrenic), expiratory (lumbar), and laryngeal neural outputs and required a mechanical model with corresponding inputs to control the abdominal, diaphragm, and laryngeal muscles. Moreover, it is well known that a given lung volume can be achieved with different configurations of the rib cage and abdomen (Konno and Mead, 1967; Younes and Riddle, 1981), and with separate neural control of the diaphragm and abdominal wall muscles as in our network model, all of these configurations can potentially be achieved. Thus, the first aim of the work reported here was to develop a model of the mechanical respiratory system that includes separate muscle models for the diaphragm, abdominal wall, and larynx, and two state variables to represent the thoracoabdominal configuration.

Our second aim was to link the resulting mechanical subsystem to an enhanced integrate and fire (IF) neural network model of the brainstem network for respiratory motor pattern generation and to assess the integrated system’s behavior with muscle activation parameters for eupneic conditions. A third related objective was to extend the simulations to include evoked coughs in order to evaluate the model’s performance in response to defined perturbations that enhance or reduce inspiratory drive. This latter goal was motivated in part by evidence that changes in the inspiratory or “operating” volume can influence airflow during the expulsive phase of the cough (Smith et al., 2012).

An additional impetus for the third aim came from results of recent model simulations, which suggested that elevated systemic arterial blood pressure – such as may occur during coughing (Sharpey-Schafer, 1953) – attenuates cough inspiratory drive, a result supported by coordinated in vivo experiments (Poliaček et al., 2011). The network model used in the present work builds upon that and other recent prior efforts (Rybak et al., 2008). The network incorporates multiple circuit paths and operations for tuning inspiratory drive that have been inferred from spike train cross-correlation feature sets (Lindsey et al., 1998; Shannon et al., 2000; Segers et al., 2008; Ott et al., 2012). These circuits include parallel channels for modulation of inspiratory phase activity in “tonic” expiratory neurons that inhibit premotor inspiratory bulbospinal neurons and drive.

In the course of sequentially eliminating sources of inspiratory drive for cough in the neuromechanical model, we also noted a contribution of tonic expiratory neuron activity to modulation of inspiratory phase drive during cough. This disinhibitory regulation predicted from the modeling results was subsequently supported by an analysis of in vivo data as described in a companion report (Segers et al., 2012).

Materials and Methods

Neural circuit components were derived from previously described respiratory network models of discrete “IF” populations after MacGregor (1987) and a “hybrid IF burster” population with Hodgkin–Huxley style equations after Breen et al. (2003). These models were developed iteratively with in vivo experiments that both guided model development and tested model predictions, as detailed in Rybak et al. (2008) and Poliaček et al. (2011). The enhanced network model used herein is described further in the Results.

Biomechanical model elements were developed using parameters derived from published work as described in Results. Of particular importance was the work of Grassino et al. (1978), who measured transdiaphragmatic pressure and diaphragm activation while controlling the thoracoabdominal configuration, making it possible to estimate the effect of rib cage motion on the abdominal volume. Our abdominal wall model is based on measurements of the curvature of the abdomen by Song et al. (2006) taken during insufflation for laparoscopic surgery. The rib cage, lung, and diaphragm volumes are derived from the measurements of Cluzel et al. (2000), who measured them from MRI’s. The thoracoabdominal configuration at extreme and resting supine lung volumes are from Konno and Mead (1967).

Models were implemented using a program package written in the C language for the UNIX environment. Simulations were run on 64-bit Intel multiprocessor-based computers under the Linux operating system. The GNU Scientific Library was used to solve the differential equations of the biomechanical model, to find the roots of the implicit model equations, to do the abdominal volume integration, and for a spline approximation of the abdominal volume function.

For each condition of the linked neural network and biomechanical model, four trials were run with different random number seeds for the stochastic network model. A pairwise two-sided t-test with non-pooled SD was used for each variable, and the p-values were adjusted for multiple testing (Holm, 1979). A difference was considered significant if the adjusted p-value was less than 0.05.

Results

The results are presented in two main parts. Section “Mechanical Model Implementation: Respiratory Muscles, Chest Wall, and Lungs” details the biomechanical model. Section “Brainstem Network Model Architecture and System Performance When Linked to the Biomechanical Model” describes the linkage between the biomechanical and neural network models and neuromechanical system behavior during various perturbations of the network.

Mechanical Model Implementation: Respiratory Muscles, Chest Wall, and Lungs

The biomechanical model described below converts respiratory neural outputs in the form of spike trains representing lumbar, phrenic, expiratory laryngeal, and inspiratory laryngeal motor neuron activity generated by a stochastic model of the brainstem respiratory network deterministically into mechanical outputs such as lung volume, tracheal flow, and alveolar pressure for a supine male human. Lung volume is fed back to the network model to simulate pulmonary stretch receptors. The mechanical model components include (i) three-element Hill muscle models of the diaphragm and abdominal muscles (Hill, 1938), (ii) a model of the larynx based on the results of Tully et al. (1990, 1991) and Rohrer’s (1915) equation, and (iii) lung/diaphragm/ribcage/abdomen volume relationships modeled on the data of Grassino et al. (1978) and the analysis of Mead and Loring (1982).

The first two equations represent the entire mechanical model. Each term is a function of the motor outputs of the network model (u_di, u_ab, and u_lm), and the diaphragm and abdominal wall volumes (V_di and V_ab) and their time derivatives ($\dot{V}$ _di and $\dot{V}$ _ab), and is defined by the subsequent equations. The parameters referenced in the model equations are listed in Table 1.

TABLE 1

Table 1. Parameters used in the biomechanical model.

Pressure balance on the rib cage

In Eq. 1, P_pl is the pleural pressure seen by the interior surface of the rib cage, P_ab−pl corrects for the fact that the rib cage sees abdominal pressure in the zone of apposition, F_diσ_di is the equivalent pressure due to the insertional forces of the diaphragm on the lower ribs, P_ica is the equivalent pressure due to the intercostal and accessory muscles, and σ_rc is the recoil pressure of the rib cage that balances the sum of the other pressures.

$P_{\text{pl}}+P_{\text{ab}-\text{pl}}+F_{\text{di}}{\sigma }_{\text{di}}+P_{\text{ica}}={\sigma }_{\text{rc}}$ (1)

Pressure balance on the diaphragm

In Eq. 2, σ_ab is the abdominal pressure, the excess of which over the pleural pressure must be balanced by σ_di, the recoil pressure of the diaphragm.

${\sigma }_{\text{ab}}-P_{\text{pl}}={\sigma }_{\text{di}}$ (2)

Pleural pressure

In Eq. 3, tracheal flow ($\dot{V}$ _L) is the derivative of lung volume. R_rs is the resistance of the airway, which when multiplied by $\dot{V}$ _L gives the pressure drop, by the hydraulic analog of Ohm’s law. σ_L is the recoil pressure of the lung.

$P_{\text{pl}}=-R_{\text{rs}}{\dot{V}}_{\text{L}}-{\sigma }_{\text{L}}$ (3)

Lung volume

Lung volume is a function of V_di and V_ab. Diaphragm volume, V_di, is the volume above the level of the diaphragm insertions on the rib cage and below the dome of the diaphragm. Abdominal wall volume, V_ab, is the volume between the abdominal wall and a frontal plane that coincides with the contracted position of the abdominal wall.

It has been commonly assumed that the abdominal contents are incompressible and that the abdominal cavity has only two movable walls – the diaphragm and the abdominal wall – and that therefore a displacement in one must be met by an equal and opposite displacement of the other (Grimby et al., 1976; Grassino et al., 1978; Macklem et al., 1978; Loring and Mead, 1982; Lichtenstein et al., 1992; Fitz-Clarke, 2007); in other words, that V_di + V_ab = V_sum (V_sum constant). Under this assumption, the abdominal volume completely determines the volume under the diaphragm, V_di. Following Lichtenstein et al. (1992), V_di determines the static contractile pressure generated by the diaphragm at a given activation. However, experimental data (Grassino et al., 1978, Figure 4) show that both rib cage volume and abdominal volume affect this pressure. Therefore, in our model, we added a term, C₁ V_rc, to the equation, allowing the rib cage and the abdominal wall to independently alter the volume under the diaphragm, effectively adding a third movable wall to the abdominal container, as discussed by Mead and Loring (1982).

The value of C₁ was determined by fitting our model to published data (Grassino et al., 1978, Figure 4) giving diaphragm pressure as a function of rib cage and abdominal volume at a fixed diaphragm activation. Given a value of C₁, our model equations will calculate the diaphragm pressure from the volumes. We then found the value of C₁ that was the best fit of the calculated pressure values to the measured values. The resulting equation is V_di + C₁ V_rc + V_ab = V_sum, which together with Eq. 20, gives us Eq. 4 for lung volume (V_L).

$V_{\text{L}}=\frac{V_{\text{sum}}-\left(1+C_1\right)V_{\text{di}}-V_{\text{ab}}-C_1{V}_{\text{c}}}{C_{\text{1}}}$ (4)

Airway resistance

Rohrer’s equation was used to calculate airway resistance (Rohrer, 1915; Hey and Price, 1982). The equation Pressure = K₁ · Flow + K₂ · Flow² or, dividing through by flow, Resistance = K₁ + K₂ |Flow| was applied twice, once for laryngeal resistance (k₁ + k₂|$\dot{V}$ _L|) and once for the resistance of the oropharynx and lower airway (0.72 + 0.44|$\dot{V}$ _L|), to give Eq. 5. We used a value of K₂ = 0.44 for the oropharynx and lower airway based on the assumption that K₂ is 0 for the oropharynx (Renotte et al., 1998; Eq. 8), and 0.44 for the lower airway (Renotte et al., 1998, Table 2). Assuming K₁ for the larynx is negligible during quiet breathing (see below), we used 0.34 for the lower airway and 0.38 for the oropharynx (Renotte et al., 1998, Table 2), for a total of 0.72.

The values of k₁ and k₂ for the larynx depend on the diameter of the glottis, as shown in Eqs 6–8. The parameter D is the diameter of the human trachea. The variable u_lm is the net laryngeal muscle activation, ranging from −1 (closed glottis) to 1 (open glottis). The variable d is the diameter of the glottis, or more precisely, the diameter of a circle with the same area as the opening of the glottis. The resting diameter (when u_lm = 0) is taken to be 10.9 mm (Baier et al., 1977; Brancatisano et al., 1983; D’Urzo et al., 1988), and is assumed (see Eq. 8) to change in proportion to u_lm (Tully et al., 1990, 1991).

The coefficient k₂ given by Eq. 7 is calculated using the equation for flow through an orifice (Simpson, 1968, Table II). The value of k₂ for the upper airway is different for inspiration and expiration (Renotte et al., 1998, Table 2), which we assume is due to changes in d. Assuming equal excursions from the resting diameter, we solved for the coefficient that would give us the reported values of k₂, and got 0.167. This approach resulted in a resting value for k₂ of 0.681.

We calculated the ratio of the mean resting value of k₁ to the mean resting value of k₂ (Tully et al., 1990, Table 2), and multiplied the ratio by 0.681 to get a resting value of k₁ of 0.0035 for the larynx, which is small relative to K₁ for the rest of the airway, justifying our assumption above that K₁ for the larynx is negligible during quiet breathing.

The Darcy–Weisbach equation for pressure loss in a pipe and the Darcy friction factor for laminar flow tell us that the resistance is proportional to 1/d⁴ (Kreith et al., 2004), which we use in Eq. 6 for k₁. Plugging in the resting value of k₁ determined above and the resting value of d, we solved for the coefficient, which gives us 49.6.

$\begin{array}{rcll}{R}_{\text{rs}}& =k_1+0.72+\left(k_2+0.44\right){\dot{V}}_{\text{L}}& & \text{(5)}\\ k_1& =\frac{49.6}{d^4}& & \text{(6)}\\ k_2& =0.167\left(\frac{D^2}{d^4}+\frac{d^2}{D^2}-1\right)& & \text{(7)}\\ d& =\left\{\begin{array}{cc}D,\hfill & u_{\text{lm}}>\frac{71}{109}\hfill \\ 10.9\left(1+u_{\text{lm}}\right),\hfill & -1\le u_{\text{lm}}\le \frac{71}{109}\hfill \\ 0,\hfill & u_{\text{lm}}<-1\hfill \\ \hfill \end{array}\right.& & \text{(8)}\end{array}$

Abdominal pressure on the rib cage

A fraction f_a of the rib cage is exposed to abdominal pressure rather than pleural pressure. The recoil pressure of the diaphragm, σ_di, is the difference between abdominal pressure and pleural pressure, so P_{ab − pl} adjusts the pressure seen by the rib cage for this difference.

$P_{\text{ab}-\text{pl}}=f_{\text{a}}{\sigma }_{\text{di}}$ (9)

Abdominal fraction of the rib cage

At total lung capacity (TLC), none of the diaphragm is apposed to the rib cage (Mead and Loring, 1982), so we assumed that at all lung volumes, a portion of V_di equal to V_di at TLC $\left(V_{\text{di}}^{\text{TLC}}\right)$ does not contribute to the zone of apposition. The remainder $\left(V_{\text{di}}-V_{\text{di}}^{\text{TLC}}\right)$ is divided by the remainder plus the lung volume $\left(V_{\text{di}}-V_{\text{di}}^{\text{TLC}}-V_{\text{L}}\right)$ to give an estimate of the fraction of the rib cage surface above the diaphragm insertions that is exposed to abdominal pressure.

The “obligatory ring” below the diaphragm insertions, which is about 15% of the rib cage surface (Mead and Loring, 1982), is always exposed to abdominal pressure, and is represented by $f_{\text{a}}^{\text{TLC}}$ in Eq. 10. Our rib cage volume, V_rc, is the volume above a plane through the diaphragm attachments (Cluzel et al., 2000, Figure 3). The estimate of 15% of the rib cage surface was in the context of a different definition of rib cage volume in which a change in rib cage volume is equal to the change in lung volume with the abdominal wall held still (Konno and Mead, 1967). This alternative definition implies a larger volume for the rib cage because it includes parts below the diaphragm insertions. From our volume equations, this means that rib cage volume is larger than ours by a factor of (1 + C₁), so we divided the previously calculated fraction of the smaller rib cage by (1 + C₁) to turn it into a fraction of the larger rib cage before adding it to $f_{\text{a}}^{\text{TLC}}.$

$f_{\text{a}}=\frac{1}{1+C_1}\left(\frac{V_{\text{di}}-V_{\text{di}}^{\text{TLC}}}{V_{\text{di}}-V_{\text{di}}^{\text{TLC}}+V_{\text{L}}}\right)+f_{\text{a}}^{\text{TLC}}$ (10)

Recoil pressure of the lung

Equation 11 assumes a linear relationship between lung volume and recoil pressure. C_L is lung compliance. V_L0 is the lung volume at zero recoil pressure, which we took to be equal to the residual volume (RV) after a maximal exhalation, the small recoil pressure remaining at RV being close enough to zero for our purposes (Permutt and Martin, 1960; Harris, 2005).

${\sigma }_{\text{L}}=\frac{V_{\text{L}}-V_{\text{L0}}}{C_{\text{L}}}$ (11)

Recoil pressure of the diaphragm

In Eq. 12, the term $u_{\text{di}}{\sigma }_{\text{di}}^{\text{max}}{F}_{\text{fl}}^{\text{di}}{F}_{\text{fv}}^{\text{di}}$ corresponds to the Hill muscle model (Ratnovsky et al., 2003, Eq. A6), except that ${\sigma }_{\text{di}}^{\text{max}}$ is a pressure rather than a force and $F_{\text{fl}}^{\text{di}}$ and $F_{\text{fv}}^{\text{di}}$ are functions of volume and its derivative (flow), respectively, rather than length and velocity.

We substituted pressure for force because, by Laplace’s (1808) Law, the pressure is proportional to the force when the curvature is constant, and the curvature of the human diaphragm dome does not change significantly with volume (Braun et al., (1982). Moreover, there is a constant ratio between diaphragm force and pressure in the dog (Kim et al., 1976).

Our substitution of volume for length (with an offset) is supported by the observation that the relationship between diaphragm pressure and length is not clearly different from linear when measured at RV, functional residual capacity (FRC), and TLC (Cluzel et al., 2000). To the extent that the action of the diaphragm resembles that of a piston (Kim et al., 1976), this linearity is expected. There is precedent for a Hill-style model in terms of pressure and flow for the respiratory system (Younes and Riddle, 1981).

In Eq. 12, u_di is phrenic activation of the diaphragm; ${\sigma }_{\text{di}}^{\text{max}}$ is static diaphragm recoil pressure at optimum length and maximum activation; R_di $\dot{V}$ _di is the pressure due to the passive resistance of the diaphragm.

${\sigma }_{\text{di}}=u_{\text{di}}{\sigma }_{\text{di}}^{\text{max}}{F}_{\text{fl}}^{\text{di}}{F}_{\text{fv}}^{\text{di}}+{\sigma }_{\text{di}}^{\text{psv}}+R_{\text{di}}{\dot{V}}_{\text{di}}$ (12)

Passive recoil pressure of the diaphragm

In Eq. 13, ${\sigma }_{\text{di}}^{\text{psv}}$ is the passive transdiaphragmatic pressure as a function of diaphragm volume. This pressure is taken to be zero at resting diaphragm volume $V_{\text{di}}^{\text{FRC}}$ (Agostoni and Rahn, 1960) and below, and quadratic above (Reid et al., 1987).

${\sigma }_{\text{di}}^{\text{psv}}=\left\{\begin{array}{cc}{K}_{\text{di}}^{\text{psv}}.{\left(V_{\text{di}}-V_{\text{di}}^{\text{FRC}}\right)}^2,\hfill & V_{\text{di}}>V_{\text{di}}^{\text{FRC}}\hfill \\ 0,\hfill & V_{\text{di}}\le V_{\text{di}}^{\text{FRC}}\hfill \\ \hfill \end{array}\right.$ (13)

Volume-pressure relationship of the diaphragm

In Eq. 14, $F_{\text{fl}}^{\text{di}}$ is the static pressure-volume relationship of the diaphragm (corresponding to Ratnovsky et al., 2003, Eq. A7 with the relative length replaced by a linear function of volume as described above). The parameter V_di0 is the volume under the diaphragm at the resting length. This is taken to be equal to the diaphragm volume at RV, based on the observation that the “highest Pdi twitch amplitude was recorded at RV” (Smith and Bellemare, 1987). The parameter $L_{\text{di}}^{min}$ is the length of the diaphragm at zero volume (i.e., when the diaphragm is flat) divided by the resting length, and is calculated by assuming that the length of the diaphragm at TLC is 65% of that at RV (Smith and Bellemare, 1987).

$F_{\text{fl}}^{\text{di}}=exp\left(-0.5{\left(\frac{\frac{1-L_{\text{di}}^{\text{min}}}{V_{\text{di0}}}{V}_{\text{di}}+L_{\text{di}}^{\text{min}}-1.05}{0.19}\right)}^2\right)$ (14)

Pressure-flow relationship of the diaphragm

In Eq. 15, $F_{\text{fv}}^{\text{di}}$ is the pressure-flow relationship of the diaphragm, with the velocity replaced by flow as discussed above (Hatze, 1981; Rosen et al., 1999; Artemiadis and Kyriakopoulos, 2005). The variable $\dot{V}$ _di is the time derivative of the volume under the diaphragm. The maximum rate of change of diaphragm volume, ${\dot{V}}_{\text{di}}^{\text{max}},$ was derived from data which gives transdiaphragmatic pressure as a function of flow at several levels of diaphragm activation up to 45% (Goldman et al., 1978, Figure 6). Because the rib cage was held still, the flow represents the rate of change of diaphragm volume. Fitting the curves to a Hill-style relation between pressure and flow (Younes and Riddle, 1981), there is a maximum flow (where the pressure goes to zero) for each level of diaphragm activation. Experimental results suggest that the maximum flow increases somewhat linearly to 80% activation and then levels off (Chow and Darling, 1999). We used that assumption together with the results for the maximum flow at lower activations to compute a maximum flow at 100% activation, which we use for ${\dot{V}}_{\text{di}}^{\text{max}}.$

$F_{\text{fv}}^{\text{di}}=\frac{0.1433}{0.1074+exp\left(-1.409\text{sinh}\left(3.2\frac{{\dot{V}}_{\text{di}}}{{\dot{V}}_{\text{di}}^{max}}+1.6\right)\right)}$ (15)

Pressure of the intercostal and accessory muscles

In Eq. 16, P_ica is the effective pressure generated by the intercostal and accessory muscles; positive values act to expand the rib cage.

In Eq. 17, $P_{\text{ica}}^{\text{di}}$ is the pressure due to the action of the inspiratory intercostals, which are assumed to be inactive when the diaphragm volume is above $V_{\text{di}}^{\text{FRC}}$ (low lung volumes). Below $V_{\text{di}}^{\text{FRC}},$ the pressure exerted by the inspiratory intercostals is assumed to be proportional to the activation of the diaphragm (u_di), and the proportionality constant itself is assumed to scale linearly from 0 at $V_{\text{di}}^{\text{FRC}}$ to its maximum value of $P_{\text{ica}}^{\text{diTLC}}$ at $V_{\text{di}}^{\text{TLC}}.$ The parameter $P_{\text{ica}}^{\text{diTLC}}$ was calculated as the pressure necessary to complete the pressure balance on the rib cage at TLC.

In Eq. 18, $P_{\text{ica}}^{\text{ab}}$ is the pressure due to the action of the expiratory intercostals, which is assumed to be proportional to the abdominal muscle activation (u_ab); the proportionality constant itself is assumed to scale linearly with the rib cage volume, changing from $P_{\text{ica}}^{\text{abRV}}$ at RV to $P_{\text{ica}}^{\text{abTLC}}$ at TLC. In Eq. 18, $P_{\text{ica}}^{\text{abTLC}}$ is the pressure due to the expiratory intercostals at TLC and maximal abdominal activation. This parameter’s value was calculated by taking the mean male maximal mouth pressure at TLC (from Ratnovsky et al., 2008, Table 1) and subtracting it from the rib cage recoil pressure at TLC. $P_{\text{ica}}^{\text{abRV}}$ is the pressure due to the expiratory intercostals at RV and maximal abdominal activation. We calculated this parameter by solving the model equations for P_ica while assuming RV conditions. This gives us the intercostal pressure necessary to reach RV.

$\begin{array}{rcll}{P}_{\text{ica}}& =P_{\text{ica}}^{\text{di}}+P_{\text{ica}}^{\text{ab}}& & \text{(16)}\\ P_{\text{ica}}^{\text{di}}& =\left\{\begin{array}{cc}{u}_{\text{di}}{P}_{\text{ica}}^{\text{diTLC}}\frac{V_{\text{di}}-V_{\text{di}}^{\text{FRC}}}{V_{\text{di}}^{\text{TLC}}-V_{\text{di}}^{\text{FRC}}},\hfill & V_{\text{di}}<V_{\text{di}}^{\text{FRC}}\hfill \\ \text{0,}\hfill & V_{\text{di}}\ge V_{\text{di}}^{\text{FRC}}\hfill \\ \hfill \end{array}\right.& & \text{(17)}\\ P_{\text{ica}}^{\text{ab}}& =u_{\text{ab}}\left(P_{\text{ica}}^{\text{abRV}}+\frac{V_{\text{rc}}-V_{\text{rc}}^{\text{RV}}}{V_{\text{rc}}^{\text{TLC}}-V_{\text{rc}}^{\text{RV}}}\left(P_{\text{ica}}^{\text{abTLC}}-P_{\text{ica}}^{\text{abRV}}\right)\right)& & \text{(18)}\end{array}$

Recoil pressure of the rib cage

The volume of the rib cage is assumed to be a sigmoid function of the recoil pressure of the rib cage, σ_rc. With increasing pressure the volume asymptotically approaches a maximum $\left(V_{\text{rc}}^{\text{max}}\right),$ and with decreasing pressure it asymptotically approaches a minimum $\left(V_{\text{rc}}^{\text{min}}\right).$ A generalized logistic function is used for the sigmoid, giving V_rc as a function of σ_rc; that equation is solved for σ_rc to give the first part of Eq. 19. The final term of Eq. 19 is the pressure due to the passive resistance of the rib cage (R_rc) and the rate of change of its volume ($\dot{V}$ _rc) The parameter ${\sigma }_{\text{rc}}^{\text{mul}}$ controls the maximum slope of the sigmoid; the slope is the compliance of the rib cage. It is calculated to make the compliance equal to C_rc/(1 + C₁). The factor of (1 + C₁) appears because C_rc is for a rib cage volume defined differently than V_rc. C_rc is the compliance of the rib cage, an average of values for three sitting subjects (Gilroy et al., 1985, Table 1).

${\sigma }_{\text{rc}}={\sigma }_{\text{rc}}^{\text{mul}}log\left(\frac{V_{{}_{\text{rc}}}^{max}-V_{\text{rc}}}{V_{\text{rc}}-V_{\text{rc}}^{min}}\right)+{\sigma }_{\text{rc}}^{\text{add}}+R_{\text{rc}}{\dot{V}}_{\text{rc}}$ (19)

Volume of the rib cage

The rib cage volume (V_rc) is the sum of the lung volume (V_L), the volume under the diaphragm (V_di), and the volume of the mediastinum and the lung blood and tissue (V_c).

$V_{\text{rc}}=V_{\text{L}}+V_{\text{di}}+V_{\text{c}}$ (20)

Recoil pressure of the abdominal wall

The abdominal wall is modeled as a surface with a circular segment cross-section in each transverse plane, all with the same radius, and a circular segment cross-section in each sagittal plane, all with another radius. The volume behind the abdominal wall, V_ab, is bounded by this surface and by a frontal plane. The values for the sagittal (r_s) and transverse (r_t) radii were derived from measurements taken during insufflation for laparoscopic surgery in humans (see Figure 3, Song et al., 2006). We fit exponential curves to the data points and the resulting relationship between the fitted sagittal and transverse radii was found to be well approximated by a linear function: r_s = 8.00 r_t − 1.10. The length of the longest transverse chord in the bounding frontal plane (c_t) was found which gave the volume change stated in the paper.

In Eq. 21, $u_{\text{ab}}{F}_{\text{CEmax}}{F}_{\text{fl}}^{\text{ab}}{F}_{\text{fv}}^{\text{ab}}$ is the Hill muscle model equation (Ratnovsky et al., 2003, Eq. A6); u_ab is the activation of the diaphragm by the lumbar motor neurons; F_CEmax is the maximal force capacity of the contractile element for a 1.5 cm² cross-section of canine external oblique muscle (Ratnovsky et al., 2003, Table 1). The constant k converts from a force to a surface tension, and (1/r_s + 1/r_t) converts the surface tension to a pressure, using the Law of Laplace (Laplace, 1808). The second term on the right, (V_ab − V_ab0)/C_ab, is the passive recoil pressure of the abdominal wall. V_ab0 is the volume behind the abdominal wall at which the recoil pressure is 0. This was taken to be equal to $V_{\text{ab}}^{\text{FRC}},$ since we assume a supine position. C_ab is the compliance of the abdominal wall. The final term is the pressure due to the passive resistance of the abdominal wall (R_ab) and the rate of change of abdominal volume ($\dot{V}$ _ab).

${\sigma }_{\text{ab}}=u_{\text{ab}}{F}_{\text{CEmax}}{F}_{\text{fl}}^{\text{ab}}{F}_{\text{fv}}^{\text{ab}}.\left(\frac{k}{r_{\text{t}}}+\frac{k}{r_{\text{s}}}\right)+\frac{V_{\text{ab}}-V_{\text{ab0}}}{C_{\text{ab}}}+R_{\text{ab}}{\dot{V}}_{\text{ab}}$ (21)

In Eq. 22, $F_{\text{fl}}^{\text{ab}}$ is the static force-length relationship of the abdominal wall (Ratnovsky et al., 2003, Eq. A7); L_CE is the length of the transversus abdominis; L_CE0 is its resting length.

$F_{\text{fl}}^{\text{ab}}=exp\left(-0.5{\left(\frac{\frac{L_{\text{CE}}}{L_{\text{CE0}}}-1.05}{0.19}\right)}^2\right)$ (22)

In Eq. 23, $F_{\text{fv}}^{\text{ab}}$ is the force-velocity relationship of the abdominal wall muscles (Hatze, 1981; Rosen et al., 1999; Artemiadis and Kyriakopoulos, 2005). The variable ${\dot{L}}_{\text{CE}}$ is the velocity of the contractile element (the time derivative of L_CE) and the parameter V_CEmax is the maximal contractile velocity of canine external oblique muscle (Ratnovsky et al., 2003, Table 1).

$F_{\text{fv}}^{\text{ab}}=\frac{0.1433}{0.1074+exp\left(-1.409sinh\left(3.2\frac{{\dot{L}}_{\text{CE}}}{V_{\text{CEmax}}}+1.6\right)\right)}$ (23)

Length-volume relationship of the abdominal wall

Equations 24 through 31 calculate the length of the abdominal muscle (L_CE) from the volume behind the abdominal wall (V_ab). Eqs 26 through 31 calculate V_ab as a function of r_t (the transverse radius); Eq. 25 uses the inverse of the resulting function to calculate r_t from V_ab; Eq. 24 calculates L_CE from r_t. In practice, the function V_ab(r_t) is pre-calculated, and approximated and inverted with a spline, and the spline is used to evaluate r_t(V_ab) during simulation.

$\begin{array}{rcll}{L}_{\text{CE}}& =100r_{\text{t}}{sin}^{-1}\left(\frac{c_{\text{t}}}{2r_{\text{t}}}\right)& & \text{(24)}\\ r_{\text{t}}\left(V_{\text{ab}}\right)& =\left\{\begin{array}{cc}{V}^{-1}\left(V_{\text{ab}}\right),\hfill & V_{\text{ab}}>0\hfill \\ 0.5c_{\text{t}},\hfill & V_{\text{ab}}\le 0\hfill \\ \hfill \end{array}\right.& & \text{(25)}\\ V\left(r_{\text{t}}\right)& =\underset{-c_{\text{s}}∕2}{\overset{c_{\text{s}}∕2}{\int }}A\left(r_{\text{t}},y\right)dy& & \text{(26)}\end{array}$

In Eq. 27, the Pythagorean Theorem is applied in the midsagittal plane to get c_s, the length of the chord that connects the ends of the abdominal wall in the frontal plane that serves as a boundary of the abdominal wall volume.

$c_{\text{s}}=2\sqrt{{r_{\text{s}}}^2-{\left(h_0-r_{\text{s}}\right)}^2}$ (27)

In Eq. 28, the Pythagorean Theorem is applied in a transverse plane to get h₀, the distance from the peak of the abdominal wall to the frontal plane that serves as a boundary of the abdominal wall volume.

$h_{\text{0}}=r_{\text{t}}-\sqrt{r_{\text{t}}^2-{\left(\frac{c_{\text{t}}}{2}\right)}^2}$ (28)

In Eq. 29, the formula for the area of a circular segment is applied in the transverse plane at a distance y from the peak of the abdominal wall to get the area between the abdominal wall and the boundary frontal plane.

$A=\frac{r_{\text{t}}^2}{2}\left(2{cos}^{-1}\left(1-\frac{h}{r_{\text{t}}}\right)+sin\left(2{cos}^{-1}\left(\frac{h}{r_{\text{t}}}-1\right)\right)\right)$ (29)

In Eq. 30, the Pythagorean theorem is applied in the midsagittal plane to get h, the distance from the abdominal wall to the boundary frontal plane at a distance y in the craniocaudal direction from the peak of the abdominal wall.

$h=\sqrt{r_{\text{s}}^2-y^2}-\sqrt{r_s^2-{\left(\frac{c_{\text{s}}}{2}\right)}^2}$ (30)

Equation 31 is the relation between the sagittal radius (r_s) and the transverse radius (r_t) derived from the results in Song et al. (2006).

$r_{\text{s}}=8.00479r_{\text{t}}-1.10158$ (31)

Equations for calculated parameters

$\begin{array}{rcll}& {\sigma }_{\text{ab}}\left(u_{\text{ab}}=1,V_{\text{ab}}^{\text{RV}},{\dot{V}}_{\text{ab}}=0,V_{\text{ab0}}\right)={\sigma }_{\text{ab}}^{\text{RV}}& & \text{(32)}\end{array}$

$\begin{array}{rcll}& V_{\text{ab0}}=V_{\text{ab}}^{\text{RV}}+V_{\text{ab}}^{\text{kmFRC}}\cdot \text{VC}& & \text{(33)}\end{array}$

$\begin{array}{rcll}& V_{\text{L}}^{\text{RV}}=V_{\text{L}}^{\text{FRC}}-\left(V_{\text{rc}}^{\text{kmFRC}}+V_{\text{ab}}^{\text{kmFRC}}\right)\cdot \text{VC}& & \text{(34)}\end{array}$

$\begin{array}{rcll}& {\sigma }_{\text{L}}^{\text{RV}}=\frac{V_{\text{L}}^{\text{RV}}-V_{\text{L0}}}{C_{\text{L}}}& & \text{(35)}\end{array}$

$\begin{array}{rcll}& {\sigma }_{\text{ab}}^{\text{RV}}={\sigma }_{\text{di}}^{\text{RV}}-{\sigma }_{\text{L}}^{\text{RV}}& & \text{(36)}\end{array}$

$\begin{array}{rcll}& V_{\text{ab}}^{\text{FRC}}=V_{\text{ab}}^{\text{RV}}+V_{\text{ab}}^{\text{kmFRC}}\cdot \text{VC}& & \text{(37)}\end{array}$

$\begin{array}{rcll}& V_{\text{sum}}=V_{\text{di}}^{\text{FRC}}+C_{\text{1}}{V}_{\text{rc}}^{\text{FRC}}+V_{\text{ab}}^{\text{FRC}}& & \text{(38)}\end{array}$

$\begin{array}{rcll}& V_{\text{di}}^{\text{RV}}=V_{\text{sum}}-V_{\text{ab}}^{\text{RV}}-C_{\text{1}}{V}_{\text{rc}}^{\text{RV}}& & \text{(39)}\end{array}$

$\begin{array}{rcll}& V_{\text{ab}}^{\text{TLC}}=V_{\text{ab}}^{\text{FRC}}+\left(V_{\text{ab}}^{\text{kmTLC}}-V_{\text{ab}}^{\text{kmFRC}}\right)\cdot \text{VC}& & \text{(40)}\end{array}$

$\begin{array}{rcll}& V_{\text{L}}^{\text{TLC}}=V_{\text{rc}}^{\text{TLC}}-V_{\text{di}}^{\text{TLC}}-V_{\text{c}}& & \text{(41)}\end{array}$

$\begin{array}{rcll}& {\sigma }_{\text{L}}^{\text{TLC}}=\frac{V_{\text{L}}^{\text{TLC}}-V_{\text{L0}}}{C_{\text{L}}}& & \text{(42)}\end{array}$

$\begin{array}{rcll}& {\sigma }_{\text{ab}}^{\text{TLC}}=\frac{V_{\text{ab}}^{\text{TLC}}-V_{\text{ab0}}}{C_{\text{ab}}}& & \text{(43)}\end{array}$

$\begin{array}{rcll}& {\sigma }_{\text{di}}^{\text{TLC}}={\sigma }_{\text{L}}^{\text{TLC}}-{\sigma }_{\text{ab}}^{\text{TLC}}& & \text{(44)}\end{array}$

$\begin{array}{rcll}& {\sigma }_{\text{rc}}^{\text{TLC}}={\sigma }_{\text{rc}}^{\text{mul}}log\left(\frac{V_{\text{rc}}^{max}-V_{\text{rc}}^{\text{TLC}}}{V_{\text{rc}}^{\text{TLC}}-V_{\text{rc}}^{\text{min}}}\right)+{\sigma }_{\text{rc}}^{\text{add}}& & \text{(45)}\end{array}$

$\begin{array}{rcll}& {\sigma }_{\text{rc}}^{\text{RV}}={\sigma }_{\text{rc}}^{\text{mul}}log\left(\frac{V_{\text{rc}}^{\text{max}}-V_{\text{rc}}^{\text{RV}}}{V_{\text{rc}}^{\text{RV}}-V_{\text{rc}}^{\text{min}}}\right)+{\sigma }_{\text{rc}}^{\text{add}}& & \text{(46)}\end{array}$

$\begin{array}{rcll}& f_{\text{a}}^{\text{RV}}=\frac{1}{1+C_{\text{1}}}\left(\frac{V_{\text{di}}^{\text{RV}}-V_{\text{di}}^{\text{TLC}}}{V_{\text{di}}^{\text{RV}}-V_{\text{di}}^{\text{TLC}}+V_{\text{L}}^{\text{RV}}}\right)+f_{\text{a}}^{\text{TLC}}& & \text{(47)}\end{array}$

$\begin{array}{rcll}& V_{\text{L0}}=V_{\text{L}}^{\text{RV}}& & \text{(48)}\end{array}$

$\begin{array}{rcll}& K_{\text{di}}^{\text{psv}}=\frac{{\sigma }_{\text{di}}^{\text{RV}}}{{\left(V_{\text{di}}^{\text{RV}}-V_{\text{di}}^{\text{FRC}}\right)}^2}& & \text{(49)}\end{array}$

$\begin{array}{rcll}& L_{\text{di}}^{\text{min}}=\frac{V_{\text{di}}^{\text{TLC}}-f_{\text{di}}{V}_{\text{di}}^{\text{RV}}}{V_{\text{di}}^{\text{TLC}}-V_{\text{di}}^{\text{RV}}∕1.05}& & \text{(50)}\end{array}$

$\begin{array}{rcll}& P_{\text{ica}}^{\text{abRV}}={\sigma }_{\text{L}}^{\text{RV}}+{\sigma }_{\text{rc}}^{\text{RV}}-\left(f_{\text{a}}^{\text{RV}}+F_{\text{di}}\right){\sigma }_{\text{di}}^{\text{RV}}& & \text{(51)}\end{array}$

$\begin{array}{rcll}& P_{\text{ica}}^{\text{diTLC}}={\sigma }_{\text{L}}^{\text{TLC}}+{\sigma }_{\text{rc}}^{\text{TLC}}-\left(f_{\text{a}}^{\text{TLC}}+F_{\text{di}}\right){\sigma }_{\text{di}}^{\text{TLC}}& & \text{(52)}\end{array}$

$\begin{array}{rcll}& {\sigma }_{\text{di}}\left(u_{\text{di}}=1,V_{\text{di}}^{\text{TLC}},{\dot{V}}_{\text{di}}=0,{\sigma }_{\text{di}}^{\text{max}}\right)={\sigma }_{\text{di}}^{\text{TLC}}& & \text{(53)}\end{array}$

$\begin{array}{rcll}& {\sigma }_{\text{rc}}^{\text{add}}={\sigma }_{\text{rc}}^{\text{mul}}log\left(\frac{V_{\text{rc0}}-V_{\text{rc}}^{\text{min}}}{V_{\text{rc}}^{\text{max}}-V_{\text{rc0}}}\right)& & \text{(54)}\end{array}$

$\begin{array}{rcll}& {\sigma }_{\text{rc}}^{\text{mul}}=\frac{\left(V_{\text{rc}}^{\text{max}}-V_{\text{rc}}^{\text{min}}\right)\left(1+C_{\text{1}}\right)}{4C_{\text{rc}}}& & \text{(55)}\end{array}$

$\begin{array}{rcll}& V_{\text{di0}}=V_{\text{di}}^{\text{RV}}& & \text{(56)}\end{array}$

$\begin{array}{rcll}& V_{\text{di}}^{\text{TLC}}=V_{\text{sum}}-V_{\text{ab}}^{\text{TLC}}-C_{\text{1}}{V}_{\text{rc}}^{\text{TLC}}& & \text{(57)}\end{array}$

$\begin{array}{rcll}& V_{\text{rc}}^{\text{max}}=V_{\text{rc}}^{\text{TLC}}+0.05\left(V_{\text{rc}}^{\text{TLC}}-V_{\text{rc}}^{\text{RV}}\right)& & \text{(58)}\end{array}$

$\begin{array}{rcll}& V_{\text{rc}}^{\text{min}}=V_{\text{rc}}^{\text{RV}}-0.99\left(V_{\text{rc}}^{\text{TLC}}-V_{\text{rc}}^{\text{RV}}\right)& & \text{(59)}\end{array}$

$\begin{array}{rcll}& V_{\text{rc}}^{\text{RV}}=V_{\text{rc}}^{\text{FRC}}-\frac{V_{\text{rc}}^{\text{kmFRC}}\cdot \text{VC}}{1+C_{\text{1}}}& & \text{(60)}\end{array}$

$\begin{array}{rcll}& V_{\text{rc}}^{\text{TLC}}=V_{\text{rc}}^{\text{FRC}}+\frac{\left(V_{\text{rc}}^{\text{kmTLC}}-{\text{V}}_{\text{rc}}^{\text{kmFRC}}\right)\cdot \text{VC}}{1+C_{\text{1}}}& & \text{(61)}\end{array}$

$\begin{array}{rcll}& V_{\text{rc0}}=V_{\text{rc}}^{\text{FRC}}& & \text{(62)}\end{array}$

Brainstem Network Model Architecture and System Performance when Linked to the Biomechanical Model

The computational model of the pontomedullary respiratory network (Figure 1) instantiated the hypothesis (Shannon et al., 1998; Rybak et al., 2008) that airway cough receptors affect several neuron populations in the ventral respiratory column (VRC) and pontine respiratory group (PRG) via cough 2nd order neurons. Evoked changes reconfigured the respiratory network to produce the cough motor pattern through the same VRC neurons involved in providing drive to respiratory muscles during normal breathing. The model incorporated recent enhancements (Poliaček et al., 2011) and additional neuron populations and other changes as detailed in Tables 2–4.

FIGURE 1

Figure 1. Schematic of the raphé-pontomedullary respiratory network model used in this study. The model extends that used in Poliaček et al. (2011) and Rybak et al. (2008) and follows labeling conventions enumerated therein. Parameters for the represented cell populations (large colored circles) and connections (see Key) are listed in Tables 2 and 4. Parameters for the I-Driver population with conditional bursting pacemaker properties were as described previously (Rybak et al., 2008; Poliaček et al., 2011). Abbreviations of brainstem regions or “compartments”: pre-BötC, pre-Bötzinger complex; VRC or VRG, ventral respiratory column or group; PRG, pontine respiratory group. Abbreviations of most populations were as enumerated in Table 1 of Rybak et al. (2008): Aug and Dec: neurons with augmenting or decrementing activity patterns, respectively, during the indicated phase (I-inspiratory; E-expiratory) of maximum firing rate. BS, bulbo-spinal; ELM, expiratory laryngeal motoneurons; EI, neurons with a peak firing rate during the E-to-I phase transition; IE, neurons with a peak firing rate during the I-to-E phase transition; ILM, inspiratory laryngeal motoneurons; NRM, non-respiratory-modulated neurons. Two phrenic motor neuron populations with different threshold ranges innervated the diaphragm (PHR, PHR-HT: high threshold); two lumbar motor neuron populations activated the abdominal muscle (LUM, LUM-HT: high threshold). I-Dec_2, second inhibitory I population in the VRC (e.g., see Ott et al., 2012); Lung Def_1s, Lung Dis_1s, Deflation 2nd order: lung deflation-sensing neuronal circuit elements. See text for further discussion.

TABLE 2

Table 2. Population parameters for network model with adjusted (Δ) and additional (+) neuronal populations modified from Poliaček et al. (2011).

TABLE 3

Table 3. Synaptic parameters for the network model.

TABLE 4

Table 4. Connectivity for the network model modified from Poliaček et al. (2011).

Linking the neural network and biomechanical models

The diaphragm received input from two phrenic motor neuron populations (PHR, PHR-HT) with different threshold ranges to generate motor unit diversity and facilitate an ordered recruitment during increased inspiratory drive. Diaphragm activation was based on the mean instantaneous firing rates (P, P₁), of the two populations by the expression (0.3P + 0.7P₁)/X where X is the firing rate for maximum diaphragm activation; values of 50–200 spikes/s were used (Nail et al., 1972) to approximate the plot for diaphragm activation in Figure 1 of Mantilla and Sieck (2011). Similarly, two lumbar motor neuron populations (LUM, LUM-HT) activated the abdominal muscle with X set to 80. Inspiratory laryngeal motor (ILM) and expiratory laryngeal motor (ELM) neuron populations regulated laryngeal resistance over a range between fully open (+1) and fully closed (−1), inclusively.

Lung afferent populations were regulated by injected currents defined at each simulation step by evaluation of expressions that included model lung volume. Pulmonary stretch receptors (PSRs) became more active with increasing lung volume V (membrane bias = 0.5V mV/%VC) and mediated the Hering–Breuer reflex. Deflation-sensitive lung receptors were also implemented (Paintal, 1955; Luck, 1970; Wei and Shen, 1985; Iscoe and Gordon, 1992; Bergren and Peterson, 1993; Matsumoto et al., 2002; Yu et al., 2003). A low threshold population (Def_1, membrane bias = − 0.225(V − 70) mV/%VC) and its afferent pathway introduced in Poliaček et al. (2011) was used to represent a class of possible network mechanisms for generating an inhibitory bias on E-Dec neurons. Simulated “vagotomy” (elimination of the effects of lung afferents) in the present model removed this inhibition, contributing to the observed prolongation of the expiratory phase (Te increased from 2.76 s to 3.15 s (p-value = 0.0004, two-sided t-test). Vagotomy also removed the influence of the PSRs and increased inspiratory phase duration (Ti) from 1.94 s to 2.61 s, (p-value = 4 × 10⁻⁷; see references and discussion in Dick et al., 2008). A higher threshold “distortion” (Dis_1) receptor population (cf. Iscoe and Gordon, 1992) excited the ELM population when lung volume was below FRC (membrane bias = −1.75(V − 10) mV/%VC if V < 10, 0 otherwise). We note that the synaptic strength and firing rates of this speculative Dis_1 population, added for development purposes in other work (Hutchison and Lindsey, 2009), resulted in negligible modulation of ELM population activity under the conditions of the present study (see “Influence of Some Added Network Connections”).

Additional enhancements to the current network model

The “I-Dec_2” population was added to provide a second inhibitory VRC inspiratory neuron population for tuning inspiratory drive as proposed in Ott et al. (2012) for central chemoreceptor modulation of breathing. In some previous models (Rybak et al., 2008; Poliaček et al., 2011), the “E-Aug-late” population inhibited numerous target populations, but also served to excite the VRG bulbospinal E-Aug-BS (+) population that drives expiratory lumbar motor neurons. A new “E-Aug (+)” population was added to facilitate differential regulation of the lumbar motor neurons and expiratory drive modulation as proposed in the literature (Iscoe, 1998; Shannon et al., 2004; Molkov et al., 2010). Other parameters were adjusted and populations added in anticipation of linking the network model to the biomechanical model derived from data from human subjects. In the antecedent model (Poliaček et al., 2011), the I-Aug-BS population output served both a premotor function and represented the “phrenic” output. The inhibitory connections from the VRC-IE population to the I-Aug-BS population were eliminated in the new model; E-Dec-P inhibition of the I-Aug-BS population was retained. The resulting eupneic respiratory cycle frequency (12.7 cycles/min) was within the range for resting breathing in the human adult (A.D.A.M. Medical Encyclopedia, 2012).

Eupneic motor pattern and “baseline” cough

The joint neuromechanical model generated a eupneic motor pattern and an evoked cough. Figure 2 shows membrane potential records from simulated neurons in representative PRG, raphé, and VRC neuron populations and the six types of motor neurons. The “IF” neuron populations do not generate action potential-like spikes; instances of threshold crossings are indicated graphically by corresponding vertical spike-like lines. Additional traces include integrated population activity of the three lung afferent populations and biomechanical system metrics, including lung volume, tracheal flow, and alveolar pressure. The three phases of the cough cycle (Bolser et al., 2003) are highlighted.

FIGURE 2

Figure 2. Simulated eupneic respiratory cycles and an evoked cough motor pattern (inspiratory, compressive, and expulsive phases respectively labeled and highlighted by colored columns). The top 29 traces show membrane potentials and discharge patterns of individual respiratory neurons from a subset of the simulated populations as indicated by the labels on the left, arranged by region (PRG, raphé, VRC) or type (Motor neuron populations). The “integrated” phrenic trace represents the threshold crossing activity of the “PHR” population summed over 60 ms windows and indicates the inspiratory and expiratory phases of the respiratory cycle. Similarly, integrated traces from three lung afferent populations are plotted below the motor neuron records. (PSR, pulmonary stretch receptors) The 13 traces below those from the afferents show, in order from top to bottom: 1: lung volume (%VC, relative to RV); 2: tracheal flow (%VC/s, expiration positive (up)); 3: alveolar pressure (cmH₂O); 4–6: diaphragm activation, abdominal muscle activation, and net laryngeal muscle activation (dimensionless ratios to maximums); 7: diaphragm volume (L); 8: abdominal volume (L); 9: derivative of diaphragm volume (L/s); 10: derivative of abdominal volume (L/s); 11–13: transdiaphragmatic, abdominal, and transpulmonary pressures (cmH₂O). The bottom trace indicates the duration of a simulated cough stimulus. A fiber population composed of 100 fibers, each with a firing probability of 0.05 at each simulation time step and 100 type Ex_1 excitatory synaptic terminals (synaptic strength 0.03), represented cough receptor excitation. These fibers excited the Cough 2nd order neuron population (Figure 1); see Table 2 for properties of this population and Table 4 for details of connections with other populations. See text for further discussion.

The inspiratory phase of the cough was characterized by increased activation of the diaphragm and enlargement of upper airway via activation of the ILM (abductor) motor neuron population, resulting in an increased lung volume (43% VC), inspiratory flow, and transdiaphragmatic and transpulmonary pressures. The subsequent compressive phase included activation of the ELM (adductor) motor neurons with transient laryngeal closure, together with activation of lumbar motor neurons and abdominal expiratory muscles. During this phase, tracheal airflow stopped and there was an increase in alveolar and abdominal pressure. In the following expulsive phase, high air flow velocity (72.2% VC/s) resulted from the opening of the larynx during continued abdominal muscle activation.

Cough behavior with changes in inspiratory drive

Two series of simulations with complementary perturbations of cough inspiratory drive were made to assess model behavior during the phases of cough. Different sets of random number seeds were used for each simulation to generate variability in model output by altering the thresholds of individual neurons in each population and the convergent and divergent connectivity patterns among populations within ranges defined by the initial baseline parameter settings.

In the first series, the activation strengths for the connections between phrenic motor neuron populations and the diaphragm were increased by factors of 2 and 4. Figure 3 shows the integrated activity of each motor neuron population together with lung volume, tracheal flow, alveolar pressure, and abdominal pressure for baseline conditions (left) and a trial with four times the baseline activation strength (right).

FIGURE 3

Figure 3. Integrated traces of motor neuron population activities (top 6 traces, population labels on left) and biomechanical model outputs (labels on left; see legends of Figures 1 and 2 for definitions) during eupneic and “baseline” cough motor patterns before (A) and after (B) increasing the gain of diaphragm activation by a factor of 4 (see center schematic inset); otherwise, cough stimulus parameters as in Figure 2. Scales on right are for (A,B). See text for further details.

Outputs from four trials for each amplified condition were compared with each other and with the baseline results. Figure 4 shows the means (±95% confidence limits) of selected biomechanical outputs measured during baseline cough (1×) and the two conditions of increased activation gain (2×, 4×). Pairs of symbols connected by a line indicate no significant difference. Successively larger peak expiratory flow rates and abdominal pressures were respectively associated with greater lung volumes during preceding inspiratory phases of the evoked coughs, even though abdominal drive did not change. This result established that differences in flow with the generated changes in inspiratory (operating) volumes were the consequence of the modeled biomechanical system.

FIGURE 4

Figure 4. Means of peak values (±95% confidence limits) of (from top) lung volume, expiratory tracheal flow, alveolar pressure, abdominal pressure, and abdominal drive together with respiratory cycle frequencies during pre-cough eupneic intervals (bottom) measured during model simulations of baseline cough (1×) and two conditions of increased phrenic-to-diaphragm activation gain (2×, 4×). Pairs of symbols connected by a line indicate no significant difference.

Mean respiratory cycle frequencies measured during pre-cough eupneic intervals for each condition were also evaluated. The respiratory frequency increased with the highest (4×) inspiratory drive, a change associated with changes in feedback from lung afferents under the “closed loop” conditions evaluated.

In the second series of simulations, inspiratory drive was altered only during the cough cycle by changing synaptic strengths of Cough 2nd order neuron inputs to selected model populations. The top panels in Figures 5A_1–C1 show schematics of a subset of the model network and sites where synaptic strengths were changed relative to the baseline conditions represented in and described for Figure 3A. Corresponding panels in Figures 5A_2–C2 show integrated traces of motor neuron population activities and biomechanical model outputs for the respective perturbations; arrows highlight significant changes in the indicated metrics (further enumerated in Figure 6).

FIGURE 5

Figure 5. Network perturbations that selectively alter cough inspiratory drive change lung volume and also influence tracheal flow and abdominal pressure during the subsequent expiratory phase of cough. Each panel (A₁–C₁) shows a schematic of a subset of the model network and changes in synaptic strengths during simulated cough stimulation relative to the baseline conditions represented in and described for Figure 3A. Corresponding panels (A₂–C₂) show integrated traces of motor neuron population activities and biomechanical model outputs (left labels) for the respective perturbations of inspiratory drive. Arrows mark changes. See text for further details.

FIGURE 6

Figure 6. Means of peak values (±95% confidence limits) of (from top) lung volume, expiratory tracheal flow, alveolar pressure, abdominal drive, and abdominal pressure measured under conditions of pre-cough eupnea (Eupnea), baseline cough (Base), doubling of the synaptic strength between Cough 2nd order and the I-Aug and I-Aug-BS populations (A), elimination of Cough 2nd order excitation of the I-Aug and I-Aug-BS populations (B), and additionally, the subsequent elimination of I-Dec_2 neuron inhibition of E-Dec-Tonic neurons during the inspiratory phase of the evoked cough cycle (C). Pairs of symbols connected by a line indicate no significant difference.

First (Figure 5A), synaptic strengths from the Cough 2nd order population to the I-Aug and I-Aug-BS populations were doubled. The highlighted segment of the inset (Figure 5A₁) shows integrated traces for the I-Aug-BS, I-Dec_2, and E-Dec-Tonic populations during a eupneic cycle and the subsequent evoked cough. I-Aug-BS activity increased under this condition (asterisk).

Next (Figure 5B), cough inspiratory drive was decreased relative to baseline by deletion (synaptic strength = 0.0) of the excitatory connections between the Cough 2nd order population and both the I-Aug and I-Aug-BS populations. The excitation of the I-Dec_2 population by 2nd order cough neurons remained, partially suppressing the recurrent inhibition of the I-Dec_2 and two I-Aug populations. The highlighted segment of the inset (Figure 5B₁) shows reduced I-Aug-BS activity during the evoked cough under this condition (asterisk).

The third perturbation further reduced cough inspiratory drive by also transiently blocking I-Dec_2 neuron inhibition of the E-Dec-Tonic population during the cough cycle. The elimination of this influence resulted in increased E-Dec-Tonic inhibition of the I-Aug and I-Aug-BS populations during the cough (asterisk in highlighted segment of Figure 5C₁).

Figure 6 plots (from the top) the means of peak values (±95% confidence limits) for lung volume, expiratory tracheal flow, alveolar pressure, abdominal drive, and abdominal pressure measured under conditions of pre-cough eupnea (Eupnea), baseline cough (Base), and the three conditions represented in Figure 5. The differences in peak lung volumes during cough under the three conditions (A–C) confirm functional roles for both the excitatory and disinhibitory influences on inspiratory drive in the model. Deletion of the excitatory component of cough inspiratory drive (B) caused peak expiratory tracheal flow to decrease relative to the previous baseline and enhanced coughs (Base, A). However, peak abdominal drive and abdominal pressure did not change. A further reduction in peak lung volume to levels below eupneic control due to transiently increased E-Dec-Tonic inhibition of inspiratory drive (C) during cough resulted in no further change in expiratory flow, although peak abdominal pressure was reduced relative to baseline cough values.

A third series of simulations was done with the isolated biomechanical model. Figure 7 plots the peak expiratory flow in four coughs simulated at different operating volumes but equal peak abdominal pressure of 26.5 cmH₂O. In each cough, the diaphragm and abdominal activations were first controlled to produce the desired operating volume, then the laryngeal muscles were controlled to close the airway, then the abdominal activation was controlled to produce an abdominal pressure of 26.5 cmH₂O, and finally the laryngeal muscles were controlled to open the airway. Note that no statistics were done on these runs because the biomechanical model is deterministic. As in the first series of simulations, successively larger peak expiratory flow rates were associated with greater lung volumes during preceding inspiratory phases of the simulated cough, but unlike the first series, the peak abdominal pressure was the same in each cough. This result established that differences in flow with the generated changes in inspiratory (operating) volumes were not entirely due to the differences in pressure seen in the first series.

FIGURE 7

Figure 7. Peak expiratory flow in four coughs simulated with the isolated biomechanical model at different operating volumes but equal peak abdominal pressure of 26.5 cmH₂O. There are no error bars because these are runs of the deterministic model.

Comparisons with behaviors of antecedent models

Table 5 shows means of inspiratory and expiratory phase durations during eupnea and cough and peak firing rates of motor neuron populations common to the present model and those described in Rybak et al. (2008) and Poliaček et al. (2011). The current model has lower firing rates, similar to those observed in vivo (Iscoe, 1998; Mantilla and Sieck, 2011), and longer respiratory phase durations; inspiratory phase durations are also more variable (see coefficients of variation, Table 5).

TABLE 5

Table 5. Comparison with previous models.

These antecedent variants of the present neuronal network model were designed without a linked biomechanical system. The new joint neuromechanical model aids tuning of phase-timing relationships and the scaling of model motor outputs. To illustrate this model feature, we linked the earlier network models to the new biomechanical model. We note that the phrenic and lumbar motor neuron activities from the Rybak and Poliaček models are not strictly comparable because the current model has a second population of each type of motor neuron to model recruitment with increased drive. Lung stretch receptor inputs to the previous network models remained filtered versions of the phrenic motor output (i.e., there was no feedback from the mechanical models to the network). Figure 8 shows records of lung volume, alveolar pressure, tracheal flow, and laryngeal muscle activation from the current neuromechanical model (Figure 8A) and for the two earlier models when connected to the biomechanical system (Figures 8B,C). The scaling and activation of the laryngeal muscles caused airway closure prior to each eupneic expiration when using the older models’ outputs (lma = −1, flow flattens at 0). During cough in the previous models, the next inspiration started before the previous expiration was complete, resulting in a progressive increase in lung volume from cough to cough.

FIGURE 8

Figure 8. Outputs of the mechanical model when linked to the current network model (A) and the networks in (B) (Rybak et al., 2008) and (C) (Poliaček et al., 2011). The earlier networks were designed without a mechanical model, but were connected to the current mechanical model to generate these plots. The plots are lung volume, alveolar pressure, tracheal flow, and laryngeal muscle activation (lma). The value of lma is 1 for a maximally open glottis, 0 for the resting diameter, and −1 for a closed glottis. The first few cycles are eupneic cycles, and the rest are coughs. The time scale is different for the current model because it was designed with slower cycles to match human respiration.

Influence of some added network connections

As noted in Sections “Linking the Neural Network and Biomechanical Models” and “Additional Enhancements to the Current Network Model,” the current network model includes lung afferents responsive to lung deflation and presynaptic inhibition by E-Dec-Tonic neurons of excitatory inputs from the I-Aug population to I-Dec_2 neurons. The latter feature was added to prolong I-Dec_2 neuron activity when E-Dec-Tonic neuron I-phase activity is reduced. Removal of these three speculative model elements resulted in shorter inspiratory phase durations (Table 5, “No speculative” and “Current” columns).

Discussion

The new biomechanical model system detailed in the Results incorporates several features developed using measures from human subjects. These include a model of the abdominal volume that captures the interaction of the diaphragm, rib cage, and abdominal wall based on Grassino et al. (1978), an abdominal wall model based on measurements of the curvature of the abdomen by Song et al. (2006) taken during insufflation for laparoscopic surgery, and a model of the larynx using results from Tully et al. (1990, 1991). The mechanical model was linked to an enhanced version of a previously described computational network model (Rybak et al., 2008; Poliaček et al., 2011) with IF neuron populations, connections, and other properties measured or inferred from in vivo and in vitro studies of mammalian brainstem circuits for breathing and cough (Shannon et al., 2000; Segers et al., 2008; Lindsey et al., 2012).

The joint neuromechanical system is a prototype for study of the neural control of breathing and airway defensive behaviors. To our knowledge, computational neural network models of cough generation have been evaluated previously using measures of motor neuron burst sequences, phase durations, and the time varying firing rates of neuron populations that only indirectly reflect possible muscle activation patterns and their attendant biomechanical consequences. The new model generated eupneic breathing and cough motor patterns together with corresponding alterations in lung volume, tracheal air flow, and other relevant metrics of cough mechanics. The present results also show the utility of the model as an aid for tuning the motor pattern and as a tool to evaluate the efficacy of phase-timing relationships.

A specific goal of this project was to assess model output during cough under conditions of altered inspiratory drive. We were motivated in part by the recent observation that lung operating volume at the onset of the compressive phase of cough influences subsequent air flow velocities during the expulsive phase (Smith et al., 2012). Inspiratory drive was altered by two distinct approaches: (i) increased gain of phrenic motor neuron activation of the diaphragm, and (ii) sequential modulation or deletion of synaptic inputs to inspiratory premotor populations. Both perturbations altered cough inspiratory volume. We also found changes in expulsive phase air flow associated with corresponding changes in peak abdominal pressure attributable to cough mechanics, results that could not have been achieved by measures of the motor pattern output alone. In the first protocol, higher end inspiratory volumes resulted in greater tracheal air flow during the subsequent expulsive phase under the same abdominal expiratory motor drives. Under the second protocol, the difference in operating volumes between the enhanced drive and reduced excitatory drive states was associated with corresponding reductions in expiratory flow and peak abdominal pressure.

Discrepancies with Experimental Results and Model Limitations

The discrepancy between the present results and those of Smith et al. (2012) is noteworthy. The latter study found changes in expulsive flow rates during voluntary coughs from different operating volumes in the absence of significant alterations in thoracic or abdominal pressures, whereas we found changes in flow associated with changes in abdominal pressure, despite no change in abdominal drive. The change in expiratory pressure in the model is due to the action of the intercostal and accessory muscles; the expiratory pressure increases because the pressure from those muscles in the model increases with rib cage volume at constant activation. Our model calculates the expiratory pressure generated by the intercostal and accessory muscles at TLC and full abdominal activation necessary to produce maximal expiratory pressure, and at RV to complete the pressure balance on the rib cage (a much smaller number; see “Additional Enhancements to the Current Network Model”). The model assumes that the expiratory pressure generated by the intercostal and accessory muscles scales linearly with rib cage volume between RV and TLC, and linearly with abdominal activation, leading to higher pressures at higher rib cage volumes with equal activation. The experimentally observed increase in maximal expiratory pressure with rib cage volume could be due to increased activation of the intercostals or improved mechanical advantage at larger volumes or a combination of the two. If improved mechanical advantage is a factor, the brainstem would have to reduce drive at higher volumes during cough to avoid higher pressures, suggesting that it may be sensing the generated pressures and adjusting drive accordingly (see e.g., discussion in “Tonic Expiratory Neurons: Model Results and Predictions”). A refined configuration to accommodate separate intercostal muscles, intercostal motor neuron populations, and muscle afferents (Shannon, 1986) and their control of the chest wall would be useful in this regard.

When we ran the biomechanical model in isolation with reduced abdominal drive at higher volumes in order to maintain an equal peak pressure (see Figure 7), we saw peak flow rate changes comparable in magnitude to those seen by Smith et al. (2012, Figure 3B), due to the increasing recoil pressure of the lung at higher volumes. The peak flow rates were comparable despite the fact that the peak abdominal pressure in our simulation was less than half that observed by Smith et al. This lower resistance is likely due to the fact that we did not model airway collapse, which is known to limit peak flow rates (Knudson et al., 1974).

We found that increased flow during cough at higher lung volume is primarily due to increased lung recoil pressure. The lung recoil pressure certainly increases with lung volume, but the accuracy of the resulting flow in the model may be affected by certain known limitations of the model: (i) airway collapse during cough is not modeled, resulting in an underestimate of airway resistance; (ii) the lung compliance is assumed to be constant in the model, whereas it is thought to vary with lung volume in vivo; (iii) the model does not take into account hysteresis in the lung flow-volume curve; (iv) volume changes due to blood shift out of the trunk during cough are not modeled; and (v) volume changes due to gas compression are not modeled (see Smith et al., 2012 for data on volume changes due to blood shift and gas compression). Nevertheless, our model suggests the hypothesis that the increased flow during cough is primarily due to increased lung recoil pressure.

Tonic Expiratory Neurons: Model Results and Predictions

The model incorporated multiple target sites for cough drive modulation, a feature of the network architecture based on correlational linkage maps of functional connectivity and associated neuronal responses to stimuli that either enhance or suppress inspiratory drive in vivo (Lindsey et al., 1998; Shannon et al., 2000; Poliaček et al., 2011; Ott et al., 2012). Deletion of excitatory mechanisms for cough inspiratory drive resulted in reductions in peak lung volume and a subsequently diminished peak air flow relative to baseline during the expulsive phase (Figure 6B). Although removal of the disinhibitory component of the drive enhancement mechanism mediated by the E-Dec-Tonic population did not further reduce expulsive phase air flow velocity, it did lead to both an additional decrease in inspiratory phase lung operating volume and a reduced expulsive phase peak abdominal pressure relative to baseline values, even though peak abdominal drive did not change.

We have previously proposed the hypothesis that tonic expiratory neurons provide a reservoir for inspiratory drive modulation. Suppression of their inspiratory phase activity during central chemoreceptor-mediated drive and spike train cross-correlation analyses both suggest that VRC tonic E neuron inhibition of premotor inspiratory neurons is reduced in high drive states, at least in part, by increased I-Dec neuron inhibition (Ott et al., 2012).

The present model included a network “module” previously introduced (Poliaček et al., 2011) for baroreceptor modulation of breathing. That circuit, inferred from spike train correlational linkages and neuron responses to baroreceptor stimulation (Lindsey et al., 1998), operated via excitatory and disinhibitory raphé neuron influences acting upon VRC E-Dec-Phasic and E-Dec-Tonic neuron populations. Simulations of baroreceptor activation using prior models (Poliaček et al., 2011; Lindsey et al., 2012) with circuits inferred from in vivo observations (see references in Lindsey et al., 1998; Poliaček et al., 2011; Ott et al., 2012) generated prolongation of expiration and reduced inspiratory drive during both eupneic respiratory cycles and evoked cough.

Collectively, these data support the hypothesis that inhibition of the E-Dec-Tonic population in the cough network amplifies inspiratory drive via disinhibition. Experimental data consistent with this hypothesis is presented in a companion paper (Segers et al., 2012). Modulation of tonic expiratory neuron activity could also operate in a push-pull mechanism in which cough drive is balanced against the potentially suppressive influences of blood pressure changes caused by cough mechanics.

Future Directions

The present model provides a framework for integrating respiratory network model development with respiratory mechanics and will guide and facilitate scaling and timing of motor neuron activity patterns and functionally antecedent connectivity for the generation of breathing, cough, and swallow. The simulations of cough and breathing suggest that an important area of focus for future modeling efforts will be reconciliation of known differential effects of pulmonary volume-related feedback on breathing and airway defensive behaviors such as coughing and the expiratory reflex. Specific components of the model that are proposed to have the greatest effect on its potential for prediction are the gain of pulmonary volume-related feedback and the interaction of this feedback with cough-related sensory information. Future models should also guide experiments targeting the control of behavior that must be tightly coordinated with breathing, such as sniffing, swallowing, and vocalization.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by NIH grants HL89104, HL89071, and HL103415 from NHLBI, and grant NS19814 from NINDS.

References

A.D.A.M. Medical Encyclopedia. (2012). Rapid Shallow Breathing. Available at: http://www.nlm.nih.gov/medlineplus/ency/article/007198.htm [accessed May 9, 2012].

Agostoni, E., Gurtner, G., Torri, G., and Rahn, H. (1966). Respiratory mechanics during submersion and negative-pressure breathing. J. Appl. Physiol. 21, 251–258.

Pubmed Abstract | Pubmed Full Text

Agostoni, E., and Rahn, H. (1960). Abdominal and thoracic pressures at different lung volumes. J. Appl. Physiol. 15, 1087–1092.

Pubmed Abstract | Pubmed Full Text

Artemiadis, P. K., and Kyriakopoulos, K. J. (2005). “Teleoperation of a robot manipulator using EMG signals and a position tracker,” in IEEE/RSJ International Conference on Intelligent Robots and Systems, Edmonton, 2005.

Baekey, D. M., Morris, K. F., Nuding, S. C., Segers, L. S., Li, Z., Lindsey, B. G., and Shannon, R. (2001). Involvement of ventral respiratory group neurons in the fictive expiration reflex. FASEB J. abstr. 15, A798.

Baier, H., Wanner, A., Zarzecki, S., and Sackner, M. A. (1977). Relationships among glottis opening, respiratory flow, and upper airway resistance in humans. J. Appl. Physiol. 43, 603–611.

Pubmed Abstract | Pubmed Full Text

Barnas, G. M., Yoshino, K., Stamenovic, D., Kikuchi, Y., Loring, S. H., and Mead, J. (1989). Chest wall impedance partitioned into rib cage and diaphragm-abdominal pathways. J. Appl. Physiol. 66, 350–359.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bergren, D. R., and Peterson, D. F. (1993). Identification of vagal sensory receptors in the rat lung: are there subtypes of slowly adapting receptors? J. Physiol. 464, 681–698.

Pubmed Abstract | Pubmed Full Text

Bolser, D. C., Davenport, P. W., Golder, F. J., Baekey, D. M., Morris, K. F., Lindsey, B. G., and Shannon, R. (2003). “Neurogenesis of cough,” in Cough: Causes, Mechanisms and Therapy, 1st Edn, eds H. Boushey, K. Chung, and J. G. Widdicombe (Malden: Blackwell Publishing Ltd.), 173–180.

Brancatisano, T., Collett, P. W., and Engel, L. A. (1983). Respiratory movements of the vocal cords. J. Appl. Physiol. 54, 1269–1276.

Pubmed Abstract | Pubmed Full Text

Braun, N. M., Arora, N. S., and Rochester, D. F. (1982). Force-length relationship of the normal human diaphragm. J. Appl. Physiol. 53, 405–412.

Pubmed Abstract | Pubmed Full Text

Breen, B. J., Gerken, W. C., and Butera, R. J. (2003). Hybrid integrate-and-fire model of a bursting neuron. Neural Comput. 15, 2843–2862.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cheng, L., Ivanova, O., Fan, H.-H., and Khoo, M. C. K. (2010). An integrative model of respiratory and cardiovascular control in sleep-disordered breathing. Respir. Physiol. Neurobiol. 174, 4–28.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cheng, L., and Khoo, M. C. (2012). Modeling the autonomic and metabolic effects of obstructive sleep apnea: a simulation study. Front. Physiol. 2:111. doi:10.3389/fphys.2011.00111

CrossRef Full Text

Chiel, H. J., Ting, L. H., Ekeberg, Ö., and Hartmann, M. J. Z. (2009). The brain in its body: motor control and sensing in a biomechanical context. J. Neurosci. 29, 12807–12814.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Chow, J. W., and Darling, W. G. (1999). The maximum shortening velocity of muscle should be scaled with activation. J. Appl. Physiol. 86, 1025–1031.

Pubmed Abstract | Pubmed Full Text

Cluzel, P., Similowski, T., Chartrand-Lefebvre, C., Zelter, M., Derenne, J.-P., and Grenier, P. A. (2000). Diaphragm and chest wall: assessment of the inspiratory pump with MR imaging – preliminary observations. Radiology 215, 574–583.

Pubmed Abstract | Pubmed Full Text

De Troyer, A., Estenne, M., Ninane, V., Van Gansbeke, D., and Gorini, M. (1990). Transversus abdominis muscle function in humans. J. Appl. Physiol. 68, 1010–1016.

Pubmed Abstract | Pubmed Full Text

Dick, T. E., Shannon, R., Lindsey, B. G., Nuding, S. C., Segers, L. S., Baekey, D. M., and Morris, K. F. (2008). Pontine respiratory-modulated activity before and after vagotomy in decerebrate cats. J. Physiol. 586, 4265–4282.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

D’Urzo, A. D., Rubinstein, I., Lawson, V. G., Vassal, K. P., Rebuck, A. S., Slutsky, A. S., and Hoffstein, V. (1988). Comparison of glottic areas measured by acoustic reflections vs. computerized tomography. J. Appl. Physiol. 64, 367–370.

Estenne, M., Yernault, J. C., and De Troyer, A. (1985). Rib cage and diaphragm-abdomen compliance in humans: effects of age and posture. J. Appl. Physiol. 59, 1842–1848.

Pubmed Abstract | Pubmed Full Text

Fitz-Clarke, J. R. (2007). Computer simulation of human breath-hold diving: cardiovascular adjustments. Eur. J. Appl. Physiol. Occup. Physiol. 100, 207–224.

CrossRef Full Text

Gaumann, A., Hoeckel, M., and Konerding, M. (1999). The anatomic basis of the transversus and rectus abdominis musculoperitoneal (TRAMP) composite flap. Hernia 3, 39–41.

CrossRef Full Text

Gilroy, R. J., Lavietes, M. H., Loring, S. H., Mangura, B. T., and Mead, J. (1985). Respiratory mechanical effects of abdominal distension. J. Appl. Physiol. 58, 1997–2003.

Pubmed Abstract | Pubmed Full Text

Goldman, M. D., Grassino, A., Mead, J., and Sears, T. A. (1978). Mechanics of the human diaphragm during voluntary contraction: dynamics. J. Appl. Physiol. 44, 840–848.

Pubmed Abstract | Pubmed Full Text

Grassino, A., Goldman, M. D., Mead, J., and Sears, T. A. (1978). Mechanics of the human diaphragm during voluntary contraction: statics. J. Appl. Physiol. 44, 829–839.

Pubmed Abstract | Pubmed Full Text

Grimby, G., Goldman, M., and Mead, J. (1976). Respiratory muscle action inferred from rib cage and abdominal V-P partitioning. J. Appl. Physiol. 41, 739–751.

Pubmed Abstract | Pubmed Full Text

Harris, R. S. (2005). Pressure-volume curves of the respiratory system. Respir. Care 50, 78–98.

Pubmed Abstract | Pubmed Full Text

Hatze, H. (1981). Myocybernetic Control Models of Skeletal Muscle: Characteristics and Applications. Pretoria: University of South Africa.

Hey, E. N., and Price, J. F. (1982). Nasal conductance and effective airway diameter. J. Physiol. (Lond.) 330, 429–437.

Pubmed Abstract | Pubmed Full Text

Hill, A. V. (1938). The heat of shortening and the dynamic constants of muscle. Proc. R. Soc. Lond. B Biol. Sci. 126, 136–195.

CrossRef Full Text

Holm, S. (1979). A simple sequentially rejective multiple test procedure. Scand. J. Statist. 6, 65–70.

Hutchison, A. A., and Lindsey, B. G. (2009). “Modeling central brainstem mechanisms in eupnea, grunting, incremental breathing and vagotomy,” in American Journal of Respiratory and Critical Care Medicine (Abstract) from American Thoracic Society International Conference, San Diego, CA.

Iscoe, S. D. (1998). Control of abdominal muscles. Prog. Neurobiol. 56, 433–506.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Iscoe, S. D., and Gordon, S. P. (1992). Chest wall distortion and discharge of pulmonary slowly adapting receptors. J. Appl. Physiol. 73, 1619–1625.

Pubmed Abstract | Pubmed Full Text

Kamel, K. S., Lau, G., and Stringer, M. D. (2009). In vivo and in vitro morphometry of the human trachea. Clin. Anat. 22, 571–579.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kim, M. J., Druz, W. S., Danon, J., Machnach, W., and Sharp, J. T. (1976). Mechanics of the canine diaphragm. J. Appl. Physiol. 41, 369–382.

Pubmed Abstract | Pubmed Full Text

Knudson, R. J., Mead, J., and Knudson, D. E. (1974). Contribution of airway collapse to supramaximal expiratory flows. J. Appl. Physiol. 36, 653–667.

Pubmed Abstract | Pubmed Full Text

Konno, K., and Mead, J. (1967). Measurement of the separate volume changes of rib cage and abdomen during breathing. J. Appl. Physiol. 22, 407–422.

Pubmed Abstract | Pubmed Full Text

Kreith, F., Berger, S., Churchill, S., Tullis, J. P., Tullis, B., White, F., Kumar, A., Todd, J., Chen, J., Irvine, T., Capobianchi, M., Kennedy, F., Booser, E. R., Wilcock, D., Boehm, R., Reitz, R., Kim, J., McDonald, A., Sherif, S., and Bhushan, B. (2004). “Fluid mechanics,” in The CRC Handbook of Mechanical Engineering, 2nd Edn, eds F. Kreith and D. Y. Goswami (Boca Raton: CRC Press), 3.1–3.232.

Laplace, P. S. (1808). Traité de mécanique céleste: supplément au dixième livre du traité de méchanique céleste sur l’action capillaire. Paris: Duprat.

Lichtenstein, O., Ben-Haim, S. A., Saidel, G. M., and Dinnar, U. (1992). Role of the diaphragm in chest wall mechanics. J. Appl. Physiol. 72, 568–574.

Pubmed Abstract | Pubmed Full Text

Lindsey, B. G., Arata, A., Morris, K. F., Hernandez, Y. M., and Shannon, R. (1998). Medullary raphé neurones and baroreceptor modulation of the respiratory motor pattern in the cat. J. Physiol. 512, 863–882.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lindsey, B. G., Rybak, I. A., and Smith, J. C. (2012). Computational models and emergent properties of respiratory neural networks. Compr. Physiol. 2, 1619–1670.

Loring, S. H., and Mead, J. (1982). Action of the diaphragm on the rib cage inferred from a force-balance analysis. J. Appl. Physiol. 53, 756–760.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Luck, J. C. (1970). Afferent vagal fibres with an expiratory discharge in the rabbit. J. Physiol. (Lond.) 211, 63–71.

Pubmed Abstract | Pubmed Full Text

MacGregor, R. J. (1987). Neural and Brain Modeling. New York: Academic Press.

Macklem, P. T., Gross, D., Grassino, G. A., and Roussos, C. (1978). Partitioning of inspiratory pressure swings between diaphragm and intercostal/accessory muscles. J. Appl. Physiol. 44, 200–208.

Pubmed Abstract | Pubmed Full Text

Mantilla, C. B., and Sieck, G. C. (2011). Phrenic motor unit recruitment during ventilatory and non-ventilatory behaviors. Respir. Physiol. Neurobiol. 179, 57–63.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Matsumoto, S., Ikeda, M., Nishikawa, T., Yoshida, S., Tanimoto, T., Ito, M., Saiki, C., and Takeda, M. (2002). Excitatory mechanism of deflationary slowly adapting pulmonary stretch receptors in the rat lung. J. Pharmacol. Exp. Ther. 300, 597–604.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

McCool, F. D. (2006). Global physiology and pathophysiology of cough. Chest 129, 48S–53S.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mead, J., and Loring, S. H. (1982). Analysis of volume displacement and length changes of the diaphragm during breathing. J. Appl. Physiol. 53, 750–755.

Pubmed Abstract | Pubmed Full Text

Molkov, Y. I., Abdala, A. P. L., Bacak, B. J., Smith, J. C., Paton, J. F. R., and Rybak, I. A. (2010). Late-expiratory activity: emergence and interactions with the respiratory CPG. J. Neurophysiol. 104, 2713–2729.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Nail, B. S., Sterling, G. M., and Widdicombe, J. G. (1972). Patterns of spontaneous and reflexly-induced activity in phrenic and intercostal motoneurons. Exp. Brain Res. 15, 318–332.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ott, M. M., Nuding, S. C., Segers, L. S., O’Connor, R., Morris, K. F., and Lindsey, B. G. (2012). Central chemoreceptor modulation of breathing via multipath tuning in medullary ventrolateral respiratory column circuits. J. Neurophysiol. 107, 603–617.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Paintal, A. S. (1955). Impulses in vagal afferent fibres from specific pulmonary deflation receptors: the response of these receptors to phenyl diguanide, potato starch, 5-hydroxytryptamine and nicotine, and their rôle in respiratory and cardiovascular reflexes. Exp. Physiol. 40, 89–111.

Permutt, S., and Martin, H. B. (1960). Static pressure-volume characteristics of lungs in normal males. J. Appl. Physiol. 15, 819–825.

Pubmed Abstract | Pubmed Full Text

Poliaček, I., Morris, K. F., Lindsey, B. G., Segers, L. S., Rose, M. J., Corrie, L. W.-C., Wang, C., Pitts, T. E., Davenport, P. W., and Bolser, D. C. (2011). Blood pressure changes alter tracheobronchial cough: computational model of the respiratory-cough network and in vivo experiments in anesthetized cats. J. Appl. Physiol. 111, 861–873.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ratnovsky, A., Elad, D., and Halpern, P. (2008). Mechanics of respiratory muscles. Respir. Physiol. Neurobiol. 163, 82–89.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ratnovsky, A., Zaretsky, U., Shiner, R. J., and Elad, D. (2003). Integrated approach for in vivo evaluation of respiratory muscles mechanics. J. Biomech. 36, 1771–1784.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Reid, M. B., Feldman, H. A., and Miller, M. J. (1987). Isometric contractile properties of diaphragm strips from alcoholic rats. J. Appl. Physiol. 63, 1156–1164.

Pubmed Abstract | Pubmed Full Text

Renotte, C., Remy, M., and Saucez, P. (1998). Dynamic model of airway pressure drop. Med. Biol. Eng. Comput. 36, 101–106.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Riddle, W., and Younes, M. (1981). A model for the relation between respiratory neural and mechanical outputs. II. Methods. J. Appl. Physiol. 51, 979–989.

Pubmed Abstract | Pubmed Full Text

Roca, J., Burgos, F., Sunyer, J., Saez, M., Chinn, S., Anto, J. M., Rodriguez-Roisin, R., Quanjer, P. H., Nowak, D., and Burney, P. (1998). References values for forced spirometry. Group of the European community respiratory health survey. Eur. Respir. J. 11, 1354–1362.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rohrer, F. (1915). Der stromungswiderstand in den menschlichen atemwegen und der einfluss der unregelmassigen verzweigung des bronchialsystems auf den atmungsverlauf in verschiedenen lungenbezirken. Archiv fur die Gesamte Physiologie des Menschen und der Tiere 162, 225–299.

Rosen, J., Fuchs, M. B., and Arcan, M. (1999). Performances of Hill-type and neural network muscle models: toward a myosignal-based exoskeleton. Comput. Biomed. Res. 32, 415–439.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rybak, I. A., O’Connor, R., Ross, A., Shevtsova, N. A., Nuding, S. C., Segers, L. S., Shannon, R., Dick, T. E., Dunin-Barkowski, W. L., Orem, J. M., Solomon, I. C., Morris, K. F., and Lindsey, B. G. (2008). Reconfiguration of the pontomedullary respiratory network: a computational modeling study with coordinated in vivo experiments. J. Neurophysiol. 100, 1770–1799.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Segers, L. S., Nuding, S. C., Dick, T. E., Shannon, R., Baekey, D. M., Solomon, I. C., Morris, K. F., and Lindsey, B. G. (2008). Functional connectivity in the pontomedullary respiratory network. J. Neurophysiol. 100, 1749–1769.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Segers, L. S., Nuding, S. C., Vovk, A., Baekey, D. M., O’Connor, R., Morris, K. F., Lindsey, B. G., Shannon, R., and Bolser, D. C. (2012). Discharge identity of medullary inspiratory neurons is altered during repetitive fictive cough. Front. Physiol. 3:223. doi:10.3389/fphys.2012.00223

CrossRef Full Text

Shannon, R. (1986). “Reflexes from the respiratory muscles and controvertebral joints,” in Handbook of Physiology. The Respiratory System. Control of Breathing, Sect. 3, Part 1, eds N. S. Cherniack and J. G. Widdicombe (Washington, DC: American Physiological Society), 431–447.

Shannon, R., Baekey, D. M., Morris, K. F., Li, Z., and Lindsey, B. G. (2000). Functional connectivity among ventrolateral medullary respiratory neurones and responses during fictive cough in the cat. J. Physiol. 525, 207–224.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Shannon, R., Baekey, D. M., Morris, K. F., and Lindsey, B. G. (1998). Ventrolateral medullary respiratory network and a model of cough motor pattern generation. J. Appl. Physiol. 84, 2020–2035.

Pubmed Abstract | Pubmed Full Text

Shannon, R., Baekey, D. M., Morris, K. F., Nuding, S. C., Segers, L. S., and Lindsey, B. G. (2004). Production of reflex cough by brainstem respiratory networks. Pulm. Pharmacol. Ther. 17, 369–376.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sharpey-Schafer, E. P. (1953). Effects of coughing on intra-thoracic pressure, arterial pressure and peripheral blood flow. J. Physiol. (Lond.) 122, 351–357.

Pubmed Abstract | Pubmed Full Text

Siafakas, N. M., Morris, A. J., and Green, M. (1979). Thoracoabdominal mechanics during relaxed and forced vital capacity. J. Appl. Physiol. 47, 38–42.

Pubmed Abstract | Pubmed Full Text

Simpson, L. L. (1968). Sizing piping for process plants. Chem. Eng. 75, 192–214.

Smith, J., and Bellemare, F. (1987). Effect of lung volume on in vivo contraction characteristics of human diaphragm. J. Appl. Physiol. 62, 1893–1900.

Pubmed Abstract | Pubmed Full Text

Smith, J. A., Aliverti, A., Quaranta, M., McGuinness, K., Kelsall, A., Earis, J., and Calverley, P. M. (2012). Chest wall dynamics during voluntary and induced cough in healthy volunteers. J. Physiol. (Lond.) 590, 563–574.

Pubmed Abstract | Pubmed Full Text

Song, C., Alijani, A., Frank, T., Hanna, G. B., and Cuschieri, A. (2006). Mechanical properties of the human abdominal wall measured in vivo during insufflation for laparoscopic surgery. Surg. Endosc. 20, 987–990.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Suárez, A. A., Pessolano, F. A., Monteiro, S. G., Ferreyra, G., Capria, M. E., Mesa, L., Dubrovsky, A., and De Vito, E. L. (2002). Peak flow and peak cough flow in the evaluation of expiratory muscle weakness and bulbar impairment in patients with neuromuscular disease. Am. J. Phys. Med. Rehabil. 81, 506–511.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tully, A., Brancatisano, A., Loring, S. H., and Engel, L. A. (1990). Relationship between thyroarytenoid activity and laryngeal resistance. J. Appl. Physiol. 68, 1988–1996.

Pubmed Abstract | Pubmed Full Text

Tully, A., Brancatisano, A., Loring, S. H., and Engel, L. A. (1991). Influence of posterior cricoarytenoid muscle activity on pressure-flow relationship of the larynx. J. Appl. Physiol. 70, 2252–2258.

Pubmed Abstract | Pubmed Full Text

Wei, J. Y., and Shen, E. (1985). Vagal expiratory afferent discharges during spontaneous breathing. Brain Res. 335, 213–219.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Younes, M., and Riddle, W. (1981). A model for the relation between respiratory neural and mechanical outputs. I. Theory. J. Appl. Physiol. 51, 963–978.

Younes, M., Riddle, W., and Polacheck, J. (1981). A model for the relation between respiratory neural and mechanical outputs. III Validation . J. Appl. Physiol. 51, 990–1001.

Yu, J., Wang, Y. F., and Zhang, J. W. (2003). Structure of slowly adapting pulmonary stretch receptors in the lung periphery. J. Appl. Physiol. 95, 385–393.

Pubmed Abstract | Pubmed Full Text