Why Firing Rate Distributions Are Important for Understanding Spinal Central Pattern Generators

Networks in the spinal cord, which are responsible for the generation of rhythmic movements, commonly known as central pattern generators (CPGs), have remained elusive for decades. Although it is well-known that many spinal neurons are rhythmically active, little attention has been given to the distribution of firing rates across the population. Here, we argue that firing rate distributions can provide an important clue to the organization of the CPGs. The data that can be gleaned from the sparse literature indicate a firing rate distribution, which is skewed toward zero with a long tail, akin to a normal distribution on a log-scale, i.e., a “log-normal” distribution. Importantly, such a shape is difficult to unite with the widespread assumption of modules composed of recurrently connected excitatory neurons. Spinal modules with recurrent excitation has the propensity to quickly escalate their firing rate and reach the maximum, hence equalizing the spiking activity across the population. The population distribution of firing rates hence would consist of a narrow peak near the maximum. This is incompatible with experiments, that show wide distributions and a peak close to zero. A way to resolve this puzzle is to include recurrent inhibition internally in each CPG modules. Hence, we investigate the impact of recurrent inhibition in a model and find that the firing rate distributions are closer to the experimentally observed. We therefore propose that recurrent inhibition is a crucial element in motor circuits, and suggest that future models of motor circuits should include recurrent inhibition as a mandatory element.

1. Introduction

Although it is known that the core neural elements of rhythmic movement, the central pattern generators (CPGs), are located in the spinal cord and the medulla, the neuronal architecture of these networks has remained perplexing. Several working hypotheses for the principle behind generation of movements have been proposed, e.g., muscle synergy and traveling wave (Cuellar et al., 2009; Saltiel et al., 2016, 2017; Yokoyama et al., 2017), multiple unit burst generators (Grillner, 1981), and multilayered half-center organization CPG (Ivanenko et al., 2006; McCrea and Rybak, 2008). A common theme in the literature is the half-center organization inspired by Brown (1914), where two rhythm generating modules, which have recurrent excitation, are coupled reciprocally via inhibitory populations to ensure an alternating flexor and extensor activity (McLean and Dougherty, 2015; Kiehn, 2016; Grillner and El Manira, 2020). Using optogenetics and light activation or inhibition of spatially restricted regions in the spinal cord it was possible to exclusively activate either flexor or extensor rhythms suggesting a modular organization (Hägglund et al., 2013). A modular organization has also been suggested based on analysis of rodent gaits and cellular ablation (Bellardita and Kiehn, 2015) and gaits of human infants (Sylos-Labini et al., 2020). It has been widely suggested that the rhythm generating modules are composed of a recurrently connected excitatory networks, where the rhythm is generated by a subset of neurons with pacemaker properties (Grillner and El Manira, 2020). Nevertheless, modular structures has been difficult to further isolate experimentally (Auyong et al., 2011a,b), and many circuit elements seem to be involved (McLean et al., 2008; Pham et al., 2020) and therefore the details of their architecture, e.g., whether they indeed are composed of recurrent excitation, have been difficult to substantiate. From a stability point of view, however, it is well-known that a population of recurrently connected excitatory neurons without inhibitory interneurons have a propensity to increase their activity in a catastrophic runaway manner (Hennequin et al., 2017; Berg et al., 2019). So, what decides if a network is stable and would a network composed purely of recurrent excitation be sensible for generation of motor activity?

2. Recurrent Networks: Stable or Unstable?

There are many different types of network structures, and the topology of a network determines its stability. This is especially true for recurrent excitatory networks, since activity can create more activity in reverberation and runaway activity. However, recurrent excitatory networks can also become completely silent. So, what decides the stability of an recurrent excitatory network? We may gain some intuition by viewing the recurrent network as a tree-like structure with feed-forward motifs (Figure 1). The overall topology can be approximated as a tree-like network, also termed a branch–process, which is simpler to understand. This tree-like assumption has been shown to be a good approximation for many types of real recurrent networks (Melnik et al., 2011).

FIGURE 1

Figure 1. Dynamics of a recurrent excitatory network. (Left) A recurrent excitatory network often has embedded tree-like structures, e.g., the dark nodes. (Right) The stability of networks with tree-like structure are more simple to characterize. They are either subcritical (left), critical (middle), or supercritical (right). Activity in supercritical networks are unstable and will activate the entire network as exponential growth, whereas the activity in subcritical networks will decay exponentially. Adapted with permission from Larremore et al. (2014).

Branching processes was first studied by Galton and Watson without relation to neuroscience, but in regards to extinction of aristocratic families, passing down family names from one generation to the next (Watson and Galton, 1875). Another example of branching processes is nuclear chain reactions. Fission is utilized in nuclear power plants, and this process also has a tree-like structure. Keeping the reaction going requires careful regulation between activity and curbing the fission process.

The activity that propagates in a tree-like network can either be subcritical (left, Figure 1), critical (middle) or supercritical (right). Whether a network is sub-critical or supercritical depends on the so-called the “branching ratio,” i.e., the expected number of action potentials in the receiver population that an action potential induces. If the branching ratio is below 1, the neuronal activity will rapidly die out (subcritical). If the branching ratio is above 1, it will rapidly increase with exponential growth until the whole population is active (supercritical). If the network is critical, i.e., between sub- and supercritical, there is a simple power-law relation between the number of activated neurons, n, i.e., the size of an avalanche, and the probability (Beggs and Plenz, 2003; Larremore et al., 2014):

