Circuit function is more robust to changes in synaptic than intrinsic conductances

Circuit function results from both intrinsic conductances of network neurons and the synaptic conductances that connect them. In models of neural circuits, different combinations of maximal conductances can give rise to similar activity. We compared the robustness of a neural circuit to changes in their intrinsic versus synaptic conductances. To address this, we performed a sensitivity analysis on a population of conductance-based models of the pyloric network from the crustacean stomatogastric ganglion (STG). The model network consists of three neurons with nine currents: a sodium current (Na), three potassium currents (Kd, KCa, A-type), two calcium currents (CaS and CaT), a hyperpolarization-activated current (H), a non-voltage-gated leak current (leak), and a neuromodulatory current (MI). The model cells are connected by seven synapses of two types, glutamatergic and cholinergic. We produced one hundred models of the pyloric network that displayed similar activities with values of maximal conductances distributed over wide ranges. We evaluated the robustness of each model to changes in their maximal conductances. We found that individual models have different sensitivities to changes in their maximal conductances, both in their intrinsic and synaptic conductances. As expected the models become less robust as the extent of the changes increase. Despite quantitative differences in their robustness, we found that in all cases, the model networks are more sensitive to the perturbation of their intrinsic conductances than their synaptic conductances.


I. Introduction
Since the early studies of Hopfield and Hebb (Hebb, 1949;Hopfield, 1982;Sejnowski and Tesauro, 1989) it has been often assumed that changes in synaptic strength provide the salient mechanism by which circuit dynamics are altered in learning (Siegelbaum and Kandel, 1991;Dayan and Abbott, 2001;Tsodyks et al., 1998;Abbott and Nelson, 2000;Sjöström et al., 2001).This tradition was reinforced by years of extracellular recordings of spikes in brain circuits, and by a myriad of computational studies using simplified model neurons such as rate models or integrate and fire neurons (Kempter et al., 1999;Fusi, 2002;Caporale and Dan, 2008;Fiebig and Lansner, 2017).While all of these studies have been extremely informative, they have led to an underestimate of the role of intrinsic membrane currents in network dynamics in larger networks.
At the same time, it is now patently clear that most, if not all, brain circuits display many voltage and time-dependent currents that give their neurons interesting and complex cellular dynamics (McCormick and Huguenard, 1992;Markram, 1997;Ranjan et al., 2020;Grienberger and Magee, 2022).Despite this, there are still prevailing beliefs that changes in synaptic strength that occur as a consequence of activity, represent the dominant mechanism for altering circuit dynamics.In contrast, those who have been studying the role of synaptic and intrinsic currents in circuit dynamics in small central pattern generating circuits have been acutely aware of the importance of modulating intrinsic currents (Calabrese and De Schutter, 1992;Swensen and Marder, 2000;Marder and Bucher, 2001;Harris-Warrick, 2011;Yang et al., 2022) and in some cases the relative insignificance of modifications of synaptic strength (Prinz et al., 2003b).Nonetheless, to the best of our knowledge, there have not been systematic sensitivity analyses in circuits with easily measurable outputs to compare the relative importance of modification of synaptic and intrinsic conductances for network performance.
In this study, we take advantage of the stereotyped motor patterns of the crustacean pyloric rhythm (Maynard, 1972;Selverston, 1976) to assess its stability to changes in its conductances.Because we did not want these results to be an idiosyncratic consequence of the set of starting parameters, we generated 100 degenerate models (different sets of parameters with similar outputs) to explore the generality of the results in this basic circuit, and then evaluated the robustness of each of these models to changes in circuit parameters.

