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 (udi, uab, and ulm), and the diaphragm and abdominal wall volumes (Vdi and Vab) and their time derivatives (di and ab), and is defined by the subsequent equations. The parameters referenced in the model equations are listed in Table 1.

Table 1

Parameters used in the biomechanical model.

Parameter	Definition	Value	Units	Source
FREE PARAMETERS USED IN THE MODEL
C1	Rib cage contribution to abdominal volume	0.369	Dimensionless	Derived from Grassino et al. (1978)
Cab	Compliance of the abdominal wall	0.108	L/cmH2O	Derived from Estenne et al. (1985)
CL	Compliance of the lung	0.201	L/cmH2O	Derived from Permutt and Martin (1960)
ct	Transverse chord of the abdominal wall	0.32	m	Derived from Song et al. (2006)
D	Diameter of the trachea	18	mm	Derived from Baier et al. (1977) and Kamel et al. (2009)
faTLC	Obligatory ring fraction	0.15	Dimensionless	Mead and Loring (1982)
FCEmax	Maximal force capacity of the external oblique	33	N	Ratnovsky et al. (2003)
Fdi	Fraction of the diaphragm pressure expanding the rib cage via insertional forces	0.15	Dimensionless	Derived from Loring and Mead (1982)
k	Conversion factor from force to surface tension in abdominal muscle	0.68	m cmH2O/N	Derived from De Troyer et al. (1990) and Ratnovsky et al. (2003)
LCEO	Length of human transversus abdominis	19.1	cm	Gaumann et al. (1999)
PicaabTLC	Maximal expiratory pressure due to intercostal and accessory muscles at total lung capacity (TLC)	−135	cmH2O	Derived from Ratnovsky et al. (2008)
Rab	Passive resistance of the abdominal wall	1.5	cmH2O/(L/s)	Derived from Barnas et al. (1989)
Rdi	Passive resistance of the diaphragm	6	cmH2O/(L/s)	Derived from Barnas et al. (1989)
Rrc	Passive resistance of the rib cage	2.7	cmH2O/(L/s)	Derived from Barnas et al. (1989)
Vc	Mediastinal plus lung blood and tissue volume	1.756	L	Derived from Cluzel et al. (2000)
VCEmax	Maximal contractile velocity of the external oblique	34.7	cm/s	Ratnovsky et al. (2003)
VdiFRC	Volume under diaphragm at functional residual capacity (FRC)	2.967	L	Derived from Cluzel et al. (2000)
V˙dimax	Maximal rate of change of volume under diaphragm	2.449	L/s	Derived from Goldman et al. (1978) and Chow and Darling (1999)
FREE PARAMETERS USED IN CALCULATED PARAMETERS
VLFRC	Volume of the lung at FRC	2.29	L	Cluzel et al. (2000)
VrckmFRC	Volume of the rib cage at FRC as a fraction of VC relative to residual volume (RV)	0.1282	Dimensionless	Konno and Mead (1967), Figure 14
VabkmFRC	Volume of the abdominal wall at FRC as a fraction of VC relative to RV	0.0400	Dimensionless	Konno and Mead (1967), Figure 14
VabkmTLC	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)
σdiRV	Passive recoil pressure of the diaphragm at RV	20	cmH2O	Derived from Agostoni et al. (1966), Grassino et al. (1978) and Siafakas et al. (1979)
VrcFRC	Volume of the rib cage at FRC	7.013	L	Derived from Cluzel et al. (2000)
VrckmTLC	Volume of the rib cage at TLC as a fraction of VC relative to RV	0.6609	Dimensionless	Konno and Mead (1967), Figure 14
fdi	Ratio of diaphragm length at TLC to RV	0.65	Dimensionless	Smith and Bellemare (1987)
Crc	Compliance of the rib cage	0.110	L/cmH2O	Derived from Gilroy et al. (1985)
CALCULATED PARAMETERS USED IN CALCULATED PARAMETERS
VabRV	Volume behind abdominal wall at RV		L	Eqs 32 and 33
VLRV	Volume of the lung at RV		L	Eq. 34
σLRV	Passive recoil pressure of the lung at RV		cmH2O	Eq. 35
σabRV	Passive recoil pressure of the abdominal wall at RV		cmH2O	Eq. 36
VabFRC	Volume behind abdominal wall at FRC		L	Eq. 37
VdiRV	Volume under diaphragm at RV		L	Eq. 38
VabTLC	Volume behind abdominal wall at TLC		L	Eq. 39
VLTLC	Volume of the lung at TLC		L	Eq. 40
σLTLC	Passive recoil pressure of the lung at TLC		cmH2O	Eq. 41
σabTLC	Passive recoil pressure of the abdominal wall at TLC		cmH2O	Eq. 42
σdiTLC	Passive recoil pressure of the diaphragm at TLC		cmH2O	Eq. 43
σrcTLC	Passive recoil pressure of the rib cage at TLC		cmH2O	Eq. 44
σrcRV	Passive recoil pressure of the rib cage at RV		cmH2O	Eq. 45
faRV	Fraction of the rib cage exposed to abdominal pressure at RV		Dimensionless	Eq. 46
CALCULATED PARAMETERS USED IN THE MODEL
Vsum	Sum of diaphragm, rib cage, and abdominal wall volume contributions		L	Eq. 47
Vab0	Volume behind abdominal wall at zero passive tension		L	Eqs 32 and 33
VL0	Volume of the lung at zero passive tension		L	Eq. 48
Kdipsv	Coefficient of passive diaphragm recoil pressure		cmH2O/L2	Eq. 49
Ldimin	Ratio of diaphragm length at zero volume to resting length		Dimensionless	Eq. 50
PicaabRV	Maximal expiratory pressure due to intercostal and accessory muscles at TLC		cmH2O	Eq. 51
PicadiTLC	Maximal inspiratory pressure due to intercostal and accessory muscles at TLC		cmH2O	Eq. 52
σdimax	Maximum active recoil pressure of the diaphragm		cmH2O	Eq. 53
σrcadd	Passive rib cage recoil pressure midway between volume limits		cmH2O	Eq. 54
σrcmul	Rib cage sigmoid compliance coefficient		cmH2O	Eq. 55
Vdi0	Volume under diaphragm at zero tension		L	Eq. 56
VdiTLC	Volume under diaphragm at TLC		L	Eq. 57
Vrcmax	Upper limit of rib cage volume		L	Eq. 58
Vrcmin	Lower limit of rib cage volume		L	Eq. 59
VrcRV	Volume of the rib cage at RV		L	Eq. 60
VrcTLC	Volume of the rib cage at TLC		L	Eq. 61
Vrc0	Volume of the rib cage at zero tension		L	Eq. 62