$\begin{array}{l}p(n)\propto n^{-3/2}\end{array}$

For activity to be between silence and runaway activity, i.e., the critical propagation, the branching ratio has to be exactly 1 (Beggs and Plenz, 2003), which makes it unlikely that the network is in a critical state by accident. Rather it likely has some self-organization toward a critical states in order to possess this property (Hesse and Gross, 2014; Rybarsch and Bornholdt, 2014). If there is no such self-organization, the network is expected to be either subcritical (quiescent) or supercritical i.e., activity is rapidly propagating throughout the whole population.

So, are CPG modules subcritical or supercritical? Since the neural activity in subcritical networks rapidly dyes out and remain silent, supercritical networks must constitute the basis for rhythm generation in spinal circuits. Spinal motor circuits may involve multiple modules and layers. To keep things simple, let us imagine a simple half-center model that is composed of two supercritical modules with reciprocal inhibition (Figure 2A). When one module becomes active, the network first undergoes a runaway escalation of firing rates where excitation rapidly propagation throughout the network. We consider this “the exponential phase” (Figures 2B,C). Once the whole population is reverberating, the neurons reach their maximum firing rates, and the distribution of firing rates become narrow at the highest possible value. In this saturated phase, adaptation starts to be engaged, e.g., via accumulation of intracellular calcium+, and the firing rates drops back down to zero. After some time, the neurons recover and the cycle can start over. Now we have gained an intuition of the basic principle of recurrent excitatory networks, let us examine how the peak firing rates would be distributed in a network-based model of a half-center module.

FIGURE 2

Figure 2. The dynamics of a burst-generator module has 3 phases. (A) The half-center organization consists of modules: flexor (F) and extensor (E) modules, which generate rhythmic bursting activity and these are reciprocally connected via inhibitory modules to ensure F/E alternation. (B,C) The F/E modules are composed of excitatory neurons with pacemaker properties, which are recurrently connected. The module has three phases: (a) When a neuron fire action potentials (arrows) these quickly activate other neurons, since it is supercritical, resulting in an exponential increase in the network activity. More activity leads to more spikes activating more neurons. Once the activity has propagated to the whole network, the neurons reach their max firing rate the network reach saturation (phase b). Shortly there after the adaptation/pacemaker property turns off the firing in the adaptation phase c, and the activity decay back toward quiescence. During the adaptation phase the quiescence also cease the inhibition to the opposite module, allowing a similar cycle to take place there with a delay. Hence, the population alternates between quiescence and maximal firing rates for all the cells in the module.

3. Dynamics of a Recurrent Excitatory Network

In the following we consider a minimal network of recurrently connected excitatory neurons. Let us assume that the excitatory neurons are recurrently connected by synapses which have a wide distribution (and no negative strengths; Figure 3A). The synaptic connection strengths are randomly drawn from this distribution resulting in a connectivity matrix (Figure 3B). If the overall strength of recurrent connections is strong enough the network becomes super-critical, so that a very small amount of spontaneous activity or external activation will result in a burst of neuronal activity. In such a burst all neurons fire at the maximum firing rate in a synchronized manner when provided by different external drives (Figures 3C,D). After a short time intrinsic adaptation in the neurons becomes active to suppress and ultimately terminate the population burst (Figure 3D, gray traces). The adaptation is also the reason that the network response is shorter than the input pulse (cf. Figures 3C,D). An important aspect with this cyclic activity is that the firing rate of the constituent neurons reach their peak firing rates simultaneously, and therefore the distribution of rates across the population is narrow and clustered around the maximum (Figures 3E,F). Hence, if this model is meaningful, it should be relatively easy to verify whether a rhythm—generating module is composed of excitatory recurrent connections by observing the distribution of firing rates across the population and identify those neurons who display a high peak firing rate. Surprisingly, however, the distribution of firing rates have rarely been reported neither in experimental studies nor in theoretical investigations of spinal motor circuits. As we will see next, however, the few experimental observations of firing rate distributions do not seem to be in accord with a high and narrow distribution of firing rates.

FIGURE 3

Figure 3. All-or-none behavior of recurrent excitatory networks: model. (A) Distribution of synaptic strengths in the model is a truncated Gaussian, with many having zero (no connection). (B) Connectivity matrix of the non-zero connections in the network composed of 200 excitatory neurons. (C) The external drive to the network has to be kept negative (inhibitory) to keep it silent because the network is unstable (supercritical). Positive pulses was given with increasing amplitude imitating a rhythmic drive. (D) The firing rate of the neuronal population (gray overlay) has a prompt increase to maximal firing once the external drive becomes positive. Intrinsic adaptation (green traces) curbs the activity and turns the firing off after some time. Note the population firing is at maximum regardless of input size. (E,F) The linear and log distributions of peak firing rates across population have a single mode at maximal firing (arrows) regardless of the input.

4. Experimentally Observed Firing Rate Distributions

