Central pattern generators are circuits generating rhythmic movements, such as walking. The majority of existing computational models of these circuits produce antagonistic output where all neurons within a population spike with a broad burst at about the same neuronal phase with respect to network output. However, experimental recordings reveal that many neurons within these circuits fire sparsely, sometimes as rarely as once within a cycle. Here we address the sparse neuronal firing and develop a model to replicate the behavior of individual neurons within rhythm-generating populations to increase biological plausibility and facilitate new insights into the underlying mechanisms of rhythm generation. The developed network architecture is able to produce sparse firing of individual neurons, creating a novel implementation for exploring the contribution of network architecture on rhythmic output. Furthermore, the introduction of sparse firing of individual neurons within the rhythm-generating circuits is one of the factors that allows for a broad neuronal phase representation of firing at the population level. This moves the model toward recent experimental findings of evenly distributed neuronal firing across phases among individual spinal neurons. The network is tested by methodically iterating select parameters to gain an understanding of how connectivity and the interplay of excitation and inhibition influence the output. This knowledge can be applied in future studies to implement a biologically plausible rhythm-generating circuit for testing biological hypotheses.

The current scientific literature suggests three main network architectures for the rhythm-generating spinal circuits, commonly called central pattern generators (CPGs) (Grillner & El Manira, 2020; Marder & Bucher, 2001; Rancic & Gosgnach, 2021). Graham Brown (1911) suggested the existence of CPGs based on his research with decerebrate cats in the early 1900s. He also proposed the first viable architecture for these circuits, the half-center oscillator where two neuronal populations mutually inhibit each other to produce an antiphasic population output for controlling an antagonistic muscle pair (e.g.; flexor/extensor; Graham Brown, 1911, 1914), which was later supported experimentally (Jankowska et al., 1967). Since then, the unit burst generator and two-layer network architectures have been proposed to account for more complex firing patterns observed during rhythmic movement (Ausborn et al., 2021; Grillner & Kozlov, 2021; Rancic & Gosgnach, 2021; Rybak et al., 2015). However, these architectures replicate alternating firing where the neurons in the active population are all firing either at the same or opposite phase and each burst consists of many spikes, covering most of the related active period during locomotion. We term this “neuronal phase” in contrast to the “population phase.” We define neuronal phase as the phase angle relationship between the local maximums of the firing rates of individual neurons firing in the network and the population output from one of the rhythm-generating populations. This provides us with the neuronal phase of any given neuron with respect to the network output. Population phase is used to refer to comparing the phase difference between the output signals of the network.

In contrast to the above models, a study by Kozlov et al. (2007) uses a mathematical model to mimic single spiking and tonic firing of neurons over a particular frequency range. This moves closer to biological recordings that show sparse neural activity (Dougherty et al., 2013; Hart & Giszter, 2010; Musienko et al., 2022; Petersen & Berg, 2016; Zhong et al., 2011). Sparse firing among neurons has further been shown to manifest as a skewed distribution across the population, which is similar to a normal distribution on a logarithmic x-axis, that is, a log-normal firing rate distribution (Berg, 2017; Petersen & Berg, 2016). Although rarely reported experimentally, the distribution of firing rates across the neuronal population provides important clues for the organization of the CPG network topology (Lindén & Berg, 2021). Based on these findings, our study implements a spiking CPG containing both sparsely and abundantly firing neurons within the rhythm-generating populations. The rhythm generators (RGs) are designed with random, sparse connectivity (Radosevic et al., 2019), which is inspired by the balanced sequence generator (BSG) (Lindén et al., 2022). A spiking version of the BSG model using leaky integrate-and-fire neurons was first developed by Najarro et al. (under preparation). We increase the complexity of this scheme by implementing a more detailed neuron model, which is capable of producing bursting and tonic firing patterns. Currents that support bursting activity in individual neurons have been demonstrated in multiple core RG populations (Brocard et al., 2013; Dougherty & Ha, 2019; Dougherty et al., 2013; Shevtsova et al., 2020; Song et al., 2020; Wilson et al., 2005). We therefore decided to explore their inclusion in our models. We implement two spiking RG networks connected in a conventional CPG architecture with reciprocal inhibition between them. The resulting network bridges recent work that implements a random, sparse architecture (Lindén et al., 2022; Radosevic et al., 2019) and the population-based architectures reviewed in Dougherty and Ha (2019) and Rancic and Gosgnach (2021). As this network merges two types of architectures previously modeled distinctly from each other, we call it a hybrid architecture. Our hybrid network places two RGs that are implemented using balanced excitation and inhibition (Berg et al., 2019) in a half-center oscillator architecture using mutual inhibition to produce antiphasic signals. The output of this architecture is able to produce alternating signals from a network with sparse firing. This means an antagonistic pair, such as the flexor and extensor, can be controlled by this network. Notably, the output also shows that the neural activity from the RG populations resembles rotational dynamics, evenly distributed neuronal phases of individual spinal neurons, observed in the lumbar spinal cord (Lindén et al., 2022), ventral respiratory column in the medulla (Ramirez & Bush, 2022), and the primate motor cortex (Churchland et al., 2012; Kalidindi et al., 2021). However, even though our model is developed using biological constraints, it only encompasses the core RGs of a spinal CPG network. Our results therefore do not exhibit perfect rotational dynamics, a neuronal phase representation that covers all phases nearly evenly, as seen in experimental recordings from a larger number of spinal neurons and explained by a new theory (Lindén et al., 2022). Regardless, our work indicates that introducing sparse firing into the RG populations combined with a topology of balance between excitation and inhibition allows for a broader neuronal phase representation. This is in contrast to previously proposed modular architectures for intralimb coordination where strict alternation is observed and all neurons burst with abundant spikes within a single discrete neuronal phase (Ausborn et al., 2018, 2021; Grillner & Kozlov, 2021; Rybak et al., 2015).

The spinal circuit is implemented as a spiking neural network. The hybrid network architecture comprises two RG populations with recurrent excitation and inhibition, a balanced network (Berg et al., 2019; Petersen et al., 2014), arranged in a CPG half-center oscillator configuration (see Figure 1).

Figure 1:

Model schematics. (A) The basic structure of a single RG population. Each RG contains two subpopulations, one excitatory (E) and one inhibitory (I) subpopulation that randomly and sparsely connect both to each other and to themselves (i.e., a balanced network). The inhibitory population is pictured as smaller because it contains fewer neurons based on the selected inhibitory/excitatory ratio (see Table 3, P1, in section 2.5). (B) Schematic showing the basic configuration of the hybrid CPG architecture. The two RG populations are connected via mutual inhibition using antagonistic inhibitory interneuron populations (AInh) to inhibit antagonistic output. Each AInh population is excited by the excitatory populations in the RG and then inhibits all neuron types in the antagonistic RG (pictured as gray lines in panel A). The amount of connectivity is set by the synaptic sparsity of the network.

