When spinal circuits generate rhythmic movements it is important that the neuronal activity remains within stable bounds to avoid saturation and to preserve responsiveness. Here, we simultaneously record from hundreds of neurons in lumbar spinal circuits of turtles and establish the neuronal fraction that operates within either a ‘mean-driven’ or a ‘fluctuation–driven’ regime. Fluctuation-driven neurons have a ‘supralinear’ input-output curve, which enhances sensitivity, whereas the mean-driven regime reduces sensitivity. We find a rich diversity of firing rates across the neuronal population as reflected in a lognormal distribution and demonstrate that half of the neurons spend at least 50 of the time in the ‘fluctuation–driven’ regime regardless of behavior. Because of the disparity in input–output properties for these two regimes, this fraction may reflect a fine trade–off between stability and sensitivity in order to maintain flexibility across behaviors.https://doi.org/10.7554/eLife.18805.001
Where and how are rhythmic movements, such as walking, produced? Many neurons, primarily in the spinal cord, are responsible for the movements, but it is not known how the activity is distributed across this group of cells and what type of activity the neurons use. Some neurons produce regular patterns of “spiking” activity, while others produce spikes at more irregular intervals. These two types of activity have different origins and represent different states of the neural network. It is not clear whether they participate equally in a movement, or if there is a hierarchy among the neurons, such that some neurons have more influence than others.
Petersen and Berg studied neurons in the lower spines of turtles during rhythmic movements. The experiments show that during rhythmic scratching some neurons are very active while most aren’t particularly active at all. This is known as a lognormal distribution and is seen in many other situations, such as the levels of income of people in a society.
Petersen and Berg also found that neurons can move between two regimes of activity, called the mean-driven and fluctuation-driven spiking regimes. During rhythmic scratching, the neurons are almost equally divided between the two regimes, and this division is also found in other types of rhythmic movement. This even division between the two regimes is likely to be important for maintaining a balance between the sensitivity and stability of the neural network. The next steps following on from this work are to reveal the mechanisms behind the two regimes and to find out what causes these differences in activity.https://doi.org/10.7554/eLife.18805.002
Rhythmic movements, such as walking, scratching, chewing and breathing, consist of a recurrent sequence of activity, which is generated by neuronal networks primarily in the spinal cord and medulla. Although, this sequential activity is formed by collective communication among the neurons, it is unknown how the participation is shared versus divided within the population. Distinct motor tasks have been reported to be divided among dedicated microcircuits in zebrafish (Ampatzis et al., 2014; Bagnall and McLean, 2014; Fetcho and McLean, 2010). Nevertheless, do all neurons, which are dedicated to a particular motor activity, spike at approximately the same rate? Or do only some neurons spike at high rate, while most others spike at lower rates? An arrangement with a spectrum of different firing rates could be beneficial by adding the possibility of increasing the overall activity, for instance during uphill walking where a stronger force is needed. In this way the spinal circuit could enhance flexibility by adopting a diversity of firing rates across the population. Other networks in the central nervous system face a similar challenge of how to distribute the activity across the population in order to collectively increase the dynamic range (Wohrer et al., 2013). In sensory processing, neural circuits must be able to retain sensitivity both to weak and strong input. Weak stimuli are amplified whereas strong stimuli are attenuated in order to reduce saturation. If there is too much activity, the circuit reaches saturation and therefore loses the ability to resolve differences in sensory input. Furthermore, amplification of weak signals by recurrent excitation pose the risk of unstable activity, which can spin out of control (Vogels et al., 2005). This computational challenge of how networks maintain both stability and sensitivity is an open question especially for spinal networks.
Stability has primarily been investigated in cortical networks and much evidence suggest that local excitation is carefully balanced by inhibition to assure stability and to widen the range of operation (Galarreta and Hestrin, 1998; Shu et al., 2003). It is well–established that unstable states such as epileptiform activity can easily be achieved by shifting the balance in favor of excitation, e.g. by blocking inhibition (Dichter and Ayala, 1987; Bazhenov et al., 2008). The concept of balanced excitation (E) and inhibition (I) (balanced networks in short) was introduced two decades ago (Shadlen and Newsome, 1994; van Vreeswijk and Sompolinsky, 1996) and has sparked numerous studies both theoretical (Amit and Brunel, 1997; Ozeki et al., 2009; van Vreeswijk and Sompolinsky, 1998; Kumar et al., 2008) as well as experimental (Berg et al., 2007; Okun and Lampl, 2008; Higley and Contreras, 2006; Wehr and Zador, 2003; Kishore et al., 2014). The primary purpose of theoretical models of balanced networks was initially to understand irregular spiking, which was widely observed in experiments (Bell et al., 1995; Shadlen and Newsome, 1994). Irregular spiking was puzzling because it could not be explained by random arrival of excitatory input alone, since this randomness was effectively regularized by temporal integration (Denève and Machens, 2016; Softky and Koch, 1993). Models of balanced networks not only were able to explain irregular spiking, but also revealed other interesting phenomena, such as emergent linearity (van Vreeswijk and Sompolinsky, 1996), multifunctionalism (Sussillo and Abbott, 2009; Hennequin et al., 2014) and self–sustained stable network activity (Amit and Brunel, 1997; Hansel and Mato, 2001; Ikegaya et al., 2013).
The consensus view thus became that irregular spiking results from a mean membrane potential, which is lurking just below threshold, where it is restrained by inhibition concurrent with excitation (Shadlen and Newsome, 1998; Bell et al., 1995; Salinas and Sejnowski, 2000), although synchrony of random excitation is sometimes needed when individual synaptic potentials are small (Stevens and Zador, 1998). This view was essentially predicted much earlier in random walk models (Gerstein and Mandelbrot, 1964). The concept of balanced E/I is now an integrated part of understanding network processing in cortex and elsewhere, but for some reason it has been forgotten in understanding spinal motor networks, with the exception of a few isolated studies (Berg et al., 2007; Petersen et al., 2014).
The balanced E/I allow a subthreshold fluctuating membrane potential, where the spikes are evoked by synaptic transients and therefore belong to the fluctuation–driven regime (Kuhn et al., 2004; Tiesinga et al., 2000). This is in contrast to the more traditional mean–driven spiking (Figure 1), where the mean membrane potential () is well above threshold and spike timing is controlled by after–hyperpolarization (Gerstner et al., 2014; Renart et al., 2007). These two regimes have contrasting manifestations (Table 1): The fluctuation–driven regime has a skewed/lognormal firing rate distribution whereas the mean–driven regime has regular spiking and a symmetric distribution. A simple mechanism has been proposed to explain the lognormal firing in the fluctuation–driven regime by Roxin et al. (2011): The skewness in distribution arises out of a supralinear transformation of the synaptic input, which is Gaussian by virtue of the central limit theorem (Figure 1A). A response to multiple input, which is larger than the sum of their individual responses (i.e. supralinear), will enhance sensitivity (Rubin et al., 2015) and therefore this mechanism may constitute an important physiological purpose.
This is in contrast to the mean–driven regime where the summation is linear or even sublinear, which will transform a normally distributed input to a normally (as opposed to lognormally) distributed firing rate (Figure 1B). Such linear (or sublinear) transformation will reduce rather than enhance sensitivity and therefore the mean–driven regime will curb network activity (Ahmadian et al., 2013). These two transformations work together into an S-shaped IO-curve, where weak input are amplified yet the network is kept stable for strong activity (Figure 1C). Sample neurons in the two regimes are shown (Figure 1D–E). If this mechanism is true, then the shape of the firing rate distribution will reveal the spiking regime of a given neuron. The degree to which neurons operate in one versus the other regime may hold the key to understanding stability, dynamic range and other important properties of network operations. Yet this still remains to be investigated, especially in spinal networks.
Here, we investigate the regimes of operation of spinal neurons during different rhythmic motor behaviors, which are generated in the lumbar spinal circuits of turtles. We test the theoretical scheme put forward by Roxin et al. (2011), by assessing the synaptic input, the spike response function in subthreshold domain, and determine the shape of the firing rate distribution. The mechanical stability of the turtle preparation allows electrophysiological recordings of unprecedented quality, such that we can combine intracellular recording with multi–electrode arrays, and thus determine the fraction of the population in the two regimes at all times. The high resistance to anoxia of turtles allows using adult animals with fully developed spinal circuitry, which have healthy network activity and which can perform multiple complex motor behaviors (Stein, 2005). Thus, we can investigate the population activity during, not just one behavior, but multiple motor behaviors. Custom designed high–density silicon electrodes recorded the population activity from hundreds of cells in the dorsoventral and rostrocaudal axes along with the intracellular of single neurons and multiple relevant motor nerves (Figure 2). This is a unique experimental investigation, because it explores the link between neuronal ensemble data, which in itself is rare in spinal motor research, and the forefront of theoretical neuroscience.
The parallel spiking activity of 200–300 single units were recorded in the medial to ventral horns of lumbar spinal segments involved in motor rhythm generation (Figure 2A). The location of the electrode arrays in the ventral area of the lumbar enlargement was verified by histology (Figure 2B–C and Figure 2—figure supplement 1). The array recordings were performed simultaneously with recording of the intracellular activity of a single neuron in parallel with electroneurograms (ENGs) from relevant motor nerves (Figure 2D). Site–specific rhythmic hindlimb scratching was induced by tactile touch of the carapace (Berkowitz et al., 2010; Stein, 2005) and could be reproduced reliably over multiple trials (Petersen et al., 2014; Vestergaard and Berg, 2015). The extracellular multielectrode arrays, which were used, were custom–designed for the spinal cord (Berg64-probe, Neuronexus inc.) to enable efficient polytrode spike sorting (Figure 2E and Figure 2—figure supplement 2). The distribution of spike count firing rates across the population was skewned (Figure 2F), but resembled a normal distribution on logarithmic x-axis (inset), i.e. a lognormal distribution. This lognormal distribution indicates a wide degree of participation in the motor activity across the population. In the following, we will investigate the participation of neurons within the mean– and fluctuation–driven regimes and how this is linked to the lognormal firing rate distribution, both across the population and within individual cells. We start by addressing the mechanism behind the lognormal firing rate distribution in intracellular recorded data, before addressing the concurrent population activity.
Two mechanisms have previously been proposed to explain the skewned lognormal firing rate distribution, which is also observed in other parts of the nervous system (Buzsáki and Mizuseki, 2014). Lognormal distributions could either arise from a nonlinear transformation of normally distributed inputs (Roxin et al., 2011) (Figure 1A) or from a linear transformation of a lognormally distributed synaptic input (Wohrer et al., 2013). The latter mechanism was considered in connection with the sparse spiking activity in auditory cortex (Koulakov et al., 2009; Hromádka et al., 2008) and since synaptic weights within neocortex have a heavy tail lognormal distribution rather than a Gaussian distribution (Ikegaya et al., 2013; Song et al., 2005). Models also show that the distribution can be either skewed or Gaussian depending on the synaptic input intensity (Ostojic, 2011). Therefore, to distinguish between the proposed mechanisms, it is important to first assess whether the synaptic current is normally versus lognormally distributed. Secondly, to test whether the transformation of the synaptic input to spiking output is linear versus supralinear. We started by addressing the first requirement by investigating the synaptic input in intracellular recordings. The most relevant part of the data was found during the peak of a locomotor cycle where the was in vicinity of and was dominated by synaptic potentials (Figures 1D and 3A). The motor activity was clearly non–stationary, which means that the spike activity was likely to move between the fluctuation– and mean–regime. Nevertheless, the rhythmic activity possessed a separation of timescales in the sense that the activity between cycles (1 s) contained much larger excursions in than within cycles (2-400 ms). Here, the mean did not change much and for practical purposes it could be considered constant within the cycle. In the following analysis of the intracellular data we regarded the dynamics in as stationary within a cycle – well aware that the comparison to theoretical models, which are based on assumption of stationarity, should be taken with a grain of salt. We intended to investigate the symmetry of the distribution of synaptic current using this assumption. The synaptic current within a cycle is difficult to assess, but rather than the mean current, we were primarily interested in the fluctuations in current, which we could approximate from via Ohm’s law under the following conditions. Within a cycle, the mean was just below threshold and did not change its value much. Therefore the voltage–activated conductances were approximately constant such that there was an Ohmic relationship between synaptic current and . This is likely justified for neurons in fluctuation–driven regime, since the conductance is often high and dominated by balanced E/I synaptic input (Destexhe et al., 2003; Kumar et al., 2008). The high conductance suppresses the coupling between and intrinsic conductance in a divisive manner (Kolind et al., 2012; Tiesinga et al., 2000). Thus, in the fluctuation–driven regime the non–Ohmic contributions were likely smaller and the -relationship more linear than in the mean–driven regime.
We intended to test the hypothesis of normally distributed input, but since the approximation of using the variability in as a proxy for the variability in synaptic current is most valid for the neurons in fluctuation–driven regime, we needed a way to distinguish neurons that were primarily in the fluctuation–driven regime. We therefore propose a novel metric, the return map ratio , which quantifies the degree of fluctuations leading up to a spike (Figure 3—figure supplement 1). The return map ratio (RMR) quantifies how direct the subthreshold –trajectory is between spikes and this forms the basis for selecting neurons in our analysis. An RMR close to 0.5 has fluctuation–driven spiking whereas a value close to 1 has mean–driven spiking (Figure 3—figure supplement 1A,B). Therefore, we defined a neuron as fluctuation-driven if its RMR ; in our sample of intracellular recordings we found 50/68 neurons in this regime. A sample neuron, which was found in the fluctuation–driven regime based on this metric illustrates how we obtained the distribution of sub–threshold (Figure 3A). The distribution was estimated both by selecting the in between spikes (temporal distribution) and by collecting instances of prior to spike peak in a spike triggered overlay (‘sigma’ in Figure 3B). These two estimates are in agreement with one another for the sample cell (Figure 3C). This agreement is also found across the population as quantified by the mean and SD (Figure 3D). The skewness for the distributions across the population is small and scattered around zero as expected for normal (symmetric) distributions (Figure 3E). From these data we conclude that the subthreshold –distributions are not skewed, but rather symmetrical and Gaussian–like (cf. inset distributions, Figure 3E). Nevertheless, the minimal requirement for confirming the two–regime hypothesis for the single neuron is that the synaptic current (not the synaptic potentials) is Gaussian (Figure 1). As we argued earlier, if there is an Ohmic relationship between current and potential, which is likely during high–conductance states, then this requirement would be granted. More importantly, now that we do find a Gaussian –distribution, it is difficult to contemplate a non-linear -relationship, which would result in such a symmetric distribution. The synaptic input current would have to have a finely matched inverse distribution to cancel out this non–linearity in order to achieve a symmetric –distribution. A more parsimonious explanation therefore is that, since the synaptic potentials are normally distributed, they are a result of a linear transformation of synaptic currents, which are also normally distributed.
So far, we have only looked at –distributions of single neurons, which operate primarily in the subthreshold domain, and found that the synaptic input is most likely normally distributed. We do not know whether the synaptic input is also normally distributed in the mean–driven regime, but since the synaptic input is normally–distributed in the subthreshold region, it is likely also normally–distributed in the suprathreshold region. Otherwise, the input statistics from the presynaptic neurons would have to depend on the threshold of the post–synaptic neuron, which is unlikely.
Above, we established that the synaptic input to a given neuron is likely normally distributed, and if this input is transformed in a supralinear fashion, the output firing rate distribution will be skewed. Nevertheless, the foundation of the skewness in population rate distribution (Figure 2F) is not necessarily directly linked to the skewness of the instantaneous rate distribution of single neurons. In principle, it is possible to have a population with a normal distribution of mean firing rates, where the cells themselves have lognormally distributed firing rates and vice versa. Therefore, we needed to address the distribution of mean across the population and test whether this was skewed or normal. Further, since the sub–threshold IO-curve is linked to threshold, it is important to establish the distance of mean from threshold with respect to the size of synaptic fluctuations, i.e. standard deviation of (). This distribution, i.e. , turns out to also be normally distributed with a mean around 3 from threshold (Figure 3—figure supplement 2, plotted for all neurons). The value used for here is the mean of the estimated thresholds for all spikes (see below). If we assume, when normalizing this way, the IO-curve has approximately the same nonlinearity across all neurons, the population distribution of firing rates will also be skewed due to the nonlinear transformation of the normally–distributed input (Figure 3—figure supplement 2F) to a lognormally–distributed output. These results are in qualitative accordance with the scheme proposed previously (Roxin et al., 2011). As another piece of the puzzle, we need to establish the shape of the neuronal response function, which rarely has been done in the subthreshold domain.
The link between a normally distributed input and a lognormally distributed output is a supralinear transformation. To test whether this is a hallmark of the fluctuation–driven regime, we needed to estimate the input–output (IO)–function for the subthreshold domain. The IO–function of neurons is a fundamental property of the nervous system, and therefore it is well-characterized both theoretically (Gerstner et al., 2014) and experimentally (Silver, 2010). Nevertheless, it has rarely been established for fluctuation–driven spiking. Here, we estimated the IO-function for subthreshold spiking via the probability of eliciting a spike as a function of in the following way. First, we collected instances of shortly before the spike–onset, where is depolarized yet still not part of the deterministic spike trajectory. The probability that a given value of will cause a spike was estimated as the histogram of –instances (gray histogram, Figure 4A) divided by the total time spent at all values of (green histogram). This gives the empirical relationship between and the firing rate (Jahn et al., 2011; Vestergaard and Berg, 2015). The IO–function had a strong non–linear shape (Figure 4B). To capture the curvature we fitted both a power–law and an exponential for all neurons and the curvature had a weak negative correlation with the SD of the –fluctuation (Figure 4C–D) as demonstrated previously (Vestergaard and Berg, 2015). Similar expansive nonlinearity has previously been characterized in sensory–driven neurons (Anderson et al., 2000; Hansel and van Vreeswijk, 2002; Miller and Troyer, 2002). It will transform the normally–distributed synaptic potentials into a lognormally–distributed spiking output in the fluctuation-driven regime (Figure 1A). For mean–driven spiking the IO-function is not supralinear, but rather linear (or even sublinear), and the normally–distributed synaptic input will therefore be transformed to a normally distributed spiking output (Figure 1B). In conclusion, neurons that have fluctuation–driven spiking also have a non–linear IO-transformation of synaptic potentials to spiking output.
The normally distributed input combined with the nonlinear IO–transformation should result in a skewed lognormal firing rate in the single neuron. To confirm this, we measured the distribution of the instantaneous firing rate, i.e. the inverse of ISIs. The quiet period in between burst cycles were not included in the analysis (Figure 1D–E), since in these periods was far from and therefore in an irrelevant part of the IO–function. The firing rate distribution of many cells was positively skewed and resembled a normal distribution with near zero skewness on a log-scale (sample cell shown in Figure 5A). This is expected for poisson–like spiking in the fluctuation–driven regime (Ostojic, 2011). Nevertheless, distributions for all the intracellularly recorded neurons () were skewed to a varying degree from strong positive to zero skewness on a linear axis and similarly shifted downwards on log axis (cf. gray and green histograms, Figure 5B). This suggests that neurons were found in a spectrum between fluctuation– and mean–driven spiking. More negative log–skewness were associated with higher mean rates (Figure 5C). This is probably due to a larger presence in the mean–regime at higher firing rates, where the distribution skewness is expected to be negative on a log–scale, i.e. Gaussian on a linear scale. Note that the spectrum of skewness was substantially larger than it was for the distributions above (Figure 3E). Skewed Gaussian distributions are shown to illustrate the range of skewness in the data (Figure 5D). In conclusion, these results suggest that the skewness in firing rates is an indicator of the degree of participation in the fluctuation–driven regime.
A neuron is not just spiking in either the fluctuation– or the mean–driven regime, rather, it likely spends time in both regimes during motor activity. To estimate the amount of time a given neuron spends in either of the two regimes we calculated the fraction of time that the smoothed was above versus below threshold. We first look at two heuristic neurons, one in the fluctuation–driven regime and one in the mean–driven regime. The fluctuation–driven neuron spent most of the time below threshold (Figure 6A) and had more irregular spiking as quantified by a local measure of irregularity, the (green line). is the difference of two adjacent ISIs divided by their mean (Holt et al., 1996; Bruno et al., 2015). In contrast, the mean–driven neuron spent most time above threshold and had more regular spiking, i.e. closer to zero (Figure 6B). Since the threshold was firing rate–dependent due to the inactivation of the –conductance (Figure 6—figure supplement 1) we used the most hyperpolarized value of threshold (broken line). The distribution of for all trials had higher mean for the fluctuation–driven cell than the mean–driven (cf. arrows, Figure 6C). Also, the cumulative time spent below threshold was higher for the fluctuation–driven cell (96%) than the mean–driven cell (35%, Figure 6D). This fraction of time spent below threshold was quantified for every neuron () and the population distribution had a strong mode at 1 (top, Figure 6E) suggesting many neurons spent much time in the fluctuation–driven regime. To compress the diversity within the population into a simpler representation, we used the reverse cumulative distribution of neurons versus time spent below threshold (bottom, Figure 6E). This indicates how many neurons (y-axis) spent at least a given fraction of time (x-axis) below threshold. The intercept with the 50%–line (broken line) indicates what fraction of time half the population at least spent below threshold. This fraction is remarkably high (84%) suggesting a prominent presence within the fluctuation–driven regime.
Mean- and fluctuation-driven spiking can be distinguished by important traits such as degree of irregularity and log-skewness of the firing–rate distribution. To verify these traits, we used another sample neuron as a heuristic illustration. We injected different levels of either positive or negative bias currents in different trials while keeping all else constant. A negative constant current injection (−1.0 nA) caused a decrease in firing rate and a slight increase in irregularity (green line) compared with zero injected current (Figure 7A–B). Similarly, a positive current injection (1.7 nA) caused more spikes and a decrease in irregularity (Figure 7C) consistent with a movement between regimes (inset in Figure 7A). The decrease in irregularity with increasing input was further quantified as a negative correlation between mean and injected current (, p0.001) over multiple trials (n=18, Figure 7D). This is qualitatively in agreement with previous reports (Prut and Perlmutter, 2003; Powers and Binder, 2000; Wohrer et al., 2013). The instantaneous firing rate in the control condition (0 nA) was lognormal as expected for the fluctuation–driven regime (top, Figure 7E). When adding input current the distribution was shifted to the right and enriched with a negative skewness as expected for mean-driven spiking (bottom, Figure 7E). This relation between input and shape of rate distribution was further confirmed by a negative correlation between multiple current injections and skewness both on linear scale (gray dots) and log–scale (red dots, Figure 7F). Hence, skewness and irregularity are indicators of the spiking regime.
An alternative to injecting electrode current is to manipulate the balance of excitation and inhibition (E/I) by pharmacological means. This is important for understanding the cause of irregularity and the fluctuation–driven regime. Hence, we manipulated the synaptic input in a reduced preparation with micro–superfusion of strychnine, a glycinergic blocker, over the transverse cut surface of the spinal cord (described in [Berg et al., 2007; Vestergaard and Berg, 2015]). This affected only neurons at the surface (300 μm) without affecting the rest of the network, which was verified by careful monitoring of flow and the network activity via the nerve recordings. Comparing the spiking during control condition (Figure 8A) with that during blockade of inhibition (Figure 8B), we noticed a strong increase in spiking. This is consistent with a depolarization due to disinhibition, thus ‘unbalancing’ the excitatory and inhibitory input. Reducing inhibition tipped the balance of E/I toward larger inward synaptic current, which resulted in a more depolarized (blue line) well above threshold (arrows, Figure A–B). It also resulted in higher firing rates and lower irregularity on the peak (cf. green lines). Generally, the irregularity () was higher in the control case than in the unbalanced case (Figure 8C) similar to the results observed with current injection (Figure 7A–D). The irregularity was also negatively correlated with depolarization of the mean when unbalancing the E/I although it was uncorrelated in the control condition, where the spiking occurred in the fluctuation–driven regime (Figure 8—figure supplement 1). The instantaneous firing rate was skewed and lognormal in the control case (top, Figure 8D), similar to the above sample cell (top, Figure 7E). This distribution became negatively skewed when adding inward current (bottom, Figure 7E). Similar effect was seen when ‘unbalancing’ the synaptic input, which also result in larger inward current. The firing rate increased (cf. broken lines, Figure 8D) and the distribution became negatively skewed (cf. −0.2 and −1.5) as expected in the mean–driven regime (bottom). To quantify the increase in time spent in the mean–driven regime, we performed an analysis similar to the analysis in the above section (Figure 6D). The cumulative time spent below threshold was larger in the control condition (78%) compared with the unbalanced case (56%, Figure 8E). These observations are largely consistent with the consensus view that irregular fluctuation–driven spiking is due to a balance between excitation and inhibition (Table 1).
In the above intracellular analyses we reported the spiking irregularity in terms of along with the mean , current injection and pharmacological manipulation of the balance of excitation and inhibition. The measure is convenient to use as an indicator of the mean– versus the fluctuation–driven regimes observed in the extracellular spiking data, since it only requires spike times. Therefore it is important to validate as an indicator of spiking regime. In the above sample cell analyses we note first, that when spent a larger fraction of time above threshold, i.e. in mean–driven regime, the was lower (Figure 6). Second, when depolarizing a neuron artificially either with constant positive current (Figure 7D), or by blocking inhibition (Figure 8C), such that more spikes were in mean–driven regime, the was decreased.
To further substantiate as an indicator of spiking regimes we looked again at the return map ratio, which is an independent metric of fluctuations during inter-spike intervals. If is an indicator of the spiking regime, it should be anti-correlated with the return map ratio. This was confirmed by plotting the mean for all cells () against the mean return map ratio, which indeed demonstrated a significant anti–correlation (, p=0.005) (Figure 3—figure supplement 1E).
A second independent indicator of fluctuation regime is the cumulative time below threshold of (Figure 6D), which should be correlated with the mean . We tested this using the most hyperpolarized value of theshold, since it was the most conservative, but there was no significant correlation between the cumulative time below threshold and the mean . Perhaps the lack of linear relationship is due to a bias from the reset voltage and after-hyperpolarization, which is different from cell to cell and therefore randomly may introduce a large fraction of time spent below threshold. Also, intense synaptic activity is known to quench the after–hyperpolarization (Berg et al., 2008) and therefore this bias may be particularly strong when the synaptic input is not balanced as in the mean–driven regime.
A third indicator of spiking regime is the skewness of the instantaneous firing rate distribution (Figure 7E and 8D). We estimated the skewness of the individual firing rate distributions for all neurons () and plotted it against the mean (data not shown). There was a significant positive correlation between the two, regardless of whether the firing rate distribution was plotted on log or linear scale (, p=0.0003, and , p=0.0006), which suggest as a valid measure for spiking regimes.
A last indicator is the local mean membrane potential depolarization, which should be anti-correlated with the instantaneous , if the is above threshold (Figure 8, Figure 8—figure supplement 1D). Here, there was a lack of correlation between and before blocking inhibition, in the fluctuation–driven regime. However, after removal of inhibition, was in supra–threshold domain, which introduced an anti-correlation between and . Hence, if the neuron is in the mean-driven regime the is an indicator for the depolarization above threshold. To further verify this we performed a similar test of the relationship between instantaneous and local depolarization for all neurons (without pharmacology). We found that all the cells with significant relationships (p<0.05, ) had anti-correlation between and (data not shown). In conclusion, the measure is correlated with other measures and indicators of spiking regimes (except the cumulative time below threshold) and therefore is a useful indicator in itself.
The irregularity in spiking could be caused by a noisy threshold rather than fluctuations in synaptic potentials. Nevertheless, a noisy threshold can only explain a small part (if any) of the spiking irregularity. First of all, if the irregularity, that we observed in spike times, was due to a noisy threshold mechanism, we should see the same irregularity regardless of the depolarization, i.e. regardless of whether the neuron was in the sub–threshold or supra–threshold domain. Yet, the spiking irregularity was strongly dependent on depolarization (Figures 6–8). There was an adaptation in threshold (Figure 6—figure supplement 1). This was not random, but rather due to a gradual inactivation of Na –channels throughout the burst (Henze and Buzsáki, 2001). The threshold of a given spike strongly depended on the threshold of the previous spike (panel F) as well as the mean firing rate (panel G). The same mechanism is behind spike–frequency adaptation, which is a well–described phenomenon (Grigonis et al., 2016). The adaptation in threshold is likely to make the IO-function more sublinear in the mean–driven regime, which will generally curb network activity.
In order to verify the extent of the threshold variance beyond the contribution from inactivation of Na+–channels, we looked at the threshold of only the first spike of each cycle, such that the neuron had ample time for recovery. The variance of the first–spike threshold () in a sample neuron was mV2 whereas the variance in synaptic potentials was more than 17–fold higher ( mV2). Therefore a randomness in the threshold had little of no effect on the irregularity of spiking compared with the randomness in synaptic input. In some recordings the threshold may appear as uncorrelated with the membrane potential prior to the spike onset. However, rather than a noisy threshold this is likely attributed to cellular morphology. If the cell is not electrically compact, the axon initial segment, where the spike is initiated, will have a different potential than what is recorded with the electrode. If this was the case, these observations would still be compatible with the two–regime hypothesis, since spikes would still be driven either by fluctuations or a large mean current, despite the disguise of a long electrotonic distance to the recording site.
So far the analysis has been performed on serially acquired intracellular recordings across trials and animals. This demonstrates that some neurons spiked primarily in the fluctuation–driven regime while others spiked in the mean–driven regime. Nevertheless, it is still unclear what the parallel population activity was during a behavior and across behaviors. How many neurons were in one versus the other regime and for how long? First, we assessed the neuronal participation in the motor patterns by their degree of spiking during motor behavior. Neurons were active during both ipsi– and contralateral scratching behaviors (Figure 9A–D). Most units had a rhythmic relationship with the nerve signals and a higher firing rate for the ipsilateral scratching compared with contralateral scratching behavior (cf. Figure 9C and D; Videos 1 and 2), which indicates participation of neurons in a hemicord to a smaller degree in the contralateral movement than the ipsilateral movement.
The distribution of firing rates across the neuronal population over several trials was strongly skewed, which indicate that most neurons spike relatively infrequently with a ‘fat-tail’ of higher spiking (Figure 9E). The distribution covered two orders of magnitudes from 0.1–10 Hz and was akin to a lognormal distribution (inset and green lines, Figure 9E). Similar lognormal–like distributions have been observed in other parts of the nervous system (Buzsáki and Mizuseki, 2014). The implication of the skewed distribution is that most neurons spiked at low rates, but there was relatively many neurons spiking at higher rates indicating an overall rich diversity of firing rates.
Although multi–functional spinal units have been reported previously (Berkowitz et al., 2010) it is unclear how their participation is distributed and whether the asymmetry in distribution is linked to different behaviors. To address this issue we analyzed the population spiking for multiple motor behaviors. The induction of a distinct scratch behavior is location–specific (Stein, 2005). Multiple behaviors can be evoked depending on exact location and which side of the body is touched. This allowed us to induce two distinct behaviors: ipsi– and contralateral hindlimb scratching, while recording from the same neuronal ensemble (Videos 1 and 2). These behaviors were reproducible over multiple trials (9 trials). Both behaviors had similar phase relationships between the muscle synergists, although ipsilateral scratching had stronger activity (cf. Figure 9A and B). The firing rate distribution was positively skewed in both behaviors with the similar qualitative shape (Figure 9E–F). This skewness was also found across animals (green bars, Figure 9G, n=5) and close to zero on log–scale, i.e. lognormal (black lower bars). To further quantify the uneven neuronal participation we used Lorenz statistics and the Gini-coefficient (O’Connor et al., 2010; Ikegaya et al., 2013). The Lorenz curve characterizes the share of cumulative participation of individual neurons of the population (Figure 9H). The diagonal corresponds to the case where all neurons have the same firing rate. The deviation from equality is quantified by the Gini–coefficient, which is the fraction of area to the total area (Figure 9H). The higher the coefficient, the more unequal the participation across neurons is. Both scratch behaviors had a Gini–coefficient of 0.5 (Figure 9I). Although the mean firing rate could change between behaviors and between animals (Figure 9J), the skewness was qualitatively similar (Figure 9K). This suggests that the skewed lognormal–like firing rate distribution, and hence a presence of the fluctuation–driven regime, was preserved across behaviors and animals.
Neurons do not occupy either the fluctuation– or the mean– driven regime all the time. Individual neurons can move back and forth between regimes depending on the synaptic current they receive. Neurons that spike predominantly in the mean–regime will have their mean firing rates closer together and more normally distributed compared with those spiking in the fluctuation–regime. Hence, we expected the skewness of the distribution of mean firing rates across the population to become more negative (on log–scale) as the general network activity increases. To address this, we analyzed the spiking across neurons in parallel. First, we estimated the time–dependent firing rate of each neuron in the population using optimal Gaussian kernel (Shimazaki and Shinomoto, 2010) and measured skewness of the population distribution. The time–dependent population distribution was achieved by binning the rates in 10 ms windows (Videos 1 and 2). The mean population rate and its SD are indicated as black gray lines (Figure 10A). As the mean firing rate increased, the skewness of the distribution (log–scale) became negative, which is a sign of more neurons were occupying the mean–driven regime (cf. inset histograms, Figure 10A). This was further confirmed by a negative correlation between the mean firing rate (black line in A) and the log–skewness for all time points (Figure 10B). Hence, as the general activity increased, the population distribution became less lognormal and more Gaussian, which suggests more neurons occupied the mean–driven regime during a higher general activity.
To further gauge the division of neurons in the two regimes we compared the irregularity of the spiking using . This metric was verified above as a reliable indicator of spiking regimes. The distribution of the mean across the population of neurons was clustered around 1 if all ISIs were included (gray histogram, Figure 10C). However, measuring the irregularity in the motor cycles alone i.e. excluding the inter–burst intervals (here, s) the mean irregularity across neurons was lower and clustered around 0.6 (red histogram). Both distributions had substantial spread around the mean, which suggests a rich diversity spiking patterns.
To get a compound measure of the behavior of the entire population across time, we considered the amount of time each neuron spent in the fluctuation–driven regime. We demarcated the fluctuation–regime as having irregularity in spiking above a critical value, i.e. . Choosing is not entirely objective. Complete Poisson–type irregularity has , but the spiking is still irregular for lower values (Feng and Brown, 1999). Based on our data, even when the 0.5, the spent as much as 96% of the time below threshold (Figure 6C–D) indicating fluctuation–driven spiking. Further, neurons that had 0.5, also had lognormal firing rate distributions (Figure 7), which also indicates the fluctuation–driven regime. For these reasons, we suggest choosing for distinguishing regular vs. irregular spiking. A similar value was previously chosen to distinguish between regular vs. irregular ‘choppers’ in the cochlear nucleus (Young et al., 1988). Thus, the population of spinal neurons had a large diversity in time spent in the fluctuation–driven regime. Some neurons spent as little as 20% in the fluctuation–driven regime while other spent as much as 80%. To get a quantitative handle on the occupation of neurons in the fluctuation–driven regime across the population, we considered the distribution of time spent with . This was formally quantified using the reverse cumulative distribution of neurons that spend a given fraction of time in the fluctuation–driven regime (Figure 10D). The reverse cumulative distribution is plotted for 3 values of (0.4, 0.5 and 0.6) to indicate the sensitivity to parameter choice. Obviously, choosing a lower results in a larger fraction of time in the fluctuation–driven regime, i.e. the curve is shifted to the right. Choosing larger has the opposite effect. This inverted –shaped curve gives the fraction of neurons (y–axis), which spend at least a given time in the fluctuation–driven regime normalized to 100% (x–axis). Hence, half of the population spent at least 58% of time in the fluctuation regime during ipsilateral scratching (intercept of curve with the broken line, Figure 10D). We refer to this metric as the time in the fluctuation–regime for half the neurons (). Similar –values were obtained for all five animals (inset histogram). Qualitatively similar results were achieved for a different motor behavior, namely contralateral scratching (Figure 10E). The metric is a time–weighted analysis of irregularity of spike trains. In addition to measuring the time in regimes, we measured how many spikes were in one regime vs. the other. Hence, we calculated the reverse cumulative distribution of neurons that had a given fraction of spikes in the fluctuation–driven regime (Figure 10—figure supplement 1). Similar to , we defined a spike–weighted metric as the spikes in fluctuation regime for half the neurons (). Both the – and –values were relatively conserved across animals as well as behaviors (Figure 10D–E, Figure 10—figure supplement 1). The large values of and indicate that the fluctuation–driven regime had a strong presence during motor behaviors, and the high similarity suggests that it may represent a conserved fundamental property of network activity.
In the data analyses presented so far we have not addressed the neuronal identity of the recorded units. Nevertheless, there is a spatial division subtypes of spinal neurons, which we could take advantage of. During development, a distinct laminar organization of different cellular subtypes is formed in the dorsoventral axis (Arber, 2012; Jessell, 2000). In particular, motoneurons are primarily found in the most ventral part of the horn whereas interneurons are found in more medial to dorsal territory. Since this is the same axis that our electrode arrays were located along, it was possible to infer a likelihood of cellular identity based on location. The electrode shanks have multiple distributed electrodes (Figure 11A), which made it possible to approximate the soma location using trilateration. Trilateration is the geometrical process of determining the location of a source in space using multiple recording sites combined with the fact that signals decay in the extracellular space (Manolakis, 1996). Thus, the node locations were approximated based on the amplitude of spike waveforms, which clearly decayed with distance (Figure 11B). Node locations were combined for all shanks, probes and animals to form a scatter (Figure 11C). Combining these locations with depth of individual shanks with respect to the surface of the spinal cord, we were able to investigate the spike patterns with respect to the absolute neuronal location. The irregularity in spiking was quantified (mean ) with respect to dorsoventral depth (Figure 11D). The distributions of mean firing rates (not shown) and the mean (Figure 11E) had no obvious dependence on depth. In particular, the spread in means was much smaller than the SD of the distributions themselves. The most parsimonious interpretation of these data is that the fluctuation–driven spiking regime was both present and equally prominent in all the neurons, regardless of whether the cell body was in the ventral horn or in the medial horn, i.e. equally present in motoneurons and interneurons.
In neuronal networks, spikes are generated in either the mean– or the fluctuation–driven regime (Brunel, 2000; Gerstner et al., 2014; Kuhn et al., 2004; Tiesinga et al., 2000). In this report we present evidence for the existence of both regimes during motor pattern generation in the spinal cord. We consistently found normally distributed synaptic input combined with the supralinear shape of the IO–function in the subthreshold region, and suggest this as a compelling mechanism behind the lognormal population firing rate distribution (Roxin et al., 2011). Using spiking irregularity across the neuronal population as a hallmark of the fluctuation regime, we found that half of the neurons spent at least 50% of the time in this regime. Thus, the fluctuation–regime was not a rarity, but rather had a prominent presence both across behaviors and across animals (Figure 10). To our knowledge this is the first report, which quantifies occupation within spiking regimes of a neuronal population, not just in the spinal cord, but also in the nervous system in general.
The fact that the relative time during which a subset of neurons occupied one of the two regimes was conserved across both behaviors and animals could indicate a key principle of neuronal processing. A fundamental challenge for neuronal networks is to perform functions while keeping the population activity from falling into either of the two extreme states: (1) the quiescent state where the neuronal spiking activity cannot remain self–sustained and (2) the unstable state of run–away recurrent spiking activity (Vogels et al., 2005; Kumar et al., 2008). It is well known that recurrent inhibition is important for maintaining stability, but other mechanisms may participate as well, e.g. synaptic depression or active adjustment of the shape of the neuronal response function by adaptation of spiking threshold. A nonlinear response function, as we observed in the fluctuation–regime (Figure 4B), will amplify input via supralinear summation (Rubin et al., 2015). The upward curvature will enhance synaptic fluctuations, which then accelerates the recurrent excitatory activity causing a potentially unstable state. In contrast, the response function in the mean–driven regime, is linear or even sublinear, which is likely to curb strong input. We therefore propose that the close proximity of the –value to 50% is an indication of a self–organizing trade–off between sensitivity and stability in order to preserve at once both network homeostasis and dynamical functionality. This conjecture remains to be further substantiated in future studies. Furthermore, the – and –values remain to be determined for other part of the nervous system and in other species.
The distinction between fluctuation– and mean–driven spiking is interesting because the two types of spiking may have radically different causes, and this may hold an important clue to understanding the enigmatic motor rhythm generation. The fluctuation–driven spiking is believed to be caused by concurrent and random arrival of excitatory and inhibitory potentials resulting in a fluctuating subthreshold (Table 1). In the mean–driven regime, on the other hand, the net membrane current is so large that the mean is above threshold, and the ISIs are therefore determined by the recharging of the membrane capacitance following the refractory period of the previous spike (Powers and Binder, 2000). This results in a deterministic trajectory of and regular ISIs. More importantly, for the mean–driven spiking the membrane current can be caused by intrinsically electrical properties as well as synaptic input, whereas the fluctuation–driven spiking is exclusively caused by synaptic input. An intrinsic property, which is commonly believed to be involved in rhythm–generation, is the pacemaker property that can autonomously generate neuronal bursting in the absence of synaptic input (Brocard et al., 2010; Ramirez et al., 2004; 2011). The prominent presence of the fluctuation–regime therefore implies that the majority of neuronal spikes were not driven primarily by intrinsic properties such as pacemaker potentials, but rather by synaptic communication. This can be interpreted in two ways: (1) if there is a pacemaker–driven rhythmogenic core of oscillatory neurons responsible for the motor activity (Huckstepp et al., 2016), the core only represents a small fraction of the network, or (2) since the fluctuation–regime is prominent and pacemaker neurons are difficult to find, the motor–rhythm may be generated by other means such as emergent collective processes in the network (Yuste, 2015). Generation of movements without the need of pacemaker neurons have been predicted theoretically in central pattern generators (Kleinfeld and Sompolinsky, 1988) as well as more complex sequence generation (Hennequin et al., 2014). Even in the respiratory system, which has the most stereotypic motor rhythm, pacemaker cells appear not to be essential for generation of the rhythmic breathing, although this topic is still debated (Feldman et al., 2013; Ramirez et al., 2011; Carroll and Ramirez, 2013; Chevalier et al., 2016). It remains to be understood how a distributed emergent processes can generate motor rhythms on a network level if, in fact, the pacemaker bursting is not an essential component.
In spinal research, neuronal identification has improved over the last decades with the development of genetic knockouts and molecular markers (Bikoff et al., 2016; Goulding, 2009; Britz et al., 2015; Kiehn, 2006). Pinning down cellular identity improves the search for a potential specialization in the circuit. However, the sole focus on cellular identity to address questions in spinal research carries a weakness as well as a strength. It contains the risk of missing the collective dynamics and the delicate interaction among neuronal cell types. Neural circuits operate to perform functions by collective interaction between all neurons, where it is difficult, if not impossible, to link a particular function to the individual neuron. Functional activity may very well arise on circuit level as opposed to cellular level. This caveat is known as the neuron doctrine versus emergent network phenomena (Yuste, 2015; Grillner, 2006), and the neuron doctrine has almost exclusively been adopted in previous investigations of spinal motor circuits. To the best of our knowledge this report is the first investigation of spinal motor circuits from an ensemble viewpoint.
Nevertheless, since motoneurons are fundamentally different from the rest of spinal neurons it would be helpful to distinguish them from interneurons. In our experiments we sampled from neurons, which were likely to be primarily interneurons since they are more numerous than motoneurons. The fraction of motoneurons to interneurons is 1:8 (Walloe et al., 2011), but we were also likely to sample motorneurons, since they have large somata. To explore this further, we investigated the population activity and its relation to cellular identity by taking advantage of their spatial segregation in the dorsoventral axis (Arber, 2012; Jessell, 2000). We were able to associate an absolute location of the cellular somata (using trilateration), and thus test for differences in spiking activity (Figure 11). The distribution of firing rates as well as the spiking irregularity did not have any dependence on location. This suggests that the fluctuation–driven spiking regime was both present and equally prominent in all the neurons, regardless of whether the cell bodies were in the ventral or medial horn, i.e. regardless of whether they were motoneurons or premotor interneurons.
Common features of network activity for different parts of the central nervous system may provide hints towards fundamental principles of neuronal operations. In the present study we identified the following features of population motor activity: (1) synaptic input to individual spinal neurons was normally distributed (Figure 3), (2) the means of these normal distributions were also normally distributed across the population. In particular, the distance to threshold in terms of fluctuations, i.e. had a normal distribution and a distance from mean to threshold of 3 on average (Figure 3—figure supplement 2F). (3) The neuronal response function was supralinear when the mean input was in the subthreshold region (Figure 4). (4) There was a rich diversity of regular to irregular spiking patterns. (5) The population firing rate was skewed and lognormal–like.
Many of these features have been identified before in other parts of CNS. The of individual neurons is often normally distributed in cortical neurons when considering either the up– or down–state (Destexhe et al., 2003; Stern et al., 1997) and the spiking is irregular with a clustered around 1 (Softky and Koch, 1993; Stevens and Zador, 1998). Similar irregularity is observed in invertebrates (Bruno et al., 2015). The distribution of mean values in our experiments was clustered around 0.6 when ignoring the inter–burst intervals (Figure 10C). This is more regular than what is observed for typical cortical neurons (although see Feng and Brown, 1999), but similar to cervical interneurons in monkeys performing isometric wrist flexion–extensions (Prut and Perlmutter, 2003).
We observed a skewed and lognormal–like population distribution across behaviors (Figure 9, Videos 1 and 2). Similar lognormal distributions have been reported in other parts of CNS (Buzsáki and Mizuseki, 2014; Hromádka et al., 2008; O’Connor et al., 2010; Wohrer et al., 2013) and it remains an open question how the skewness arises out of neuronal ensembles. Roxin et al proposed the mechanism where the skewness arises from a nonlinear transformation of Gaussian input (Roxin et al., 2011). Our data supports this hypothesis. First, we observed a normally distributed for individual cells, which is a proxy for the requirement of normally distributed input currents (Figure 3). Second, a supralinear IO–function covering most of this input (Figure 4). Third, a firing rate distribution of individual cells which was typically highly skewed and lognormal–like although some did not have lognormal firing (Figure 5). Nevertheless, there is a distinction between the lognormal firing of individual neurons and the lognormal distribution of mean rates across the population. Whereas the lognormal population firing rate remains to be fully understood, the skewed firing rate distribution of individual neurons is fairly well understood. Here, the skewness is due to the fluctuating input and irregularity of spiking (Ostojic, 2011). Nevertheless, we argue the mechanism for the lognormal population firing is the same as that for the individual neuron. If the subthreshold IO-function is approximately similar across the population, which our data implies (Figure 4), we can explain the lognormal population firing by a supralinear transformation, if the mean across the population is also Gaussian. We did in fact find the distribution of mean to be Gaussian (Figure 3—figure supplement 2F).
Classical studies of spinal motoneurons indicate two regimes of spiking: a primary and a secondary range (Kernell, 2006; Meehan et al., 2010), which corresponds to different parts of the mean–driven spiking regime. This characterization was associated with the intrinsic properties without synaptic input being present. Nevertheless, a different type of fluctuation–driven spiking was discovered in experiments where synaptic input were present, in what was referred to the subprimary range in mice (Manuel and Heckman, 2011) and humans (Kudina, 1999; Matthews, 1996). This subprimary range conforms to the fluctuation–regime though under a different terminology. As the name indicates, the primary range has been considered to represent the dominant mode of spiking whereas the subprimary range is a peculiarity. Nevertheless, a recent study recorded for the first time the motoneuron discharge and muscle force and found that the subprimary range accounts for 90% of the contraction force (Manuel and Heckman, 2011). This indicates that the fluctuation–regime may have a more noteworthy role in covering the dynamical range in motor control than previously assumed, which is in agreement with the observations of the present study.
We used the integrated turtle preparation with the spinal motor network intact ( for the multi–electrode recordings and for the serially aqquired intracellular recordings), in order to address how the neuronal firing rates are distributed across the population of interneurons and motoneurons in the spinal cord (Petersen et al., 2014). These sample sizes where assumed to be large enough in the experimental design and because of a consistency in results, although a specific power analysis was not conducted. To avoid the confounding factors of supraspinal input, we spinalized the turtle. The transection was performed at the spinal cord at segments (D3-4) caudal to the cervical segments, where the local circuitry has only little or no involvement in generation of motor patterns (Mortin and Stein, 1989; Hao et al., 2014; Mui et al., 2012). The adult turtle preparation is capable of producing elaborate motor patterns such as scratching. We used the semi-intact spinal cord of adult turtles (Keifer and Stein, 1983; Petersen et al., 2014) and recorded from the segments D8-D10. These segments contain the essential CPG circuits (Mortin and Stein, 1989). Most of the spinal cord including the sensory periphery is left intact. The blood is replaced and the spinal column is provided with oxygenated Ringer's solution so that the neurons and the network have optimal conditions. In this experimental situation the motor behavior is as close to in vivo situation as possible, and is indistinguishable from the intact condition (Keifer and Stein, 1983). The turtle preparation allow for mechanical stability and the turtle’s resistance to anoxia allow for remarkable durability of both the recording conditions and the motor pattern reproducibility (Vestergaard and Berg, 2015).
Adult red-eared turtles (Trachemys scripta elegans) of either sex were placed on crushed ice for 2 hr to ensure hypothermic anesthesia. The turtles were killed by decapitation and the blood was substituted by the perfusion with a Ringer’s solution containing (mM): 100 NaCl; 5 KCl; 30 NaHCO3; 2MgCl2; 3CaCl2; and 10 glucose, saturated with 95% O2 and 5% CO2 to obtain pH 7.6, to remove the blood from the nervous system. We isolated the carapace containing the spinal cord segments D4-Ca2 by transverse cuts (Keifer and Stein, 1983; Petersen et al., 2014) and perfused the cord with Ringer’s solution through the vertebral foramen , using a steel tube and gasket pressing against the D4 vertebra. We opened the spinal column on the ventral side along D8-D10 and gently removed the dura mater with a fine scalpel and forceps. For each insertion site for the silicon probed, we opened the pia mater with longitudinal cuts along the spinal cord with the tip of a bend syringe needle tip (BD Microlance 3: 27G3/4", 0.4x 19 mm). We performed the cuts parallel to the ventral horn between the ventral roots. The surgical procedures comply with Danish legislation and were approved by the controlling body under the Ministry of Justice.
We used a fire polished tip of a bent glass rod for mechanical stimulation, that was mounted linear actuator. The actuator was controlled with a function generator: frequency, amplitude and duration of the stimulus.
Extracellular recordings were performed in parallel at 40 KHz using a 256 channel multiplexed Amplipex amplifier (KJE-1001, Amplipex). Up to four 64-channel silicon probes were inserted in the incisions perpendicular to the spinal cord from the ventral side. We used the 64-channel Berg silicon probes (Berg64 from NeuroNexus Inc., Ann Arbor, MI, USA) with 8 shanks, and 8 recording sites on each shank arranged in a staggered configuration with 30 μm vertical distance. The shanks are distanced 200 μm apart. Recordings were performed at depths in the range of 400-1000 μm.
The intracellular recordings were performed in current-clamp mode with an Axon Multiclamp 700B amplifier (Molecular devices). Glass pipettes were pulled with a P-1000 puller (Sutter instruments) and filled with a mixture of 0.9 M potassium acetate and 0.1 M KCl. Data were sampled at about 20 kHz with a 12-bit analog-to-digital converter (Axon Digidata 1440a, Molecular devices). We inserted the glass electrodes from the ventral side of D8-D10 perpendicularly to the spinal cord. Neurons were located at depths ranging from about 300–800 μm. Typically we had stable intracellular recordings for multiple trials.
Electroneurogram (ENG) recordings were performed with suction electrodes. The scratch behavior was measured by the activity of the nerves: Hip Flexor, Knee Extensor, dD8 and HR-KF. The nerve activities were recorded with a differential amplifier Iso-DAM8 (World Precision Instruments) with bandwidth of 100 Hz–1 kHz.
For histological verification, we combined several staining techniques: The silicon probes were painted with DiI (1–2% diluted in ethanol) before insertion into the spinal cord (Blanche et al., 2005; Vandecasteele et al., 2011). Following successful experiments, we performed Nissl– and ChAT–staining of the tissue, to determine the location of respectively neurons and motoneurons.
The histological processing is detailed in (Petersen et al., 2014). We carefully removed the tissue, perfused it and put it in phosphate buffered saline (PBS) with 4% paraformaldehyde for 24–48 hrs and further stored it in PBS. The tissue was mounted in an agar solution and sliced into 100 μm slices using a microtome (Leica, VT1000 S). The slices were washed with PBS and incubated overnight at 5°C with primary choline acetyltransferase antibodies goat anti-ChAT antibodies (1:500, Milipore, USA) in blocking buffer, which is PBS with 5% donkey serum and 0.3% Triton X-100. The slices were washed three times with PBS and incubated for 1 hr at room temperature with the secondary antibody Alexa488 conjugated to donkey anti-goat antibodies (1:1000 Jackson) in blocking buffer. After three washes with PBS, the slice was mounted on cover slit with a drop of ProLong Gold antifade reagent (Invitrogen Molecular Probes, USA) and cured overnight at room temperature before microscopy. The slice was viewed using a confocal microscope, Zeiss LSM 700 with diode lasers, on a Zeiss Axiolmager M2 using 10x/0.30 EC Plan-Neofluar, 40x/0.6 Corr LD Plan-Neofluar, and 63x/1.40 oil DIC Plan-Apochromat objectives (Zeiss).
The data analysis was primarily done in the programming languages Matlab and Python. The correlation coefficient was calculated as the Pearson product-moment correlation coefficient.
We use skewness (Press et al., 1992) or the third moment as a measure of asymmetry in the distribution around the mean, sometimes referred to as Pearson’s moment coefficient of skewness. It can be estimated using the method of moment estimator as
where are all the observations ( or firing rate) and and are the sample standard deviation and sample mean of the distribution. The skewness is a unitless number and a value of zero indicates perfect symmetry. A positive skew has a tale pointing in the positive direction of the axis and a negative value points in the opposite direction.
Spike sorting was performed in the Klustakwik-suite: SpikeDetekt, KlusterKwik v.3.0 and KlustaViewa (Kadir et al., 2014). Raw extracellular signals were bandpass filtered from 400–9000 Hz, and spikes were detected by a median based amplitude threshold with SpikeDetekt (Takekawa et al., 2012; Kadir et al., 2014; Quiroga et al., 2004). An automatic clustering of the spikes was performed in KlustaKwik, followed by manual cluster-cutting and cluster verification in KlustaViewa. The cluster quality was evaluated by several measures: The shape of the autocorrelation function, the amount of contamination in the refractory period, the Isolation distance (Harris et al., 2001) and the (Schmitzer-Torbert and Redish, 2004) (Figure 2—figure supplement 2). Only well isolated units was used in the further data analysis.
The time-dependent firing rate was estimated by a gaussian kernel by convolving the spike times, , with a Gaussian kernel :
where is defined as
with the bandwidth optimized for each spike train with the sskernel method (Shimazaki and Shinomoto, 2010). The estimated width was in the range of 100–500 μs.
The coefficient is a measure of statistical dispersion and it is defined as a ratio of the areas on the Lorenz curve diagram
where is the area below the line of no dispersion (the diagonal, i.e. ), and is the Lorenz curve, i.e. the cumulative distribution of firing rates (Figure 9H).
The irregularity of the spiking of individual neurons can be described by several measures. The most common measures are the coefficient of variation () and the Fano factor (), but both measures easily overestimate the irregularity when the firing rate is non-stationary (Holt et al., 1996; Ponce-Alvarez et al., 2010; Softky and Koch, 1993). More advanced methods of estimating the time dependent variations in the irregularity have been developed (Shinomoto et al., 2009; Shimokawa and Shinomoto, 2009; Miura et al., 2006), and here we use the widely used metric , which has been suggested to be the most robust measure of local spiking irregularity (Wohrer et al., 2013; Ponce-Alvarez et al., 2010). The time dependent is defined by pairs of adjacent inter-spike intervals and :
where for a Poisson process and for a regular process. can take values in the range from zero to two.
We noticed a small difference in the distribution of irregularity among the neurons recorded with intracellular versus extracellular electrodes (data not shown). The neurons were recorded with intracellular electrodes had more regular spiking than those recorded with extracellular electrodes. This may be caused by a systematic bias in the way the intracellularly recorded neurons were collected, as there is an experimental bias towards high firing rates. Spike sorting processing of the extracellular recordings, on the other hand, is likely to both miss spikes and contain false positives, which will cause overestimation of spiking irregularity.
To get a quantitative handle on the fraction of neurons found in the fluctuation–regime across the population, we consider the distribution of neurons, , which spends a given amount of normalized time in the fluctuation regime, i.e. with . We consider three values of , 0.4, 0.5 and 0.6, as indicators for when the neurons are in the fluctuation–regime. Formally we quantify the time in fluctuation–regime for the population using the reverse cumulative distribution of neurons (Figure 10D–E and Figure 10—figure supplement 1). The reverse cumulative fraction of neurons in the fluctuation regime for a given fraction of normalized time is
This fraction is the fraction of neurons, which spend at least amount of normalized time in the fluctuation regime. To compress the distribution into a single number we use the fraction of time in fluctuation regime of half of the population, , which is the value of for which (arrows and broken lines, Figure 10D–E).
Since the firing rate is rarely constant, one may want to know how many spikes are elicited in the mean– versus fluctuation regime. This is calculated in similar way, using the distribution of neurons having a normalized fraction of spikes in the fluctuation regime, i.e. spikes with , . The reverse cumulative of again gives the fraction of neurons which have at least spikes in fluctuation regime, normalized to 100%,
Again we compress the distribution into a single number and use the fraction of spikes, which occur in fluctuation regime of half of the population, , which is the value of for which (arrows and broken lines Figure 10—figure supplement 1).
We use a definition of the action potential threshold, which is based on the phase plot of versus the derivative . This is the second method reported in Sekerli et al. (2004). The threshold is found as the point in the trajectory in phase space, where there is a strong departure from rest prior to the cycle. Since is proportional to the membrane current, this point represents a strong initiation of the inward current. Defining the slope of in time, , the threshold is defined as the largest peak in second derivative with respect to in phase space, i.e. the maximum of (red dots, Figure 6—figure supplement 1B–C). This is the point with the largest acceleration from baseline prior to the peak of the action potential. The trace was low–pass filtering at 5000 Hz to reduce the vulnerable to electrical noise of the estimates of derivatives.
The method for estimating the response rate as a function of has been described previously (Vestergaard and Berg, 2015). The relationship between firing rate, , and membrane depolarization is based on the assumption that spikes occur as a random renewal point–process i.e. a Poisson process. The rate is directly related to the probability, , of a spike occurring in a small time window at a certain time :
The window has to be small such that the chance of getting more than one spike in the window is negligeble. The firing rate can thus be defined in terms of the probability of achieving a spike in an infinitesimally small time window (Gerstner et al., 2014):
This definition of is also called the ‘stochastic intensity’. Since the probability is strongly dependent on the depolarization of the membrane potential, the firing rate will be similarly dependent. To determine as a function of we have to empirically determine the probability, , for the smallest possible value of , which is the sampling interval of the intracellular recordings. To get as a function of membrane potential, , we first empirically determine the stochastic distribution of prior to the spike (1.5-1.7 ms prior), which we know will cause a spike. Then we normalize this distribution with the amount of time spent at each -level at all time. This is the estimated probability of getting a spike, , within a small time window for a given , i.e. the firing rate as a function of . This empirical method of relating firing rate and was relatively recently invented (Jahn et al., 2011) and used in determining IO properties of e.g. motoneurons (Vestergaard and Berg, 2015). The shape of the spike response function is highly non-linear with upward curvature. This has been observed in previous experiments (using a different method) and has often been referred to as expansive non-linearity (Hansel and van Vreeswijk, 2002; Miller and Troyer, 2002; Murphy and Miller, 2003; Priebe and Ferster, 2005, 2008). An exponential
was fitted to capture the curvature, where the curvature is represented in the exponent , which have units of , and is a constant of units . Such expansive non-linearities have also been investigated in the visual cortex where they are often characterized as a power-law relationship, i.e.
where is a constant and is the power , i.e. supralinear, and often ranging from 2-5 (Hansel and van Vreeswijk, 2002; Miller and Troyer, 2002). This exponent is also a measure of the expansive curvature of the non-linearity. represent a subthreshold level of , where the spiking probability is zero, such that the values in the sampled traces are always larger than , i.e. . The curvature dependence on synaptic fluctuations was assessed by the standard deviation of the distribution of traces prior to the spike in the diffusion regime, i.e. where there was no link to the and the spike occurrence. This distribution was chosen 18 ms prior to the spike (Figure 3B). The analysis and fits were performed in Matlab with generic fitting functions.
In order to distinguish neurons in fluctuation– versus mean–regime, we employ a new metric for quantifying the degree of fluctuations in in between action potentials. We plot the values of in a return map, which is a plot of versus . If the inter–spike has a direct trajectory from the reset potential to the next spike, will smoothly increase and thus will always be larger than . Therefore each point will be above the line of unity (Figure 3—figure supplement 1A–B). On the other hand, if has fluctuations, it will have an indirect and convolved trajectory from the reset value to the threshold. This will manifest as containing values of which are actually smaller than . Thus we use the ratio of points above versus below the unity line as a metric for how convolved and fluctuating the path of is from reset to threshold. If the ratio is 0.5 then is highly fluctuating, whereas if the ratio is approaching 1 the path is straight without any fluctuations. We choose a mean value of the histogram of all values to 0.7 to classify neurons as fluctuation– or mean–driven (Figure 3—figure supplement 1C). This metric of straight versus convolved trajectory had significant negative correlation with other measures of fluctuation– regime, e.g. spike rate skewness, spike irregularity () and least time below threhold (LTBT, Figure 3—figure supplement 1D–F). The choice of is not important as long as it is larger than the timescale of electronic fluctuations of the amplifiers and smaller than the timescale of synaptic fluctuations in . We consistently used ms for all neurons. The return map ratio is intended as a metric to analyze sub-threshold activity and therefore spikes were removed from the traces, including a 6 ms window before and after the peak. Also, the containing the interburst (defined as having ISIs > 300 ms) intervals was removed.
Trilateration is a geometrical process of determining the location of a source in 2D–space using multiple recording sites scattered in space. We adapted the method to take advantage of a distance–dependent decay of the electrical signal from the action potential in the extracellular space. In this way, the amplitudes of the waveforms, which were simultaneously recorded on multiple electrodes, revealed the location of the source in space relative to the position of the electrodes. We assumed that the electrical signal decayed as , where is the distance.
In Figure 2, the following trials were used: for ipsilateral pocket scratch and for contralateral pocket scratch. Data used in Figure 7 has already been published in a different context (Berg et al., 2007). A small subset of the neurons used in Figure 3D–E ( out of 68) has been acquired in a reduced preparation (Petersen et al., 2014) and published for an investigation of a different matter (Berg et al., 2007; Berg and Ditlevsen, 2013). The data from experiments of blockade of inhibition using superfusion of strychnine has also been published previously in the investigation of a different matter (Vestergaard and Berg, 2015). Regarding excluding spikes from the analysis in Figure 3C–E: For the temporal distribution (panel C), only ISIs > 6 ms was included and for the spike triggered -distribution only ISIs > 20 ms was included, all having ISIs < 300 ms. Estimating the FV-curve (Figure 4) all spikes having ISIs > 1.7 ms was included.
Neuronal spiking has traditionally been considered to occur when the mean inward current of the cellular membrane is large enough to cross the rheobase such that the mean membrane potential () is above threshold (). In practice, the mean will not exceed by very much due to the active spiking and after–hyperpolarization, but if this mechanism was turned off the mean membrane current () would drive across threshold, formally written as where is the membrane resistance. Spikes elicited in this manner are in the mean–driven regime (Gerstner et al., 2014; Renart et al., 2007). They have shorter inter–spike intervals (ISIs) because of the large and regular spiking due to the after–hyperpolarization. In contrast, when the mean is below threshold, i.e. , spikes are elicited by temporary fluctuations in due to synaptic bombardment. Such spiking is in the fluctuation–driven regime (Kuhn et al., 2004; Tiesinga et al., 2000; Gerstner et al., 2014; Roxin et al., 2011). The random synaptic fluctuations cause the spiking to be more irregular, which results in a higher coefficient of variation (CV, defined as the standard deviation () divided by the mean of ISIs), than for the mean–driven regime (cf. Figure 1D–E). Therefore irregularity is an indicator of the spiking regime. Another indicator of the fluctuation–driven regime is positive skewness of the firing rate distribution (Figure 1A–B). These indicators are used to quantify the fraction of the population that is in one versus the other regime.
Cellular and network mechanisms of electrographic seizuresDrug Discovery Today. Disease Models 5:45–57.https://doi.org/10.1016/j.ddmod.2008.07.005
'Balancing' of Conductances May Explain Irregular Cortical Firing (Technical Report INC-9502)La Jolla: Institute for Neural Computation, University of California, San Diego.
Synaptic inhibition and excitation estimated via the time constant of membrane potential fluctuationsJournal of Neurophysiology 110:1021–1034.https://doi.org/10.1152/jn.00006.2013
Polytrodes: high-density silicon electrode arrays for large-scale multiunit recordingJournal of Neurophysiology 93:2987–3000.https://doi.org/10.1152/jn.01023.2004
Do pacemakers drive the central pattern generator for locomotion in mammals?The Neuroscientist : A Review Journal Bringing Neurobiology, Neurology and Psychiatry 16:139–155.https://doi.org/10.1177/1073858409346339
Dynamics of networks of randomly connected excitatory and inhibitory spiking neuronsJournal of Physiology-Paris 94:445–463.https://doi.org/10.1016/S0928-4257(00)01084-6
The log-dynamic brain: how skewed distributions affect network operationsNature Reviews. Neuroscience 15:1–15.https://doi.org/10.1038/nrn3687
Cycle-by-cycle assembly of respiratory network activity is dynamic and stochasticJournal of Neurophysiology 109:296–305.https://doi.org/10.1152/jn.00830.2011
The high-conductance state of neocortical neurons in vivoNature Reviews. Neuroscience 4:739–751.https://doi.org/10.1038/nrn1198
Understanding the rhythm of breathing: so near, yet so farAnnual Review of Physiology 75:423–452.https://doi.org/10.1146/annurev-physiol-040510-130049
Coefficient of variation of interspike intervals greater than 0.5. How and when?Biological Cybernetics 80:291–297.https://doi.org/10.1007/s004220050526
Some principles of organization of spinal neurons underlying locomotion in zebrafish and their implicationsAnnals of the New York Academy of Sciences 1198:94–104.https://doi.org/10.1111/j.1749-6632.2010.05539.x
Random walk models for the spike activity of a single neuronBiophysical Journal 4:41–68.https://doi.org/10.1016/S0006-3495(64)86768-0
Neuronal Dynamics: From Single Neurons to Networks and Models of CognitionNew York, NY, USA: Cambridge University Press.
Circuits controlling vertebrate locomotion: moving in a new directionNature Reviews. Neuroscience 10:507–518.https://doi.org/10.1038/nrn2608
Existence and stability of persistent states in large neuronal networksPhysical Review Letters 86:4175–4178.https://doi.org/10.1103/PhysRevLett.86.4175
Balanced excitation and inhibition determine spike timing during frequency adaptationThe Journal of Neuroscience : The Official Journal of the Society for Neuroscience 26:448–457.https://doi.org/10.1523/JNEUROSCI.3506-05.2006
Motoneuron membrane potentials follow a time inhomogeneous jump diffusion processJournal of Computational Neuroscience 31:563–579.https://doi.org/10.1007/s10827-011-0326-z
Neuronal specification in the spinal cord: inductive signals and transcriptional codesNature Reviews. Genetics 1:20–29.https://doi.org/10.1038/35049541
High-dimensional cluster analysis with the masked EM algorithmNeural Computation 26:2379–2394.https://doi.org/10.1162/NECO_a_00661
The Motoneurone and Its Muscle FibresNew York: Oxford University Press.
Locomotor circuits in the mammalian spinal cordAnnual Review of Neuroscience 29:279–306.https://doi.org/10.1146/annurev.neuro.29.051605.112910
Opposing Effects of Intrinsic Conductance and Correlated Synaptic Input on V-Fluctuations during Network ActivityFrontiers in Computational Neuroscience 6:40.https://doi.org/10.3389/fncom.2012.00040
Correlated connectivity and the distribution of firing rates in the neocortexJournal of Neuroscience 29:3685–3694.https://doi.org/10.1523/JNEUROSCI.4500-08.2009
Analysis of firing behaviour of human motoneurones within 'subprimary range'Journal of Physiology-Paris 93:115–123.https://doi.org/10.1016/S0928-4257(99)80142-9
Neuronal integration of synaptic input in the fluctuation-driven regimeThe Journal of Neuroscience : The Official Journal of the Society for Neuroscience 24:2345–2356.https://doi.org/10.1523/JNEUROSCI.3349-03.2004
Efficient solution and performance analysis of 3-D position estimation by trilaterationIEEE Transactions on Aerospace and Electronic Systems 32:1239–1248.https://doi.org/10.1109/7.543845
Intrinsic properties of mouse lumbar motoneurons revealed by intracellular recording in vivoJournal of Neurophysiology 103:2599–2610.https://doi.org/10.1152/jn.00668.2009
Neural noise can explain expansive, power-law nonlinearities in neural response functionsJournal of Neurophysiology 87:653–659.https://doi.org/10.1152/jn.00425.2001
Estimating spiking irregularities under changing environmentsNeural Computation 18:2359–2386.https://doi.org/10.1162/neco.2006.18.10.2359
Distributions of active spinal cord neurons during swimming and scratching motor patternsJournal of Comparative Physiology. A, Neuroethology, Sensory, Neural, and Behavioral Physiology 198:877–889.https://doi.org/10.1007/s00359-012-0758-6
Interspike interval distributions of spiking neurons driven by fluctuating inputsJournal of Neurophysiology 106:361–373.https://doi.org/10.1152/jn.00830.2010
Premotor spinal network with balanced excitation and inhibition during motor patterns has high resilience to structural divisionThe Journal of Neuroscience : The Official Journal of the Society for Neuroscience 34:2774–2784.https://doi.org/10.1523/JNEUROSCI.3349-13.2014
Comparison of local measures of spike time irregularity and relating variability to firing rate in motor cortical neuronsJournal of Computational Neuroscience 29:351–365.https://doi.org/10.1007/s10827-009-0158-2
Numerical Recipes in FORTRAN: The Art of Scientific ComputingCambridge University Press.
The role of spiking and bursting pacemakers in the neuronal control of breathingJournal of Biological Physics 37:241–261.https://doi.org/10.1007/s10867-011-9214-z
Pacemaker neurons and neuronal networks: an integrative viewCurrent Opinion in Neurobiology 14:665–674.https://doi.org/10.1016/j.conb.2004.10.011
On the distribution of firing rates in networks of cortical neuronsThe Journal of Neuroscience : The Official Journal of the Society for Neuroscience 31:16217–16226.https://doi.org/10.1523/JNEUROSCI.1677-11.2011
Estimating action potential thresholds from neuronal time-series: new metrics and evaluation of methodologiesIEEE Transactions on Bio-Medical Engineering 51:1665–1672.https://doi.org/10.1109/TBME.2004.827531
Kernel bandwidth optimization in spike rate estimationJournal of Computational Neuroscience 29:171–182.https://doi.org/10.1007/s10827-009-0180-4
Estimating instantaneous irregularity of neuronal firingNeural Computation 21:1931–1951.https://doi.org/10.1162/neco.2009.08-08-841
Relating neuronal firing patterns to functional differentiation of cerebral cortexPLoS Computational Biology 5:e1000433.https://doi.org/10.1371/journal.pcbi.1000433
Input synchrony and the irregular firing of cortical neuronsNature Neuroscience 1:210–217.https://doi.org/10.1038/659
Chaotic balanced state in a model of cortical circuitsNeural Computation 10:1321–1371.https://doi.org/10.1162/089976698300017214
Divisive gain modulation of motoneurons by inhibition optimizes muscular controlThe Journal of Neuroscience : The Official Journal of the Society for Neuroscience 35:3711–3723.https://doi.org/10.1523/JNEUROSCI.3899-14.2015
Neural network dynamicsAnnual Review of Neuroscience 28:357–376.https://doi.org/10.1146/annurev.neuro.28.061604.135637
Stereological estimate of the total number of neurons in spinal segment D9 of the red-eared turtleThe Journal of Neuroscience : The Official Journal of the Society for Neuroscience 31:2431–2435.https://doi.org/10.1523/JNEUROSCI.3938-10.2011
Population-wide distributions of neural activity during perceptual decision-makingProgress in Neurobiology 103:156–193.https://doi.org/10.1016/j.pneurobio.2012.09.004
Jan-Marino RamirezReviewing Editor; Seattle Children's Research Institute and University of Washington, United States
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Lognormal firing rate distribution reveals prominent fluctuation-driven regime in spinal motor networks" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom, Jan-Marino Ramirez (Reviewer #1), is a member of our Board of Reviewing Editors and the evaluation has been overseen by Eve Marder as the Senior Editor. The other two reviewers involved in the review of your submission have agreed to reveal their identity: Mark Humphries (Reviewer #2) and Alexander Roxin (Reviewer #3).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
This study characterizes the firing rate distribution of neurons within spinal motor networks. The authors have beautifully combined intracellular and large-scale extracellular recordings to reveal the wide-spread existence of "fluctuation-driven" activity in spinal motor networks. These neurons (50% or more) discharge in a "subthreshold" manner. The authors suggest that the existence of these neurons is a strategy to increase stability and sensitivity in the locomotor network. The return map ratio as used here is a metric that has the potential for wide application in a variety of networks. Indeed, the authors nicely discuss their work in the context of a variety of networks. The study elegantly links spiking dynamics to subthreshold activity and therefore shows that the fluctuation-driven versus mean-driven regimes as described in theoretical work appears to be consistent with real networks. The strength of the study is the use of electrophysiology to test a theoretical framework on neuronal dynamics. There are not many studies, specially in spinal circuits that explore this in great detail. Another strength is the use of the turtle model, which allows the characterization of neuronal discharge in a relatively intact animal.
1) The manuscript requires serious editorial revisions to more clearly describe the results. This will be critical to make the text more interesting for a general readership.
– The discussion is of broad interest and compares the data obtained in the turtle spinal cord with data obtained in numerous networks. However, the Introduction is too specialized, and does not provide a general perspective. The reader is immediately confronted with the fluctuation/mean-driven regimes, without an explanation why this is interesting in a general context. It is not clear what the authors mean e.g. with networks being in unresponsive quiescence and saturation. A thorough revision of the Introduction is critical. There is a body of literature on "fluctuation-driven regimes in which excitatory and inhibitory inputs are balanced (e.g. van Vreeswijk and Sompolinsky, Science 1996), and recent in-vivo electrophysiology in cortical circuits supports this idea, e.g. Renart et al. Science 2010. The Introduction fails to make this general point, and also fails to discuss how this study is different from the more "traditional way" to study spinal circuits.
As also indicated by one of the reviewers: "The writing obscures the logic of the work".
– The same need for editorial revision applies to the Results section. Careful proof-reading is needed: see e.g. the first paragraph of the Results makes little sense; Figure 3—figure supplement 2 legend; Figure 3 legend confusing the grey and green histograms; and so on. Most seriously: An entire sub-section of the results is missing, explaining and describing Figure 4A-G
Present the data in a logical manner.i) Establish necessary and sufficient conditions for the existence of "fluctuation-driven" regime:a) Is the subthreshold input approximately Gaussian? [Figure 5]b) Do the neurons have supralinear input-output functions in the subthreshold regime? [Figure 3]ii) How often are neurons in the fluctuation-driven regime? [Figure 6]iii) Does the mode of operation (fluctuation- vs mean-driven) depend on E/I balance? [Figure 4]iv) Can the fluctuation-driven regime be detected across the spinal motor networks? [Figure 2]v) Is it stable between animals and behaviors? [Figure 2]vi) How often and how much of the population is in the fluctuation-driven regime? [Figure 7]
2) Following this logic one notes that a key analysis step is missing. Figure 7 uses the duration of time spent in high irregularity (CV2 > CVcrit) to detect the fluctuation-driven regime, as a proxy for the direct measurements of threshold possible in the intracellular data (Figure 6). In some respects, this is their main result – that 50% of the neurons spend more than 50% of their time in this regime. To show that using CV2 can detect this regime, the authors should apply this same analysis of CV2 to their intracellular data: they need to show it is indeed a proxy for the fluctuation-driven regime.
The authors note that there is different irregularity of firing between the intracellular and extracellular recordings (Consequently, it is not possible to specifically validate the CV2 threshold chosen for the extracellular data. Nonetheless, the authors need to show that in principle this CV2 analysis can recover the fluctuation-driven regime in the intracellular data.
3) In the theoretical work on networks, the network state is considered to be stationary. This means that measures of the spiking activity such as firing rate distributions and CVs of inter-spike-intervals represent very long-time averages and, in fact, in a simulation will converge to the theoretical predictions if ever-longer time intervals are used to evaluate them. In short, there is a proper self-consistent network state for which such statistical measures can be estimated. In the case of the present work the activity is highly non-stationary. It seems the authors take advantage of a sort of separation of time scales in that the bursts, which drive the scratching behavior, as shown in Figure 1D and E for example, are long compared to the inter-spike interval within the burst. However, it seems this non-stationarity should introduce additional variability to any measures of spiking or subthreshold activity beyond the effects of a pure ‘fluctuation-driven' or ‘mean-driven' regime. The authors should find a way to better characterize (quantify) this, since it is one of the important messages of the paper.
The theory on fluctuation-driven versus mean-driven network states does not take into account variable or adaptive thresholds. The authors should also discuss how this might affect such states, e.g. make fluctuation-driven more robust, or less robust? Specifically, a strongly variable threshold would create the appearance of a fluctuation-driven regime, even given approximately constant input (because the timing of ISIs then depends on the noise in the threshold, not the noise in the input). To solve this, the authors could perhaps estimate the degree of non-stationarity by estimating the local variance in the threshold within the burst e.g. computing Var(threshold) in overlapping windows of 10 ISIs. In that way, they could determine if there is only strong variation at the start and end of the bursts (as suggested in some of the example figures), and so would little affect their conclusions.
4) Unfortunately, the authors have not differentiated between interneurons and motoneurons. Yet, they talk about the implications for the rhythm-generating network. Without identifying whether these neurons were either interneurons or motoneurons we don't know whether the difference in the mean/fluctuation driven neuron population is explained by the different types of neurons. The authors should discuss this caveat or even better show for a subset of interneurons (identified for example by staining) that the different types of discharge patterns are present among these neurons. From the data presented here, we don't know whether e.g. the spike distribution characterized here is only representative for motoneurons.
5) The Gaussian distribution of synaptic inputs is not directly demonstrated by actually characterizing the distribution of synaptic events, which is a caveat. The authors deduce this from the somatic membrane potential (Vm = RI) which has a Gaussian distribution – and they show that it does [Figure 3 panels E-G; and Figure 3—figure supplement 2]. Given that this is such an important aspect the authors shouldn't show this key result as an inset of panel A of Figure 3—figure supplement 2. This should be a main figure to illustrate that all neurons have an approximately symmetrical distribution of membrane potentials between spikes (as the skewness is close to 0). And the caveat should be clearly stated that the authors did not perform a characterization of synaptic inputs.
6) Similarly hidden is the definition of the threshold, which is another key point. The authors define the threshold at the maximum slope of dVm/Vm when dVm > 0 [Figure 4—figure supplement 1, panel C], but do not explain it well. The authors need to make this clearer.
Also, in Figure 4—figure supplement 1G, the authors show that the threshold increases for increasing firing rate. In Figure 4A on the other hand, the dashed lines seem to indicate that the threshold is constant during the duration of the burst, despite the fact that the firing rate is clearly changing in time. This apparent discrepancy should be addressed.
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "Lognormal firing rate distribution reveals prominent fluctuation-driven regime in spinal motor networks" for further consideration at eLife. Your revised article has been favorably evaluated by Eve Marder (Senior editor), and three reviewers, one of whom is a member of our Board of Reviewing Editors. The reviewers noted considerable improvement, in both prose and organisation of the content. In terms of the prose, the Introduction now does a much better job of laying out the issues in understanding the dynamical regimes of neuronal networks; and the Discussion nicely links the findings to a range of issues. In particular, with the new results on the physical distribution of the "fluctuating" regime neurons, the Discussion now makes clear the biggest contribution of this study: that we have to think of all neuronal networks, even spinal ones, as acting en masse to generate dynamics, not as a collection of arbitrarily labelled individuals.
The Reviewing editor notes however, there are some remaining issues that need to be addressed before acceptance, as outlined below:
1) Please consider changing the flow of the Abstract:
"When spinal circuits generate rhythmic movements it is important that the neuronal activity remains within stable bounds to avoid saturation and to preserve responsiveness. Here, we simultaneously record from hundreds of neurons in lumbar spinal circuits and establish the neuronal fraction that operates within either a ‘mean-driven’ or a ‘fluctuation-driven’ regime. Fluctuation-driven neurons have a ‘supralinear’ input-output curve, which enhances sensitivity, whereas the mean-driven regime reduces sensitivity. We find…”.
The Introduction has greatly improved! I have some minor comments, though:
2a) In "Which is generated by large circuits primarily in the spinal cord…” replace large circuits with "neuronal network". (Recent studies indicate that locomotion is generated by coupled "microcircuits" – hence the term "large circuits" is very misleading). State "in the spinal cord and medulla" too, since you refer to breathing, which is generated in the medulla.
2b) In the Introduction I would also like to see reference and discussion of the concept that different speeds of locomotion seem to be generated by different microcircuits (see studies by El Manira, Fetcho, etc.). And then differentiate your work from this concept by stating that there must also be a control of different intensities… and then go to your concept of a pool of mean-driven neurons, etc.
3) "This view was essentially predicted much earlier in random walk models [Gerstein and Mandelbrot, 1964]". Follow with a statement like: "However, this concept has been forgotten in explaining locomotion. Yet, it has been adapted to explain cortical processing". Please state something along these lines. I believe it is important to emphasize that this paper tries to apply lessons and approaches that are now commonly used in the cortex to better understand locomotion. In other words, the paper should further crystalize the novelty approach for the field of locomotion, and the opportunities by using approaches that are now commonly used for understanding neuronal circuits in areas other than the spinal cord and brainstem.
4) Small typo: third sentence down in subsection "Normally distributed synaptic input": "in fluctuation-driven regime" should be "in the fluctuation-driven regime".
5) Section "Normally distributed synaptic input": the key definition of fluctuation-driven neurons is here – you go on to relate everything else to this RMR-based definition. But, it is not clear; nor is it said anywhere how many neurons are defined as such. Clearer would be: "…and this forms the basis for selecting neurons in our analysis. An RMR close to 0.5 has fluctuation-driven spiking whereas a value close to 1 has mean-driven spiking (Figure 3—figure supplement 1A,B). Therefore, we defined a neuron as fluctuation-driven if its RMR < 0.7; in our sample of intracellular recordings we found x/68 neurons in this regime."
6) Section "Mean Vm across the population is normally distributed". Three issues here:i) You make use of the threshold here; but do not define it until a few pages later. Is this the same definition of threshold? If so, note it here.ii) "The IO-curve has approximately the same non-linearity across all neurons (Figure 3E)". Figure 3E doesn't address the IO-curve; was this meant to refer to the theoretical scheme (i.e. Figure 1)?iii) Is this analysis only for the fluctuation-driven neurons, or all 68 neurons?
7) Section "Neuronal response-function in subthreshold domain is nonlinear". Again: is this analysis for only the fluctuation-driven neurons, or all 68 neurons?
8) Section "CV2 as an indicator of spiking regime". Reported here is no linear correlation between the time spent below threshold, and the mean CV2 of a neuron. This would seem to be an issue with later using CV2 to diagnose regimes in the population recording. A clear explanation of why this is not an issue would be good – I think the authors are saying that the time below threshold need have no linear relationship with CV2, because different neurons, all in the fluctuation-driven regime, will have different relationships between their CV2 and time below threshold (because of their specific AHP behavior, etc.). But still, Figure 2C implies there should be some relationship between CV2 and time-below-threshold.
That is, we should still be able to see the two broad classes of fluctuating/mean-driven neurons. One wonders if the issue is the same as in Figure 3—figure supplement 1F: the data are not particularly suited for linear correlation. In that figure, correlating LTBT versus RMR, a linear correlation is poor because most of the data-points are clustered close to LTBT=1, and so points inside that cluster dominate the correlation, obscuring any relationship over the whole range of LTBT. So presumably LTBT vs CV2 has the same issue. Perhaps try estimating the distribution of RMR (or CV2) values per LTBT bin e.g. for the bin of, say, LTBT in [0.8 0.9], take the median RMR (or CV2) value. That way, you can correlate (or regress) the medians of the distribution per bin, estimating the relationship between the LTBT and the centre of mass of the values of RMR and CV2.
9) Section "Noisy threshold has no effect". It is stated that "not random, but rather due to a gradual inactivation of Na+-channels throughout the burst." How do we know this is the mechanism?
10) Section "Skewness preserved across behaviours". It is stated that "the ipsilateral behaviour had a slightly higher Gini-coefficient". In Figure 9I can only see 2 out of the 5 animals for which the ipsilateral behaviour has a higher Gini coefficient?
11) Figure 4. Two issues:
i) I think the Results text and/or Figure legend needs a simple explanation for how a firing rate output (y-axis) is derived from the ratio of two histograms (Figure 4A) [as some of the details are already in the Methods]. I understand that this is essentially P(spike) transformed into a firing rate?
ii) Define more carefully what you mean by "sub threshold": Figure 4B shows that "sub threshold Vm" far exceeds the threshold (star)! I think you mean that all these neurons are fluctuation-driven neurons, and so spend the majority of their time sub-threshold.
"In neuronal networks, spikes are generated in either in the mean- or the…” (omit "in").
13) "An intrinsic property, which is commonly believed to be involved in rhythm-generation, is the pacemaker property that can autonomously generate neuronal bursting in the absence of synaptic input [Brocard et al., 2010]". This paper is a good example. Since the authors also discuss their paper in the context of breathing it would be appropriate to also cite: Ramirez et al. 2011, PMID: 22654176 and Ramirez et al. 2004: Pacemaker neurons and neuronal networks: an integrative view (Curr Opin Neurobiol). The authors may also want to discuss Carroll and Ramirez 2013, a paper which discusses the role of pacemaker neurons, and discharge pattern in respiratory rhythm generation using a population approach as similarly applied here for the locomotor network.https://doi.org/10.7554/eLife.18805.028
- Rune W Berg
- Rune W Berg
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Thanks to Gyrgy Buzsáki, Daniel F English and Henrik Lindén for reading and commenting on an earlier version of the manuscript. Funded by the Novo Nordisk Foundation (RB), the Danish Council for Independent Research Medical Sciences (RB and PP) and the Dynamical Systems Interdisciplinary Network, University of Copenhagen.
Animal experimentation: The surgical procedures comply with Danish legislation and were approved by the controlling body under the Ministry of Justice.
- Jan-Marino Ramirez, Reviewing Editor, Seattle Children's Research Institute and University of Washington, United States
© 2016, Petersen 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.