Lognormal firing rate distribution reveals prominent fluctuation–driven regime in spinal motor networks

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. DOI: http://dx.doi.org/10.7554/eLife.18805.001


Introduction
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 (V m ) is well above threshold and spike timing is controlled by after-hyperpolarization 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 eLife digest 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. 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.  Figure 1. Skewness of the rate distribution reveals two regimes of neuronal spiking. (A) In the fluctuation-driven regime the mean input is below the spiking threshold and the IO-curve has a nonlinear shape. A normally distributed input current (shown below x-axis) is transformed into a skewed firing rate distribution (y-axis). (B) In contrast, if the mean input is above threshold, the transformation is linear and the firing rate distribution is symmetric. (C) IO-function for both regimes: Linear for suprathreshold region and nonlinear for subthreshold region. The noise level affects the curvature of the nonlinearity (3 curves illustrate different levels of noise). (D) Sample recordings during motor activity from two spinal neurons in the subthreshold region, where the spiking is irregular and driven by fluctuations, and the supra-threshold region (E), where the mean input is above threshold and spiking is regular. Highlighted area shown at bottom. Spikes in bottom panel are clipped. Tick marks: À50 mV, scale bars: 5 mV. (A-B) adapted from (Roxin et al., 2011). DOI: 10.7554/eLife.18805.003 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 V m 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.

Results
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 Table 1. Two regimes of neuronal spiking and their definition, properties and causes.

Fluctuation-driven Mean-driven Key references
Definition R m I total < V thres R m I total > V thres Brunel, 2000) Properties Lower firing rates Higher firing rates Irregular spiking Regular spiking (Amit and Brunel, 1997;Shadlen and Newsome, 1998;van Vreeswijk and Sompolinsky, 1998) Lognormal/Skewed distribution Symmetric distribution (Buzsáki and Mizuseki, 2014) (Roxin et al., 2011;Mizuseki and Buzsáki, 2013) Cause Balanced E/I Intrinsic currents, unbalanced E/I (Bell et al., 1995;Shadlen and Newsome, 1994;Softky and Koch, 1993) Synchronized excitation (Stevens and Zador, 1998)  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. Mechanisms behind lognormal distribution and the fluctuation-regime

Nissl
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 V m 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 V m was in vicinity of V thres 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 meanregime. Nevertheless, the rhythmic activity possessed a separation of timescales in the sense that the activity between cycles (~1 s) contained much larger excursions in V m than within cycles (~2-400 ms). Here, the mean V m 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 V m 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 V m via Ohm's law under the following conditions. Within a cycle, the mean V m 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 V m . This is likely justified for neurons in fluctuationdriven 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 V m 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 IV m -relationship more linear than in the mean-driven regime.