Our survey of the literature provides two observations: First, a distribution of firing rates across neuronal population has rarely been reported (Petersen and Berg, 2016; Cuellar et al., 2018). Second, those experimental reports that do provide the firing rate distribution, they seem in impressive agreement. The precise firing rates vary from experiment to experiment and region to region, but their distributions all seem to have the same long-tail and skewness toward zero. Spinal interneuronal recordings from cervical CPG of the mudpuppy (Figure 4A), show rhythmic (black) and nonrhythmic discharging units during locomotion, both types have distributions lopsided toward the origin. Recordings from the spinal cord in awake macaque monkeys, which were trained to perform visually guided flexion and extension of their wrist in an active ramp provided similar insight (Figure 4B). During both flexion (left) and extension (inset) the distribution of firing rates of interneurons in the cervical spine had a strong skewness toward origin qualitatively similar to the mudpuppy data, in spite of their differences in species and motor function. Interneuronal recording from respiratory motor circuitry in the thoracic spinal cord of cats also provided insight to firing rate distributions in motor circuitry (Figure 4C). Left column: inspiratory interneurons, right column: expiratory interneurons. Top: identified by antidromic stimulation, bottom both antidromic and by location, white and hatched, respectively. Last, recordings from the lumbar spinal cord of turtle during scratching had skewed firing rate distribution, that closely resembles a log-normal distribution (blue fit, Figure 4D), i.e., normally distributed on a log-scale (inset). It is unknown whether the exact shape of the distribution is best described by a log-normal distribution or other skewed distributions, e.g., the gamma-distribution.

FIGURE 4

Figure 4. Firing rate distributions across experimental preparations are skewed toward the origin. (A) The cervical CPG of the mudpuppy while performing locomotion, show rhythmic (black) and nonrhythmic discharge. (B) Cervical interneurons in macaque monkeys performing wrist flexion and extension. (C) Interneuronal recording from respiratory motor circuitry in the thoracic cord of cats. Left: inspiratory interneurons, right: expiratory interneurons. Top: identified by antidromic stimulation, bottom inclusion identification by location, white and hatched, respectively. (D) Lumbar spinal cord of turtle during scratching has skewed firing rate distribution (left) that is bell-shaped on log scale (blue fit). Reproduced with permission (A) Cheng et al. (2002), (B) Prut and Perlmutter (2003), (C) Kirkwood et al. (1988), and (D) Petersen and Berg (2016).

Due to the rareness of experimental reports, we include and analyze a data set from zebrafish larvae during active locomotion in response to visual input, that was kindly made available (Severi et al., 2018). It is known that inhibitory interneuron V1 and V2b provide feedback inhibition (Callahan et al., 2019; Sengupta et al., 2021) most notably via Renshaw cells. In this preparation it is possible to record the fluorescent signal from a genetically encoded calcium sensor (GCaMP5G), which is expressed in glycinergic interneurons by the domain for expression of the transcription factor “engrailed1b,” which is highly correlated with the onset of locomotion. The fluorescent signal in these inhibitory neurons is an indication of the degree of increase in spiking activity. Hence, using the maximal fluorescent signal across the population, we found the distribution to be skewed, remarkably close to a log-normal (Figure 5). The distribution of fluorescent signal could be interpreted in various ways. First, the activity across the cells is uniform, but the distribution of soma size—hence the fluorescent signal—could be log-normal, although the cells have equal firing rates. Second, the cell size if normally distributed, but their firing rates are log-normally distribution, hence giving a fluorescent signal with a skewed distribution. The increase in fluorescent signal may not ad linearly, as the spiking activity increases. Nevertheless, the observation of a skewed distribution (Figure 5) appear in accord with the previous observation and the most parsimonious explanation is that the neurons has a log-normal distribution of firing rates (blue fit, Figure 4).

FIGURE 5

Figure 5. Neuronal imaging of calcium-dependent fluorescent signal in zebrafish larvae hindbrain during visuomotor response has a skewed lognormal-like distribution in peak signal. (A) In-vivo tail-bending in response to visual input. (B) Fluorescent signal (ΔF/F) at neuronal resolution of a glycinergic subpopulation which provide inhibitory feedback during escape response (50 cells shown). (C) The maximal calcium response across cells in the population is skewed and lognormal-like (lognormal fit, gray line). (D) Same distribution on log scale resembles a normal distribution (gray fit). (E) High skewness across the cohort (n = 10 animals). (F) Skewness on a log scale is closer to zero, i.e., normal distribution. Data kindly provided by Severi et al. (2018) and reanalyzed with permission.

In summary, the meta-analysis indicates a profile of activity across the population, which is skewed toward zero, with a long tail toward higher firing rates. This is at odds with the prediction from a recurrent excitatory network as a driver, where the firing rates should cluster around the maximal firing rate and zero. What could explain this apparent lack of congruence between the conventional model and experimental data?

5. Could Reciprocal Inhibition Help?

One of the cornerstones of the half-center organization is reciprocal inhibition. Could reciprocal inhibition help stabilize the firing rates so the population does not switch between silence and maximal firing? The short answer is “no.” Reciprocal inhibition exerts its effect on the antagonist module, which is either already silent, or will become silent once the reciprocal inhibition performs its action. When the activity of the antagonistic module is reduced or completely silent, the returning reciprocal inhibition is also reduced or silent. Hence, reciprocal inhibition actually causes mutual dis-inhibition, which is, in effect, positive feedback. Thus, the reciprocal inhibition in a half-center organization actually provides no help in dampening the firing rates, but rather makes the situation worse, by removing inhibition when it is needed.

6. Inclusion of Recurrent Inhibition