Pressure balance on the rib cage

In Eq. 1, Ppl is the pleural pressure seen by the interior surface of the rib cage, Pab−pl corrects for the fact that the rib cage sees abdominal pressure in the zone of apposition, Fdiσdi is the equivalent pressure due to the insertional forces of the diaphragm on the lower ribs, Pica 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.

Pleural pressure

In Eq. 3, tracheal flow (L) is the derivative of lung volume. Rrs is the resistance of the airway, which when multiplied by 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 Vdi and Vab. Diaphragm volume, Vdi, is the volume above the level of the diaphragm insertions on the rib cage and below the dome of the diaphragm. Abdominal wall volume, Vab, 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 Vdi + Vab = Vsum (Vsum constant). Under this assumption, the abdominal volume completely determines the volume under the diaphragm, Vdi. Following Lichtenstein et al. (1992), Vdi 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, C1Vrc, 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).

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.

The value of C1 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 C1, our model equations will calculate the diaphragm pressure from the volumes. We then found the value of C1 that was the best fit of the calculated pressure values to the measured values. The resulting equation is Vdi + C1Vrc + Vab = Vsum, which together with Eq. 20, gives us Eq. 4 for lung volume (VL).

Airway resistance

Rohrer’s equation was used to calculate airway resistance (Rohrer, 1915; Hey and Price, 1982). The equation Pressure = K1 · Flow + K2 · Flow2 or, dividing through by flow, Resistance = K1 + K2 |Flow| was applied twice, once for laryngeal resistance (k1 + k|L|) and once for the resistance of the oropharynx and lower airway (0.72 + 0.44|L|), to give Eq. 5. We used a value of K2 = 0.44 for the oropharynx and lower airway based on the assumption that K2 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 K1 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.

Table 2

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

Population name	Size	Resting threshold (mV)	THO variability (mV)	Membrane time constant	Post-spike increase in GK+	Post-spike GK+ time constant (ms)	Adaptation threshold increase	Adaptation (ms)	Noise amplitude	DC (mV)
	N	THO		TMEM	B	TGK	C	TTH
I-pons	100	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.03	2.0
rostral IE-pons	100	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.3	5.0
caudal IE-pons	100	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.3	5.0
E-pons	100	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.3	13.0
EI-pons	100	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.3	20.0
NRM-pons	100	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.03	25.0
NRM-BötC	300	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.03	25.0
ΔI-DRIVER	300	See Table A2 in Rybak et al. (2008) for I-Driver neuron properties
ΔI-Dec	300	10.0	1.0	6.0	75.0	8.5	0.9	1500.0	0.2	20.0
+I-Dec_2	300	10.0	1.0	6.0	75.0	8.5	0.8	1200.0	0.2	5.0
ΔI-Aug	300	12.0	3.0	6.0	75.0	5.0	0.0	5000.0	0.6	18.0
ΔVRC IE	99	18.0	2.0	9.0	50.0	5.0	0.0	1000.0	0.075	0.0
ΔE-Dec-Phasic	300	9.0	1.0	6.0	75.0	8.5	0.9	1500.0	0.1	4.0
ΔE-Dec-Tonic	300	8.0	1.0	9.0	50.0	3.8	0.8	2000.0	0.3	0.0
E-Aug-early	300	10.0	1.0	6.0	27.0	2.5	0.0	500.0	0.3	30.0
E-Aug-late	300	10.0	1.0	9.0	27.0	2.5	0.0	500.0	0.1	27.0
ΔE-Aug-Cough	300	12.0	2.0	9.0	75.0	7.5	0.0	500.0	0.2	0.0
+E-Aug (+)	300	10.0	2.0	6.0	75.0	7.5	0.1	1000.0	0.5	20.0
ΔE-Aug-BS (+)	300	12.0	3.0	6.0	100.0	6.0	0.1	1000.0	0.5	0.0
ΔPump (+)	300	5.0	0.5	6.0	25.0	3.8	0.08	500.0	0.1	0.0
ΔPump (-)	300	5.0	0.0	6.0	25.0	3.8	0.08	500.0	0.1	0.0
ΔI-Aug-BS	300	12.0	3.0	6.0	75.0	5.0	0.0	5000.0	0.5	0.0
+Phrenic	210	10.0	2.0	5.0	200.0	6.0	0.08	500.0	0.5	0.0
+Phrenic-HT	70	16.0	2.0	60.0	200.0	5.0	0.08	500.0	0.5	0.0
ΔLumbar	210	15.0	2.0	6.0	75.0	7.5	0.08	500.0	0.5	0.0
+Lumbar-HT	70	18.0	2.0	30.0	200.0	7.5	0.1	500.0	0.5	0.0
ΔILM	300	20.0	1.0	6.0	25.0	3.8	0.08	500.0	0.1	2.0
ΔE-Dec-pre-ELM	300	11.0	0.0	6.0	100.0	6.0	0.8	500.0	0.5	1.0
ΔELM	300	18.0	2.0	6.0	100.0	6.0	0.9	100.0	0.5	0.0
ΔLUNG PSRs	300	11.0	1.0	9.0	20.0	7.0	0.0	500.0	0.5	0.0
+Lung deflation receptors (Def_1)	300	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.5	0.0
+Lung distortion receptors (Dis_1)	300	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.5	0.0
>Cough 2nd order	100	10.0	1.0	9.0	20.0	7.0	0.3	500.0	0.1	0.0
ΔDeflation 2nd order	300	8.0	1.0	9.0	27.0	2.5	0.5	1000.0	0.3	0.0
Raphé 8	100	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.01	0.0
Raphé 28	100	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.1	0.0
Raphé 29	100	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.5	0.0
Raphé 30	100	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.5	0.0
Raphé 31	100	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.1	10.0
Raphé 32	100	10.0	1.0	9.0	20.0	7.0	0.0	500.0	0.1	10.0