Figure 1:

Model schematics. (A) The basic structure of a single RG population. Each RG contains two subpopulations, one excitatory (E) and one inhibitory (I) subpopulation that randomly and sparsely connect both to each other and to themselves (i.e., a balanced network). The inhibitory population is pictured as smaller because it contains fewer neurons based on the selected inhibitory/excitatory ratio (see Table 3, P1, in section 2.5). (B) Schematic showing the basic configuration of the hybrid CPG architecture. The two RG populations are connected via mutual inhibition using antagonistic inhibitory interneuron populations (AInh) to inhibit antagonistic output. Each AInh population is excited by the excitatory populations in the RG and then inhibits all neuron types in the antagonistic RG (pictured as gray lines in panel A). The amount of connectivity is set by the synaptic sparsity of the network.

Close modal

2.1  Rhythm Generator

The implementation of each RG is inspired by the balanced sequence generator, a network with recurrent excitation and inhibition where the activity is stable yet close to instability in relation to equilibrium points (Lindén et al., 2022). Each RG contains 2000 neurons divided into two subpopulations, one excitatory and one inhibitory, that connect both to each other and to themselves (see Figure 1A). The balance is indicated as the ratio of excitatory to inhibitory neurons, one of the parameters to be tested (P1) and is further explained in section 2.5. In alignment with theoretical findings, the subpopulations are designed to be random and sparsely connected (Berg et al., 2019; Radosevic et al., 2019; Zhou & Yu, 2018). A sparsely connected network is defined here as having a low probability of connections between neurons; in our final model, this is 3% and 9% depending on the populations.

2.2  Central Pattern Generator Network

The central pattern generator network (CPG) network architecture is based on the flexor/extensor half-center RG architecture (Ausborn et al., 2021, 2018; Kiehn, 2016; Rybak et al., 2015) and comprises two RG populations mutually inhibiting each other via two inhibitory interneuron populations (see Figure 1B). The population of excitatory neurons (E, n=1667) within the RG excites a population of antagonistic inhibitory interneurons (AInh, n=100). The neurons within the AInh populations extend inhibitory connections to all neurons in the antagonistic RG (Rybak et al., 2015). The connectivity between the RG and AInh populations is set using the synaptic sparsity parameter and is explained in section 2.5, parameter P2. The inhibitory subpopulation (I, n=333) in the RG only projects locally within the RG.

2.3  Neuron Model

The spiking neurons are modeled as adaptive exponential integrate-and-fire (AdEx) neurons (Naud et al., 2008) with two different sets of parameters (see Table 1). The majority of neurons fire tonically when activated, while the remainder of the neurons produce bursting behaviors; the percentage of tonic versus bursting neurons is one of the tested parameters (see P3 in section 2.5). Some neuronal parameters are initialized as normal distributions in order to create a variance in characteristics while maintaining overall behavior—tonic firing or bursting. These distributed parameters are indicated by providing their mean and standard deviation (example: 200 ± 40 (mean ± standard deviation)). Additionally, the membrane potential of individual neurons is assigned an initial value of -60 ± 10 mV. Finally, the capacitance, leak conductance, and spiking threshold have been manually tuned to move closer to the physiological output frequency range of 0.27 to 1.84 Hz as reported in Danner et al. (2015).

Table 1:

Parameter Values of Two Neuronal Firing Behaviors (Tonic Firing and Bursting).

ParameterTonic FiringBursting
Vth (mV) -50 -52 ±
tref (ms) 1 ± 0.2 1 ± 0.2 
C (pF) 200 ± 40 500 ± 80 
Ibias (pA) 320 ± 80 160 ± 40 
gL (nS) 10 26 
EL (mV) -70 -60 
ΔT (mV) 
τw (ms) 30 130 
a (nS) -11 
b (pA) 30 
Vreset (mV) -58 -48 
ParameterTonic FiringBursting
Vth (mV) -50 -52 ±
tref (ms) 1 ± 0.2 1 ± 0.2 
C (pF) 200 ± 40 500 ± 80 
Ibias (pA) 320 ± 80 160 ± 40 
gL (nS) 10 26 
EL (mV) -70 -60 
ΔT (mV) 
τw (ms) 30 130 
a (nS) -11 
b (pA) 30 
Vreset (mV) -58 -48 

In Table 1, Vth is the voltage threshold potential, tref is the refractory period, C is the membrane capacitance, Ibias is the bias current, Vm is the membrane potential, gL is the leak conductance, EL is the resting potential, ΔT is the sharpness factor, τw is the adaptation time constant, a is the subthreshold adaptation conductance, b is the spike-triggered adaptation, and Vreset is the reset potential (Naud et al., 2008). The noise is created as gaussian white noise; it is updated every 0.1 ms, which is the designated time step of the simulation. At each time step, the noise is added to the bias current of each neuron. Noise is one of the tested parameters (see P5, Table 3) so it changes for different tests. However, the mean is always 0 pA; this means that bias current can be reduced by the noise if it receives a negative noise value.

Table 2:

Metrics for Measuring the Performance of the Developed Network Architecture.

MetricDesired ResultMethod
Antagonistic output Output activity shows antiphasic population output; peaks should be 180° apart when comparing the output from each RG population. Plot population firing rate output and comparison of time of local maximums of each population’s output (see equation 2.7). 
Firing sparsity A percentage of neurons fire two times per period or less often. Plot a line graph tallying the number of neurons that spike a certain number of times throughout the simulation (example: if 100 neurons spike 10 times, the point is placed at (10,100), showing a bias toward zero). Calculate the percentage of sparsely firing neurons in the RGs. 
MetricDesired ResultMethod
Antagonistic output Output activity shows antiphasic population output; peaks should be 180° apart when comparing the output from each RG population. Plot population firing rate output and comparison of time of local maximums of each population’s output (see equation 2.7). 
Firing sparsity A percentage of neurons fire two times per period or less often. Plot a line graph tallying the number of neurons that spike a certain number of times throughout the simulation (example: if 100 neurons spike 10 times, the point is placed at (10,100), showing a bias toward zero). Calculate the percentage of sparsely firing neurons in the RGs. 
Table 3:

Parameters Tested to Determine Their Effect on Network Output.