Although there are indications of a sparse connectivity in the CPG structure (Carroll and Ramirez, 2013; Radosevic et al., 2019) the architecture and connectomics of circuits generating motor programs are largely unknown. We know that reciprocal inhibition is indirectly enhancing the instability via dis-inhibition and therefore it does not provide a solution to the problem. So, what is the simplest means to remedy the divergence between these states? As it turns out, recurrent inhibition has an important element in stabilizing the spinal activity. It is known that firing rate distributions are broad for recurrent networks with a balanced between inhibition and excitation (Vogels et al., 2005; Hennequin et al., 2017) and strongly skewed (van Vreeswijk and Sompolinsky, 1996) as seen in experiments. Recurrent inhibition pulls the membrane potential to be less depolarized and in this way the neuron will spike at a lower rate. This also has the effect that more neurons fire action potentials in the fluctuation–driven regime rather than in the mean–driven regime. It is known that at least half of the neurons are active in the fluctuation driven regime during rhythmic scratching (Petersen and Berg, 2016; Berg, 2017). At the same time balanced inhibition and excitation has been observed in at least a subset of neurons in various motor networks (Berg et al., 2007; Petersen et al., 2014; Vestergaard and Berg, 2015; Ramirez and Baertsch, 2018). To illuminate this further, we expanded the model of a recurrent excitatory network presented earlier (Figure 3) to include recurrent inhibition, which is also known as “a balanced network.” While the pure excitatory network had an all-or-none firing activity, the inclusion of recurrent inhibition to replace the role of intrinsic adaptation enriches the network with the capacity to have a graded response to an external drive (Figure 6). Furthermore, the firing rate distribution turns from having a single mode at the maximal firing rate (Figures 3E,F) to having a widely distributed peak firing across the population (Figures 6E,F). The distributions is also similar to those observed in experiments, i.e., skewed toward zeros with a log-normal appearance. The difference between recurrent excitatory networks and networks where recurrent inhibition is included can be summarized as: (1) the firing rate distributions are very different both in width an location of peaks. (2) the all-or-none vs. graded response to external drive (Figure 7).

FIGURE 6

Figure 6. Graded activity in a recurrent excitatory network balanced by inhibition: model. (A) Distribution of synaptic strengths in the model is a Gaussian, with a balance around zero. It also have many connections with zero strength (no connection). (B) Connectivity matrix of the non-zero connections in the network composed of 100 excitatory (blue) and 100 inhibitory (red) randomly connected neurons. (C) The external drive is negative (inhibitory) with positive pulses imitating a rhythmogenic input with increasing amplitude. (D) the firing rate of a subset of the neuronal population has diverse peak firing. (E,F) The linear- and log-scale distributions of peak firing rates are broad and skewed toward zero far from the maximal firing (arrow). The colors represent increasing external input, light blue: lowest, blue: middle, dark blue: largest input.

FIGURE 7

Figure 7. All-or-none vs. graded population response to external drive in two types of model networks: Purely recurrent excitatory network and one including recurrent inhibition, i.e., a balanced network. (A) Firing rates of individual neurons in a recurrent excitatory network (gray) and a recurrent network with both excitation and inhibition (orange), i.e., a balanced network, in response to external drive (bottom trace). All neurons in the excitatory network has the same firing rate, i.e., the maximum or minimum, in response to input, whereas the balanced net has a graded mean response with a wide distribution. (B) The distribution of firing rates across the population in the balanced net (orange) is wide and has a graded response to input, whereas the excitatory network has an all-or-none response that is narrow.

7. Blocking Inhibition Makes the Network Unstable

If recurrent inhibition is indeed a vital element of motor networks, the effect of reducing inhibition should be to increase the overall firing rates and to synchronize the spinal network. Several experiments have been performed where the effects of a systemic block of inhibition by application of either glycinergic antagonists (e.g., strychnine) or GABAergic antagonists (e.g., picrotoxin, bicuculline; Cowley and Schmidt, 1995; Beato and Nistri, 1999; Talpalar et al., 2011, 2013). The general effect is a widespread barrage of intense activity across multiple segments of the spinal cord, which was recorded via the motor nerve output. When reducing inhibition the effective branching ratio grows dramatically, i.e., the expected number of post-synaptic action potential that a presynaptic neuron will increase, and hence the network becomes highly supercritical (Figure 1). A seizure-like barrage of activity quickly spread throughout within milliseconds and reverberates the spinal network until some form of adaptation or fatigue turns it off. This is similar to the dynamics of a recurrent excitatory network (Figure 3). The activity slowly alternates between intense activity with maximal firing of individual neurons and quiescence, which represents a recovery phase.

8. Discussion

In spite of the striking insight into the CPG architecture that firing rates distributions can provide, they have received remarkably little attention in the motor control literature. Other parts of neuroscience have identified distributions as an important element in theoretical neuroscience (Vegué and Roxin, 2019) and e.g., in the processing of sleep regulation (Levenstein et al., 2017) and decision making (Wohrer et al., 2013). The idea of the “log-normal” brain is inspired by the skewed shape of the firing rate distributions (Mizuseki and Buzsáki, 2013; Buzsáki and Mizuseki, 2014).

Here, we inspect the dynamics of excitatory recurrent networks driven by an external drive. Despite the absence of experimental evidence, this type of network architecture has often been proposed to constitute a module in a half-center model or unit burst generator (McCrea and Rybak, 2008; Grillner and El Manira, 2020). Nevertheless, for a recurrent network to not be silent, it has to be supercritical. In such networks, propagation of activity is rapid due to the tree-like structure (Figure 1), which makes it unstable and prone to acceleration and reverberation until all cells fire at maximal rates. This also causes recurrent excitatory networks to essentially work by all-or-none activity: Either it is silent or it has maximal activity (Figure 7). A consequence of this is that it does not respond in a graded fashion to various external drives and has no ability to modulate the output amplitude. This lack of flexibility seems sub-optimal given the large variety of motor task that the body needs to execute, e.g. fine motor skills vs. heavy lifting, and walking peacefully vs. carrying a large load or walking uphill. Hence there is a strong functional argument against an modular architecture of pure recurrent excitatory connectivity.

