Sigh and Eupnea Rhythmogenesis Involve Distinct Interconnected Subpopulations: A Combined Computational and Experimental Study1,2,3

How a single neural network can generate several rhythmic activities at different time scales remains an open question. Here, in addition to the already described reconfiguring process, we propose a new mechanism by which the respiratory network can generate simultaneously two distinct inspiration-related activities (eupnea and sigh) at different frequencies.


Introduction
In several motor systems, the generation of multiple rhythms by a single central network is often attributed to circuit reconfiguring processes. For example, in the crustacean stomatogastric nervous system, the two interacting pyloric and gastric networks under certain neuromodulatory conditions can generate distinctive motor patterns through functional rearrangements of these circuits (Dickinson et al., 1990;Meyrand et al., 1991;Blitz and Nusbaum, 1997;Combes et al., 1999;Faumont et al., 2005). A similar reconfiguration process has also been described in the mouse respiratory network of the pre-Bötzinger complex (preBötC) for gasping motor activity produced under anoxic conditions. Indeed, a decrease in synaptic inhibition enables all post-inspiratory and expiratory neurons (interneurons with functionality related to other phases of the respiratory cycle) to become active during inspiration (Lieske et al., 2000;Paton et al., 2006). Interestingly, the respiratory neural network is also able to generate a third type of respiration-related activity, namely sighs (or augmented breaths). Sighs are biphasic large-amplitude inspiratory events generated intermittently with eupneic activity but at a lower frequency and higher amplitude (Cherniack et al., 1981;Orem and Trotter, 1993;Takeda and Matsumoto, 1997;Lieske et al., 2000;Tryba et al., 2008;Koch et al., 2013;Chapuis et al., 2014). The shape of sigh motor events depends on inhibitory connections among inspiratory neurons, since disturbances of glycinergic inhibition convert sigh bursts from biphasic to monophasic (Lieske et al., 2000;Chapuis et al., 2014).
In contrast to gasps, the generation of sigh activity appears not to derive from a network reconfiguring process. Indeed, the vast majority of the neurons active during eupnea are also active during sighs and no functional remodeling within the network has been identified so far. Intriguingly, two previous in vitro studies have described a very small subpopulation of inspiratory neurons that are active exclusively during sigh-like activity (Tryba et al., 2008;Chapuis et al., 2014). However, the specific role of these cells in the production of sigh bursts remains unclear, largely because the so called sigh-only cells represent a very limited neuronal population (the existing recordings suggest that these neurons constitute less than 5% of the inspiratory neurons of the respiratory network), and thereby are difficult to target for investigation. Furthermore, the excitatory glutamatergic connections are widespread within the respiratory network and likely to be involved in the production of eupneic and sigh bursts. Such a reciprocally and intensively connected network makes distinguishing between sigh and eupnea using classical electrophysiological approaches problematic. Consequently, whether two distinct subpopulations within the respiratory network are responsible for the generation of sigh and eupnea remains an open question. Here, we addressed this issue by combining computational modeling and in vitro experimentation. We designed a two-compartment computational model for sigh and eupnea subpopulations of neurons with several different parameters reflecting distinct burst-generating mechanisms. The sigh subpopulation generates low-frequency rhythm based on slow intracellular Ca 2ϩ oscillations and the eupnea subnetwork generates fast oscillations mainly driven by activation/inactivation of the persistent Na ϩ current. We also tested connections between these two subpopulations and found that a model in which eupnea compartment inhibits sigh compartment and sigh compartment excites eupnea compartment is consistent with experimental data obtained in vitro. Furthermore, we used this model to make several predictions that we tested using in vitro electrophysiological recordings of fictive inspirationrelated activities generated in brainstem slices.

