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.

II. Results

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_MI, 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_Na; 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 hyperpolarization- activated inward current, I_H ; a leak current, I_leak; and a neuromodulatory current I_MI.

FIG. 1.

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 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 . 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 and were changed, but all intrinsic conductances were changed by a similar percent. In the LP cell, the sodium conductance decreased while the delayed-rectifier potassium conductance increased, leading to a loss of spiking activity. The sodium maximal conductance in the PY cell also decreased but decreased even further, thus maintaining spiking activity. These effects are due to the relative magnitude of and 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 , 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, 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_MI 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).

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 , 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.

FIG. 2.

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 contribution 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. 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.

FIG. 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_MI 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

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 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 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 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.

FIG. 4.

To quantify the observation that the intrinsic conductances are more sensitive to perturbations than the synaptic conductances, we calculated the average sensitivity curves in 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), , that best approximates it. The parameters are a, which controls the width, and b, which determines the sharpness. Parameter c controls the overall amplitude of the sigmoid and was kept fixed at c = 100, since by definition it must be S(0) = 100.

FIG. 5.

The sigmoid functions approximate the sensitivity curves well, with an average correlation coefficient of r² = 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² = 0.98), whereas it was 1.54% (r² = 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 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 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, and . 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).

FIG. 6.

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.

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²⁺-dependent spike, before eventually acquiring a rapid Na⁺-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 neuro-modulation. 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).

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.

IV. Methods

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_Na; 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_MI (Golowasch and Marder, 1992b; Swensen and Marder, 2000). This model of the pyloric network is identical to previous ones, except for the inclusion of the neuromodulatory current I_MI.

The membrane potential of the cells evolves according to a Hodgkin-Huxley equation for a single compartment,

The terms in the sum are the ionic currents and I_e is injected current. The maximal conductance of a channel is given by . 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_MI currents were kept fixed at E_Na = 30mV, E_K = −80mV, E_H = −20mV, E_leak = −50mV, and E_IMI = −10mV. The reversal potential of the calcium currents E_Ca was updated using the Nernst equation assuming an extracellular calcium concentration of 3×10³µM. The kinetic equations describing the seven voltage-gated conductances were modeled as in Liu et al. (1998),

The functions and 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⁺²] (Liu et al., 1998). The dynamics of the intracellular calcium concentration are given by,

Here is a current-to-concentration factor and Ca₀ = 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 , where 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

with,

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 and Marder (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 < Δ_PD−LP > and < Δ_PD−PY >. 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,

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 t_sw = −50mV (if #_sw ≠#_b then E_sw is large), E_sw2 is the same as E_sw but for a different threshold value (t_sw2 = −45mV), E_ph compares the lags between bursts (in units of the bursting period to a target lag Δ_tg. These measures are discussed in more detail in Alonso and Marder (2019, 2020).

For each set of maximal conductances we can then define an objective function,

The coefficients (α, β, γ, γ₂, η) determine which features are more important. In this work we used α = 10, β = 1000, γ = 1, γ₂ = 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 PD and 25% for the LP and PY cells . The lag between bursts was targeted to be Δ_{PD−LP tg} = 0.5 and Δ_{PD−PY 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 . All synaptic conductances were searched in the range . We discretized the search space by taking 10³ 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 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 . 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 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 PD_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 classifier on a database of 1200 model solutions or traces featuring a wide variety of wave-forms. 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 curves (N_δ = 11), and y(δ_i) are the values of a sensitivity curve for each value of δ_i. We evaluated all possible combinations of the a and b parameters over a 100 × 100 2D grid of possible pairs with (−0.5 to −0.01 for a, 0 to 2 for b), and kept the values of a and b that minimized the error function E.

F. Conductance Ranges

There are regions in conductance space where the activity is a normal pyloric rhythm. Visualizing and characterizing these regions is challenging because they are high dimensional objects, and we are limited to visualizing their two-dimensional projections. We evaluated the output of the classifier over 2-dimensional ranges in conductance space. For this we chose a pair of conductances and picked 100 values between 0 and for each pair. We then evaluated network activity using the classifier on the 10⁴ possible pairs of over this range. The values of that produce pyloric activity were plotted in color (blue, red) and the values of that result in non-pyloric activity were plotted in grey. For the intrinsic conductances, we chose to visualize the and maximal conductances projection for each cell. For the synaptic conductances, we visualized three pairs: the cholinergic inputs to LP and PY from PD, the glutamatergic inputs to LP from PD and PY, and the glutamatergic inputs from PD to LP and LP to PD. We performed this procedure on two different models. For Model A, for and was 250nS. In the case of the synaptic conductances was set to 0.45nS on the x-axis and 1.75nS on the y-axis. For Model B, was 200nS for and 175nS for for the synaptic conductances was 0.15nS on the x-axis and 1.75nS on the y-axis.

Funding acknowledgments

NIH R35 NS097343, R01 MH046742 (E.M.)