Control of breathing by interacting pontine and pulmonary feedback loops

The medullary respiratory network generates respiratory rhythm via sequential phase switching, which in turn is controlled by multiple feedbacks including those from the pons and nucleus tractus solitarii; the latter mediates pulmonary afferent feedback to the medullary circuits. It is hypothesized that both pontine and pulmonary feedback pathways operate via activation of medullary respiratory neurons that are critically involved in phase switching. Moreover, the pontine and pulmonary control loops interact, so that pulmonary afferents control the gain of pontine influence of the respiratory pattern. We used an established computational model of the respiratory network (Smith et al., 2007) and extended it by incorporating pontine circuits and pulmonary feedback. In the extended model, the pontine neurons receive phasic excitatory activation from, and provide feedback to, medullary respiratory neurons responsible for the onset and termination of inspiration. The model was used to study the effects of: (1) “vagotomy” (removal of pulmonary feedback), (2) suppression of pontine activity attenuating pontine feedback, and (3) these perturbations applied together on the respiratory pattern and durations of inspiration (TI) and expiration (TE). In our model: (a) the simulated vagotomy resulted in increases of both TI and TE, (b) the suppression of pontine-medullary interactions led to the prolongation of TI at relatively constant, but variable TE, and (c) these perturbations applied together resulted in “apneusis,” characterized by a significantly prolonged TI. The results of modeling were compared with, and provided a reasonable explanation for, multiple experimental data. The characteristic changes in TI and TE demonstrated with the model may represent characteristic changes in the balance between the pontine and pulmonary feedback control mechanisms that may reflect specific cardio-respiratory disorders and diseases.


