## Abstract

Neurons and circuits each with a distinct balance of intrinsic and synaptic conductances can generate similar behavior but sometimes respond very differently to perturbation. Examining a large family of circuit models with non-identical neurons and synapses underlying rhythmic behavior, we analyzed the circuits' response to modifications in single and multiple intrinsic conductances in the individual neurons. To summarize these changes over the entire range of perturbed parameters, we quantified circuit output by defining a global stability measure. Using this measure, we identified specific subsets of conductances that when perturbed generate similar behavior in diverse individuals of the population. Our unbiased clustering analysis enabled us to quantify circuit stability when simultaneously perturbing multiple conductances as a nonlinear combination of single conductance perturbations. This revealed surprising conductance combinations that can predict the response to specific perturbations, even when the remaining intrinsic and synaptic conductances are unknown. Therefore, our approach can expose hidden variability in the balance of intrinsic and synaptic conductances of the same neurons across different versions of the same circuit solely from the circuit response to perturbations. Developed for a specific family of model circuits, our quantitative approach to characterizing high-dimensional degenerate systems provides a conceptual and analytic framework to guide future theoretical and experimental studies on degeneracy and robustness.

**SIGNIFICANCE STATEMENT** Neural circuits can generate nearly identical behavior despite neuronal and synaptic parameters varying several-fold between individual instantiations. Yet, when these parameters are perturbed through channel deletions and mutations or environmental disturbances, seemingly identical circuits can respond very differently. What distinguishes inconsequential perturbations that barely alter circuit behavior from disruptive perturbations that drastically disturb circuit output remains unclear. Focusing on a family of rhythmic circuits, we propose a computational approach to reveal hidden variability in the intrinsic and synaptic conductances in seemingly identical circuits based solely on circuit output to different perturbations. We uncover specific conductance combinations that work similarly to maintain stability and predict the effect of changing multiple conductances simultaneously, which often results from neuromodulation or injury.

## Introduction

The intrinsic and synaptic conductances in diverse neural circuits, ranging from central pattern generating circuits in invertebrates (Marder and Goaillard, 2006; Schulz et al., 2006; Tobin et al., 2009), to central circuits in mammalian brains (Swensen and Bean, 2005; Günay et al., 2008), can vary several-fold and still generate nearly identical activity patterns. This degeneracy is due to the extensive overlap that ion channels display in their biophysical properties and how they shape neural activity (Marder and Tang, 2010). Although multiple circuits generate similar activity patterns under normal conditions, they can respond very differently when perturbed. One example is the pyloric circuit in the stomatogastric ganglion of crustaceans, which normally generates a stereotypical rhythmic pattern of activity (Tang et al., 2010; Haddad and Marder, 2018). Following extreme temperature perturbations, the circuits in each preparation display so-called “crashes” in different ways, affecting all circuit components (Robertson and Money, 2012; Tang et al., 2012; Rinberg et al., 2013). Yet, in such circuits, some types of perturbations affect individual aspects of circuit output (Ransdell et al., 2013; Sakurai et al., 2014) and can be used to reveal differences in the circuit's intrinsic and synaptic conductances. Beyond rhythmic circuits, degeneracy is also observed during neuropathic pain in primary afferents of nerve-injured rats (Ratté et al., 2014; Ratté and Prescott, 2016).

Given this degeneracy in a circuit's intrinsic and synaptic conductances, many studies have examined how to reliably modulate circuit output while maintaining robustness (Marder et al., 2014). However, experimentally, it is more difficult to record all circuit components over time, compared with a single measure of circuit output. To address this challenge, we used computational modeling where we can read out the values of all biological components and examine how they interact. Unlike previous modeling work, which studied few circuits as representative examples with a desired response (Grashow et al., 2009, 2010; Dethier et al., 2015), we incorporated intrinsic and synaptic biophysical variability in an ensemble modeling approach (Prinz et al., 2004; Lamb and Calabrese, 2013). Specifically, we studied a family of degenerate models of half-center oscillators, small central pattern generating circuits of two neurons coupled with recurrent inhibition (Brown, 1911), a common circuit motif underlying rhythmic behaviors such as breathing and locomotion (Perkel and Mulloney, 1974; Wang and Rinzel, 1992; Sharp et al., 1996). Although modeling studies typically assume identical neurons and synapses, experiments highlight the impact of their variability on circuit output and the circuits' ability to cope with perturbations (Sharp et al., 1996; Grashow et al., 2009, 2010). Therefore, we studied a large population of degenerate half-center circuits composed of non-identical neurons and synapses with properties from the stomatogastric nervous system (Golowasch and Marder, 1992; Liu et al., 1998; Goldman et al., 2001).

We took a novel approach opposite from studies that assumed experimentally observed relationships between the circuits' building blocks. By examining circuit output in response to perturbations, fully or partially blocking or enhancing a given ion channel, we identified surprising conductance combinations that can predict the response to specific perturbations, even when the remaining intrinsic and synaptic conductances are unknown. We first defined a concise measure of circuit stability when altering conductances in biologically plausible ranges based on how significantly circuit output was modified for the entire population of circuits. Unbiased clustering of stability exposed specific subsets of conductances whose perturbation had a similar effect on the output of many degenerate circuits with similar behavior. Unlike classical sensitivity analysis (Goldman et al., 2001; Olypher and Calabrese, 2007; Weaver and Wearne, 2008; Drion et al., 2015), our measure is non-local and goes beyond a linear approximation of the covariation of interacting conductances. We demonstrate this by explaining the stability to perturbations in two conductances simultaneously from the combination of single perturbations, which is relevant for many neuromodulators that alter more than one conductance. Although implemented on a particular family of circuits, our approach provides a broader framework for studying degeneracy and robustness quantitatively, which has the potential to inspire future experimental investigation of these phenomena.

## Materials and Methods

### Single neuron model

For the dynamics of single neurons, we used Hodgkin–Huxley models with seven different channel types (Fig. 1): fast Na^{+} (Na), fast transient K^{+} (A), Ca^{2+}-dependent K^{+} (KCa), a delayed rectifier K^{+} (Kd), a hyperpolarization-activated inward cation current (H) and fast transient Ca^{2+} (CaT) and slow Ca^{2+} (CaS). The transient calcium conductance, CaT, has a faster time constant than the slow CaS, especially for the closing gating variable. Therefore, this conductance increases rapidly when the cell becomes depolarized, and afterward vanishes rapidly. In contrast, during depolarization CaS increases more slowly than CaT, and afterward persists for longer.

Parameter values for the channel dynamics, including activation and inactivation time constants and voltage dependencies, were chosen to capture rhythmic activity of the lateral pyloric neuron in the stomatogastric ganglion of *Cancer borealis* (Liu et al., 1998) based on experimental work (Golowasch and Marder, 1992):