Variable names used by MacGregor (.

The values of k1 and k2 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 ulm 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 ulm = 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 ulm (Tully et al., 1990, 1991).

The coefficient k2 given by Eq. 7 is calculated using the equation for flow through an orifice (Simpson, 1968, Table II). The value of k2 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 k2, and got 0.167. This approach resulted in a resting value for k2 of 0.681.

We calculated the ratio of the mean resting value of k1 to the mean resting value of k2 (Tully et al., 1990, Table 2), and multiplied the ratio by 0.681 to get a resting value of k1 of 0.0035 for the larynx, which is small relative to K1 for the rest of the airway, justifying our assumption above that K1 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/d4 (Kreith et al., 2004), which we use in Eq. 6 for k1. Plugging in the resting value of k1 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 fa 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 Pab − 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 (Mead and Loring, 1982), so we assumed that at all lung volumes, a portion of Vdi equal to Vdi at TLC does not contribute to the zone of apposition. The remainder is divided by the remainder plus the lung volume 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 in Eq. 10. Our rib cage volume, Vrc, 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 + C1), so we divided the previously calculated fraction of the smaller rib cage by (1 + C1) to turn it into a fraction of the larger rib cage before adding it to

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 . Scales on right are for (A,B). See text for further details.

Recoil pressure of the lung

Equation 11 assumes a linear relationship between lung volume and recoil pressure. CL is lung compliance. VL0 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 corresponds to the Hill muscle model (Ratnovsky et al., 2003, Eq. A6), except that is a pressure rather than a force and and 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, udi is phrenic activation of the diaphragm; is static diaphragm recoil pressure at optimum length and maximum activation; Rdidi is the pressure due to the passive resistance of the diaphragm.

Passive recoil pressure of the diaphragm

In Eq. 13, is the passive transdiaphragmatic pressure as a function of diaphragm volume. This pressure is taken to be zero at resting diaphragm volume (Agostoni and Rahn, 1960) and below, and quadratic above (Reid et al., 1987).

Volume-pressure relationship of the diaphragm

In Eq. 14, 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 Vdi0 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 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, 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 di is the time derivative of the volume under the diaphragm. The maximum rate of change of diaphragm volume, 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

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 . Pairs of symbols connected by a line indicate no significant difference.

Pressure of the intercostal and accessory muscles

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

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

In Eq. 18, is the pressure due to the action of the expiratory intercostals, which is assumed to be proportional to the abdominal muscle activation (uab); the proportionality constant itself is assumed to scale linearly with the rib cage volume, changing from at RV to at TLC. In Eq. 18, 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. is the pressure due to the expiratory intercostals at RV and maximal abdominal activation. We calculated this parameter by solving the model equations for Pica 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 and with decreasing pressure it asymptotically approaches a minimum A generalized logistic function is used for the sigmoid, giving Vrc 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 (Rrc) and the rate of change of its volume (rc) The parameter 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 Crc/(1 + C1). The factor of (1 + C1) appears because Crc is for a rib cage volume defined differently than Vrc. Crc 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 (Vrc) is the sum of the lung volume (VL), the volume under the diaphragm (Vdi), and the volume of the mediastinum and the lung blood and tissue (Vc).

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, Vab, is bounded by this surface and by a frontal plane. The values for the sagittal (rs) and transverse (rt) 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: rs = 8.00 rt − 1.10. The length of the longest transverse chord in the bounding frontal plane (ct) was found which gave the volume change stated in the paper.

In Eq. 21, is the Hill muscle model equation (Ratnovsky et al., 2003, Eq. A6); uab is the activation of the diaphragm by the lumbar motor neurons; FCEmax is the maximal force capacity of the contractile element for a 1.5 cm2 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/rs + 1/rt) converts the surface tension to a pressure, using the Law of Laplace (Laplace, 1808). The second term on the right, (Vab − Vab0)/Cab, is the passive recoil pressure of the abdominal wall. Vab0 is the volume behind the abdominal wall at which the recoil pressure is 0. This was taken to be equal to since we assume a supine position. Cab is the compliance of the abdominal wall. The final term is the pressure due to the passive resistance of the abdominal wall (Rab) and the rate of change of abdominal volume (ab).

In Eq. 22, is the static force-length relationship of the abdominal wall (Ratnovsky et al., 2003, Eq. A7); LCE is the length of the transversus abdominis; LCE0 is its resting length.

In Eq. 23, is the force-velocity relationship of the abdominal wall muscles (Hatze, 1981; Rosen et al., 1999; Artemiadis and Kyriakopoulos, 2005). The variable is the velocity of the contractile element (the time derivative of LCE) and the parameter VCEmax is the maximal contractile velocity of canine external oblique muscle (Ratnovsky et al., 2003, Table 1).

Length-volume relationship of the abdominal wall

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

In Eq. 27, the Pythagorean Theorem is applied in the midsagittal plane to get cs, 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 h0, 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 (rs) and the transverse radius (rt) derived from the results in Song et al. (2006).

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

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 4

