A Positive Feedback at the Cellular Level Promotes Robustness and Modulation at the Circuit Level

Julie Dethier, Guillaume Drion, Alessio Franci, Rodolphe Sepulchre

1 Department of Electrical Engineering and Computer Science, University of Liège, Liège B-4000, Belgium

2 Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ 08544, USA

3 Laboratory of Pharmacology and GIGA Neurosciences, University of Liège, Liège B-4000, Belgium

4 Volen Center for Complex Systems, Brandeis University, Waltham, MA 02454, USA

5 Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, UK

* These authors contributed equally to this work

## Abstract

The paper highlights the role of a positive feedback gating mechanism at the cellular level in the robustness and modulation properties of rhythmic activities at the circuit level. The results are presented in the context of half-center oscillators, which are simple rhythmic circuits composed of two reciprocally connected inhibitory neuronal populations. Specifically, we focus on rhythms that rely on a particular excitability property, the post-inhibitory rebound, an intrinsic cellular property that elicits transient membrane depolarization when released from hyperpolarization. Two distinct ionic currents can evoke this transient depolarization: a hyperpolarization-activated cation current and a low-threshold T-type calcium current. The presence of a slow activation is specific to the T-type calcium current and provides a slow-positive feedback at the cellular level that is absent in the cation current. We show that this slow-positive feedback is necessary and sufficient to endow the network rhythm with physiological modulation and robustness properties. This study thereby identifies an essential cellular property to be retained at the network level in modeling network robustness and modulation.

## Introduction

Biological rhythms play a major role in the functioning of the brain but much remains to be understood regarding their control, regulation, and function. Many advances in this important question have come from experimental and computational studies of central pattern generators (CPGs), which endogenously produce precise rhythmic outputs directly related to motor functions [1, 2, 3, 4, 5]. In this effort, experimental work benefits from computational models but models at the circuit level usually rely on mathematical simplifications at the component level. The question of which cellular details must be retained at the network level is largely open [6].

Motivated by this general question, we highlight a simple feedback mechanism at the cellular level that has a key influence on circuit robustness and modulation. We illustrate this property via the computational study of an archetype model of CPG circuits, the half-center oscillator (HCO): two neuronal populations that do not oscillate in isolation, but oscillate in an antiphase rhythm when reciprocally connected [7, 8, 9]. Because of the widespread occurrence of this circuit motif, the mechanisms have been extensively studied, both computationally and experimentally [10, 3, 5, 11, 12, 13, 14]. A specific cellular excitability property, the post-inhibitory rebound (PIR) [8], and its two specific ionic currents, and [15], are well-known key players in circuit oscillations.

Previous studies have focused on distinguishing those two currents from their contribution to the escape or release mechanism [13, 12, 14] and their modulation of rhythmic activity [16, 17, 18, 19]. This paper highlights that those two currents differ in another simple yet fundamental aspect: : both generate a PIR, but only one of them acts, through its slow activation, as a source of positive feedback in the network rhythm timescale. The present paper demonstrates through a computational study that this particular feedback is fundamental for the robustness and modulation properties of the circuit rhythm, and that its absence is detrimental both to robustness and modulation at the circuit level.

Our results predict that PIR per se is not a sufficient cellular excitability property to be retained at the network level. In addition, the current slow-regenerative regenerative (slow-positive-feedback) nature must be retained as an important dynamical parameter. We emphasized in previous work the importance of this positive feedback at the cellular level in defining bursting excitability [20] and its widespread regulation in different neuron types [21]. The present paper prolongs this work, moving from the role of regenerativity at the cellular level to its importance at the circuit level.

## Results

### Slow activation of T-type calcium channels is critical to robustness of network rhythmic activity

To assess the role of cellular properties in network rhythms, we consider one of the simplest and most studied networks: the half-center oscillator (HCO). The network rhythm results from the mutual inhibition (I) of two neurons that do not oscillate endogenously in isolation [7, 8, 9]. HCOs have been identified at the core of most endogenous rhythmic circuits, such as CPGs governing locomotion [7, 8, 9, 10, 14] or respiration [1, 22, 14]. Oscillations in HCOs are triggered by an external pulse of hyperpolarizing current. When released from hyperpolarization, the cell generates a burst-like transient depolarization with one or more spikes. This activity hyperpolarizes the other cell via the inhibitory synaptic connection, which in turn triggers a transient burst. The cycle repeats leading to an antiphase rhythm between the two neurons.