IDDescriptionTest Strategy
P1 Ratio of excitatory versus inhibitory neurons in RG Test with 5:1, 4:1, 3:1 excitatory-to-inhibitory neurons (ratio based on biological findings in Sukenik et al., 2021). 
P2 Synaptic sparsity Test with network connectivity up to 10% (RG) and 30% (AInh). 
P3 Neuronal subtype (tonic, bursting) Update percentage of bursting neurons to find the minimal number of bursting neurons necessary. 
P4 Network balance Update the balance by skewing the RGs to either excitation or inhibition through synaptic weights. 
P5 Noise Update the amount of current noise in the network to see how much noise can be added before the output characteristics are compromised. 
IDDescriptionTest Strategy
P1 Ratio of excitatory versus inhibitory neurons in RG Test with 5:1, 4:1, 3:1 excitatory-to-inhibitory neurons (ratio based on biological findings in Sukenik et al., 2021). 
P2 Synaptic sparsity Test with network connectivity up to 10% (RG) and 30% (AInh). 
P3 Neuronal subtype (tonic, bursting) Update percentage of bursting neurons to find the minimal number of bursting neurons necessary. 
P4 Network balance Update the balance by skewing the RGs to either excitation or inhibition through synaptic weights. 
P5 Noise Update the amount of current noise in the network to see how much noise can be added before the output characteristics are compromised. 
Equations 2.1 and 2.2 describe the AdEx neuron’s dynamics:
CdVmdt=-gL(Vm-EL)+gL(ΔT)eVm-VthΔT-ge(t)(Vm-Ee)-gi(t)(Vm-Ei)-w+Ibias,
(2.1)
τwdwdt=a(Vm-EL)-w.
(2.2)
The reset equation for each is
IfVm>0mV,thenVmVresetww+b.
(2.3)
Equation 2.1 defines the change in membrane potential per time step, whereas equation 2.2 outlines the current adaptation (w) in pA and the reset functions for each are in equation 2.3. The parameters are the same as those defined in Table 1 with the addition of the parameters for synaptic conductance. The neural simulation tool (NEST) used in this study (Sinha et al., 2023) models the synaptic conductance within the neuron model equation by adding a calculation for both excitatory (ge(t)(Vm-Ee)) and inhibitory (gi(t)(Vm-Ei)) synaptic conductances into equation 2.1. The term ge(t) is the excitatory conductance, Ee is the excitatory reversal potential, gi(t) is the inhibitory conductance, and Ei is the inhibitory reversal potential. The rise and decay times of the synapses are modeled as alpha functions in the selected NEST neuron model, aeif_cond_alpha. Both excitatory and inhibitory kinetics are provided in equation 2.4:
ge(t)=Wetτsyn_exe-tτsyn_ex,gi(t)=Witτsyn_ine-tτsyn_in,
(2.4)
where We is the excitatory weight, t is the time of the presynaptic spike, τsyn_ex is the excitatory rise time, Wi is the inhibitory weight, and τsyn_in is the inhibitory rise time. The synaptic parameters are left at the default values defined by NEST. Specifically, the excitatory rise time is 0.2 ms and the inhibitory rise time is 2 ms. The synaptic weights are initialized as a distribution and tuned so that each RG population is balanced with a slight excitatory bias up to 10%. The synaptic connections from the inhibitory populations to the postsynaptic RG neurons are two times stronger than the inhibitory weights within the RG itself. This trade-off produces antiphasic population output from the RGs while allowing individual RGs to maintain sparse firing. Balance within each RG and in the overall network is calculated using equations 2.5 and 2.6:
balance=sum_of_excitatory_weights-sum_of_inhibitory_weightstotal_weight,
(2.5)
individual_weight=synaptic_weight*(g(t))*#_of_spikes.
(2.6)
In equation 2.5, the balance is approximated by taking the difference between the sum of all excitatory (excitatory_weight) versus inhibitory weights (inhibitory_weight) and dividing this by the sum of all synaptic weights (total_weight) in the populations being measured. If the balance value is positive, the network skews excitatory; if it is negative, the network skews inhibitory. Equation 2.6 describes the calculation of an individual excitatory_weight or inhibitory_weight. In order to account for the actual impact of each synapse, equation 2.6 multiplies the number of spikes (#_of_spikes) sent on a particular synapse with the synaptic weight (synaptic_weights) of that synapse. Furthermore, the integral of the conductance ((g(t))) is also multiplied in order to account for the different rise times of excitatory and inhibitory synapses. The conductance, g(t), is the same that is referred to in equation 2.4.

2.4  Network Performance

The performance of the network is measured using two metrics in order to make sure the developed model is operating as desired. Table 2 outlines each metric, the desired result, and the method.

The neural network is run using NEST and written in Python, code that is also available on GitHub (Strohmer, 2023). The network data analysis is performed within the Python script. The first metric, antagonistic output, is found by analyzing the output activity from each RG population. This activity is analyzed by detecting spikes for each rhythm-generating population and applying a sliding time window, counting spikes per window to produce a plot of population activity. The time window is 5 ms in length and moves at an interval of 0.1 ms to create a smooth population firing rate approximation from each RG. The signals are then compared by recording the timing of their local peaks. This is used to confirm antiphasic population output by applying equation 2.7,
avg_phase_diff=peaksrg1-peaksrg2avg_period,
(2.7)
where avg_phase_diff is the average population phase difference between the two output signals, peaksrg1 represent the local maximums of the first output signal, peaksrg2 are the local maximums of the second output signal, and avg_period is the average period of both signals over the complete simulation.
The frequency of bursting should be close to the interval 0.27 to 1.84 Hz to align with surface electromyographic (EMG) signals recorded after stimulating a functionally isolated human spinal cord (Danner et al., 2015). The most realistic comparison to the developed network is an isolated spinal cord because the model does not incorporate sensory feedback. The frequency is calculated by subtracting two consecutive peak times of an individual RG output signal and averaging them as described in equations 2.8 and 2.9:
freq=1peak2-peak1,
(2.8)
avg_freq=Σfreq#_of_data_points,
(2.9)
where freq is the frequency of a single population, and peak2 and peak1 are the times of two consecutive peaks recorded from a single population. avg_freq is the average frequency of a single population over the complete simulation time, and #_of_data_points is the total number of peak times recorded for a single population.

The second metric of firing sparsity is measured by counting the number of spikes produced by each neuron throughout the complete simulation. The total number of spikes per neuron is plotted as a line graph with the number of neurons on the y-axis and the number of spikes on the x-axis to confirm that the plot skews toward zero. Additionally, the number of neurons spiking four times or less per second (two times or less per period based on an average population firing rate output frequency of 2 Hz) is counted and divided by the total number of spiking neurons to provide a percentage of sparsely firing neurons.

The network is manually tuned using the operational performance metrics. Antagonistic output is ensured by plotting the population firing rate output of each RG and calculating the average population phase difference (see equation 2.7). Firing sparsity is checked using the percentage of sparsely firing neurons. By tuning the network parameters using these metrics, we are able to define a functional network for testing. The parameters of the control network after tuning are found in section 3.2.

2.5  Network Testing

In order to understand the network, a methodology is constructed to step through changes in the network and neuronal parameters. The parameters are selected in order to test connectivity, excitatory/inhibitory balance, and robustness. One parameter is tested at a time in an iterative process to understand the relationship between the specific parameter and network output. All of the tests are performed using the same random seed to ensure a direct causality between changing a parameter and the change in output (see Table 3).

The ratio of excitatory-versus-inhibitory neurons (P1) is defined by the number of excitatory and inhibitory neurons within each RG. Therefore, an RG with a 5:1 ratio of excitatory-to-inhibitory neurons has five excitatory neurons for every one inhibitory neuron, or 1/6 of neurons are inhibitory. Synaptic sparsity (P2) refers to the network sparsity or the amount of connectivity in the network. A connectivity percentage of 10% means that there is a 10% chance that any presynaptic neuron is connected to any postsynaptic neuron. There are two values used for connectivity—one for connections within the RG and from the RG neurons to the AInh population, the other for connections from the AInh population to the postsynaptic RG neurons. The RG connectivity is three times as sparse, meaning if AInh connectivity is 3%, RG connectivity is 1%. This ratio was found through manual testing where AInh needed greater connectivity to ensure antiphasic population output. Neuronal subtypes (P3) specifically relate to firing behavior when isolated from the network, that is, whether the neurons are innately tonically firing or bursting. Network balance (P4) is updated by changing the weights of the excitatory and inhibitory synapses to bias the RG populations to be overall excitatory, inhibitory, or balanced. This is a valid approach according to biological findings that confirm the human brain regulates balance through changing excitatory/inhibitory connections (Sukenik et al., 2021). Finally, robustness of the network is tested through the addition of noise (P5). This is added directly to individual neurons as current noise: the mean of the noise is always 0 pA, but the standard deviation is updated based on the test. For example, if initial current bias for a neuron is 320 pA and it receives a noise value of -100 pA, the current bias at that time step will be 220 pA. This will be updated again at the next time step when the noise is recalculated. The noise is always applied to the initial current bias, which remains constant throughout a simulation.

In order to further test robustness, another round of tests is performed using the control parameters but changing the initial value (seed) for the random number generator each time. This is completed 20 times to confirm that output characteristics remain consistent across randomly initialized trials. The random initialization seed value determines which neurons are connected and changes the overall initialization of noise in the network. Therefore, changing the random seed tests that the network is not overfitted to produce the desired output for only a single network configuration.

There are two main results of this study: (1) the implementation of a network architecture that is able to produce rhythmic, antiphasic population activity while exhibiting sparse firing and (2) how specific parameters affect network output.

3.1  Network Performance

The behavior of individual neurons as tonically firing or bursting is confirmed prior to connecting them within a network (see Table 1 and Figure 2). For this initial behavior test, the neurons are provided an increasing external current in order to find a suitable current input range based on firing pattern.

Figure 2:

Spiking and membrane potentials of sample neurons. (a) Top: Raster plot of spiking activity for neurons initialized as tonically firing. Bottom: Membrane potential of a representative tonic neuron, corresponding to neuron number 6 in the raster plot. (B) Top: Raster plot for neurons initialized as bursting. Bottom: Membrane potential of a representative bursting neuron, corresponding to neuron number 6 in the raster plot. Initialization parameters of both are found in Table 1. The external input current is represented by the blue steps in all plots.

Figure 2:

Spiking and membrane potentials of sample neurons. (a) Top: Raster plot of spiking activity for neurons initialized as tonically firing. Bottom: Membrane potential of a representative tonic neuron, corresponding to neuron number 6 in the raster plot. (B) Top: Raster plot for neurons initialized as bursting. Bottom: Membrane potential of a representative bursting neuron, corresponding to neuron number 6 in the raster plot. Initialization parameters of both are found in Table 1. The external input current is represented by the blue steps in all plots.

Close modal

The plots show that using a normal distribution to initialize some parameters allows for individual neurons to start firing at different times and fire at different frequencies while keeping their overall tonic or bursting behavior. The external current ramp confirms that tonic firing of most neurons begins around 250 pA and around 110 pA for bursting.

The hybrid CPG architecture is constructed in Figure 1B with the following parameter values: P1: 5:1 excitatory/inhibitory ratio; P2: 3% (RG), 9% (AInh) connectivity; P3: 30% bursting neurons; P4: balanced with a slight excitatory bias (average of RGs: 9.79%); P5: 320 pA (tonic), 160 pA (bursting) current noise (see Table 3). The parameter values are manually tuned until an antagonistic output is produced from a network with sparse firing (see Figure 3, top row).

Figure 3:

Hybrid network with (top row) and without (bottom row) inhibition within the RGs and consequent sparse versus dense activity. (A) Schematics showing the configuration of the network. (B) Normalized population firing rate output per RG showing antiphasic population firing from both network architectures. At each time step, the number of spikes within a 5 ms window is counted, producing an analog output signal from each RG. (C) The plot compares the number of spikes on the x-axis to the number of neurons on the y-axis. The peak of the plot shows a bias toward zero in the network containing sparsely firing neurons. (D) Neuronal phase sorted firing rate per RG. The output of RG1 is on top and RG2 is below. The color represents the firing rate and is normalized between 0 and 1, where 0 denotes no firing. The antagonistic inhibitory interneuron population (AInh) firing is not included. (E) Neuronal phase distribution histogram using all populations of the hybrid CPG model, including the AInh populations. Each bin contains 0.25 radians. Plots D and E indicate that all neuronal phases are represented when sparsely firing neurons are present. (F) PCA of the firing rate activity from all populations, including the AInh populations. The network allowing sparse firing shows a circular shape in comparison to the U-shape of the abundantly firing network (seed = 76,25,853).

Figure 3:

Hybrid network with (top row) and without (bottom row) inhibition within the RGs and consequent sparse versus dense activity. (A) Schematics showing the configuration of the network. (B) Normalized population firing rate output per RG showing antiphasic population firing from both network architectures. At each time step, the number of spikes within a 5 ms window is counted, producing an analog output signal from each RG. (C) The plot compares the number of spikes on the x-axis to the number of neurons on the y-axis. The peak of the plot shows a bias toward zero in the network containing sparsely firing neurons. (D) Neuronal phase sorted firing rate per RG. The output of RG1 is on top and RG2 is below. The color represents the firing rate and is normalized between 0 and 1, where 0 denotes no firing. The antagonistic inhibitory interneuron population (AInh) firing is not included. (E) Neuronal phase distribution histogram using all populations of the hybrid CPG model, including the AInh populations. Each bin contains 0.25 radians. Plots D and E indicate that all neuronal phases are represented when sparsely firing neurons are present. (F) PCA of the firing rate activity from all populations, including the AInh populations. The network allowing sparse firing shows a circular shape in comparison to the U-shape of the abundantly firing network (seed = 76,25,853).

Close modal

The interburst interval (see Figure 3B) and firing sparsity (see Figure 3C) validate the network with sparsely firing neurons, producing the desired output. The population firing rate shows that the RG populations are firing in an antiphasic pattern to allow for control of an antagonistic muscle pair. This is evidenced by comparing the local maximums from each population and ensuring they are antiphasic (see Figure 3B). When equation 2.7 is used, the average population phase difference is determined to be 188.44°. The average frequency of the output from the RGs is 2.08 Hz. The firing is sparse across the population, and the distribution of spike count is low and biased toward zero (see Figure 3C). Dividing the number of sparsely firing neurons by the total number of firing neurons reveals that more than one-fifth of the neurons (24.36%) fired twice or less per cycle. Importantly, some neurons only fire once or twice throughout the whole simulation, in accordance with the experimental observation of firing rate distribution, which is skewed toward zero (Berg, 2017; Lindén & Berg, 2021; Petersen & Berg, 2016). The presented results indicate that the manually tuned network produces antagonistic output while maintaining sparse firing. Furthermore, we observe that the neuronal phase-sorted plot of firing rates (see Figure 3D) along with the neuronal phase distribution plot (E) have a somewhat wide distribution. The low-dimensional activity is revealed by a simple structure of the population firing rate trajectory, illustrated by the principal component analysis (PCA) in time (F). Even though the neuronal phase distribution has a distinct peak at π and 2π (see Figure 3E) and limited representation at intermediate neuronal phases, it is not as limited in neuronal phase as the abundantly firing network. When removing inhibition within the RG, thus increasing firing and removing balance, the distinct mode in the neuronal phase distribution is distinctly pronounced around π and 2π (see Figure 3, bottom row). The individual neuronal firing rates for all neurons are compared to the output from RG1; therefore, neurons in RG1 will fire mostly in-phase (0 and 2π), and neurons in RG2 will fire mostly in antiphase (π).

The “unbalance” of excitation and inhibition within the RGs resulted in a clear alternation between the two RG populations (see Figures 3B to 3D), and a noncircular, U-shaped PCA representation (F). The firing rates across the population were very high with 0% sparsely firing neurons, and the distribution had a clear mode (see Figure 3E) due to the removal of inhibition and the uncontrolled escalation of activity, as described previously (Lindén & Berg, 2021). In such an unbalanced RG network, all the neurons of a particular population tend to fire together; hence, the neuronal phase of their activity peaks around the same time during each period and the distribution has strong peaks at zero/2π and π (see Figure 3E). These results show that while alternating activity is still produced by unbalanced abundantly firing networks, fewer neuronal phases are represented.

3.2  Parameter Testing

Table 4 shows the values or ranges per parameter producing alternating output with sparse firing.

Table 4:

Range of Select Parameters Able to Maintain Antagonistic Output with Sparse Firing.

IDDescriptionValue Range
P1 Ratio of excitatory versus inhibitory neurons in RG 4:1 or 5:1 
P2 Connectivity (RG) 3% 
P2 Connectivity (AInh) 9% 
P3 Neuronal subtype (bursting) Minimum of 25% 
P4 RG balance Balanced with a slight bias toward excitation (average of RGs: 9.5%) 
P4 Network balance Inhibitory bias 
P5 Noise 90% to 500% 
IDDescriptionValue Range
P1 Ratio of excitatory versus inhibitory neurons in RG 4:1 or 5:1 
P2 Connectivity (RG) 3% 
P2 Connectivity (AInh) 9% 
P3 Neuronal subtype (bursting) Minimum of 25% 
P4 RG balance Balanced with a slight bias toward excitation (average of RGs: 9.5%) 
P4 Network balance Inhibitory bias 
P5 Noise 90% to 500% 

Note: The corresponding plots are found in the supplementary material (Figures 4–8) and specific figure references can be found in the text following the table.

Additionally, running the simulation with the control parameters but initializing with a new random seed gives an average frequency over 20 trials of 1.92 Hz, an average population phase difference between signals of 174.69° and an average firing sparsity of 23.69%. Based on the collected data, the network can be described as robust because even with a noise-to-signal ratio of 5:1, some sparse firing (4.18%) and alternating activity (191.73°) are still recorded (see Figure 8, top row, in the supplementary material). In other words, the network still produces desired results when the standard deviation of the noise parameter is five times larger than the input bias current. The output of the network with 6:1 noise depicts degradation in population firing rate output (see Figure 8, bottom row, in the supplementary material). It is important to note that the network requires 90% noise to produce enough random spiking to promote continuous oscillations throughout the simulation. Furthermore, it can be concluded that the RG populations must be balanced with a slight excitatory bias to produce an alternating output for the duration of the simulation. If the individual RG populations have an inhibitory bias (see Figure 7, first row, in the supplementary material), the network does not have persistent oscillating output. Alternatively, if the RG populations are too excitatory in their bias, the neural activity becomes alternating, showing neuronal phase representation at close to zero/2π and π, and the spiking per neuron increases (see Figure 7, second row, in the supplementary material). When inhibition is removed from the individual RG populations, there is an anti-phasic population output from the network but neuron spiking increases and neuronal phase representation is not distributed (see Figure 3, bottom row). Furthermore, when the mutual inhibition is removed so the RG populations fire as individual populations, output from each RG is independent so the population phase is not consistent. There is an increase in output frequency and neuron spiking as well as a preference for neuronal in-phase firing (see Figure 7, fourth row, in the supplementary material). In-phase firing is determined by a large single peak at 0/2π radians in the neuronal phase distribution histogram. Finally, when all inhibition is removed from the network, neurons fire abundantly and fire in-phase within each RG without a consistent population phase difference between RG populations (see Figure 7, third row, in the supplementary material). Additionally, the connectivity must be weak as decreasing synaptic sparsity forces neurons to spike at the same time, reducing neuronal phase representation (see Figure 5, middle row, in the supplementary material) and eventually leads to oversuppression of network output (see Figure 5, bottom row, in the supplementary material). However, if connectivity is too weak, firing will be random (see Figure 5, top row, in the supplementary material). Finally, the percentage of bursting neurons must be at least 25% to produce continuous oscillations (see Figure 6, top row, in the supplementary material). Based on biological findings, the percentage of neurons that are bursting varies based on extracellular ion levels (Brocard et al., 2013), so we attempt to find the lowest percentage of bursting neurons for our network that still provides alternating output.

The connectivity of the network is shown to have a significant impact on the firing dynamics of the population. If there is strong connectivity, oversuppression can occur (see Figure 5, bottom row, in the supplementary material) whereas weak connectivity can devolve into random firing (see Figure 5, top row, in the supplementary material), which does not seem to depend on other neurons in the network. Even a smaller increase in connectivity can change the behavior of the network, where the RGs can still produce antagonistic output and sparse firing is recorded but neuronal phase representation is not broad. This is due to more neurons firing together and the sparse neurons also firing within these bursts of activity. The initial test for 5% connectivity within the RG and 15% connectivity from the inhibitory populations shows an overrepresentation of sparsely firing neurons (see Figure 5, middle row, in the supplementary material). When the test is rerun with a different seed, to create different random connections and noise, the sparse firing is reduced to 8.91% of neurons, a significant reduction from the original trial showing 39.93%. Based on these observations, the connectivity must be carefully considered when implementing a sparsely spiking CPG. Optimization methods could be used to tune the connectivity parameters to different desired outputs. During tuning, the synaptic weights were also observed to have a significant impact on the output of the network. Increasing the mean synaptic weight from 0.6 nS to 1.2 nS was enough to reduce firing sparsity and promote neurons firing together.

The specific connectivity of the network is also shown to affect how quickly alternating activity emerges in the population dynamics. The random seed trials show different initialization times (see Figure 9, in the supplementary material). This reinforces the expectation that the network is producing alternating output and it is not a product of neurons being initialized in a specific way. For example, if all neurons were given the same characteristics at start-up, they would begin spiking at the same time. On the other hand, some of our tests show that firing is random at first, and as the neurons are allowed to fire, oscillations can emerge based on the connectivity within the network.

The interplay of excitation and inhibition in the network also proves to be a sensitive parameter. The postsynaptic connections from the AInh populations were less sparse and stronger than those from the RG populations in order to promote antiphasic population output. The AInh populations do not have any self-recurrent connections so all inhibitory connections were with neurons in the RG populations. The balance within the RG populations had to be balanced with a slight excitatory bias to ensure the output did not devolve into random firing or complete suppression. This means that the sum of the weights of the excitatory synapses within the RG was larger than that of the inhibitory synapses (as defined by equations 2.5 and 2.6). If the RG was inhibitory, the population firing rate output was noisy (see Figure 7B, first row, in the supplementary material). However, inhibition is required within the network; otherwise, the population output from the RG populations will not be antiphasic (see Figure 7, third row, in the supplementary material). When inhibition is removed from different parts of the network, there is a strong preference for neuronal in-phase firing, and spiking activity increases significantly (see Figure 3, bottom row, and Figure 7, third and fourth row, in the supplementary material). This is an intuitive result and confirms the important role of inhibition within the network (Lindén & Berg, 2021; Petersen & Berg, 2016). It is important to note that our calculation of balance is a simplification, using the synaptic conductances and spikes traveling across each synapse to approximate the amount of excitation versus inhibition within a simulation run. A more accurate calculation of balance would include the driving force of the synaptic currents when an individual spike is sent.

The network shows that a minimum of 25% bursting neurons within the RGs is necessary to drive continuous oscillations. Bursting neurons reinforce rhythm by increasing the number of spikes within specific time intervals. This study cannot conclude whether it is the percentage of bursting neurons or the number of spikes within a time interval that is necessary to reinforce output oscillations. However, we assert that bursting neurons are a biologically plausible method to produce a minimum threshold of spikes at regular intervals.

It is also notable that running the RGs independently by removing mutual inhibition in the network still produces a strong preference for neuronal in-phase firing (see Figure 7, fourth row, in the supplementary material). This is not the expected result as the RGs should show nearly equal neuronal phase representation when operating by themselves as found in Lindén et al. (2022). This could be an artifact of the firing behavior of the spiking neuron model since neurons are initialized with similar parameters coupled with the self-excitation of the excitatory neurons, which reinforces neuronal in-phase firing. Additionally, spiking neurons have a refractory period that reduces the ability of the neuron to fire at specific times, namely, during the refractory period.

The presented results show that random firing of neurons causes the lowest average neuronal phase distribution difference (for an example, see Figure 5E, bottom row, in the supplementary material). This makes intuitive sense as the neurons would have no preference for firing at any particular time. Therefore, this must be accounted for when developing a quantitative measure of neuronal phase distribution.

The neuronal phase distribution plots (see Figure 3E) show that the most represented phases are at 0/2π and π. This means that while neurons fire at all neuronal phases, they are more likely to fire in-phase or antiphasic to each other. This can be accounted for by the configuration of the RG populations, which are set up as RG and AInh populations. Within the RGs, the excitatory neurons will fire close to in-phase with each other due to recurrent excitatory connections. When the inhibitory neurons are firing, they will suppress excitatory neuron output, leading to fewer neurons firing between neuronal phases.

The test results for 500% noise in the network appear to contradict our conclusion that sparse firing enables broad neuronal phase representation in a spiking network (see Figure 8, top row, in the supplementary material). The population firing rate output is alternating, and the neuronal phase distribution histogram shows firing in all phases, but only 4.18% of neurons are sparsely firing. This indicates that a high amount of noise in the network increases firing per neuron. However, spikes that are elicited by peaks in the noise do not arise from network connectivity and therefore are not correlated with network activity. These noise-induced spikes thus likely do not contribute to an increase in alternation between the RGs or in-phase firing within each RG. Therefore, we maintain that this test is still consistent with our conclusion that sparse firing is one of the factors that allows for a broad neuronal phase representation of neuron firing.

In the network where all inhibition is removed (see Figure 7, third row, in the supplementary material), the PCA for each RG is still circular, although oval, which taken alone would indicate rotational activity and an expectation of broad neuronal phase representation. However, the neuronal phase distribution plot shows that only phases close to 0/2 π are present. This illustrates that circular trajectory in PCA space is necessary but not sufficient for the network to exhibit rotational activity as defined in Lindén et al. (2022). Instead, the PCA plot should be considered together with the neuronal phase distribution. Due to the method of relating data using orthogonal vectors, a PCA can show circular behavior with as little as two neuronal phases represented.

Finally, it is important to note that the simulation shows that some neurons are silent and do not fire at all. These silent neurons were removed from the firing rate neuronal phase plots. We argue that silent neurons are biological as they have been observed in fictive locomotion of the mouse spinal cord where some neurons stop firing at lower frequencies (Zhong et al., 2011) and in zebrafish where different motor neurons are recruited at increasing speeds of motion (Jha & Thirumalai, 2020). Furthermore, neurons are known to be multistable; they are able to switch between tonic, bursting, and silence depending on the level of hyperpolarization (Malashchenko et al., 2011). Therefore, the observation of silent neurons during slower frequency firing is predicted by biological studies.

A hybrid CPG network was developed and confirmed to produce antiphasic population output in the presence of sparse firing. Furthermore, it was observed that sparse firing allowed broad neuronal phase representation of neurons firing at a population level, which brings the model closer to exhibiting rotational dynamics as observed in experiments (Lindén et al., 2022). Testing the network by changing parameters independently and comparing the output to a control experiment produced the following observations:

  • The network is robust against noise.

  • The RG populations should be balanced with a slight excitatory bias to maintain oscillatory output.

  • The network should have weak connectivity (3% RG, 9% AInh) to promote coordinated output while allowing for a broad neuronal firing phase distribution.

The implementation of this network, which increases biological fidelity by reproducing a phenomenon observed in biological studies, namely, sparse firing, paves the way for further investigation into rhythmic movement. The next step will be testing the model’s behavior for different cycle frequencies and interfacing the computational model with a musculoskeletal model that provides feedback in order to adapt the network output.

B.S. and E.N. developed the main idea of the article and contributed equally to this work. E.N. implemented the spiking rhythm generator (RG) in Brian and the data analysis tools. B.S. researched the biological neural network concepts, adapted the spiking RG into a neural simulation tool (NEST), and implemented the central pattern generator architecture. The manuscript was written by B.S. and E.N. with support from J.A., R.B., and S.T., J.A., and R.B provided theoretical and experimental neuroscientific knowledge. J.A., R.B., and S.T. discussed the results. All authors contributed to the article and approved the submitted version.

We thank Simon Danner for his valuable input during the model development process. This work was funded by the Lundbeck Foundation, grant no. R370-2021-948, and the National Science Foundation, NSF CRCNS/DARE 2113069. We declare that there are no competing interests.

Ausborn
,
J.
,
Shevtsova
,
N. A.
, &
Danner
,
S. M
. (
2021
).
Computational modeling of spinal locomotor circuitry in the age of molecular genetics
.
International Journal of Molecular Sciences
,
22
(
13
),
6835
.
Ausborn
,
J.
,
Snyder
,
A. C.
,
Shevtsova
,
N. A.
,
Rybak
,
I. A.
, &
Rubin
,
J. E
. (
2018
).
State-dependent rhythmogenesis and frequency control in a half-center locomotor CPG
.
Journal of Neurophysiology
,
119
(
1
),
96
117
.
Berg
,
R. W
. (
2017
).
Neuronal population activity in spinal motor circuits: Greater than the sum of its parts
.
Frontiers in Neural Circuits
,
11
,
103
.
Berg
,
R. W.
,
Willumsen
,
A.
, &
Lindén
,
H
. (
2019
).
When networks walk a fine line: Balance of excitation and inhibition in spinal motor circuits
.
Current Opinion in Physiology
,
8
,
76
83
.
Brocard
,
F.
,
Shevtsova
,
N. A.
,
Bouhadfane
,
M.
,
Tazerart
,
S.
,
Heinemann
,
U.
,
Rybak
,
I. A.
, &
Vinay
,
L
. (
2013
).
Activity-dependent changes in extracellular Ca2+ and K+ reveal pacemakers in the spinal locomotor-related network
.
Neuron
,
77
(
6
),
1047
1054
.
Churchland
,
M. M.
,
Cunningham
,
J. P.
,
Kaufman
,
M. T.
,
Foster
,
J. D.
,
Nuyujukian
,
P.
,
Ryu
,
S. I.
, &
Shenoy
,
K. V
. (
2012
).
Neural population dynamics during reaching
.
Nature
,
487
(
7405
),
51
56
.
Danner
,
S. M.
,
Hofstoetter
,
U. S.
,
Freundl
,
B.
,
Binder
,
H.
,
Mayr
,
W.
,
Rattay
,
F.
, &
Minassian
,
K
. (
2015
).
Human spinal locomotor control is based on flexibly organized burst generators
.
Brain
,
138
(
3
),
577
588
.
Dougherty
,
K. J.
, &
Ha
,
N. T
. (
2019
).
The rhythm section: An update on spinal interneurons setting the beat for mammalian locomotion
.
Current Opinion in Physiology
,
8
,
84
93
.
Dougherty
,
K. J.
,
Zagoraiou
,
L.
,
Satoh
,
D.
,
Rozani
,
I.
,
Doobar
,
S.
,
Arber
,
S.
, . . .
Kiehn
,
O
. (
2013
).
Locomotor rhythm generation linked to the output of spinal Shox2 excitatory interneurons
.
Neuron
,
80
(
4
),
920
933
.
Graham Brown
,
T
. (
1911
).
The intrinsic factors in the act of progression in the mammal
.
Proceedings of the Royal Society of London. Series B, Containing Papers of a Biological Character
,
84
(
572
),
308
319
.
Graham Brown
,
T
. (
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
.
Journal of Physiology
,
48
(
1
),
18
.
Grillner
,
S.
, &
El Manira
,
A
. (
2020
).
Physiological reviews review article
.
Physiological Reviews
,
100
,
271
320
.
Grillner
,
S.
, &
Kozlov
,
A
. (
2021
).
The CPGs for limbed locomotion–facts and fiction
.
International Journal of Molecular Sciences
,
22
(
11
),
5882
.
Hart
,
C. B.
, &
Giszter
,
S. F
. (
2010
).
A neural basis for motor primitives in the spinal cord
.
Journal of Neuroscience
,
30
(
4
),
1322
1336
.
Jankowska
,
E.
,
Jukes
,
M. G.
,
Lund
,
S.
, &
Lundberg
,
A
. (
1967
).
Half-centre organization of interneurones transmitting effects from the flexor reflex afferents
.
Acta Physiologica Scandinavica
,
70
(
3
),
389
402
.
Jha
,
U.
, &
Thirumalai
,
V
. (
2020
).
Neuromodulatory selection of motor neuron recruitment patterns in a visuomotor behavior increases speed
.
Current Biology
,
30
(
5
),
788
801
.
Kalidindi
,
H. T.
,
Cross
,
K. P.
,
Lillicrap
,
T. P.
,
Omrani
,
M.
,
Falotico
,
E.
,
Sabes
,
P. N.
, &
Scott
,
S. H
. (
2021
).
Rotational dynamics in motor cortex are consistent with a feedback controller
.
eLife
,
10
,
e67256
.
Kiehn
,
O
. (
2016
).
Decoding the organization of spinal circuits that control locomotion
.
Nature Reviews Neuroscience
,
17
(
4
),
224
238
.
Kozlov
,
A. K.
,
Lansner
,
A.
,
Grillner
,
S.
, &
Kotaleski
,
J. H
. (
2007
).
A hemicord locomotor network of excitatory interneurons: A simulation study
.
Biological Cybernetics
,
96
,
229
243
.
Lindén
,
H.
, &
Berg
,
R. W
. (
2021
).
Why firing rate distributions are important for understanding spinal central pattern generators
.
Frontiers in Human Neuroscience
,
15
,
719388
.
Lindén
,
H.
,
Petersen
,
P. C.
,
Vestergaard
,
M.
, &
Berg
,
R. W.
(
2022
).
Movement is governed by rotational neural dynamics in spinal motor networks
.
Nature
,
610
(
7932
),
526
531
.
Malashchenko
,
T.
,
Shilnikov
,
A.
, &
Cymbalyuk
,
G
. (
2011
).
Six types of multistability in a neuronal model based on slow calcium current
.
PLOS One
,
6
(
7
),
e21782
.
Marder
,
E.
, &
Bucher
,
D
. (
2001
).
Central pattern generators and the control of rhythmic movements
.
Current Biology
,
11
(
23
),
R986–R996
.
Musienko
,
P. E.
,
Lyalka
,
V. F.
,
Gorskii
,
O. V.
,
Zelenin
,
P. V.
, &
Deliagina
,
T. G.
(
2022
).
Activity of spinal interneurons during forward and backward locomotion
.
Journal of Neuroscience
,
42
(
17
),
3570
3586
.
Naud
,
R.
,
Marcille
,
N.
,
Clopath
,
C.
, &
Gerstner
,
W
. (
2008
).
Firing patterns in the adaptive exponential integrate-and-fire model
.
Biological Cybernetics
,
99
(
4–5
),
335
.
Petersen
,
P. C.
, &
Berg
,
R. W
. (
2016
).
Lognormal firing rate distribution reveals prominent fluctuation–driven regime in spinal motor networks
.
eLife
,
5
.
Petersen
,
P. C.
,
Vestergaard
,
M.
,
Jensen
,
K. H.
, &
Berg
,
R. W
. (
2014
).
Premotor spinal network with balanced excitation and inhibition during motor patterns has high resilience to structural division
.
Journal of Neuroscience
,
34
(
8
),
2774
2784
.
Radosevic
,
M.
,
Willumsen
,
A.
,
Petersen
,
P. C.
,
Lindén
,
H.
,
Vestergaard
,
M.
, &
Berg
,
R. W
. (
2019
).
Decoupling of timescales reveals sparse convergent CPG network in the adult spinal cord
.
Nature Communications
,
10
(
1
),
2937
.
Ramirez
,
J.-M.
, &
Bush
,
N.
(
2022
).
Latent neural population dynamics underlie normal and disrupted breathing
.
bioRxiv
.
Rancic
,
V.
, &
Gosgnach
,
S
. (
2021
).
Recent insights into the rhythmogenic core of the locomotor CPG
.
International Journal of Molecular Sciences
,
22
(
3
),
1394
.
Rybak
,
I. A.
,
Dougherty
,
K. J.
, &
Shevtsova
,
N. A.
(
2015
).
Organization of the mammalian locomotor CPG: Review of computational model and circuit architectures based on genetically identified spinal interneurons
.
eNeuro
,
2
(
5
).
Shevtsova
,
N. A.
,
Ha
,
N. T.
,
Rybak
,
I. A.
, &
Dougherty
,
K. J
. (
2020
).
Neural interactions in developing rhythmogenic spinal networks: Insights from computational modeling
.
Frontiers in Neural Circuits
,
14
,
614615
.
Sinha
,
A.
,
de Schepper
,
R.
,
Pronold
,
J.
,
Mitchell
,
J.
,
Mørk
,
H.
,
Nagendra Babu
,
P.
, . . .
Plesser
,
H. E.
(
2023
,
February
).
Nest 3.4.
Zenodo
.
Song
,
J.
,
Pallucchi
,
I.
,
Ausborn
,
J.
,
Ampatzis
,
K.
,
Bertuzzi
,
M.
,
Fontanel
,
P.
, . . .
El Manira
,
A
. (
2020
).
Multiple rhythm-generating circuits act in tandem with pacemaker properties to control the start and speed of locomotion
.
Neuron
,
105
(
6
),
1048
1061
.
Strohmer
,
B.
(
2023
).
Hybrid CPG source code.
https://github.com/bstrohmer/hybrid-cpg
Sukenik
,
N.
,
Vinogradov
,
O.
,
Weinreb
,
E.
,
Segal
,
M.
,
Levina
,
A.
, &
Moses
,
E
. (
2021
).
Neuronal circuits overcome imbalance in excitation and inhibition by adjusting connection numbers
.
Proceedings of the National Academy of Sciences
,
118
(
12
),
e2018459118
.
Wilson
,
J. M.
,
Hartley
,
R.
,
Maxwell
,
D. J.
,
Todd
,
A. J.
,
Lieberam
,
I.
,
Kaltschmidt
,
J. A.
, . . .
Brownstone
,
R. M
. (
2005
).
Conditional rhythmicity of ventral spinal interneurons defined by expression of the Hb9 homeodomain protein
.
Journal of Neuroscience
,
25
(
24
),
5710
5719
.
Zhong
,
G.
,
Sharma
,
K.
, &
Harris-Warrick
,
R. M
. (
2011
).
Frequency-dependent recruitment of V2a interneurons during fictive locomotion in the mouse spinal cord
.
Nature Communications
,
2
(
1
),
274
.
Zhou
,
S.
, &
Yu
,
Y
. (
2018
).
Synaptic E-I balance underlies efficient neural coding
.
Frontiers in Neuroscience
,
12
,
46
.

Author notes

Beck Strohmer and Elias Najarro contributed equally to this work and share first authorship.

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. For a full description of the license, please visit https://creativecommons.org/licenses/by/4.0/legalcode

Supplementary data