INTRODUCTION
The respiratory rhythm and motor pattern controlling breathing in mammals are generated by a respiratory central pattern generator (CPG) located in the lower brainstem (Cohen, 1979;Bianchi et al., 1995;Richter, 1996;Richter and Spyer, 2001). The pre-Bötzinger complex (pre-BötC), located within the ventrolateral respiratory column (VRC) in the medulla, contains mostly inspiratory neurons (Smith et al., 1991;Rekling and Feldman, 1998;Koshiya and Smith, 1999). The pre-BötC, interacting with the adjacent Bötzinger complex (BötC), containing mostly expiratory neurons (Cohen, 1979;Ezure, 1990;Jiang and Lipski, 1990;Bianchi et al., 1995;Tian et al., 1999;Ezure et al., 2003), represents a core of the respiratory CPG (Bianchi et al., 1995;Tian et al., 1999;Rybak et al., 2004Rybak et al., , 2007Rybak et al., , 2008Rybak et al., , 2012Smith et al., 2007Smith et al., , 2009Rubin et al., 2009;Molkov et al., 2010Molkov et al., , 2011. This core circuitry generates primary respiratory oscillations defined by the intrinsic biophysical properties of respiratory neurons, the architecture of network interactions within and between the pre-BötC and BötC, and the inputs and drives from other brainstem compartments, including the pons, retrotrapezoid nucleus (RTN), raphé, and nucleus tractus solitarii (NTS). It has been suggested (Rybak et al., , 2008Smith et al., 2007) that these external inputs and drives may have a specific spatial mapping onto respiratory neural populations within the pre-BötC/BötC core network, so that changes in these inputs or drives can alter the balance in excitation between key populations within the core network, thereby affecting their interactions and producing specific changes in the respiratory motor patterns observed under different conditions.
Most CPGs controlling rhythmic motor behaviors in invertebrates and vertebrates operate under control of multiple afferent feedbacks and often provide feedback to the sources of their descending and afferent inputs hence allowing feedback regulation of the descending and afferent control signals (Dubuc and Grillner, 1989;Ezure and Tanaka, 1997;Blitz and Nusbaum, 2008;Buchanan and Einum, 2008), and this regulation often operates via presynaptic inhibition (Nushbaum et al., 1997;Ménard et al., 2002;Côté and Gossard, 2003;Blitz and Nusbaum, 2008).
As in other CPGs, afferent feedbacks are involved in the control of the mammalian respiratory CPG and the generation and shaping of the breathing pattern. Many peripheral mechano-and chemo-sensory afferents, including those from the lungs, tracheobronchial tree and carotid bifurcation, provide feedback signals involving in the homeodynamic control of breathing, cardiovascular function, and different types of motor behaviors coordinated with breathing, such as coughing (see Loewy and Spyer, 1990, for review). The NTS is the major integrative site of these afferent inputs. The present study focuses on the mechanoreceptor feedback mediated by pulmonary stretch receptors (PSRs). These mechanoreceptors respond to mechanical deformations of the lungs, trachea, and bronchi, and produce a burst of action potentials during each breath, thereby providing the central nervous system with feedback regarding rate and depth of breathing (see Kubin et al., 2006, for review). Activation of PSRs elicits reflex effects including inspiratory inhibition or expiratory facilitation (representing the so-called Hering-Breuer reflex), enhancement of early inspiratory effort, bronchodilatation, and tachycardia. PSR axons travel within the vagus nerve, and form excitatory synapses in NTS pump cells (Averill et al., 1984;Backman et al., 1984;Berger and Dick, 1987;Bajic et al., 1989;Anders et al., 1993;Kubin et al., 2006). Pharmacological microinjection and lesion studies (McCrimmon et al., 1987;Ezure et al., 1991Ezure et al., , 1998Ezure andTanaka, 1996, 2004;Kubin et al., 2006) suggest that NTS pump cells mediate the Hering-Breuer reflex (lunginflation induced termination of inspiration). Through pump cells, PSR-originating information alters the activity of CPG neurons in manners consistent with their proposed roles in rhythm generation.
Bilateral injections of NMDA antagonists (MK-801 and AP-5) into the rostral pons reversibly increase the duration of inspiration in vagotomized rats, and this increase is dose-dependent (Fung et al., 1994). This suggests that the rostral pons contains neurons with NMDA-receptors participating in the inspiratory off-switch mechanism. Morrison et al. (1994) showed that lesions of the parabrachial nuclei in the decerebrate, vagotomized, unanesthetized rat produced a significant (4-fold) increase in the duration of inspiration and a doubling of the duration of expiration, supporting a role for this pontine area in the regulation of the timing of the phases of respiration. This abnormal breathing pattern is known as apneusis. Administration of MK-801 into the rostral dorsolateral pons was shown to induce apneusis in vagotomized ground squirrels (Harris and Milsom, 2003). Systemic injection of MK-801 increases the inspiratory duration or results in an apneustic-like breathing in vagotomized and artificially ventilated rats (Foutz et al., 1989;Monteau et al., 1990;Connelly et al., 1992;Pierrefiche et al., 1992Pierrefiche et al., , 1998Fung et al., 1994;Ling et al., 1994;Borday et al., 1998). Similarly, Jodkowski et al. (1994) showed that electrical and chemical lesions in the ventrolateral pons produced apneustic breathing in vagotomized rats. At the same time, apneustic breathing is not usually developed if the vagi remained intact and can be reversed by vagal stimulation, suggesting that NMDA receptors are not involved in the pulmonary (vagal) feedback mechanism.  recorded cells in the rostral pons that exhibited respiratory modulation only when lung inflation, via a cycle-triggered pump, was stopped. The emergence of this respiratory-modulated activity suggests that afferent vagal input may have an inhibitory effect on the respiratory modulated cells in the pons (see also Feldman and Gautier, 1976;Cohen and Feldman, 1977). In the same work, it was noticed that this activity had no apparent influence on the tonic discharge of pontine neurons, suggesting that this inhibition might be presynaptic. Dick et al. (2008) recorded several hundred cells in the dorsolateral pons of decerebrate cats, artificially ventilated by a cycle-triggered pump before and after vagotomy. In their experiments, vagotomy led to either an emergence or facilitation of respiratory modulation in the pons. Sustained electrical stimulation of the vagus nerve elicited the classic Hering-Breuer reflex. Systemic or local blockade of NMDA receptors can result in an apneustic breathing pattern (Foutz et al., 1989;Connelly et al., 1992; Fung et al., 1994;Ling et al., 1994;Borday et al., 1998) similar to that demonstrated by pontine lesions or transections. The specifics of feedback control in the brainstem respiratory CPG is that the latter operates under control of two control loops (pulmonary and pontine ones), which both regulate key neural interactions within the CPG, thereby affecting the respiratory rate, respiratory phase durations and breathing pattern, and, at the same time, interact with each other so that each of them may dominate in the control of breathing depending on the conditions and/or the state of the system. Such feedback interactions and a state-dependent feedback control of the CPG may have broader implication in other CPGs in vertebrates and/or invertebrates.
Specifically, our study focuses on the following major feedback loops involved in the control of breathing ( Figure 1A): (1) the peripheral, pulmonary (vagal) loop that controls the medullary rhythm-generating kernel via afferent inputs from PSRs mediated by the NTS circuits, and (2) the pontine control loop, that provides pontine control of the respiratory rhythm and pattern. Our central hypothesis is that both the peripheral afferent and pontine-medullary loops control the respiratory frequency and  phase durations via key medullary circuits responsible for the respiratory phase transitions (onset of inspiration, E-I, and inspiratory off-switch, I-E, see Figure 1A). In addition, these loops interact changing, balancing, and adjusting their control gain via interaction between NTS and VRC and pontine circuits. To investigate the involvement and potential roles of these feedback loops and their interactions with the medullary respiratory circuits we simulated the effects of suppression/elimination of each and both these feedbacks on the respiratory pattern and respiratory phase durations. The results of simulations were compared with the related experimental data and showed good qualitative correspondence hence providing important insights into feedback control of breathing.

SIMULATION PACKAGE
All simulations in this study were performed using a neural simulation package NSM-3.0 developed at Drexel by Drs. Markin, Shevtsova, and Rybak and ported to the high-performance computer cluster systems running OpenMPI by Dr. Molkov. This simulation environment has been specifically developed and used for multiscale modeling and computational analysis of crosslevel integration of: (a) the intrinsic biophysical properties of single respiratory neurons (at the level of ionic channel kinetics, dynamics of ion concentrations, synaptic processes, etc.); (b) population properties (synaptic interactions between neurons within and between populations with random distributions of neuronal parameters); (c) network properties (connectivity strength and type of synaptic interactions, with user-defined or random distribution of connections), (d) morpho-physiological structure (organization of interacting modules/compartments) (see Rybak et al., 2003Rybak et al., , 2004Rybak et al., , 2007Rybak et al., , 2012Smith et al., 2007;Baekey et al., 2010;Molkov et al., 2010Molkov et al., , 2011. NSM-3.0 has special tools for simulation of various in vivo and in vitro experimental approaches, including suppression of specific ionic channels or synaptic transmission systems, various lesions/transections, application of various pharmacological, electrical and other stimuli to particular neurons or neural populations, etc.

MODELING BASIS: NEURONAL PARAMETERS AND IONIC CHANNEL KINETICS
The model presented in this paper continues a previously published series of models of neural control of respiration (Rybak et al., 2004Smith et al., 2007;Baekey et al., 2010;Molkov et al., 2010Molkov et al., , 2011 and, specifically, represents an extension of Smith et al. (2007) model. Following that model, each neuron type in the present model was represented by a population of 20-50 neurons. Each neuron was modeled as a singlecompartment neuron described in the Hodgkin-Huxley (HH) style. These neuron models incorporated the currently available data on ionic channels in the medullary neurons and their characteristics. Specifically, the kinetic and voltage-gated and characteristics of fast (Na) and persistent (NaP) sodium channels in the respiratory brainstem were based on the studies of the isolated pre-BötC neurons in rats (Rybak et al., 2003). The kinetics and steady-state characteristics of activation and inactivation of highvoltage activated (CaL) calcium channels were based on the earlier studies performed in vitro (Elsen and Ramirez, 1998) and in vivo (Pierrefiche et al., 1999). Temporal characteristics of intracellular calcium kinetics in respiratory neurons were drawn from studies of Frermann et al. (1999). Other descriptions of channel kinetics were derived from previous models Smith et al., 2007). Heterogeneity of neurons within each population was set by a random distribution of some neuronal parameters and initial conditions to produce physiological variations of baseline membrane potential levels, calcium concentrations, and channel conductances. A full description of the model and its parameters can be found in the Appendix. All simulations were performed using the simulation package NSM 3.0 (see above). Differential equations were solved using the exponential Euler integration method with a step of 0.1 ms. We utilized the high-performance computational capabilities of the Biowulf Linux cluster at the National Institutes of Health, Bethesda, MD (http://biowulf.nih.gov).

MODEL ARCHITECTURE AND OPERATION IN NORMAL CONDITIONS
The main objective of this study was to investigate the mechanisms underlying control of the mammalian breathing pattern that is generated in the respiratory CPG circuits in the medulla and modulated by two major feedback loops, one involving interactions of medullary respiratory circuits with the lungs, and the other resulting from interactions of these circuits with the pontine circuits contributing to control of breathing ( Figure 1A). We used an explicit computational modeling approach and focused on investigating the anticipated changes in the motor output (activity of the phrenic nerve, PN), specifically the changes in the duration of the inspiratory and expiratory phases under conditions of removal or suppression of the above feedback interactions ( Figure 1A). The full schematic of our model is shown in Figure 1B. While developing this model, we used as a basis and extended the well-known large-scale computational model of the brainstem respiratory network developed by Smith et al. (2007). This basic model focused on the interactions among respiratory neuron populations within the medullary VRC. Similar to that model, the medullary respiratory populations in the present model (see Figure 1B) include (right-to-left): a ramp-inspiratory (ramp-I) population of premotor bulbospinal inspiratory neurons and an inhibitory earlyinspiratory [early-I(2)] population-both in the rostral ventral respiratory group (rVRG); a pre-inspiratory/inspiratory (pre-I/I) and an inhibitory early-inspiratory [early-I(1)] populations of the pre-BötC; and an inhibitory augmenting-expiratory (aug-E) and inhibitory (post-I) and excitatory (post-Ie) post-inspiratory populations in the BötC. As suggested in the previous modeling studies (Rybak et al., 2004Smith et al., 2007), these populations interact within and between the pre-BötC and BötC compartments and form a core circuitry of the respiratory CPG. In addition, multiple inputs and drives from other brainstem components, including the pons, RTN, NTS, and raphé affect interactions within this core circuitry and regulate its dynamic behavior and the motor output expressed in the activity of phrenic nerve (PN). Respiratory oscillations in the basic and present models emerge within the BötC/pre-BötC core due to the dynamic interactions among: (1) the excitatory neural population, located in the pre-BötC and active during inspiration (pre-I/I); (2) the inhibitory population in the pre-BötC providing inspiratory inhibition within the network [early-I(1)]; and (3) the inhibitory populations in the BötC generating expiratory inhibition (post-I and aug-E). A full description of these interactions leading to the generation of the respiratory pattern can be found in previous publications (Rybak et al., 2004Smith et al., 2007). Specifically, during expiration the activity of the inhibitory post-I neurons in BötC decreases because of their intrinsic adaptation properties (defined by the high-threshold calcium and calcium-dependent potassium currents) and augmenting inhibition from the aug-E neurons (Figures 1B and  2A,B). At some moment, the pre-I/I neurons of pre-BötC release from the deceasing post-I inhibition and start firing (Figure 2) providing excitation to the inhibitory early-I(1) population of pre-BötC and the premotor excitatory ramp-I populations of rVRG ( Figure 1B). The early-I(1) population inhibits all postinspiratory and expiratory activity in the BötC leading to the disinhibition of all inspiratory populations including the ramp-I hence completing the onset of inspiration (E-I transition). During inspiration early-I(1) inhibition of BötC expiratory neurons decreases due to intrinsic adaptation properties defined by the high-threshold calcium and calcium-dependent potassium currents (Figure 2). This decrease of inspiratory inhibition leads to the onset of expiration and termination of inspiration (inspiratory off-switch) (Figure 2). In the rVRG, the premotor ramp-I neurons receive excitation from the pre-I/I neurons and drive phrenic motoneurons and PN activity. The early-I(2) population shapes augmenting pattern of ramp-I neurons and PN. The PN projects to the diaphragm ( Figure 1B) hence controlling changes in the lung volume (inflation/deflation) providing breathing.
The architecture of network interactions within the medullary VRC column (i.e., within and between the BötC, pre-BötC and rVRG compartments) in the present model is the same as in the preceding model of Smith et al. (2007). The extension of the basic model in the present study includes: (1) a more detailed simulation of the pontine compartment (in the Smith et al. model, the pontine compartment did not have neuron populations but simply provided tonic drive to medullary respiratory populations), (2) incorporation of suggested interactions between the pontine and medullary populations that form the pontine control loop in the model (Figures 1A,B), and (3) incorporation of the pulmonary (vagal) control loop that included models of the lungs and pump cells in the NTS (Figures 1A,B).
To simulate the pontine feedback loop, we incorporated in the pontine compartment of the model the following populations (see Figure 1B): the excitatory populations of neurons with inspiratory-modulated (I), inspiratory-expiratorymodulated (IEe) and expiratory-modulated (E) activities, and the inhibitory population of neurons with an inspiratory-expiratorymodulated (IEi) activity. As described above, pontine neurons with such types of modulated activity were found in both rat and cat. However, the existing experimental data on intrapontine and pontine-medullary interactions are insufficient and do not provide exact information on the specific connections between these neuron types; they only suggest general ideas and principles for organization of these interactions, such as the possible reciprocal interconnections between the pontine and medullary neurons with similar respiratory-related patterns (see references in the previous paragraph) and the existence of pontine projections to key medullary neurons involved in the respiratory phase switching (such as post-I, see references above). Therefore in the model, respiratory modulation of neuronal activity in pontine populations was provided by excitatory inputs from the medullary respiratory neurons with the corresponding phases of activity within the respiratory cycle. Specifically, the inspiratory modulation activity in the pontine I population was provided by excitatory inputs from the medullary ramp-I population, the IE modulation in the pontine IEe and IEi populations resulted from excitatory inputs from the medullary ramp-I and post-Ie populations, and the expiratory-modulation in the pontine E population was provided by inputs from the medullary post-Ie population. In addition, to simulate the presence of neurons with respiratory modulated phasic and tonic activities, each of the above four population was split into two equal subpopulations with neurons having the same properties and neuronal connections, but differed by tonic drive, which was received only by tonically active subpopulations (not shown in Figure 1B).
In turn, the pontine feedback in the model included (see Figure 1B): (1) excitatory inputs from the pontine I neurons (from both tonic and phasic subpopulations) to the medullary pre-I/I and ramp-I populations; (2) excitatory inputs from the pontine IEe neurons (both tonic and phasic subpopulations) to the medullary post-I population; (3) inhibitory inputs from the pontine IEi neurons (again both subpopulations) to the medullary early-I(1) population; and (4) excitatory inputs from the pontine E neurons (both subpopulations) to the medullary post-I, post-Ie, and aug-E populations. These neuronal connections from pons to medulla (especially pontine inputs to the medullary post-I and pre-I/I populations) allowed the pontine feedback to control operation of the respiratory network in the BötC/pre-BötC core and specifically to control the durations of the respiratory phases and phase switching. Specifically, the connection weights in the model were tuned so that (a) the durations of inspiration (T I ) and expiration (T E ) in the model without vagal feedback would be within the corresponding physiological ranges for the vagotomized rat in vivo (T I = 0.2-0.55 s and T E = 0.8-1.7 s, e.g., see Monteau et al., 1990;Connelly et al., 1992) and (b) after full suppression or removal of the pons, the value of T I would dramatically increase (3-4 times or more) to be consistent with apneusis (Jodkowski et al., 1994;Morrison et al., 1994;Fung and St. John, 1995;St. John, 1998).

PULMONARY (VAGAL) FEEDBACK LOOP
The busting activity of phrenic motoneurons produces rhythmic inflation/deflation of the lungs, which in turn causes rhythmic activation of PSRs projecting back to the medullary respiratory network within the vagus nerve and hence providing pulmonary (vagal) feedback. The activity of pulmonary afferents in the medulla is relayed by the NTS pump (P) cells. To simulate pulmonary feedback loop, we incorporated simplified models of the lungs and PSRs, so that changes in the lung volume were driven by the activity of PN (see Figures 1A,B). The resultant lung inflation activates PSRs that projected back activating the excitatory (Pe) and inhibitory (Pi) pump cells populations in the NTS. The latter finally projected to the VRC and pons (Figure 1B). Hence in the model, both Pe and Pi populations were involved in the Hering-Breuer reflex preventing over-inflation of the lungs. Specifically (Figure 1B), the Pe population excited the post-I population, which was based on the previous experimental data that both lung inflation and electrical stimulation of the vagus nerve produced an additional activation of decrementing expiratory neurons (Hayashi et al., 1996). Following the previous model (Rybak et al., 2004) we suggested that vagal feedback inhibits the early-I(1) population (in this model, via the Pi population). Both these interactions produced a premature termination of inspiration with switching to expiration and a prolongation of expiration.

INTERACTIONS BETWEEN THE LOOPS
As mentioned in the section "Introduction," the respiratorymodulated activity in the pons is usually much stronger in the absence of lung inflation and in vagotomized animals (e.g., see Dick et al., 2008). One explanation for these effects is that the respiratory-modulated activity in the pons is suppressed by vagal afferents via NTS neurons projecting to the pons. There is indirect evidence that this suppression is based on presynaptic inhibition (Feldman and Gautier, 1976;Dick et al., 2008). Therefore in our model, this presynaptic inhibition is provided by the Pi population of NTS and affects all excitatory synaptic inputs from medullary to pontine neural populations ( Figure 1B). Therefore, this presynaptic inhibition suppresses the respiratory modulation in the activities of pontine neurons and reduces the influence of pontine feedback on the medullary respiratory network operation and the respiratory pattern generated.
Because of the lack of specific data, the synaptic weighs of connections from both pump cell populations (Pe and Pi) were set so that (a) significantly reduce the respiratory nodulation in all types of pontine neurons and (b) keep the durations of inspiration and expiration in simulations with vagal feedback intact within their physiological ranges for the rat in vivo (T I = 0.17-0.3 s and T E = 0.3-0.5 s, e.g., see Connelly et al., 1992).

SIMULATION OF VAGOTOMY (PULMONARY FEEDBACK REMOVAL)
Under normal conditions the "intact" model generated the respiratory pattern with the duration of inspiration T I = 0.189 ± 0.046 s and the duration of expiration T E = 0.388 ± 0.064 s (Figures 2, 3A, 4A, and 5A). "Vagotomy" was simulated by breaking the pulmonary feedback, specifically by a removal of afferent inputs from PSRs to the pump cells in the NTS ( Figure 1A). The resultant changes in the activity of different neural populations and in the output respiratory pattern in the model after simulated vagotomy are shown in Figures 3B and 4B. As a result of vagotomy the pump cells (Pi and Pe populations) become silent (only the activity of Pi is shown in Figures 3B and 4B; the activity of Pe population is similar, i.e., it also becomes silent). This eliminates the excitatory effect of lung inflation (PSR) on the post-I population (and post-Ie, pre-I/I, and ramp-I), mediated by Pe, and its inhibitory effect on the aug-E population, provided by Pi ( Figure 1B). This also eliminates the pulmonary (vagal) control of respiratory phase switching and phase durations. However, this breaking of the pulmonary feedback also removes the presynaptic inhibition of all medullary inputs to pontine neural populations (provided in the intact case by the NTS's Pi population) hence increasing respiratory-modulated activities in the pontine neurons involved in the feedback control of the respiratory network operation (Figures 1A,B). This therefore increases the gain of pontine feedback and its role in the control of respiratory phase switching and phase durations. Figure 3 shows that the vagotomy resulted in increases in the respiratorymodulated activity of pontine populations, a prolongation of inspiration (T I = 0.277 ± 0.108 s), and a dramatic increase in the expiratory phase duration (T E = 0.938 ± 0.065 s). Figure 4 shows that the applied vagotomy produced a significant increase of inspiratory (I), inspiratory-expiratory (IE), and expiratory (E) modulation in the activity of the corresponding pontine neurons with tonic activity and releases the corresponding firing in pontine neurons with phasic I, IE, and E activities not active in the intact case.

SIMULATION OF PONTINE FEEDBACK SUPPRESSION WITH AND WITHOUT PULMONARY FEEDBACK
A complete removal of the pons (i.e., a removal of pontine feedback) in the model with an intact pulmonary feedback produced a prolongation of inspiration (T I = 0.337 ± 0.052 s) and a slightly reduced in average (in comparison to the intact model) but highly variable expiratory duration (T E = 0.353 ± 0.159 s) characterized by occasional deletions of aug-E bursts (see Figures 5B and  6A). To compare our simulations with the existing experimental data on the effects of pontine suppression by local injections of MK801, a blocker of NMDA receptors, that might not completely suppress the excitatory synaptic transmission in the pontine neurons and their activity, we also simulated a partial suppression of excitatory synaptic weights in the pontine compartment (e.g., by 25% see Figure 6A). Such partial suppression produced a visible prolongation of inspiration (T I = 0.262 ± 0.028 s with T E = 0.297 ± 0.028 s at 25% suppression, Figure 6A). In contrast to pontine suppression with the intact pulmonary feedback, the same procedures after vagotomy led to a dramatic increase in the average duration of inspiration (making the inspiratory duration highly variable) at relatively constant duration of expiration (Figures 5C and 6A). This prolongation of inspiration after vagotomy increased with the degree of pontine suppression (reducing the weights of excitatory synaptic inputs to pontine neurons) ( Figure 6A) and accompanied by a suppression or full elimination of post-I activity and reduced amplitude of integrated PN ( Figure 5C). Both these features are typical for apneusis (see Cohen, 1979;Wang et al., 1993;Jodkowski et al., 1994;Morrison et al., 1994;Fung and St. John, 1995;St. John, 1998). The durations of inspiration and expiration after vagotomy at different degrees of pontine suppression were the following: T I = 0.437 ± 0.143 s with T E = 0.433 ± 0.030 s at 25% suppression; T I = 0.885 ± 0.339 s with T E = 0.417 ± 0.004 s at 75% suppression; and T I = 571 ± 0.310 s with T E = 0.431 ± 0.003 s at 100% suppression.
The results of our simulations reflecting changes in T I and T E following different combinations of vagotomy with pontine suppression at different degrees are shown together in Figure 6A. Our general conclusions made from these simulations are the following. (1) A suppression of pontine activity with the intact pulmonary feedback leads to a moderate prolongation of inspiration, slight shortening of expiration, and an increase in variability of T E (with 100% pontine suppression). and other typical characteristics of apneusis (suppressed post-I activity and reduced PN amplitude).

COMPARISON WITH EXPERIMENTAL DATA
To test our model, we performed simulation with 25%, 75%, and 100% suppression of the pontine control loop before and after simulated vagotomy (removal of the pulmonary feedback).
The resultant changes in T I and T E are shown in Figure 6A. To compare these simulation results with the related experimental data, we built similar diagrams from the early study of Connelly et al. (1992), which examined spontaneously breathing in Wistar rats during the administration of NMDA blocker MK-801 before and after vagotomy ( Figure 6B). In this study, the experiments on Wistar rats (in contrast to the Sprague-Dawley strain) did not end with apneusis, due to (in our opinion) an insufficient suppression of the pontine feedback by the performed MK-801 injections. Nevertheless, the effects of vagotomy and MK-801 administration on T I and T E before and after vagotomy reported in Connelly et al. study are qualitatively similar to our simulations with 25% suppression of pontine feedback (see Figures 6A,B). Specifically, the 25% pontine suppression in our simulations and the administration of MK-801 in Connelly et al. experiments result in an increase of T I and slight reduction of T E before vagotomy and in a significant prolongation of inspiration after vagotomy. In addition, vagotomy alone without other perturbations in both cases results in an increase of T I and significant prolongation of T E (see Figures 6A,B). Moreover, the changes in the respiratory frequency and the shape and amplitude of integrated phrenic activity after vagotomy and/or pontine suppression in our model are similar to that in the experimental studies with MK-801 administration (Figure 7). The other comparison of our simulations was made with the experimental study of Monteau et al. (1990) performed in anaesthetized vagotomized rats by using MK-801 administration, which results are summarized in Figure 6C. This study did demonstrate that MK-801 application after vagotomy produced switching from a normal breathing pattern to the typical apneusis. The relationships between T I and T E in our simulation after vagotomy and their changes following 100% pontine suppression (apneusis) are similar to these in the Monteau et al. study (see Figures 6A,C).

DISCUSSION
The results of our simulations promote the concept that both pulmonary and pontine feedback loops contribute to the control of the respiratory pattern and, specifically, the durations of inspiration (T I ) and expiration (T E ). Furthermore, our modeling results are consistent with the previous suggestion of specific interactions between these feedback loops, in particular that the PSR afferents involved in the pulmonary control of T I and T E attenuate the gain of the pontine control of these phase durations (via the presynaptic inhibition of excitatory inputs from medullary to pontine populations) (Feldman and Gautier, 1976; Frontiers in Neural Circuits www.frontiersin.org February 2013 | Volume 7 | Article 16 | 9 FIGURE 5 | The effects of pontine suppression before and after simulated vagotomy. Activity of major medullary [post-I, aug-E, early-I(1), pre-I/I, early-I(1), early-I(2), and ramp-I], NTS (Pi) and pontine (I, IEe, and E) neural populations, lung inflation and PN activity under control conditions (A) and following the 100% suppression of pontine activity before (B) and after (C) simulated vagotomy. The activity pattern shown in (C) represents typical apneusis. Vertical dashed line indicate the inspiratory (I) and expiratory (E) phases. See text for details. Cohen and Feldman, 1977;Cohen, 1979;Mörschel and Dutschmann, 2009). Nevertheless, according to our simulations, pontine activity still plays a role in the control of inspiration and expiration even when the pulmonary feedback is intact, although the gain of this pontine control is significantly reduced by the presynaptic inhibition. This presynaptic inhibition is expected to suppress the respiratory modulation in the activity of pontine neurons expressing either tonic or phasic firing patterns (Feldman and Gautier, 1976;Cohen and Feldman, 1977;Cohen, 1979;St. John, 1987Shaw et al., 1989;Dick et al., 1994Dick et al., , 2008Song et al., 2006;Segers et al., 2008), which is reproduced by our model (Figure 4). Also, the model offers a plausible mechanistic explanation for the previous experimental findings that injection of NMDA antagonists in the dorsolateral pons (specifically in the Kölliker-Fuse area) leads to a prolongation of inspiration and to apneusis in the case of a lack of pulmonary feedback (Foutz et al., 1989;Connelly et al., 1992;Pierrefiche et al., 1992Pierrefiche et al., , 1998Fung et al., 1994;Ling et al., 1994;Bianchi et al., 1995;Borday et al., 1998;St. John, 1998). In contrast to previous suggestions and models (Okazaki et al., 2002;Cohen and Shaw, 2004;Rybak et al., 2004;Dutschmann and Herbert, 2006;Mörschel and Dutschmann, 2009;Dutschmann and Dick, 2012), the mechanisms of action of the two feedbacks considered in the current model are not exactly symmetric. Excitatory inputs from both these feedbacks (from PSRs via the NTS's Pe cells, and from the pontine I, IEe, and E populations) activate the ramp-I, pre-I/I, post-Ie, and post-I medullary populations (see Figure 1B). The majority of these excitatory connections are the ones activating the inhibitory post-I population that controls the inspiratory off-switching, i.e., the timing Frontiers in Neural Circuits www.frontiersin.org February 2013 | Volume 7 | Article 16 | 10

FIGURE 6 | Changes in the durations of inspiration (T I ) and expiration (T E ) following pontine suppression and/or vagotomy. (A) Changes in T I
and T E following the simulated pontine suppression at different degrees (25%, 75%, and 100%) before and after (vag. +) vagotomy. (B) Changes in T I and T E in the study of Connelly et al. (1992): diagrams are built for spontaneously breathing Wistar rats under control conditions and after administration of NMDA blocker MK-801 before and after vagotomy. (C) Changes in T I and T E in the study of Monteau et al. (1990) performed in anaesthetized vagotomized rats using MK-801 administration.
of inspiratory phase termination and T I , and those activating the excitatory pre-I/I population which, in a balance with the inputs to post-I, control the onset of inspiration (and T E ). However the effect of these excitatory inputs from the two feedbacks on the medullary circuitry is not identical and depends on the particular synaptic weights and the activity pattern of the inhibitory NTS's Pi cells providing presynaptic inhibition of medullary inputs to the pontine neurons ( Figure 1B). The organization of inhibitory inputs of these feedbacks to the medullary populations in the model is different. While the pulmonary feedback inhibits the aug-E population (via PSRs and Pi cells) causing a complex effect on the respiratory pattern, the pontine IEi population inhibits the early-I(1) population hence promoting expiration, which is clearly seen after vagotomy ( Figure 1B).
It is important to mention that the current model of the medullary core respiratory circuits in the VRC (including the BötC, pre-BötC, and rVRG) used in our model was derived from the model of Smith et al. (2007) without significant changes. Starting with that first publication, this basic model (with necessary additions) was able to reproduce multiple experimental results, including the characteristic changes of the respiratory pattern following a series of pontine and medullary transections and effect of riluzole (persistent sodium current blocker) on the intact and sequentially reduced in situ preparation Smith et al., 2007), the emergence of the additional late-expiratory oscillations in the RTN/parafacial respiratory group (RTN/pFRG) during hypercapnia and interactions between the BötC/pre-BötC and RTN/pFRG oscillators Molkov et al., 2010), the effects of baroreceptor stimulation and the respiratory-sympathetic coupling including this following the intermittent hypoxia (Baekey et al., 2010;Molkov et al., 2011;Rybak et al., 2012), etc. The extended model described here was also able to reproduce the above behaviors, including the biologically plausible changes of membrane potentials and firing patterns of different respiratory neurons ( Figure 2B). The ability of the extended model to reproduce the experimentally observed effects of the two feedback loops provides an additional support for the model of the core respiratory circuits used in all these previous models.
The exact mechanisms of pontine control of breathing are not well-understood and the pontine-medullary connections incorporated in the model are currently speculative. However, the general importance of the pons in the control of the respiratory pattern is well-recognized (see Dutschmann and Dick, 2012, for review). Studies utilizing the classic neurophysiological approaches of lesioning, stimulating and recording neurons have established that the lateral pons influences not only phase duration, phrenic amplitude, and response to afferent stimulation, but also the dynamic changes in respiratory pattern associated with persistent stimuli. For instance, blocking neural activity in the dorsolateral pons not only prolongs inspiration but also blocks the adaptation to vagal stimulation (Siniaia et al., 2000), and the shortening of expiration associated with repeated lung inflation . Thus, the pons is not only intimately involved in the initial response to various stimuli, but also in the complex processes of accommodation and habituation. In the cardiovascular control system, parabrachial stimulation attenuates the NTS response to carotid sinus nerve stimulation by inhibition of NTS neurons receiving these inputs (Felder and Mifflin, 1988).
With normally operating pontine-medullary interactions, the simulated vagotomy results in a prolongation of inspiration and significant increase of the expiratory duration ( Figures 3B and 6A). However, despite these changes, the breathing pattern after vagotomy remains similar to that in eupnea (Figure 3). This maintenance of the eupneic breathing pattern occurs because the control performed by the pulmonary loop is now partly mimicked by the pontine loop, whose gain is increasing after vagotomy, as the latter removes the presynaptic inhibition of medullary inputs to pontine neurons ( Figure 1B). Our model suggests that the pulmonary feedback yet performs the major function in the control of respiratory phase transitions and phase durations, and that a removal of this control loop places the full responsibility for this control on the pontine feedback loop. The complementary role of the pontine and pulmonary feedbacks in control of phase duration (especially T I ) in our model is consistent with the classical interpretation of their function in respiratory control (see Dutschmann and Dick, 2012, for review). In particular, a premature termination of inspiration and switching to expiration can be elicited by stimulation of either the rostral pons or the pulmonary afferents (Bertrand and Hugelin, 1971;Cohen, 1979;Oku and Dick, 1992;Wang et al., 1993;St. John, 1998;Haji et al., 1999;Okazaki et al., 2002;Rybak et al., 2004;Dutschmann and Herbert, 2006). This observation was explained by their common excitatory input on the post-inspiratory neurons in the medullary VRC which are critically involved in this phase transition (Okazaki et al., 2002;Rybak et al., 2004;Dutschmann and Herbert, 2006;Mörschel and Dutschmann, 2009).
Alternatively, our results suggest that the pontine-medullary feedback does not simply function as an "internal pulmonary feedback," performing a redundant function and compensating for the potential loss of vagal input. The specific increase in the variability of T E with the suppression pontine activity and the significant prolongation of T E after vagotomy ( Figure 6A) indicate that the pontine and pulmonary feedbacks differ in the control of T E . Indeed, our modeling results show that these control loops may complement each other in differential control of phase duration and breathing pattern variability. For example, an increase of T E variability with pontine suppression, as seen in Figures 5B and 6A, may be the case during various breathing disorders, such as sleep apnea or ventilator weaning (Tobin et al., 2012). In this connection, the stability of T E can be critically important and is primarily being controlled by the pons. Moreover, the Kölliker-Fuse area of the dorsolateral pons was explicitly identified to contribute to breathing disorders in a mouse model for a neurodevelopmental disease called Rett-syndrome (Stettner et al., 2007;Abdala et al., 2010). Consistent with the many earlier and recent experimental data from cats and rats (Lumsden, 1923;Cohen, 1979;Wang et al., 1993;Jodkowski et al., 1994;Morrison et al., 1994;St. John, 1998), our simulations show that a strong pontine suppression (e.g., 75%) or its removal after vagotomy leads to apneusis, characterized by a significant increase of inspiratory duration and its variability (Figures 5C and 6A). The other specific characteristics of apneusis are a lack of post-inspiratory activity and a reduction of phrenic amplitude during inspiration (Cohen, 1979;Wang et al., 1993;Jodkowski et al., 1994;Morrison et al., 1994;Fung and St. John, 1995;St. John, 1998), which were reproduced in our simulations ( Figure 5C).
Our understanding of interactions between individual components of complex systems is often insufficient to explain emergent properties of these systems. The present study elucidates the important role of two major feedback loops and interactions between them in regulation of the respiratory rate and breathing pattern allowing the brainstem respiratory network to maintain system's homeostasis and adjust breathing to various metabolic and physiologic demands.
The expressions for steady state activation and inactivation variables and time constants are shown in Table A1. The value of maximal conductances for all neuron types are shown in Table A2.
The kinetics of intracellular calcium concentration Ca is described as follows (Rybak et al., 1997): where the first term constitutes influx (with the coefficient k Ca ) and buffering (with the probability P B) , and the second term where B is the total buffer concentration and K is the rate parameter.
The calcium reversal potential is considered a variable and is a function of Ca: E Ca = 13.27 · ln(4/Ca) (at rest Ca = Ca o = 5 × 10 −5 mM and E Ca = 150 mV).
(A6) heterogeneity of neurons within neural populations, the value of E L was randomly assigned from normal distributions using average value ± SD. Leakage reversal potential for all neurons (except for the pre-I ones) was E L = −60 ± 1.2 mV; for pre-I neurons E L = −68 ± 1.36 mV.

MODELING OF LUNGS, PN, AND PSR
The phrenic motoneuron population and phrenic nerve (PN) were not modeled. Integrated activity of the ramp-I population were considered as PN motor output. An increase in lung volume (lung inflation) V was modeled as a low-pass filter of PN activity: where τ V = 100 ms is a lung time constant. The PSR output was considered proportional to the lung inflation V.