The transient depolarization following the termination of an hyperpolarizing input is an essential cellular property for the network rhythm, best known in the literature as post-inhibitory rebound (PIR) [8]. Two major ionic currents have been shown to underlie the PIR (Figure 1A): i) the hyperpolarization-activated cation current, , an hyperpolarization-activated inward current that contributes to rebound responses in a diverse array of neurons in invertebrates and vertebrates [15]; ii) the low-threshold T-type calcium current, , which is deinactivated by hyperpolarization and then activates upon release from inhibition [23]. Many studies have highlighted the distinction between these two currents in HCOs from an “escape or release” mechanism perspective, T-type calcium currents inducing the release mechanism and -like currents promoting the escape mechanism [13, 12, 14]): ): either the active cell “releases” its inhibitory effect on the silent cell (release mechanism), or the silent cell “escapes” from inhibition via the activation of an -like current (escape mechanism).

Separately, the two types of current generate similar PIR traces in single cells (Figure 1A; see Methods for cellular models and model difference between Mechanism A and B). While both mechanisms are redundant for the generation of oscillations in a network with reciprocal inhibition, we emphasize a fundamental difference between the two: in presence of physiological variability, i.e., variability in the intrinsic cellular properties and synaptic connections (see Methods for a description of variability), only the rhythm generated by Mechanism B is robust (Figure 1B).This robustness property highlights a fundamental difference between the two mechanisms.

This difference lies in the dynamical feedback loops generated by the gating variables of the two currents. Both currents generate an ultraslow (i.e., that lasts over the course of several action potentials) inward current in response to hyperpolarization, which is the foundation of the PIR. This inward current counteracts the external hyperpolarization and acts as a source of negative feedback on membrane potential variations—or restorativity in the terminology of [21]—in the ultraslow timescale (Figure 2, red feedback loops). But, in contrast to , provides a slow positive feedback on membrane potential variations—or regenerativity in the terminology of [21]—via its slow activation variable (Figure 2, green feedback loop). This slow positive feedback is absent in Mechanism A. At the cellular level, the slow positive feedback is revealed by a specific signature during hyperpolarization [21]: a transient excitatory pulse that triggers a single spike in Mechanism A (), triggers a burst in Mechanism B () (Figure 2, bottom panel). This signature reveals that bursts are endogenously generated with a PIR with slow regenerativity (Mechanism B), as opposed to a purely restorative—i.e., only —PIR (Mechanism A).

A frequent modeling simplification is to neglect the slow activation kinetics of T-type calcium channels and to consider the activation at steady-state (i.e., instantaneous). It should be noted that the slow regenerativity is lost in this approximation, which eliminates the bursting signature observed in Mechanism B, as illustrated in Figure 2, center panel.

In the rest of the paper, we investigate the impact of the difference between Mechanism A and B at the network level. More specifically, we look at a few simple quantities that characterize the network rhythm: the network—or interburst—frequency, i.e., the inverse of the time duration between two burst onsets, the duty cycle, which is the ratio between a burst duration and the time duration between two bursts, and the ratio between the duty cycle in neuron 1 and in neuron 2. In order to isolate the contribution of the slow regenerativity only, we compare the two mechanisms in a model that only includes, , , and one “PIR current”, , modeled by T-type calcium channels with instantaneous activation for Mechanism A, designated simply by “PIR”f or the rest of the paper, and modeled by T-type calcium channels with slow activation for Mechanism B, designated by “PIR + slow regenerativity” (see Methods for cellular models and details of the simulations). Accordingly, the cellular models contain an identical PIR current in the two cases except for the activation time constant of the PIR current and both neuron models possess identical I/V curves. Results do not depend on other properties such as the role of the sag brought by or the difference between release and escape, both models differing only in their slow regenerativity. We stress that all results obtained under Mechanism A can be reproduced in a model where the PIR is modeled by channels only.

### Robustness of network oscillations requires PIR with slow regenerativity

There exists extensive experimental evidence that the rhythmic activity of neuronal circuits is robust against variability in intrinsic parameters (such as ionic conductances across neurons), extrinsic parameters (such as synaptic conductances), and exogenous noise (such as synaptic currents external to the circuit) [24, 25, 26, 27, 28]. We tested the robustness of HCOs in a network with PIR without slow regenerativity against a network with PIR with slow regenerativity (see Methods for network description and details of the simulations). The results show the drastic influence of cellular slow regenerativity in the robustness of the network.

Intrinsic variability of the network was studied by introducing variability (see Methods for a description of variability) in the maximal conductance of the PIR current, , in a network with two populations (Figure 3, network connections). Variability in the cellular properties dramatically impacts the rhythmic activity of the network without slow regenerativity (Figure 3, left panel). The network rhythm becomes unstable beyond 75% of variability and is significantly perturbed for smaller values. In sharp contrast, the network oscillations with slow regenerativity are robust against intrinsic variability up to 200% (Figure 3, right panel). Remarkably, the network frequency is almost unaffected by the intrinsic variability, a consequence of the positive feedback brought by slow regenerative currents. Instead, the network frequency is strongly affected by intrinsic variability without slow-regenerative currents.

