Mapping circuit dynamics during function and dysfunction
- Volen Center and Biology Department, Brandeis University, United States
- Federated Department of Biological Sciences, New Jersey Institute of Technology and Rutgers University, United States
Abstract
Neural circuits can generate many spike patterns, but only some are functional. The study of how circuits generate and maintain functional dynamics is hindered by a poverty of description of circuit dynamics across functional and dysfunctional states. For example, although the regular oscillation of a central pattern generator is well characterized by its frequency and the phase relationships between its neurons, these metrics are ineffective descriptors of the irregular and aperiodic dynamics that circuits can generate under perturbation or in disease states. By recording the circuit dynamics of the well-studied pyloric circuit in Cancer borealis, we used statistical features of spike times from neurons in the circuit to visualize the spike patterns generated by this circuit under a variety of conditions. This approach captures both the variability of functional rhythms and the diversity of atypical dynamics in a single map. Clusters in the map identify qualitatively different spike patterns hinting at different dynamic states in the circuit. State probability and the statistics of the transitions between states varied with environmental perturbations, removal of descending neuromodulatory inputs, and the addition of exogenous neuromodulators. This analysis reveals strong mechanistically interpretable links between complex changes in the collective behavior of a neural circuit and specific experimental manipulations, and can constrain hypotheses of how circuits generate functional dynamics despite variability in circuit architecture and environmental perturbations.
Editor's evaluation
This study applies an unsupervised dimensionality reduction (t-SNE) to characterize neural spiking dynamics in the pyloric circuit in the stomatogastric ganglion of the crab, an important system for mechanistic analysis of rhythmic circuit function. The application of unsupervised methods to characterize qualitatively distinct regimes of spiking neural circuits is interesting and novel. The challenges and lessons learned in this study are of broader interest to those seeking to quantitatively characterize large sets of neural data across many subjects. The method is demonstrated across hundreds of animal subjects and used to investigate circuit responses to a variety of perturbations.
https://doi.org/10.7554/eLife.76579.sa0Introduction
Neural circuits can generate a wide variety of spiking dynamics, but must constrain their dynamics to function appropriately. Cortical circuits maintain irregular spiking patterns through a balance of excitatory and inhibitory inputs (van Vreeswijk and Sompolinsky, 1996; Mariño et al., 2005; Brunel and Wang, 2003) and the loss of canonical dynamics is associated with neural diseases like channelopathies and epilepsy (Marbán, 2002; Staley, 2015). Preserving functional dynamics can be a challenge for neural circuits for the following reasons. The same spike pattern can be generated by diverse circuits with many different topologies and broadly distributed synaptic and cellular parameters (Prinz et al., 2004; Golowasch et al., 2002; Alonso and Marder, 2019; Memmesheimer and Timme, 2006). Furthermore, neural circuits are constantly being reconfigured, with ion channel protein turnover, and homeostatic feedback mechanisms modifying conductance and synapse strengths continuously (Turrigiano et al., 1994; Turrigiano et al., 1995; O’Leary et al., 2014; Franci et al., 2020). The problem of maintaining functional activity patterns is aggravated by the fact that functional circuit dynamics tend to lie within a low-dimensional subspace within the high-dimensional state space: of the numerous possible solutions, only a few are functional and are found in animals (Cunningham and Yu, 2014; Pang et al., 2016). How do neural circuits preserve functional dynamics despite these obstacles?
Answering this question requires, as a prerequisite, a quantitative description of the dynamics of neural circuits during function and dysfunction. When rhythms are regular, this is relatively simple, but when rhythms become irregular, classifying them becomes hard (Haddad and Marder, 2018; Tang et al., 2012; Haley et al., 2018). In this article, we study the dynamics of a well-studied central pattern generator, the pyloric circuit in the stomatogastric ganglion (STG) in Cancer borealis (Marder and Bucher, 2007). The pyloric circuit is small, in crabs consisting of 13 neurons coupled by inhibitory and electrical synapses. Its topology and cellular dynamics are well understood, and the circuit generates a clearly defined ‘functional’ collective behavior where bursts of spikes from three different cell types alternate rhythmically to generate a triphasic motor pattern. The stereotypy and periodicity of the motor pattern suggest that the baseline dynamics of the pyloric circuit are fundamentally low dimensional. This has allowed for the effective parameterization of the rhythm by a small number of ad hoc descriptors such as the burst period, duty cycles, and phase of each neuron (Hartline and Maynard, 1975; Eisen and Marder, 1984; Miller and Selverston, 1982).
In response to perturbations that span many cycles, pyloric circuit dynamics are not always periodic, and descriptors that work well to characterize the canonical rhythm are inadequate to describe these atypical dynamic states. Efforts to study circuit dynamics under these regimes, and to characterize how the circuit responds to, and recovers from perturbations, have been frustrated by the inability to quantitatively describe irregular and non-stationary dynamics (Haddad and Marder, 2018; Tang et al., 2012; Haley et al., 2018).
In this article, we set out to address the problem of quantitatively describing neural circuit dynamics under a variety of conditions. We reasoned that circuit dynamics lie on some lower-dimensional set within the full high-dimensional space of possible dynamics, even when circuits exhibit atypical and nonfunctional behavior, because even circuits generating dysfunctional dynamics are still constrained by cellular parameters and network topology. We therefore set out to find and visualize this subset of spike patterns using an unsupervised machine learning approach to visualize patterns in the high-dimensional data in two dimensions. This method allows us to visualize the totality of a large and complex dataset of spike patterns, while being explicit about the assumptions and biases in the analysis. Using this method, we found nontrivial features in the distribution of the data that hinted at diverse, stereotyped responses to perturbations. Using this compact representation allowed us to efficiently manually classify these patterns and measure transitions between these patterns. We were thus able to characterize the diversity of circuit dynamics under baseline and perturbed conditions, and identify anecdotally observed atypical states within the full repertoire of spiking patterns for many hundreds of animals.
Results
Perturbations can destabilize the triphasic pyloric rhythm
Studies that measure the pyloric rhythm commonly involve recording from nerves from the STG in ex vivo preparations. Preparations typically also include the stomatogastric nerve (stn) that carries the axons of descending neuromodulatory neurons from the esophageal and commissural ganglia that project into the STG. Under baseline conditions (11°C, with the stn intact, Figure 1a), the periodic triphasic oscillation of the pyloric circuit can be measured by extracellular recordings of the lateral pyloric, pyloric dilator, and pyloric nerves (lpn, pdn, and pyn) (Figure 1a). Bursts of spikes from the pyloric dilator (PD) neurons on the pdn are followed by bursts of spikes from the lateral pyloric neuron (LP) on lpn and bursts of spikes from the pyloric neurons (PY) on pyn. Spikes from lateral posterior gastric (LPG) neurons are also found on the pyn nerve in these recordings and can be differentiated from PY spikes by their shape and timing (LPG is active during PD bursts). Under these control conditions, where the rhythm is robust and spikes from these neurons are easily identifiable both by their location on the nerve and their phase in the cycle, the dual problems of identifying spikes from raw extracellular recordings and meaningfully describing circuit dynamics are easily resolvable.
Figure 1
The triphasic pyloric rhythm can become irregular and hard to characterize under perturbation.
(a) Simplified schematic of part of the pyloric circuit (left). Filled circles indicate inhibitory synapses, solid lines are glutamatergic synapses, and dotted lines are cholinergic synapses. Resistor symbol indicates electrical coupling. The pyloric circuit is subject to descending neuromodulatory control from the stomatogastric nerve (stn). Right: simultaneous extracellular recordings from the lvn, lpn, pdn, and pyn motor nerves. Action potentials from lateral pyloric (LP), pyloric dilator (PD), and pyloric (PY) are visible on lpn, pdn, and pyn. Under these baseline conditions, PD, LP, and PY neurons burst in a triphasic pattern. The anterior burster (AB) neuron is an endogenous burster and is electrically coupled to PD neurons. (b) When the stn is cut, neuromodulatory input is removed and the circuit is ‘decentralized.’ In this case, the pyloric rhythm can become irregular and hard to characterize. In addition, spikes from multiple PY neurons can become harder to reliably identify on pyn.
In studies that characterize the changes in circuit dynamics to prolonged perturbations, spike identification and circuit dynamics characterization are less straightforward. For example, when descending neuromodulatory projections from the stn are cut (i.e., when the STG is decentralized, Figure 1b), the collective dynamics of the pyloric circuit can become less regular. This loss of regularity is concomitant with spikes being harder to reliably identify in extracellular recordings. While PD and LP neuron spikes can still be typically easily identified on the pdn and lpn nerves (Figure 1b), identifying PY on the pyn in the absence of a regular rhythm can be challenging. This problem is aggravated by the fact that spikes from the LPG neuron are frequently found on pyn, and because there are several copies of the PY neuron, whose spikes can range from perfect coincidence to slight offsets that can unpredictably change the amplitude and shape of PY spikes due to partial summation. For these reasons, some previous works studying the response of pyloric circuits to perturbations have consistently recorded from the lpn and pdn nerves, but not from the pyn (Hamood et al., 2015; Haley et al., 2018; Haddad and Marder, 2018; Rosenbaum and Marder, 2018). Therefore, in order to include the largest number of experiments in our analysis, we chose to characterize the dynamics of the LP and PD neurons.
Nonlinear dimensionality reduction allows for the visualization of diverse pyloric circuit dynamics
The regular pyloric rhythm involves out-of-phase bursts of spikes between LP and PD, and is observed under baseline conditions (Figure 2a1-3). Perturbations such as the removal of descending neuromodulatory inputs, changes in temperature, or changes in pH can qualitatively alter the rhythm, leading to a large variety of hard-to-characterize spiking patterns (Figure 2a4-6). Because these irregular states may lose the strong periodicity found in the canonical motor pattern, burst metrics such as burst period or phase offsets between bursts that work well to characterize the regular rhythm perform poorly. Efforts to characterize and quantify these atypical spike patterns must overcome the slow timescales in observed dynamics, the large quantity of data, and irregularity and variability in observed spike trains. Previous work used ad hoc categorization systems to assign observations of spike trains into one of a few groups (Haddad and Marder, 2018; Haley et al., 2018), but these categorization methods scaled poorly and relied on subjective annotations.
Figure 2 with 5 supplements see all
Visualization of diverse neural circuit dynamics.
(a) Examples of canonical (1–3) and atypical (4–6) spike patterns of pyloric dilator (PD; blue) and lateral pyloric (LP; red) neurons. Rasters show 10 s of data. (b–d) Schematic of data analysis pipeline. (b) Spike rasters in (a2) can be equivalently represented by interspike intervals (ISIs) and phases. 20 s bins shown. Each 20 s bin contains a variable number of spikes/ISIs. (c) Summary statistics of ISI and phase sets in (d), showing tenth percentiles. Using percentiles converts the variable length sets in (b) to vectors of fixed length. (d) -scored data assembled into a single vector, together with some additional measures (Materials and methods). (e) Embedding of data matrix containing all vectors such as the one shown in (d) using t-distributed stochastic neighbor embedding (t-SNE). Each dot in this image corresponds to a single 20 s spike train from both LP and PD. Example spike patterns shown in (a) are highlighted in the map. points from animals. In (a–d), features derived from LP spike times are shown in red, and features derived from PD spike times are shown in blue.
We sought instead to visualize the totality of pyloric circuit dynamics under all conditions using a method that did not rely on a priori identification of (non)canonical dynamic patterns. Such a data visualization method, while descriptive, would generate a quantitative vocabulary to catalog the diversity of spike patterns observed both when these patterns were regular and also when they were irregular and aperiodic, thus allowing for the quantitative characterization of data previously inaccessible to traditional methods (Börner et al., 2005; Nguyen and Holmes, 2019).
The visualization was generated as follows: time-binned spike trains were converted into their equivalent interspike interval (ISI) and phase representations (Figure 2b, Materials and methods). For all analyses, we consider nonoverlapping 20 s time bins. We chose this time bin following inspection of circuit dynamics across many conditions in several animals. Because there can be an arbitrary number of spikes in a bin, there are an arbitrary number of ISIs and phases. This makes it challenging to find a basis to represent the entire dataset. Ideally, we want to represent the spike pattern in each 20 s bin with a point in some space of high but fixed dimensionality. To convert this into a vector of fixed length, we measured percentiles of ISIs and phases (Figure 2c). Together with other metrics (like ratios of ISIs, measures that capture discontinuities in ISI distributions, see Materials and methods for details), these percentiles were assembled into a fixed-length vector and each dimension was -scored across the entire dataset (Figure 2d). A collection of spike trains from an arbitrary number of neurons has thus been reduced to a matrix where each row consists of -scored percentiles of ISIs and other metrics. This matrix can be visualized using a nonlinear dimensionality reduction technique such as -distributed stochastic neighbor embedding (t-SNE) (Van der Maaten and Hinton, 2008), which can generate a two-dimensional representation of the full dataset (Figure 2e).
In this representation, each dot corresponds to a single time bin of spike trains from both neurons. We found by manual inspection that spike trains that are visually similar (Figure 2a1-3) tend to occur close to each other in the embedding (Figure 2e1-3). Spike patterns that are qualitatively different from each other (Figure 2a4-6) tended to occur far from each other, often in clusters separated by regions of low data density (Figure 2e4-6, Supplementary file 1).
How useful is such a visualization and does it represent the variation in spike patterns in the data in a reasonable manner? We colored each point by classically defined features such as the burst period or the phase (Figure 2—figure supplement 1). We found that the embedding arranges data so that differences between clusters and within clusters had interpretable differences in various burst metrics. For example, clusters on the left edge of the map tended not to have defined LP phases, typically due to silent or very sparse LP firing (Figure 2—figure supplement 1b). Location of data in the largest cluster was correlated to firing rate in the PD neuron (Figure 2—figure supplement 1c). We observed that burst metrics, when they were defined, tended to vary smoothly across the map. To quantify this observation, we built a Delaunay triangulation (Materials and methods) on the embedded data and measured the triadic differences between PD burst periods and PD duty cycles (Figure 2—figure supplement 3). Triadic differences in these metrics were significantly smaller in the map than triadic differences in a projection of the first two principal components or a shuffled map (p<0.0001, Kolmogorov–Smirnv test), suggesting that the t-SNE cost function generates a useful embedding where spike features vary smoothly within clusters. Finally, to validate our approach, we generated a synthetic dataset with different classes of spike patterns (Materials and methods) and analyzed it similarly. Coloring points in the t-SNE embedding by the original class revealed that clusters in the t-SNE map corresponded to different classes in the synthetic data, suggesting that this method can identify and recover stereotyped spike patterns in neural data (Figure 2—figure supplement 5).
Visualization of circuit dynamics allows manual labeling and clustering of data
Previous studies have shown that regular oscillatory bursting activity of the pyloric circuit can qualitatively change on perturbation. Circuit dynamics can be highly variable and has been categorized into various states such as ‘atypical firing,’ ‘LP-01 spikes,’ or ‘atypical’ (Haddad and Marder, 2018; Haley et al., 2018). Both the process of constructing these categories and the process of classifying data into these categories are typically done manually, and therefore requires expert knowledge that is not explicitly captured and is impossible to reproduce. Because the embedding distributed data into clusters, we hypothesized that clusters corresponded to stereotyped dynamics that were largely similar, and different clusters represented the qualitatively different circuit dynamics identified by earlier studies.
We therefore manually inspected circuit dynamics at randomly chosen points in each apparent cluster and generated labels to describe the dynamics in that region (Figure 3). This process colored the map and segmented it into distinct regions that broadly followed, and were largely determined by, the distribution of the data in the embedding (Figure 3a). Most of the data (57%) were assigned the regular label, where both PD and LP neurons burst regularly in alternation with at least two spikes per burst, and all identified regular states occurred in a single contiguous region in the map (blue). In the LP-weak-skipped state, PD bursts regularly, but LP does not burst every cycle, or only fires a single spike per burst. Irregular-bursting states showed bursting activity on both neurons, which were interrupted or otherwise irregular. In contrast, the irregular state showed spiking that was more variable and did not show strong signs of bursting at any point. LP-silent-PD-bursting states had regular bursting on PD, with no spikes on LP, while LP-silent states also had no spikes on LP, but activity on PD was more variable, and did not show regular bursting.
Figure 3 with 3 supplements see all
Map allows identification of distinct spiking dynamics.
(a) Map of all pyloric dynamics in dataset where each point is colored by manually assigned labels. Each point corresponds to a 20 s paired spike train from lateral pyloric (LP) and pyloric dilator (PD) neurons. Each panel in (b) shows two randomly chosen points from that class. The number of points in each class is shown in parentheses above each panel. points from animals. Labels are ordered by likelihood in the data.
The time evolution of the pyloric dynamics of every preparation constitutes a trajectory in the map, and every point in the map is therefore associated with an instantaneous speed of motion in the map. We hypothesized that instantaneous speed could vary across the map, with points labeled regular moving more slowly through the map than points with labels corresponding to atypical states such as irregular because regular rhythms would vary less over time. Consistent with this, we found that points in the regular cluster tended to have smaller speeds than points in other clusters (Figure 3—figure supplement 1a). Speeds in the regular state were significantly lower than every other state except PD-silent-LP-bursting (p<0.004, permutation test), suggesting that atypical states were associated with increased variability in circuit dynamics (Figure 3—figure supplement 1b).
Do the clusters we see in the data, and the resultant categorization of the data, depend strongly on the details of the dimensionality reduction method we used (t-SNE)? We used an entirely different embedding algorithm (Uniform Manifold Approximation and Projection [uMAP], McInnes et al., 2018) to embed the feature vectors in two-dimensional space. The map generated by uMAP preserved the coarse feature of the t-SNE embedding, suggesting that the features in the map reflected the features of the distribution of the data more strongly than details of the dimensionality reduction method. Coloring points in the uMAP embedding (Figure 3—figure supplement 3) revealed a roughly similar organization of data in the embedding space, suggesting that our categorization method did not strongly depend on the details of the dimensionality reduction.
Variability in baseline circuit dynamics across a population of wild-caught animals
Work on the pyloric circuit has used a wild-caught crustacean population. This uncontrolled environmental and genetic variability serves as a window into the extant variability of a functional neural circuit in a wild population of animals. In addition, experimental and computational work has shown that similar rhythms can be generated by a wide variety of circuit architectures and cellular parameters (Prinz et al., 2003; Hamood and Marder, 2014; Alonso and Marder, 2019). We therefore set out to study the variability in baseline circuit dynamics in the 346 pyloric circuits recorded under baseline conditions in this dataset.
The burst period of the pyloric circuit in the lobster can vary two- to threefold under baseline conditions at 11°C across animals (Bucher et al., 2005). Despite this sizable variation, other burst metrics, such as the phase onset of follower neurons, or the duty cycles of individual neurons, are tightly constrained (Bucher et al., 2005), likely related to the fact that these circuits are under activity-dependent feedback regulation (Turrigiano et al., 1995; O’Leary et al., 2014; Gorur-Shandilya et al., 2020) as they develop and grow. Activity-dependent regulation of diverse pyloric circuits could constrain variability in a single circuit across time to be smaller than variability across the population.
To test this hypothesis, we measured a number of burst metrics such as burst period and the phases and duty cycles of the two neurons across these 346 preparations in baseline conditions (Figure 4) when data are labeled regular because metrics are well-defined in this state. The mean values of each of these metrics were unimodally distributed (Figure 4a) and the coefficient of variation (CV) for all metrics was approximately 0.1 (Figure 4b). Using the mean CV in each individual as a proxy for the within-animal variability, and the CV of the individual means as a proxy for the across-animal variability, we found that every metric measured was more variable across animals than within animals (Figure 4c). Shuffling experimental labels generated null distributions for excess variability across animals and showed that across-animal variability was significantly greater than within-animal variability (Figure 4d, p<0.007, permutation test, Table 1).
Figure 4 with 3 supplements see all
Variability of burst metrics under baseline conditions.
(a) Variability of burst metrics in pyloric dilator (PD) and lateral pyloric (LP) neurons across a population of wild-caught animals. Metrics are only computed under baseline conditions and in the regular cluster. (b) Distribution of coefficient of variation (CV) of metrics in each animal across all data from that animal. In (a, b), each dot is from a single animal, and distributions show variability across the entire population. (c) Across-animal variability (CV of individual means, Δ) is greater than within-animal variation (mean of CV in each animal, Ο) for every metric. (d) Difference between across-animal variability and within-animal variability (colored dots). For each metric, gray dots and distribution show differences between across-animal and within-animal variability for shuffled data. points from animals.
Table 1
ANOVA results and power analysis for Figure 4.
ANOVA results for burst metrics in baseline conditions. For each metric, each animal is treated as a group and the variability (mean square difference) is compared within and across group. is the ratio of across-animal to within-animal mean square differences. N.99 is the estimate of the sample size required to reject the null hypothesis with a probability of 0.99 when the alternative hypothesis is true. animals.
| Metric | Across-animal MS | Within-animal MS | F | N.99 |
|---|---|---|---|---|
| LP delay off (s) | 1.1391 | 0.010 956 | 103.97 | 6 |
| LP delay on (s) | 0.616 47 | 0.0111 | 55.54 | 6 |
| LP durations (s) | 0.363 86 | 0.012 366 | 29.424 | 4 |
| LP duty cycle | 0.159 86 | 0.001 309 3 | 122.09 | 10 |
| LP phase off | 0.234 06 | 0.007 227 9 | 32.383 | 11 |
| LP phase on | 0.216 55 | 0.008 811 5 | 24.576 | 9 |
| PD burst period (s) | 3.557 | 0.036 872 | 96.469 | 4 |
| PD durations (s) | 0.079 397 | 0.000 549 44 | 144.5 | 6 |
| PD duty cycle | 0.053 472 | 0.000 413 23 | 129.4 | 16 |
-
LP: lateral pyloric; PD: pyloric dilator.MS: mean square.
It is reasonable to suppose that all baseline data exist in the regular cluster. While most baseline data are confined to the regular cluster (, Figure 4—figure supplement 1a), the remaining data, nominally recorded under baseline conditions, contains atypical circuit dynamics (Figure 4—figure supplement 1b and c). What causes these atypical circuit dynamics in this large, unbiased survey of baseline pyloric activity? One possibility could be inadvertent damage to the preparation caused by dissection and preparation of the circuit for recording. Consistent with this, we found that the probability of observing regular states was significantly reduced when cells were recorded from intracellularly (Figure 4—figure supplement 2), which may be due to increase in leak currents owing to impaling cells with sharp electrodes (Cymbalyuk et al., 2002) or due to cell dialysis (Hooper et al., 2015). No significant correlation was observed between sea surface temperatures (a proxy for environmental conditions for these wild-caught animals) and burst metrics (Figure 4—figure supplement 3a–c) or the probability of observing a regular state (Figure 4—figure supplement 3d). Taken together, these results underscore the importance of verifying that baseline or control data does not include uncontrolled technical variability that could mask biological effects of interest.
Perturbation modality alters state probability
The pyloric circuit and other circuits in the crab must exhibit robustness to the environmental perturbations that these animals are likely to encounter. Previous studies have characterized the ability of crustacean circuits to be robust to environmental perturbations such as pH (Haley et al., 2018; Ratliff et al., 2021; Qadri et al., 2007), temperature (Tang et al., 2010; Tang et al., 2012; Rinberg et al., 2013; Haddad and Marder, 2018; Kushinsky et al., 2019), oxygen levels (Clemens et al., 2001), and changes in extracellular ionic concentrations (He et al., 2020). Robustness to these perturbations exists up to a limit, likely reflecting the bounds of the natural variation in these quantities that these circuits are evolved to function in. When challenged with extremes of any of these perturbation modalities, the pyloric rhythm breaks down, displaying irregular or aberrant states, and may even cease spiking entirely. Such states are commonly referred to as ‘crashes’ and can have many flavors (Haddad and Marder, 2018; Tang et al., 2010; Tang et al., 2012) and involve the loss of the characteristic antiphase activity in the LP and PD neurons.
It remains unclear if extreme perturbations of different modalities share common pathways of destabilizing and disrupting the pyloric rhythm (Ratliff et al., 2021). In principle, these environmental perturbations can disrupt neuron and circuit function in qualitatively different ways: for example, changes in extracellular potassium concentration can alter the reversal potential of potassium (He et al., 2020) vs. changes in temperature can have varied effects on the timescales and conductances of all ion channels (Tang et al., 2010; Caplan et al., 2014). Because prior work was focused on studying the limits of robustness and lacked a detailed quantitative description of irregular behavior, the fine structure of the transition between functional dynamics and silent or ‘crashed’ states remain poorly characterized (Ratliff et al., 2021). We therefore set out to measure how pH, temperature, and extracellular potassium perturbations alter circuit state probability.
Where in the map are data under extreme environmental perturbations? Circuit spike patterns under high pH (>9.5), high temperature (>25°C), and high extracellular potassium () are distributed across a wide region of the map, spanning both regions in the regular cluster and other nonregular clusters (Figure 5a). Spike patterns observed under high-temperature conditions in the regular region were clustered in the lower extremity, in the region containing high firing rates and small burst periods of PD (Figure 2—figure supplement 1), consistent with earlier studies showing that elevated temperatures tend to speed up the pyloric rhythm (Tang et al., 2010; Tang et al., 2012).
Figure 5 with 1 supplement see all
Effect of three different environmental perturbations.
(a) Map showing regions that are more likely to contain data recorded under extreme environmental perturbations. (b) Treemaps showing probability distributions of states under baseline and perturbed conditions. (c) Probability distribution of states preceding silent state under perturbation. pH perturbations: from 6 animals; perturbations: from 20 animals; temperature perturbations: from 414 animals.
Subjecting the pyloric circuit to extremes of pH, temperature, and extracellular potassium altered the distribution of observed states (Figure 5c). In all cases, the probability of observing regular was significantly reduced (p<0.001, paired permutation test), and a variety of nonregular states were observed. We observed that high pH (>9.5) did not silence the preparation, but silent states were observed in low pH (<6.5), consistent with previously published manual annotation of this data (Haley et al., 2018). Silent states were also observed in , as reported earlier by He et al., 2020. Previous work has shown that the isolated pacemaker kernel (AB and PD neurons) has a stereotyped trajectory from bursting through tonic spiking to silence when subjected to pH perturbations (Ratliff et al., 2021), but moves through a different trajectory (bursting to weak bursting to silence) during temperature perturbations. Do pathways to silent states share similarities across perturbation modality in intact circuits? To answer this, we plotted the probability of observing states conditioned on the transition to silence in low pH, high temperature, and (Figure 5d). In the transitions between states detected, we never observed a transition from regular to silent, suggesting that the timescales of silencing are slow, longer than the width of one data bin (20 s). Trajectories to silent states always transition through a few intermediate states such as sparse-irregular, LP-silent, or PD-silent (Figure 5d).
Transitions between states during environmental perturbations
Changes in temperature, pH, and have different effects on the cells in the pyloric circuit and therefore can destabilize the rhythm in different ways. Increasing the extracellular changes the reversal potential of K+ ions, altering the currents flowing through potassium channels, and typically depolarizes the neuron (He et al., 2020). Ion channels can be differentially sensitive to changes in temperature or pH, and changes in these variables can have complex effects on ionic currents in neurons (Tang et al., 2010; Tang et al., 2012; Haley et al., 2018). We therefore asked if different environmental perturbations changed the way in which regular rhythms destabilized.
Our analysis mapped a time series of spike times from PD and LP neurons to a categorical time series of labels such as regular. We therefore could measure the transitions between states during different environmental perturbations (Materials and methods). We found that transition matrices between states shared commonalities across environmental perturbations (Figure 6a), such as likely transitions between regular and LP-weak-skipped states. PD-silent-LP-bursting states tended to be followed by PD-silent states, in which the LP neuron is spiking, but not bursting regularly. The LP neuron becomes less regular in both transitions, contributing to the loss of regular rhythms. We never observed a transition from regular rhythms LP-silent or PD-silent states, suggesting slow (>20 s) timescales of rhythm collapse. In high pH, every transition away from the regular state was to the LP-weak-skipped state, hinting at increased sensitivity of the LP neuron to high pH. High pH perturbations also never silenced the circuit, as previously reported (Haley et al., 2018), and showed fewer and less varied transitions than other perturbations. Are some transitions over- or underrepresented in the transition matrix? To determine this, we constructed a null model where transitions occurred with probabilities that scaled with the marginal probability of final states (Materials and methods). Transitions that occurred significantly more often than predicted by the null model are shown with black borders and those that occurred significantly less often than predicted are shown with filled circles (Figure 6a). Transitions that never occurred but occurred at significantly nonzero rates in the null model are indicated with diamonds.
Figure 6
Effect of environmental perturbations on transitions between states.
(a) Transition matrix between states during environmental perturbations. Each matrix shows the conditional probability of observing the final state in the next time step given an observation of the initial state. Probabilities in each row sum to 1. Size of disc scales with probability. Discs with dark borders are transitions that are significantly more likely than the null model (Materials and methods). Dark solid discs are transitions with nonzero probability that are significantly less likely than in the null model. are transitions that are never observed and are significantly less likely than in the null model. States are ordered from regular to silent. (b) Coefficient of variation (CV) of burst period of pyloric dilator (PD) (purple) and lateral pyloric (LP) (red) vs. time before transition away from the regular state. are from Spearman test (on binned data, as plotted) to check if variability increases significantly before transition. Temperature perturbations: transitions in 61 animals; pH perturbations: transitions in 6 animals; perturbations: transitions in 20 animals.
Earlier work has shown that transitions from regular bursting are preceded by an increase in variability in the voltage dynamics of bursting in PD neurons pharmacologically isolated from most of the pyloric circuit (Ratliff et al., 2021). Can we detect similar signatures of destabilization before transitions from regular states in the intact circuit? We measured the CV of the burst periods of PD and LP neurons in regular states just before transitions away from regular (Figure 6b). Because we restricted our measurement of variability to regular states, we could disambiguate true cycle-to-cycle jitter in the timing of bursts from the apparent variability in cycle period due to alternations between bursting and nonbursting dynamics. We found that transitions away from regular were correlated with a steady and almost monotonic increase in variability in PD and LP burst periods for low pH and high perturbations, but not for high pH and high-temperature perturbations (Spearman rank correlation test). This suggests mechanistically different underpinnings to the pathways of destabilization between these sets of perturbations and is consistent with previous work showing that robustness to perturbations in pH only moderately affects temperature robustness in the same neuron (Ratliff et al., 2021).
Decentralization elicits variable circuit dynamics
The pyloric circuit is modulated by a large and chemically diverse family of neuromodulators that it receives via the stomatogastric (stn) nerve (Marder, 2012). Decentralization, or the removal of this neuromodulatory input via transection and/or chemical block of the stn, has been shown to affect the pyloric rhythm in a number of ways (Russell, 1976). Decentralization can stop the rhythm temporarily, which can recover after a few days (Golowasch et al., 1999; Thoby-Brisson and Simmers, 1998). Decentralization slows down the pyloric rhythm (Eisen and Marder, 1982; Rosenbaum and Marder, 2018) and makes the rhythm more variable (Hamood and Marder, 2014; Hamood et al., 2015). Decentralization can evoke variable circuit dynamics, sometimes with slow timescales (Figure 7—figure supplement 1), and can lead to changes in ion channel expression (Mizrahi et al., 2001).
The variability in circuit dynamics elicited by decentralization and the animal-to-animal variability in response to decentralization have made a quantitative analysis of the effects of decentralization difficult. We therefore set about to characterize the variable and invariant features of the changes in circuit spiking dynamics on removal of descending neuromodulation across a large () population.
We first asked where in the map decentralized data were (Figure 7a). A large fraction () of the data was found outside the regular cluster, suggesting the existence of atypical circuit dynamics on decentralization. Decentralization also changed probabilities of observing many states. The regular state was significantly less likely on decentralization, and several atypical states were significantly more likely (Figure 7b and c, Table 2, Figure 7—figure supplement 2).
Figure 7 with 4 supplements see all
Effect of decentralization.
(a) Map occupancy conditional on decentralization. Shading shows all data, and bright colored dots indicate data when preparations are decentralized. (b) State probabilities before and after decentralization. (c) Fold change in state probabilities on decentralization. States marked n.s. are not significantly more or less likely after decentralization. All other states are (paired permutation test, p<0.00016). (a, b) points from animals. (d) Transition matrix during decentralization. Probabilities in each row sum to 1. Size of disc scales with probability. Discs with dark borders are transitions that are significantly more likely than the null model (Materials and methods). Dark solid discs are transitions with nonzero probability that are significantly less likely than in the null model. are transitions that are never observed and are significantly less likely than in the null model. States are ordered from regular to silent. transitions. (e) Coefficient of variation of pyloric dilator (PD, purple) and lateral pyloric (LP; red) burst periods before transition away from regular states. from Spearman test. points from animals.
Table 2
State counts before and after decentralization for the data shown in Figure 7.
p-Values of change in probability of observing change estimated from paired permutation tests.
| State | ||||
|---|---|---|---|---|
| Regular | 7,967 | 5,791 | lt0.001 | -0.308 77 |
| LP-silent | 22 | 724 | lt0.001 | 0.030 65 |
| LP-silent-PD-bursting | 14 | 577 | lt0.001 | 0.045 926 |
| PD-silent | 11 | 140 | 4 | 0.018 51 |
| PD-silent-LP-bursting | 20 | 18 | 0.469 59 | 0.000 188 91 |
| Aberrant-spikes | 111 | 168 | 0.300 37 | 0.003 285 3 |
| LP-weak-skipped | 317 | 1,628 | lt0.001 | 0.099 875 |
| PD-weak-skipped | 142 | 118 | 0.292 19 | 0.003 453 8 |
| Sparse-irregular | 4 | 154 | lt0.001 | 0.013 263 |
| Irregular | 13 | 116 | 0.000 23 | 0.010 877 |
| Silent | 0 | 321 | lt0.001 | 0.024 825 |
| Irregular-bursting | 72 | 753 | lt0.001 | 0.057 913 |
-
LP: lateral pyloric; PD: pyloric dilator.
How do preparations switch between different states when decentralized? The transition matrix during decentralization revealed many transitions between diverse states (Figure 7d), with the most likely transitions being significantly overrepresented compared to the null model (p<0.05, Materials and methods). Transitions away from regular included significantly more likely transitions into states where one of the neurons was irregular such as LP-weak-skipped and PD-weak-skipped. Similar to rhythm destabilization in high or low pH, transitions away from regular were associated with a near-monotonic increase in the variability of PD and LP burst periods before the transitions (Figure 7e, , p<0.006, Spearman rank correlation test).
The time series of identified states on a preparation-by-preparation basis showed striking variability in the responses to decentralization (Figure 7—figure supplement 3a), with the probability of observing regular states decreasing immediately after decentralization (Figure 7—figure supplement 3b). What causes the observed animal-to-animal variability in circuit dynamics on decentralization? One possibility is that seasonal changes in environmental conditions alter the sensitivity of the pyloric circuit to neuromodulation. We tested this hypothesis by measuring the correlation between measures such as the probability of observing the regular state, the change in burst period, and the change in firing rate on decentralization and the sea surface temperature at the approximate location of these wild-caught animals (Figure 7—figure supplement 4). None of these measures was significantly correlated with sea surface temperature (p>0.07, Spearman rank correlation test).
Stereotyped effects of decentralization on burst metrics
Despite the animal-to-animal variation in responses to decentralization, are there stereotyped responses to decentralization? Decentralization removes some unknown mixture of modulators that are released by the stn, which can vary from animal to animal. Previous work has shown that decentralization typically slows down the pyloric rhythm (Eisen and Marder, 1982; Rosenbaum and Marder, 2018) and (Figure 8—figure supplement 1), but a finer-grained analysis of rhythm metrics was confounded by the irregular dynamics that can arise when preparations are decentralized. For example, alteration between regular and atypical states could bias estimates of burst metrics that are not defined in atypical states. Because our analysis allows us to identify the subset of data where pyloric circuit dynamics are regular enough that burst metrics are well-defined, we measured the changes in a number of burst metrics like the burst period, duty cycle, and phases on decentralization (Figure 8a). Every metric measured was significantly changed except the phase at which LP bursts start (p<0.007, paired permutation test). Consistent with earlier studies, we found that the CV in every metric increased following decentralization (Figure 8b).
Figure 8 with 2 supplements see all
Effects of decentralization on burst metrics.
(a) Change in mean burst metrics on decentralization. (b) Change in coefficient of variation of burst metrics on decentralization. In (a) and (b), each dot is a single preparation; * indicate distributions whose mean is significantly different from zero (p<0.007, paired permutation test). Firing rates (c), burst period (d), lateral pyloric (LP) phases (e), and duty cycles (f) vs. time since decentralization. In (c–f), thick lines indicate population means, and shading indicates the standard error of the mean. points from preparations.
What are the dynamics of changes in burst metrics on decentralization? Firing rates of both LP and PD neurons decreased immediately on decentralization, roughly halving their pre-decentralized values (Figure 8c). This occurred together with a doubling of PD burst periods (Figure 8d), suggesting that the entire rhythm is slowing down. Intriguingly, decentralization led to significant advance in the phase of LP burst ends, but not starts (Figure 8e), leading to a large decrease in the duty cycle of the LP neuron (Figure 8f) that was significantly more than the decrease in PD’s duty cycle (p<10-8, paired t-test).
The stereotyped slowing of the rhythm on decentralization can also be quantified by looking at the distribution of the data in the regular cluster before and after decentralization (Figure 8—figure supplement 2). Data are concentrated in the upper-left edge of the regular cluster when decentralized, where burst periods are large and firing rates low (Figure 2—figure supplement 1a and c), suggesting that decentralization could elicit a more stereotyped rhythm for circuits that continue to burst regularly, because circuits that do so tend to share a common, slow bursting dynamics. Counterintuitively, it may appear that regular rhythms in baseline conditions are more variable than regular rhythms after decentralization. To test this hypothesis, we measured the dispersion of each preparation in the map (Figure 8—figure supplement 2b) before and after decentralization. Dynamics before decentralization were significantly more dispersed in the regular cluster than dynamics after decentralization (Figure 8—figure supplement 2c, , paired t-test) because they then tended to be concentrated in the upper-left edge of that cluster. To first approximation, our analysis shows that there are many ways to manifest a regular rhythm under baseline conditions, but regular rhythms on decentralization are typically slow, and stereotyped in comparison.
Neuromodulators differentially affect state probabilities
The crustacean STG is modulated by more than 30 substances (Harris-Warrick and Marder, 1991; Marder, 2012) that tune neuronal properties at an intermediate timescale between feedback homeostasis and intrinsic cellular properties (Daur et al., 2016). Earlier work has focused on understanding the effect modulators have on restoring (or destabilizing) the canonical rhythm, in part because the restoration of regular oscillatory dynamics is a common feature of neuromodulator action. Other effects that neuromodulators might have on pyloric circuit dynamics are harder to investigate and are hindered by the difficulty in characterizing circuit dynamics when nonregular. Here, we set out to systematically characterize the effects of neuromodulators on dynamic states identified in the full space of circuit behaviors (Figure 3).
We focused our analysis on the effect of four neuromodulators: red pigment-concentrating hormone (RPCH), proctolin, oxotremorine, and serotonin. In the experiments analyzed, these neuromodulators were added to decentralized preparations so that endogenous effects of these (and other) neuromodulators were minimized. We therefore first characterized the distribution of states in decentralized preparations where neuromodulators were subsequently added (Figure 9a).
Figure 9 with 2 supplements see all
Effect of bath-applied modulators.
(a) State distribution in decentralized preparations. (b) State distribution in bath application of neuromodulators. Change percentages show difference in probability of regular state from decentralized to addition of neuromodulator. (c) Probability distribution of states conditional on transition to (for red pigment-concentrating hormone [RPCH], proctolin, and oxotremorine) or from (for serotonin) the regular state. is the number of data points, and is the number of animals.
RPCH is a neuropeptide that targets a number of cells in the circuit (Nusbaum and Marder, 1988; Swensen and Marder, 2001) and has been shown to increase the number of spikes per burst in PD and LP (Dickinson et al., 2001; Thirumalai and Marder, 2002), though it has little effect on the pyloric period (Thirumalai et al., 2006). RPCH increased the probability of the regular state, suggesting stabilization of the triphasic rhythm, and decreased the probability of most other atypical states (Figure 9b, Table 3, p<0.004, paired permutation test). Consistent with earlier studies that reported that RPCH can activate rhythms in silent preparations (Nusbaum and Marder, 1988), the probability of observing the silent state was driven to 0 in the presence of RPCH, together with other atypical states such as LP-silent and LP-silent-PD-bursting (Figure 9b).
Table 3
Probability distribution of states during modulator application, as shown in Figure 9.
| State | Decentralized | RPCH | Proctolin | Oxotremorine | Serotonin |
|---|---|---|---|---|---|
| Regular | 0.39 | 0.73 | 0.69 | 0.78 | 0.27 |
| LP-silent | 0.06 | 0 | 0.02 | 0 | 0.07 |
| LP-silent-PD-bursting | 0.09 | 0 | 0.07 | 0 | 0.1 |
| PD-silent | 0.07 | 0 | 0 | 0 | 0.04 |
| PD-silent-LP-bursting | 0.01 | 0 | 0 | 0 | 0.03 |
| Aberrant-spikes | 0.01 | 0.04 | 0.01 | 0.01 | 0.03 |
| LP-weak-skipped | 0.14 | 0.11 | 0.07 | 0.17 | 0.19 |
| PD-weak-skipped | 0.02 | 0.05 | 0 | 0 | 0 |
| Sparse-irregular | 0.03 | 0 | 0.01 | 0 | 0.02 |
| Irregular | 0.02 | 0.02 | 0.01 | 0 | 0.07 |
| Silent | 0.07 | 0 | 0 | 0 | 0.01 |
| Irregular-bursting | 0.1 | 0.04 | 0.11 | 0.03 | 0.17 |
-
LP: lateral pyloric; PD: pyloric dilator.
Proctolin also targets a number of cells in the circuit (Swensen and Marder, 2001) and strengthens the pyloric rhythm through various mechanisms: by increasing the amplitude of slow oscillations in AB and LP (Hooper and Marder, 1987; Nusbaum and Marder, 1989), depolarizing the LP neuron (Golowasch and Marder, 1992; Turrigiano and Marder, 1993), and increasing the number of spikes per burst in LP and PD (Hooper and Marder, 1987; Marder et al., 1986; Hooper and Marder, 1984). Oxotremorine, a muscarinic agonist, has also been shown to enhance the robustness of the pyloric rhythm (Bal et al., 1994; Haddad and Marder, 2018; Rosenbaum and Marder, 2018). Similar to RPCH, both proctolin and oxotremorine significantly increase the probability of the regular state (Figure 9b, Table 3, p<0.004, paired permutation test), and the regular state is the only one significantly more likely when the neuromodulator is added. The strengthening effects of RPCH and oxotremorine are also manifested in the significantly lower probabilities of observing atypical and dysfunctional states such as silent, LP-silent, PD-silent, and sparse-irregular (Table 3).
Serotonin can have variable effects on the pyloric circuit, varying from animal to animal, and can either speed up or slow down the rhythm (Beltz et al., 1984; Spitzer et al., 2008). In Panularis, serotonin depolarizes LP in culture, but hyperpolarizes LP in situ, unlike other neuromodulators that typically have the same effect in situ and in culture (Turrigiano and Marder, 1993). Consistent with earlier work in C. borealis showing that serotonin destabilizes the rhythm in decentralized preparations (Haddad and Marder, 2018), we found that the probability of regular states was significantly lower on addition of serotonin (Figure 9b, Table 3, p<0.004, paired permutation test), together with a significantly higher probability of atypical dysfunctional states such as LP-silent, aberrant-spikes, PD-silent-LP-bursting, and irregular, suggesting loss of coordination between the many neurons in the pyloric circuit with serotonin receptors (Clark et al., 2004).
Do these modulators share common features in how they (de)stabilize the rhythm? We computed the probability distribution of states conditional on transitions to the regular state for RPCH, proctolin, and oxotremorine, and conditional on transitions from the regular state for serotonin (Figure 9c). For all four neuromodulators, the conditional state distribution predominantly comprised these three states: LP-weak-skipped, irregular-bursting, and aberrant-spikes, suggesting that trajectories of recovery or destabilization of the regular rhythm share common features. Serotonin destabilizes the rhythm, decreasing the likelihood of observing regular states, similar to environmental perturbations (Figure 5) and decentralization (Figure 7).
Different neuromodulators activate different forms of the rhythm (Marder and Weimann, 1992; Marder et al., 1985; Marder, 2012), partly because different neuron types express different receptors to varying extents (Garcia et al., 2015). Moreover, similar rhythmic motor patterns can be produced by qualitatively different mechanisms, such as one that depends on voltage-gated sodium channel activity, and one that can persist in their absence (Harris-Warrick and Flamm, 1987; Epstein and Marder, 1990; Rosenbaum and Marder, 2018). To determine if different neuromodulators elicit regular rhythms that occupy different parts of the map, we plotted the location of data elicited by various neuromodulators in the full map (Figure 9—figure supplement 2). Regular data elicited by different neuromodulators tended to lie in clusters, whose distribution in the map was significantly different between serotonin and CCAP (Crustacean cardioactive peptide), and proctolin and every other neuromodulator tested (p<0.05, two-dimensional Kolmogorov–Smirnov test, using the method of Peacock, 1983). The differential clustering of regular states in the map with neuromodulator suggests that neuromodulators can elicit characteristic, distinct rhythms.
Neuromodulators differentially affect transition between states
RPCH, proctolin, and oxotremorine activate a common voltage-dependent modulatory current, (Swensen and Marder, 2001), but can differentially affect neurons in the STG because different cell types express receptors to these modulators to different degrees. For example, RPCH activates strongly in LP neurons, but the effects of oxotremorine and proctolin are more broadly observed in the circuit (Swensen and Marder, 2000; Swensen and Marder, 2001). Though these three modulators strengthen the slow-wave oscillations in pyloric neurons, only oscillations elicited by oxotremorine and RPCH persist in tetrodotoxin, and proctolin rhythms do not, hinting that qualitatively different mechanisms underlie the generation of these seemingly similar oscillations (Rosenbaum and Marder, 2018). We therefore measured the transition rates between states during neuromodulator application to how similar or different trajectories towards recovery were.
In RPCH, proctolin, and oxotremorine application, ≈100 transitions were observed between states (Figure 10). Transitions could not always be predicted by a null model assuming that transition probabilities scaled with the conditional probability of observing states after a transition. For example, some transitions, such as the transition from irregular to regular, were never observed in RPCH, a significant deviation from the expected number of transitions given the likelihood of observing regular states after transitions (Materials and methods). Others, such as the transition LP-silent to LP-silent-PD-bursting in proctolin and oxotremorine, were observed at rates significantly higher than expected from the null model. Transitions into regular state are distributed across aberrant-spikes, LP-weak-skipped, and irregular-bursting states for all three, but no invariant feature emerges in the rest of the transition matrix.
Figure 10
Effect of red pigment-concentrating hormone (RPCH), proctolin, oxotremorine, and serotonin on transition probabilities.
Each matrix shows the conditional probability of observing the final state in the next time step given an observation of the initial state during bath application of that neuromodulator. Probabilities in each row sum to 1. Size of disc scales with probability. Discs with dark borders are transitions that are significantly more likely than the null model (Materials and methods). Dark solid discs are transitions with nonzero probability that are significantly less likely than in the null model. are transitions that are never observed and are significantly less likely than in the null model. States are ordered from regular to silent. Bar graphics show the coefficient of variability (CV) of pyloric dilator (PD) and lateral pyloric (LP) burst periods before transition away from regular states. from Spearman rank correlation test. RPCH: transitions in animals; proctolin: transitions in animals; oxotremorine: transitions in animals; serotonin: transitions in animals. Bar graphs show the CV of burst periods of PD and LP vs. time before a transition away from regular states during serotonin application. from Spearman rank correlation test.
Serotonin destabilizes the rhythm in decentralized preparations, and the transition matrix under serotonin reveals several features of the irregularity behavior observed under serotonin (Figure 10). A number of irregular and low-firing states from silent to irregular never transition into the regular state, which is unlikely in the null model (p<0.05, Materials and methods). Transitions between pairs of states are symmetric and occur at rates significantly larger than in the null model, such as between LP-silent and LP-silent-PD-bursting. Intriguingly, destabilizing transitions from regular to LP-weak-skipped, aberrant-spikes, and irregular-bursting are observed at rates significantly higher than in the null model. These three abnormal states are also observed immediately preceding regular states in RPCH, proctolin, and oxotremorine (Figure 9c), suggesting that the mechanisms for both stabilization and destabilization of the rhythm share stereotyped trajectories.
Are transitions away from regular states also associated with increases in variability of burst periods? Similar to preparations in high and low pH, and when decentralized, transitions away from regular states in serotonin were associated with significantly rising variability in the burst periods of PD and LP neurons (Figure 10, p<0.05, Spearman rank correlation test).
Discussion
Advances in neural recording technology have made it possible to generate increasingly large datasets, and an ongoing challenge is in developing computational tools to find structure in the neural haystack (Pachitariu et al., 2016). Nonlinear dimensionality reduction algorithms such as t-SNE can create a useful representation of datasets that are too large to visualize in their entirety using traditional methods. We combined domain-specific expert knowledge with an unsupervised dimensionality reduction process (t-SNE) by manually segmenting and labeling clusters of dynamics representing biologically significant behavior. This approach conferred two advantages: it allowed for a more accurate measure of traditional metrics such as burst phases in large datasets (Figures 4 and 8), and it allowed for the analysis of irregular dynamics that are typically intractable with conventional analysis methods (e.g., Figure 9), with the disadvantage of not being fully automated, and requiring human intervention to inspect data in the embedding and draw cluster boundaries. Our work hints at a possibility to characterize nonregular spike patterns in small neural circuits and can thus provide a deeper understanding of circuit activity under baseline conditions and in response to perturbations. Our approach makes limited assumptions of the dynamics of the circuit, yet provides a formal framework based on domain-specific knowledge for characterizing circuit activity. Additionally, this way of analyzing neural spike data can readily be adapted to other circuits and systems.
Reliable identification of regular rhythms allows for accurate, interpretable analysis of rhythm metrics
Characterizing the statistics of neural oscillations has several subtle challenges. For example, variations in cycle period arising from cycle-to-cycle fluctuations are not distinguished from those arising from alteration between epochs of regular oscillations interrupted by spans of irregular activity where metrics like cycle period are undefined. One way to disambiguate the two is to construct elaborate checks to make sure that the spike pattern being measured meets certain criteria. However, edge cases abound, and this is a challenging and subjective approach. A fortuitous consequence of the embedding method we used is to reliably identify when rhythms were regular, and we found that burst metrics were well defined for this subset of data (blue region in Figure 3). We were therefore able to measure the mean and variability of various burst metrics (Figure 4) only in stretches of data where it made sense to do so, and thus the measured variability stemmed almost entirely from cycle-to-cycle variations.
Consistent with previous studies (Bucher et al., 2005; Hamood and Marder, 2014; Hamood et al., 2015), our results (Figure 4) show that within-animal variability in pyloric burst metrics is lower than across-animal variability. Our results are from an analysis of data from several experimenters from different laboratories, collected over a span of 10 years. It is therefore an ideal dataset in which to measure variability. We find that the CV of all burst metrics measured is (Figure 4b), which can now be used as a standard for regular baseline pyloric oscillations. Measuring burst metrics on decentralization (Figure 8) also allowed us to characterize how regular rhythms change, while still being recognizably regular. In addition to recapitulating well-understood phenomena such as the slowing down and increased variability in rhythms, we found that the phase of LP burst onset did not change significantly, but the phase of LP burst termination did, suggesting that features of the rhythm are differentially robust to the removal of neuromodulation.
Earlier work categorized the varied dynamics of the pyloric circuit during perturbations (Haddad and Marder, 2018; Haley et al., 2018; Ratliff et al., 2021). In those studies, categories were typically constructed by hand and were not rigorously shown to be mutually exclusive. Categories in this work, while manually chosen, emerge naturally from the distribution of the data in the reduced space (Figure 3) and no segment of data can have more than one label because it can exist only at a single point in the map. For instance, earlier work categorized rhythms that were labeled regular into two categories, ‘normal triphasic’ and ‘normal triphasic slow’ (Haddad and Marder, 2018), while we did not observe a distinctly bimodal distribution of burst periods. In contrast, the catch-all ‘atypical firing’ state was separated here into a number of states (irregular, irregular-bursting, sparse-irregular) that span several well-separated clusters in the map (Figure 3). In summary, this work recapitulates every label constructed to categorize spike patterns from PD and LP neurons in earlier work, and additionally finds new spike patterns that were either not detected or not identified as distinct because they are hard to detect by manual inspection.
Diversity and stereotypy in trajectories from functional to crash states
Are there preferred paths to go from regular rhythms to crash? Diversity in the solution space of functional circuits, and the varied effects of perturbations on these circuits, argues for an assortment of trajectories from functional dynamics to irregular or silent states. While transition matrices measured during different perturbations were varied (Figure 6), we did observe universal features in transition matrices measured during environmental perturbations, decentralization, and addition of neuromodulators (Figures 6, 7 and 10). The destabilizing transition from regular → LP-weak-skipped was overrepresented in every transition matrix, suggesting that the weakening of the LP neuron is a crucial step in the trajectories towards destabilization, perhaps because there is only one copy of LP in the circuit. Earlier work studying trajectories of destabilization of regular bursting in the isolated pacemaker kernel also found a conserved motif in trajectories towards destabilization: from regular bursting to tonic spiking to silence in response to pH perturbations, and another conserved motif (bursting to weak bursting to silence) in response to temperature perturbations (Ratliff et al., 2021). Transitions away from regular rhythms were also associated with increased variability in burst periods during all perturbations except high temperature and low pH (Figures 6, 7 and 10). An increase in variability in PD voltage dynamics before transitions from regular bursting has been observed in the isolated pacemaker kernel (Ratliff et al., 2021), similar to the effect we observed in the intact circuit.
The structure of the transitions between states also hints at features of the circuit that are critical for rhythm (de)stabilization. Unsurprisingly, PD-silent states precede silent states in low pH, high temperature, and high perturbations (Figure 6). This makes sense because PD cells are electrically coupled to the endogenous burster AB in the pacemaker kernel, and silencing the pacemaker kernel can cause the circuit to go silent. Though the states are determined purely from clusters in the embedding (Figure 2), and thus from statistical features of spike times, some states may be identified predominantly with cell-specific features (e.g., LP-weak-skipped where the LP neuron fails to burst regularly, but the PD neurons continue to burst regularly), or with circuit-level features (e.g., aberrant-spikes where one or both neurons fire spikes outside the main burst, which may be caused by incomplete inhibition). Decentralization elicits the largest number of transition types, with of all transition types observed, which could be a consequence of the complex change in the neuromodulator milieu following transection of descending nerves.
Linking circuit output to circuit mechanisms
A large body of work has shown that there is more than one way to make a neural circuit with similar patterns of activity (Prinz et al., 2003; Prinz et al., 2004; Gutierrez et al., 2013). Several combinations of circuit parameters such as synapse strengths, ion channel conductances, and network topology can be found in circuits that generate similar emergent collective dynamics (Gonçalves et al., 2020). The dimensionality of the space of neuronal and synaptic parameters in a neural circuit is much larger than the dimensionality of the circuit output (Marder and Bucher, 2007). This disparity in dimensionality leads to an inherently many-to-one mapping from the space of circuit architecture to the space of circuit dynamics. Circuits can therefore exhibit ‘cryptic’ architectural variability (Haddad and Marder, 2018), where the diversity of topologies and neuronal parameters is masked by the relatively low-dimensional nature of the observed circuit outputs. However, perturbations can reveal differences between seemingly identical circuits. For instance, current injections in an oscillator network can shift phases, thus revealing connection weights between individual neurons (Timme, 2007). This work reveals a path towards analysis that can reveal cryptic variability and build mechanistic links from circuit architecture to function. By characterizing the totality of circuit dynamics under a variety of conditions, our framework provides a way to fit biophysically detailed models of the pyloric circuit to diverse circuit dynamics under baseline conditions and perturbations. From the large diversity of neuron and circuit parameters that can reproduce a snapshot of activity, only a subset of models could potentially recapitulate the diverse irregular behavior seen under extreme perturbations. Recent work that reproduced how circuits change cycle periods with temperature (Alonso and Marder, 2020) can be extended to find parameter sets that also generate the irregular states characterized in this study at the rates observed in the data. Crucially, the characterization of the pyloric circuit dynamics in this work can be used to rule out models and parameter sets that generate irregular activity that is qualitatively dissimilar to any of the irregular states observed in the pyloric circuit. Future experimental work can pair data analysis methods such as this work with quantitative measurements of cellular and circuit parameters using emerging techniques (Schulz et al., 2006; Schulz et al., 2007; Tobin et al., 2009) to find parameter sets that generate robust rhythms and irregular states.
Applicability to other systems
The analysis method in this study is well-suited for large datasets of neural recordings from identified neurons. Data where the identity of each neuron is not or cannot be known, such as large-scale mammalian brain recordings, would require modifications to the analysis pipeline described in Figure 2. First, it would no longer be possible to construct a data vector of fixed length because ordering of the different neurons would not be meaningful. Each data point would instead be an unordered set of spike times from each neuron, and a distance function that operated on spike times (Christen et al., 2006; Victor and Purpura, 2009; Schreiber et al., 2003; van Rossum, 2001) could be used to generate a distance matrix between raw data points, which would be the input to the embedding algorithm. In our analysis, we included features such as the ‘spike phase’ (Figure 2b and c) because the neurons in this circuit interact with one another strongly in each cycle of oscillations. The analysis of neural circuits that do not show such strong intrinsically phase-controlled behavior could use other features more suitable to those systems.
Comparison with other methods
Visualization and other forms of analysis of large neural datasets rely on dimensionality reduction (Nguyen and Holmes, 2019). Here, we used the t-SNE algorithm as a core method to reduce the dimensionality of the dataset and visualize our data. t-SNE has been widely used in the unsupervised analysis of many types of biological data (Berman et al., 2014; Kollmorgen et al., 2020; Chen et al., 2021; Macosko et al., 2015; Kobak and Berens, 2019; Leelatian et al., 2020), including neural recordings (Dimitriadis et al., 2018). t-SNE is a technique that allows high-dimensional data to be visualized in a lower-dimensional space (Van der Maaten and Hinton, 2008; Linderman and Steinerberger, 2019), and works by preserving pairwise distances between points in the high-dimensional space and the low-dimensional embedding, within a certain neighborhood. This feature makes t-SNE an attractive tool to try to visualize large, structured datasets, such as those examined in this study, because it can demonstrate how similar spike patterns are to each other (Dimitriadis et al., 2018). t-SNE has been shown rigorously to be capable of recovering well-separated data clusters (Linderman and Steinerberger, 2019). In our application, t-SNE generated embeddings where spike patterns in different regions could be described as qualitatively different. For example, spike patterns in the top-most cluster (colored green in Figure 3) all had weak PD spiking, but regular and strong LP spiking. This was qualitatively different from the two closest clusters LP-weak-skipped and irregular. In regions of the map where clusters were not cleanly separated (e.g., in the connection between the regular and irregular-bursting clusters), manual inspection revealed a number of intermediate states. The clustered or not-clustered regions of the map are therefore informative of the underlying distribution of spike patterns and emerge robustly from the embedding.
t-SNE is widely used in the analysis and visualization of high-dimensional data, but is important to acknowledge its limitations. t-SNE can generate embeddings that appear to have clusters from purely randomly distributed data, can distort sizes of clusters, and can fail to preserve large-scale topological features of the data in some embeddings (Wattenberg et al., 2016). The visualization we generated was useful in that it guided manual clustering and made feasible a previously intractable task, that of classifying hundreds of hours of spike patterns from hundreds of animals.
A variety of other dimensional reduction techniques, including multidimensional scaling (Cox and Cox, 2008), convolutional non-negative matrix factorization (Mackevicius et al., 2019) and their extensions (Williams et al., 2020), tensor component analysis (Williams et al., 2018), and dynamical component analysis (Clark et al., 2019), have been developed that aid in visualizing and analysis of large neural datasets. Methods based on neural networks offer powerful tools to analyze unstructured neural data by modeling the data with a recurrent neural net and then analyzing that model (Vyas et al., 2020). Topological autoencoders are one such technique that combine autoencoders with methods from topological data analysis to produce representation in lower-dimensional spaces (Moor et al., 2019). These methods are similar in spirit to the analysis presented here, but use sophisticated neural nets whose parameters yield the lower-dimensional representation. Other analysis methods include SOM-VAE, which combines self-organizing maps (SOMs) and variational auto-encoders (VAEs) (Fortuin et al., 2018) to analyze high-dimensional time series and find transitions between states, and deep temporal clustering, which combines dimensionality reduction and temporal clustering (Madiraju et al., 2018).
Technical considerations
In this study, we have used the activity of the LP and PD neurons as a proxy for the pyloric circuit. However, the pyloric circuit contains other neurons: AB (which is electrically coupled to PD neurons), PY neurons (which are anti-phase to both PD and LP), and VD and IC neurons. A richer description of the dynamics of the pyloric circuit would include spikes from these neurons, and the methods we have described here can be scaled up to include these neurons. It is likely that we are underestimating the number of states, and thus, transitions between states, because we do not have access to the dynamics of these neurons. Datasets that contain recordings from all pyloric neurons as preparations are subjected to the perturbations studied here and are not available for large numbers of animals. We therefore chose to focus on the functional antagonists LP and PD. Additionally, neurons in the pyloric circuit are coupled using graded synapses, and the circuit can generate coordinated activity even when spiking is abolished (Rosenbaum and Marder, 2018), suggesting that subthreshold oscillations may be an important feature we are not measuring by only recording spikes. However, the data required necessitates substantially harder to perform experiments because intracellular electrodes must be used. Furthermore, the signal to the muscles – arguably the physiologically and functionally relevant signal – is the spike signal, suggesting that spike patterns from the pyloric circuit are a useful feature to measure.
The unit of data we operated on was a time series of spikes from the LP and PD neurons. In order to describe what the dynamics of these neurons is at a given point in time, we chose to look at a neighborhood in time. In this article, we chose 20 s nonoverlapping bins, based on inspection of the data by eye. Choosing a time bin imposes certain tradeoffs in the analysis of time series: changes in dynamics on timescales smaller than the bin are counted as different states, and changes in dynamics on timescales longer than the bin size are counted as transitions between states. The statistics of the transitions we measure are therefore dependent on the bin size we chose. We note that dwell times in each state are almost always in excess of null model predictions generated by shuffling states (transition matrices in Figures 6, 7 and 10), supporting the validity of our choice of 20 s bins.
Conclusion and outlook
Our work provides a way to characterize nonregular spike patterns in small neural circuits. It thus provides a bridge between experimental or simulation work grounded in the biophysical detail of ion channels and synaptic currents; and the rich body of observations of circuits under baseline and perturbed conditions. The methods we have employed can easily be adapted to other circuits and systems, make limited assumptions of the dynamics of the circuit, yet provide a robust framework on which to hang a large volume of previously ineffable expert domain knowledge.
Prior to this work, crashes in the pyloric circuit and irregular dynamics in a normally regular circuit were difficult to characterize. We present a method to tame the complexity of the distribution of irregular states by exploiting the fact that pyloric dynamics are not unbounded even in their irregularity. By using a t-SNE in conjunction with manual inspection of reduced data and manual clustering, we have made this previously intractable problem feasible and found undiscovered spike patterns and transitions. Our approach was successful because we used a dataset with long recordings from identified neurons in a circuit that can be subjected to many different perturbations, which is one of the advantages of using the STG system. It will be interesting to see if this can be applied to other systems with identified neurons in a functional circuit to characterize their function and failure modes.
In intact pyloric circuits, and in the presence of modulatory input from the stn, almost all networks are ‘normal’ and exhibit regular rhythms. Decentralization can generate more variable dynamics, presumably because the underlying differences in network structure that were compensated by modulator action now manifest as different collective dynamics in the network. Although it may appear that modulators can have similar effects when added to a decentralized network, they are in fact distinguishable when looking at how they influence the totality of circuit dynamics, not just the regular state.
A major unanswered question was whether crashes triggered by different perturbations share dynamical mechanisms and common pathways. Earlier work looking at a simpler subset of the pyloric circuit argued that different perturbations led to stereotyped but diverse transitions before crash, and we have extended this result in the intact circuit. We show that different perturbations can have different trajectories to crash, but the stereotypy observed in the simpler system was not observed, presumably due to the larger number of pathways accessible to the intact circuit. The new insight from this work stems from the fact that this is the first time transitions through multiple physiological conditions in so many modalities have been characterized and shows that there are many paths through possible circuit dynamical states from canonical states to crash. Several studies focus on one perturbation at a time. By studying a number of perturbations together, we compare responses to different kinds of perturbations on the same physiological network.
Materials and methods
Animals and experimental methods
Request a detailed protocolAdult male Jonah crabs (C. borealis) were obtained from Commercial Lobster (Boston, MA), Seabra’s Market (Newark, NJ), and Garden Farm Market (Newark, NJ). Dissections were carried out as previously described (Gutierrez and Grashow, 2009). Decentralization was carried out either by cutting the stn or by additionally constructing a well on the stn and adding sucrose and TTX (tetrodotoxin) as described in Haddad and Marder, 2018. Temperature was controlled as described in Tang et al., 2010; Tang et al., 2012; Haddad and Marder, 2018. Extracellular potassium concentrations were varied as described in He et al., 2020. pH perturbations are described in Haley et al., 2018.
Data selection and curation
Request a detailed protocolOur goal was to include as much data as possible to create as complete a description of pyloric dynamics as possible. Following our strategy of including only the LP and PD neurons, we used every available dataset that recorded from these neurons from the Marder lab. We also included available datasets from the Nadim and Bucher labs. No dataset was explicitly excluded for reasons linked to the activity of the pyloric circuit in those datasets. Data where crucial metadata was not recorded (e.g., if the temperature of the preparation was not recorded) was excluded. Data where only lvn was recorded from was only included in cases of exceptional data quality, where it was judged that PD and LP could be reliably identified.
Spike identification and sorting
Request a detailed protocolSpikes are identified from extracellular recordings of motor nerves or from intracellular recordings. LP spikes were identified from intracellular recordings, lvn, lpn, and gpn nerves (in descending order of likelihood). PD spikes were identified from pdn, intracellular recordings, and lvn. We used a custom-designed spike identification and sorting software (called ‘crabsort’) that we have made freely available at https://github.com/sg-s/crabsort (copy archived at swh:1:rev:6a67e765e90caa536e6a11f67d9d4737d059af50; Gorur-Shandilya, 2021), previously described in Powell et al., 2021. Spikes are identified using a fully connected neural network that learns spike shapes from small labeled datasets. A new network is typically initialized for every preparation. Predictions from the neural network also indicate the confidence of the network in these predictions, and uncertain predictions are inspected and labeled and the neural network learns from these using an active learning framework (Settles, 2009).
Data curation and data model
Request a detailed protocolEach file was split into 20 s nonoverlapping bins, and spike times, together with metadata, were assembled into a single immutable instance of a custom-built class (embedding.DataStore). The data store had the following attributes:
Spike times containing LP and PD spike times.
ISIs containing ISIs and spike phases
Labels categorical data containing manually generated labels from Figure 3
Metadata such as concentration of modulators, pH, temperature, whether the preparation was decentralized or not, etc.
Using an immutable data structure, reduced risks of accidental data alteration during analysis. Every attribute was defined for every data point.
Embedding
ISI and phase representation (Figure 2b)
Request a detailed protocolEach data point is a 20 s bin containing spike times from LP and PD neurons (Figure 2a). For each data point, spike times are converted into ISIs. A set of spike times uniquely identifies a set of (ordered) ISIs. The set of LP spike times generates a set of LP ISIs, and the set of PD spike times generates a set of PD ISIs (Figure 2b).
For every spike in PD or LP, a ‘spike phase’ can be calculated as follows. Spike phases are not defined when either LP or PD are silent in that data point, or for LP/PD spikes with no spikes from the other neuron before or after that spike. Thus, the ‘spike phase’ of the ith spike on neuron X w.r.t. neuron Y is given by
where is the time of the ith spike on neuron , is the time of the last spike on before , and is the time of the first spike after . Note that this definition can be generalized to neurons, though the number of spike phases grows combinatorially with .
Construction of vectorized data frame (Figure 2c and d)
Request a detailed protocolEach data point can contain an arbitrary number of spikes, and thus an arbitrary number of ISIs and spike phases. Ideally, each data point is a data frame of fixed length (a point in some fixed high-dimensional space). To do so, we computed percentiles ISIs and spike phases (Figure 2c). We chose 10 bins per ISI type (deciles). The end result is not strongly dependent on the number of bins chosen as long as there are sufficiently many bins to capture the distinctly bimodal distribution in ISIs during bursting.
We included four other features to help separate spike patterns that appeared qualitatively different. First, firing rates of LP and PD neurons. Second, the ratios of second-order to first-order ISIs, defined as
where is the nth order set of ISIs computed as the time between n spikes. is the simple set of ISIs defined between subsequent spikes. This measure is included because it captures the difference between single spike bursts and normal bursts well. Third, the ratio between the largest and second-largest ISIs for each neuron.
Finally, we also included a metric defined as follows:
where is a vector of sorted ISIs, and is the sorted ISI for which the difference between it and the previous sorted ISI is maximum. This metric was included as it captures to a first approximation how ‘burst-like’ a spike train is. Intuitively, this metric is high for spike trains with bimodal ISI distributions, as is the case during bursts.
All these features were combined into a single data frame and -scored (Figure 2d).
In some cases, these features were not defined, for example, when there are no spikes on either neuron, the concepts of spike phases or ISIs are meaningless. In these cases, ‘filler’ values were used that were located well off the extremes of the distribution of the metric when defined. For example, ISIs were filled with values of 20 s (the size of the bin) when no spikes were observed. The overall results and shape of the embedding did not depend sensitively on the value of the filler values used.
These features were chosen to capture various modes of spiking and bursting that have been previously identified by manual inspection (Haddad and Marder, 2018; Tang et al., 2012; Haley et al., 2018). Other features may be more appropriate in other systems where spike patterns span different axes of variability. However, we note that these features while being appropriate for this data were not ‘fine-tuned’ to specialize in features that are exclusively found in spike patterns from the pyloric circuit. For example, these features do not explicitly measure bursting, the dominant feature of the pyloric rhythm, but instead use distributions of ISIs that are sufficiently descriptive to capture the variability in bursting and transitions from bursting to other spiking.
Embedding using t-SNE
Request a detailed protocolSo far, we have described how we converted a 20 s snippet containing spike times from LP and PD into a data frame (a vector). We did this for every 20 s snippet in the dataset. Data that did not fit into any bin was discarded (e.g., data at the trailing end of an experiment shorter than 20 s). Thus, our entire dataset is represented by matrix, where is the number of features in the data frame and is the number of data points.
We used the t-SNE algorithm (Van der Maaten and Hinton, 2008) to visualize the vectorized data matrix in two dimensions. Our dataset contained points and was therefore too large for easy use of the original t-SNE algorithm. We used the FI-t-SNE approximate algorithm (Linderman et al., 2019) to generate these embeddings. We used a perplexity of to generate these embeddings. Varying perplexity caused the embedding to change in ways consistent with what is expected for t-SNE embeddings, and the coarse features of the embedding did not sensitively depend on this choice of perplexity (Figure 2—figure supplement 4). t-SNE is often used with random initialization, and different random initializations can lead to different embeddings with clusters located at different positions in the map. The importance of meaningful initializations has recently been highlighted (Kobak and Linderman, 2021), and we used a fixed initialization where the x-axis corresponded to the shortest ISI in each data point and the y-axis corresponded to the maximum ratio of second-order to first-order ISI ratios (described above). For completeness, we also generated embeddings using other initializations (Figure 3—figure supplement 2). For both random initializations (Figure 3—figure supplement 2a–d) and initializations based on ISIs (Figure 3—figure supplement 2e and f), we observed that regular states tended to occur in a single region, surrounded by clusters that were dominated by a single color corresponding to irregular states. Thus, the precise location of different clusters can vary with the initialization, but the overall structure of the embedding, and the identity of points that tend to co-occur in a cluster, does not vary substantially with initialization.
Manual clustering and annotation of data
Once the feature vectors were embedded using t-SNE, we manually inspected these points to get a sense of the spike patterns in each point cloud. To do so, we built an interactive tool that visualized spike patterns that corresponded to each point when clicked on. Random points within regions of high density were sampled to check that interior points had similar spike patterns. Points were assigned labels by drawing boundaries around them and labeling all points within that boundary. Finally, we generated plots of ISIs and rasters from points in clusters to ensure that patterns of spiking were visually similar.
Triangulation and triadic differences (Figure 2—figure supplement 3)
Request a detailed protocolThe output of the embedding algorithm is a set of points in two dimensions. We built a Delaunay triangulation on these points. For each triangle in the triangulation, we computed the maximum difference between some burst metric (e.g., burst period of PD neurons) across the three vertices of that triangle. These triadic differences are represented colored dots, where the dots are located at the incenters of each triangle in the triangulation.
Time-series analysis
Measuring burst metrics (Figure 4)
Request a detailed protocolBurst metrics were measured following previous definitions (Prinz et al., 2004; Bucher et al., 2005). Briefly, bursts were identified by observing that ISI distributions were bimodal, with smaller ISIs corresponding to ISIs within a burst, and longer ISIs corresponding to inter-burst intervals. This allowed us to threshold ISIs, and this identifies burst starts and burst ends. From here, burst periods could be calculated, which allowed us to measure phases and delays relative to the start of the PD burst.
Measuring transition matrices (Figures 6, 7 and 10)
Request a detailed protocolThe transition matrix is a square matrix of size that describes the probability of transitioning from one to another of possible states. The transition matrix we report is the right stochastic matrix, where rows sum to 1. Each element of the matrix corresponds to the conditional probability that we observe state given state . To compute this, we iterate over the sequence of states and compare the current state to the state in the next state. Breakpoints in the sequence are identified by discontinuities in the timestamps of that sequence and are ignored. We then zeroed the diagonal of the matrix and normalized each row by the sum.
Measuring variability before transitions away from regular states (Figures 6 and 7)
Request a detailed protocolWe first identified continuous segments that corresponded to uninterrupted recordings from the same preparation at the appropriate condition. For each segment, we found all transitions away from the regular state. We therefore computed a vector as long as the segment containing the time to the next transition. We then collected points corresponding to time to next transition ranging from to . For each time bin, we measured the CV of the burst period by dividing the standard deviation of the burst period in that datum by the mean in that datum.
Data visualization
Raincloud plots (Figure 4)
Request a detailed protocolRaincloud plots (Allen et al., 2019) are used to visualize a univariate distribution. Individual points are plotted as dots, and a shaded region indicates the overall shape of the distribution. This shape is obtained by estimating a kernel smoothing function estimate over the data. Individual points are randomly jittered along the vertical axis for visibility.
Occupancy maps (Figures 5 and 7)
Request a detailed protocolTo visualize where in the map data from a certain condition occurred, the full embedding is first plotted with colors corresponding to the state each point belongs to. The full dataset is made semi-transparent and plotted with larger dots to emphasize the data of interest. Data in the condition of interest is then plotted as usual. Each bright point in these plots corresponds to a 20 s snippet of data in the condition indicated.
Treemaps (Figures 7 and 9)
Request a detailed protocolTreemaps (Shneiderman and Wattenberg, 2001) were used to visualize state probabilities in a given experimental condition. For each preparation, the probability of each state was computed, and the mean probability of a given state was computed by averaging across all preparations. Thus, each preparation contributes equally. The area of the region in the treemap scales with the probability of that state.
Transition matrices (Figures 6, 7 and 10)
Request a detailed protocolTransition matrices were visualized as in Corver et al., 2021. Initial states are shown along the left edge, and final states are shown along the bottom edge of each matrix. Lines are colored by origin (horizontal lines) or destination (vertical) states. The size of each disc at the intersection of each line scales with the conditional probability of moving from the initial state to the final state. Note that the size of all discs is offset by a constant to make small discs visible.
Statistics
Comparing within-group to across-group variability (Figure 4)
Request a detailed protocolTo compare the variability of various burst metrics within each animal and across animals, we first measured the means and CVs of each burst metrics in every animal. We then used the mean of the coefficients of variations as a proxy for the within-animal variability and used the CV of the means as a proxy for the across-animal variability. Note that both measures are dimensionless. They can therefore be directly compared.
To test if the within-animal variability was significantly less than the across-animal variability, we performed a permutation test. We shuffled the labels identifying the animal to which each data point belonged to and measured a new ‘within-animal’ and ‘across-animal’ variability measure using these shuffled labels. We repeated this process 1000 times to obtain a null distribution of differences between within- and across-animal variability. Identifying where in the null distribution the data occurred allowed us to estimate a p-value for the measured difference. For example, if the measured difference between within- and across-animal variability in metric X was greater than 99% of the null distribution obtained by shuffling labels, we conclude that the p-value is 0.01. The significance level of 0.05 was divided by the number of burst metrics we tested to determine if any one metric was significantly more or less variable across animals.
Measuring trends in variability in regular rhythms before transitions (Figures 6b, 7f and 9d)
Request a detailed protocolTo determine if variability significantly increased in the 200 s preceding a transition away from regular, we measured the Spearman rank correlation between time before transition (x-axis) and mean variability. The Spearman rank correlation is 1 if quantities monotonically increase.
Measuring transition rate significance (Figures 6a, 7e and 10)
Request a detailed protocolIn the empirical transition matrices, certain transitions never occur, and certain transitions occur with relatively high probability. Each element of the transition matrix corresponds to the conditional probability . Our null model assumes that transitions occur at random between states, and therefore the probability of observing any transition scales with the marginal probability of observing state after transitions. We therefore built a null distribution of transition rates by sampling with replacement from the marginal counts of states after transitions. The fraction of this null distribution that was above or below the empirical transition rate was interpreted to be the p-value and thus determined significance.
Code availability
Request a detailed protocolTable 4 lists the code used in this article. The code can be downloaded by prefixing https://github.com/ to the project name.
Table 4
Code availability.
| Project | Notes |
|---|---|
| sg-s/crabsort | Interactive toolbox to sort spikes from extracellular data |
| sg-s/stg-embedding | Contains all scripts used to generate every figure in this article |
| KlugerLab/FIt-SNE | Fast interpolation-based t-distributed stochastic neighbor embedding, used to make embedding |
| sg-s/SeaSurfaceTemperature | Wrapper to scrape NOAA databases |
Data availability
All data needed to reproduce figures in this paper are available at https://zenodo.org/record/5090130.
-
ZenodoMapping circuit dynamics during function and dysfunction.https://doi.org/10.5281/zenodo.5090130
References
-
Muscarinic modulation of a pattern-generating network: control of neuronal propertiesThe Journal of Neuroscience 14:3019–3035.
-
Mapping the stereotyped behaviour of freely moving fruit fliesJournal of the Royal Society, Interface 11:20140672.https://doi.org/10.1098/rsif.2014.0672
-
Visualizing knowledge domainsAnnual Review of Information Science and Technology 37:179–255.https://doi.org/10.1002/aris.1440370106
-
Animal-to-animal variability in motor pattern production in adults and during growthThe Journal of Neuroscience 25:1611–1619.https://doi.org/10.1523/JNEUROSCI.3679-04.2005
-
Many Parameter Sets in a Multicompartment Model Oscillator Are Robust to Temperature PerturbationsThe Journal of Neuroscience 34:4963–4975.https://doi.org/10.1523/JNEUROSCI.0280-14.2014
-
Measuring spike pattern reliability with the lempel–ziv-distanceJournal of Neuroscience Methods 156:342–350.https://doi.org/10.1016/j.jneumeth.2006.02.023
-
A modulatory role for oxygen in shaping rhythmic motor output patterns of neuronal networksRespiration Physiology 128:299–315.https://doi.org/10.1016/s0034-5687(01)00308-5
-
Distinct movement patterns generate stages of spider web-buildingCurrent Biology 31:4983–4997.https://doi.org/10.1016/j.cub.2021.09.030
-
BookMultidimensional scalingIn: Cox MA, editors. Handbook of Data Visualization. Springer. pp. 315–347.https://doi.org/10.1007/978-3-540-33037-0_14
-
Dimensionality reduction for large-scale neural recordingsNature Neuroscience 17:1500–1509.https://doi.org/10.1038/nn.3776
-
Bursting in leech heart interneurons: cell-autonomous and network-based mechanismsThe Journal of Neuroscience 22:10580–10592.
-
The complexity of small circuits: the stomatogastric nervous systemCurrent Opinion in Neurobiology 41:1–7.https://doi.org/10.1016/j.conb.2016.07.005
-
Rpch modulation of a multi-oscillator network: effects on the pyloric network of the spiny lobsterJournal of Neurophysiology 85:1424–1435.https://doi.org/10.1152/jn.2001.85.4.1424
-
t-SNE Visualization of Large-Scale Neural RecordingsNeural Computation 30:1750–1774.https://doi.org/10.1162/neco_a_01097
-
A mechanism for production of phase shifts in a pattern generatorJournal of Neurophysiology 51:1375–1393.https://doi.org/10.1152/jn.1984.51.6.1375
-
Multiple modes of a conditional neural oscillatorBiological Cybernetics 63:25–34.https://doi.org/10.1007/BF00202450
-
Positive Dynamical Networks in Neuronal Regulation: How Tunable Variability Coexists With RobustnessIEEE Control Systems Letters 4:946–951.https://doi.org/10.1109/LCSYS.2020.2997214
-
Proctolin activates an inward current whose voltage dependence is modified by extracellularThe Journal of Neuroscience 12:810–817.https://doi.org/10.1523/JNEUROSCI.12-03-00810.1992
-
Network Stability from Activity-Dependent Regulation of Neuronal ConductancesNeural Computation 11:1079–1096.https://doi.org/10.1162/089976699300016359
-
Failure of Averaging in the Construction of a Conductance-Based Neuron ModelJournal of Neurophysiology 87:1129–1131.https://doi.org/10.1152/jn.00412.2001
-
Cancer borealis stomatogastric nervous system dissectionJournal of Visualized Experiments: JoVE 25:1207.https://doi.org/10.3791/1207
-
ConferenceAnimal-to-Animal Variability in Neuromodulation and Circuit FunctionCold Spring Harbor Symposia on Quantitative Biology. pp. 21–28.https://doi.org/10.1101/sqb.2014.79.024828
-
Multiple mechanisms of bursting in a conditional bursting neuronThe Journal of Neuroscience 7:2113–2128.
-
Modulation of neural networks for behaviorAnnual Review of Neuroscience 14:39–57.https://doi.org/10.1146/annurev.ne.14.030191.000351
-
Motor patterns in the stomatogastric ganglion of the lobster Panulirus argusThe Journal of Experimental Biology 62:405–420.https://doi.org/10.1242/jeb.62.2.405
-
Rapid adaptation to elevated extracellular potassium in the pyloric circuit of the crab, Cancer borealisJournal of Neurophysiology 123:2075–2089.https://doi.org/10.1152/jn.00135.2020
-
Modulation of the lobster pyloric rhythm by the peptide proctolinThe Journal of Neuroscience 7:2097–2112.https://doi.org/10.1523/JNEUROSCI.07-07-02097.1987
-
Cell dialysis by sharp electrodes can cause nonphysiological changes in neuron propertiesJournal of Neurophysiology 114:1255–1271.https://doi.org/10.1152/jn.01010.2014
-
The art of using t-SNE for single-cell transcriptomicsNature Communications 10:5416.https://doi.org/10.1038/s41467-019-13056-x
-
In vivo effects of temperature on the heart and pyloric rhythms in the crab Cancer borealisJournal of Experimental Biology 222:jeb.199190.https://doi.org/10.1242/jeb.199190
-
Clustering with t-SNE, provablySIAM Journal on Mathematics of Data Science 1:313–332.https://doi.org/10.1137/18m1216134
-
Neurotransmitter modulation of the stomatogastric ganglion of decapod crustaceansModel Neural Networks and Behavior 7:319–337.https://doi.org/10.1007/978-1-4757-5858-0_17
-
Modulatory action and distribution of the neuropeptide proctolin in the crustacean stomatogastric nervous systemThe Journal of Comparative Neurology 243:454–467.https://doi.org/10.1002/cne.902430403
-
Modulatory control of multiple task processing in the stomatogastric nervous systemNeurobiology of Motor Programme Selection 6:3–19.https://doi.org/10.1016/B978-0-08-041986-2.50006-0
-
Understanding circuit dynamics using the stomatogastric nervous system of lobsters and crabsAnnual Review of Physiology 69:291–316.https://doi.org/10.1146/annurev.physiol.69.031905.161516
-
Umap: Uniform manifold approximation and projectionJournal of Open Source Software 3:861.https://doi.org/10.21105/joss.00861
-
Designing the dynamics of spiking neural networksPhysical Review Letters 97:18.https://doi.org/10.1103/PhysRevLett.97.188101
-
Ten quick tips for effective dimensionality reductionPLOS Computational Biology 15:e1006907-19.https://doi.org/10.1371/journal.pcbi.1006907
-
A modulatory proctolin-containing neuron (MPN). I. Identification and characterizationThe Journal of Neuroscience 9:1591–1599.https://doi.org/10.1523/JNEUROSCI.09-05-01591.1989
-
Dimensionality reduction in neuroscienceCurrent Biology 26:R656–R660.https://doi.org/10.1016/j.cub.2016.05.029
-
Two-dimensional goodness-of-fit testing in astronomyMonthly Notices of the Royal Astronomical Society 202:615–627.https://doi.org/10.1093/mnras/202.3.615
-
Alternative to Hand-Tuning Conductance-Based Models: Construction and Analysis of Databases of Model NeuronsJournal of Neurophysiology 90:3998–4015.https://doi.org/10.1152/jn.00641.2003
-
Similar network activity from disparate circuit parametersNature Neuroscience 7:1345–1352.https://doi.org/10.1038/nn1352
-
Temperature and acid–base balance in the american lobster Homarus americanusThe Journal of Experimental Biology 210:1245–1254.https://doi.org/10.1242/jeb.02709
-
Neuronal oscillator robustness to multiple global perturbationsBiophysical Journal 120:1454–1468.https://doi.org/10.1016/j.bpj.2021.01.038
-
The effects of temperature on the stability of a neuronal oscillatorPLOS Computational Biology 9:e1002857.https://doi.org/10.1371/journal.pcbi.1002857
-
A new correlation-based measure of spike timing reliabilityNeurocomputing 52–54:925–931.https://doi.org/10.1016/S0925-2312(02)00838-X
-
BookActive Learning Literature Survey Doctoral DissertationUniversity of Wisconsin-Madison.
-
ConferenceOrdered treemap layoutsIEEE Symposium on Information Visualization, 2001. pp. 73–78.https://doi.org/10.1109/INFVIS.2001.963283
-
Modulators with Convergent Cellular Actions Elicit Distinct Circuit OutputsThe Journal of Neuroscience 21:4050–4058.
-
Robustness of a Rhythmic Circuit to Short- and Long-Term Temperature ChangesThe Journal of Neuroscience 32:10075–10085.https://doi.org/10.1523/JNEUROSCI.1443-12.2012
-
Colocalized Neuropeptides Activate a Central Pattern Generator by Acting on Different Circuit TargetsThe Journal of Neuroscience 22:1874–1882.
-
Neuromodulatory inputs maintain expression of a lobster motor pattern-generating network in a modulation-dependent state: evidence from long-term decentralization in vitroThe Journal of Neuroscience 18:2212–2225.
-
Revealing network connectivity from response dynamicsPhysical Review Letters 98:224101.https://doi.org/10.1103/PhysRevLett.98.224101
-
Modulation of identified stomatogastric ganglion neurons in primary cell cultureJournal of Neurophysiology 69:1993–2002.https://doi.org/10.1152/jn.1993.69.6.1993
-
Activity-dependent changes in the intrinsic properties of cultured neuronsScience (New York, N.Y.) 264:974–977.https://doi.org/10.1126/science.8178157
-
Selective regulation of current densities underlies spontaneous changes in the activity of cultured neuronsThe Journal of Neuroscience 15:3640–3652.
-
Chaos in neuronal networks with balanced excitatory and inhibitory activityScience (New York, N.Y.) 274:1724–1726.https://doi.org/10.1126/science.274.5293.1724
-
Computation Through Neural Population DynamicsAnnual Review of Neuroscience 43:249–275.https://doi.org/10.1146/annurev-neuro-092619-094115
Decision letter
-
Ronald L CalabreseSenior and Reviewing Editor; Emory University, United States
-
Alex WilliamsReviewer; Stanford University, United States
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Decision letter after peer review:
[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]
Thank you for submitting the paper "Mapping Circuit Dynamics During Function and Dysfunction" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Alex Williams (Reviewer #1).
We are sorry to say that, after consultation with the reviewers, we have decided that this work will not be considered further for publication by eLife.
The reviewers find that paper has intrinsic interest especially to specialists in motor pattern generation by small circuits but requires considerable work and refocusing to make it of broader significance.
The expert reviews are quite compatible and provide rich feedback that could lead to a major overhaul and a new submission to eLife.
1. t-SNR is not in itself a clustering algorithm so its resulting 'clusters' should be validated by application of a true clustering algorithm. t–SNE does often provide a nice 2–D embedding that can simplify subsequent clustering. Here the reviewers suggest that the feature vectors themselves should be subjected to clustering algorithm.
2. The feature vectors are not as well motivated and explained as they should be so that readers in other systems can use the strategy presented here in cases where other features than ISI and spike phase of the activity might be more relevant.
3. There was discussion as to whether the framing of the manuscript as a ML pipeline is appropriate because the pipeline described relies heavily on manual curation. The manual curation should be more explicitly described (with details), and the limitations of this approach – ML followed by manual curation – discussed. Such a pipeline probably can't be as easily downloaded and applied directly to someone else's data, but overall, the community could benefit from more examples/consideration of how to systematically and carefully interleave the use of machine–learning techniques with manual analysis by domain experts but the process should be more explicitly explained.
4. Discussion should be revamped to discuss more fully specific insights into the pyloric circuit and the limitations of the analyses presented, for example the omission of the PY neuron activity means that the map as given is incomplete: potentially there are many more states, and hence transitions, within or beyond those already found that correspond to changes in PY neuron activity.
Reviewer #1:
The authors sought to establish a standardized quantitative approach to categorize the activity patterns in a central pattern generator (specifically, the well–studied pyloric circuit in C. borealis). While it is easy to describe these patterns under "normal" conditions, this circuit displays a wide range of irregular behaviors under experimental perturbations. Characterizing and cataloguing these irregular behaviors is of interest to understand how the network avoids these dysfunctional patterns under "normal" circumstances.
The authors draw upon established machine learning tools to approach this problem. To do so, they must define a set of features that describe circuit activity at a moment in time. They use the distribution of inter–spike–intervals ISIs and spike phases of the LP and PD neuron as these features. As the authors mention in their Discussion section, these features are highly specialized and adapted to this particular circuit. This limits the applicability of their approach to other circuits with neurons that are unidentifiable or very large in number (the number of spike phase statistics grows quadratically with the number of neurons).
The main results of the paper provide evidence that ISIs and spike phase statistics provide a reasonable descriptive starting point for understanding the diversity of pyloric circuit patterns. The authors rely heavily on t–distributed stochastic neighbor embedding (tSNE), a well–known nonlinear dimensionality reduction method, to visualize activity patterns in a low–dimensional, 2D space. While effective, the outputs of tSNE have to be interpreted with great care (Wattenberg, et al., "How to Use t–SNE Effectively", Distill, 2016. http://doi.org/10.23915/distill.00002). I think the conclusions of this paper would be strengthened if additional machine learning models were applied to the ISI and spike phase features, and if those additional models validated the qualitative results shown by tSNE. For example, tSNE itself is not a clustering method, so applying clustering methods directly to the high–dimensional data features would be a useful validation of the apparent low–dimensional clusters shown in the figures.
The authors do show that the algorithmically defined clusters agree with expert–defined clusters. (Or, at least, they show that one can come up with reasonable post–hoc explanations and interpretations of each cluster). The very large cluster of "regular" patterns – shown typically in a shade of blue – actually looks like an archipelago of smaller clusters that the authors have reasoned should be lumped together. Thus, while the approach is still a useful data–driven tool, a non–trivial amount of expert knowledge is baked into the results. A central challenge in this line of research is to understand how sensitive the outcomes are to these modeling choices, and there is unlikely to be a definitive answer.
Nonetheless, the authors show results which suggest that this analysis framework may be useful for the community of researchers studying central pattern generators. They use their method to qualitatively characterize a variety of network perturbations – temperature changes, pH changes, decentralization, etc.
In some cases it is difficult to understand the level of certainty in these qualitative observations. A first look at Figure 5a suggests that three different kinds of perturbations push the circuit activity into different dysfunctional cluster regions. However, the apparent spatial differences between these three groups of perturbations might be due to animal–level differences (i.e. each preparation produces multiple points in the low–D plot, so the number of effective statistical replicates is smaller than it appears at first glance). Similarly, in Figure 9, it is somewhat hard to understand how much the state occupancy plots would change if more animals were collected –– with the exception of proctolin, there are ~25 animals and 12 circuit activity clusters which may not be a favorable ratio. It would be useful if a principled method for computing "error bars" on these occupancy diagrams could be developed. Similar "error bars" on the state transition diagrams (e.g. Figure 6a) would also be useful.
Finally, one nagging concern that I have is that the ISIs and spike phase statistics aren't the ideal features one would use to classify pyloric circuit behaviors. Sub–threshold dynamics are incredibly important for this circuit (e.g. due to electrical coupling of many neurons). A deeper discussion about what is potentially lost by only having access to the spikes would be useful.
Overall, I think this work provides a useful starting point for large–scale quantitative analysis of CPG circuit behaviors, but there are many additional hurdles to be overcome.
Reviewer #2:
This manuscript uses the t–SNE dimensionality reduction technique to capture the rich dynamics of the pyloric circuit of the crab.
Strengths:
– The integration of a rich data–set of spiking data from the pyloric circuit.
– Use of nonlinear dimension reduction (t–SNE) to visualise that data.
– Use of clusters from that t–SNE visualisation to create subsets of data that are amenable to consistent analyses (such as using the "regular" cluster as a basis for surveying the types of dynamics possible in baseline conditions).
– Innovative use of the cluster types to describe transitions between dynamics within the baseline state and within perturbed states (whether by changes to exogenous variables, cutting nerves, or applying neuromodulators).
Some interesting main results:
– Baseline variability in the spiking patterns of the pyloric circuit is greater within than between animals.
– Transitions to silent states often (always?) pass through the same intermediate state of the LP neuron skipping spikes.
Weaknesses:
– t-SNE is not, in isolation, a clustering algorithm, yet here it is treated as such. How the clusters were identified is unclear: the manuscript mentions manual curation of randomly sampled points, implying that the clusters were extrapolations from these. This would seem to rather defeat the point of using unsupervised techniques to obtain an unbiased survey of the spiking dynamics and raises the issue of how robust the clusters are.
– The main purpose and contribution of the paper is unclear, as the results are descriptive, and mostly state that dynamics in some vary between different states of the circuit; while the collated dataset is a wonderful resource, and the map is no doubt useful for the lab to place in context what they are looking at, it is not clear what we learn about the pyloric circuit, or more widely about the dynamical repertoire of neural circuits.
– In some places the contribution is noted as being the pipeline of analysis: unfortunately as the pipeline used here seems to rely in manual curation, it is of limited general use; moreover, there are already a number of previous works that use unsupervised machine–learning pipelines to characterise the complexity of spiking activity across a large data–set of neurons, using the same general approach here (quantify properties of spiking as a vector; map/cluster using dimension reduction), including Baden et al. (2016, Nature), Bruno et al. (2015, Neuron), Frady et al. (2016, Neural Computation).
Some key limitations are not considered:
– The omission of the PY neuron activity means that the map as given is incomplete: potentially there are many more states, and hence transitions, within or beyond those already found that correspond to changes in PY neuron activity.
– The use of long, non–overlapping time segments (20s) – this means, for example, that the transitions are slow and discrete, whereas in reality they may be abrupt, or continuous.
– tSNE cannot capture hierarchical structure, nor has a null model to demonstrate that the underlying data contains some clustering structure. So, for example, distances measured on the map may not be strictly meaningful if the data is hierarchical.
– The Discussion does not include enough insight and contextualisation of the results.
Recommendations:
Explain and validate the clusters:
– The paper explains only that the clusters were assigned by manual curation of randomly sampled points. How then were all points assigned to clusters given that sample? What metric was used?
– We would suggest validating the clusters by using a clustering algorithm to recover them (approximately) from the vectors
Clarify the paper's goals and contributions:
– As noted above, the use of an unsupervised pipeline to map spiking activity is not novel; the survey of data is interesting, but purely descriptive. The key contributions were then not obvious to this reader – please clarify. To give some suggestions, while reading 3 possible contributions came to mind, but none quite fit: if the paper is an announcement of the dataset, then the datasets need a detailed explanation; if the paper is the introduction of a pipeline, then the pipeline ought to be fully unsupervised, else it is of little use outside the hands of the present lab; if the paper is about insights into the pyloric circuit, then conclusions and insights ought to be drawn from the descriptive results
Improve the Discussion:
– Explicitly link the results to insights into the pyloric circuit. The abstract states "strong mechanistically interpretable links between complex changes in the collective behavior of a neural circuit and specific experimental manipulations" – but only the short section on "Diversity and stereotypy in trajectories from functional to crash states" touches on this.
– Long section on other clustering methods is a survey of some alternatives and draws no conclusions about why tSNE was chosen here: either use to justify tSNE, or omit.
– Lines 678–716 are not based on any results in the paper.
Discuss limitations:
– The omission of the PY neuron activity means that the map as given is incomplete – what might it look like with it (new states, new transitions etc)?
– The transitions are slow and discrete by design (20 s long segments), whereas in reality they may be abrupt, or continuous.
– tSNE's limitations and their implications e.g. it cannot capture hierarchical structure, nor has a null model to demonstrate that the underlying data contains some clustering structure.
– Why those features of the spike–trains, and how would the map change if features were omitted or added? (e.g. the regularity of spiking).
Reviewer #3:
Gorur–Shandilya et al. apply an unsupervised dimensionality reduction (t–SNE) to characterize neural spiking dynamics in the pyloric circuit in the stomatogastric ganglion of the crab. The application of unsupervised methods to characterize qualitatively distinct regimes of spiking neural circuits is very interesting and novel, and the manuscript provides a comprehensive demonstration of its utility by analyzing dynamical variability in function and dysfunction in an important rhythm–generating circuit. The system is highly tractable with small numbers of neurons, and the study here provides an important new characterization of the system that can be used to further understand the mapping between gene expression, circuit activity, and functional regimes. The explicit note about the importance of visualization and manual labeling was also nice, since this is often brushed under the rug in other studies.
While the specific analysis pipeline clearly identifies qualitatively distinct regimes of spike patterns in the LP/PD neurons, it is not clear how much of this is due to t–SNE itself vs the initial pre–processing and feature definition (ISI and spike phase percentiles). Analyses that would help clarify this would be to check whether the same clusters emerge after (1) applying ordinary PCA to the feature vectors and plotting the projections of the data along the first two PCs, or (2) defining input features as the concatenated binned spike rates over time of the LP & PD neurons (which would also yield a fixed–length vector per 20 s trial), and then passing these inputs to PCA or t–SNE. As the significance of this work is largely motivated by using unsupervised vs ad hoc descriptors of circuit dynamics, it will be important to clarify how much of the results derive from the use of ISI and phase representation percentiles, etc. as input features, vs how much emerge from the dimensionality reduction.
Please define all acronyms when they are introduced.
Why is 20 seconds chosen as the time window to compute ISIs/phase representations before feeding the data into t–SNE? Is this roughly the window one needs to look at to visually identify differences across dynamical regimes and can the authors explain why?
It would be useful to note in the manuscript that t–SNE is tuned to preserve local over global similarities, which gives it its utility in clustering. It would also be useful to discuss the relationship of t–SNE to UMAP (McInnes, Healy, Melville 2018 arXiv: https://arxiv.org/abs/1802.03426), another popular dimensionality reduction technique in neuroscience.
Most uses of t–SNE in biology in the past have been in the context of classifying behavior (e.g. Berman et al. 2014 J R Soc Interface: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4233753/ or Clemens et al. Current Biology 2018: https://www.sciencedirect.com/science/article/pii/S0960982218307735). Including these references would be helpful as a point of comparison and to emphasize the novelty of applying t–SNE to neural spiking data.
Can the authors comment on how much pyloric rhythms can deviate from the standard triphasic pattern before behavior (or digestion?) is significantly disrupted?
In Figure 6 it would be helpful to show example time–series going into the transition matrix analysis, with identified state labels at the different timepoints. It wasn't clear whether timepoints were e.g. moment–to–moment, 20 sec chunks, a 20 sec sliding window, etc.
In Figure 8 it would be helpful to show examples of what bursts look like with and without decentralization, in order to contextualize the quantitative metrics.
It would also be useful to connect this work in the discussion to modern approaches for identifying model parameter regimes underlying distinct neural activity patterns (see e.g. Bittner et al. 2021 bioRxiv https://www.biorxiv.org/content/10.1101/837567v3)
[Editors’ note: further revisions were suggested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "Mapping circuit dynamics during function and dysfunction" for further consideration by eLife. Your revised article has been evaluated by Ronald Calabrese (Senior Editor) and a Reviewing Editor.
The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:
Two reviewer concerns remain, which are explained more fully in the reviews below.
1. The analysis pipeline outlined in this paper is rather ad hoc. This isn't a general "tool" as it has only been demonstrated in one particular biological system.
2. It is acceptable for the authors to call the core method (tSNE) unsupervised. The "manual curation" is part of any exploratory data analysis, though it is generally preferable to for this manual curation / interpretation step to be as objective and reproducible as possible (which is arguably not entirely achieved in this paper). The paper can be seen as an investigation into how to apply an existing tool (tSNE) to a very well-studied and simple neural system.
In revision, we ask the authors to take a final pass through the Introduction and Discussion and to tone down any comments about this being a general, off-the-shelf tool and to note that while the core algorithmic method (tSNE) is unsupervised, there is still the need for expert interpretation of the output.
Reviewer #1:
The manuscript has been revised along the lines suggested in both the consensus suggestions and the individual reports. The rewritten Discussion is considerably improved. There is a lot of hard work here, and I've no doubt the results are useful to the PI's lab, and presumably other workers on this CPG.
To my mind the authors did not really address two of the main concerns of the reviewers: (1) what the purpose of this paper is and (2) the validity of the manual clustering of tSNE. In more detail:
Purpose of the paper: the paper prominently refers to it introducing an "unsupervised" method/approach (e.g. lines 69-70) and a "tool" (lines 537-541, 711-714). Neither are true: the approach here is by definition not unsupervised as it uses manual curation to identify the clusters upon which all further analyses are based; such manual curation means that the only "tool" parts are using tSNE on a feature vector, which doesn't merit the term. Rather, it seems that the purpose of the paper is instead to provide a cohesive overview of the dynamics of this specific CPG system, by finding a way to systematically combine a large number of separate recordings from baseline and perturbed preparations, thus allowing knowledge discovery (e.g. the types of transitions; the changes in the regular state that precede transitions) on those combined data.
tSNE and validating the clustering: the authors did considerable work related to this issue, by using their approach on synthetic data, replacing tSNE with UMAP, and by testing simple hierarchical and k-means clustering of their feature vectors. But none of these spoke directly to the validity of the clustering: the synthetic data and UMAP were still manually clustered, so could be made to conform to the desired results; the unsupervised clustering approach used shows the well-known problem that directly clustering a high-dimensional dataset (a 48-dimensional feature vector) is often impossible using classical techniques simply due to the curse of dimensionality. The normal approach here would be to use dimension reduction then unsupervised clustering (e.g. in the classic combination of PCA then k-means), matching the way the two steps are done in the present paper.
As the paper shows, the clusters found in tSNE by manual curation are arbitrary. For example, changing the perplexity parameter changes the fragmentation of the tSNE map, so would likely lead to different manual clusterings. In another example, the clusters found by manual curation are not all contiguous in the 2D space e.g. the LP-silent cluster is in 3 groups, one of which is on the opposite side of the 2D space to the others – so distance in this space is not interpreted as meaningful for all clusters. In another example, the dividing line between some "clusters" is arbitrarily put in a contiguous run of points – e.g. the LP-weak-skipped (pink) and irregular-bursting (brown).
So to me it seems we're still left with the original main issues raised, and consequently am still unclear who this paper is aimed at – this is a useful way for the lab who owns that data to organise it, but how useful is it in revealing the structure of dynamics in the CPG, and hence to others? Moreover, though the Discussion is much improved, the paper is still rather descriptive throughout.
These are clearly subjective concerns, not issues with the quality of the analyses done based on the clustering, which are solid, nor with their interpretation, which has been considerably improved. Nonetheless, I think the paper would still be considerably improved by a concerted redraft, and perhaps reorganisation, to deal with these issues.
Reviewer #2:
I appreciate that the authors carefully considered my review and tried new analyses (e.g. hierarchical clustering and UMAP) and have expanded the Discussion section to comment on the issues raised during review. I also agree that "a rigorous error analysis would be useful but is not trivially done" -- and since I do not have a concrete suggestion here, I think the manuscript should be published without one. Overall, I feel that my concerns have been satisfactorily addressed.
Reviewer #3:
One small remaining edit that would be helpful is to show the projection of the data onto the first two PCs (re: the PCA analysis referenced in Figure 2-Figure supp 2), so that the reader can visualize what structure there is prior to computing the triadic difference statistics. If not too computationally intensive, it would also be useful to show what emerges from applying a naive clustering or visualization algorithm directly on the spike-train time-series [see (2) in Reviewer #3 major concern], even if only to show that it is quite messy.
https://doi.org/10.7554/eLife.76579.sa1Author response
[Editors’ note: the authors resubmitted a revised version of the paper for consideration. What follows is the authors’ response to the first round of review.]
1. t-SNR is not in itself a clustering algorithm so its resulting 'clusters' should be validated by application of a true clustering algorithm. t–SNE does often provide a nice 2–D embedding that can simplify subsequent clustering. Here the reviewers suggest that the feature vectors themselves should be subjected to clustering algorithm.
We agree with the reviewers is t-SNE is not in itself a clustering algorithm. We have been careful to use it as a visualization and data reduction tool, whose strengths are uniquely suited to the task we applied it to. The only clustering we performed in the original manuscript was entirely manual, using the t-SNE embedding to visualize the data and a tool to quickly inspect the raw spike patterns to guide our clustering. As the reviewers point out, “t-SNE does often provide a nice 2-D embedding that can simplify subsequent clustering”. This is exactly what we did. The strength of our approach is that it made feasible a previously intractable task: manually classifying hundreds of hours of spike patterns from hundreds of animals.
We agree with the reviewers that this paper can be strengthened by validation of this technique and by also comparing to alternative techniques. We therefore did the following:
1. We created a synthetic data set with spikes from two neurons drawn from different patterns. We then subjected this data set to the workflow described in the manuscript and found we could recover distinct classes of spike patterns that existed in the synthetic data. This new analysis is shown in Figure 2—figure supplement 4.
2. In addition, we used another nonlinear dimensionality reduction technique (uMAP) to also embed the feature vectors in two-dimensional space. We do this to emphasize that t-SNE is only a tool to visualize patterns in high dimensional data and is not the only one that can be useful, and that clustering of this data comes not from t-SNE but is a distinct step following embedding (Figure 3— figure supplement 3).
3. Finally, we use hierarchical clustering to directly cluster the feature vectors. In Author response image 1, (a) shows how the number of clusters varies with the cutoff chosen. Allowing the cutoff to vary to yield a maximum of 20 clusters, we color the points in the t-SNE embedding by cluster identity (b). The vast majority of the data are grouped into a single cluster, with disparate spike patterns grouped together. This is likely because we used Euclidean distances to determine linkage between clusters, and clustering may be dominated by global rather than local structure in the feature vectors. A similar picture emerges when we use kmeans clustering on the feature vectors (c) Setting the number of clusters k = 20. Cluster boundaries do not correspond to sharp discontinuities in spike patterns, as observed by manual inspection. (d) shows the distribution of points in each cluster. These results demonstrate that the methods we propose offer utility beyond naïve hierarchical or k-means clustering and are better suited to recover patterns in neural data.
Author response image 1
2. The feature vectors are not as well motivated and explained as they should be so that readers in other systems can use the strategy presented here in cases where other features than ISI and spike phase of the activity might be more relevant.
We have expanded the description of the feature vectors in the Methods:
Construction of vectorized data frame to go into greater detail. We have included a discussion of why we chose these feature vectors, and proposed alternatives to ease applicability in other systems.
3. There was discussion as to whether the framing of the manuscript as a ML pipeline is appropriate, because the pipeline described relies heavily on manual curation. The manual curation should be more explicitly described (with details), and the limitations of this approach – ML followed by manual curation – discussed. Such a pipeline probably can't be as easily downloaded and applied directly to someone else's data, but overall, the community could benefit from more examples/consideration of how to systematically and carefully interleave the use of machine–learning techniques with manual analysis by domain experts but the process should be more explicitly explained.
The method described here uses a nonlinear dimensionality reduction algorithm to generate a two-dimensional visualization that allows a domain expert to efficiently label and cluster large amounts of data. We agree that this is a manual process that requires human intervention. We have extensively rewritten the Discussion to take this into account and included a section “Applicability to other systems” that goes into detail what needs to be done to apply this method to other circuits. We have included a statement in the Discussion making clear that a drawback of this method is that it is not fully automated and requires a human to inspect data and draw cluster boundaries (lines 532-541).
4. Discussion should be revamped to discuss more fully specific insights into the pyloric circuit and the limitations of the analyses presented, for example the omission of the PY neuron activity means that the map as given is incomplete: potentially there are many more states, and hence transitions, within or beyond those already found that correspond to changes in PY neuron activity.
We agree with the reviewers that the omission of the PY neurons’ activity means that the map is incomplete. There are likely many more states, and hence many more transitions, than the ones we have identified. In addition, we note that there are other pyloric neurons whose activity is also missing (AB, IC, LPG, VD). However, measuring just LP and PD allows us to monitor the activity of the most important functional antagonists in the system (because they are effectively in a half-center oscillator because PD is electrically coupled to AB). In general, the more neurons one measures, the richer the description of the circuit dynamics will be, and the more transitions one will observe. Collecting datasets at this scale (~500 animals) from all pyloric neurons is challenging, and we have revised the manuscript to make this important point (see Discussion: Technical considerations).
Reviewer #1:
The authors sought to establish a standardized quantitative approach to categorize the activity patterns in a central pattern generator (specifically, the well–studied pyloric circuit in C. borealis). While it is easy to describe these patterns under "normal" conditions, this circuit displays a wide range of irregular behaviors under experimental perturbations. Characterizing and cataloguing these irregular behaviors is of interest to understand how the network avoids these dysfunctional patterns under "normal" circumstances.
The authors draw upon established machine learning tools to approach this problem. To do so, they must define a set of features that describe circuit activity at a moment in time. They use the distribution of inter–spike–intervals ISIs and spike phases of the LP and PD neuron as these features. As the authors mention in their Discussion section, these features are highly specialized and adapted to this particular circuit. This limits the applicability of their approach to other circuits with neurons that are unidentifiable or very large in number (the number of spike phase statistics grows quadratically with the number of neurons).
We agree with the reviewer that the size of the feature vectors as described grows quadratically with the number of neurons. The feature sets we describe are most suited for “identified” neurons – neurons whose identity and connectivity are known and can be reliably recorded from multiple animals. The method described here is best suited for systems with small numbers of identified neurons. For other systems, other feature vectors may be chosen, as we have suggested in the Discussion: Applicability to other systems.
The main results of the paper provide evidence that ISIs and spike phase statistics provide a reasonable descriptive starting point for understanding the diversity of pyloric circuit patterns. The authors rely heavily on t–distributed stochastic neighbor embedding (tSNE), a well–known nonlinear dimensionality reduction method, to visualize activity patterns in a low–dimensional, 2D space. While effective, the outputs of tSNE have to be interpreted with great care (Wattenberg, et al., "How to Use t–SNE Effectively", Distill, 2016. http://doi.org/10.23915/distill.00002). I think the conclusions of this paper would be strengthened if additional machine learning models were applied to the ISI and spike phase features, and if those additional models validated the qualitative results shown by tSNE. For example, tSNE itself is not a clustering method, so applying clustering methods directly to the high–dimensional data features would be a useful validation of the apparent low–dimensional clusters shown in the figures.
We thank the reviewer for these suggestions and agree with the reviewer that t-SNE is not a clustering method, and directly clustering on t-SNE embeddings is rife with complexities. Instead, we have used t-SNE to generate a visualization that allows domain experts to quickly label and cluster large quantities of data. This makes a previously intractable task feasible, and offers some basic guarantees on quality (e.g., no one data point can have two labels, because labels derive from position of data points in two-dimensional space). In addition:
– We used uMAP, another dimensionality reduction algorithm, to perform the embedding step, and colored points by the original t-SNE embedding. (Figure 3—figure supplement 3). Large sections of the map are still strikingly colored in single colors, suggesting that the manual clustering did not depend on the details of the t-SNE algorithm, but is rather informed by the statistics of the data.
We validated our method using synthetic data. We generated synthetic spike trains from different “classes” and embedded the resultant feature vectors using t-SNE. Data from different classes are not intermingled and form tight “clusters” (Figure 2 —figure supplement 4).
– Finally, we attempted to use hierarchical clustering to cluster the raw feature vectors and were not able to find a reasonable portioning of the linkage tree that separated qualitatively different spike patterns (Figure at the top of this document). We speculate that this is because feature vectors may contain outliers that bias clustering algorithms that attempt to preserve global distance to lump the majority of the data into a single cluster, in order to differentiate outliers from the bulk of the data.
The authors do show that the algorithmically defined clusters agree with expert–defined clusters. (Or, at least, they show that one can come up with reasonable post–hoc explanations and interpretations of each cluster). The very large cluster of "regular" patterns – shown typically in a shade of blue – actually looks like an archipelago of smaller clusters that the authors have reasoned should be lumped together. Thus, while the approach is still a useful data–driven tool, a non–trivial amount of expert knowledge is baked into the results. A central challenge in this line of research is to understand how sensitive the outcomes are to these modeling choices, and there is unlikely to be a definitive answer.
We agree with the reviewer entirely.
Nonetheless, the authors show results which suggest that this analysis framework may be useful for the community of researchers studying central pattern generators. They use their method to qualitatively characterize a variety of network perturbations – temperature changes, pH changes, decentralization, etc.
In some cases it is difficult to understand the level of certainty in these qualitative observations. A first look at Figure 5a suggests that three different kinds of perturbations push the circuit activity into different dysfunctional cluster regions. However, the apparent spatial differences between these three groups of perturbations might be due to animal–level differences (i.e. each preparation produces multiple points in the low–D plot, so the number of effective statistical replicates is smaller than it appears at first glance). Similarly, in Figure 9, it is somewhat hard to understand how much the state occupancy plots would change if more animals were collected –– with the exception of proctolin, there are ~25 animals and 12 circuit activity clusters which may not be a favorable ratio. It would be useful if a principled method for computing "error bars" on these occupancy diagrams could be developed. Similar "error bars" on the state transition diagrams (e.g. Figure 6a) would also be useful.
We agree with the reviewer. Despite this paper containing data from hundreds of animals, the dataset may not be sufficiently large to perform some necessary statistical checks. We agree with the reviewer that a more rigorous error analysis would be useful but is not trivially done.
Finally, one nagging concern that I have is that the ISIs and spike phase statistics aren't the ideal features one would use to classify pyloric circuit behaviors. Sub–threshold dynamics are incredibly important for this circuit (e.g. due to electrical coupling of many neurons). A deeper discussion about what is potentially lost by only having access to the spikes would be useful.
We agree with the reviewer that spike times aren’t the ideal feature to use to describe circuit dynamics. This is especially true in the STG, where synapses are graded, and coupling between cells can persist without spiking. However, the data required simply do not exist, as it requires intracellular recordings, which are substantially harder to perform (and maintain over challenging perturbations) than extracellular recordings.
Finally, the signal to the muscles – arguably the physiologically and functionally relevant signal – is the spike signal, suggesting that spike patterns from the pyloric circuit are a useful feature to measure. Nevertheless, this is an important point, and we thank the reviewer for raising it, and we have included it in the section titled Discussion: Technical considerations.
Overall, I think this work provides a useful starting point for large–scale quantitative analysis of CPG circuit behaviors, but there are many additional hurdles to be overcome.
Reviewer #2:
This manuscript uses the t–SNE dimensionality reduction technique to capture the rich dynamics of the pyloric circuit of the crab.
Strengths:
– The integration of a rich data–set of spiking data from the pyloric circuit.
– Use of nonlinear dimension reduction (t–SNE) to visualise that data.
– Use of clusters from that t–SNE visualisation to create subsets of data that are amenable to consistent analyses (such as using the "regular" cluster as a basis for surveying the types of dynamics possible in baseline conditions).
– Innovative use of the cluster types to describe transitions between dynamics within the baseline state and within perturbed states (whether by changes to exogenous variables, cutting nerves, or applying neuromodulators).
Some interesting main results:
– Baseline variability in the spiking patterns of the pyloric circuit is greater within than between animals.
– Transitions to silent states often (always?) pass through the same intermediate state of the LP neuron skipping spikes
Weaknesses:
– t-SNE is not, in isolation, a clustering algorithm, yet here it is treated as such. How the clusters were identified is unclear: the manuscript mentions manual curation of randomly sampled points, implying that the clusters were extrapolations from these. This would seem to rather defeat the point of using unsupervised techniques to obtain an unbiased survey of the spiking dynamics and raises the issue of how robust the clusters are.
We have used t-SNE to visualize the circuit dynamics in a two-dimensional map. We have exploited t-SNE’s ability to preserve local structure to generate an embedding where a domain expert can efficiently manually identify and label stereotyped clusters of activity. As the author points out, this is a manual step, and we have emphasized this in the manuscript. The strength of our approach is to combine the power of a nonlinear dimensionality reduction technique such as t-SNE with human curation to make a task that was previously impossible (identifying and labelling very large datasets of neural activity) feasible.
To address the question of how robust the manually identified clusters are, we have:
1. Used another dimensionality reduction technique, uMAP, to generate an embedding and colored points by the original t-SNE map (Figure 3 —figure supplement 3). To rough approximation, the coloring reveals that a similar clustering exists in this uMAP embedding.
2. We generated synthetic spike trains from pre-determined spike pattern classes and used the feature vector extraction and t-SNE embedding procedure as described in the paper. We found that this generated a map (Figure 2—figure supplement 4) where classes of spike patterns were well separated in the t-SNE space.
– The main purpose and contribution of the paper is unclear, as the results are descriptive, and mostly state that dynamics in some vary between different states of the circuit; while the collated dataset is a wonderful resource, and the map is no doubt useful for the lab to place in context what they are looking at, it is not clear what we learn about the pyloric circuit, or more widely about the dynamical repertoire of neural circuits.
– In some places the contribution is noted as being the pipeline of analysis: unfortunately as the pipeline used here seems to rely in manual curation, it is of limited general use; moreover, there are already a number of previous works that use unsupervised machine–learning pipelines to characterise the complexity of spiking activity across a large data–set of neurons, using the same general approach here (quantify properties of spiking as a vector; map/cluster using dimension reduction), including Baden et al. (2016, Nature), Bruno et al. (2015, Neuron), Frady et al. (2016, Neural Computation).
Some key limitations are not considered:
– The omission of the PY neuron activity means that the map as given is incomplete: potentially there are many more states, and hence transitions, within or beyond those already found that correspond to changes in PY neuron activity.
We agree with the reviewer that the omission of the PY neurons’ activity means that the map is incomplete. There are likely many more states, and hence many more transitions, than the ones we have identified. In addition, we note that there are other pyloric neurons whose activity is also missing (AB, IC, LPG, VD). However, measuring just LP and PD allows us to monitor the activity of the most important functional antagonists in the system (because they are effectively in a half-center oscillator because PD is electrically coupled to AB). In general, the more neurons one measures, the richer the description of the circuit dynamics will be. Collecting datasets at this scale (~500 animals) from all pyloric neurons is challenging, and we have revised the manuscript to make this important point (see Discussion: Technical considerations).
– The use of long, non–overlapping time segments (20s) – this means, for example, that the transitions are slow and discrete, whereas in reality they may be abrupt, or continuous.
We agree with the reviewer. There are tradeoffs in choosing a bin size in analyzing time series – choosing longer bins can increase the number of “states” and choosing shorter bins can increase the number of transitions. We chose 20s bins because it is long enough to include several cycles of the pyloric rhythm, even when decentralized, yet was short enough to resolve slow changes in spiking. We have included a statement clarifying this (see Discussion: Technical considerations).
– tSNE cannot capture hierarchical structure, nor has a null model to demonstrate that the underlying data contains some clustering structure. So, for example, distances measured on the map may not be strictly meaningful if the data is hierarchical.
We agree with the reviewer. t-SNE can manifest clusters when none exist (Section 4 of https://distill.pub/2016/misread-tsne/) and can obscure or merge true clusters. We have restricted analyses that rely on distances measured in the map to cases where there are qualitative differences in behavior (e.g., with decentralization, Figure 7) or have compared distances within subsets of data where a single parameter is changed (e.g., pH or temperature, Figure 5). The only conclusion we draw from these distance measures is that data are more (or less) spread out in the map, which we use as a proxy for variability. We have included a statement, discussion limitations of using t-SNE (Discussion: Comparison with other methods).
– The Discussion does not include enough insight and contextualisation of the results.
We have completely rewritten the discussion to address this.
Recommendations:
Explain and validate the clusters:
– The paper explains only that the clusters were assigned by manual curation of randomly sampled points. How then were all points assigned to clusters given that sample? What metric was used?
We have added a new section in the Methods (Manual clustering and annotation of data) that describes this process.
– We would suggest validating the clusters by using a clustering algorithm to recover them (approximately) from the vectors.
We have validated the clusters and the embedding algorithm by:
1. Creating a synthetic data set with spikes from two neurons drawn from different patterns. We then subjected this data set to the workflow described in the manuscript and found we could recover distinct classes of spike patterns that existed in the synthetic data. This new analysis is shown in Figure 2—figure supplement 4.
2. In addition, we used another nonlinear dimensionality reduction technique (uMAP) to also embed the feature vectors in two-dimensional space. We do this to emphasize that t-SNE is only a tool to visualize patterns in high dimensional data and is not the only one that can be useful, and that clustering of this data comes not from t-SNE but is a distinct step following embedding. (Figure 3— figure supplement 3)
3. Finally, we use hierarchical clustering to directly cluster the feature vectors. Author response image 1, (a) shows how the number of clusters varies with the cutoff chosen. Allowing the cutoff to vary to yield a maximum of 20 clusters, we color the points in the t-SNE embedding by cluster identity (b). The vast majority of the data is grouped into a single cluster, with disparate spike patterns grouped together. This is likely because we used Euclidean distances to determine linkage between clusters, and clustering may be dominated by global rather than local structure in the feature vectors. A similar picture emerges when we use kmeans clustering on the feature vectors (c) setting the number of clusters k = 20. Cluster boundaries do not correspond to sharp discontinuities in spike patterns, as observed by manual inspection. (d) shows the distribution of points in each cluster. These results demonstrate that the methods we propose offer utility beyond naïve hierarchical or k-means clustering and are better suited to recover patterns in neural data.
Clarify the paper's goals and contributions:
– As noted above, the use of an unsupervised pipeline to map spiking activity is not novel; the survey of data is interesting, but purely descriptive. The key contributions were then not obvious to this reader – please clarify. To give some suggestions, while reading 3 possible contributions came to mind, but none quite fit: if the paper is an announcement of the dataset, then the datasets need a detailed explanation; if the paper is the introduction of a pipeline, then the pipeline ought to be fully unsupervised, else it is of little use outside the hands of the present lab; if the paper is about insights into the pyloric circuit, then conclusions and insights ought to be drawn from the descriptive results
Improve the Discussion:
– Explicitly link the results to insights into the pyloric circuit. The abstract states "strong mechanistically interpretable links between complex changes in the collective behavior of a neural circuit and specific experimental manipulations" – but only the short section on "Diversity and stereotypy in trajectories from functional to crash states" touches on this.
We have rewritten the entire Discussion. A new section (Linking circuit outputs to circuit mechanisms) reviewers links between our results and what we learn about the pyloric circuit.
– Long section on other clustering methods is a survey of some alternatives and draws no conclusions about why tSNE was chosen here: either use to justify tSNE, or omit.
We have removed this section.
– Lines 678–716 are not based on any results in the paper
We have removed this section.
Discuss limitations:
– The omission of the PY neuron activity means that the map as given is incomplete – what might it look like with it (new states, new transitions etc)?
We agree with the reviewer that the omission of the PY neurons’ activity means that the map is incomplete. There are likely many more states, and hence many more transitions, than the ones we have identified. In addition, we note that there are other pyloric neurons whose activity is also missing (AB, IC, LPG, VD). However, measuring just LP and PD allows us to monitor the activity of the most important functional antagonists in the system (because they are effectively in a half-center oscillator because PD is electrically coupled to AB). In general, the more neurons one measures, the richer the description of the circuit dynamics will be. Collecting datasets at this scale (~500 animals) from all pyloric neurons is challenging, and we have revised the manuscript to make this important point (see Discussion: Technical considerations).
– The transitions are slow and discrete by design (20 s long segments), whereas in reality they may be abrupt, or continuous.
We agree with the reviewer. There are tradeoffs in choosing a bin size in analyzing time series – choosing longer bins can increase the number of “states” and choosing shorter bins can increase the number of transitions. We chose 20s bins because it is long enough to include several cycles of the pyloric rhythm, even when decentralized, yet was short enough to resolve slow changes in spiking. We have added a new section called “Technical considerations” to make this important point.
– tSNE's limitations and their implications e.g. it cannot capture hierarchical structure, nor has a null model to demonstrate that the underlying data contains some clustering structure.
This is an important point. We have added this to a new section in the Discussion called “Technical considerations”.
– Why those features of the spike–trains, and how would the map change if features were omitted or added? (e.g. the regularity of spiking).
We have repeated the embedding deleting randomly chosen columns in our matrix of feature vectors. This gives us some insight into whether our results are strongly dependent on any one feature. In the following figure, each panel shows the result of embedding following the deletion of a randomly chosen column in the feature matrix. All coloring is from the original map. All maps look similar, and gross features of the map are preserved, suggesting that the map does not depend sensitively on any one feature.
Reviewer #3:
Gorur–Shandilya et al. apply an unsupervised dimensionality reduction (t–SNE) to characterize neural spiking dynamics in the pyloric circuit in the stomatogastric ganglion of the crab. The application of unsupervised methods to characterize qualitatively distinct regimes of spiking neural circuits is very interesting and novel, and the manuscript provides a comprehensive demonstration of its utility by analyzing dynamical variability in function and dysfunction in an important rhythm–generating circuit. The system is highly tractable with small numbers of neurons, and the study here provides an important new characterization of the system that can be used to further understand the mapping between gene expression, circuit activity, and functional regimes. The explicit note about the importance of visualization and manual labeling was also nice, since this is often brushed under the rug in other studies.
While the specific analysis pipeline clearly identifies qualitatively distinct regimes of spike patterns in the LP/PD neurons, it is not clear how much of this is due to t–SNE itself vs the initial pre–processing and feature definition (ISI and spike phase percentiles). Analyses that would help clarify this would be to check whether the same clusters emerge after (1) applying ordinary PCA to the feature vectors and plotting the projections of the data along the first two PCs, or (2) defining input features as the concatenated binned spike rates over time of the LP & PD neurons (which would also yield a fixed–length vector per 20 s trial), and then passing these inputs to PCA or t–SNE. As the significance of this work is largely motivated by using unsupervised vs ad hoc descriptors of circuit dynamics, it will be important to clarify how much of the results derive from the use of ISI and phase representation percentiles, etc. as input features, vs how much emerge from the dimensionality reduction.
We agree with the reviewer that is important to clarify how much of our results come from the data itself, and how we parameterize them using ISIs and phases, and how much comes from the choice of t-SNE as a dimensionality reduction algorithm. We have addressed this concern in the following ways:
1. We used principal components analysis on the feature vectors and measured triadic differences in features such as the period and duty cycle of the PD neuron. We found that triadic differences were lower in the t-SNE embedding than in the first two PCA features, or in shuffled t-SNE embeddings (Figure 2– Figure supplement 2), suggesting that the embedding is creating a useful representation that captures key features of the data.
2. We have used uMAP to reduce the dimensionality of the feature matrix to two dimensions and found that it too preserved the coarse features of the embedding that we observe with t-SNE. Coloring the uMAP embedding by the t-SNE labels revealed that the overall classification scheme was intact (Fig 3 – figure supplement 3).
3. We generated a synthetic dataset and applied the unsupervised part of our algorithm to it (conversion to ISIs, phases, etc., then t-SNE). We colored the points in the t-SNE embedding by the category in the synthetic dataset. We found that categories were well separated in the t-SNE plot, and each cluster tended to have a single color. This validates the overall power of our approach and shows that it can recover clustering information in large spike sets (Figure 2—figure supplement 4).
4. We have run k-means and hierarchical clustering on the feature vectors directly and shown that our method is superior to these naïve clustering algorithms running on the feature vectors. We speculate that this is because these clustering methods attempt to partition the full space using global distances, at the expense of distance along the manifold on which the data is located. Algorithms like t-SNE are biased towards local distances, and discount global distances between points outside a neighborhood, and are this better suited here.
Please define all acronyms when they are introduced.
Fixed.
Why is 20 seconds chosen as the time window to compute ISIs/phase representations before feeding the data into t–SNE? Is this roughly the window one needs to look at to visually identify differences across dynamical regimes and can the authors explain why?
Binning and analyzing data in bins is routine in any time-series analysis. 20s bins were chosen to be large enough that there are rich spike patterns within the bin, and small enough so that long-term changes could be visualized as changes in states. Smaller bins increase the number of points that need to be embedded. We have added a statement about this to the section titled Discussion: Technical considerations.
It would be useful to note in the manuscript that t–SNE is tuned to preserve local over global similarities, which gives it its utility in clustering. It would also be useful to discuss the relationship of t–SNE to UMAP (McInnes, Healy, Melville 2018 arXiv: https://arxiv.org/abs/1802.03426), another popular dimensionality reduction technique in neuroscience.
We have revised the manuscript to include using uMAP to perform the dimensionality reduction instead of t-SNE (Figure 3 – figure supplement 3). As we pointed out in the original manuscript (lines 885), recent work (Kobak and Linderman 2021) has found that t-SNE and uMAP are more similar than previously thought, with differences arising substantially from differences in initialization.
Most uses of t–SNE in biology in the past have been in the context of classifying behavior (e.g. Berman et al. 2014 J R Soc Interface: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4233753/ or Clemens et al. Current Biology 2018: https://www.sciencedirect.com/science/article/pii/S0960982218307735). Including these references would be helpful as a point of comparison and to emphasize the novelty of applying t–SNE to neural spiking data.
We cite Berman et al. and other relevant papers for these reasons (line 551 in original manuscript).
Can the authors comment on how much pyloric rhythms can deviate from the standard triphasic pattern before behavior (or digestion?) is significantly disrupted?
Unfortunately, there are very few studies of motor patterns in intact animals, and it would be very hard to know how much real disruption of food movement would occur, even if the pyloric rhythm were somewhat disrupted. So sadly, we have no idea how to answer this question. That said, the pyloric rhythm is remarkably reliable in vivo under normal conditions.
In Figure 6 it would be helpful to show example time–series going into the transition matrix analysis, with identified state labels at the different timepoints. It wasn't clear whether timepoints were e.g. moment–to–moment, 20 sec chunks, a 20 sec sliding window, etc.
Time series of states are shown in Fig 5—figure supplement 1. Transition analysis was performed taking care only to measure transitions in contiguous chunks of data. All analysis in the manuscript uses 20s non-overlapping bins. We have added a sentence clarifying this in the Results.
In Figure 8 it would be helpful to show examples of what bursts look like with and without decentralization, in order to contextualize the quantitative metrics.
We have included a new figure supplement (Figure 8 – figure supplement 1) showing example traces before and after decentralization.
It would also be useful to connect this work in the discussion to modern approaches for identifying model parameter regimes underlying distinct neural activity patterns (see e.g. Bittner et al. 2021 bioRxiv https://www.biorxiv.org/content/10.1101/837567v3)
This is an interesting topic for further discussion, and is beyond the scope of our paper. In this manuscript we do not try to infer parameters of a mechanistic model that can generate these data; instead we focus on trying to understand the structure of the data and what stereotyped patterns we observe in it.
[Editors’ note: what follows is the authors’ response to the second round of review.]
The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:
Two reviewer concerns remain, which are explained more fully in the reviews below.
1. The analysis pipeline outlined in this paper is rather ad hoc. This isn't a general "tool" as it has only been demonstrated in one particular biological system.
We agree with the reviewers that this isn’t a general “tool” and we have demonstrated the utility of this approach only in one biological system. We have revised all descriptions of our approach to remove the word “tool”. Even though we have limited ourselves to a single circuit type, we have shown that our approach allows us to ask many interesting questions of the system and refer the reader to the Discussion section where we talk about its applicability to other systems.
2. It is acceptable for the authors to call the core method (tSNE) unsupervised. The "manual curation" is part of any exploratory data analysis, though it is generally preferable to for this manual curation / interpretation step to be as objective and reproducible as possible (which is arguably not entirely achieved in this paper). The paper can be seen as an investigation into how to apply an existing tool (tSNE) to a very well-studied and simple neural system.
In revision, we ask the authors to take a final pass through the Introduction and Discussion and to tone down any comments about this being a general, off-the-shelf tool and to note that while the core algorithmic method (tSNE) is unsupervised, there is still the need for expert interpretation of the output.
We have carefully revised the Introduction and Discussion and removed any statements suggesting that our approach is an off-the-shelf tool and have revised its description to be explicit about what we did. We have explicitly stated that our approach requires manual annotation and classification of spike patterns (Introduction lines 68-71; Results – Visualization of circuit dynamics allows manual labelling and clustering of data; Discussion lines 494, 531, 623-631, 676-684).
Reviewer #1:
The manuscript has been revised along the lines suggested in both the consensus suggestions and the individual reports. The rewritten Discussion is considerably improved. There is a lot of hard work here, and I've no doubt the results are useful to the PI's lab, and presumably other workers on this CPG.
To my mind the authors did not really address two of the main concerns of the reviewers: (1) what the purpose of this paper is and (2) the validity of the manual clustering of tSNE. In more detail:
Purpose of the paper: the paper prominently refers to it introducing an "unsupervised" method/approach (e.g. lines 69-70) and a "tool" (lines 537-541, 711-714). Neither are true: the approach here is by definition not unsupervised as it uses manual curation to identify the clusters upon which all further analyses are based; such manual curation means that the only "tool" parts are using tSNE on a feature vector, which doesn't merit the term.
We agree with the reviewer. We have removed all references to “unsupervised” in isolation and made it clear that t-SNE is an unsupervised method, but what we have done requires manual identification and labelling of clusters.
We further agree that the use of “tool” isn’t warranted. We have removed “tool” from all references to what we did and rephrased to make it clear that what we did uses both an off-the-shelf tool (t-SNE) and expert annotation.
Rather, it seems that the purpose of the paper is instead to provide a cohesive overview of the dynamics of this specific CPG system, by finding a way to systematically combine a large number of separate recordings from baseline and perturbed preparations, thus allowing knowledge discovery (e.g. the types of transitions; the changes in the regular state that precede transitions) on those combined data.
We agree with the reviewers’ assessment of our paper. This is indeed the purpose of this paper, and we hope it is clearer to the reader now that we have focused our attention on this.
tSNE and validating the clustering: the authors did considerable work related to this issue, by using their approach on synthetic data, replacing tSNE with UMAP, and by testing simple hierarchical and k-means clustering of their feature vectors. But none of these spoke directly to the validity of the clustering: the synthetic data and UMAP were still manually clustered, so could be made to conform to the desired results;
We apologize for not being clear – we did not manually re-cluster either the synthetic data or the UMAP embedding. We used the labelling generated in the original t-SNE map to color all other maps (t-SNE embeddings with different parameters, initializations, UMAP embeddings). The fact that these other maps appear “clustered” even when we use the original labelling offers (qualitative) evidence that the heavy lifting in the labelling isn’t being done by the human but by the dimensionality reduction technique we used. Similarly, we did not re-color or label clusters in the synthetic data embedding (Figure 2—figure supplement 5). The fact that clusters are almost entirely of a single color suggests that the first part of our method (feature reduction + t-SNE) is powerful enough that separating clusters in the synthetic data is largely trivial by a human.
The unsupervised clustering approach used shows the well-known problem that directly clustering a high-dimensional dataset (a 48-dimensional feature vector) is often impossible using classical techniques simply due to the curse of dimensionality. The normal approach here would be to use dimension reduction then unsupervised clustering (e.g. in the classic combination of PCA then k-means), matching the way the two steps are done in the present paper.
We thank the reviewer for this suggestion. We have added a new figure (Figure 2—figure supplement 2) showing the first two principal components of the feature vectors. We applied k-means to the first ten principal vectors and colored the dots in both the projection of the first two principal components and the t-SNE embedding using this k-means clustering, demonstrating the point the reviewer makes.
As the paper shows, the clusters found in tSNE by manual curation are arbitrary. For example, changing the perplexity parameter changes the fragmentation of the tSNE map, so would likely lead to different manual clusterings.
Figure 2—figure supplement 4 shows the opposite. We use the same coloring for all perplexities, and we see that as long as the perplexity is sufficiently high, we get qualitatively the same result. If one were to manually re-cluster embeddings at higher perplexities, the result would be similar.
In another example, the clusters found by manual curation are not all contiguous in the 2D space e.g. the LP-silent cluster is in 3 groups, one of which is on the opposite side of the 2D space to the others – so distance in this space is not interpreted as meaningful for all clusters.
This is correct. We previously reasoned that most of the distance measurements were between pairs of points within the same cluster, but we agree with the reviewer that we did not rigorously show this. We have removed all distance-based metrics and plots from the paper, because it is confounded by this effect.
In another example, the dividing line between some "clusters" is arbitrarily put in a contiguous run of points – e.g. the LP-weak-skipped (pink) and irregular-bursting (brown).
This is correct. This stems from a smooth variation in observed spike patterns between one canonical form to another, so there truly is no discontinuity when the raw spike patterns are observed. This can be seen clearly in our interactive visualization of the data (a preview of which is available at https://srinivas.gs/stg). When putting smoothly varying data into categorical bins, this is unavoidable. A strength of our visualization method is that it is clear when there exist clear distinctions between categories (e.g., between PD-weak-skipped and LP-weak-skipped) and when there does not (e.g., between the LP-weak-skipped and irregular bursting, as the reviewer points out).
Reviewer #3:
One remaining edit that would be helpful is to show the projection of the data onto the first two PCs (re: the PCA analysis referenced in Figure 2-Figure supp 2), so that the reader can visualize what structure there is prior to computing the triadic difference statistics.
We show the first two principal components in the new Figure 2 —figure supplement 2.
If not too computationally intensive, it would also be useful to show what emerges from applying a naive clustering or visualization algorithm directly on the spike-train time-series [see (2) in Reviewer #3 major concern], even if only to show that it is quite messy.
There is no straightforward way to do this. If we treat the spike trains as binary series, small shifts in the time of first spike can lead to spurious large differences between similar or identical spike trains. This was our motivation for converting spike trains to ISI sets. Even then, because ISI sets are of arbitrary length, it is not trivial to directly apply a naïve clustering or visualization algorithm.
https://doi.org/10.7554/eLife.76579.sa2Article and author information
Author details
Srinivas Gorur-Shandilya
Elizabeth M Cronin
Funding
National Institutes of Health (T32 NS007292)
- Srinivas Gorur-Shandilya
National Institutes of Health (R35 NS097343)
- Srinivas Gorur-Shandilya
- Eve Marder
National Institutes of Health (MH060605)
- Dirk M Bucher
- Farzan Nadim
Deutsche Forschungsgemeinschaft (DFG SCHN 1594/1-1)
- Anna C Schneider
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This article includes data collected by Lamont Tang, Lily He, Mara Rue, Jessica Haley, Daniel Powell, Anatoly Rinberg, and Ekaterina Morozova. We gratefully acknowledge helpful conversations with Paul Miller, Mark Zielinksi, Sriram Sampath, and Alec Hoyland. This work was funded by NIH grants T32 NS007292 (SGS) and R35 NS097343 (EM and SGS), NIH MH060605 (FN and DB), and DFG SCHN 1594/1-1 (ACS).
Senior and Reviewing Editor
- Ronald L Calabrese, Emory University, United States
Reviewer
- Alex Williams, Stanford University, United States
Version history
- Preprint posted: July 7, 2021 (view preprint)
- Received: December 21, 2021
- Accepted: March 18, 2022
- Accepted Manuscript published: March 18, 2022 (version 1)
- Version of Record published: April 11, 2022 (version 2)
Copyright
© 2022, Gorur-Shandilya et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 1,239
- Page views
-
- 227
- Downloads
-
- 5
- Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
A two-part list of links to download the article, or parts of the article, in various formats.
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Mapping circuit dynamics during function and dysfunction
eLife 11:e76579.
https://doi.org/10.7554/eLife.76579
Further reading
-
- Computational and Systems Biology
A deep analysis of multiple genomic datasets reveals which genetic pathways associated with atherosclerosis and coronary artery disease are shared between mice and humans.
-
- Computational and Systems Biology
Mouse models have been used extensively to study human coronary artery disease (CAD) or atherosclerosis and to test therapeutic targets. However, whether mouse and human share similar genetic factors and pathogenic mechanisms of atherosclerosis has not been thoroughly investigated in a data-driven manner. We conducted a cross-species comparison study to better understand atherosclerosis pathogenesis between species by leveraging multiomics data. Specifically, we compared genetically driven and thus CAD-causal gene networks and pathways, by using human GWAS of CAD from the CARDIoGRAMplusC4D consortium and mouse GWAS of atherosclerosis from the Hybrid Mouse Diversity Panel (HMDP) followed by integration with functional multiomics human (STARNET and GTEx) and mouse (HMDP) databases. We found that mouse and human shared >75% of CAD causal pathways. Based on network topology, we then predicted key regulatory genes for both the shared pathways and species-specific pathways, which were further validated through the use of single cell data and the latest CAD GWAS. In sum, our results should serve as a much-needed guidance for which human CAD-causal pathways can or cannot be further evaluated for novel CAD therapies using mouse models.