Connectivity for the network model modified from Poliaček et al. (.

Source population	Target population	Synaptic type	Conduction times		No. of terminals	Synaptic strength	Source pop. N	Target pop. N	Divergence	Mean no. of terminals	Convergence
			Min	Max
I-Driver	I-Dec	ex_1	2	6	100	0.006	300	300	84.99 ± 3.14	1.18	84.99 ± 7.54
I-Driver	I-Aug	ex_1	2	6	100	0.01	300	300	84.93 ± 3.00	1.18	84.93 ± 7.73
I-Driver	I-Driver	ex_1	0	4	50	0.003	300	300	46.34 ± 1.76	1.08	46.34 ± 5.84
E-Dec-Phasic	I-Driver	inh_22	2	6	50	0.03	300	300	46.16 ± 1.77	1.08	46.16 ± 9.24
E-Dec-Phasic	E-Aug-early	inh_4	0	2	150	0.012	300	300	118.24 ± 4.01	1.27	118.24 ± 8.86
E-Dec-Phasic	E-Aug-late	inh_4	2	4	150	0.04	300	300	118.10 ± 3.94	1.27	118.10 ± 9.61
E-Dec-Phasic	VRC-IE	inh_4	0	2	50	0.1	300	99	39.33 ± 2.31	1.27	119.18 ± 7.22
E-Dec-Phasic	I-Dec	inh_4	0	2	200	0.2	300	300	146.07 ± 4.42	1.37	146.07 ± 8.71
+E-Dec-Phasic	I-Aug-BS	inh_4	2	6	100	0.15	300	300	85.21 ± 3.08	1.17	85.21 ± 7.27
E-Dec-Phasic	I-Aug	inh_4	2	6	50	0.1	300	300	46.19 ± 1.72	1.08	46.19 ± 5.43
E-Dec-Phasic	rostral IE-pons	ex_13	2	4	100	0.001	300	100	63.13 ± 2.98	1.58	189.39 ± 6.74
E-Dec-Phasic	caudal IE-pons	ex_13	2	4	100	0.001	300	100	63.23 ± 3.19	1.58	189.70 ± 9.56
+E-Dec-Phasic	I-Dec_2	inh_4	2	6	150	0.1	300	300	118.39 ± 4.42	1.27	118.39 ± 8.69
+E-Dec-Phasic	E-Aug (+)	inh_4	2	6	100	0.03	300	300	85.18 ± 2.99	1.17	85.18 ± 10.46
+E-Dec-Phasic	E-Aug-Cough (−)	inh_22	2	6	100	0.025	300	300	85.20 ± 3.05	1.17	85.20 ± 9.59
I-Dec	E-Aug-early	inh_7	2	6	115	0.5	300	300	95.67 ± 3.39	1.20	95.67 ± 7.14
I-Dec	E-Dec-Phasic	inh_25	2	6	200	0.2	300	300	145.80 ± 5.20	1.37	145.80 ± 9.31
I-Dec	I-Aug	inh_25	2	6	120	0.025	300	300	99.27 ± 3.48	1.21	99.27 ± 8.70
I-Dec	E-Aug-late	inh_7	2	6	115	0.5	300	300	95.71 ± 3.51	1.20	95.71 ± 7.84
I-Dec	VRC-IE	inh_25	0	4	33	0.025	300	99	28.15 ± 1.85	1.17	85.29 ± 6.01
I-Dec	E-Dec-Tonic	inh_25	2	6	100	0.001	300	300	84.90 ± 2.96	1.18	84.90 ± 9.43
I-Dec	ILM	ex_1	0	3	50	0.002	300	300	46.35 ± 1.71	1.08	46.35 ± 6.00
I-Dec	E-Aug-BS (+)	inh_7	0	4	100	0.05	300	300	85.44 ± 3.02	1.17	85.44 ± 9.96
I-Dec	EI-pons	ex_13	2	4	100	0.001	300	100	63.33 ± 3.16	1.58	190.00 ± 8.78
I-Dec	I-pons	ex_13	2	4	100	0.0005	300	100	63.50 ± 3.23	1.57	190.51 ± 8.03
I-Dec	I-pons	inh_7	2	4	100	0.0005	300	100	63.66 ± 3.21	1.57	190.98 ± 8.42
I-Dec	rostral IE-pons	inh_4	2	4	100	0.0001	300	100	63.55 ± 3.04	1.57	190.65 ± 8.05
I-Dec	caudal IE-pons	inh_7	2	4	100	0.0001	300	100	63.27 ± 3.06	1.58	189.80 ± 8.43
I-Dec	Lumbar	inh_7	0	4	100	0.1	300	210	79.70 ± 3.14	1.25	113.86 ± 8.97
I-Dec	E-Dec-pre-ELM	inh_4	2	6	200	0.06	300	300	146.12 ± 4.70	1.37	146.12 ± 8.41
+I-Dec	E-Aug (+)	inh_4	0	2	130	1.0	300	300	105.84 ± 3.65	1.23	105.84 ± 7.43
+I-Dec	ELM	inh_4	0	5	200	0.06	300	300	146.22 ± 4.96	1.37	146.22 ± 9.62
+I-Dec	Lumbar-HT	inh_7	2	6	100	0.1	300	70	53.32 ± 2.84	1.88	228.53 ± 6.63
+I-Dec	I-Dec	inh_25	2	6	140	0.0125	300	300	112.11 ± 3.97	1.25	112.11 ± 9.09
+I-Dec	E-Aug-Cough (−)	inh_7	2	6	115	0.16	300	300	95.71 ± 3.51	1.20	95.71 ± 7.84
I-Aug	I-Aug	ex_1	0	5	50	0.025	300	300	45.99 ± 1.77	1.09	45.99 ± 5.31
I-Aug	caudal IE-pons	inh_7	2	4	100	0.0001	300	100	63.43 ± 3.19	1.58	190.28 ± 9.05
I-Aug	I-Aug-BS	ex_1	2	6	100	0.06	300	300	85.20 ± 3.01	1.17	85.20 ± 7.27
I-Aug	ILM	ex_1	2	6	70	0.035	300	300	62.72 ± 2.32	1.12	62.72 ± 6.93
I-Aug	VRC-IE	ex_1	2	6	16	0.0015	300	99	14.82 ± 0.97	1.08	44.92 ± 6.80
I-Aug	I-pons	ex_13	2	4	100	0.0025	300	100	63.59 ± 2.94	1.57	190.76 ± 7.71
I-Aug	rostral IE-pons	inh_4	2	4	100	0.0001	300	100	63.43 ± 2.95	1.58	190.29 ± 8.57
+I-Aug	I-Dec_2	ex_28	2	6	200	0.01	300	300	146.05 ± 4.72	1.37	146.05 ± 9.65
E-Aug-early	E-Dec-Phasic	inh_4	2	6	110	0.007	300	300	92.32 ± 3.45	1.19	92.32 ± 7.95
E-Aug-early	I-Dec	inh_4	0	5	100	0.06	300	300	85.13 ± 2.98	1.17	85.13 ± 7.70
E-Aug-early	I-Aug	inh_4	2	6	100	0.135	300	300	85.44 ± 3.23	1.17	85.44 ± 7.76
E-Aug-early	VRC-IE	inh_4	0	2	24	0.05	300	99	21.46 ± 1.29	1.12	65.02 ± 6.96
E-Aug-early	I-Aug-BS	inh_4	0	2	150	0.001	300	300	118.31 ± 4.21	1.27	118.31 ± 7.57
E-Aug-early	E-Aug-late	inh_10	0	2	50	0.001	300	300	46.01 ± 1.81	1.09	46.01 ± 6.57
E-Aug-early	E-pons	ex_13	2	4	100	0.002	300	100	63.08 ± 2.85	1.59	189.25 ± 7.68
E-Aug-early	I-pons	inh_7	2	4	100	0.0005	300	100	63.08 ± 2.85	1.59	189.25 ± 7.68
+E-Aug-early	I-Dec_2	inh_4	2	6	150	0.08	300	300	118.31 ± 4.21	1.27	118.31 ± 7.57
+E-Aug-early	Phrenic	inh_4	0	2	150	0.001	300	210	107.57 ± 4.12	1.39	153.68 ± 8.34
+E-Aug-early	Phrenic-HT	inh_4	0	2	150	0.001	300	70	61.97 ± 2.40	2.42	265.57 ± 5.37
+E-Aug-early	ELM	inh_4	2	6	175	0.007	300	300	132.52 ± 4.21	1.32	132.52 ± 9.21
E-Aug-late	E-Aug-early	inh_10	0	2	200	0.04	300	300	145.91 ± 4.71	1.37	145.91 ± 9.07
E-Aug-late	I-Dec	inh_4	2	6	55	0.2	300	300	50.33 ± 1.87	1.09	50.33 ± 6.68
E-Aug-late	I-Aug	inh_4	2	6	100	0.1	300	300	85.27 ± 3.04	1.17	85.27 ± 7.41
E-Aug-late	E-Dec-Phasic	inh_4	2	6	120	0.015	300	300	98.86 ± 3.48	1.21	98.86 ± 9.67
E-Aug-late	I-Aug-BS	inh_4	2	6	150	0.06	300	300	118.14 ± 4.18	1.27	118.14 ± 8.44
E-Aug-late	VRC-IE	inh_4	0	2	24	0.02	300	99	21.32 ± 1.39	1.13	64.62 ± 7.03
E-Aug-late	E-Dec-Tonic	inh_4	0	2	100	0.05	300	300	85.18 ± 3.01	1.17	85.18 ± 6.53
E-Aug-late	E-Dec-pre-ELM	inh_22	2	6	115	0.05	300	300	95.69 ± 3.35	1.20	95.69 ± 7.07
E-Aug-late	ELM	inh_7	2	6	200	0.1	300	300	145.91 ± 4.71	1.37	145.91 ± 9.07
+E-Aug-late	I-Dec_2	inh_4	2	6	150	0.075	300	300	118.14 ± 4.18	1.27	118.14 ± 8.44
+E-Aug-late	Phrenic	inh_4	4	8	150	0.06	300	210	107.09 ± 3.99	1.40	152.98 ± 7.22
+E-Aug-late	Phrenic-HT	inh_4	0	2	150	0.06	300	70	61.83 ± 2.37	2.43	265.00 ± 5.23
+E-Aug-late	ILM	inh_4	0	2	115	0.04	300	300	95.59 ± 3.32	1.20	95.59 ± 8.06
VRC-IE	I-Dec	inh_7	0	4	200	0.035	99	300	146.65 ± 4.99	1.36	48.39 ± 5.40
Raphé 8	Raphé 29	ex_1	0	3	50	0.0125	100	100	39.73 ± 2.28	1.26	39.73 ± 5.01
Raphé 28	Raphé 30	ex_1	0	3	50	0.0125	100	100	39.51 ± 2.47	1.27	39.51 ± 4.93
+I-Aug-BS	Phrenic-HT	ex_1	3	6	18	0.05	300	70	15.99 ± 1.24	1.13	68.51 ± 7.54
+I-Aug-BS	Phrenic	ex_1	3	6	50	0.05	300	210	44.95 ± 2.01	1.11	64.21 ± 6.33
NRM-BötC	rostral IE-pons	inh_7	0	1	100	0.002	300	100	63.41 ± 3.22	1.58	190.23 ± 7.92
NRM-BötC	caudal IE-pons	inh_7	0	1	100	0.002	300	100	63.37 ± 3.02	1.58	190.12 ± 7.75
NRM-BötC	I-pons	ex_13	0	1	100	0.002	300	100	63.25 ± 3.00	1.58	189.74 ± 8.31
+E-Aug-Cough (−)	E-Aug (+)	inh_22	2	6	100	0.05	300	300	85.39 ± 3.06	1.17	85.39 ± 7.75
E-Aug-Cough (−)	E-Aug-BS (+)	inh_22	2	6	100	0.5	300	300	85.34 ± 3.04	1.17	85.34 ± 7.22
+E-Aug-Cough (−)	I-Dec_2	inh_22	0	3	200	0.05	300	300	145.65 ± 4.67	1.37	145.65 ± 8.84
E-Aug-Cough (−)	E-Dec-pre-ELM	inh_4	2	6	100	0.1	300	300	85.22 ± 3.02	1.17	85.22 ± 7.02
+E-Aug-Cough (−)	I-Aug-BS	inh_4	0	3	200	0.025	300	300	145.65 ± 4.67	1.37	145.65 ± 8.84
+E-Aug-Cough (−)	ILM	inh_4	0	3	200	0.025	300	300	146.45 ± 4.63	1.37	146.45 ± 8.96
Lung PSRs	Pump (+)	ex_1	0	3	75	0.015	300	300	66.50 ± 2.43	1.13	66.50 ± 6.72
Lung PSRs	Pump (−)	ex_1	0	3	50	0.015	300	300	46.23 ± 1.78	1.08	46.23 ± 9.43
caudal IE-pons	I-Driver	ex_1	0	5	100	0.001	100	300	85.68 ± 2.78	1.17	28.56 ± 4.45
Pump (–)	E-pons	pre-ex_13	0	4	100	0.99	300	100	63.19 ± 2.88	1.58	189.58 ± 7.28
+Pump (–)	I-Dec_2	inh_4	0	2	25	0.0035	300	300	24.04 ± 0.88	1.04	24.04 ± 5.11
Pump (–)	I-Dec	inh_4	0	2	25	0.0035	300	300	23.98 ± 0.90	1.04	23.98 ± 5.97
Pump (–)	I-pons	pre-ex_13	0	4	100	0.99	300	100	63.53 ± 2.94	1.57	190.58 ± 7.02
Pump (–)	EI-pons	pre-ex_13	2	4	100	0.99	300	100	63.55 ± 2.95	1.57	190.64 ± 7.44
Pump (–)	Lung Def_1s	inh_4	0	4	100	0.02	300	300	85.22 ± 3.09	1.17	85.22 ± 7.46
Pump (–)	rostral IE-pons	pre-ex_13	0	4	100	0.99	300	100	63.63 ± 3.02	1.57	190.88 ± 9.09
Pump (–)	caudal IE-pons	pre-ex_13	0	4	100	0.99	300	100	63.63 ± 3.02	1.57	190.88 ± 9.09
Pump (+)	E-Dec-Phasic	ex_1	0	2	100	0.01	300	300	85.47 ± 2.95	1.17	85.47 ± 8.14
Pump (+)	VRC-IE	ex_1	2	6	100	0.01	300	99	63.10 ± 2.99	1.58	191.20 ± 11.20
Pump (+)	I-Aug	ex_1	0	2	25	0.0	300	300	24.07 ± 0.91	1.04	24.07 ± 4.15
Pump (+)	E-Dec-T	ex_1	0	2	100	0.002	300	300	85.12 ± 3.16	1.17	85.12 ± 6.94
+Pump (+)	I-Dec_2	ex_1	2	6	100	0.005	300	300	85.20 ± 3.05	1.17	85.20 ± 9.59
E-Dec-T	Raphé 29	inh_4	0	3	100	0.001	300	100	63.59 ± 3.05	1.57	190.76 ± 7.20
E-Dec-T	I-Aug-BS	inh_4	2	6	100	0.03	300	300	84.99 ± 3.14	1.18	84.99 ± 7.54
E-Dec-T	rostral IE-pons	ex_13	2	4	100	0.001	300	100	63.26 ± 3.11	1.58	189.77 ± 11.78
E-Dec-T	I-pons	ex_13	2	4	100	0.0005	300	100	63.33 ± 3.02	1.58	189.98 ± 9.72
E-Dec-T	I-pons	inh_7	2	4	100	0.0005	300	100	63.20 ± 3.13	1.58	189.60 ± 7.79
E-Dec-T	rostral IE-pons	inh_4	2	4	100	0.0005	300	100	63.57 ± 3.16	1.57	190.70 ± 7.57
E-Dec-T	caudal IE-pons	ex_13	2	4	100	0.001	300	100	63.20 ± 3.18	1.58	189.61 ± 7.72
E-Dec-T	caudal IE-pons	inh_7	2	4	100	0.0005	300	100	63.40 ± 3.17	1.58	190.20 ± 9.74
E-Dec-T	ELM	inh_4	0	4	100	0.04	300	300	85.04 ± 3.30	1.18	85.04 ± 8.24
E-Dec-T	I-Aug	inh_4	2	6	100	0.0075	300	300	85.12 ± 3.35	1.17	85.12 ± 8.17
+E-Dec-T	I-Dec_2	pre-ex_28	2	6	100	0.2	300	300	84.99 ± 3.14	1.18	84.99 ± 7.54
rostral IE-pons	EI-pons	inh_4	2	4	100	0.03	100	100	63.79 ± 3.24	1.57	63.79 ± 4.98
rostral IE-pons	VRC-IE	ex_1	0	1	100	0.001	100	99	62.90 ± 3.28	1.59	63.54 ± 4.65
rostral IE-pons	E-Dec-Phasic	ex_1	0	5	100	0.02	100	300	85.06 ± 2.76	1.18	28.35 ± 4.15
EI-pons	rostral IE-pons	ex_1	2	4	100	0.002	100	100	63.47 ± 3.20	1.58	63.47 ± 4.78
EI-pons	caudal IE-pons	ex_1	2	4	100	0.002	100	100	63.36 ± 3.47	1.58	63.36 ± 4.41
EI-pons	VRC-IE	ex_1	0	4	50	0.0003	100	99	39.46 ± 2.35	1.27	39.86 ± 5.11
EI-pons	E-Dec-T	ex_1	0	4	100	0.01	100	300	85.16 ± 3.33	1.17	28.39 ± 5.10
E-Dec-pre-ELM	ELM	ex_19	2	6	250	0.0125	300	300	169.56 ± 5.08	1.47	169.56 ± 7.81
Def 2nd (−)	E-Dec-Phasic	inh_4	2	6	100	0.04	300	300	85.17 ± 3.11	1.17	85.17 ± 9.16
+E-Aug (+)	E-Aug-BS (+)	ex_19	2	6	100	0.02	300	300	85.16 ± 2.97	1.17	85.16 ± 7.43
Raphé 8	Raphé 31	inh_4	0	3	50	0.005	100	100	39.38 ± 2.09	1.27	39.38 ± 6.21
Raphé 8	Raphé 32	inh_4	0	3	50	0.005	100	100	39.51 ± 2.47	1.27	39.51 ± 4.93
Raphé 8	E-Aug-BS (+)	inh_22	0	3	400	0.0	100	300	221.45 ± 4.77	1.81	73.82 ± 4.47
Raphé 29	Raphé 30	ex_1	0	3	50	0.01	100	100	39.51 ± 2.47	1.27	39.51 ± 4.93
Raphé 29	E-Dec-T	ex_19	0	3	100	0.15	100	300	84.74 ± 3.12	1.18	28.25 ± 4.23
Raphé 29	E-Dec-Phasic	ex_19	0	3	100	0.2	100	300	84.74 ± 3.12	1.18	28.25 ± 4.23
Raphé 30	Raphé 29	inh_4	0	3	50	0.01	100	100	39.51 ± 2.47	1.27	39.51 ± 4.93
Raphé 32	Raphé 31	inh_4	0	3	50	0.005	100	100	39.51 ± 2.47	1.27	39.51 ± 4.93
Raphé 32	E-Dec-Tonic	inh_22	0	3	100	0.01	100	300	85.28 ± 3.13	1.17	28.43 ± 4.65
Raphé 32	E-Dec-Phasic	inh_22	0	3	100	0.01	100	300	84.74 ± 3.12	1.18	28.25 ± 4.23
Raphé 32	E-Dec-pre-ELM	inh_22	0	3	100	0.01	100	300	84.74 ± 2.97	1.18	28.25 ± 4.67
>+Cough 2nd order (+)	I-Aug-BS	ex_1	2	6	100	0.02	100	300	85.25 ± 2.83	1.17	28.42 ± 5.15
>Cough 2nd order (+)	I-Aug	ex_1	2	6	100	0.0045	100	300	85.25 ± 2.83	1.17	28.42 ± 5.15
>Cough 2nd order (+)	I-Dec	ex_1	2	6	100	0.0045	100	300	85.27 ± 2.89	1.17	28.42 ± 4.68
>+Cough 2nd order (+)	I-Dec_2	ex_1	2	6	100	0.05	100	300	85.54 ± 3.07	1.17	28.51 ± 4.72
>Cough 2nd order (+)	E-Aug-late	ex_1	2	6	100	0.005	100	300	85.25 ± 2.83	1.17	28.42 ± 5.15
>Cough 2nd order (+)	E-Aug-early	ex_1	0	3	100	0.01	100	300	85.07 ± 3.05	1.18	28.36 ± 4.19
>Cough 2nd order (+)	VRC-IE	inh_4	0	3	100	0.2	100	99	63.13 ± 3.05	1.58	63.77 ± 5.20
>Cough 2nd order (+)	caudal IE-pons	ex_1	0	3	100	0.001	100	100	63.59 ± 3.21	1.57	63.59 ± 5.84
>Cough 2nd order (+)	rostral IE-pons	ex_1	0	3	100	0.001	100	100	63.59 ± 3.21	1.57	63.59 ± 5.84
>Cough 2nd order (+)	I-pons	ex_1	0	3	100	0.001	100	100	63.59 ± 3.21	1.57	63.59 ± 5.84
>Cough 2nd order (+)	E-pons	ex_1	2	6	100	0.001	100	100	63.59 ± 3.21	1.57	63.59 ± 5.84
>Cough 2nd order (+)	EI-pons	ex_1	0	3	100	0.001	100	100	63.59 ± 3.21	1.57	63.59 ± 5.84
>Cough 2nd order (+)	E-Dec-pre-ELM	ex_19	2	6	100	0.004	100	300	85.14 ± 3.06	1.17	28.38 ± 3.97
>+Cough 2nd order (+)	E-Aug (+)	ex_1	2	6	100	0.05	100	300	85.25 ± 2.83	1.17	28.42 ± 5.15
>+Cough 2nd order (+)	E-Aug-Cough (−)	ex_1	2	6	100	0.04	100	300	85.25 ± 2.83	1.17	28.42 ± 5.15
>+Cough 2nd order (+)	ILM	ex_1	2	6	100	0.001	100	300	84.92 ± 3.23	1.18	28.31 ± 4.63
E-pons	rostral IE-pons	inh_4	2	4	100	0.0001	100	100	63.17 ± 3.15	1.58	63.17 ± 5.22
E-pons	caudal IE-pons	inh_4	2	4	100	0.0001	100	100	63.47 ± 3.13	1.58	63.47 ± 5.60
E-pons	I-Dec	inh_4	0	1	100	0.008	100	300	85.14 ± 3.03	1.17	28.38 ± 4.11
NRM-pons	I-pons	ex_1	0	4	100	0.015	100	100	63.26 ± 3.28	1.58	63.26 ± 4.40
NRM-pons	I-pons	inh_4	0	4	100	0.05	100	100	63.62 ± 3.11	1.57	63.62 ± 4.61
NRM-pons	I-Driver	ex_1	2	6	100	0.11	100	300	85.10 ± 3.00	1.18	28.37 ± 4.51
NRM-pons	VRC-IE	ex_1	0	1	100	0.01	100	99	63.02 ± 2.53	1.59	63.66 ± 4.67
NRM-pons	I-Aug	ex_1	0	1	100	0.01	100	300	85.21 ± 2.94	1.17	28.40 ± 4.81
NRM-pons	E-Aug-early	ex_1	0	4	100	0.025	100	300	85.82 ± 3.10	1.17	28.61 ± 4.20
NRM-pons	E-Aug-late	ex_1	0	4	50	0.003	100	300	45.82 ± 1.89	1.09	15.27 ± 3.67
NRM-pons	E-Dec-Phasic	ex_1	0	1	100	0.01	100	300	84.87 ± 3.22	1.18	28.29 ± 4.18
NRM-pons	E-Dec-Tonic	ex_1	0	1	100	0.1	100	300	85.35 ± 3.04	1.17	28.45 ± 4.06
NRM-pons	NRM-BötC	inh_4	0	1	100	0.001	100	300	85.11 ± 2.96	1.17	28.37 ± 5.07
E-Aug-BS (+)	Lumbar	ex_1	6	10	25	0.03	300	210	23.59 ± 1.14	1.06	33.70 ± 5.45
+E-Aug-BS (+)	Lumbar-HT	ex_1	3	6	10	0.05	300	70	9.34 ± 0.75	1.07	40.03 ± 3.18
I-pons	rostral IE-pons	ex_1	0	4	100	0.005	100	100	62.93 ± 2.89	1.59	62.93 ± 5.56
I-pons	VRC-IE	ex_1	0	5	100	0.005	100	99	63.61 ± 3.41	1.57	64.25 ± 4.84
I-pons	I-Aug	ex_1	0	4	50	0.005	100	300	46.17 ± 1.67	1.08	15.39 ± 3.39
I-pons	caudal IE-pons	ex_1	0	4	100	0.005	100	100	63.67 ± 2.85	1.57	63.67 ± 4.57
+I-Dec_2	E-Dec-Tonic	inh_4	2	6	125	0.1	300	300	102.25 ± 3.60	1.22	102.25 ± 7.81
+I-Dec_2	E-Dec-pre-ELM	inh_4	2	6	100	0.01	300	300	85.11 ± 3.10	1.17	85.11 ± 7.52
+I-Dec_2	E-Aug-BS (+)	inh_4	2	6	50	0.00005	300	300	46.04 ± 1.73	1.09	46.04 ± 8.13
+I-Dec_2	ELM	inh_4	2	6	100	0.02	300	300	85.61 ± 3.37	1.17	85.61 ± 7.12
+I-Dec_2	E-Aug (+)	inh_4	2	6	100	0.02	300	300	85.19 ± 3.09	1.17	85.19 ± 8.08
+Lung Def_1s	Def 2nd (−)	ex_1	0	3	35	0.03	300	300	33.04 ± 1.38	1.06	33.04 ± 6.68
+Lung DIS_1s	ELM	ex_1	0	3	100	0.09	300	300	84.97 ± 2.93	1.18	84.97 ± 7.40

+, Connection added to the network in Poliaček et al. (.

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

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 (cmH2O); 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 (cmH2O). 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).

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 5A1–C 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 5A2–C 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

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

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 5A1) 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 5B1) 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 5C1).

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

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