The robustness of the network oscillations against variability in extrinsic parameters was studied by varying the maximal synaptic conductance parameters, , in a two-neuron network with reciprocal connections (see Methods for a description of variability). Without slow regenerativity, a small variability in the synaptic conductances affects dramatically the network activity (Figure 4, left panel): identical maximal synaptic conductances generate oscillations but oscillations become unstable when the maximal synaptic conductances differ between the two cells. Oscillations with a PIR without slow regenerativity are fragile to network variability. In sharp contrast, variability in the synaptic conductances is possible for a much larger range with slow regenerativity and the network frequency is also almost independent of the synaptic variability (Figure 4, right panel). Oscillations persist up to a variability higher than 80%. A source of slow-positive feedback in the PIR mechanism is therefore essential to robustness of network oscillations against network variability.

The robustness of the network oscillations against exogenous disturbances was investigated by adding a Gaussian white noise in the equation that models membrane potential variations (see Methods for a description of noise). This emulates the external perturbations—spike train inputs from surrounding neurons—received by a network when studied in a noisy environment rather than in isolation [29]. We simulated a sixteen-neuron network with two populations, with a different noise source for each neuron (Figure 5).

The results are consistent with the robustness against parameter variability. Without slow regenerativity, oscillations are sensitive to noise and completely disappear with a noise level greater than 0.15. With slow regenerativity, oscillations are robust to noise up to a level of 0.225. Similarly to the introduction of variability in intrinsic and extrinsic parameters, the network frequency is also less affected with in the presence of slow regenerativity.

### Robust modulation of network properties requires slow regenerativity

Neuromodulators can tune and reconfigure the network dynamics, affecting both the frequency and phasing of neurons [30, 31]. For instance, in the Tritonia swim CPG, intrinsic modulation can produce an enhanced level of excitability, triggering circuit activity that is maintained after the initial signal, generating an escape swimming response to particular aversive stimulus [32]. Neuromodulation can also switch the circuit between rhythms: in the crustacean stomatogastric ganglion, neuromodulators can switch the circuit activity from the fast pyloric rhythm, to two slower rhythms, the gastric mill rhythm and the cardiac sac rhythm [33]. In addition, neuromodulators can determine the active neuronal elements in the circuit or combine elements from different circuits into one [33, 5].

Experimentally, the network properties, i.e., network frequency and duty cycle—or phase relation—can be modulated via both intrinsic neuron parameters and synaptic parameters on multiple timescales [34, 35, 36, 31, 37, 5, 33, 38, 39]. In this section, we investigate how the network rhythmic activity responded to these modulations, both with PIR without slow regenerativity and with PIR with slow regenerativity (see Methods for network description and details of the simulations).

Extrinsic parameters, i.e., the synaptic parameters and , given intrinsic (cellular) characteristics, modulate the network frequency. Synaptic coupling is very plastic [40, 41] and synapses are a primary target of modulators [39]. Synaptic currents can be generated by the cooperation of several ion channel subtypes which can have slightly different kinetics. Variation of the synaptic parameters results from a variation of the contribution of all the subtypes. Absolute variation of the different ion channels influences the maximal conductance whereas their relative variation can modulate the time constant of the synaptic current that aggregates all the different subtypes in a model. Therefore, both the synaptic magnitude, , and the synaptic kinetics, , can be sources of modulation.

Oscillations with cellular slow regenerativity can be modulated over a large range by extrinsic parameters (Figure 6, right panal). Variation of the and parameters generates a 150% increase in network frequency (see Methods for a description of mean frequency). Such a span of modulation is observed in physiological rhythms: for instance, there is a 150% increase in frequency from slow-wave sleep () to sleep spindles () and a 250% increase from beta-band oscillations () to gamma-band oscillations (). In addition, with slow regenerativity, the network frequency is only weakly sensitive to the variability in parameters as shown with the standard deviation plot and the highly similar network frequencies in the membrane voltage traces (see Methods for a description of variability and standard deviation). In opposition, variations of without slow regenerativity allow for network frequency modulation for a much smaller parameter range (Figure 6, left panel). Moreover, this modulation is very fragile and very sensitive to variability: the standard deviation reaches higher values than with slow regenerativity and membrane potential traces, for a same set of parameters but different simulations, are drastically different (Figure 6, bottom left panel). Variation of is almost impossible: must lie in a very specific timescale for the oscillations to develop in the network. The modulation requires a tight coupling between intrinsic and extrinsic parameters: the network oscillations are a direct reflection of the unicellular activity. The oscillation frequency is set by the neuron intrinsic dynamics and almost no variation can be induced by the synaptic dynamics.

