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 rhythmgenerating 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, 1911Brown, , 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 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.

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(Ausborn et al., , 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).

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.

Generator Network. The central pattern generator network (CPG) network architecture is based on the flexor/extensor halfcenter RG architecture (Ausborn et al., 2021(Ausborn et al., , 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.

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

In Table 1, V th is the voltage threshold potential, t re f is the refractory period, C is the membrane capacitance, I bias is the bias current, V m is the membrane potential, g L is the leak conductance, E L 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 V reset 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.

Equations 2.1 and 2.2 describe the AdEx neuron's dynamics:

The reset equation for each is

(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 (g e (t)(V m -E e )) and inhibitory (g i (t)(V m -E i )) synaptic conductances into equation 2.1. The term g e (t) is the excitatory conductance, E e is the excitatory reversal potential, g i (t) is the inhibitory conductance, and E i 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:

where W e is the excitatory weight, t is the time of the presynaptic spike, τ syn_ex is the excitatory rise time, W i 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_weights total_weight ,

(2.5)

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 excita-tory_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.

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,

where avg_phase_diff is the average population phase difference between the two output signals, peaks rg1 represent the local maximums of the first output signal, peaks rg2 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:

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.

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.

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.

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

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 lowdimensional 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. 4 shows the values or ranges per parameter producing alternating output with sparse firing.

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.