Furthermore, comparing the firing rate distributions of recurrent excitatory networks with those from experimental reports, they are incompatible: Experiments show broad distributions skewed toward very low rates, akin to a log-normal (Figures 4, 5), whereas recurrent excitatory networks have a single mode at maximal firing, or at zero (Figures 3E,F). A single mode at maximal firing is the exact opposite of the experimentally observed low firing with a wide skewed distribution. Hence, there are functional arguments against the hypothesis of a recurrent excitatory network as well as experimental evidence against it.

So, if there are both experimental and functional arguments against recurrent excitatory modules, what type of network architecture could constitute modules? As the most parsimonious solution, we suggest including recurrent inhibition to curb the activity. Networks that encompass both recurrent excitation and inhibition has previously been considered “heterodox” and “anarchistic” (Grillner and Jessell, 2009), but so far there are no clear experimental evidence against the presence of recurrent inhibition in motor networks. On the contrary, indications of concurrent inhibition have been observed in various medullary and spinal experiments (Ramirez and Baertsch, 2018; Berg et al., 2019). Such balanced networks are also known to have skewed firing rate distributions with low mean firing (van Vreeswijk and Sompolinsky, 1996; Petersen and Berg, 2016). Further, recurrent inhibition is certainly provided via a type of inhibitory interneuron, the Renshaw cell, which receives collateral excitatory fibers from alpha motor neurons, which it also inhibits. Although, the role of these cells is unknown, they may represent just one example of a general and widespread element of recurrent inhibition in spinal motor circuits. Due to the close proximity of Renshaw cells to motor neurons they have been easy to identify, and therefore well-described. Other types of recurrent inhibition in motor network may be less easily identified and therefore likely to receive less attention. There are many other types of inhibitory interneurons (Bikoff et al., 2016), that have remained uncharacterized, because their circuit motifs are more synaptic layers removed from motoneurons.

To test the impact of recurrent inhibition in our simple model, we substituted half of the excitatory neurons in our model with inhibitory neurons, and this changes the dynamics qualitatively in two ways (Figure 6). First, the population response to a different external drive resulted in a graded output (Figure 7), hence allowing a more functional dynamics. Second, the firing rate distribution was wide and skewed toward zero in accord with experimental observations (Figures 6E,F). It is remarkable that such a simple alteration of recurrent excitatory network can change the population dynamics to something more functional and closer to the experimentally observed. Hence, we propose that recurrent inhibition is an indispensable element in motor networks.

9. Methods

To understand the effect of recurrent connectivity of the firing rate distribution of spinal networks we consider a population of N=200 neurons that we model in terms of their firing rates. The firing rate r_i of an example neuron i is determined by a static non-linear firing rate function ϕ(x) that relates that activity variable x_i (analogous to a membrane potential) to the output firing rate: r_i(t) = ϕ[x_i(t)]. We used