Intrinsic parameters, i.e., the cellular parameters and , given extrinsic (synaptic) characteristics, modulate the duty cycle and duty cycle ratio (see Methods for a description of duty cycle and duty cycle ratio). Many neuromodulators act on the neuron intrinsic properties by altering the balance of conductances, modifying their excitability properties [31]. The maximal conductance of the PIR current is a natural candidate for modulation by intrinsic parameters.

The high robustness brought by cellular slow regenerativity allows also for the modulation by intrinsic parameters even in presence of variability in the network (Figure 7; see Methods for a description of variability). Covariation of the maximal PIR conductances, and , leads to an increase in duty cycle ratio of 150% (Figure 7, top right panel). Independent variation of the same parameters, i.e., varying and independently, modulates the duty cycle ratio up to a factor two (Figure 7, bottom right panel). Variation in phase relation have been observed for instance in cats, during normal locomotion, where the shortening, by a factor two or three, of one of the phase (the extensor phase) leads to faster walking [42]. In contrast, our computational model suggests that modulation with PIR without slow regenerativity is so fragile that it is unrealistic. Stable oscillations with variation of intrinsic parameters do not cover a large parameter range (Figure 7, left panel).

In brief, the high robustness brought by cellular slow regenerativity allows for the modulation by both extrinsic and intrinsic parameters. In addition, the network frequency and duty cycle can be modulated independently and in presence of physiological variability in the network.

## Discussion

### Cellular slow regenerativity is essential to robustness and modulation of network rhythmic activity

The main message of this paper is to highlight the role of slow regenerativity, a cellular excitability property, in endowing network oscillations with robustness and modulation properties that seem ubiquitous in physiological neuronal networks. An ionic current is slowly regenerative if it provides a source of positive feedback around resting potential in the slow timescale of repolarization [21]. The importance of this cellular property was assessed in one of the simplest and best understood network oscillation mechanisms, the antiphasic rhythm observed between two populations of neurons reciprocally connected by inhibitory synaptic connections. Many earlier studies have emphasized the role of post-inhibitory rebound (PIR) at the cellular level as a core mechanism for the network oscillation, and have identified and as two distinct ionic currents that can participate in the PIR. Our novel contribution is to observe that the cellular PIR will enable a robust and subject to modulation network oscillation only in the presence of a slow-regenerative ionic current. Because both and are sources of PIR currents but only is slow regenerative, our paper suggests a novel and somewhat fundamental complementarity between T-type calcium and channels in PIR mechanisms.

As a source of positive feedback, regenerative currents make the PIR endogenous, that is, robust to intrinsic and extrinsic sources of variability. As a consequence, a PIR with slow-regenerative currents allows for network oscillations that are robust and subject to modulation. The network oscillation is robust because it can sustain large variability across the neuronal population both in intrinsic (cellular) and extrinsic (synaptic) parameters. It is also subject to modulation because the frequency and phase properties of the oscillation can be controlled over a broad range by a relative modulation of extrinsic or intrinsic conductances. Our computational investigation illustrated that this robustness and modulation properties are lost when the PIR is purely ultraslow restorative.

### Complementarity between the two types of PIR currents

In the context of HCOs, many neurons possess both and , the two main currents that contribute to PIR [10, 43, 44, 16, 17, 18, 19]. When those two currents are present, both and can be a source of modulation. In this case, the presence of T-type calcium currents as a source of slow regenerativity is sufficient to guarantee network oscillation robustness. On the other hand, the hyperpolarization-activated cation current can modulate drastically the network frequency and duty cycle [16, 17]. However, our computational model suggests that this is only the case if a slow-regenerative current, and therefore cellular endogenous characteristics, is supplied by another mechanism. It is noteworthy that the necessity for slow restorativity can be achieved by other means, such as the presence of high-threshold calcium channels. This necessary condition for slow regenerativity reveals a somewhat fundamental complementarity, distinct from the release or escape view, between the two channels: allows for stable rhythmic oscillations to emerge and enlarges the modulation possibilities.

### Positive feedback as a source of endogenous activity

