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., , 2000Baekey 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., 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 . 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 . 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  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;, and with separate neural control of the diaphragm and abdominal wall muscles www.frontiersin.org 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 .
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 Shannon et al., 2000;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 .

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. (1990Tully et al. ( , 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 (V di andV ab ), and is defined by the subsequent equations. The parameters referenced in the model equations are listed in Table 1.

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 Frontiers in Physiology | Computational Physiology and Medicine  Goldman et al. (1978) and Chow and Darling (1999) Konno and Mead (1967), Figure 14 V kmFRC ab Volume of the abdominal wall at FRC as a fraction of VC relative to RV 0.0400 Dimensionless Konno and Mead (1967), Figure 14 V kmTLC ab Volume of the abdominal wall at TLC as a fraction of VC relative to RV 0.3391 Dimensionless Konno and Mead (1967), Figure 14 VC Vital capacity 5.370 L Roca et al. (1998) Konno and Mead (1967), Figure 14 f di Ratio of diaphragm length at TLC to RV 0.65 Dimensionless Smith and Bellemare (1987) C rc Compliance of the rib cage 0.110 L/cmH 2 O Derived from Gilroy et al. (1985) (Continued) www.frontiersin.org 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.

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.
Frontiers in Physiology | Computational Physiology and Medicine

Pleural pressure
In Eq. 3, tracheal flow (V L ) is the derivative of lung volume. R rs is the resistance of the airway, which when multiplied byV L gives the pressure drop, by the hydraulic analog of Ohm's law. σ L is the recoil pressure of the lung.

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 1 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 1 was determined by fitting our model to published data   Figure 4) giving diaphragm pressure as a function of rib cage and abdominal volume at a fixed diaphragm activation. Given a value of C 1 , our model equations will calculate the diaphragm pressure from the volumes. We then found the value of C 1 that was the best fit of the calculated pressure values to the measured values. The resulting equation is V di + C 1 V rc + V ab = V sum , which together with Eq. 20, gives us Eq. 4 for lung volume (V L ).

Airway resistance
Rohrer's equation was used to calculate airway resistance (Rohrer, 1915;Hey and Price, 1982). The equation Pressure = K 1 · Flow + K 2 · Flow 2 or, dividing through by flow, Resistance = K 1 + K 2 |Flow| was applied twice, once for laryngeal resistance (k 1 + k 2 |V L |) and once for the resistance of the oropharynx and lower airway (0.72 + 0.44|V L |), to give Eq. 5. We used a value of K 2 = 0.44 for the oropharynx and lower airway based on the assumption that K 2 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 1 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 1 and k 2 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(Tully et al., , 1991. The coefficient k 2 given by Eq. 7 is calculated using the equation for flow through an orifice (Simpson, 1968, Table II). The value of k 2 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 2 , and got 0.167. This approach resulted in a resting value for k 2 of 0.681.
We calculated the ratio of the mean resting value of k 1 to the mean resting value of k 2 (Tully et al., 1990, Table 2), and multiplied the ratio by 0.681 to get a resting value of k 1 of 0.0035 for the larynx, which is small relative to K 1 for the rest of the airway, justifying our assumption above that K 1 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 4 (Kreith et al., 2004), which we use in Eq. 6 for k 1 . Plugging in the resting value of k 1 determined above and the resting value of d, we solved for the coefficient, which gives us 49.6.

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.

Abdominal fraction of the rib cage
At total lung capacity (TLC), none of the diaphragm is apposed to the rib cage , so we assumed that at all www.frontiersin.org lung volumes, a portion of V di equal to V di at TLC (V TLC di ) does not contribute to the zone of apposition. The remainder (V di − V TLC di ) is divided by the remainder plus the lung volume (V di −V TLC di +V L ) 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 , is always exposed to abdominal pressure, and is represented by f TLC a 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 1 ), so we divided the previously calculated fraction of the smaller rib cage by (1 + C 1 ) to turn it into a fraction of the larger rib cage before adding it to f TLC a .

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

Recoil pressure of the diaphragm
In Eq. 12, the term u di σ max di F di fl F di fv corresponds to the Hill muscle model (Ratnovsky et al., 2003, Eq. A6), except that σ max di is a pressure rather than a force and F di fl and F di fv 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 .
In Eq. 12, u di is phrenic activation of the diaphragm; σ max di is static diaphragm recoil pressure at optimum length and maximum activation; R diVdi is the pressure due to the passive resistance of the diaphragm.
Passive recoil pressure of the diaphragm In Eq. 13, σ psv di is the passive transdiaphragmatic pressure as a function of diaphragm volume. This pressure is taken to be zero at resting diaphragm volume V FRC di (Agostoni and Rahn, 1960) and below, and quadratic above (Reid et al., 1987).

Volume-pressure relationship of the diaphragm
In Eq. 14, F di fl 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 min di 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).

Pressure-flow relationship of the diaphragm
In Eq. 15, F di fv 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 vari-ableV di is the time derivative of the volume under the diaphragm. The maximum rate of change of diaphragm volume,V max di , was derived from data which gives transdiaphragmatic pressure as a function of flow at several levels of diaphragm activation up to 45%   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 , 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 forV max di .

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 di ica is the pressure due to the action of the inspiratory intercostals, which are assumed to be inactive when the diaphragm volume is above V FRC di (low lung volumes). Below V FRC di , 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 FRC di to its maximum value of P diTLC ica at V TLC di . The parameter P diTLC ica was calculated as the pressure necessary to complete the pressure balance on the rib cage at TLC.
In Eq. 18, P ab ica 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 abRV ica at RV to P abTLC ica at TLC. In Eq. 18, P abTLC ica 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 abRV ica 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.

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 (V max rc ), and with decreasing pressure it asymptotically approaches a minimum (V min rc ). 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 (V rc ) The parameter σ mul rc 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 1 ). The factor of (1 + C 1 ) 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).

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

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 ab F CEmax F ab fl F ab fv 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 2 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 toV FRC ab , 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 (V ab ). (21) In Eq. 22, F ab fl 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.
In Eq. 23, F ab fv is the force-velocity relationship of the abdominal wall muscles (Hatze, 1981;Rosen et al., 1999;Artemiadis and Kyriakopoulos, 2005). The variableL 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,

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.
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.
In Eq. 28, the Pythagorean Theorem is applied in a transverse plane to get h 0 , the distance from the peak of the abdominal wall to the frontal plane that serves as a boundary of the abdominal wall volume.
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.
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.
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).
Equations for calculated parameters Frontiers in Physiology | Computational Physiology and Medicine

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

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 1 ), of the two populations by the expression (0.3P + 0.7P 1 )/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 −7 ; 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,    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.   including lung volume, tracheal flow, and alveolar pressure. The three phases of the cough cycle (Bolser et al., 2003) are highlighted.

Frontiers in Physiology | Computational Physiology and Medicine
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 Frontiers in Physiology | Computational Physiology and Medicine     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). 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.
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 -C 1 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 -C 2 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).
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 1 ) 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 1 ) 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 1 ). 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.   ) www.frontiersin.org

FIGURE 2 | Continued
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 2 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 2 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   for (A,B). See text for further details.
(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 2 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 2 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. 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).

Comparisons with behaviors of antecedent 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.

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. (1990Tully et al. ( , 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;Lindsey et al., 2012). 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 . 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.
www.frontiersin.org The behaviors of the networks in Rybak et al. (2008) and Poliaček et al. (2011) (Holm, 1979); significant values (at the 0.05 level) are marked with an asterisk.

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 flowvolume 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  ., 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.
www.frontiersin.org 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 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 , 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 . 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.