Comparison with previous models.

Variable	Unit	Rybak			Poliaček					No speculative				Current
		Mean	p	CV	p	Mean	p	CV	p	Mean	p	CV	p	Mean	CV
EUPNEA
ELM	spk/s	170	0.00*	0.10	0.00*	231	0.00*	0.08	0.00*	19	1.00	0.23	1.00	19	0.23
ILM	spk/s	83	0.00*	0.06	0.00*	49	0.00*	0.07	0.00*	26	0.62	0.29	1.00	19	0.47
LUMBAR	spk/s	24	0.00*	0.05	0.07	135	0.00*	0.08	1.00	10	1.00	0.08	1.00	11	0.08
PHRENIC	spk/s	109	0.00*	0.10	1.00	129	0.00*	0.11	1.00	62	0.33	0.08	1.00	56	0.13
Ti	s	0.663	0.00*	0.04	0.00*	1.506	0.00*	0.07	1.00	1.744	0.03*	0.07	1.00	1.939	0.10
Te	s	1.053	0.00*	0.12	1.00	1.396	0.00*	0.10	1.00	2.905	1.00	0.12	1.00	2.760	0.11
COUGH
ELM	spk/s	302	0.00*	0.33	0.04*	521	0.00*	0.22	0.86	55	1.00	0.07	1.00	55	0.09
ILM	spk/s	238	0.00*	0.11	1.00	154	0.00*	0.14	1.00	37	1.00	0.10	1.00	36	0.10
LUMBAR	spk/s	235	0.00*	0.07	1.00	536	0.00*	0.10	1.00	76	1.00	0.11	1.00	76	0.06
PHRENIC	spk/s	348	0.00*	0.14	0.16	705	0.00*	0.16	0.07	98	1.00	0.04	1.00	97	0.05
Ti	s	0.494	0.00*	0.20	1.00	0.503	0.00*	0.23	1.00	2.471	1.00	0.09	1.00	2.302	0.12
Te	s	0.564	0.00*	0.24	1.00	0.476	0.00*	0.19	1.00	3.490	1.00	0.19	1.00	3.242	0.15

The behaviors of the networks in Rybak et al. (.

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

Outputs of the mechanical model when linked to the current network model (A) and the networks in (B) (Rybak et al., . 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.

Table 3

Synaptic parameters for the network model.

Synapse name	Synapse type	Synapse equilibrium potential (mV)	Synapse time constant (ms)
Ex_1	Excitatory	115.0	1.5
Inh_4	Inhibitory	−25.0	1.5
Inh_7	Inhibitory	−25.0	2.0
Inh_10	Inhibitory	−25.0	1.5
Ex_13	Excitatory	115.0	1.5
Pre-ex_13	Inhibitory (pre-synaptic to Ex_13)	0.0	1.5
Ex_19	Excitatory	115.0	5.0
Inh_22	Inhibitory	−25.0	5.0
Inh_25	Inhibitory	−25.0	4.5
Ex_28	Excitatory	115.0	1.5
Pre-ex_28	Inhibitory (pre-synaptic to Ex_28)	−25.0	3.5