A. Feature extraction to classify model network activity
The pyloric rhythm in crustaceans is produced by the sequential bursting of the PD (pyloric dilator), LP (lateral pyloric), and PY (pyloric) neurons.The AB (anterior burster) neuron is electrically coupled to the PD neurons and together they form a pacemaker kernel.
In these simulations the AB and PD neurons are represented by a single compartment cell, which we call PD.We modified the model of the pyloric network in Prinz et al. (2004) to include a neuromodulatory current, I M I , in each cell (Golowasch and Marder, 1992b;Swensen and Marder, 2000).Figure 1A shows a schematic of the network studied in this work.Each cell is modeled as a single compartment with multiple currents (Golowasch and Marder, 1992a;Buchholtz et al., 1992;Goldman et al., 2001).The interactions in the network are mediated by seven chemical synapses of two types: glutamatergic and cholinergic (Prinz et al., 2004).Each neuron has a sodium current, I N a ; transient and slow calcium currents, I CaT and I CaS ; a transient potassium current, I A ; a calcium-dependent potassium current, I KCa ; a delayed rectifier potassium current, I Kd ; a hyperpolarizationactivated inward current, I H ; a leak current, I leak ; and a neuromodulatory current I M I .
The voltage traces in Figure 1B are an example of normal triphasic pyloric activity produced by one of the models in this study.To evaluate the activity of the network, we performed a number of measurements on the voltage traces of each cell.The colored lines correspond to the quantities used to calculate features of the network's activity.Duty cycles are calculated as the burst duration divided by the bursting period, τ b , of the cell.The phases of the rhythm at which the cells begin and end their bursts, are measured using the start of the PD burst as a reference (Figure 1B).Because the network activity is approximately periodic, it is possible to normalize the times of these events by the network period τ n .In Fournier Z, Alonso LM, and Eve Marder this way, we computed a number between 0 and 1 that corresponds to the ON and OFF phases for each cell.By computing these quantities for several periods of the activity, we obtained a distribution of values of duty cycles, burst periods, and ON/OFF phases, for each cell.The mean values and standard deviations of these distributions were used as features for classification purposes.For example, a cell that produces irregular bursts will yield a distribution of duty cycles with a standard deviation that would be larger than in the case of a cell that bursts regularly in a normal triphasic rhythm.
The activity of the network depends on the values of the maximal conductances ḡi .In this study we compare the sensitivity of the network to changes in intrinsic versus synaptic conductances.Figure 1C shows a normal (control) triphasic rhythm in black.The blue trace in Figure 1C corresponds to the activity of the same model when each of the intrinsic conductances are changed by a random amount.In this case, altering the intrinsic conductances results in a pattern of activity that is different from the control activity and is not "pyloric" because the cells fail to produce bursts of action potentials.Figure 1D shows control activity in black, and activity when the synaptic conductances are changed in red.
In this case the pattern is similar to control activity, but there are some differences.Despite only changing the synaptic conductances, the network's activity becomes slower and the phases of the cycle at which cell fires are slightly different.Figure 1E shows the control values of the intrinsic conductances in black, and the changed values of the conductances in blue used in the simulations in Figure 1C.For simplicity, we only show how the g N a and g Kd were changed, but all intrinsic conductances were changed by a similar percent.In the LP cell, the sodium g N a conductance decreased while the delayed-rectifier potassium ḡKd conductance increased, leading to a loss of spiking activity.The sodium g N a maximal conductance in the PY cell also decreased but g Kd decreased even further, thus maintaining spiking activity.These effects are due to the relative magnitude of g N a and g Kd in the cell, as these are the two largest maximal conductances.Finally Figure 1F shows the control values of the synaptic maximal conductances in black and the perturbed values in red.