The capacitance was set to *c* = 1, the membrane potential was *V* and each current (*I*_{i}) was described by the maximum conductance value , the reversal potential *E*_{i}, and the opening and closing variables *m*_{i} and *h*_{i} and the integers *p*_{i} and *q*_{i}. For the leak current, *m* and *h* were set to 1 and *g*_{leak} = 0.01 µS. The reversal potentials *E*_{leak} = –50 mV, *E*_{Na} = 50 mV, *E*_{K} = –80 mV and *E*_{H} = –20 mV were kept fixed, whereas the reversal potential for the calcium-dependent channels was dynamically computed by using the Nernst equation with an extracellular Ca^{2+} concentration of 3 mM following Liu et al. (1998) for the sake of consistency with the other parameters, although the extracellular Ca^{2+} concentration in the crab is 13 mM (Golowasch et al., 1999; Zhang and Golowasch, 2007). The intracellular Ca^{2+} concentration depended on the influx of Ca^{2+} through the Ca^{2+}-dependent channels, as well as on a linear buffer rate:
with a time constant ms. The opening and closing parameters were modeled by the following equations (where we omitted the subscript *i* on *m* and *h* for clarity):

To build single neuron models we selected the maximum conductance values (in µS) from the following ranges: .

### Half-center circuit model

To generate half-center circuits (Fig. 2), two cells were coupled via a graded synaptic connection (Graubard et al., 1980; Johnson et al., 1995) where the synaptic input current *I*_{syn} adds an additional input to Equation 1 for the postsynaptic cell (Cell 1) and depends on the membrane voltage of the presynaptic cell (Cell 2) (Golomb et al., 1994):
with *V*_{half} = 45 mV, *V*_{slope} = –2 mV, ms and *E*_{syn} = –78 mV. Here, denotes the connection from Cell 1 to Cell 2. Similarly, the reciprocal connection from Cell 2 to Cell 1 is . To determine synaptic strengths that generate a specific output, we scanned the space of all possible synaptic connections in a Monte Carlo fashion as follows. For each randomly selected pair of neurons, we chose a pair of synaptic conductances from a grid in the range [0,1] µS with a resolution of 0.025 µS. To introduce variability in an unbiased fashion, we further added jitter to each conductance from the range [–0.0125, 0.0125] µS with a resolution of 0.0025 µS. Each potential circuit was simulated for 4.5 s. We efficiently sampled the entire space of synaptic conductances using a batch approach; performing 32 circuit simulations per batch, we evaluated the output of each circuit. If Cell 1 suppressed Cell 2 for a specific pair of synaptic conductances, we observed that Cell 1 continued suppressing Cell 2 when further increasing the synaptic connection from Cell 1 to Cell 2 (). Thus, we excluded all bigger synaptic conductances from the remaining selection process after each batch. We call this Monte Carlo process “self-refining”.

During this self-refining random sampling, we observed three different outcomes of circuit output: (1) the circuit is considered “functional” when both cells are active, fire regularly and in alternation with a phase difference, period and burst exclusion metric as defined in the main text (Fig. 2*A*,*B*, green circles), (2) both cells are active but they do not fulfill the requirements for the circuit to be considered functional: although they fire in alternation, they do not have the phase difference, period and burst exclusion matrix that we require (Fig. 2*A*,*B*, yellow circles), or (3) one of the cells is silent (Fig. 2*A*,*B*, red crosses). The random sampling was terminated when the circuit output evaluation generated at least three functional circuits (though in Fig. 2*B* we did not terminate the scan for the purpose of visualization), or if the whole grid of synaptic conductances generated no functional circuit output. In our analysis, we used only one of the three functional circuits for each pair of neurons and synaptic conductances.

In summary, our population generation protocol can be described as follows. We started with 750 single neurons. These could be combined into a total of 281,625 pairs, allowing us to form a total of 450,600,000 possible circuits because each pair could be combined into a circuit with 1600 different possible synaptic connections. We scanned all pairs and found a surprisingly large number, approximately half of them (137,028), which could form circuits that satisfy our conditions. In particular, we scanned a total of ∼300 million possible circuits made by these pairs, based on the self-refining procedure. Because of time and computational constraints, we finally randomly selected 7690 circuits for our detailed analysis. The functional circuits were further simulated for a longer period of 120 s to ensure that their regular firing properties persisted.

### Burst detection

During each simulation, an event was counted as a spike when the membrane voltage of the neuron crossed the threshold of −30 mV from below. To extract bursts, we applied a threshold to the interspike intervals (ISIs) that was determined by one-half the sum of 90th percentile of the whole spike train ISI distribution and the minimal ISI. When this threshold was too close to the maximum or minimum values of the ISI (<10 ms), the cell was considered to spike. For coupled cells, all bursts of one cell that occurred between bursts or spikes of the other cell in the circuit were considered one long burst.

In addition, we used the burst exclusion metric *x*_{network} (Grashow et al., 2009) to determine whether the bursts fired in an alternating fashion (*x*_{network} = 1). The metric is computed by calculating the total active time of both cells (*t*_{cell1} and *t*_{cell2}) and the overlap of both cells being active at the same time in the circuit (*O*_{network}). The latter is compared with the overlap time that would happen randomly for uncorrelated cells (*O*_{random}) and the minimum possible overlap time (*O*_{min}) during the total simulation time T_{trial}:

### Stability value

For each circuit we determined the initial unperturbed phase difference ϕ_{0}, defined as the time between the onset of a burst from one cell to the onset of the next burst of the other cell divided by the period of the circuit. This also included bursts consisting of just a single spike. We perturbed the circuit in *N* = 15 steps (in each direction for increases and decreases) and for each step *i* extracted the resulting phase difference ϕ* _{i}*. Then we can calculate the proximity of the phase difference ϕ

*to the initial phase difference ϕ*

_{i}_{0}:

The stability value θ is defined as the average proximity over all perturbation steps (Fig. 4):

### Influence of the unperturbed cell on stability

To determine how specific circuit attributes correlate with circuit stability in response to conductance perturbations in circuits with an identical perturbed cell, we chose a total of 100 subsets, each consisting of 18–29 circuits, for a total of 2002 circuits (Fig. 5). For every subset, we selected one perturbation (e.g., an increase in ) to calculate the Spearman correlation coefficient between the stability value (e.g., ) and an attribute of either the unperturbed cell when uncoupled or an attribute of the unperturbed circuit (e.g., the synaptic connection to the perturbed cell; Fig. 5). We then evaluated the statistical significance of these correlations.

### Statistical analysis

The significance of the Spearman correlation was determined using the scipy function *scipy.stats.spearmanr*. Because the analysis is repeated 100 times, for a significance value of *p* < 0.05 we expect to see significant correlations by chance, even with uncorrelated datasets.

### Hierarchical clustering algorithm

