Robustness of respiratory rhythm generation across dynamic regimes

A central issue in the study of the neural generation of respiratory rhythms is the role of the intrinsic pacemaking capabilities that some respiratory neurons exhibit. The debate on this issue has occurred in parallel to investigations of interactions among respiratory network neurons and how these contribute to respiratory behavior. In this computational study, we demonstrate how these two issues are inextricably linked. We use simulations and dynamical systems analysis to show that once a conditional respiratory pacemaker, which can be tuned across oscillatory and non-oscillatory dynamic regimes in isolation, is embedded into a respiratory network, its dynamics become masked: the network exhibits similar dynamical properties regardless of the conditional pacemaker node’s tuning, and that node’s outputs are dominated by network influences. Furthermore, the outputs of the respiratory central pattern generator as a whole are invariant to these changes of dynamical properties, which ensures flexible and robust performance over a wide dynamic range.

Introduction A variety of neuronal circuits, including a range of brainstem and spinal cord central pattern generators (CPGs) in many species, exhibit rhythmic activity patterns. In many CPGs, these patterns consist of sequential activations of different neuronal populations that interact through synaptic connections [1,2]. Significant effort has gone into exploring, using experimental and theoretical methods, the extent to which the intrinsic bursting or pacemaking capabilities of neurons within these populations are responsible for the existence of the network rhythms in which they participate. For example, experimental studies have established the existence of intrinsically bursting neurons in the respiratory pre-Bötzinger complex (preBötC), the inspiratory oscillator within the mammalian brainstem [3], and certain experimental manipulations of burst-supporting conductances in these neurons have eliminated respiratory rhythms [4][5][6][7].
Although these and a wide range of other investigations have explored burst generation mechanisms in respiratory preBötC neurons and have debated the role of this burst generation capability in respiratory rhythms across a variety of conditions [8,9], these studies have not answered a key question: What happens to this intrinsic bursting when the intrinsically burstcapable neurons are embedded within the full network with which they interact synaptically? In fact, it remains unknown whether the intrinsic bursting capabilities of subsets of preBötC neurons affect the dynamics that these neurons actually exhibit once they are embedded within a synaptically interconnected respiratory circuit and how this bursting capability contributes to the properties of the circuit's rhythmic outputs.
These issues relate to the broader concept of robustness. Respiration is critical for survival in mammalian species; respiratory rhythms within individuals must be produced across variations in environmental conditions, metabolic demands, and physiological states, and this persistence must occur across individuals despite their phenotypic heterogeneity. Robustness of a dynamical output such as a respiratory rhythm refers to the maintenance of that pattern of output across conditions or, in a modeling context, across changes in the model parameters or features that represent these conditions. An even stronger form of robustness can arise in which not only the output pattern but also the dynamical mechanisms that give rise to that pattern are preserved. Past experimental and computational work has revealed significant robustness of rhythmicity in the stomatogastric ganglion CPG, which generates certain digestive rhythms in crustacean species, with respect to changes in extrinsic factors such as temperature as well as intrinsic factors such as conductances of ionic currents [10][11][12].
In this study, we use highly reduced models of the respiratory CPG, composed of a small number of units representing important neuron populations [13], to highlight some key ideas relating to robustness of rhythmic outputs across parameter variations, including changes that would induce a switch in intrinsic dynamics if applied to isolated preBötC neurons. The model network generates a three-phase respiratory rhythm and includes inspiratory neurons that abruptly activate at the onset of inspiration, with gradually declining input; a post-inspiratory component that strongly activates at the inspiratory-to-expiratory transition and then gradually diminishes its activity; and a component with activity that represents a later expiratory phase, gradually augmenting after the main post-inspiratory activity surge [14,15]. We model the burst-capable neurons in the respiratory network using differential equations that incorporate a representation of the persistent sodium current, which supports an established oscillatory burst generating mechanism and has been well documented in these neurons experimentally [6,[16][17][18]. Basic principles of dynamical systems theory imply, however, that any alternative mechanisms providing a similar voltage-dependence of activity and autorhythmic capability would yield similar results [19], and hence the ideas that we illustrate are not dependent on a specific bursting mechanism. In particular, our key finding is that network rhythmicity extends smoothly and robustly across different dynamical regimes of the intrinsically burst-capable neurons in the network. Within the multi-phase rhythmic outputs that emerge, these neurons' outputs can become dominated by network influences. While network rhythms exhibit similar features whether the tuning of these neurons renders them intrinsically burst-capable or not, the maintenance of output function across different tunings of intrinsic dynamics yields a broadened dynamic range of network output properties.

Properties of pre-I neuronal unit dynamics are masked within network rhythms
We simulated a highly reduced rhythm-generating network model consisting of ordinary differential equations representing coupled neuronal units. We based our network connectivity on published models for the generation of rhythmic activity patterns in respiratory circuits [5] and we tuned the network to produce rhythmic outputs that are qualitatively similar to the experimentally observed tri-phasic rhythmic respiratory pattern, including the inspiratory (I), the post-I expiratory, and the subsequent augmenting-expiratory (aug-E) stages. The circuit depicted (Fig 1) consists of a preBötC excitatory (glutamatergic) unit (often referred to as pre- I/I because its firing occurs before and during the I phase, but abbreviated as pre-I throughout this work) representing rhythmogenic neurons in the preBötC, a brainstem region that can exhibit intrinsic oscillatory dynamics and drives I activity, along with inhibitory units named after the timing within the respiratory cycle at which their activity occurs (early-I, post-I and aug-E). Among these, it is thought that critical subpopulations of post-I and aug-E inhibitory neurons in respiratory CPG circuits reside in the BötC region adjacent to the preBötC [5,8], and these inhibitory neurons have been shown to provide strong inhibition to the preBötC to regulate inspiratory rhythmicity [20]. In some simulations, as noted in the corresponding text, we also included a unit representing the recently identified glutamatergic post-inhibitory complex (PiCo) [21]. Within this framework, we tuned the parameters in the equations for each unit to represent a distinct activity profile experimentally identified to arise in a corresponding subpopulation of respiratory neurons during the respiratory cycle [5]. In each of these populations, neurons become active and fall silent in a synchronized way but need not spike in synchrony, and hence we omitted spiking currents in our model units.
Within the respiratory rhythm-generating circuit, some members of the pre-I population in the preBötC have been found to be capable of producing rhythmic bursting activity even when isolated from synaptic inputs, at least over some range of constant input currents or extracellular potassium concentrations [3,22,23]. We simulate the conditional bursting property of pre-I neurons by endowing the pre-I unit's voltage equation with a persistent sodium current, I NaP , which is an established feature of preBötC excitatory neurons [6,16,18] and provides a voltage-controlled relaxation-oscillator type of rhythmic oscillation. In the reduced modeling framework, the dynamics produced by the pre-I unit can be understood by using a phase plane representation. The phase plane shows values of voltage, V pre−I , and I NaP inactivation, h pre−I , along with these variables' nullclines, which are curves along which dV pre−I /dt = 0 or dh pre−I /dt = 0, respectively (Fig 2A). Nullcline intersections correspond to fixed points, or combinations of V pre−I , h pre−I for which the unit is at equilibrium. Fixed points may be asymptotically stable, stable, or unstable; in the last case, if we choose a The V pre−I -nullcline structure for the pre-I unit in the preBötC depends on its tonic drive parameter, c 11 . Blue curves denote the V pre−I -nullcline for several values of c 11 , with different line patterns used for contrast enhancement. The red curve is the h pre−Inullcline, which is independent of c 11 . A fixed point occurs where the two nullclines intersect. (B): Bifurcation diagram for the pre-I unit with respect to c 11 . Solid (dashed) curve denotes stable (unstable) equilibria. Red curves are max and min voltages along a family of stable periodic orbits; these grow extremely rapidly in amplitude when they first appear, near c 11 = −0.060. Pre-I is oscillatory where these orbits exist, on c 11 2 (−0.060, −0.011). Inset shows dependence of periods on c 11 ; period jumps sharply from 0 at the onset of oscillations, drops abruptly, and then decays more gradually as c 11 increases.
https://doi.org/10.1371/journal.pcbi.1006860.g002 starting point for (V pre−I , h pre−I ) at random, then we expect these quantities to stay away from the fixed point as time evolves.
Because the tonic excitatory drive to the pre-I unit affects dV pre−I /dt, the V pre−I -nullcline's position and shape depend on the strength of this drive, denoted by c 11 . Generally, increases in c 11 cause nullcline alterations that promote elevations in voltage (Fig 2A). More specifically, for relatively low values of c 11 , the unit's V pre−I -nullcline is cubic or N-shaped and features an asymptotically stable fixed point on its hyperpolarized (or left) branch, corresponding to a quiescent state. As c 11 is increased, the V pre−I -nullcline changes, causing the location of its intersection with the h pre I -nullcline to change as well. As a result, the fixed point crosses to the middle branch and a bifurcation occurs (an Andronov-Hopf, or A-H, bifurcation with a canard explosion [24]) that destabilizes the fixed point, after which the unit will, in isolation from synaptic inputs from other units, produce relaxation oscillations. In this regime, the unit is intrinsically rhythmic, and this rhythmicity depends on I NaP . Finally, with additional increases in c 11 , the V pre−I nullcline switches from cubic to monotonically increasing and another (A-H) bifurcation occurs, producing the loss of oscillations and the restabilization of the fixed point, now at a depolarized voltage, corresponding to a tonically activated state ( Fig  2). In this regime, the excitatory drive to the unit keeps voltage high enough that I NaP deinactivation is limited, including in the tonically activated state. Note that the actual numerical values of c 11 used here are given in arbitrary units.
Once this pre-I unit is embedded within the full respiratory network (as shown in Fig 1, but without the PiCo), its intrinsic properties are masked, or superseded by the network dynamics. There are multiple manifestations of this form of masking. First, at the low end of the c 11 range over which the isolated pre-I unit generates oscillations, the coupled network does not produce oscillatory dynamics but rather settles into a tonic activity regime in which the post-I unit maintains a steady elevated activation level while the pre-I and early-I units exhibit low activation levels. Second, the ranges of periods and amplitudes that the pre-I unit exhibits during its intrinsic oscillations differ from those that it produces when operating within the oscillatory dynamics of the full network (Fig 3). The sensitivity of the pre-I oscillation period and amplitude to changes in c 11 is also altered by the embedding.
Interestingly, the activity pattern of the synaptically connected respiratory network remains remarkably similar regardless of whether c 11 is tuned to yield oscillatory or tonic pre-I dynamics in isolation. These network activity patterns are shown in Fig 4A and 4B for representative c 11 values, one per regime of intrinsic pre-I dynamics, in terms of nondimensional variables (scaled between 0 and 1, with the same scaling used across all c 11 ) representing the output of the pre-I, early-I, post-I, and aug-E units that comprise the network, when PiCo is silent. In both cases, pre-I and to a lesser extent early-I output slowly ramps up before pre-I and early-I activity levels surge together, corresponding to the onset of inspiration. Pre-I and early-I then gradually decline before an abrupt take-over by post-I unit activity and corresponding inhibition of the two I units, marking the phase switch to expiration. After this switch, post-I output declines, disinhibiting aug-E, the output of which ramps up, yielding what is known as the lateexpiratory or E-2 phase, with aug-E output exceeding and inhibiting but not completely suppressing post-I. This similarity in the pattern of network dynamics across different pre-I unit tunings represents another form of masking, since the intrinsic tuning of the pre-I unit dynamics becomes hidden and cannot be deduced from observation of the network activity pattern.
A comparison between the examples shown reveals that the periods of the outputs generated for these parameter values do differ. We will next explore frequency responses to parameter changes across the two regimes of pre-I intrinsic dynamics, and later we will consider network dynamics when PiCo activity is included (Fig 4C and 4D).

Period and amplitude responses to reductions of g NaP
In the model that we are considering, we use a persistent sodium current, I NaP , to allow the possibility of intrinsic rhythmicity in our pre-I unit. In experiments, addition of the I NaP blocker riluzole to the perfusate of an in situ rat preparation that preserved the intact respiratory network yielded little change in respiratory period but caused the output of the phrenic nerve, which corresponds to the excitatory output from the inspiratory neurons in the respiratory circuit, to decrease in amplitude [5]. This decrease occurred roughly linearly with riluzole concentration, reaching about 1/3 of the original amplitude with the maximal riluzole dose applied.
We tested the dependence of the period of our network rhythm (T), the duration of its inspiratory (T I ) and expiratory (T E ) phases, and the amplitude of the pre-I unit's output on the strength of I NaP , represented by its maximal conductance g NaP,exc . Moreover, we evaluated this dependence over a range of values of the pre-I tonic drive parameter c 11 spanning across regimes of both oscillatory and tonic pre-I intrinsic dynamics. Note that the value of c 11 at which the pre-I unit's intrinsic dynamics changes from oscillatory to tonic depends on g NaP,exc . Since this transition occurs at an A-H bifurcation as shown in Fig 2B, we used the continuation capabilities of XPPAUT [25] to follow this bifurcation in the two parameters (g NaP,exc , c 11 ), thereby generating a transition curve. This curve appears in red in all panels of Fig 5. We found that the quantitative properties of the network rhythm did depend on both g NaP,exc and c 11 (Fig 5). For example, T was sensitive to g NaP,exc over a range of c 11 values but less so for sufficiently large c 11 , possibly because units other than pre-I exerted more control over network timing for large c 11 values associated with tonic intrinsic pre-I dynamics (Fig 2). A sufficient reduction of g NaP,exc prevented network output of a rhythmic pattern for most of the range of c 11 considered because pre-I became unable to activate, although for c 11 sufficiently large, rhythmic outputs persisted down to g NaP,exc = 0 nS, despite the lack of intrinsic pre-I rhythmicity at this c 11 level. Wider ranges of T, T I , and T E values could occur for relatively small (i.e., relatively negative) c 11 than for large c 11 (Fig 5A, 5C and 5D), while amplitude changed more with g NaP,exc when c 11 was large (i.e., near 0 or positive, Fig 5B), due to the wider range of g NaP,exc values preserving rhythmicity. Importantly, however, when g NaP,exc and c 11 were varied from one side of the oscillatory-tonic transition curve (red curves in Fig 5) to the other, the transition in pre-I intrinsic dynamics caused absolutely no associated jump or shift in the quantitative properties of the network rhythm ( Fig 5). In other words, the patterns of changes in period, amplitude, T I , and T E with respect to g NaP,exc and c 11 continue uninterrupted across this transition curve.
Finally, we compared the dependence of the model on g NaP,exc to that found experimentally [5]. We found that over a range of c 11 values straddling 0 (c 11 2 [−0.006, 0.004]), we could select a range of g NaP,exc values (g NaP,exc 2 [1.7 nS, 4.1 nS]) such that when g NaP,exc was reduced

Dynamical mechanisms of network rhythms do not depend on intrinsic tuning of preBötC excitatory unit
Although the network's rhythmic activity pattern is similar for both oscillatory and tonic intrinsic preBötC unit dynamics, it remains possible that this intrinsic tuning could affect the detailed dynamical mechanisms by which this pattern is produced, potentially leading to different dynamics in response to certain inputs or perturbations. Since we use model units with dynamics that can be visualized in the phase plane to model each population in the network, we can represent the network oscillation mechanisms geometrically and compare them between the two cases. To do so, although the full network model consists of an eight-dimensional dynamical system, we consider the network outputs projected to various two-dimensional phase planes, each defined using the two variables for one specific unit (Figs 6 and 7). Within these phase planes, we plot the nullclines, or curves of zero rate of change, for each of these two variables. Since each unit's voltage equation depends on the outputs of the other units that are coupled to it synaptically, the location of each projected voltage nullcline varies with the activity of each presynaptic unit. We plot the nullclines at several key moments within a network activity cycle that are particularly helpful for understanding the network dynamics. We also plot the projection of the network output trajectory to each phase plane. Because voltage changes much more rapidly than the persistent sodium inactivation variable, these trajectory projections generally lie on appropriate voltage nullclines except during brief excursions when the network output rapidly changes, which occur in the transitions between inspiratory Black curve is the model output projected to this plane. Solid blue curves are V pre−I -nullclines. The lower right nullcline corresponds to 0 inhibition to pre-I, the upper left nullcline (which continues out of the image) to the maximal inhibition received by pre-I (0.0965), and the middle nullcline to the inhibition level received by pre-I (0.09) when it begins its escape to the active phase by hitting the curve of knees (LK, grey dashed). The red curve is the h pre−I -nullcline, which is independent of inhibition level. The vertical black curve denotes the V pre−I value at which the output of pre-I is half of its maximum value. Solid black circles and squares correspond to time points when inhibition to the pre-I unit is 0.09 and 0.0965, respectively. (B) Same solution projected to the plane in which the inhibition to pre-I is plotted versus h pre−I (black solid). As inhibition wears off, the trajectory first crosses the curve of pre-I fixed points (FP, green dashed) and then crosses the curve of V pre−I -nullcline left knees (LK, blue dashed), allowing escape (black circle) and initiation of inspiration. Trajectories for reduced/increased inhibition from post-I to pre-I (b 31  . From right to left, the V aug−E -nullclines shown correspond to minimal inhibition to aug-E (inhibition is 0.05); to inhibition to aug-E at the moment when pre-I initiates its escape (inhibition is 0.0542, also marked by black circle at same time point as in other plots); to inhibition to aug-E when V pre−I is at its peak (0.283); and to inhibition to aug-E when V post−I is at its peak (0.35, marked by black square at same time point as in other plots).
https://doi.org/10.1371/journal.pcbi.1006860.g006 and expiratory phases. Note that we omit the early-I unit from this discussion, because it always settles to a stable fixed point corresponding to the level of input that it receives, and it lacks the auto-rhythmic capability needed to initiate a transition to the I phase.
The phase plane for the pre-I unit shows that even if the pre-I unit is tuned to be intrinsically oscillatory, at the moment when it receives its maximal level of inhibition, which occurs at the onset of post-I activity (black square in Fig 6A), the unit cannot activate through I NaP deinactivation alone, because its voltage-or V pre−I -nullcline (leftmost blue curve, Fig 6A) intersects the h pre−I -nullcline (red curve) in a stable fixed point at large h pre−I , corresponding to stable quiescence. When inhibition to the pre-I is at its peak, V post−I is at its maximum (black square in Fig 6C), while V aug−E is hyperpolarized due to the strong inhibition that the aug-E unit receives from post-I (black square in Fig 6D). As can be seen in Fig 6C and 6D, during the ensuing post-I phase, V post−I gradually decreases, which occurs due to an inhibitory input from the aug-E unit (which is already present-although very weak-near the onset of the expiratory phase in our model tuning; see Fig 4), and V aug−E gradually increases (Fig 6C and 6D arrows); together, these changes form a self-reinforcing feedback loop. As a result, the inhibitory output from the post-I unit drops, the V pre−I -nullcline position changes, and the value of V pre−I increases correspondingly while h pre−I continues to deinactivate (Fig 6A, arrow).
A key moment in the respiratory cycle is the transition from expiration to inspiration. From the phase plane in Fig 2, recall that when the pre-I unit is tuned to be intrinsically oscillatory, the V pre−I -nullcline has a cubic shape whether it receives any inhibition or not. In a phase plane with a cubic voltage nullcline, a unit transitions from the hyperpolarized phase to the depolarized phase when its voltage rises above the left fold or knee of its voltage nullcline, corresponding to sufficient I NaP deinactivation relative to the level and time course of input that the unit receives. Because the inhibition from post-I to pre-I is stronger than the inhibition from aug-E to pre-I, it follows that the inhibition to pre-I declines over the course of expiration. This decline allows the pre-I unit to reach its curve of knees once V post−I is sufficiently reduced. This event is illustrated in Fig 6A, where the V pre−I -nullcline is shown at the moment (marked with a black circle) when V pre−I reaches this curve of knees (grey curve), allowing the transition to inspiration to commence (see also black circles and corresponding voltage nullclines in Fig 6C and 6D). Another view of this transition is provided by plotting the total level of inhibition to the pre-I unit versus h pre−I (Fig 6B). The (V pre−I , h pre−I ) coordinates of both the knee itself (dashed blue curve) and the fixed point (dashed green curve), where the V pre−I and h pre−I nullclines intersect, depend on this inhibition level. The projection of the trajectory would tend to the curve of fixed points if inhibition were held fixed, preventing the onset of inspiration. But the gradual decrease of total inhibition to pre-I carries the trajectory across the curve of knees, allowing V pre−I to rapidly increase and inspiration to begin. During the subsequent inspiratory phase, the pre-I unit progresses along the right branch of its uninhibited voltage nullcline as I NaP inactivates (rightmost blue curve, Fig 6A; see also sharp rise in pre-I output in Fig 4), while the post-I and aug-E units become strongly inhibited and jump to corresponding inhibited voltage nullclines at hyperpolarized voltages (Fig 6C and 6D). Finally, once V pre−I becomes close enough to its synaptic threshold (black vertical line, Fig 6A), the inhibition to post-I and aug-E weakens, allowing both of their nullclines to shift to depolarized voltages. Both units attempt to activate but post-I wins, yielding a transient spike and subsequent hyperpolarization of the aug-E unit and the transition to expiration, which completes the respiratory cycle.
Note that the only possible moment in this process when the intrinsic oscillation capability of the pre-I unit could have mattered was at the onset of inspiration, which depends on the cubic shape of the V pre−I -nullcline. Yet examining similar phase plane representations when the pre-I unit is tuned to be intrinsically tonic shows that even in this tuning, the V pre−I -nullcline is still cubic at the moment of onset of inspiration and the mechanism underlying the transition to inspiration is the same as in the oscillatory tuning case (Fig 7A and 7B). That is, even though in the absence of inhibition the V pre−I -nullcline becomes monotonic in this tuning, sufficiently strong inhibition can yield a prolonged quiescent period with enhanced I NaP deinactivation, corresponding to a cubic V pre−I -nullcline. When the inhibition to the pre-I unit starts to decay after a period of strong inhibition during expiration, the deinactivation of I NaP can allow the pre-I unit to cross its inhibition-induced knee curve and activate. This mechanism is robust; the V pre−I -nullcline has a cubic shape, with a corresponding left knee, over a wide range of inhibition levels to the pre-I unit (dashed blue curve, Fig 7B). Not surprisingly, we observed that the phase plane projections for the post-I and aug-E units are qualitatively identical between the two pre-I unit tunings as well.

Parameter dependence of network rhythms does not depend on intrinsic tuning of preBötC excitatory unit
The fact that the same dynamical mechanisms underlie network rhythmicity across parameter tunings explains the insensitivity of output period, phase durations, and amplitude to the switch in intrinsic preBötC dynamics as g nap and c 11 are varied (Fig 5). Similarly, we can expect the changes of the network rhythm properties with variations in other network parameters to be independent of intrinsic preBötC tuning. Although we already explored variations in c 11 , we first plotted period (T) and inspiratory (T I ) and expiratory (T E ) phase durations versus c 11 for comparison with effects of other parameter variations ( Fig 8A). As seen previously (Fig 5), there was no abrupt change as the preBötC tuning switched from intrinsically oscillatory to intrinsically tonic. The period T decreased as c 11 increased, via a shortening of T E as also observed experimentally with preBötC stimulation (e.g., [26]), from a maximum of just over 7 seconds to a minimum of about 3 seconds.
Next, since the role of inhibition in respiratory rhythmicity has received significant recent attention [8,20,[26][27][28][29], we considered the effects of varying the strength of the inhibition from the post-I to the pre-I unit, b 31 . This parameter variation produced similar effects whether the pre-I unit was tuned to be intrinsically oscillatory (Fig 8B) or tonic (Fig 8C): increases in b 31 produced increases in overall period, largely via increases in T E but with a small accompanying increase in T I , again yielding a range of periods between about 3 and 7 seconds. The stronger maximal inhibition allowed faster deinactivation of I NaP , represented by a faster increase in h pre−I , than with weaker inhibition and therefore allowed the pre-I unit to cross its curve of left knees and activate at a larger inhibition level than previously (dashed black curves, Figs 6B and 7B). Nonetheless, the amount of time for inhibition to decay from its larger maximum value enough to allow the trajectory to hit the curve of left knees and initiate inspiration was still longer than the escape time with weaker inhibition, leading to an overall longer T E . The stronger I NaP deinactivaton represented by the larger value of h pre−I at the onset of inspiration in this case yielded the small increase in inspiratory duration T I . Similarly, reducing b 31 allowed an earlier crossing of the pre-I curve of left knees to initiate inspiration and a shorter T E (dotted black curve, Figs 6B and 7B).
After these manipulations of parameters associated with the pre-I unit, we varied the drive levels to the other units in the network, one by one, to analyze their abilities to tune network oscillation frequency. In all three cases, these parameter changes produced qualitatively similar results across both tunings of the pre-I unit (Fig 9). Increases in the drive to the early-I unit, c 12 , caused decreases in network period via decreases in T E . In contrast, increases in c 13 (post-I drive) and c 14 (aug-E drive) both increased T, again via changes in T E . Interestingly, the period was much more sensitive to c 14 than to c 13 (also consistent with previous modeling [13]). This sensitivity emerged because in the baseline tuning, the aug-E unit activation was much more modest than the post-I activation, with the aug-E output much farther from maximal. Thus, changes in drive to aug-E influenced the inhibition level to the pre-I unit more than changes in drive to post-I. Finally, because the baseline period was shorter for the intrinsically tonic than the intrinsically oscillatory pre-I tuning, the range of periods in all cases extended to larger values, producing a wider overall range, for the oscillatory regime.

Network responses to changes in synaptic inhibition
Recent experiments considered perturbations of respiratory rhythms resulting from microinjection of GABA A receptor and glycine receptor antagonists into the preBötC or the BötC of rats to disrupt post-synaptic inhibition [29]. The former manipulation led to a decrease in respiratory period and amplitude, as measured by integrated phrenic nerve activity. In an in situ rat brainstem-spinal cord perfused preparation, sufficient blockade of inhibition to the preBötC completely disrupted the respiratory rhythm. Sufficient blockade of inhibition to the BötC, where the subpopulations of post-I and aug-E inhibitory neurons providing input to preBötC circuits are thought to reside [8,20,29], also perturbed respiratory rhythms, in both in vivo and in situ experiments. In contrast to the experiments on the preBötC, however, this manipulation increased respiratory period while decreasing the amplitude of phrenic nerve outputs. In some cases, rhythmic activity was completely suppressed with block of inhibition in the BötC. Responses of our model to reductions in inhibition were similar, but not identical, to the experimental findings.
Our results agreed completely with the experimental findings [29] associated with reduction of inhibition to the pre-I and early-I units of the preBötC, both for intrinsically oscillatory and for intrinsically tonically active pre-I tunings. Indeed, in both cases, period was approximately halved and amplitude, measured in terms of inspiratory output, dropped by about 40% (Fig 10A and 10B), quantitatively consistent with the experimental results. As observed in vivo, the drop in period included shortening of both T I and, especially, T E . The reduction in T E was clearly expected based on the lessened ability of the expiratory units to suppress the inspiratory units, while the decrease of T I resulted from the fact that inspiratory onset occurred with less deinactivation of I NaP , yielding a shorter time needed for I NaP inactivation to occur (Fig 11). The production of a rhythmic output failed at a higher remaining fraction of inhibition when the pre-I unit was intrinsically tonic than when it was oscillatory, because the pre-I and early-I units became tonically active within the network with less of a loss of inhibition than in the former case.
When we decreased inhibition to the post-I and aug-E units, corresponding to experimental perturbations of synaptic inhibition in the BötC, the network did not exhibit an appreciable change in inspiratory output amplitude. When the pre-I unit was intrinsically oscillatory, output period was relatively insensitive to the inhibition level to the BötC within the range over which rhythmicity persisted, although it did increase by about 20% (Fig 10C). With a reduction of inhibition of only about 10%, output rhythmicity was lost. An intrinsically tonic tuning of the pre-I unit provided more robustness to the stronger inhibition from expiratory units that resulted from a decrease in inhibition to the BötC, such that network rhythmicity persisted down to about 65% of normal inhibition levels. Moreover, network period increased by almost 50% over this range of inhibition levels, due to the same effect that was observed to be dominant experimentally, namely an increase in T E , consistent with the mechanism of network rhythmogenesis (Fig 10D). This comparison with experimental results supports the idea that an intrinsically tonic tuning of pre-I neurons may be more consistent with experimental observations than an intrinsically pacemaking tuning, since the tonic tuning provides more robustness against decreased inhibition in the BötC [29].

Network rhythmicity remains similar when an excitatory post-inspiratory complex (PiCo) unit is included
Recently, a glutamatergic post-inspiratory complex (PiCo) including neurons with intrinsic oscillatory capabilities has been identified [21]. As its name suggests, PiCo neurons activate during the post-inspiratory phase. The timing of PiCo activation shifts to inspiration when the inhibition that it normally receives during the I phase is blocked, suggesting that the PiCo receives excitatory inputs from glutamatergic inspiratory neurons. To explore its possible role, we incorporated an intrinsically oscillatory PiCo unit into our model network, sending excitation to the post-I unit, excited by the pre-I unit, and, like the post-I inhibitory unit, inhibited by the early-I and aug-E units. We modeled the PiCo as an intrinsically oscillatory unit [21], with a voltage equation including the same formulation of I NaP used in the pre-I unit and, for simplicity, with the same persistent sodium and leak conductances as the pre-I unit as well.
The inclusion of the PiCo did not qualitatively alter network activity patterns: the PiCo remained suppressed by inhibition from early-I during inspiration and activated together with post-I during the post-inspiratory phase (Fig 4C and 4D). Excitatory output from the PiCo promoted post-I activation, causing the post-I unit to exhibit less adaptation during its active phase and, correspondingly, making the activation of the aug-E unit slightly less gradual (compare Fig 4C and 4D vs. Fig 4A and 4B). Because the network was already rhythmic without the presence of the PiCo, turning on and strengthening the excitation from the PiCo to the post-I unit had little effect on expiratory phase duration, regardless of whether the pre-I unit was intrinsically oscillatory or tonic (Fig 12A and 12B). Examining this effect in more detail, we observe that excitation from PiCo to the post-I unit induces a stronger post-I inhibition of its postsynaptic targets, and of pre-I in particular, while PiCo is active. As a result of this stronger inhibition, the persistent sodium current for the pre-I unit deinactivates more than it would otherwise during the expiratory phase. After the post-I unit falls silent, the pre-I is left with approximately the same inhibitory input (from the aug-E unit) as it would have been without the PiCo, but with slightly larger h pre−I (Fig 12C and 12D; compare grey and magenta circles). With similar inhibition levels, a similar level of persistent sodium deinactivation is needed to allow pre-I activation; with PiCo, since the value of h pre−I at the start of aug-E is larger, the augmenting-expiratory phase ends up being slightly shorter, resulting in a mildly shorter overall cycle period.
Given that the PiCo unit oscillates and drives post-I, it seems theoretically possible that if the pre-I unit is tuned to produce tonic inspiratory activity, then the PiCo could conceivably rescue the network rhythm by activating and recruiting post-I, which could in turn shut down the pre-I and early-I units. In our network configuration, however, extensive parameter explorations with networks tuned to exhibit either sustained pre-I/early-I activation or sustained post-I activation did not reveal evidence that adding a rhythmic PiCo unit to the network without other parameter alterations could induce oscillations. Thus, although this issue needs to be explored further, our findings support the conjecture that a post-inspiratory oscillatory unit cannot on its own induce network oscillations in an otherwise non-rhythmic respiratory network.
Finally, we examined how tuning the drive to the PiCo unit (c 15 ) impacted network outputs and how the presence of the PiCo affected network responses to tuning of drives to other units (c 1i , i < 5). Varying c 15 switched intrinsic PiCo dynamics from complete quiescence (stable fixed point at hyperpolarized voltage) to oscillations to tonic activation (stable fixed point at depolarized voltage; Fig 13A and 13B), as is typical for the voltage-dependent I NaP -mediated oscillation mechanism. For small enough c 15 , even excitation from pre-I could not induce PiCo recruitment at or near the onset of the post-I phase on all rhythm cycles, while for large enough c 15 , PiCo excitatory output allows post-I to suppress aug-E throughout expiration. Between these extremes, c 15 had little effect on network period (Fig 13A and 13B). Although variation of c 15 altered the projection of the network output trajectory to the PiCo phase plane, we found that the PiCo active phase duration remained essentially invariant. Larger c 15 produced slightly stronger PiCo output and hence affected post-I output quantitatively, but PiCo deactivation remained associated with a significant drop in post-I output and corresponding emergence of aug-E activity, with little change in overall timing (*200 msec change over the full range of c 15 ). Similarly, because PiCo activation remained stereotyped, its inclusion produced little impact on network responses to other drive parameters (cf. Fig 13C and 13D, and compare to Fig 9, middle row). One subtle effect that emerged from inclusion of the PiCo is that its output could support post-I activation for lower levels of c 13 than otherwise possible, and, interestingly, could preserve rhythmic network outputs for slightly higher levels of c 13 as well (Figs 13C, 13D and 9, middle row, compare range of c 13 values), through effects of inhibition on I NaP .

Discussion
In this work, we considered a computational network of neuronal units connected synaptically and receiving tonic excitatory synaptic drive, with network parameters tuned to produce a multi-phase rhythmic activity pattern, as a model of components of the mammalian respiratory CPG. Our most fundamental result was the observation that network outputs persist across regimes of tuning of a conditionally rhythmic unit in the network, representing the inspiratory pre-I oscillator in the preBötC. This robustness and insensitivity of qualitative network output pattern to pre-I intrinsic dynamics was not inevitable but also did not require special tuning. These properties emerged because synaptic interactions shaped the network dynamics, altering the pre-I unit's outputs compared to what they would be in the isolated pre-BötC and resulting in similar network dynamical mechanisms across intrinsic tuning regimes for pre-I. While qualitatively similar network outputs were maintained, the ability of the network to function across a broad range of pre-I excitability resulted in an extended dynamic range of output amplitude and period, which could be valuable for function across diverse conditions.

Implications for the respiratory CPG
Our findings can be viewed as exposing general properties of networks incorporating conditional oscillators under certain conditions and are also of specific relevance to an ongoing debate about rhythmicity in the respiratory CPG, which motivated our choices of network connectivity and unit dynamics to study. In the respiratory context, our results point out that a CPG network structure that has been hypothesized based on experimental data [5,13] yields network dynamics that are robust across a range of tunings of the key glutamatergic respiratory neurons in the preBötC. Our findings also imply that experimental manipulations that would alter dynamics or excitability in the isolated preBötC [4,6,26,30,31] cannot be expected to yield analogous results when applied to the connected respiratory network, and that the failure of such manipulations to compromise respiratory outputs does not imply a lack of a role for pre-I dynamics in respiratory rhythm generation. Moreover, our modeling analyses extend our previous results [13] indicating that all of the units participating in rhythmic respiratory activity that are proposed to be core elements of the respiratory CPG can represent important nodes for the control of rhythm features. This control appears to occur seamlessly without altering the basic three-phase organization of activity within the CPG, consistent with experimental observations [20,31] and the notion that input tuning of excitability of different circuit elements is robust in terms of maintaining the coordinated activity necessary for proper behavioral function.
More specifically, our work predicts that increases in tonic drive to pre-I and early-I respiratory neurons should shorten the respiratory cycle period via decreases in the duration of expiration, which is consistent with experimental results [20,26,31]. Increases in tonic drive to post-I and aug-E respiratory neurons are predicted to increase the period, again via changes in expiratory duration, also consistent with experimental observations [20]. We found that tuning of the aug-E population produces a broader range of periods than tuning of post-I; aug-E output in baseline conditions is far from maximal levels, due to lingering inhibition from post-I, and hence increased drive to aug-E can have a significant impact on the timing of inspiratory onset. We also predict that the maximal amplitude of pre-I outputs will be larger in the connected network than in the isolated preBötC, while the minimal network amplitude will also exceed the minimal preBötC amplitude (Fig 3). The inclusion of I NaP in our pre-I model yields oscillations with small amplitude and short period when pre-I excitability is high, as well as very long oscillations just before the transition to quiescence as excitability drops, in the isolated pre-I population. While the connected network cannot capture this full range of periods, the fact that it continues to produce functional outputs across different regimes of pre-I intrinsic dynamics significantly extends its dynamic range relative to what it could produce if restricted to a specific pre-I tuning. Interestingly, our model produces some evidence in support of the idea that in conditions associated with respiratory outputs in the in situ perfused rat preparation, pre-I respiratory neurons are tuned to be tonically active, based on model responses to decreases in g NaP used to simulate partial I NaP blockade (Fig 5). We also found that the tonic pre-I tuning gives more robustness of rhythmic network activity to decreases in inhibition to BötC neurons than does the oscillatory tuning (Fig 10), because the tonic tuning allows inspiratory neurons to activate despite receiving stronger inhibitory outputs from the BötC. Note, however, that a tonic pre-I tuning does not imply that I NaP does not play any role in pre-I activity nor that it is completely inactivated; indeed, it is evident in Fig 7  that h pre−I does not drop to zero during oscillations in this case. On the other hand, we observed that the oscillatory pre-I tuning generally could yield longer respiratory cycle periods (hence extended dynamic range) and associated expiratory phase durations over changes in various drive and synaptic parameters (Fig 9), because with this tuning, the escape of the pre-I unit from inhibition, needed to initiate inspiration, could become strongly delayed. The oscillatory tuning also gave more robustness to decreases in inhibition to preBötC neurons because this tuning yielded more protection against tonic activation of pre-I and early-I within the full network (Fig 10). It is therefore possible that neuromodulatory adjustment of pre-I excitability could be used to transition the network between these regimes, each allowing for particular network tuning mechanisms to become accessible.
The ideas that we explore are consistent with arguments made by Richter and Smith [8] that in the intact respiratory circuit, endogenous pacemaker currents of preBötC neurons need not be autonomously active to drive respiration and indeed may not be active during rhythm generation by the full respiratory CPG circuit due to control of membrane potentials by inhibitory synaptic currents from network connections from the BötC; furthermore, the frequency of the full rhythm may not match the intrinsic frequency of pacemaker neurons in isolation. Diekman et al. have also considered the related ideas that distinct mechanisms may underlie network bursting under different feedback conditions and that an intrinsic tuning of preBötC neurons that does not favor bursting may enhance the robustness of network bursting across normoxic and hypoxic regimes in the presence of chemosensation [31].

Relation to studies of other CPGs and general implications
Related results to ours have recently been demonstrated in a half-center oscillator network, using a similar reduced model that was shown to capture many properties found in simulations of a network composed of a large collection of spiking Hodgkin-Huxley style model neurons [33]. A key point in that work was that rather than the intrinsic dynamics of individual units in the network, the dynamical transition mechanisms by which switches in active populations occurred, involving synaptic effects known as escape and release [13,[34][35][36], determined network output properties. This finding is consistent with our proposal that masking of pre-I properties arises because the dynamic regimes set up by network interactions are similar across pre-I unit tunings.
Our results are also in the same spirit as many years of experimental and modeling results from the CPG in the crustacean stomatogastric ganglion (STG), which generates digestive rhythms. These experiments and simulations have shown that network patterning in the STG is similar across wide ranges of values of parameters associated with the neurons in the circuit [10,12]. Our work carries this idea to the respiratory network and explores additional issues related to the two forms of masking that we have discussed. In particular, we demonstrate that changes in parameter values that alter the intrinsic dynamics of a critical node in a CPG if applied in isolation, for example by inducing a bifurcation, need not produce any noticeable impact on the dynamics of a network in which that node is embedded. Moreover, once that node is embedded in the network, the influence of the network can cause its output tuning to change, effectively masking or overriding its intrinsic dynamics due to circuit interactions.
These principles could be important for the robustness of multi-phase rhythmic outputs of a variety of CPG systems. That is, interactions within a network can yield a preservation of function in the face of changes in conditions that could significantly alter the dynamics of the network components in isolation, as long as these interactions are sufficiently strong to contribute to phase transitions without allowing one node to take over and suppress the others. Furthermore, network interactions can enhance a network's flexibility by allowing for similar rhythms to result from different dynamic mechanisms in different parameter regimes, often leading to a broader range of quantitative output features than would be attainable through reliance on a single rhythm generation mechanism.

Multiple respiratory oscillators?
We performed some simulations with the newly discovered oscillatory post-I population, the PiCo [21], included in our model network, receiving excitatory input from pre-I and inhibition from early-I and aug-E, while providing excitation to the post-I inhibitory unit, based on the proposal that PiCo is critical for generation of post-I activity [21]. Interestingly, our simulations predict that changes in drive to the PiCo have no effect on respiratory period and increases in the strength of excitatory drive from the PiCo to post-I only slightly alter the period, causing it to decrease by a small amount (Fig 12). We modeled the intrinsic dynamics of the PiCo identically to the dynamics of the pre-I unit, based on observations of its intrinsic rhythmicity [21], although it remains for future experiments to test the mechanisms underlying this rhythmicity. Our results raise the question: in what way does the presence of PiCo enhance respiratory network function? According to our simulations, PiCo has limited effectiveness in restoring network rhythmicity if rhythmicity is not present without PiCo, but perhaps PiCo could serve a rescue function in some other circumstances that we have not explored. For example, Anderson et al. did demonstrate an impact on network dynamics by showing that PiCo photostimulation leads to vagal bursts and delays subsequent inspiration [21]. These results are consistent with the idea that stimulating PiCo excitatory neurons drives post-I motor output and can also drive post-I inhibitory neurons to effect the dynamical transition from inspiration to expiration as proposed in our model. The recently proposed triple oscillator hypothesis includes PiCo as a critical component, linked via excitatory coupling with preBötC rhythm-generating neurons and the parafacial respiratory group (pFRG) conditional oscillator, which generates late expiratory activity [37]. We have not included in the present analysis the pFRG oscillator, which is not active and hence not involved in generating the aug-E unit activity during normal eupneic breathing. We have previously (e.g., [38]) analyzed in detail network dynamics when the pFRG and pre-I intrinsic oscillators are coupled and their activity is coordinated by a more extensive inhibitory connectome than represented by the present model, which is intended to represent generation of the normal three-phase rhythmic pattern during eupneic breathing. An interesting future direction would be to analyze the network dynamics in an extended model with coupling between pre-I, PiCo, and pFRG intrinsic oscillators and to evaluate how these interactions impact network performance in conditions (e.g., elevated CO 2 ) under which the pFRG is functionally active.

Modeling assumptions and limitations
Our results were based on simulations of a network of units with two-dimensional dynamics, connected via synaptic coupling. This framework is perfectly valid for illustrating the fundamental point that intrinsic dynamics can become masked by network interactions. It makes the link with the respiratory CPG more abstract, however. To make this connection, we think of each unit within the network as representing the average activity of a population with synchronized activation and inactivation but with asynchronous spiking. The assumed synchrony within the preBötC is consistent with observed eupneic dynamics, presumably reflecting a lack of strong, destabilizing inhibition intrinsic to the preBötC [39]. The model does include the key features of conditional autorhythmicity of the excitatory unit, along with adaptation of the inhibitory units. Furthermore, past work has established that it can produce similar output patterns and even similar responses to changes in drives and coupling to a much more detailed model involving many coupled Hodgkin-Huxley neurons [13,33,40]. Agreement of pre-I outputs with experimental recordings of inspiratory (pre)motor nerve outputs also supports the use of this framework.
We based the conditional intrinsic rhythmicity in our model pre-I unit on the persistent sodium current, I NaP , a well studied current known to exist and to support the possibility of rhythmogenesis in preBötC excitatory neurons. Our results immediately generalize, however, if the pre-I unit instead features other currents or combinations of currents that include one amplifying factor and one resonant (or depolarization-resisting) factor acting on appropriate timescales [19], such that increasing drive yields transitions from quiescence to oscillations to sustained activation via A-H bifurcations (Fig 2B). They may also apply to scenarios where individual pre-I neurons cannot become rhythmic but instead the collective pre-I population can be tuned to produce rhythmicity, either through expansion of dynamic range due to synaptic coupling (e.g., [41][42][43]) or through inclusion of various ionic currents (e.g., [8,44,45]), if the bifurcation structure associated with the activity transitions is similar to what we have considered. Nonetheless, specific predictions about the respiratory network, especially at the quantitative level, may suffer from the abstraction in the model, as well as from the omission of other currents, including the CAN current, which have been shown to contribute to respiratory outputs [46][47][48]. Our model also cannot represent the heterogeneity that is naturally present within biological populations of neurons, which has been predicted to add robustness to respiratory outputs [41][42][43]. Indeed, with heterogeneity, it is likely that different neurons have different ionic conductances leading to different intrinsic dynamics [43]. In this case, the results of our study lead us to conjecture that the maintenance of functional outputs from the full respiratory network should be robust to some variations in these intrinsic dynamical properties, although we recognize that significant heterogeneties can lead to important dynamical effects (e.g., [40,[49][50][51]). Another caveat of our results is that we saw loss of network rhythmicity with sufficient block of inhibition in all regimes considered, whereas outcomes from experiments involving inhibitory block have been more diverse [27,29,52]. In reality, experiments can show variable perturbations from disruption of inhibition (e.g., whether or not rhythm is disrupted) across different experimental trials, possibly corresponding to different levels of inhibition block achieved in different experiments.

Conclusion
CPGs subserve many functions critical to survival, including mammalian respiration. It is critical for functional outputs to persist across wide ranges of metabolic and environmental conditions including diverse developmental and disease states. Maintenance of function as well as control of outputs to allow effective behavior across conditions likely involve significant alterations in feedback and top-down drives to network components and in neuromodulatory influences [53]. Hence, even if particular components within a CPG circuit exhibit intrinsic rhythmicity for certain drive and neuromodulator levels, their dynamics may change as these levels vary. Therefore, the invariance of network outputs and tuning under changes in intrinsic dynamics that we highlight is likely a critical feature of maintaining stable CPG performance across conditions. In the case of the respiratory network, this invariance ensures flexible and robust performance over a wide dynamic range of rhythm generation, which is likely a critical property for such a vital homeostatic CPG.

Methods
All results in this paper were derived from simulations of highly reduced ordinary differential equation models for the excitatory kernel of the preBötC and for a respiratory CPG network that includes this preBötC component as well as inhibitory early-I, post-I, and aug-E components and, in some simulations, an additional excitatory PiCo component (Fig 1). The models were based on those used in several previous respiratory modeling studies [13,38,40,54]. In this framework, each neuronal population is assumed to undergo synchronized transitions between active and silent states, but without spike-level synchrony. Thus, each is represented by a single non-spiking neuron model unit, the voltage of which corresponds to the average of the membrane potentials of the neurons in that population.
For purposes of model specification, we number the model units as follows: 1-pre-I, 2early-I, 3-post-I, 4-aug-E, 5-PiCo. All units obeyed the equations dV i =dt ¼ ðÀ I NaP;i À I K;i À I L;i À I synÀ I;i À I synÀ E;i Þ=C; where i 2 {1, . . ., 5}. In system (1), I NaP,i represents the current through persistent sodium channels, I K,i is a potassium current, I L,i = g L,α � (V i − E L,α ) is the leakage current with α = exc for i 2 {1, 5} and α = inh for i 2 {2, 3, 4}, respectively. Moreover, I syn−I,i and I syn−E,i are inhibitory and excitatory synaptic currents, respectively.
For the excitatory units, with i 2 {1, 5}, the intrinsic currents in (1) take the following forms: For the inhibitory units, with i 2 {2, 3, 4}, the intrinsic currents in (1) are given by: with That is, the persistent sodium current takes the same form for all units but has a stronger conductance for excitatory than inhibitory units. On the other hand, as in previous models (e.g., [13]), the potassium current for the excitatory units is a standard delayed rectifier that is set to be quite weak, representing potassium current activated by the depolarization attributable to the persistent sodium current, while for the inhibitory units it is a much stronger adaptation current with a dedicated gating variable, p i . The voltage-dependent activation functions and time constants in Eqs (1), (2) and (3) In Eq (4) for the adaptation current gating variable, the function f out,inh (V) appears. This function represents a nonlinear filtering of the output of a unit and takes the sigmoidal form f out;inh ¼ 1=ð1 þ exp ½ðV À y out;inh Þ=s out;inh �Þ: ð6Þ Thus, Eq (4) specifies how the conductance of an inhibitory unit's adaptation current grows while that unit is active.
To define the synaptic currents in Eq (1), we use the function given in Eq (6) along with similar sigmoidal functions f out;preÀ I ¼ 1=ð1 þ exp ½ðV À y out;preÀ I Þ=s out;preÀ I �Þ for pre À I; f out;pico ¼ 1=ð1 þ exp ½ðV À y out;pico Þ=s out;pico �Þ for PiCo: The synaptic currents to unit i are given by I synÀ I;i ðV i Þ ¼ g synÀ I � ðV i À E synÀ I Þ � P j2N inh i b ji f out;inh ðV j Þ; i a ji f out;exc ðV j Þ þ where N inh i ; N exc i denote the sets of inhibitory and excitatory units, respectively, presynaptic to unit i, the constants b ji , a ji give the strengths of inhibitory and excitatory connections, and the terms c ji , j 2 {1, 2}, are parameters representing the strength of the tonic drive to each unit, with c 1i tunable and c 2i fixed for each i.
Some figures in this paper include nullclines or bifurcation diagrams. Note that system (1) for a single unit in the excitatory case (i.e., without a p variable) consists of a pair of coupled differential equations with dependent variables V, h, subject to time-dependent inputs from the synaptic currents I syn−I , I syn−E . If these inputs are held fixed or are set to zero to represent a unit in isolation, then the behavior of solutions of the system can be visualized in the (V, h) plane. Within this plane, the set of (V, h) for which dV/dt = 0 holds is the V-nullcline, with an analogous definition for the h-nullcline. For system (1), each of these sets is a curve (Fig 2A). In particular, for a certain range of parameters, the V-nullcline is a cubic-shaped curve with two folds, or knees, the left knee (LK) and right knee (RK).
The position of the each nullcline will vary with changes in model parameters that appear in the corresponding ODE. Intersection points of the nullclines are equilibria of the system; if (V, h) are initially chosen to correspond to their values at an equilibrium point, then (V, h) will remain constant for all time. For sufficiently small c 11 , the equilibrium point (or fixed point, FP) lies on the left branch of the V-nullcline. Yet since the V-nullcline for the pre-I unit moves under modulation of c 11 (Fig 2A), the position of the unit's equilibrium point also changes. We use XPPAUT to generate a bifurcation diagram for the pre-I unit by tracking both the location and the stability of the equilibrium point (Fig 2B), along with information about periodic orbits that emerge (disappear) when the equilibrium loses (gains) stability, as c 11 is increased.
As c 11 is increased, the c 11 value where the pre-I equilibrium loses stability represents the transition from quiescent to oscillatory intrinsic dynamics for our baseline parameter set, while the c 11 value where the equilibrium regains stability corresponds to the transition from oscillatory to tonic intrinsic dynamics. If any other parameter from the pre-I system is altered, then this variation may change these transitional c 11 values. Thus, for each transition, we can use XPPAUT to obtain a curve of values in the two-dimensional (c 11 , g NaP,exc ) parameter space where this transition occurs; the curve corresponding to the oscillatory-to-tonic transition appears in each panel of Fig 5. In the description above, we set the inhibition level to the pre-I unit equal to 0. We can also treat that inhibition level as a parameter of the pre-I model. Variation in the inhibition level will alter the position of the V-nullcline and hence of the LK and FP. We also use XPPAUT to compute how these features depend on the inhibition level, with results appearing in Figs 6, 7, 12.