Slow regenerativity is nothing but a source of positive feedback in the slow timescale of repolarization. It is a slow analog of the positive feedback brought by sodium activation in the fast timescale of spike upstroke. In previous work [21], we showed that this positive feedback is essential for the robust coexistence of hyperpolarized and spiking states at the cellular level. We subsequently showed in [20] that this positive feedback is essential for modulation and robustness of bursting activities. Here we show that the same positive feedback at the cellular level is also essential for robustness and modulation at the network level. The common feature of the positive feedback in those three phenomena is that it makes the neuronal excitability in the slow timescale an endogenous property, robust to intrinsic and extrinsic variability. Making an activity endogenous is the very nature of positive feedback and has been emphasized in a number of contexts. The results presented in this paper are in line with the discussion of the role of positive feedback in other biological models, such as for instance the biochemical mechanisms underlying the mitotic oscillator [45, 46, 47]: the oscillator is endogenous and robust in the presence of positive feedback, whereas it becomes exogenous and entrainable when the source of positive feedback disappears. The importance at the network level of positive feedback at the cellular level is thought to be general and not specific to the case study of HCOs chosen in this paper for its simplicity and physiological relevance.

### Slow regenerativity in half-center oscillator models

There is a rich literature on computational models of oscillations generated by reciprocal inhibition. HCOs have been used to model rhythmic motor outputs in many invertebrates and vertebrates [7, 8, 9, 10, 3, 5]. In a different context, models of thalamocortical spindle oscillations suggest that the rhythm originates from the thalamic reticular nucleus, which consists in interacting inhibitory nonoscillatory neurons [13, 48, 49, 50].

It is of interest to observe the varying degree of cellular regenerativity in published models of HCOs. Early models are conductance-based and usually include at the cellular level both and , the two main physiological currents eliciting the PIR [51, 10, 43, 44, 49]. However, network computational studies often lead to a subsequent mathematical simplifications of the cellular details and the cellular slow-positive feedback is often lost in this reduction process. A frequent simplification in the literature (see e.g. [13, 48, 52, 53, 14]) is to resort to a steady-state approximation of the calcium activation in the same way as it is normally done for sodium activation. But this approximation rests on neglecting fast dynamics, which amounts to consider calcium channels as a source of fast rather than slow positive feedback. The resulting reduced models have therefore lost their source of slow regenerativity, which makes them unsuitable for robustness and modulation studies at the network level.

The alternative model reduction consists to model the cellular level as Morris-Lecar type of neurons, retaining the slow calcium currents but neglecting the fast sodium currents [11, 12]. Those models do retain the slow-positive feedback source necessary for robustness but they lose the modulation capabilities illustrated in the present paper because the network interconnection properties are spike-dependent. This prevents exogenous modulation of the rhythm. In addition, if sodium spikes were added to a Morris-Lecar neuron with the addition of the spike currents (as suggested in [54]) while keeping the calcium activation at steady-state, the slow regenerativity would be destroyed.

It should be highlighted that it is possible to derive reduced neuronal models that do retain the balance of slow positive and negative feedbacks as an explicit parameter, see e.g. the recent models in [55, 20]. The results of the present paper suggest that it is an important feature to retain in a simplified model aimed at network computational studies.

## Methods

All the numerical simulations and analyses were performed with MATLAB, MathWorks. The models were implemented in a MATLAB code and simulated using a forward Euler method with a time step of .

### Cellular model

The cellular model is inspired from the crab stomatogastric ganglion (STG) conductance-based neuron model [56, 24]. The model contains the standard Hodgkin-Huxley (HH) currents (Hodgkin and Huxley, 1952): the transient sodium current, , a fast depolarizing current, and the delayed-rectifier potassium current, , a slower hyperpolarizing current, plus a leak current, . The two currents responsible for PIR are: a low threshold T-type calcium current, , and a hyperpolarization activated cation current, . The membrane potential dynamics writes as follows:

where is the membrane capacitance and is the applied current.

Each ionic current takes the standard HH form:

where is the maximal conductance for current , and are integers, and is the reversal potential of the ion . Table 1 lists the values of , , , and for the different currents.

60 | 3 | 1 | 50 | |

40 | 4 | 0 | -70 | |

0.035 | 0 | 0 | -49 | |

0.3 | 3 | 1 | 120 | |

0.04 | 1 | 0 | -20 |

Notation is explained in the text. All conductances are in and membrane potentials in .

Activation and inactivation variable dynamics follow the classical formalism:

where the functions for , , , and are given in Table 2. Note that, if and/or , and and/or and are not listed, respectively.

Notation is explained in the text. All time constants are in .

In Figures 3-7, the only current responsible for the PIR is , either with instantaneous activation, i.e., or with slow activation, i.e., of Table 2 (). The former case corresponds to Mechanism A or “PIR”, whereas the latter entails Mechanism B or “PIR + slow regenerativity” (see text for details). In both cases, we denote this current and its maximal conductance as and .

### Network model

The inhibitory synaptic connections are gamma-aminobutyric acid (GABA) and are made with exponential synapses of the GABA type. The synaptic current between two neurons takes the form [50]:

(4) | ||||

where is the presynaptic membrane potential and the number of presynaptic neurons. If not stated otherwise, , , , and , , .

### Variability