$\begin{array}{ll}\varphi (x_i)=\{\begin{array}{cc}{r}_0(1+\text{tanh}\left[(x_i-r_0)/r_0\right]),& \text{for }{x}_i\le r_0\\ r_0+r_{\text{max}}\text{tanh}\left[(x_i-r_0)/r_{\text{max}}\right],& \text{for }{x}_i>r_0\\ \end{array}& (1)\end{array}$

where r₀ represents the rate at the inflection point and r_max = 80 is the maximum deviation of the firing rate from this point. The dynamics of the network is given by

$\begin{array}{ll}\{\begin{array}{c}{\tau }_m{ẋ}_i(t)=-x_i+{\displaystyle {\sum }_j}{J}_{ij}\varphi \left[x_j(t)\right]-g_w{w}_i(t)+I_e(t)\\ {\tau }_a{ẇ}_i(t)=-w_i+\varphi \left[x_i(t)\right]\end{array}& (2)\end{array}$

where J_ij represents the recurrent connections in the network (see below) and I_e(t) is a time-varying external drive. The variable w_i represents an intrinsic adaptation that depends on the firing rate and contributes as a negative input to the activity variable x_i with a strength set by the parameter g_w. The time constants τ_m and τ_a determines the time scale of the activity variable x_i and adaptation w_i, respectively. Here we set τ_m = 50 ms and τ_a = 300 ms. Below we outline two different scenarios used in this study, either a network with recurrent excitation and intrinsic adaptation (Figure 3) or a network with recurrent excitation and inhibition (but lacking intrinsic adaptation) (Figure 6).

Recurrent Excitation and Intrinsic Adaptation

For the case of a network with only recurrent excitation that is counteracted by an intrinsic adaptation, we first generate connection weights J_ij from Gaussian distribution with zero mean and a width (standard deviation) set to σ = 0.05. Negative weights are then set to zero, resulting in a truncated distribution with positive only values, where approximately half of the possible connections are zero (Figure 3A). We set the strength of adaptation g_w to match the average sum of the incoming connection weights to each neuron: $g_w=1/N{\sum }_{ij}{J}_{ij}$.

Recurrent Excitation and Inhibition

In the model with recurrent excitation and inhibition we first generate the connections as for the case above. Here, however, we consider half of the neurons to be excitatory and half of them to be inhibitory and adjust all connections from the inhibitory connections by a factor g_i=-1.5. This results in a network where recurrent connections are dominated by inhibition (Figure 6A). For simplicity we exclude the effect of the intrinsic adaption, i.e., we set g_w = 0.

Software and Code Availability

Numerical simulations of the network model were done using the forward Euler method implemented in custom-written software using Python 3.8. Code for reproducing Figures 3, 6, 7 is available at Berg Lab website (https://berg-lab.net).

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.berg-lab.net.

Author Contributions

Both authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Funding

This work was supported by The Independent Research Fund Denmark (RB), European Research Council (Mobilex, HL), and the Lundbeck Foundation (RB).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Thanks to Claire Wyart and her colleagues for making their zebra fish data set available for further analysis.

References

Auyong, N., Ollivier-lanvin, K., and Lemay, M. A. (2011a). Population spatiotemporal dynamics of spinal intermediate zone interneurons during air-stepping in adult spinal cats. J. Neurophysiol. 106, 1943–1953. doi: 10.1152/jn.00258.2011

CrossRef Full Text | Google Scholar

Auyong, N., Ollivier-lanvin, K., and Lemay, M. A. (2011b). Preferred locomotor phase of activity of lumbar interneurons during air-stepping in subchronic spinal cats. J. Neurophysiol. 105, 1011–1022. doi: 10.1152/jn.00523.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Beato, M., and Nistri, A. (1999). Interaction between disinhibited bursting and fictive locomotor patterns in the rat isolated spinal cord. J. Neurophysiol. 82, 2029–2038.

PubMed Abstract | Google Scholar

Beggs, J. M., and Plenz, D. (2003). Neuronal avalanches in neocortical circuits. J. Neurosci. 23, 11167–11177. doi: 10.1523/JNEUROSCI.23-35-11167.2003

CrossRef Full Text | Google Scholar

Bellardita, C., and Kiehn, O. (2015). Phenotypic characterization of speed-associated gait changes in mice reveals modular organization of locomotor networks. Curr. Biol. 25, 1426–1436. doi: 10.1016/j.cub.2015.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Berg, R. W. (2017). Neuronal population activity in spinal motor circuits: greater than the sum of its parts. Front. Neural Circuits 11:103. doi: 10.3389/fncir.2017.00103

PubMed Abstract | CrossRef Full Text | Google Scholar

Berg, R. W., Alaburda, A., and Hounsgaard, J. (2007). Balanced inhibition and excitation drive spike activity in spinal half-centers. Science 315, 390–393. doi: 10.1126/science.1134960

PubMed Abstract | CrossRef Full Text | Google Scholar

Berg, R. W., Willumsen, A., and Lindén, H. (2019). When networks walk a fine line: balance of excitation and inhibition in spinal motor circuits. Curr. Opin. Physiol. 8, 76–83. doi: 10.1016/j.cophys.2019.01.006

CrossRef Full Text | Google Scholar

Bikoff, J. B., Gabitto, M. I., Rivard, A. F., Drobac, E., Machado, T. A., Miri, A., et al. (2016). Spinal inhibitory interneuron diversity delineates variant motor microcircuits. Cell 165, 207–219. doi: 10.1016/j.cell.2016.01.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, T. G. (1914). On the nature of the fundamental activity of the nervous centres; together with an analysis of the conditioning of rhythmic activity in progression, and a theory of the evolution of function in the nervous system. J. Physiol. 48, 18–46. doi: 10.1113/jphysiol.1914.sp001646

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., and Mizuseki, K. (2014). The log-dynamic brain: how skewed distributions affect network operations. Nat. Rev. Neurosci. 15, 264–278. doi: 10.1038/nrn3687

PubMed Abstract | CrossRef Full Text | Google Scholar

Callahan, R. A., Roberts, R., Sengupta, M., Kimura, Y., Higashijima, S. I., and Bagnall, M. W. (2019). Spinal V2b neurons reveal a role for ipsilateral inhibition in speed control. eLife 8, 1–27. doi: 10.7554/eLife.47837

PubMed Abstract | CrossRef Full Text | Google Scholar

Carroll, M. S., and Ramirez, J.-M. (2013). Cycle-by-cycle assembly of respiratory network activity is dynamic and stochastic. J. Neurophysiol. 109, 296–305. doi: 10.1152/jn.00830.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, J., Jovanovic, K., Aoyagi, Y., Bennett, D., Han, Y., and Stein, R. (2002). Differential distribution of interneurons in the neural networks that control walking in the mudpuppy (Necturus maculatus) spinal cord. Exp. Brain Res. 145, 190–198. doi: 10.1007/s00221-002-1102-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Cowley, K. C., and Schmidt, B. J. (1995). Effects of inhibitory amino acid antagonists on reciprocal inhibitory interactions during rhythmic motor activity in the in vitro neonatal rat spinal cord. J. Neurophysiol. 74, 1109–1117.

PubMed Abstract | Google Scholar

Cuellar, C. A., De La Torre Valdovinos, B., Huidobro, N., Delgado-Lezama, R., Ornelas-Kobayashi, R., and Manjarrez, E. (2018). The spinal neurons exhibit an ON-OFF and OFF-ON firing activity around the onset of fictive scratching episodes in the cat. Front. Cell. Neurosci. 12:68. doi: 10.3389/fncel.2018.00068

PubMed Abstract | CrossRef Full Text | Google Scholar

Cuellar, C. A., Tapia, J. A., Juarez, V., Quevedo, J., Linares, P., Martinez, L., et al. (2009). Propagation of sinusoidal electrical waves along the spinal cord during a fictive motor task. J. Neurosci. 29, 798–810. doi: 10.1523/JNEUROSCI.3408-08.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Grillner, S. (1981). “Control of locomotion in bipeds, tetrapods and fish,” in Handbook of Physiology, The Nervous System II. Motor Control, ed V. B. Brooks (Bethesda, MD: American Physiological Society; Waverly Press), 1179–1236.

Grillner, S., and El Manira, A. (2020). Current principles of motor control, with special reference to vertebrate locomotion. Physiol. Rev. 100, 271–320. doi: 10.1152/physrev.00015.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

Grillner, S., and Jessell, T. M. (2009). Measured motion: searching for simplicity in spinal locomotor networks. Curr. Opin. Neurobiol. 19, 572–586. doi: 10.1016/j.conb.2009.10.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Hägglund, M., Dougherty, K. J., Borgius, L., Itohara, S., Iwasato, T., and Kiehn, O. (2013). Optogenetic dissection reveals multiple rhythmogenic modules underlying locomotion. Proc. Natl. Acad. Sci. U.S.A. 110, 11589–11594. doi: 10.1073/pnas.1304365110

PubMed Abstract | CrossRef Full Text | Google Scholar

Hennequin, G., Agnes, E. J., and Vogels, T. P. (2017). Inhibitory plasticity: balance, control, and codependence. Annu. Rev. Neurosci. 40, 557–79. doi: 10.1146/annurev-neuro-072116-031005

PubMed Abstract | CrossRef Full Text | Google Scholar

Hesse, J., and Gross, T. (2014). Self-organized criticality as a fundamental property of neural systems. Front. Syst. Neurosci. 8:166. doi: 10.3389/fnsys.2014.00166

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanenko, Y. P., Poppele, R. E., and Lacquaniti, F. (2006). Spinal cord maps of spatiotemporal alpha-motoneuron activation in humans walking at different speeds. J. Neurophysiol. 95, 602–618. doi: 10.1152/jn.00767.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Kiehn, O. (2016). Decoding the organization of spinal circuits that control locomotion. Nat. Rev. Neurosci. 17, 224–238. doi: 10.1038/nrn.2016.9

PubMed Abstract | CrossRef Full Text | Google Scholar

Kirkwood, P. A., Munson, J. B., Sears, T. A., and Westgaard, R. H. (1988). Respiratory interneurones in the thoracic spinal cord of the cat. J. Physiol. 395, 161–192. doi: 10.1113/jphysiol.1988.sp016913

PubMed Abstract | CrossRef Full Text | Google Scholar

Larremore, D. B., Shew, W. L., and Restrepo, J. G. (2014). “Critical dynamics in complex networks,” in Criticality in Neural Systems, eds D. Plenz, and E. Niebur (Weinheim: Wiley-VCH Verlag GmbH & Co. KGaA).

Google Scholar

Levenstein, D., Watson, B. O., Rinzel, J., and Buzsaki, G. (2017). Sleep regulation of the distribution of cortical firing rates. Curr. Opin. Neurobiol. 44, 34–42. doi: 10.1016/j.conb.2017.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

McCrea, D. A., and Rybak, I. A. (2008). Organization of mammalian locomotor rhythm and pattern generation. Brain Res. Rev. 57, 134–146. doi: 10.1016/j.brainresrev.2007.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

McLean, D. L., and Dougherty, K. J. (2015). Peeling back the layers of locomotor control in the spinal cord. Curr. Opin. Neurobiol. 33, 63–70. doi: 10.1016/j.conb.2015.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

McLean, D. L., Masino, M. A., Koh, I. Y., Lindquist, W. B., and Fetcho, J. R. (2008). Continuous shifts in the active set of spinal interneurons during changes in locomotor speed. Nat. Neurosci. 11, 1419–1429. doi: 10.1038/nn.2225

PubMed Abstract | CrossRef Full Text | Google Scholar

Melnik, S., Hackett, A., Porter, M. A., Mucha, P. J., and Gleeson, J. P. (2011). The unreasonable effectiveness of tree-based theory for networks with clustering. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 83, 1–12. doi: 10.1103/PhysRevE.83.036112

PubMed Abstract | CrossRef Full Text | Google Scholar

Mizuseki, K., and Buzsáki, G. (2013). PreconfiLgured, skewed distribution of firing rates in the hippocampus and entorhinal cortex. Cell Rep. 4, 1010–1021. doi: 10.1016/j.celrep.2013.07.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Petersen, P. C., and Berg, R. W. (2016). Lognormal firing rate distribution reveals prominent fluctuation-driven regime in spinal motor networks. eLife 5:e18805. doi: 10.7554/eLife.18805

PubMed Abstract | CrossRef Full Text | Google Scholar

Petersen, P. C., Vestergaard, M., Jensen, K. H. R., and Berg, R. W. (2014). Premotor spinal network with balanced excitation and inhibition during motor patterns has high resilience to structural division. J. Neurosci. 34, 2774–2784. doi: 10.1523/JNEUROSCI.3349-13.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Pham, B. N., Luo, J., Anand, H., Kola, O., Salcedo, P., Nguyen, C., et al. (2020). Redundancy and multifunctionality among spinal locomotor networks. J. Neurophysiol. 124, 1469–1479. doi: 10.1152/JN.00338.2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Prut, Y., and Perlmutter, S. I. (2003). Firing properties of spinal interneurons during voluntary movement. I. State-dependent regularity of firing. J. Neurosci. 23, 9600–9610. doi: 10.1523/JNEUROSCI.23-29-09600.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Radosevic, M., Willumsen, A., Petersen, P. C., Lindén, H., Vestergaard, M., and Berg, R. W. (2019). Decoupling of timescales reveals sparse convergent CPG network in the adult spinal cord. Nat. Commun. 10:2937. doi: 10.1038/s41467-019-10822-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramirez, J.-M., and Baertsch, N. A. (2018). The dynamic basis of respiratory rhythm generation : one breath at a time. Annu. Rev. Neurosci. 41, 475–499. doi: 10.1146/annurev-neuro-080317-061756

PubMed Abstract | CrossRef Full Text | Google Scholar

Rybarsch, M., and Bornholdt, S. (2014). Avalanches in self-organized critical neural networks: a minimal model for the neural SOC universality class. PLoS ONE 9:e93090. doi: 10.1371/journal.pone.0093090

PubMed Abstract | CrossRef Full Text | Google Scholar

Saltiel, P., D'Avella, A., Tresch, M. C., Wyler, K., and Bizzi, E. (2017). Critical points and traveling wave in locomotion: experimental evidence and some theoretical considerations. Front. Neural Circuits 11:98. doi: 10.3389/fncir.2017.00098

PubMed Abstract | CrossRef Full Text | Google Scholar

Saltiel, P., d'Avella, A., Wyler-Duda, K., and Bizzi, E. (2016). Synergy temporal sequences and topography in the spinal cord: evidence for a traveling wave in frog locomotion. Brain Struct. Funct. 221, 3869–3890. doi: 10.1007/s00429-015-1133-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Sengupta, M., Daliparthi, V., Roussel, Y., Bui, T. V., and Bagnall, M. W. (2021). Spinal V1 neurons inhibit motor targets locally and sensory targets distally to coordinate locomotion. bioRxiv [Preprint]. 016413. doi: 10.1101/2021.01.22.427836

CrossRef Full Text | Google Scholar

Severi, K. E., Böhm, U. L., and Wyart, C. (2018). Investigation of hindbrain activity during active locomotion reveals inhibitory neurons involved in sensorimotor processing. Sci. Rep. 8, 1–11. doi: 10.1038/s41598-018-31968-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Sylos-Labini, F., Scaleia, V. L., Cappellini, G., Fabiano, A., Picone, S., Keshishian, E. S., et al. (2020). Distinct locomotor precursors in newborn babies. Proc. Natl. Acad. Sci. U.S.A. 117, 9604–9612. doi: 10.1073/pnas.1920984117

PubMed Abstract | CrossRef Full Text | Google Scholar

Talpalar, A. E., Bouvier, J., Borgius, L., Fortin, G., Pierani, A., and Kiehn, O. (2013). Dual-mode operation of neuronal networks involved in left-right alternation. Nature 500, 85–88. doi: 10.1038/nature12286

PubMed Abstract | CrossRef Full Text | Google Scholar

Talpalar, A. E., Endo, T., Low, P., Borgius, L., Hagglund, M., Dougherty, K. J., et al. (2011). Identification of minimal neuronal networks involved in flexor-extensor alternation in the mammalian spinal cord. Neuron 71, 1071–1084. doi: 10.1016/j.neuron.2011.07.011

PubMed Abstract | CrossRef Full Text | Google Scholar

van Vreeswijk, C., and Sompolinsky, H. (1996). Chaos in neuronal networks with balanced excitatory and inhibitory activity. Science 274, 1724–1726.

PubMed Abstract | Google Scholar

Vegué, M., and Roxin, A. (2019). Firing rate distributions in spiking networks with heterogeneous connectivity. Phys. Rev. E 100, 1–15. doi: 10.1103/PhysRevE.100.022208

PubMed Abstract | CrossRef Full Text | Google Scholar

Vestergaard, M., and Berg, R. W. (2015). Divisive gain modulation of motoneurons by inhibition optimizes muscular control. J. Neurosci. 35, 3711–3723. doi: 10.1523/JNEUROSCI.3899-14.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Vogels, T. P., Rajan, K., and Abbott, L. F. (2005). Neural network dynamics. Annu. Rev. Neurosci. 28, 357–76. doi: 10.1146/annurev.neuro.28.061604.135637

PubMed Abstract | CrossRef Full Text | Google Scholar

Watson, H. W., and Galton, F. (1875). On the probability of the extinction of families. J. Anthropol. Inst. Great Britain Ireland 4, 138–144.

Google Scholar

Wohrer, A., Humphries, M. D., and Machens, C. K. (2013). Population-wide distributions of neural activity during perceptual decision-making. Prog. Neurobiol. 103, 156–193. doi: 10.1016/j.pneurobio.2012.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Yokoyama, H., Hagio, K., Ogawa, T., and Nakazawa, K. (2017). Motor module activation sequence and topography in the spinal cord during air-stepping in human: insights into the traveling wave in spinal locomotor circuits. Physiol. Rep. 5, 1–15. doi: 10.14814/phy2.13504

PubMed Abstract | CrossRef Full Text | Google Scholar