The hierarchical clustering algorithm (Fig. 6) was implemented by the scipy function *scipy.cluster.hierarchy.linkage*. The distance between two columns was calculated as the Euclidean distance (L2 norm). The clustering mechanism uses the linkage method “single”, which assigns a distance d between the clusters *u* and *v* given by the following:
for all points *i* in cluster *u* and all points *j* in cluster *v*.

### Logistic regression classifier

We trained a binary logistic regression classifier for each of the 14 different perturbations of the seven intrinsic conductances in each cell in a circuit (Fig. 7). To prepare the data we assigned the labels “stable” and “unstable” to the 1500 circuits (top and bottom 20%) with the highest and lowest stability value for the chosen perturbation, respectively. As features, we used all 14 maximum conductances (7 for each cell) in the circuit, as well as the strength of the synaptic connection between them. Because of the different ranges of each conductance, each was scaled to the range [0,1]. These data were then used to train a standard logistic regression classifier from the *sklearn.linear_model* Python library. We used 10-fold cross-validation and the L2 norm with a regularization parameter λ = 1.

### Prediction of double perturbation by single perturbation

We fit the function using the scipy function *scipy.optimize.curve_fit* (Fig. 10). For each stability value of the double perturbation (), the two corresponding single stability values ( and θ_{2}) were used to determine the coefficients α, β, and γ.

### Code/software accessibility

All code is provided on github at https://github.com/comp-neural-circuits/circuit-stability.

## Results

Despite the ubiquity of degeneracy in biological systems, studying degenerate circuits that produce similar output from different intrinsic and synaptic conductances is a challenge because of the many timescales of the underlying parameters and the high dimension of the parameter space. First, we describe some of this complexity for a family of degenerate circuit models; two neurons coupled with recurrent inhibition that form a half-center circuit (HCC) with rhythmic output. Then, we propose a computational approach to expose hidden variability in a circuit's intrinsic and synaptic conductances. We achieve this by quantifying circuit stability to single and double conductance perturbations for a large population of circuits using unbiased statistical analysis.

### Constructing a population of degenerate rhythmic circuits with similar output