Physiological variability is modeled by randomly selecting the values for the parameter subjected to variability in an interval, called variation range, centered on a given parameter value. The variability level quantifies the width of this interval, in percentage of the given value. For instance, a 200% variability in means that, for a given parameter value of , the variation range has a width of and is centered on . The parameters are therefore randomly selected out of the interval .

### Analyses

A network is categorized as having a stable rhythmic activity (rhythm ON) if all the neurons in the network are still bursting in the stationary state. Due to the specific network structure—connection all-to-all from one population to the other, and vice versa—all the neurons in one population receive the same input (coming from all the neurons in the other population). Therefore, if all the neurons are bursting, the bursts have been elicited by the same transient hyperpolarizing input—bursting cannot happen without this hyperpolarization—and the bursts are synchronous, i.e., all the bursts overlap but not necessarily the spikes. This provokes HCO antiphasic oscillations. Bursting in neurons is detected by having two consecutive spikes less than appart. Practically, we detect busting in all neurons after the transient phase: due to the time constants in play, analyzing the last of data was sufficient to have bursting in the two populations in the stationary state. If at least one neuron was not bursting during this time period, the network is categorized as having no stable rhythmic activity (rhythm OFF).

The frequency is the inverse of the time duration between the beginning of two bursts, or period, averaged over the two neurons. The duty cycle is the ratio between the burst duration and the period, averaged over the two neurons. The duty cycle ratio is the ratio between the duty cycle in neuron 1 and in neuron 2. Each value for the mean frequency, mean duty cycle and mean duty cycle ratio was computed from 10 simulations with the same set of parameters but with 40% variability in and 20% variability in . Only simulations endowed with rhythmic activity in the sense defined previously were considered to compute the means and standard deviations. If no rhythmic activity was detected, the computed means and standard deviations were set to 0. The proportion of oscillatory HCO quantifies the percentage of simulations that showed rhythmic activity out of the 10 simulations.

### Simulation details

Cell simulations of Figure 1, are performed using the cellular model described above, with for Mechanism A and for Mechanism B. The applied current, , on both neurons takes a value of . During the hyperpolarization, the applied current drops to .

The network simulations of Figure 1 were obtained with the cellular model with , and with instantaneous activation for the left panel, i.e., from Table 2, and slow activation with time constant (see Table 2) for the right panel. The top panel does not present any parameter variability, the two neurons and synaptic connections are identical. Physiological variability is simulated with 40% variability in and 20% variability in .

The cell models in Figure 2 were performed using the cellular model described above, with in the left panel, and instantaneous activation in the center panel, and and slow activation in the right panel. The applied current, , takes a value of . During the hyperpolarization, the applied current drops to and takes a value of for for the fast depolarizing input.

The two populations of Figure 3, each composed of 8 neurons, are connected all-to-all with the synaptic current described in the previous section. The intrinsic neuron parameter is with a variability level from 0% to a 200%.

The simulations of Figure 4 were done for a two-neuron network. The extrinsic neuron parameter is with a variability level from 0% to a 100%.

The two populations considered in Figure 5, each composed of 8 neurons, are connected all-to-all with the synaptic current described in the previous section. A Gaussian white noise is added in the voltage equation to model the typical spike train input received from the many other unmodeled neurons [29]. The noise is modeled by , where is the noise intensity and varies from 0 to 0.25, and is drawn from a normal distribution with zero mean and unitary standard deviation and is different for each neuron.

In Figure 6, the maximal synaptic conductance, (from to ), and the synaptic time constant, proportional to ( varies from to ), vary simultaneously for the two neurons. The variability level is 40% for and 20% for . Membrane potential plots ([], []): top left panel–(5, 0.05) and (0.5, 0.05); top right panel–(5, 0.05) and (6.5, 0.5); bottom panel–(6, 0.2) with variability.

In Figure 7, the maximal PIR conductances, and (from to ), vary independently for the two neurons. The variability level is 40% for and 20% for . The zoom for the duty cycle and duty cycle ratio covers a range from to .

## References