B. A model database of degenerate solutions
Conductance based models of neuronal activity are often employed to test and generate hypotheses about the functioning of individual neurons and/or circuits.One difficulty in working with these models is that tuning them to perform specific tasks can be challenging.
Finding values of the parameters, such as the maximal conductances ḡi , so that a given circuit will display a specific set of behaviors has been done in the past by painstakingly hand-tuning the models' parameters (Traub et al., 1991;Nadim et al., 1995;Buchholtz et al., 1992).One alternative to hand-tuning consists of evaluating all possible combinations of parameter values within a specified domain, building a model database, and filtering the solutions of the models that satisfy a set of criteria (Prinz et al., 2003a(Prinz et al., , 2004;;Crasto et al., 2007).Another alternative was introduced recently in Alonso and Marder (2019).
For any given set of parameters, it is possible to simulate the model, and score the solutions depending on their features.The approach consists of defining a cost function that takes values of conductances as inputs, and returns a score.The goal then is to define a cost function such that low scores correspond to solutions of the model that display the target or desired behavior.Thus, training the model corresponds to optimizing -finding local minima -of the cost function.This approach was employed successfully to produce models of single neurons and pyloric networks that can operate robustly over a range of temperatures (Alonso and Marder, 2020).In this study we followed the approach in Alonso and Marder (2020).
We modified the model in Prinz et al. (2004) to include a modulatory current I M I and we also modified the cost function in Alonso and Marder (2020) to enforce a greater voltage separation between slow oscillations and spikes (see methods).In this work, optimization of the cost function yields values of conductances such that the circuit displays a pyloric rhythm as in Figure 1B.We used a genetic algorithm to generate a database of N = 100 models with different values of maximal conductances (Goldberg and Holland, 1988).All models in the database display pyloric rhythms whose features -such as frequency and burst duration -are consistent with experimental measures performed in the pyloric circuit in the crustacean STG (Hamood et al., 2015).
Fournier Z, Alonso LM, and Eve Marder The values of the conductances in our model database are widely distributed over manyfold ranges.For example, Figure 2A and 2B shows the distributions of the PD A-current conductances and the PD to PY glutamatergic synapse conductances.These distributions span a broad range and are not normally distributed (p < 0.05 Shapiro-Wilk test (Shapiro and Wilk, 1965)).In this study, each model is fully specified by its set of maximal conductances ḡi , while all other parameters such as reversal potentials, activation functions, membrane capacitance, were kept fixed and are the same across models.The total number of conductances in a model is N g = 34.To try to visualize the conductances in the database we employed a popular technique (t-SNE) to visualize high-dimensional vectors.Figure 2C shows a t-SNE plot of the model conductances, where we reduced each models' conductance set for visualization in a two-dimensional space (Hinton and Roweis, 2002).The clusters in this plot indicate that there are groups of models with similar underlying conductance relationships in the model database.
The models in the database were constrained to produce solutions whose features were similar to the experimental measurements in the STG, but their activities are not identical.
For example, Figures 2D and 2E show the distributions of two features: the duty cycle of the LP cell, and the phase of the cycle at which the LP cell begins to fire.Interestingly these features are normally distributed (p > 0.05).Similarly to the values of the conductances, the features we extracted from the activity also can be matched to a high dimensional vector with a total number of features of N F = 11.Figure 2F shows a t-SNE plot of the models' features.In this case it is not possible to identify clusters, suggesting that there are no underlying relationships between the features in the model database.Even though the conductances are broadly distributed over many-fold ranges, the output of the circuits result in tight yet uncorrelated distributions.
C. Model networks with similar behavior and different underlying currents.
Computational and experimental studies suggest that neurons can produce similar activities employing different contributions of each ionic current.This observation is known as degeneracy and one way to illustrate this was introduced recently in Alonso and Marder (2019).A currentscape is a representation of neural dynamics that displays the percent con-Fournier Z, Alonso LM, and Eve Marder tribution of each ionic current to the total inward and outward currents over time.Figure 3 shows the currentscapes for the PD, LP, and PY neurons in two different models.We compared the robustness of intrinsic and synaptic conductances by modifying these two sets separately.For this, we independently varied the values of each maximal conductance ḡi by some random amount within a range (Giovannini et al., 2013;Zang and Marder, 2023).For each conductance, this variation range is defined as a percent deviation from its starting or control value.The larger this percentage, the larger the change in maximal conductances.The values of the modified conductances are obtained by sampling from a uniform distribution over each range.When the range is small and there is a small percent deviation from the control values, the expectation is that the activity will change very little and stay pyloric in most cases.As the variation range is increased, corresponding to larger Fournier Z, Alonso LM, and Eve Marder percent deviations from control values, it is reasonable to expect that the activity will be disrupted in a greater proportion of cases.
We explored changes to maximal conductances between 0 to 100 percent.The change size or variation range was defined as δ ∈ [0, 1].For each conductance we perform the change ḡi → ḡi * γ with γ = 1 + δU [−1, 1], a random number drawn from a uniform distribution and scaled by δ.In this way values of δ ≈ 0 correspond to small percent changes and the limit case δ = 1 allows for changes of 100% of the value of a given conductance.For each value of δ we sampled N = 1000 sets of conductances.For each set of conductances we simulated the model, extracted features, and classified the activity as pyloric or not pyloric using a random forest classifier (Ho, 1995).We then computed the percentage of traces that displayed pyloric activity, and evaluated how this percentage changes as the variation range δ increases.This procedure was performed independently for both intrinsic and synaptic conductances, allowing us to compare the sensitivity to changes of each set.We repeated the procedure for all one hundred models in our database.
Figure 4 shows the result of this analysis on ten models.In all panels, the x-axis corresponds to the variation range δ.The y-axis indicates the fraction of traces that were classified as pyloric activity.By definition, when the variation range is zero (δ = 0), all "perturbations" result in the same model, and therefore, the activity is the same as in control.As expected, the fraction of traces that display pyloric activity decreases as the variation range (δ) increases.This is true for both intrinsic (blue) and synaptic conductances (red).The synaptic sensitivity curve is above the intrinsic sensitivity curve for all values of the variation range δ.This means that the models are more resilient to changes in their synaptic conductances.
To quantify the observation that the intrinsic conductances are more sensitive to perturbations than the synaptic conductances, we calculated the average sensitivity curves in The sigmoid functions approximate the sensitivity curves well, with an average correlation coefficient of r 2 = 0.98, and an average error of approximately 1% for each value of δ.The fits to the synaptic sensitivity curves were marginally better than in the intrinsic case as the mean error for the synaptic was 0.55% (r 2 = 0.98), whereas it was 1.54% (r 2 = 0.98) for the intrinsic.
To compare the intrinsic and synaptic sensitivity curves, we inspected the distributions of the sigmoidal fit parameters (Figure 5C, 5D).These sets of parameters were not normally distributed (p < 0.05) except for the midpoint parameter in the intrinsic curves (p = 0.13).
To evaluate if two distributions were statistically different, we used Student's t-Test (Student, 1908).The midpoint parameter distributions are different, and there is little to no overlap between them (p < 0.001).In addition, we found that the intrinsic conductances have much smaller midpoint parameters, consistent with the observation that their sensitivity curves decrease faster as the size of the conductance changes δ increases.
Although the distributions of midpoint parameters are visibly different, the width parameter distributions largely overlap.However, these distributions are significantly different (p < 0.001) because the red distribution is long-tailed.A smaller midpoint parameter and a smaller width parameter both indicate that a curve decreases more quickly.This difference in the midpoint parameters indicates that the synaptic conductances must be varied further from their starting values to disrupt pyloric activity.
To further compare the intrinsic and synaptic sensitivity curves, we calculated the areas under the curves (Figure 5E).The areas under the synaptic curves are significantly larger (p < 0.001), consistent with the previous observation that the network is more robust to perturbations in their synaptic conductances.
To characterize the variability in the sensitivity curves, we performed principal component analysis (PCA).The first component captures 97.6% of the variance, which indicates that the curves are similar and predominantly vary along a single dimension.The distribution of the eigenvalues required to reconstruct the curves from the principal component are shown Fournier Z, Alonso LM, and Eve Marder in Figure 5F.The eigenvalues for the intrinsic and synaptic curves result in visibly different distributions (p < 0.05) with means of opposite signs.These eigenvalues are normally distributed in the case of the intrinsic curves (p = 0.82), and not normally distributed in the synaptic case (p < 0.05).The intrinsic and synaptic curves can be reasonably approximated by the same sigmoid-looking curve, but flipping the sign in each case.