Normally distributed synaptic input
We intended to test the hypothesis of normally distributed input, but since the approximation of using the variability in V m 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 V m -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 <0:7; 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 V m ( Figure 3A). The distribution was estimated both by selecting the V m in between spikes (temporal distribution) and by collecting instances of V m 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 V m -distributions are not skewed, but rather symmetrical and Gaussianlike (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 V m -distribution, it is difficult to contemplate a non-linear IV m -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 V m -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 V m -distributions of single neurons, which operate primarily in the subthreshold domain, and found that the synaptic input is most likely normally distributed. We do

Membrane potential [mV]
Spike-triggered (18 ms) 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 normallydistributed 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.

Mean V m across the population is normally distributed
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 V m 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 V m from threshold with respect to the size of synaptic fluctuations, i.e. standard deviation of V m (s). This distribution, i.e. ðV m À V thres Þ=s, turns out to also be normally distributed with a mean around 3 s from threshold ( Figure 3-figure supplement 2, plotted for all n ¼ 68 neurons). The value used for V thres here is the mean of the estimated thresholds for all spikes (see below). If we assume, when normalizing V m 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.
Neuronal response-function in subthreshold domain is nonlinear.
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  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 V m in the following way. First, we collected instances of V m shortly before the spike-onset, where V m is depolarized yet still not part of the deterministic spike trajectory. The probability that a given value of V m will cause a spike was estimated as the histogram of V m -instances (gray histogram, Figure 4A) divided by the total time spent at all values of V m (green histogram). This gives the empirical relationship between V m 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 n ¼ 68 neurons and the curvature had a weak negative correlation with the SD of the V m -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.

Lognormal firing rate distribution in single neurons
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 V m was far from V thres 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 (n ¼ 68) 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 logskewness 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 V m 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. The empirical probability of evoking a spike in a small window as a function of V m is determined using spike-triggered overlays. The probability distribution is estimated as the V m -distribution of trajectories prior to spike-onset (gray histogram, 1.7 ms prior to peak) normalized with the total (temporal) V m -distribution (green histogram). Dividing this probability by the sampling interval gives the firing rate (see Materials and methods). (B) The firing rate versus V m for a sample neuron is strongly nonlinear. A power-law (broken line) and an exponential (blue line) are fitted to capture the nonlinearity. Note that the mean threshold ($) is below the largest subthreshold fluctuation (Å), likely due to a depolarization of threshold associated with a higher firing rate (see also

Time spent in regimes: intracellular data
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 V m 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 CV 2 (green line). CV 2 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. CV 2 closer to zero ( Figure 6B). Since the threshold was firing rate-dependent due to the inactivation of the Na þ -conductance ( Figure 6-figure supplement 1) we used the most hyperpolarized value of threshold (broken line). The distribution of CV 2 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 (n ¼ 68) 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.

Transition between regimes by current injection
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 CV 2 and injected current (R ¼ À0:84, p(0.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.

Blocking inhibition causes change in 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 mm) 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 V m (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 (CV 2 ) 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 V m when unbalancing the E/I although it was uncorrelated in the control condition, where the spiking occurred in the fluctuationdriven 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).
CV 2 as an indicator of spiking regime In the above intracellular analyses we reported the spiking irregularity in terms of CV 2 along with the mean V m , current injection and pharmacological manipulation of the balance of excitation and inhibition. The CV 2 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 CV 2 as an indicator of spiking regime. In the above sample cell analyses we note first, that when V m spent a larger fraction of time above threshold, i.e. in mean-driven regime, the CV 2 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 CV 2 was decreased.
To further substantiate CV 2 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 CV 2 is an indicator of the spiking regime, it should be anti-correlated with the return map ratio. This was confirmed by plotting the mean CV 2 for all cells (n ¼ 68) against the mean return map ratio, which indeed demonstrated a significant anti-correlation (R ¼ À0:34, p=0.005) (Figure 3-figure supplement 1E).
A second independent indicator of fluctuation regime is the cumulative time below threshold of V m (Figure 6D), which should be correlated with the mean CV 2 . 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 CV 2 . 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 (n ¼ 68) and plotted it against the mean CV 2 (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 (R log ¼ 0:43, p=0.0003, and R lin ¼ 0:41, p=0.0006), which suggest CV 2 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 CV 2 , if the V m is above threshold (Figure 8, Figure 8-figure supplement 1D). Here, there was a lack of correlation between CV 2 and V m before blocking inhibition, in the fluctuation-driven regime. However, after removal of inhibition, V m was in supra-threshold domain, which introduced an anti-correlation between CV 2 and V m . Hence, if the neuron is in the mean-driven regime the CV 2 is an indicator for the depolarization above threshold. To further verify this we performed a similar test of the relationship between instantaneous CV 2 and local depolarization for all neurons (without pharmacology). We found that all the cells with significant relationships (p<0.05, n ¼ 16=68) had anti-correlation between V m and CV 2 (data not shown). In conclusion, the CV 2 measure is correlated with other measures and indicators of spiking regimes (except the cumulative time below threshold) and therefore CV 2 is a useful indicator in itself.
Noisy threshold has no effect 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 (n ¼ 51) in a sample neuron was s 2 thres ¼ 0:8 mV 2 whereas the variance in synaptic potentials was more than 17-fold higher (s 2 Vm ¼ 14:0 mV 2 ). 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.

Rich diversity in population firing rates
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.

Skewness preserved across behaviors
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 a to the total area a þ b ( 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.

Skewness in firing rate distribution is activity-dependent
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 b Gini = a / (a + b)  achieved by binning the rates in 10 ms windows (Videos 1 and 2). The mean population rate and its SD are indicated as black AE 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 logskewness 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.

Occupancy within regimes across population and time
To further gauge the division of neurons in the two regimes we compared the irregularity of the spiking using CV 2 . This metric was verified above as a reliable indicator of spiking regimes. The distribution of the mean CV 2 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, ISI < 0:5 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. CV 2 > i crit . Choosing i crit is not entirely objective. Complete Poisson-type irregularity has CV 2 ¼ 1, but the spiking is still irregular for lower values (Feng and Brown, 1999). Based on our data, even when the CV 2 » 0.5, the V m spent as much as 96% of the time below threshold ( Figure 6C-D) indicating fluctuation-driven spiking. Further, neurons that had CV 2 » 0.5, also had lognormal firing rate distributions (Figure 7), which also indicates the fluctuation-driven regime. For these reasons, we suggest choosing i crit ¼ 0:5 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   Figure 10 continued on next page 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 CV 2 > i crit . 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 i crit (0.4, 0.5 and 0.6) to indicate the sensitivity to parameter choice. Obviously, choosing a lower i crit results in a larger fraction of time in the fluctuation-driven regime, i.e. the curve is shifted to the right. Choosing i crit larger has the opposite effect. This inverted S-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 (TIF 50 ). Similar TIF 50 -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 TIF 50 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 TIF 50 , we defined a spike-weighted metric as the spikes in fluctuation regime for half the neurons (SIF 50 ). Both the SIF 50 -and TIF 50 -values were relatively conserved across animals as well as behaviors ( Figure 10D-E, Figure 10-figure supplement 1). The large values of TIF 50 and SIF 50 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.

Cell types and spiking 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 CV 2 ) with respect to dorsoventral depth ( Figure 11D). The distributions of mean firing rates (not shown) and the mean CV 2 ( 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.

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

Stability and the two regimes
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 TIF 50 -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 TIF 50 -and SIF 50 -values remain to be determined for other part of the nervous system and in other species.

Rhythm generation and regimes
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 V m ( Table 1). In the mean-driven regime, on the other hand, the net membrane current is so large that the mean V m AEs 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 V m 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.

Cell identity and circuit function
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.
Comparison with other parts of the CNS 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. ðV m À V thres Þ=s had a normal distribution and a distance from mean to threshold of 3s 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 V m of individual neurons is often normally distributed in cortical neurons when considering either the up-or downstate (Destexhe et al., 2003;Stern et al., 1997) and the spiking is irregular with a CV 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 CV 2 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).

Lognormal population firing
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 V m 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 V m across the population is also Gaussian. We did in fact find the distribution of mean V m to be Gaussian (Figure 3-figure  supplement 2F).

Fluctuation-driven regime as a subprimary range in motoneurons?
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 meandriven 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.

Experimental procedures
The experimental procedures are described in more detail at Bio-protocol (Petersen and Berg, 2017). We used the integrated turtle preparation with the spinal motor network intact (n ¼ 5 for the multi-electrode recordings and n ¼ 60 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).

Integrated preparation
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 NaHCO 3 ; 2MgCl 2 ; 3CaCl 2 ; and 10 glucose, saturated with 95% O 2 and 5% CO 2 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.

Network activation
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
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 mm vertical distance. The shanks are distanced 200 mm apart. Recordings were performed at depths in the range of 400-1000 mm.

Intracellular recordings
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 mm. Typically we had stable intracellular recordings for multiple trials.

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

Histology
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 ChATstaining 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 mm 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).

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

Skewness of distribution
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 x 1 ; :::; x N are all the observations (V m or firing rate) and s and x 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
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 L ratio (Schmitzer-Torbert and Redish, 2004) (Figure 2-figure supplement 2). Only well isolated units was used in the further data analysis.

Time-dependent firing rates
The time-dependent firing rate n was estimated by a gaussian kernel by convolving the spike times, sðtÞ, with a Gaussian kernel kðtÞ: where kðtÞ 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 ms.

Gini coefficient
The Gini coefficient is a measure of statistical dispersion and it is defined as a ratio of the areas on the Lorenz curve diagram where a þ b is the area below the line of no dispersion (the diagonal, i.e. a þ b ¼ 1=2), and b is the Lorenz curve, i.e. the cumulative distribution of firing rates ( Figure 9H).

Irregularity of the spiking activity
The irregularity of the spiking of individual neurons can be described by several measures. The most common measures are the coefficient of variation (CV ¼ s=) and the Fano factor (F ¼ s 2 =), 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 Shimokawa and Shinomoto, 2009;Miura et al., 2006), and here we use the widely used metric CV 2 , 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 CV 2 is defined by pairs of adjacent inter-spike intervals ISI i and ISI iþ1 : where CV 2 ¼ 1 for a Poisson process and CV 2 ¼ 0 for a regular process. CV 2 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.
TIF 50 and SIF 50 : time and spikes in fluctuation regime based on 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, f ðtÞ, which spends a given amount of normalized time t in the fluctuation regime, i.e. with CV 2 > i crit . We consider three values of i crit , 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 FðtÞ for a given fraction of normalized time t is This fraction FðtÞ is the fraction of neurons, which spend at least t 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, TIF 50 , which is the value of t for which FðtÞ ¼ 50% (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 CV 2 > i crit , f ðsÞ. The reverse cumulative of f ðsÞ again gives the fraction of neurons which have at least s 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, SIF 50 , which is the value of s for which FðsÞ ¼ 50% (arrows and broken lines Figure 10-figure supplement 1).

Estimating threshold
We use a definition of the action potential threshold, which is based on the phase plot of V m versus the derivative dV m =dt. 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 dV m =dt is proportional to the membrane current, this point represents a strong initiation of the inward current. Defining the slope of V m in time, f ¼ dVm dt , the threshold is defined as the largest peak in second derivative with respect to V m in phase space, i.e. the maximum of d 2 f dV 2 m (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 V m trace was low-pass filtering at 5000 Hz to reduce the vulnerable to electrical noise of the estimates of derivatives.

Spike rate versus V m (FV-curve)
The method for estimating the response rate as a function of V m has been described previously (Vestergaard and Berg, 2015). The relationship between firing rate, n, 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, P, of a spike occurring in a small time window at a certain time t: Pðt; t þ DtÞ ¼ nDt The window Dt 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 : This definition of n is also called the 'stochastic intensity'. Since the probability P is strongly dependent on the depolarization of the membrane potential, the firing rate will be similarly dependent. To determine n as a function of V m we have to empirically determine the probability, P, for the smallest possible value of Dt, which is the sampling interval of the intracellular recordings. To get P as a function of membrane potential, PðV m Þ, we first empirically determine the stochastic distribution of V m 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 V m -level at all time. This is the estimated probability of getting a spike, P, within a small time window Dt for a given V m , i.e. the firing rate as a function of V m . This empirical method of relating firing rate and V m 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;Ferster, 2005, 2008). An exponential nðV m Þ ¼ ce bVm was fitted to capture the curvature, where the curvature is represented in the exponent b, which have units of 1=mV, and c is a constant of units 1=s. Such expansive non-linearities have also been investigated in the visual cortex where they are often characterized as a power-law relationship, i.e. nðV m Þ ¼ k½V m À E a a where k is a constant and a is the power >1, 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. E a represent a subthreshold level of V m , where the spiking probability is zero, such that the values in the sampled traces are always larger than E a , i.e. V m > E a . The curvature dependence on synaptic fluctuations was assessed by the standard deviation of the distribution of V m traces prior to the spike in the diffusion regime, i.e. where there was no link to the V m 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.

Return map ratio: Intracellular metric for mean-vs. fluctuation-regime
In order to distinguish neurons in fluctuation-versus mean-regime, we employ a new metric for quantifying the degree of fluctuations in V m in between action potentials. We plot the values of V m in a return map, which is a plot of V m ðtÞ versus V m ðt þ DtÞ. If the inter-spike V m has a direct trajectory from the reset potential to the next spike, V m will smoothly increase and thus V m ðt þ DtÞ will always be larger than V m ðtÞ. Therefore each point will be above the line of unity (Figure 3-figure supplement 1A-B). On the other hand, if V m has fluctuations, it will have an indirect and convolved trajectory from the reset value to the threshold. This will manifest as containing values of V m ðt þ DtÞ which are actually smaller than V m ðtÞ. Thus we use the ratio of points above versus below the unity line as a metric for how convolved and fluctuating the path of V m is from reset to threshold. If the ratio is 0.5 then V m 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 (CV 2 ) and least time below threhold (LTBT, Figure 3-figure supplement 1D-F). The choice of Dt 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 V m . We consistently used Dt ¼ 1:5 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 V m containing the interburst (defined as having ISIs > 300 ms) intervals was removed.

Determining cellular location using trilateration
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 distancedependent 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 1=r 2 , where r is the distance.

Data selection
In Figure 2, the following trials were used: n ¼ ½6; 4; 9; 5; 6 for ipsilateral pocket scratch and n ¼ ½6; 3; 10; 5; 6 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 (n ¼ 10 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 V m -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.

Definition of fluctuation-and mean-driven spiking
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 (V m ) is above threshold (V thres ). In practice, the mean V m will not exceed V thres by very much due to the active spiking and after-hyperpolarization, but if this mechanism was turned off the mean membrane current (I m ) would drive V m across threshold, formally written as I m thres=R m where R m is the membrane resistance. Spikes elicited in this manner are in the mean-driven regime Renart et al., 2007). They have shorter inter-spike intervals (ISIs) because of the large I m and regular spiking due to the after-hyperpolarization. In contrast, when the mean V m is below threshold, i.e. I m < V thres =R m , spikes are elicited by temporary fluctuations in V m 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 (s) 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.