- 1. Grillner S, Wallén P (1985) Central pattern generators for locomotion, with special reference to vertebrates. Annu Rev Neurosci 8: 233-61.
- 2. Getting PA (1989) Emerging principles governing the operation of neural networks. Annu Rev Neurosci 12: 185-204.
- 3. Marder E, Calabrese RL (1996) Principles of rhythmic motor pattern generation. Physiol Rev 76: 687-717.
- 4. Marder E, Bucher D (2001) Central pattern generators and the control of rhythmic movements. Curr Biol 11: R986-96.
- 5. Harris-Warrick RM (2011) Neuromodulation and flexibility in central pattern generator networks. Curr Opin Neurobiol 21: 685-92.
- 6. Marder E, Goaillard JM (2006) Variability, compensation and homeostasis in neuron and network function. Nat Rev Neurosci 7: 563-74.
- 7. Brown TG (1911) The intrinsic factors in the act of progression in the mammal. Proc R Soc Lond B 84: 308-319.
- 8. Perkel DH, Mulloney B (1974) Motor pattern production in reciprocally inhibitory neurons exhibiting postinhibitory rebound. Science 185: 181-3.
- 9. Satterlie RA (1985) Reciprocal inhibition and postinhibitory rebound produce reverberation in a locomotor pattern generator. Science 229: 402-4.
- 10. Calabrese RL, De Schutter E (1992) Motor-pattern-generating networks in invertebrates: modeling our way toward understanding. Trends Neurosci 15: 439-45.
- 11. Skinner FK, Turrigiano GG, Marder E (1993) Frequency and burst duration in oscillating neurons and two-cell networks. Biol Cybern 69: 375-83.
- 12. Skinner FK, Kopell N, Marder E (1994) Mechanisms for oscillation and frequency control in reciprocally inhibitory model neural networks. J Comput Neurosci 1: 69-87.
- 13. Wang XJ, Rinzel J (1992) Alternating and synchronous rhythms in reciprocally inhibitory model neurons. Neural Comput 4: 84-97.
- 14. Daun S, Rubin JE, Rybak IA (2009) Control of oscillation periods and phase durations in half-center central pattern generators: a comparative mechanistic analysis. J Comput Neurosci 27: 3-36.
- 15. Angstadt JD, Grassmann JL, Theriault KM, Levasseur SM (2005) Mechanisms of postinhibitory rebound and its modulation by serotonin in excitatory swim motor neurons of the medicinal leech. J Comp Physiol A Neuroethol Sens Neural Behav Physiol 191: 715-32.
- 16. Cymbalyuk GS, Gaudry Q, Masino MA, Calabrese RL (2002) Bursting in leech heart interneurons: cell-autonomous and network-based mechanisms. J Neurosci 22: 10580-92.
- 17. Sorensen M, DeWeerth S, Cymbalyuk G, Calabrese RL (2004) Using a hybrid neural system to reveal regulation of neuronal network activity by an intrinsic current. J Neurosci 24: 5427-38.
- 18. Olypher A, Cymbalyuk G, Calabrese RL (2006) Hybrid systems analysis of the control of burst duration by low-voltage-activated calcium current in leech heart interneurons. J Neurophysiol 96: 2857-67.
- 19. Doloc-Mihu A, Calabrese RL (2011) A database of computational models of a half-center oscillator for analyzing how neuronal parameters influence network activity. J Biol Phys 37: 263-83.
- 20. Franci A, Drion G, Sepulchre R (2014) Modeling the modulation of neuronal bursting: a singularity theory approach. SIAM J Appl Dyn Syst 13: 798-829.
- 21. Franci A, Drion G, Seutin V, Sepulchre R (2013) A balance equation determines a switch in neuronal excitability. PLoS Comput Biol 9: e1003040.
- 22. Butera RJ Jr, Rinzel J, Smith JC (1999) Models of respiratory rhythm generation in the pre-bötzinger complex. ii. populations of coupled pacemaker neurons. J Neurophysiol 82: 398-415.
- 23. Steriade M, Gloor P, Llinás RR, Lopes da Silva FH, Mesulam MM (1990) Basic mechanisms of cerebral rhythmic activities. Electroen Clin Neuro 76: 481-508.
- 24. Liu Z, Golowasch J, Marder E, Abbott LF (1998) A model neuron with activity-dependent conductances regulated by multiple calcium sensors. J Neurosci 18: 2309-20.
- 25. Golowasch J, Abbott LF, Marder E (1999) Activity-dependent regulation of potassium currents in an identified neuron of the stomatogastric ganglion of the crab Cancer borealis. J Neurosci 19: RC33+.
- 26. Golowasch J, Casey M, Abbott LF, Marder E (1999) Network stability from activity-dependent regulation of neuronal conductances. Neural Comput 11: 1079-96.
- 27. Goldman MS, Golowasch J, Marder E, Abbott LF (2001) Global structure, robustness, and modulation of neuronal models. J Neurosci 21: 5229-38.
- 28. Destexhe A, Rudolph-Lilith M (2012) Neuronal noise. In: Destexhe A, Brette R, editors, Springer Series in Computational Neuroscience, Springer US, volume 8.
- 29. Lindner B, Longtin A, Bulsara AR (2003) Analytic expressions for rate and CV of a type I neuron driven by white gaussian noise. Neural Comput 15: 1761-88.
- 30. Marder E, Weimann JM (1992) Modulatory control of multiple task processing in the stomatogastric nervous system. In: Kien J, Mccrohan C, Winlow B, editors, Neurobiology of motor progamme selection, New York: Pergamon Press. pp. 3-19.
- 31. Marder E, Thirumalai V (2002) Cellular, synaptic and network effects of neuromodulation. Neural Netw 15: 479-93.
- 32. Katz PS, Getting PA, Frost WN (1994) Dynamic neuromodulation of synaptic strength intrinsic to a central pattern generator circuit. Nature 367: 729-31.
- 33. Marder E (2012) Neuromodulation of neuronal circuits: back to the future. Neuron 76: 1-11.
- 34. Harris-Warrick RM, Marder E (1991) Modulation of neural networks for behavior. Annu Rev Neurosci 14: 39-57.
- 35. Fellous JM, Linster C (1998) Computational models of neuromodulation. Neural Comput 10: 771-805.
- 36. El Manira A, Wallén P (2000) Mechanisms of modulation of a neural network. News Physiol Sci 15: 186-91.
- 37. Dickinson PS (2006) Neuromodulation of central pattern generators in invertebrates and vertebrates. Curr Opin Neurobiol 16: 604-14.
- 38. Dayan P (2012) Twenty-five lessons from computational neuromodulation. Neuron 76: 240-56.
- 39. Nadim F, Bucher D (2014) Neuromodulation of neurons and synapses. Curr Opin Neurobiol 29C: 48-56.
- 40. Gerstner W, Kistler W (2002) Spiking neuron models: Single neurons, populations, plasticity. Cambridge: Cambridge University Press.
- 41. Dayan P, Abbott LF (2005) Theoretical neuroscience: Computational and mathematical modeling of neural systems. Cambridge, Mass: The MIT Press.
- 42. Halbertsma JM (1983) The stride cycle of the cat: the modelling of locomotion by computerized analysis of automatic recordings. Acta Physiol Scand Suppl 521: 1-75.
- 43. De Schutter E, Angstadt JD, Calabrese RL (1993) A model of graded synaptic transmission for use in dynamic network simulations. J Neurophysiol 69: 1225-35.
- 44. Nadim F, Olsen OH, De Schutter E, Calabrese RL (1995) Modeling the leech heartbeat elemental oscillator. i. interactions of intrinsic and synaptic currents. J Comput Neurosci 2: 215-35.
- 45. Novak B, Tyson JJ (1993) Numerical analysis of a comprehensive model of M-phase control in Xenopus oocyte extracts and intact embryos. J Cell Sci 106: 1153-68.
- 46. Pomerening JR, Sontag ED, Ferrell JE Jr (2003) Building a cell cycle oscillator: hysteresis and bistability in the activation of Cdc2. Nat Cell Biol 5: 346-51.
- 47. Tyson JJ, Novak B (2008) Temporal organization of the cell cycle. Curr Biol 18: R759-68.
- 48. Wang XJ, Rinzel J (1993) Spindle rhythmicity in the reticularis thalami nucleus: synchronization among mutually inhibitory neurons. Neuroscience 53: 899-904.
- 49. Destexhe A, Contreras D, Sejnowski TJ, Steriade M (1994) A model of spindle rhythmicity in the isolated thalamic reticular nucleus. J Neurophysiol 72: 803-18.
- 50. Golomb D, Wang XJ, Rinzel J (1994) Synchronization properties of spindle oscillations in a thalamic reticular nucleus model. J Neurophysiol 72: 1109-26.
- 51. Huguenard JR, Prince DA (1992) A novel T-type current underlies prolonged Ca(2+)-dependent burst firing in gabaergic neurons of rat thalamic reticular nucleus. J Neurosci 12: 3804-17.
- 52. Golomb D, Rinzel J (1993) Dynamics of globally coupled inhibitory neurons with heterogeneity. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 48: 4810-4.
- 53. Rinzel J, Terman D, Wang X, Ermentrout B (1998) Propagating activity patterns in large-scale inhibitory neuronal networks. Science 279: 1351-5.
- 54. Kopell N, Ermentrout G (2002) Mechanisms of phase-locking and frequency control in pairs of coupled neural oscillators. In: Fielder B, Iooss G, Kopell N, editors, Handbook of Dynamical Systems II: Towards Applications, Amsterdam: Elsevier B.V., volume 2, chapter 1. pp. 3-54.
- 55. Drion G, Franci A, Seutin V, Sepulchre R (2012) A novel phase portrait for neuronal excitability. PLoS ONE 7: e41806.
- 56. Turrigiano G, LeMasson G, Marder E (1995) Selective regulation of current densities underlies spontaneous changes in the activity of cultured neurons. J Neurosci 15: 3640-52.