E. Relative robustness of models visualized in conductance space
These 3 neuron models require the specification of N i = 27 intrinsic conductances and N s = 7 synaptic conductances.There are multiple combinations of ḡi that result in a proper pyloric rhythm.This means that there are regions in conductance space G for which the activity is pyloric, G P .The shape and properties of G P have been explored and studied in the past (Taylor et al., 2006).To gain some intuition about the properties of these regions we exploited the trained classifier, and evaluated its output over 2D ranges of conductance space.
Figure 6 shows 2D projections of the regions of conductance space for which the activity is pyloric G P , for two different models.The top rows show the projection of G P over pairs of intrinsic conductances, ḡNa and ḡKd .The blue regions correspond to values of the conductances for which the activity is pyloric (G P ), and the white cross indicates the "control" or initial values of the conductances for that model.The white box indicates the maximum range of the perturbations used in the sensitivity analyses (δ = 1).
The shape of G P is intricate and complex.Despite this, some basic properties of G P are revealed.The top panels of Figure 6 show the projection of G P over the same pair of intrinsic conductances, across cell types.The blue regions across panels are different, indicating that the same perturbation in different cell types may result in drastically different outputs: cell identity matters.The white cross is closer to the edge in the case of the LP cell, meaning that this particular model is more sensitive to perturbations in the LP cell than to perturbations in the other cells.
The second row in Figure 6, model A, shows the projection of G P over pairs of synaptic conductances.The red areas correspond to values of the conductances for which the activity is pyloric (G P ), and the white cross indicates the starting or "control" values of this particular model.As before the shape of these projections are intricate and complex.They also are different across pairs of conductances.Notice that with the exception of the cholinergic synapses (first panel), the initial value of the conductances (white cross) is farther away from the edge than in the intrinsic cases, consistent with the finding that the model is more robust to perturbations in its synaptic conductances.In the same vein, note that the proportion of area of the white box (δ = 1) occupied by the red (synaptic) region is visibly larger than the proportion of area occupied by the blue (intrinsic) region.
Figure 6 (bottom) shows the same calculation in a different model (model B).As before, the blue regions and the red regions are different across cells and pairs of synapses.The shape of the projections of G P are different across models.In the case of the intrinsic conductances the shapes are similar overall, but there are visible differences, such as the white cross in the LP cell being farther away from the edge as in model A. The projection of G P on pairs of synaptic conductances is very different in this model.In this case, almost the totality of the area enclosed by the white boxes (δ = 1) corresponds to G P , suggesting that model B is more robust than model A to changes in synaptic conductances.This is consistent with the finding that the synaptic sensitivity curve is steeper in model A (Model #1 in Figure 4) than in model B (Model #2 in Figure 4).These results provide an intuition of the complexities of G P .Not only are these regions hard-to-impossible to characterize in one circuit, but they are also different across circuits.

III. Discussion
A major goal of much of neuroscience is to understand how experience alters the function of brain circuits, both in the context of development and learning.At the same time, we are aware that damage to the brain can produce altered brain structure.It is therefore important to understand the resilience of brain circuits to damage, and to understand the plasticity mechanisms that foster learning and recovery from damage.Thus, it is critical to understand which changes in circuit parameters are most influential in altering circuit function.
Fournier Z, Alonso LM, and Eve Marder Despite the enormous body of work on the importance of alterations of synaptic strength in numerous preparations, there are also many well-documented cases of neurons that change their intrinsic excitability during development or as a function of neuromodulation.For example the Rohan-Beard cells first show a slow Ca 2+ -dependent spike, before eventually acquiring a rapid N a + -dependent action potential (Baccaglini and Spitzer, 1977) and ACh and biogenic amine neuromodulators directly alter the intrinsic properties of thalamic and other neurons (McCormick and Pape, 1990;McCormick et al., 1991;McCormick and Prince, 1986, 1987, 1988;McCormick and Wang, 1991).There are myriad other examples in which neuromodulators completely transform the dynamics of neuronal activity (Bargmann, 2012;Marder, 2012;Nusbaum and Blitz, 2012).A neuromodulator that changes a neuron's activity profile from silent, to tonic firing, to producing plateau potentials or oscillatory bursts of activity is likely to appreciably alter the activity of the circuits in which these neurons are found.Because altering an intrinsic current is often accompanied by changes in the time scales of neuronal activity, even modest changes in an intrinsic conductance can result in major alterations in the dynamics of a target neuron and the circuits in which it functions.In a sense, we biased the results of this study in manner that underestimates the influence of intrinsic activity by selecting for neurons that were within a restricted region of the degenerate range of solutions.
We expect that results similar to these will be found in many brain regions that show strong propensities for plateau and oscillatory properties, especially in response to neuromodulation.It remains to be seen whether cortical or hippocampal regions, typically thought to be modified by changes in synaptic strength (Hasselmo et al., 1995), are less sensitive to modulation of intrinsic properties.
In this study we evaluated the role of changes in inhibitory synaptic strengths, as in the STG synaptic connections are inhibitory.Because increases in inhibitory synaptic conductances push the membrane potential closer to reversal potential, there can be a diminishing effect of increasing the inhibitory synaptic strength, although the increase in the post-synaptic conductance still decreases the efficacy of all other currents, and changes in inhibitory synaptic currents can saturate (Prinz et al., 2003b).
Fournier Z, Alonso LM, and Eve Marder Although increasing the strength of excitatory synaptic currents also is somewhat diminished as the membrane potential of neurons moves closer to the reversal potential of an EPSP, the effect of increasing the EPSP may be amplified as voltage dependent currents are activated.Therefore, it is possible that a circuit dominated by excitation may be more dominated by changes in synaptic than intrinsic currents.Nonetheless, a large body of work now argues that much cortical function depends on a balance of excitatory/inhibitory function (Atallah and Scanziani, 2009;Isaacson and Scanziani, 2011;Lim and Goldman, 2013;Sadeh and Clopath, 2021), so a similar study such as this needs to be done in a modeled cortical or hippocampal circuit.That said, this study benefits greatly from a simple measure of circuit performance that can be used to determine circuit robustness, and this is less obvious in many other circuits.

A. The Model
We employed a model circuit inspired by the STG of crustaceans that was developed and studied previously by multiple authors (Golowasch et al., 1999;Prinz et al., 2004;Alonso and Marder, 2020).The model network is composed of three cells named PD, LP and PY, each modeled by a single compartment with nine currents: all cells have a sodium current, I N a ; a delayed rectifier potassium current, I Kd ; transient and slower calcium currents, I CaT and I CaS ; a transient potassium current, I A ; a calcium-dependent potassium current, I KCa ; a hyperpolarization-activated inward current, I H ; a leak current, I leak ; and a neuromodulatory current I M I (Golowasch and Marder, 1992b;Swensen and Marder, 2000).This model of the (1) The terms in the sum are the ionic currents and I e is injected current.
The maximal conductance of a channel is given by ḡi .The variables m i and h i correspond to the activation and inactivation processes.The integers p i and q i are the number of gates in each channel.The reversal potential of the ion associated with the i − th current is given by E i .
The reversal potential of the Na, K, H, leak, and I M I currents were kept fixed at The functions τ m i (V ), m ∞ i (V ), τ h i (V ) and h ∞ i (V ) are based on the experimental work of Turrigiano et al. (1995).The activation functions of the K Ca current require a measure of the internal calcium concentration [Ca +2 ] ( Liu et al., 1998).The dynamics of the intracellular calcium concentration are given by, Here Ca F = 0.94 µM nA is a current-to-concentration factor and Ca 0 = 0.05 µM .The total capacitance of the cell is C = 1nF .
The interactions in the network consist of 7 chemical synapses and are similar to Prinz et al. (2004).The synaptic current is given by I s = ḡs s(V post − E s ), where ḡs is the synaptic strength, V post is the membrane potential of the postsynaptic neuron and E s is the reversal potential of the synapse.The activation of a synapse s(t) is given by and These equations are identical to Prinz et al. (2004) except for the inclusion of a bound for the timescale of activation τ r = 20msec (Alonso and Marder, 2020).
B. Generating the model database.
We generated 100 different models of the pyloric network following the approach outlined in Alonso andMarder (2019, 2020), which we briefly review here.Each model is specified by a set of maximal conductances, while all other parameters (such as reversal potentials) are the same across models.Our approach consists of defining a target or objective function such that upon minimization, it will yield sets of conductances that result in normal or control pyloric activity.
The target function takes values of conductances as its input.The model is integrated for 20 seconds and the first 10 seconds are dropped to minimize transient activity.We then compute several quantities: the average (<>) burst frequency < f b >, the average duty cycle < dc >, the number of crossings with a slow wave threshold # sw = −50mV , the number of bursts # b , and the average lags between bursts < ∆ P D−LP > and < ∆ P D−P Y >.In this study we also sought solutions that would display a large voltage separation between the slow wave and the spike troughs.This was done by defining a second slow wave threshold at # sw = −45mV and checking that the number of crossings is the same as the number of bursts.Finally, to penalize unstable solutions we compute the standard deviation of the burst frequency and duty cycle distributions; a solution was penalized if std({f b }) ≥< f b > ×0.1 or std({dc}) ≥< dc > ×0.2.If a solution is not discarded we can use these quantities to score its similarity to the target, pyloric behavior.We define the errors as follows, Fournier Z, Alonso LM, and Eve Marder E f measures how far the bursting frequency of each cell is to a target frequency f tg = 3, E dc accounts for the duty cycle, E sw accounts for the difference between the number of bursts and the number of crossings with the slow wave threshold at E ph compares the lags between bursts (in units of the bursting period τ b = 1 <f b > ) to a target lag ∆ tg .These measures are discussed in more detail in Alonso andMarder (2019, 2020).
For each set of maximal conductances ḡi we can then define an objective function, The coefficients (α, β, γ, γ 2 , η) determine which features are more important.In this work we used α = 10, β = 1000, γ = 1 , γ 2 = 1, η = 10.We found that using these weights results in a target function that produces useful sets of conductances (that produce normal pyloric triphasic activity).All penalties E i were calculated using T = 10 secs.The simulations were done using an RK4 integration routine with dt = 0.05 msecs.
The target control behavior is defined as all cells bursting with 20% duty cycle for P D (dc tg P D = 0.2) and 25% for the LP and PY cells (dc tg LP,P Y = 0.25).The lag between bursts was targeted to be ∆ P D−LP tg = 0.5 and ∆ P D−P Y tg = 0.75.The target burst frequency of all cells was set to f tg = 1Hz (Bucher et al., 2006;Tang et al., 2010;Hamood et al., 2015).
Minimization of the target function successfully produces sets of maximal conductances g i that produce normal or control triphasic rhythms.
We defined a search space of allowed values listed here: for each cell we searched for ḡNa ∈ ḡH ∈ [0, 10 1 ], ḡL ∈ [0, 1], ḡMI ∈ [0, 10].All synaptic conductances were searched in the range ḡsyn ∈ [0, 1] ([µS]).We discretized the search space by taking 10 3 equally spaced values for each parameter.The objective function (9) was minimized using a custom genetic algorithm (Holland, 1992) on a desktop computer with 32 cores.We ran the genetic algorithm using 1000 initial individuals taken at random in the search space for 300 generations.Finding a successful population of models takes several hours.We repeated this process to produce Fournier Z, Alonso LM, and Eve Marder 100 different models with values of conductances distributed over large ranges.
C. Extracting features from model traces.
The models were simulated for 10 seconds, and the final 5 seconds of the simulation were utilized.We calculated spike times by finding the local maxima of the time series of the membrane potential of each cell.Some subthreshold depolarizations resulted in local maxima that were not spikes.To isolate the spikes from these subthreshold depolarizations, local maxima were filtered to only include spike times within 5ms of the local maxima of the derivative of the voltage dV dt .This is because local maxima of the voltage derivative will typically be within 5ms of the local maxima of the voltage if the rise in voltage is very rapid, which is the case for spikes produced by rapid sodium activation.In addition, we computed the maximum value of voltage in the 5 seconds of simulation that we analyzed (V max ), and we only kept spikes that were less than 25mV below this value.To prevent spike times from being detected in non-firing cells with fast but small oscillations, we imposed that if V max < −25mV , then the spike threshold was set to −25mV .Models with cells that fire spikes beneath −25mV do exist but were not included in this work.
To calculate the bursts' start and end times, we set a threshold to distinguish inter-burst intervals (IBIs) from inter-spike intervals (ISIs) ISI th .In the case of a regularly bursting cell, the ISI distribution is bimodal, with one peak corresponding to the intraburst ISI, and the other peak corresponding to the IBIs.We calculated the ISI threshold by first detecting the peaks in the histograms of ISIs and then chose the value halfway between these peaks.
In some cases the distribution of ISIs was not bimodal due to irregular bursting or tonic spiking.This made it difficult to determine an ISI threshold to correctly identify bursts.
In cases where the threshold could not be estimated, or it was estimated to be larger than 300mS or smaller than 100mS, we set ISI th = 100mS.This allowed us to reliably calculate phases by counting tonic spikes in the triphasic rhythm as burst starts and burst ends while also detecting spikes with large ISIs at the end of a burst.Using this method led to errors in burst identification in some cases, and this is a source of error and a limitation of the study.
The features extracted from each model include the mean and standard deviations of the distributions of duty cycles that each cell produces over the period of 5 seconds.We also Fournier Z, Alonso LM, and Eve Marder calculated the means of the distributions of the phases of the rhythm at which each cell begins and ends a burst.We calculated the duty cycles for each cell as the duration of the burst B d divided by the bursting period (τ b ), The phases were calculated as the difference between the cell's on/off times and the start of burst in the PD cell P D on (taken as a reference) divided by the period of the rhythm (τ N ), We defined τ n as the bursting period τ b of the PD cell, thus when PD either did not fire or fired irregularly phases were not calculated and these quantities were set to 0.
D. Classification of model traces.
The activity of each trace was evaluated by first building a vector containing features of its activity.The features are: the means of the distributions of duty cycles for each cell, the standard deviations of the distributions of duty cycles for each cell, and the means of the distributions of phases for each cell.We attempted classification of the vectors using a number of algorithms including random forest classifiers and gradient boosting machines (Ho, 1995;Friedman, 2001), and we chose to use random forest classifiers to avoid the possibility that gradient boosting machines might over fit to the training data for our purposes (Natekin and Knoll, 2013).A random forest classifier is a machine learning method used to classify items based on a list of the item's features (Ho, 1995).The classifier we created consisted of 500 decision trees (max depth = 4), where each tree was fed a different subset of the features for classification.A winner-takes-all vote between the trees chooses the output label of the classifier.The classifier must first be trained with a set of data that was manually labeled, with the help of experimentalists that study the pyloric network in the crustacean STG.
When the cells were bursting in a triphasic pattern in the order PD-LP-PY, with at least two spikes per burst, we considered the activity to be pyloric (see Figure 1).We trained the Fournier Z, Alonso LM, and Eve Marder classifier on a database of 1200 model solutions or traces featuring a wide variety of waveforms.Any other type of activity was considered 'not pyloric'.In some cases we excluded rhythms that displayed pyloric characteristics but also exhibited irregular features, such as abnormally large duty cycles, cells with low amplitude spikes, and inconsistent firing and/or bursting patterns.
We found that the trained classifier generalized with ≥ 95% accuracy for many, but not all models.To verify that the classifier properly classified the models we chose, we manually checked thousands of traces.Only models that were classified with this level of accuracy after perturbation were used in this study.On rare occasions, the classifier misclassified models that are clearly pyloric as not pyloric even though the list of features fed to the classifier were calculated properly.The percentage of misclassifications was approximately 2 − 3%, which is consistent with the accuracy of the classifier.

E. Sigmoidal fits
To analyze the sensitivity curves quantitatively, we fit them with sigmoid functions.We optimized for the sigmoidal midpoint and width parameters by minimizing an error function using a brute force approach.The sigmoid function is given by, where the a parameter controls the width of the sigmoid, the b parameter determines the midpoint of the sigmoid, and the c parameter corresponds to the maximum asymptote of the sigmoid function.Variable δ represents the variation range parameter.Parameter c was fixed at 100% because this is the maximum percentage of models that can be pyloric, and by definition all curves start at 100% when the size of the perturbation is δ = 0 (i.e.we impose S(δ = 0) = 100).To find the values of a and b that best approximate the sensitivity curves we minimize an error function defined as follows, where delta i = [0, 0.1, ....1] are the values of δ for which we computed the sensitivity The voltage traces of each network are plotted above the currentscapes.The total outward and inward current of each cell are shown as black-filled curves above and below each currentscape plot on a logarithmic scale.The dashed lines correspond to reference current values.The percent contributions of each ionic current are represented in colors.The outward currents are on top and inward currents on the bottom.

Figure 3
Figure 3 shows the currentscapes of each cell in two model networks.The currentscapes show visible differences in each cell within a model as well as across the two different models, yet the activity of the networks is similar.For example, there is a visibly larger contribution of the I M I current to the PD cell in Model B as compared to Model A. Likewise, the contribution of the I Kd current is more prominent in Model B's LP cell than Model A's LP cell.In addition, there is little to no contribution of the KCa current in the Model A PY cell, while contributions of this current are more visible in the Model B PY cell.Furthermore, there is a small contribution of the I CaS current in the Model A PD cell, but a larger contribution in the Model B PD cell.Because the activity is driven by different combinations of currents, the excitability properties of the cells are different, and this in turn results in differential sensitivities to perturbations.D. Model networks are more sensitive to perturbation of intrinsic than synaptic conductances.

Figure 5A .
Figure 5A.The plot shows the average fraction of pyloric traces across all models for all values of the variation range δ.To further quantify these observations we fitted each sensitivity curve with a sigmoid function (Figure 5B), S(δ) = c • 1 + pyloric network is identical to previous ones, except for the inclusion of the neuromodulatory current I M I .The membrane potential of the cells evolves according to a Hodgkin-Huxley equation for a single compartment, , and E IM I = −10mV .The reversal potential of the calcium currents E Ca was updated using the Nernst equation assuming an extracellular calcium concentration of 3×10 3 µM .The kinetic equations describing the seven voltage-gated conductances were modeled as in Liu et al. (1998),

FIG. 1 .
FIG. 1. Feature extraction to classify model network activity (A) Schematic of the model used in Prinz et al. (2004).The three cells are connected by seven inhibitory chemical synapses of two types: the red curves are cholinergic synapses and the grey curves are glutamatergic synapses.(B) Pyloric network activity.Each cell displays periodic bursting activity and the cells fire in a sequence PD-LP-PY.The color lines indicate the measures used for calculating duty cycle and firing phases.(C) Membrane voltages of control network activity (in black), overlaid with the network activity that results from changing the intrinsic conductances (blue).(D) Membrane voltages of control network activity (in black) overlaid with network activity that results from changing the synaptic conductances (in red).(E) Changes in intrinsic maximal conductance values.The black dots show the unperturbed or control conductances and the blue dots correspond to the changed conductances that produced the blue voltage traces in Figure 1C.Only 2 maximal conductances per cell are shown for simplicity, but all 9 × 3 = 27 intrinsic conductances were changed.(F) Synaptic maximal conductance values for the control network (in black) and the perturbed network (in red) in Figure 1D.

FIG. 2 .
FIG. 2. A model database of degenerate solutions.(A, B) Distributions of gA maximal conductance in PD cells and P D − P Y glut glutamatergic synapse maximal conductance (N = 100).(C) Dimensionality reduction (t-SNE) and visualization of models' maximal conductances.(D, E) Distributions of features: LP duty cycle and LP-ON phase.(F)Dimensionality reduction (t-SNE) and visualization of models' conductances models' features