We modeled individual neurons in the half-center circuit with intrinsic and synaptic conductances based on an established model for the rhythmic activity of neurons in the stomatogastric ganglion of *Cancer borealis* (Golowasch and Marder, 1992; Liu et al., 1998). In particular, we used Hodgkin–Huxley single neuron models with seven different channel types (see Materials and Methods). In contrast to previous mathematical reductions of HCCs (Perkel and Mulloney, 1974; Wang and Rinzel, 1992; Skinner et al., 1993, 1994; Daun et al., 2009), we allowed the intrinsic and synaptic conductances in the HCCs to be different and investigated the robustness of circuit output when perturbing the intrinsic conductances. To build the population of circuits, we first generated single neurons with variable maximum conductances for each channel type, which produced a range of activity patterns, from single spikes to bursts with different periods (Fig. 1*A*). These conductances were randomly and independently chosen from a uniform range of biologically plausible values (Fig. 1*A*; Liu et al., 1998; Goldman et al., 2001; Prinz et al., 2004; O'Leary et al., 2014). The resulting population of 750 neurons had a significant degree of degeneracy, where different sets of maximum conductances produced similar output patterns (Fig. 1*B*).

Following the construction of single neuron models, we combined randomly selected pairs of these neurons into HCCs by coupling them with inhibitory synaptic conductances, where we allowed the synaptic conductances between the two cells to be different (see Materials and Methods). To generate a population of HCCs with similar behavior but distinct intrinsic and synaptic conductances, we required that HCC output satisfied three conditions: (1) non-overlapping activity between the two cells, quantified by a burst exclusion metric close to 1 (see Materials and Methods; Grashow et al., 2009); (2) a phase difference of 0.5 ± 0.03 (measured as a fraction of the circuit period), and (3) a period between 100 ms and 800 ms. Several biological systems maintain such a precise phase difference of 0.5, including the leech heartbeat (Calabrese, 1977) and motor nerve activity generated during swimming in newly hatched *Xenopus* tadpoles (Roberts et al., 2010). A precise phase difference was chosen to give us a population of circuits that despite highly variable intrinsic and synaptic conductances have a similar output. This eventually enabled us to compare the deviation from this chosen phase difference following a perturbation in all simulated circuits. We allowed the period to vary over a range based on experimental results; for instance, the leech heartbeat period can be regulated by changes in temperature and during swimming (Masino and Calabrese, 2002).

A randomly chosen pair from our population of 750 single neurons, when synaptically coupled, could generate one of three possible behaviors: one of the cells was silent, both cells were active but the circuit did not satisfy all three of the above criteria, and both cells were active and the circuit satisfied all three of the above criteria; in which case we called the circuit “functional” (Fig. 2*A*). For each randomly selected pair of neurons, we identified synaptic strengths that led to a functional circuit. For this purpose, we scanned the space of possible synaptic connections in a *self-refining* Monte-Carlo fashion (see Materials and Methods; Fig. 2*B*). This approach enabled us to efficiently sample the entire space of synaptic conductances in a biologically plausible range, without simulating all possible circuits (see Materials and Methods). The process was terminated early when our output evaluation found at least three functional circuits, from which we randomly selected one for our HCC population (Fig. 2*B*, green symbols). For some randomly chosen neuron pairs, the construction of functional circuits required that one synaptic conductance was larger than the other (Fig. 2*B*, left) and for others that the two synaptic conductances had similar strength, yielding a much larger set of functional solutions for the given pair (Fig. 2*B*, middle). Other neuron pairs never generated a functional circuit when synaptically coupled (Fig. 2*B*, right). Overall, of all possible neuron pairs (281,625) we found that ∼50% (137,028) produced functional circuits when scanning through different synaptic conductances (see Materials and Methods). As for the single cells, we found many examples of degeneracy within those functional circuits (Fig. 2*C*). In summary, we have developed an unbiased and efficient approach to generate degenerate HCCs, which generate similar behavior despite different intrinsic and synaptic conductances.

### Intrinsic conductance perturbations affect circuit output differently

We next asked how these circuits with variable intrinsic and synaptic conductances respond to perturbations that fully or partially block or enhance a given ion channel in one of the constituent neurons. Many biological systems are equipped with activity-dependent compensatory mechanisms whereby an activity sensor actively regulates and reconfigures other conductances when one conductance is perturbed to return circuit activity back to a target level (Marder, 2011; O'Leary, 2018). Here, we consider a different form of compensation where following a perturbation of a given conductance, the remaining conductances can together compensate for it, due to the overlap in their timescales of operation of the opening and closing gating variables of the different ion channels (Grashow et al., 2010; Tang et al., 2010, 2012; Soofi et al., 2014). This type of compensation is much faster and as such would be applicable to circuits shortly after a perturbation is applied and when slower activity-dependent mechanisms have not yet had the chance to regulate other channels in response to the perturbation.

We were specifically interested in channel deletions, and therefore incrementally decreased the maximum conductance of a single ion channel to 0. Because there is no upper bound for conductance increases, we incrementally increased the maximum conductance to 200% of its original value. For both increases and decreases, we changed each conductance in 15 discrete steps, which produced 30 perturbed circuits. We simulated these perturbed circuits and scored their functionality (Fig. 3). For some circuits, modifying a specific conductance did not significantly affect circuit output, whereas for others it led to a complete failure to generate a half-center rhythmic pattern or rendered one cell completely silent. We termed this most severe failure a “crash” similar to the experimental nomenclature of activity changes following temperature perturbations (Tang et al., 2010; Haddad and Marder, 2018). For instance, decreasing typically led to silencing the perturbed cell because the hyperpolarization-activated current provides the only source of depolarization in the perturbed cell. This agrees with previous studies, which have shown that the H current plays an important role in generating stable bursting behavior in HCCs (Angstadt et al., 2005; Daun et al., 2009; Grashow et al., 2009), specifically in the context of the 'escape and release' mechanism (Wang and Rinzel, 1992; Skinner et al., 1994; Sharp et al., 1996): When two cells are coupled in a half-center oscillator, one cell can either release the other by reducing its inhibitory influence on the other, or one cell can escape the suppression from the other by hyperpolarization-activated currents, like the H current.

Surprisingly, we found that most of the perturbations silenced the unperturbed, as opposed to the perturbed cell. Decreases in the conductances and , as well as increases in , were the most obvious conductances for which this happened, where >50% of the tested circuits failed to generate a stable output on a perturbation (Fig. 3). This suggests that the effects of each perturbation on circuit output result from an interplay of the different timescales of each conductance in the two constituent cells and are not trivially caused by failures in the perturbed cell. Although the specific way by which circuits crashed following perturbations differed, we observed that when one cell silenced the other from one perturbation step to the next, the number of spikes per burst increased in 77% of the cases compared with the corresponding unperturbed circuits. This increase in spiking was accompanied by a shorter ISI within a burst in 70% of the cases, and occurred when either the perturbed or the unperturbed cell was the silenced cell.

Interestingly, we found that deleting specific ion channels was more detrimental to circuit stability than increasing them, leading to more dysfunctional circuits. This suggests that some conductances may need to exceed a threshold value to ensure circuit stability, i.e., they do not need to be fine-tuned as long as they are sufficiently large. Together, our analysis argues that perturbing any single intrinsic conductance in one of the constituent cells can have drastically different effects on circuit output that cannot be trivially predicted, because of its interactions with other conductances in the same (perturbed) or in the unperturbed cell of the circuit.

### A new measure for quantifying circuit stability in response to intrinsic conductance perturbations

A concise quantification of the response to a specific perturbation is difficult because degenerate circuit models respond to perturbation in different and nonlinear ways (Tang et al., 2010, 2012; Ransdell et al., 2013; Sakurai et al., 2014; Haddad and Marder, 2018; Ratliff et al., 2018; Alonso and Marder, 2019). Classical sensitivity analysis typically estimates how a subtle change in one or multiple system parameters changes the output of the system. This method has been used to demonstrate how much one parameter should change from its original value to compensate a deviation of another parameter from its original value. This typically yields a linear and local approximation of circuit output as a function of parameter change (Olypher and Calabrese, 2007). To quantify changes in circuit output and compare them across circuits and perturbations, we developed a novel quantitative measure of stability based on the properties of circuit output following a perturbation that is global, taking into account the entire range of conductance change.

For each circuit we generated the output in response to individual conductance perturbations ranging from 0 to 200% of the original conductance values (Fig. 4*A*, Step 1). Because the circuits generated very different outputs in response to increased or decreased conductances, we computed the stability values for each direction of change separately; thus, denoted the stability value when increasing the maximum conductance *X* from 100 to 200% and denoted the stability value when decreasing the conductance *X* from 100 to 0%. For each perturbation, we computed the mean phase difference over the range of functional outputs (Fig. 4*A*, Step 2). To quantify circuit stability, we focused on the phase difference only, because this property seems to be preserved on a perturbation in rhythmic circuits. For example, in the case of temperature perturbations in the stomatogastric ganglion, despite period changes, the pyloric circuit remains functional as long as phase relationships are preserved (Tang et al., 2010). Thus, for each perturbation, we calculated the proximity of each resulting phase difference to the initial phase difference of the unperturbed circuit: A value equal to the initial phase difference results in a proximity of 1, and a dysfunctional circuit results in a proximity value of 0. This procedure generated proximity values for each simulated circuit (Fig. 4*A*, Step 3). Finally, we averaged all proximity values for each conductance increase or decrease to get a stability value between 0 and 1. Circuits that maintained a functional output with the same phase difference as the original circuit (0.5 ± 0.03) for the entire range of the perturbation were quantified as stable, with a stability value of 1, whereas circuits that crashed for the smallest perturbation were quantified as unstable, with a stability value of 0. For the majority of circuits, the stability was somewhere in between (2 examples given in Fig. 4*B*).

This method provides us with a concise and quantitative measure of circuit stability in response to perturbations of a single conductance in a given constituent cell of an HCC. We highlight that computing a given stability value does not require knowledge of the specific conductance composition of the neurons, only the phase difference, and thus, has the potential to be used for quantifying stability in other rhythmic circuits with a clearly defined output (see Discussion about the applicability of our approach to other circuits).

### A population of degenerate circuits shows a range of stability to intrinsic conductance perturbations

Equipped with a concise measure of circuit stability, we aimed to quantify the response of the population of degenerate HCCs with similar output but distinct intrinsic and synaptic conductances following a perturbation. Similarly to single neuron models (Alonso and Marder, 2019), HCCs responded to intrinsic conductance perturbations in different ways (Fig. 5). HCC output was particularly vulnerable to some, but not other, conductance perturbations. For instance, the majority of circuits were largely unaffected by increasing , as reflected in the highly skewed distribution of stability values to 1 (Fig. 5*B*–*D*, bars, left). In contrast, decreasing had a strong effect on the majority of circuits (Fig. 5*B*–*D*, bars, right). Yet other perturbations, for instance decreasing showed a range of effects on the circuits as seen in the relatively uniform distribution of stability values (Fig. 5*B*–*D*, bars, middle).

Some of the studied circuits shared one of the neurons, i.e., they had either an identical perturbed or unperturbed neuron. One might conceive that when the circuits shared the perturbed neuron, the identity of the second, unperturbed neuron did not matter for the stability of circuit output. To determine this, we compared a subset of 20 circuits where the perturbed cell was identical (Fig. 5*B*–*D*, colored dots). The synaptic conductances differed across the different circuits to ensure that the different unperturbed cell could still generate functional HCC output based on our criteria when coupled with the perturbed cell. Particularly, we wondered whether the stability values in response to a given perturbation correlated with specific attributes of the circuit. We searched for correlations between or and different circuit attributes, including: the synaptic conductance from the unperturbed to the perturbed cell (Fig. 5*B*, blue dots), the synaptic conductance from the perturbed to the unperturbed cell (Fig. 5*C*, purple dots), and the period of the unperturbed cell when uncoupled from the circuit (Fig. 5*D*, green dots). The subset of circuits that had the same perturbed cell showed diverse correlation relationships between the stability values and the circuit attributes, with a range of statistical significance (Fig. 5*B*–*D*. We found that some significant correlations have interesting implications; for instance, circuits with high in the perturbed cell seemed to receive a strong synaptic connection from the unperturbed cell (Fig. 5*B*, left), but provide a weaker synaptic connection to the unperturbed cell (Fig. 5*C*, left). These same circuits were also characterized with longer periods of the unperturbed cell (Fig. 5*D*, left). Furthermore, circuits with low , seemed to receive a stronger synaptic connection from the unperturbed cell (Fig. 5*B*, right), but had a variable synaptic connection to the perturbed cell (Fig. 5*C*, right).

To summarize these findings in an unbiased manner across many simulated circuits, we found a total of 100 such distinct subsets each consisting of ∼20 (18–29) circuits that have the same perturbed cell but different unperturbed cells, giving us a total of 2002 circuits. For the circuits within these subsets, we computed the correlations between each of and with each of and the period of the unperturbed cell when uncoupled. Considering the statistically significant correlations (*p* < 0.05) revealed some general relationships (Fig. 5*E*), was generally positively correlated with and but negatively correlated with (Fig. 5*E*). This suggests that having a high , the synaptic conductance from the unperturbed to the perturbed cell, makes circuits ultra-sensitive to decreases in , consistent with a classical “escape” mechanism whereby one cell escapes inhibition from the other during rhythm generation (Wang and Rinzel, 1992; Skinner et al., 1994; Sharp et al., 1996). The special role of in rhythm generation is further highlighted in the lack of correlation between and the period of the unperturbed cell when uncoupled, in contrast to the other positive correlations with and (Fig. 5*E*). We also found that was negatively correlated with all and in most circuits (Fig. 5*E*). This implies that a stronger synaptic conductance from the perturbed to the unperturbed cell will result in a lower circuit stability. Interestingly, a similar relationship was found experimentally in the mollusk, *Tritonia diomedea*, where the extent of motor impairment during swimming following injury correlates with the inhibitory synapse strength converging onto a specific neuron (Sakurai et al., 2014).

Together, these results demonstrate that the relationships between circuit stability when perturbing a single conductance and distinct circuit attributes can be quite complex. For some circuits, having a different unperturbed cell has a huge impact on circuit stability after a perturbation and the resulting stability exhibits strong correlations with specific circuit attributes. For other circuits, the different unperturbed cell continues to have a big impact on circuit stability, but the stability fails to correlate with specific circuit attributes. For yet other circuits, having a different unperturbed cell has minimal influence on circuit stability. These results highlight that circuit stability depends on the entire circuit, rather than the individual cell that is perturbed, and demonstrate the need for an unbiased analysis of the effects of perturbations on circuit stability, which is what we sought to accomplish next.

### Perturbations reveal subsets of conductances that co-regulate circuit stability

Applying a perturbation to a specific conductance can be interpreted as moving through the parameter space of all intrinsic conductances along one direction (the perturbed parameter) to find parameter combinations that yield a similar circuit output. Our newly defined stability measure allows us to quantify this movement at different points in this parameter space for a large population of circuits. More specifically, the measure describes whether, following a perturbation, the circuit retains its output (high stability) or not (low stability). In this view, we postulated that any correlations between two stability values should locally describe the shape of the parameter space in which functional circuits reside. Therefore, we examined the correlation between all possible pairs of stability values in the population of degenerate circuits (Fig. 6*A*). Surprisingly, the majority of the stability values were uncorrelated (e.g., vs and vs ; Fig. 6*B*, top). Additionally, no strongly negative correlations were observed, which suggests that high circuit stability when perturbing one conductance does not imply low circuit stability when perturbing another conductance. Several pairs of stability values exhibited strong correlations (e.g., vs and vs ; Fig. 6*B*, bottom). Applying a hierarchical clustering algorithm (see Materials and Methods) to reorganize the columns of the correlation matrix revealed four clusters of strong positive correlations (Fig. 6*A*, yellow boxes), uncovering subsets of conductances that when perturbed produced a similar effect on diverse individuals of the population, thus co-regulating circuit stability. Two of the clusters ( vs and vs ) imply co-regulation of stability by opposite changes in and . This relationship predicts that a neuromodulator that increases (decreases) the A current has similar effects on circuit stability as a neuromodulator that decreases (increases) the CaS current. A plausible explanation for this correlation in stability may be that the activation functions of the corresponding currents operate on a similar timescale, but in opposite directions; the A current has a hyperpolarizing action, whereas CaS current has a depolarizing action (Table 1). Similarly, the Kd current operates in the same direction and with a similar activation timescale as A, consistent with being a member of the stability cluster containing and .

Particularly interesting is the large cluster consisting of the stability values to five different perturbations, which contains and . The H and KCa conductances are critical for the robustness of circuit output for several reasons. First, even weakly increasing or decreasing affects the membrane potential of the perturbed cell right after its burst ends. During normal rhythm generation, the H current gradually vanishes over the course of the burst, whereas the KCa current builds up. During a perturbation, however, slightly decreasing leads to a smaller hyperpolarizing current after the burst, while slightly increasing leads to a larger depolarizing current in the same time window. Applying a stronger perturbation that more significantly decreases , or increases , effectively prevents the intrinsic mechanism of the perturbed cell to terminate its bursting, leading to a severe crash where the unperturbed cell is silenced (Fig. 3). Although the other stability values in the large cluster correspond to perturbations of conductances that act on different timescales (Table 1), their ultimate effect on circuit output is likely similar; increasing or and decreasing lead to a longer ISI and thus shorten the time window where the membrane voltage is depolarized. For instance, a bigger Na^{+} current speeds up the depolarization during each spike in a burst, preventing other depolarizing currents with much slower timescales to build up. Immediately following a spike, the Na^{+} current no longer contributes to the membrane potential because the Na^{+} channels are completely closed. At this time, the lack of other depolarizing currents becomes apparent, leading to a stronger after-spike hyperpolarization and eventually a longer ISI within a burst. Because of the longer ISI and more hyperpolarized membrane potential, the KCa current fails to build up (or the H current is not reduced as strongly) as before the perturbations, resulting in similar effects as decreasing (or increasing ) directly.

We also identified isolated clusters of stability that did not correlate with any other, and (Fig. 6*A*). Perturbations in these conductances had unique effects on circuit output stability independent of the other conductances. Increasing did not significantly impact circuit output. In contrast, decreasing had severe impact on circuit output, because a threshold amount of is necessary for the circuit to generate the desired bursting output (Fig. 3; Wang and Rinzel, 1992; Skinner et al., 1994; Sharp et al., 1996; Angstadt et al., 2005; Daun et al., 2009; Grashow et al., 2009).

Together, our analysis uncovers subsets of intrinsic conductances that co-regulate circuit stability by examining only the response to perturbations of these same conductances in a population of degenerate circuits with entirely uncorrelated intrinsic and synaptic conductances. By identifying subsets of perturbations that act alike, the correlations between pairs of stability values at different locations within the parameter space of intrinsic conductances provide insights into the shape of the parameter space where functional circuits live.

### Features predictive of circuit stability to single conductance perturbations

We next wondered whether certain intrinsic or synaptic conductance combinations can predict circuit stability to specific conductance perturbations. To extract predictive circuit features in an unbiased fashion in the entire population of degenerate circuits, we used a logistic regression classifier with the intrinsic conductances of each cell and the two synaptic conductances as features (total of 16). We trained a different classifier for each perturbation (see Materials and Methods; Fig. 7*A*). We refer to the classes as stable and unstable, with each containing the 20% most stable or 20% most unstable circuits, respectively. While the classification rate was >75% for most of the classifiers (Table 2), each classifier weighted the features differently to predict circuit stability. This information enabled us to look at the features with the largest weights as those that most strongly predict whether a circuit is stable given a specific perturbation. For example, the feature “ in the perturbed cell” for the classifiers of (Figs. 6*A*, green box, 7*B*) and (Fig. 6*A*, orange box, 7*B*) shows that a circuit with a high value of in the perturbed cell will result in low stability values and because of the negative weights for each classifier (Fig. 7*B*, green dotted box), but in a high stability value in due to the positive weight (Fig. 7*B*, orange dotted box). Therefore, a high value of in the perturbed cell makes the circuit more likely to be sensitive to decreases in and and increases in , but less sensitive to decreases in . In the parameter space of all intrinsic conductances, this analysis reveals that the initial position in parameter space, defining the unperturbed circuit, determines the stability of that circuit with respect to the different perturbations, which move the circuit along different directions.

To compare general differences between the weights of the classifiers without looking at each weight individually, we computed the angle between the weight vectors for each classifier. To visualize the angles between the weight vectors for the four classifiers highlighted in Figure 7*B*, we projected the vectors into the space spanned by the first three principal components of the weights of all 14 classifiers (Fig. 7*C*). The angles capture the similarity between the projected weight vectors. Here, the orange vector is nearly orthogonal to the green vectors. Calculating the angles between the weight vectors in the full 16-dimensional space for all classifiers produces a matrix similar to the correlation matrix between stability values (compare Figs. 6*A*, 7*D*). This suggests that the specific subsets of conductances that co-regulate circuit stability, identified with unbiased clustering, correspond to circuit features that are predictive of circuit stability to perturbations.

These two results can also be reconciled in the parameter space of all intrinsic conductance in which a perturbation can be interpreted as a movement through parameter space. The previously identified correlations in the stability to a perturbation are independent from the initial position in the parameter space (which defines the unperturbed circuit), because they were determined for the entire population of degenerate circuits. Complementary to that, the classifier developed here identifies the direction along the weight vector within the parameter space, which leads to the best classification. Thus, when the two stability values in response to perturbations are perfectly correlated, the weight vectors completely overlap. Indeed, we find that the perturbations with highly correlated stability values (Fig. 6*A*) also have small angles between the weight vectors of their classifiers (Fig. 7*B*).

### Quantifying circuit stability in response to double conductance perturbations

Considering that neuromodulators often impact more than one intrinsic conductance simultaneously (Ratté et al., 2014; Ratté and Prescott, 2016; Ratliff et al., 2018), it is important to understand how to determine circuit stability in response to multiple conductance perturbations. Few studies have explored the effects of the interaction of multiple conductance perturbations at the circuit level, especially providing quantitative descriptions of co-modulation (Harris-Warrick, 2011; Marder, 2012; Nadim and Bucher, 2014; Li et al., 2018). This co-modulation can directly result from a perturbation that simultaneously affects two conductances or from a secondary effect where a perturbation in one conductance triggers a change in a second conductance. To address this, we next extended our measure of stability to simultaneously perturbing two intrinsic conductances in one cell of the HCCs (Fig. 8*A*). These double perturbations resulted in an order of magnitude more (600 vs 30) perturbed circuits. As before, we calculated the phase difference between the two cells for each perturbed circuit (for an example circuit, see Fig. 8*B*). This generated a set of directions in the two-dimensional plane of conductance changes. The horizontal and vertical directions correspond to the single perturbations considered previously (compare Figs. 8*A*, 4*A*). Other directions than the cardinal correspond to perturbations in which both conductances were simultaneously perturbed at a given ratio. For each of these directions, we extracted a single stability value using the same procedure as for the single perturbations (compare to Fig. 8*B*, 4*A*). These stability values were then examined in a two-dimensional space (Fig. 8*C*).

To quantitatively describe the stability of all circuits, we generated a two-dimensional histogram where each quadrant provides important information about how the two conductance perturbations interact (Fig. 9*A*–*C*, left). For instance, the top quadrant corresponds to perturbations where the H conductance was increased to 200% from its reference value (Fig. 9*A*), whereas at the same time the A conductance was either decreased (left half) or increased (right half). One naive expectation is that these double perturbations affect circuit output more than perturbations in the single conductances, thus leading to lower stability values along perturbation directions different from the cardinal. However, we found this not to be the case. In fact, simultaneous perturbations of two conductances often improved circuit stability relative to the single perturbations. For instance, when decreased, decreasing rather than increasing improved the stability for most circuits (Fig. 9*A*, left quadrant, a peak in the histogram near 1 in the bottom half). In contrast, when increased, increasing rather than decreasing improved stability (Fig. 9*A*, right quadrant, a peak in the histogram near 1 in the top half). Therefore, we find that changing in the same direction as increases circuit robustness. This finding is consistent with previous experimental results that and are co-regulated (MacLean et al., 2003; Amendola et al., 2012) and actively kept at a fixed ratio to maintain a given output (Krenz et al., 2013, 2014).

Extracting a similar relationship when examining the effect of changes in after perturbing appears to be more challenging (Fig. 9*A*, top and bottom quadrants). In particular, we find a peak in the histogram near 0.2 when decreasing and increasing , which agrees with our previous finding that stability is decreased when and change in opposite directions (Fig. 9*A*, bottom quadrants). To better illustrate the impact of perturbations in on the circuits already perturbed in , we plotted circuit stability for the double perturbations, and , as a function of circuit stability for the single perturbation, (Fig. 9*A*, right). As for the single perturbations, perturbing had a strong effect on circuit stability because decreasing removes a cell's ability to burst (Fig. 3; Grashow et al., 2010). This is a minimal requirement for which changes in cannot compensate. Furthermore, circuit stability when changing in addition to decreasing is stereotypic: Nearly all circuits became more stable when also decreased (Fig. 9*A*, right, red) and more unstable when increased (Fig. 9*A*, right, green). This relationship was surprisingly robust across the vast majority of tested circuits, although the intrinsic conductances in the population of degenerate circuits were uncorrelated by construction, thus confirming that and need to be co-regulated for circuit stability.

Previously, the single conductance perturbations revealed that changes in and , when applied in opposite directions, have similar effects on circuit stability (Fig. 6). This relationship persisted when considering the double perturbation (Fig. 9*A*): Simultaneously increasing and decreasing decreased stability in the majority of circuits compared with only increasing or only decreasing (Fig. 9*B*, right). In contrast, changing and in the same direction enhances circuit stability (Fig. 9*B*). This opposing action of A and CaS is further evident when comparing the double perturbations involving and with those involving and . Circuit stability when comparing changes in in addition to decreasing ( and ) is mirror opposite to circuit stability with changes in in addition to decreasing ( and ; Fig. 9, compare *A*, *C*, right). These findings confirm that even more than for single perturbations, the circuit response to double perturbations is non-intuitive, highlighting the strength of our newly derived stability measure to summarize this diversity in responses.

### Explaining circuit stability to double perturbations by combining single perturbations

Next, we asked whether circuit stability in response to double perturbations can be predicted from the stability to single perturbations, and , by fitting the nonlinear function with constants α, β and γ: (Fig. 10*A*). The fitted exponent γ tells us whether the double perturbation results in a lower (γ > 1) or higher (γ < 1) stability than the weighted combination of stability to single perturbations. In addition, we also considered the goodness of fit, *R*^{2}. We found that a fit withγ > 1 (supralinear) explains the double perturbation better (has a higher *R*^{2}) than a fit with γ < 1 (sublinear; Fig. 10*A*, middle). This can tell us more about the impact of the interaction between the two simultaneously perturbed conductances on circuit stability.

We first investigated fits with high γ and *R*^{2}, which indicate that simultaneously perturbing both conductances more strongly destabilizes circuit output than when perturbing the conductances one at a time. For instance, the stability to the double perturbation can be well explained by the supralinear combination of and with (Fig. 10*B*, left). We previously showed that and are positively correlated (Fig. 6*A*). Interestingly, we found this to be a general result: the well correlated stability values to single perturbations (Fig. 6*A*) could be supralinearly combined to yield excellent predictions for the stability to double perturbations (Fig. 10*A*, right). The high correlation between the two stability values to the single perturbations indicates that these conductances have the same effect on circuit output when perturbed. Thus, simultaneously perturbing those same conductances yields even more perturbed output, suggesting that their interaction is weak. When the two stability values in response to single conductance perturbations are weakly correlated (e.g., and ), the double perturbation leads to better circuit stability than the nonlinear combination of the stability to the single perturbations (Figs. 10*B*, middle, 9*B*). Therefore, the interaction between the two perturbed conductances is likely strong, and the double perturbation cannot be well-explained by the nonlinear combination of single perturbations (low *R*^{2}).

An exception to the relationship between γ and *R*^{2} are stability values involving perturbations that decrease H (Fig. 10*A*, black symbols), due to the powerful effect of this conductance on bursting in HCCs as previously discussed. Fitting the stability to the double perturbation involving a decrease in typically results in high *R*^{2} values because the single perturbation involving tends to dominate the double perturbation. This is reflected in the large weight assigned to compared with the stability when perturbing the second conductance (e.g., in Fig. 10*B*, right).

When considered together, our results reveal a set of conditions under which a double perturbation can be explained by a nonlinear combination of single perturbations, and when such a double perturbation improves or worsens circuit stability, depending on the interaction strength of individual conductances in a circuit. Beyond the specifics of the family of half-center circuits we consider, this finding cautions against studying single conductance perturbations one at a time to determine the effect of neuromodulators or external perturbations that simultaneously disrupt more than conductance.

## Discussion

We developed quantitative measures to reveal hidden variability in a circuit's intrinsic and synaptic conductances from the responses of a population of model circuits to distinct perturbations (increases and decreases) in the intrinsic conductances. Based on the simulated output to such perturbations and unbiased statistical analysis, we showed how these models can be used to develop intuitions about which conductance combinations predict circuit stability to particular perturbations when all other conductances are unknown. We focused on the smallest circuit of two reciprocally coupled neurons with inhibition, a common motif found in rhythmic networks (Brown, 1911; Perkel and Mulloney, 1974; Calabrese and De Schutter, 1992; Marder and Calabrese, 1996; Daun et al., 2009; Harris-Warrick, 2011; Doloc-Mihu and Calabrese, 2014; Dethier et al., 2015). Our circuits were composed of non-identical conductance-based neurons with multiple conductances preserving essential dynamics for bursting (Dethier et al., 2015; Franci et al., 2013, 2018). This enabled us to achieve a balance between computational tractability and biological realism.

Understanding the cooperativity rules of intrinsic conductances is a challenging task due to the surprising amount of degeneracy in achieving a particular circuit output. The parameters that govern intrinsic excitability, synaptic strength, and neuronal architecture can vary several-fold, despite nearly identical circuit output (Goldman et al., 2001; Prinz et al., 2003, 2004; Swensen and Bean, 2005; Marder and Goaillard, 2006; Schulz et al., 2006, 2007; Tobin et al., 2009; Doloc-Mihu and Calabrese, 2011; Roffman et al., 2012; Temporal et al., 2012; Srikanth and Narayanan, 2015). One reason for having so many degrees of freedom to tune circuit function may be the circuits' ability to compensate perturbations. There are at least two complementary ways to achieve such compensation. One way relies on an activity sensor, which actively regulates the unperturbed conductances to return the circuit to a target activity level. The mechanisms implementing this form of compensation, including synaptic and homeostatic plasticity, function over long timescales after perturbation (Marder, 2011; O'Leary, 2018). The other way to achieve compensation, pursued in our study, does not rely on active mechanisms but on the overlap of different timescales in the interacting neuronal conductances (Grashow et al., 2010; Tang et al., 2010, 2012; Soofi et al., 2014), and thus depends on the context provided by all other conductance identities and expression levels. This compensation occurs on a much faster timescale and before the onset of activity-dependent mechanisms. A common approach to understand conductance interactions has been to build on identified correlation relationships between individual parameters (Schulz et al., 2006, 2007; Khorkova and Golowasch, 2007; Goaillard et al., 2009; Calabrese et al., 2011), arguing that reliable circuit output and robustness to perturbation is encoded in the correlation rules, rather than the value of any one parameter (Olypher and Calabrese, 2007; Tobin et al., 2009; Zhao and Golowasch, 2012). Our approach addressed the problem in reverse: starting from a measured output to a given perturbation in a specific two-cell circuit, we revealed conductance subsets of the circuit neurons that similarly affect circuit stability when perturbed, as well as conductances that predict circuit stability after a perturbation.

Single neurons in biological circuits have many conductances, producing numerous parameter combinations that maintain appropriately tuned electrical phenotypes; this problem intensifies at the circuit level, making it difficult to evaluate all possible changes to these circuits. To surmount this curse of dimensionality, we developed a novel measure to quantify circuit stability when modifying intrinsic conductances in biologically plausible ranges based on how significantly circuit output was altered (Fig. 4). This measure provided a concise description across the entire population of model half-center circuits to the same perturbation, as opposed to qualitatively describing a few representative examples (Grashow et al., 2009, 2010; Lamb and Calabrese, 2013; Dethier et al., 2015). Compared with classical sensitivity analysis (Olypher and Calabrese, 2007), our measure does not assume a linear change of circuit output for small and local changes in the modulated parameter but embodies strong nonlinearities especially when the circuits fail to generate functional output. This allowed us to scan the entire space of model instances without being restricted to local sensitivity effects (Weaver and Wearne, 2008) and consequently probe interactions between parameters. Previous studies have also attempted to define more extensive robustness measures (Doloc-Mihu and Calabrese, 2014, 2016) but in a more qualitative sense.

The quantitative nature of our stability measure enabled us to evaluate stability correlations to different perturbations in the family of HCC we studied. Unbiased clustering of the stability to single conductance perturbations led to the discovery of specific subsets of intrinsic neuronal conductances that, when modified, co-regulate the circuits' response to perturbations in these conductances (Fig. 6). Interestingly, the identified subsets of conductances resulted from generalizing over all the circuits in our population of circuits with entirely uncorrelated sets of conductances. For instance, altering the A and CaS currents in opposite directions resulted in similar circuit stability after a perturbation, consistent with the similar activation timescales of these currents. Decreases in the H current had a severe impact on circuit function, in agreement with the important role that this current plays in the generation of stable bursting (Wang and Rinzel, 1992; Skinner et al., 1994; Grashow et al., 2009). The strong correlations of stability in response to single perturbations were directly related to subsets of conductances, which were used as features to predict circuit stability using nonlinear classifiers (Fig. 7), significantly extending linear classification techniques like principal component analysis (Grashow et al., 2010; Doloc-Mihu and Calabrese, 2014).

Many neuromodulators either target the same ion channel or can have distinct targets within the same neuron or circuit (Harris-Warrick, 2011; Marder, 2012). A wealth of experimental data has shown that the co-modulatory actions of converging neuromodulators can be similar or opposing, resulting in additive, synergistic, antagonistic, or other nonlinear co-modulatory effects (Nadim and Bucher, 2014). Therefore, we investigated when the stability to single perturbations can be combined to predict circuit stability to perturbations in two conductances simultaneously. Strong correlations between stability to the single perturbations resulted in much worse circuit stability when those same conductances were simultaneously perturbed, as in the case of increasing CaS and decreasing A (Fig. 9). For conductances that exhibited weak correlation when perturbed independently, the double perturbation often led to higher stability compared with the single perturbations applied alone, as in the case of increasing CaS and increasing A (Fig. 9). It is important to note that we only considered perturbations to the maximal conductances in individual neurons in the circuits. Many neuromodulators, however, also act on the level of synaptic transmission, and may follow different rules at different subcellular targets (Nadim and Bucher, 2014; Li et al., 2018). Moreover, we did not explicitly model neuromodulatory action, only the effect that it could have on intrinsic conductance changes. Still, our approach is an important advance toward understanding how distinct conductances might interact to shape circuit output by quantitatively evaluating the interaction rules when modifying several intrinsic neuronal conductances simultaneously. To our knowledge, no previous analysis of stability has been applied in such a nonlocal and nonlinear manner.

Some of the conductance relationships we identified are already known. This includes simultaneously modifying A and H in the same direction (Fig. 9*A*), which improved circuit stability consistent with experimental data, which shows that these two currents are co-regulated and actively maintained in the same ratio (MacLean et al., 2003; Amendola et al., 2012; Krenz et al., 2013, 2014). Importantly, our analysis also revealed novel relationships. For instance, modifying A and CaS in opposite directions one at a time had similar effects on circuit stability, possibly due to the similar timescales of operation. We also found that increasing the H current and decreasing the KCa current affected the perturbed cell's ability to release the unperturbed cell from inhibition. We identified similar effects on circuit stability from increasing the Na or Kd currents and decreasing the CaT current. Perturbing any one of a larger subset of conductances (H↓, A↑, CaS↓, Na↓) silenced the perturbed cell (Fig. 3) preventing it from escaping the suppression of the unperturbed cell. To fully understand how the “escape and release” mechanisms are combined by the different circuits will require a concise classification into each of these categories before and after a perturbation (Skinner et al., 1994).

Although our approach was developed for a specific family of circuits, it can be applied to other circuits for which a well-defined output exists (in our case, the phase difference between the two neurons). One would need to know the individual conductances in the circuit neurons and the connectivity diagram to generate a population of circuits, and perform a similar large-scale perturbation analysis. However, once this analysis has been performed, it can be used to identify combinations of conductances that co-regulate and are predictive of stability for those circuits solely from the circuit's response to perturbation. The outcome of the analysis will depend on the overlap of timescales of the opening and closing gating variables of each conductance, which determine whether specific conductances can compensate for others on a perturbation. Going beyond the details of the circuit model we considered, our quantitative and unbiased approach to characterizing high-dimensional degenerate systems, thus, provides a novel framework to understand circuit robustness and to guide future experimental and theoretical studies.

## Synthesis

Demb, Jonathan Yale University

## Footnotes

This work was supported by the Max Planck Society. We thank Eve Marder and Marion Silies for valuable comments and critical feedback on the paper; members of the Gjorgjieva group for discussions; and Christoph Miehl, Leonidas Richter, Luis Riquelme, and Marina Wosniack for reading drafts of the paper.

The authors declare no competing financial interests.

- Correspondence should be addressed to Julijana Gjorgjieva at gjorgjieva{at}brain.mpg.de

This is an open-access article distributed under the terms of the Creative Commons Attribution License Creative Commons Attribution 4.0 International, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.