Computational model
The sighϪeupnea computational network combines two single-compartment models representing sigh and eupnea subpopulations, connected through excitatory and inhibitory synapses. We assumed that the two subpopulations are identical except for several parameters ( Table 1). The model for each subpopulation is based on a previous model of inspiratory preBötC neurons (Toporikova and Butera, 2011), where the fast-spiking component has been removed (Rubin and Terman, 2002) and voltage-gated Ca 2ϩ current has been added. Each compartment of the model represents a subpopulation of such preBötC neurons that are interconnected with predominantly excitatory synaptic coupling.
The voltage of each compartment (V) is determined by using a balance of essential burst-generating currents (Toporikova and Butera, 2011): a voltage-gated persistent Na ϩ current (I NaP ) and a K ϩ -dominated passive leak current (I leak ), a non-inactivating voltage-gated calcium current (I Ca ), a voltage-gated hyperpolarization-activated current (I h ), calcium-activated nonspecific cation current (I CaN ), and a synaptic current (I syn ). Therefore, the voltage component of the model is described by the following equation: where the currents are given by: (2) Here, the parameter g represents the maximal conductance of each current (NaP, K, CaN, Ca, H, or synaptic), while m and n are activation functions and h is the inactivation variable. To reduce the number of parameters in the model, the activation of I Ca was approximated with the activation function for I NaP . This fast, non-inactivated current is intended to track average changes in Ca 2ϩ during the burst. Specific parameters for I h were fitted to experimental activation curve in Figure 7A, except for reversal potential, which was chosen to fit experimental data on sigh shape in Figure 7B. The synaptic activation variable s can be defined as: where k syn is synaptic dissociation constant. Inactivation variable h is described by the following equation: where h ϱ (V) is the steady state activation/inactivation curve and h is the voltage-dependent time constant. The steady state activation/inactivation curves and synaptic activation function H(V) are modeled as a sigmoid, with x representing m, n, s, and h: while the time constant is modeled as follows We used a two-pool model to account for Ca 2ϩ fluxes through the cells' plasma membrane and endoplasmic reticulum (ER). Therefore, Ca 2ϩ kinetics are modeled by three equations representing the intracellular Ca 2ϩ ([Ca] i ) balance, ER Ca 2ϩ ([Ca] ER ) balance and IP 3 channel gating variable (l) (Li and Rinzel, 1994).
where is the ratio of ER to plasma membrane surfaces, ϭ A ER / A pm , f i is a constant reflecting the fraction of bound-to-free Ca 2ϩ concentration normalized to effective area (Wagner and Keizer, 1994), A is a scaling constant, K d is the dissociation constant for IP 3 receptor inactivation by Ca 2ϩ , and is the ratio of cytosolic to ER volumes. The flux into the cytosol from the ER (J ER IN ) is regulated by the activity of IP 3 receptors and is defined as: Synaptic g syn ϭ 9 nS; V syn ϭ 0 mV; V ss ϭ Ϫ10 mV; S S ϭ Ϫ5 mV; s ϭ 5 ms; k syn ϭ 1 g syn ϭ 3 nS; V syn ϭ Ϫ70 mV Other Cm ϭ 21 pF New Research where P IP 3 R is the maximum total permeability of IP 3 channels, L IP 3 R is the ER leak permeability, [IP 3 ] is the IP 3 concentration, K I and K a are the dissociation constants for IP 3 receptor activation by IP 3 and Ca 2ϩ , respectively. The flux from the cytosol back to the ER (J ER OUT ) is controlled by the activity of sarco/endoplasmic reticulum Ca 2ϩ ATPase (SERCA) pumps: where V SERCA is the maximal SERCA pump rate, K SERCA is the coefficient for the SERCA pumps. The influx of Ca 2ϩ through the plasma membrane is proportional to the voltage-gated Ca current where ␣ is a proportionality constant.
The outflux of Ca 2ϩ through the plasma membrane is controlled by the activity of PMCA pumps as follows: where V PMCA is the maximum activity of PMCA pumps and K PMCA is the kinetic constant for PMCA pumps.
The values for sigh and eupneic cell populations' parameters are given in Table 1.
The differential equations were solved numerically using Python programming language. Numerical integration was carried out using odeint function of Scypy library with an adaptable time step. This method solves system of ordinary differential equations using lsoda method for stiff or nonstiff systems. To ensure that the model was in a stable oscillatory regime, the first 6 s of simulations were systematically discarded. (The source code for the model is available at http://modeldb.yale.edu/181962).

Experimental methods
All animal procedures were performed in accordance with the authors' university's and the national Committee of Animal Care regulations. Experiments were performed on mouse embryos (either sex) that were obtained from pregnant OF1 females raised in our laboratory's breeding facility. To remain consistent, we performed our physiological experiments at the same developmental stage at which data used to build our model were obtained; i.e., embryonic day (E) 18.5.

Rhythmic slice preparation
Embryonic brainstem transverse slices isolating the preBötC network were obtained using the following procedure: pregnant mice were killed by cervical dislocation on E18.5, the day of the plug being considered to be E0.5. Embryos were excised from their uterine bags, manually stimulated to trigger spontaneous breathing behavior and then kept under a slight heating light until their experimental use. Slice preparations were dissected in cold oxygenated artificial CSF (aCSF) composed of (in mM): 120 NaCl, 8 KCl, 1.26 CaCl 2 , 1.5 MgCl 2 , 21 NaHCO 3 , 0.58 NaH 2 PO 4 , 30 glucose, pH 7.4. First, the hindbrain was isolated from the embryo's body by a rostral section made at the level of the rhombencephalon and a caudal section made at the level of the first cervical roots (between C2 and C4). Second, the isolated hindbrain was placed in a low melting point agar block and carefully oriented in order to proceed to serial transverse sectioning from rostral to caudal, using a vibratome (Leica). A 450-m-thick slice, with its anterior limit set 250-300 m caudal to the more caudal extension of the facial nucleus, was isolated. Other anatomical landmarks such as the wide opening of the fourth ventricle, the presence of the inferior olive, the nucleus ambiguous, and the hypoglossal nucleus were also used to select the appropriate sectioning axis (as also referred to in newborn mice by Ruangkittisakul et al., 2011;see also Fig. 1A). Such slice preparations encompass the region containing a significant portion of the preBötC network capable of spontaneously generating rhythmic fictive inspiratory activities. Slices were then transferred to the recording chamber, rostral surface up, and continuously alimented with oxygenated aCSF at a temperature of 30°C. Preparations were allowed to recover from the slicing procedure over a period of 20 min before any recording sessions were started.

Recording procedures
Recordings of local population activity in slices were performed using glass micropipettes (with a tip diameter of 80-100 m) positioned at the surface of the slice on top of the preBötC network that was located ventral to the nucleus ambiguus (a structure visible when the slice is exposed to direct light from below; Fig. 1A). The micropipettes used as suction electrodes were fabricated from borosilicate glass tubes broken at the tip (Harvard Apparatus), filled with aCSF and connected to a high-gain amplifier (AM Systems). The signals were filtered (bandwidth between 3 Hz and 3 kHz), integrated (time constant 100 ms; Neurolog, Digitimer), and recorded on a computer through a Digidata 1440 interface and with the PClamp10 software (Molecular Devices). The stored files were analyzed offline.
Whole cell patch-clamp recordings were performed under visual control using differential interference contrast. Patch pipettes were fabricated with borosilicate glass tubes using a puller (Sutter Instrument) and filled with a solution composed of the following (in mM): 140 K-gluconate acid, 1 CaCl 2 ·6H 2 O, 10 EGTA, 2 MgCl 2 , 4 Na 2 ATP, 10 HEPES (pH 7.2); and had tip resistances of 5-7 M⍀ when filled with this solution. Electrophysiological signals were recorded using an Axoclamp 2A amplifier (Molecular Devices) and the same digitizing interface and software as described above. Recorded neurons were all located close to the extracellular pipette and were selected on the basis that their discharge was in phase with the population activity (see Fig. 7A). The I h current was activated by applying 2 s duration 10 mV incrementing hyperpolarizing voltage steps from a holding potential of Ϫ50 mV. To trace the corresponding I/V curve, the amplitude of I h was obtained by subtracting the current amplitude measured at the end of the voltage step from the current measured at the beginning of the voltage step.

Calcium imaging
The activity of multiple preBötC neurons was simultaneously monitored using calcium imaging. Slice preparations were first incubated in the dark for 45 min at room temperature in a solution of oxygenated aCSF containing the cell-permeable calcium indicator dye Calcium Green-1 AM (10 M; Life Technologies). After loading, preparations were positioned in the recording chamber, rostral side up. A subsequent 30 min delay was observed to wash out the excess dye and for the preparation to stabilize in oxygenated aCSF at 30°C. Fluorescent signals were captured through a FN1 upright microscope (Nikon) equipped with an epifluorescent illumination system and a fluorescein filter coupled to an Exiblue camera (QImaging). Images (100 ms exposure time) were acquired over periods lasting 120 s and analyzed using software developed by Dr N. Mellen (Mellen and Tuong, 2009).

Pharmacology
Drugs, obtained from Sigma or Tocris, were dissolved in aCSF or DMSO and bath-applied for 15-30 min at their final concentration of 1 M for strychnine; 4 M for cadmium; 10-20 M for Riluzole; 50 M for ZD7288, an I h blocker; 1-10 M for thapsigargin; 50-100 M for cyclopiazonic acid (CPA); and 5-10 M for ryanodine. The drug effects were measured at the end of the exposure period and frequency values are given as mean Ϯ SEM. Statistical significance was assessed by Student's t test, MannϪWhitney, or one-way ANOVA, as appropriate. Mean values were considered to be significantly different at p Ͻ 0.05.

A two-compartment computational model of the preBötC network generates eupneic-and sigh-like activities
In embryonic brainstem slice preparations (Fig. 1A, top), the preBötC respiratory network generates two distinct patterns of rhythmic activity ( Fig 1A; see also Chapuis et al., 2014), as also observed in newborn slices (Lieske et al., 2000). Extracellular population recordings showed the expression of sigh-like bursts at a low frequency (1.3 Ϯ 0.1 burst/min, n ϭ 15) and with a large amplitude ( Fig. 1A, traces, orange stars) and eupnea-like bursts with higher frequency (16.3 Ϯ 0.5 burst/min, n ϭ 15) and smaller amplitude (Fig. 1A, traces, bursts without stars). Each sigh episode was followed by an extended period of eupnea (mean duration 8.4 Ϯ 0.3 s, n ϭ 7), a delay referred to as post-sigh apnea. Both fictive activities were generated spontaneously and concomitantly in control conditions.
We developed a two-compartment mathematical model of the preBötC network that reproduced the pattern of in vitro fictive sigh and eupneic activities. Each compartment of the model represents a group of preBötC neurons connected predominantly through excitatory synapses. Then the rhythmic activity in each compartment reflects a combination of cellular properties and intercellular connections. The sigh and eupnea compartments of the model are identical except for conductance of a persistent sodium current (g NaP ), ER capacity expressed as a ratio of ER membrane to plasma membrane areas (), and activation threshold of a hyperpolarizationactivated current (I h ), (Fig. 1B1, diagrams). All other parameters of the model can be found in Table 1. Based on inspiratory-related activity features observed in in vitro recordings, we assumed that each compartment possesses distinctly different rhythm-generating mechanisms. We hypothesized that the slow sigh rhythm results mainly from slow Ca 2ϩ oscillations, whereas the fast eupnea rhythm follows mainly the kinetics of the persistent Na ϩ current (I NaP ), since it has been shown that the pre-BötC activity relies heavily upon I NaP (Del Negro et al., 2002;Peña et al., 2004;Del Negro et al., 2005;Koizumi and Smith, 2008). To demonstrate the independent oscillations in each compartment, we first plotted the voltage and intracellular Ca 2ϩ ([Ca] i ) of each compartment in an uncoupled model (g sym INH ϭ g syn STIM ϭ 0) as depicted in Figure 1B1. The voltage oscillations in the sigh compartment followed changes in intracellular Ca 2ϩ periodically released out of ER (Fig. 1B1, left two panels). The slow activation of I h led to a depolarization and activation of voltage-gated Ca 2ϩ current, I Ca , which eventually activated inositol 1,4,5-triphosphate receptors (IP 3 R) on the ER membrane, thereby triggering Ca 2ϩ release into the cytosol. The slow inactivation of IP 3 R by Ca 2ϩ terminated the Ca 2ϩ oscillation and induced a prolonged period of hyperpolarization (post-sigh apnea) during which SERCA pumps transported Ca 2ϩ back into the ER. The voltage oscillations in the eupnea compartment followed the kinetics of I NaP (Fig. 1B1, diagram). This current exhibits a fast activation and slow inactivation, which can induce bursting activity in a preBötC model neuron (Butera et al., 1999;Del Negro et al., 2005;Koizumi and Smith, 2008). In contrast to the oscillations in the sigh compartment, Ca 2ϩ influx in eupnea compartment did not play a rhythmogenic role. Although depolarization-activated I Ca led to an influx of Ca 2ϩ during each eupnea oscillation, the ER capacity in continued of sigh (left) and eupnea (right) network models. Except for parameter values for the ER capacity (), I NaP , and I h , the sigh and eupnea models are identical. Intracellular Ca 2ϩ (top) and voltage (bottom) obtained for individual uncoupled compartments (g syn ϭ 0). Left and right panels show outputs for sigh and eupnea compartments, respectively. B 2 , Model with two compartments connected through an inhibitory synapse (blue) from eupnea to sigh and an excitatory synapse (red) from sigh to eupnea subpopulations. The difference in symbol thickness represents the difference in synaptic strength. Bottom trace shows the average voltage trace of the coupled model that is concomitantly generating sigh and eupnea bursts. Orange stars indicate sigh events. XII, Hypoglossal nucleus; io, inferior olive; na, nucleus ambiguus. New Research this model was not large enough to induce Ca 2ϩ release from the ER. Thus, each compartment produced its own intrinsic rhythmic activity: the eupnea compartment generated high-frequency, low-amplitude oscillations while the sigh compartment generated low-frequency, high-amplitude oscillations (Fig. 1B1, bottom panels). We chose relatively short cycle periods for sigh rhythmicity (1.2 min) to better illustrate the shape of the global output. However, the cycle period of sighs could be extended up to 5 min by varying the combination of g NaP and g Ca (data not shown).
Next, we established synaptic connections between the two compartments and found that with an inhibitory synapse from the eupneic subpopulation to the sigh subpopulation and an excitatory synapse from the sigh to eupneic subpopulations (Fig. 1B2, diagram), the output of the model became very similar to the one recorded in vitro (compare bottom traces in Fig. 1A, and B2). We also found that the weight of the inhibitory synaptic influence from the eupneic to sigh subpopulations has to be set considerably higher than the reciprocal excitatory one (see Discussion, below). With these parameters, largeamplitude low-frequency bursts (Fig. 1B2, orange stars) were intermingled with small-amplitude, high-frequency eupneic-like bursts (Fig. 1B2). The frequencies of the sigh and eupnea rhythms were 0.85 and 13.9 bursts/min, respectively, which presents a temporal relationship that is similar to that seen in the experimental data (Fig. 1A).
Although shorter than in in vitro recordings, the model expressed a theoretical post-sigh apnea, consisting of a 23% increase in the duration of the eupneic period immediately after the sigh (5.34 s post-sigh burst compared to 4.32 s post-eupnea burst). In contrast, when we tried to connect the two compartments through other combinations of synaptic interactions (i.e., reciprocally inhibitory, reciprocally excitatory, and inhibitory from sigh to eupnea subpopulations) the model failed to reproduce the typical eupnea-sigh discharge pattern found experimentally (data not shown and see below). Therefore, with two compartments possessing distinct endogenous cellular properties and connected through a specific combination of inhibitory and excitatory synapses, our model can mimic the bimodal motor burst rhythms of the biological preBötC network.

The model reproduces the changes in sigh burst shape after the blockade of inhibitory synapses
In addition to frequency and amplitude differences, sigh motor output events differ from eupneic bursts in terms of shape (Figs. 1, 2A), (Lieske et al., 2000;Tryba et al., 2008;Chapuis et al., 2014). Whereas eupneic activity has a monophasic shape, sighs are defined as biphasic bursts with the initial phase being comparable to a eupneic event and the second phase having larger amplitude ( 2A, the two components of a sigh burst can be separated by blockade of inhibitory glycinergic synaptic connections within the preBötC network. We confirmed this finding in embryonic mouse brain stem preparation in vitro by bath application of strychnine (1 M), which rendered the highamplitude, low-frequency events monophasic in shape (Fig. 2C, blue trace). Blockade of glycinergic synapses also induced an increase in amplitude of both eupneic and sigh bursts (Fig. 2B,C), suggesting that inhibitory glycinergic inputs are able to limit the overall magnitude of bursting activity throughout the inspiratory neuronal population.
We then used the outcome of the strychnine experiments to narrow options for the connections between the sigh and eupnea compartments in our model. By testing in silico four possible connections between the two compartments (all results not shown), we found that only one of these configurations reproduced the observations made in the in vitro strychnine experiments (Fig. 2D): when eupnea inhibited sigh and sigh excited eupnea, the sigh burst shape was biphasic, similar to the experimental in vitro recordings (Fig. 2D, top trace). When the inhibitory connection between the sigh and eupnea compartments was blocked (g syn INH ϭ 0), the sigh shape became monophasic (Fig. 2D, bottom trace; 2E, right), while the average shape of eupnea oscillations remained unchanged (Fig.  2E, left). It should be noted that the data from the strychnine experiment could only be reproduced for a specific ratio of inhibition to excitation of 1:3 (g syn STIM ϭ 9nS, g syn INH ϭ 3nS). Therefore, the biphasic sigh shape in our model reflected the strength of coupling between the two compartments and required both a weaker inhibitory synapse from the eupnea subpopulation to the sigh subpopulation and a stronger reciprocal excitatory connection.

Sighs but not eupnea are abolished by a reduction of calcium conductances and disruption of ER activity
In our model, activity in the eupnea compartment depends mainly upon the activation of I NaP while sigh compartment activity relies on a specific intracellular Ca 2ϩ mechanism that combines the activation of membrane conductances (I CaN and I Ca ) and second messenger pathways including the ER machinery. Because the intracellular source of calcium does not contribute significantly to eupneic rhythm generation in the preBötC (Beltran-Parrazal et al., 2012) and given the design of our sigh compartment, we expected that partially blocking voltage-gated Ca 2ϩ input would affect rhythmogenesis of sigh but not eupnea. We simulated this scenario in silico by reducing the voltage-gated Ca 2ϩ conductance (g Ca ) to 50% of its original value. The resulting population activity is shown in Figure 3. Compared to control (Fig. 3A, top), a reduction in g Ca removed sighs completely, but did not significantly affect eupnea oscillations (Fig. 3A, bottom). To test the physiological correlate, we conducted equivalent in vitro experiments by applying 4 M cad- New Research mium (Cd), which at this low concentration allows partial blockade of a large spectra of unspecified calcium currents (Elsen and Ramirez, 1998). Cd application selectively prevented fictive sigh burst occurrence while inducing a slight but significant increase in the frequency of ongoing fictive eupneic bursts (15.5 Ϯ 0.2 burst/min in control vs 18.4 Ϯ 0.2 burst/min in Cd; n ϭ 8, p Ͻ 0.005; Fig. 3B). These results are consistent with our model and with previous findings showing that sigh generation in the newborn mouse is critically dependent upon the activation voltage-gated calcium channels and more specifically the P/Q-type (Lieske et al., 2000;Lieske and Ramirez, 2006a;Koch et al., 2013).
In a second series of experiments, we investigated whether disrupting Ca 2ϩ sequestration and release from the ER would more specifically affect sigh than eupnea activity. We tested this both in silico and in vitro and results obtained are presented in Figure 3. In order to significantly affect ER activity in vitro, we simultaneously targeted the calcium fluxes coming from and going to the ER by co-applying thapsigargin (10 M; n ϭ 2) or CPA (100 M; n ϭ 3) and ryanodine (10 M). We observed that reducing SERCA activity down to 10% in the model (Fig.  3C) or blocking ER activity using pharmacological agents in vitro ( Fig. 3D; n ϭ 5) prevented sigh from being generated. In such conditions, eupnea was the only inspiration-related activity detectable and was generated at a higher frequency compared to control (11.3 Ϯ 0.8 burst/min in control vs 17.5 Ϯ 0.6 burst/min in SERCA blocked; n ϭ 5, p Ͻ 0.001). Note that when the drugs were applied at lower concentrations (Ͻ10 M for thapsigargin and ryanodine and Ͻ100 M for CPA; n ϭ 6) or when SERCA activity was kept higher than 10% of its maximal value in the model, sighs were still generated (data not shown), suggesting that a minimal ER activity can sustain sigh generation. Altogether, the outcome of our combined in vitro and in silico experiments supports our assumption that Ca 2ϩ entry through voltage-gated channels, together with the intracellular calcium handling machinery, are major components of the sigh-generation mechanisms.

Sigh and eupnea are sensitive to I NaP blockade
We acknowledge that both I CaN and I NaP are important in respiratory rhythmogenesis (Del Negro et al., 2002;Peña et al., 2004;Del Negro et al., 2005;Pace et al., 2007a,b); however, in the present version of our model, we mainly examined the potential role of I NaP in sigh generation only. Although I NaP is expressed in the sigh compartment of our model, its value was set lower than the one required for stable oscillations. In contrast, the generation of rhythmic activity in the eupnea compartment relies mainly on the persistent sodium conductance (g NaP ). Consequently, a reduction of g NaP should affect the eupnea rhythm more than the sigh rhythm. We therefore tested the involvement of I NaP in rhythmic activity in silico by progressively reducing g NaP . As expected, a decrease in g NaP (g NaP ϭ 80% of control) reduced the frequency of eupnea activity without much effect on the sigh frequency (Fig. 4A,C). Complete removal of I NaP abolished oscillations in the eupnea compartment. Since in the sigh compartment, g NaP is activated only during Ca 2ϩ oscillations, its removal decreased the size of sigh bursts that became indistinguishable in amplitude from typical eupneic bursts.
We conducted comparable in vitro experiments by blocking I NaP using Riluzole (Ril; 10-20 M; n ϭ 6). In five preparations, exposure to 10 M Ril had little effect on either sigh or eupnea burst frequencies (p ϭ 0.6). However, when applied at 20 M, Ril completely blocked both fictive rhythms (Fig. 4B,C). In the remaining preparation, both eupnea and sigh activities were blocked in the presence of 10 M Ril. These findings therefore suggest that, in the embryo, in vitro sigh and eupneic rhythm generation rely to a qualitatively similar extent on I NaP expressed within inspiratory preBötC neurons. These results also show that, in the embryo, I NaP plays a major role in respiratory rhythmogenesis, while the I CaN current, being unable to sustain the rhythm when I NaP is blocked, very probably does not play a critical role in rhythm generation, in contrast to what is considered at postnatal stages (Del Negro et al., 2002;Peña et al., 2004;Del Negro et al., 2005;Pace et al., 2007a,b). Thus, our combined in vitro and in silico experiments suggest that, in the embryo, both sigh and eupnea rhythm generation relies on I NaP , with the distinction that g NaP is directly involved in the rhythmogenic mechanisms of the eupnea subpopulation, whereas in the sigh subpopulation, it serves to amplify already existing depolarizing potentials in each cycle.

Sigh frequency is insensitive to changes in extracellular K ؉ concentration
We hypothesized that sighs are generated mainly through ultra-slow (order of minutes) [Ca 2ϩ ] i oscillations. In our model, such metabolic oscillations rely on periodic Ca 2ϩ release from ER store. Although depolarization-activated Ca 2ϩ influx through voltage-gated Ca 2ϩ channels might participate to the generation of these oscillations, the release from ER is only weekly voltage-dependent and thus depolarization should not significantly affect the frequency of such metabolic oscillations. In contrast, oscillations in the eupnea compartment reflect depolarizationdependent kinetics of I NaP and its frequency should directly reflect the changes in depolarization. Therefore, we expected that oscillations in the eupnea compartment of our model would be more sensitive to external K ϩ concentration ([K ϩ ] o ) than oscillations in the sigh compartment. We tested the effect of changing [K ϩ ] o on the sigh and eupnea rhythm frequencies in silico by progressively increasing the reversal potential for the K ϩ dominated leak current (V K ). Figure 5A shows average voltage profiles of the eupneaϪsigh model for four different values of V K . For V K ϭ Ϫ64 mV, the eupnea compartment was below oscillatory threshold but the sigh compartment was still able to generate slow Ca 2ϩ oscillations (0.4 burst/ min). Thus, the model produced only the ultra-slow sigh oscillations and no eupneic activity (Fig. 5A, top panel). Note that this feature might be specific of the developmental stage examined here, since this is not observed at postnatal ages (Ruangkittisakul et al., 2011). As V K increased (V K ϭ Ϫ62 mV), the eupnea compartment reached its threshold and the two-compartment model produced both sigh and eupnea oscillations (Fig. 5A, second panel). The eupnea frequency was 8.1 burst/min and the sigh frequency increased slightly to 0.6 burst/min. A further increase in V K to its control value (V K ϭ Ϫ60 mV) doubled the eupnea frequency (14.8 burst/min), but the frequency of sighs increased only slightly to 0.8 burst/min (Fig. 5A, third panel). For a relatively depolarized value of V K ϭ Ϫ58 mV, we observed a much greater increase in the eupnea frequency compared to sighs, with the former increasing geometrically to 25 burst/min, and the latter increasing linearly to 1.1 burst/min (Fig. 5A, bottom panel). These data clearly demonstrate that for different levels of V K , the range of frequency increase in the eupnea compartment is significantly greater than that of the sigh compartment (Fig. 5B).
To test our modeling prediction in vitro, we applied different concentrations of [K ϩ ] o to brainstem slices in order to change overall levels of neuronal excitability while recording the integrated network output (Fig. 5C). We monitored six slice preparations, which spontaneously produced both fictive eupnea and sigh rhythms under the standard 8 mM [K ϩ ] o (Fig. 5C). In all six preparations, decreasing [K ϩ ] o to 4 mM caused eupneic bursting to cease completely while sigh bursts continued to be generated without any significant decrease in their cycle frequency (0.20 burst/min in control vs 0.14 burst/min in 4 mM [K ϩ ] _o ; p ϭ 0.6; Fig. 5C,D). At 6 mM [K ϩ ] o , however, weak eupneic activity was also present (Fig. 5C) and thereafter showed a significant increase in burst frequency with increasing [K ϩ ] o (1.27 burst/min in 6 mM [K ϩ ] _o , 10.32 burst/min in 8 mM [K ϩ ] _o , and 19.6 burst/min in 10 mM [K ϩ ] _o ; p Ͻ 0.05; Fig. 5C,D). In contrast, sigh bursts were expressed under all [K ϩ ] o tested (Fig. 5C,D) and at frequencies that remained statistically unchanged (0.21 burst/min in 6 mM [K ϩ ] _o and 0.32 burst/min in 10 mM [K ϩ ] _o ; p ϭ 0.9 and 0.5, respectively; Fig. 5D). These data therefore indicate that excitability levels within the preBötC network have differential consequences for the generation of the two types of inspiratory activity, with the production of the slower sigh rhythm requiring less neuronal excitability than the expression of eupneic bursting. Accordingly, our in vitro experiments confirm our modeling prediction that frequency of eupnea is more sensitive to [K ϩ ] o than sighs. Moreover, these results support the modeling assumption that sighs are generated through slow, voltage-insensitive Ca 2ϩ oscillations, whereas eupnea relies mainly on voltage-sensitive ionic conductances.

PreBötC activity in low [K ؉ ] o involves both neuronal subpopulations
Our simulations indicated that in low [K ϩ ] o , the network operated under a sigh-only regime. On this basis, we reasoned that the low [K ϩ ] o constitutes an experimental condition that allows the specific investigation of sighonly neuron behavior. To test this, we examine the activity of the sigh and eupnea compartments of the model individually under low [K ϩ ] o ; i.e., V K ϭ Ϫ64 mV. Unexpectedly, both subpopulations were active in phase in low extracellular [K ϩ ] (Fig. 6A). This effect could be explained by the presence of the model's excitatory synaptic connection from the sigh to eupnea compartments. To test whether this was merely an artifact of our model, we ran in vitro experiments (n ϭ 5) in which we examined the activity of several neurons simultaneously in control conditions ([K ϩ ] o ϭ 8 mM) and in low [K ϩ ] o (4 mM) using calcium imaging. Consistent with our modeling results, we observed that all neurons active during fictive eupnea and sigh activities in control conditions were also active in low [K ϩ ] o (160 cells from five preparations, a mean of 32 cells per slice; Fig. 6B). Together, these experiments demonstrate that the synapse from the sigh subpopulation onto the eupnea subpopulation is strong enough to excite the entire network in low-excitability conditions. These results also indicate how intricately intermingled the two subpopulations are and hence the difficulty with which one can be functionally distinguished from the other.

Sigh generation requires I h current activation
The hyperpolarization-activated inward I h is a conductance reported to play an important role in regulating respiratory cycle frequency (Thoby-Brisson et al., 2000), but its specific role in sigh generation has never been examined. To address this, we first tested whether I h was expressed in embryonic neurons generating both eupnea and sigh. We recorded 44 inspiratory neurons in voltage clamp conditions and found that the vast majority (n ϭ 38; i.e., 86.4%) did indeed express the I h current. The conductance started to activate around Ϫ65 mV and was effectively blocked by ZD 7288 (50 M) (Fig. 7A). Next, we determined in vitro whether I h was involved in fictive sigh generation by blocking its conductance in slices where the preBötC network exhibited a bimodal (eupnea and sigh) discharge pattern. In all preparations examined (n ϭ 10), blocking I h with 50 M ZD 7288 selectively prevented the generation of fictive sighs (Fig. 7B), while eupneic activity persisted, although at a lower   ). B, Schematic representation of the in vitro slice preparation from which an electrophysiological recording of preBötC activity was made simultaneously with calcium imaging performed on the contralateral respiratory network. Right, Image of a Ca 2ϩ -dye loaded slice obtained with a 40ϫ objective. Yellow circles indicate rhythmically active cells. Bottom, Paired recordings of population electrical frequency (13.8 Ϯ 0.4 burst/min in control vs 12.3 Ϯ 0.4 burst/min in ZD 7288; n ϭ 10, p Ͻ 0.05). Comparable results were also found using the less specific I h blocker, caesium (5 mM; n ϭ 4; data not shown). Together, these data indicate that, unlike eupneic activity, the induction of sigh bursts is critically dependent on I h .  Since I h blockade preferentially affects sighs in vitro, it is possible that this current is expressed only in the sigh subpopulation. To determine the distribution of I h between the two compartments of the model, we tested three scenarios in silico: I h expressed in (1) the eupnea subpopulation only, (2) the sigh subpopulation only, or (3) both. The simulation where I h was expressed solely in the eupnea compartment yielded no sigh activity (Fig. 7C, top  trace). The simulation where I h was expressed in the sigh compartment only produced sighs but with a monophasic shape (Fig. 7C, middle trace). This is probably the result of weaker synchronization between eupnea and sigh, due to lower eupneic depolarization and reduced inhibition of sigh. Only when I h was expressed in both compartments were sigh bursts biphasic (Fig. 7C, bottom trace). In accordance with our patch clamp recordings, these results demonstrate that a high percentage of inspiratory neurons in the preBötC network express I h (see above) and that it is not a unique property of sigh cells. In addition, when we simulated I h blockade in silico by setting g h ϭ 0, sigh bursting in our model was abolished and the model generated only eupnea activity (Fig. 7D). Therefore, our modeling data are consistent with the expression of I h in both eupnea and sigh subpopulations and confirms that I h activation is essential for sigh generation.

Voltage-insensitivity of sigh activity is determined by I h activation threshold
The results of our V K experiments clearly indicate that sigh oscillations are only weakly sensitive to [K ϩ ] o . However, the Ca 2ϩ oscillations in the sigh compartment rely partially on voltage-gated Ca 2ϩ influx for the initiation of Ca 2ϩ release from the ER, and thereby should be somewhat sensitive to depolarization. Therefore, we suspected that some additional mechanism(s) is (are) responsible for the voltage-insensitivity of the sigh subpopulation. Since I h is involved in the slow depolarization of the sigh compartment, it is possible that this current was responsible for the compartment's excitability insensitivity. To test this hypothesis, we ran in silico simulations to identify the consequences of changes in [K ϩ ] o in sigh generation in relation to I h activation threshold in the sigh compartment (Fig. 8A). When I h half-activation in sigh compartment was the same as in eupnea compartment (V n sigh ϭ V n eupnea ϭ Ϫ 90 mV), the period of sighs decreased sharply with a V K increase (Fig. 8B, green curve). When I h half activation in sigh compartment became more depolarized (Fig. 8A, red curve), the sigh period still responded to increase in V K , but this dependence was less steep (Fig. 8B, red curve). When the activation curve in sigh compartment was moved even further to the right (V n ϭ Ϫ70 mV; Fig.  8A, black curve), the sigh compartment became insensitive to increase in V K (Fig. 8B, black curve). In contrast, changes in I h activation threshold did not affect the period of eupnea oscillations (Fig. 8C). Thus, our model predicted that the differential sensitivity to overall excitability between the sigh and eupnea might be linked to distinct activation properties of the I h current in the two subpopulations.

Discussion
In this work, we combined experimental and modeling approaches to test the hypothesis that sigh and eupnea are generated by different subpopulations of neurons in the embryonic preBötC network. We developed a twocompartment mathematical model of the respiratory network with eupnea and sigh subpopulations being identical except for three parameters. In the model, we assumed that the activity of a persistent sodium current drives fast oscillations in the eupnea compartment while the release of Ca 2ϩ from intracellular stores generates slow oscillations in the sigh compartment. We established that only the combination of an inhibitory connection from the eupnea to sigh compartments and a reciprocal excitatory connection correctly reproduced the in vitro experimental recordings. Our model also predicted that sighs are less sensitive than eupnea to changes in [K ϩ ] o concentration; this prediction was confirmed experimentally. Finally, we found that the hyperpolarization-activated current, I h , had to be present in both sigh and eupnea subpopulations in Figure 8 The [K ϩ ] o insensitivity of sighs is linked to I h activation characteristics. A, Model tested with three activation curves for I h in the sigh compartment (V h 1/2 ϭ Ϫ 90mV, green;V h 1/2 ϭ Ϫ 80mV, red; V h 1/2 ϭ Ϫ 70mV, black). The I h activation curve for eupnea was kept constant at Ϫ90 mV. B, The cycle period of sighs in response to different levels of V K was measured for three levels of I h half-activation. The lowest half-activation is less sensitive to changes in V K . C, The period of eupnea does not depend on I h half-activation. order to produce typical biphasic sigh bursts. Another prediction of our model is that I h had a higher activation threshold in the sigh subpopulation, suggesting a stronger activation at physiological resting membrane potentials. Overall, our combined in vitro and in silico findings support the conclusion that sigh and eupnea were generated by two distinct subpopulations and provided the basis for continuing to understand the cellular and network mechanisms by which the sigh rhythm is produced.
It has been previously shown that the respiratory network of the preBötC starts to generate sighs at embryonic ages (Chapuis et al., 2014). We mainly used these recent data to build our model, and therefore also performed our in vitro experiments on preparations obtained at embryonic stages. However, the preBötC continues to generate sighs in vitro later in development. It is very likely that changes/maturation/adaptation of the mechanisms underlying sigh generation occur at later developmental stages. For example, it is known that chloride-mediated signaling undergoes significant functional changes during the perinatal period. Because it has been shown that inhibition plays an important role in coordinating the two components of a sigh burst, maturation of glycine signaling might be one of the parameters undergoing developmental changes with potential developmental consequences on sigh rhythmogenesis. In rat, Ren and colleagues (2006) have shown that the transition from excitatory to inhibitory effects of chloride-mediated conductances occurs right prior birth. Here in the mouse, experiments were performed at the corresponding time (the day directly preceding birth), suggesting that chloride-mediated events might already be inhibitory. However, Koch and colleagues (2013) recently showed that sigh generation requires the maturation of some membrane and synaptic properties. Performing physiological experiments at postnatal stages would be informative in testing the relevance of our findings at more mature states.
Although our present results are consistent with a twosubpopulation mechanism underlying fictive eupnea and sigh generation, the alternative hypothesis that both rhythms are generated by a single network cannot be completely excluded. Indeed, if we increase the ER Ca 2ϩ storage in the eupnea compartment (parameter ), the model can produce both large amplitude Ca 2ϩ oscillations and fast I NaP oscillations within the same compartment. Such a mechanism has been reported previously (Jasinski et al., 2013). However, the voltage profile of such single-compartment model is missing several crucial aspects of in vitro sigh recordings. First, the high value of g NaP needed for Ca 2ϩ oscillations reduces the period of sigh-like oscillations to the order of seconds instead of minutes. We found it is impossible to produce minutelong Ca 2ϩ oscillations in this configuration: the excessive depolarization provided by the large value of g NaP activates I Ca too fast and triggers Ca 2ϩ release from the ER a few times per second. Second, we were unable to reproduce the change in the sigh shape from biphasic to monophasic after blocking inhibitory synaptic connectivity in a single-compartment model. And finally, as reported previously (Jasinski et al., 2013), this latter model is able to generate sigh in a stable manner only within a very narrow parameter regime. Altogether, these observations favor the hypothesis of two subpopulations involved in the generation of the two different types of inspiratory activity. Another possibility for sigh/eupnea generation is the involvement of complex network connections among preBötC cells. Indeed it has been shown that network connectivity can lead to self-organized activity on different time scales (Eytan and Marom, 2006;Bullmore and Sporns, 2009). In addition, the existence of specific metabotropic connections that distinguish sighs from eupneic activity has been proposed (Lieske and Ramirez, 2006b). However, in the absence of clear experimental evidence on connectivity among the cells in pre-BötC, we have not made any explicit assumptions on the connectivity within the network.
Although including a number of essential ionic channels and cellular pathways, our model remains simplistic and can only reproduce some of the features of respiratory activity generated by the preBötC network. For instance, although our model displays a prolonged period of inactivity after each sigh, it is still much shorter than the post-sigh apnea observed in our electrophysiological recordings. It is thought that the post-sigh pause in vivo is associated with chemomodulation and/or with HeringϪBreuer reflex inhibition of inspiration via vagal stimulation of stretch receptors. However, because this short pause is also observed in reduced preparations, such as transverse brainstem slice and even island preparations, it is also likely to derive from specific characteristics of the intrinsic membrane properties of inspiratory neurons (Lieske et al., 2000;Tryba et al., 2008;Viemari et al., 2013;Chapuis et al., 2014). Indeed, intracellular recordings have shown that during sighs, inspiratory neurons are activated to depolarized levels that are higher than those occurring during eupneic bursts. This stronger activation of bursting mechanisms thus results in shifting the membrane potential further from triggering thresholds, leading to an increased refractory period that causes the subsequent eupneic burst to be delayed. It is also possible that some mechanisms essential for post-sigh apnea are not explicitly modeled in our simulations and are required to more accurately reproduce this feature of sigh activity. Another result of model simplification is that partial blockade of Ca 2ϩ conductance decreased eupnea frequency in our model, while the slice recordings showed instead an increase. Again, this is probably due to an incomplete list of ionic currents incorporated into our model. One possibility is that Ca 2ϩ -activated K ϩ channels could lead to an overall hyperpolarization and decrease in eupnea frequency. It is also very likely that applying cadmium in vitro does not only affect calcium-dependent conductances in the slice but might also cause changes in other calcium-dependent mechanisms, including in synaptic interactions (Lieske and Ramirez, 2006a;Koch et al., 2013), and is therefore not as specific as changing a single conductance as we were able to do in the model. Third, another difference of the model is that reducing the synaptic inhibitory conductance to zero does not induce an increase in eupneic burst amplitude as observed experimentally. This can be explained by the fact that the synaptic wiring between the two compartments in our model is thus far very simplistic. Indeed, our present data and previously published studies (Winter et al., 2009;Morgado-Valle et al., 2010) suggest that inspiratory neurons receive inhibitory inputs during eupneic bursts that participate in limiting their burst size. In addition, it has to be noted that our data assigned a new functional role to glycinergic neurons within the preBötC network in the coordination of eupnea and sigh rhythms. Fourth, discrepancies also exist regarding the effect of blocking I NaP . While applying riluzole in vitro (at concentrations currently used by others) prevents both eupnea and sigh activities to be generated at comparable concentrations (a feature that might be specific to embryonic ages), blocking the persistent sodium current in silico leads first to a complete blockade of eupnea activity while sighs are still generated but with a smaller amplitude. Once again, this outcome probably reflected an incomplete set of membrane conductances in the different compartments of our model. Finally, it is possible that individual cell properties are more important for sigh generation than our model has assumed. A study by Tryba et al. (2008) indicated that there are some preBötC neurons that generate two different types of bursts at two different frequencies under the conditions of pharmacological blockade of synaptic transmission. These experimental findings suggested that a single preBötC neuron can generate both eupnea and sigh rhythms. However the same study demonstrated that such neurons are extremely rare (less than 1% of experimental cells), which make it highly unlikely that they generate the sigh burst in the network.
The present study also explored the potential role of I h in sigh generation. Our patch-clamp recordings performed on embryonic inspiratory neurons revealed that I h was expressed in a high proportion (ϳ86%) of the neuronal population sampled here. This was corroborated by the finding that I h must be expressed in both compartments of our model (and would therefore be expressed in the majority of network constitutive elements). However several previous studies performed in newborn rodents found that around half of the rhythmically active pre-Bötzinger interneurons are devoid of I h (Rekling et al., 1996;Thoby-Brisson et al., 2000;Picardo et al., 2013). This experimental disagreement might be due to the fact that our experiments were conducted at earlier developmental stage. This suggests that the prevalence for I h might change during development, being highly expressed at more immature states. Such difference also has functional consequences: I h blockade induces opposite effects on eupneic frequency in the embryo versus newborn, triggering a frequency decrease at postnatal stages and an increase before birth. Despite the fact that our data clearly show that I h activation is required for sigh generation in the embryo, it remains possible that this specific role also evolves during later stages of development. Testing this hypothesis would require additional experiment, which is beyond the scope of the study.
Another very important characteristic to be considered is intra-network synaptic connectivity. One primary limitation of our model is that it simulates cell subpopulations and thus ignores individual cell-to-cell interactions. Although we have not modeled explicitly the connections within each subpopulation, we assumed that excitatory synapses among cells are widespread and that blocking these connections would inevitably result in the absence of rhythmic output in our model. Second, our combined in vitro and in silico approaches point to a role of synaptic inhibition in sigh generation. The possible explanations for the monophasic shape of the sigh during strychnine application could be the absence of phase synchronization between sigh and eupnea oscillations, as shown in the newborn (Lieske et al., 2000) and in the embryo (Chapuis et al., 2014) or, in contrast, a complete synchronization of the sigh's two components. The first depolarizing phase of the sigh is a presumed eupnea, which is followed by Ca 2ϩ release from the ER. Such a mechanism would in turn require a precise balance between excitatory and inhibitory synapses. Thus, our model predicted that a mild decrease in synaptic inhibition would extend the time between the two phases of sigh events, thus resulting in separate monophasic bursts. But our model also predicts that a mild increase in synaptic inhibition would lead to a tighter synchronization and loss of the sigh's biphasic shape. A full synchronization is theoretically possible, but we have no experimental evidence supporting this in the embryo.
The initial indication for the possible existence of two subpopulations involved in sigh and eupnea generation was the report of so-called "sigh-only" neurons in two separate studies (Tryba et al., 2008;Chapuis et al., 2014). These neurons, active exclusively during sigh events, were identified either by individual patch-clamp recordings or through calcium imaging, and in each case they constituted a very small population (suggested to be Ͻ5% of the neurons active during respiration-related activities). This population size imbalance was not directly taken into account in our model: the two compartments represented subpopulation entities without simulation of individual constituent elements. Consequently, we do not know how many cells are required to mimic the output of each subpopulation. However, the fact that we can fit our data only when the connection from the eupnea to sigh compartments is stronger than the reciprocal one is consistent with the assumption that eupnea subpopulation is much larger than sigh subpopulation. It is also possible that the two subpopulations have different morphologies. As previously shown, adding a dendritic compartment to a model of persistent sodium pacemaker enables the generation of slow Ca 2ϩ oscillations (Toporikova and Butera, 2011). Thus, the sigh-only cells might have a larger dendritic tree that enables Ca 2ϩ accumulation and slow calcium-dependent intracellular processes. Testing this possibility would obviously require additional experiments and simulations.
One important question is: could a small neuronal population cause synchronization among all neurons of the network? Inspiratory neurons voltage-clamp recordings previously published (Chapuis et al., 2014) and some of our current-clamp recordings (unpublished data) show that during a sigh burst, both discharge frequency of inspiratory neurons and synaptic drive potentials are more intense. It is impossible to determine which one is the consequence of the other, but it illustrates the fact that all neurons of the pre-Bötzinger network are strongly and probably recipro-cally connected through glutamatergic synapses . Therefore, it could be reasonable to consider that even if the sigh-only neuron population might be limited, it could significantly influence the whole network through an efficient recurrent excitatory loop targeting a large proportion of the constitutive neuronal elements of the network.
To our knowledge, our model is the first to accurately reproduce the bimodal discharge pattern of the pre-BötC network generating simultaneously, but with very different frequencies, eupnea and sigh rhythmic burst activities. Although our findings did not completely resolve the question of the mechanism of distinct rhythmic breathing activities, our results definitively support the hypothesis that sighs and eupnea are generated by separate networks within the preBötC. By continuing to combine physiological data with computational simulations, we will gain additional insights on mechanisms underlying sigh generation.