Neural Interactions in Developing Rhythmogenic Spinal Networks: Insights From Computational Modeling

The mechanisms involved in generation of rhythmic locomotor activity in the mammalian spinal cord remain poorly understood. These mechanisms supposedly rely on both intrinsic properties of constituting neurons and interactions between them. A subset of Shox2 neurons was suggested to contribute to generation of spinal locomotor activity, but the possible cellular basis for rhythmic bursting in these neurons remains unknown. Ha and Dougherty (2018) recently revealed the presence of bidirectional electrical coupling between Shox2 neurons in neonatal spinal cords, which can be critically involved in neuronal synchronization and generation of populational bursting. Gap junctional connections found between functionally-related Shox2 interneurons decrease with age, possibly being replaced by increasing interactions through chemical synapses. Here, we developed a computational model of a heterogeneous population of neurons sparsely connected by electrical or/and chemical synapses and investigated the dependence of frequency of populational bursting on the type and strength of neuronal interconnections. The model proposes a mechanistic explanation that can account for the emergence of a synchronized rhythmic activity in the neuronal population and provides insights into the possible role of gap junctional coupling between Shox2 neurons in the spinal mechanisms for locomotor rhythm generation.


INTRODUCTION
The mammalian spinal cord contains neuronal circuits that can generate locomotor-like oscillations in the absence of supraspinal and afferent inputs. These circuits include rhythm-generating kernels capable of producing, maintaining, and coordinating populational rhythmic activity. A rhythmogenic kernel on each side of the cord is thought to be comprised of mutually connected excitatory neurons with ipsilateral projections (Kjaerulff and Kiehn, 1996;Kiehn, 2006;Hägglund et al., 2010). Although no single genetically identified neuron type has been found to be solely responsible for rhythm generation in the spinal cord, the genetically identified Shox2 interneurons have been suggested to be involved in locomotor rhythmogenesis (Dougherty et al., 2013). Later investigations provided evidence that a subpopulation of Shox2 neurons, specifically the Shox2 non-V2a neuron type, may belong to the rhythm-generating kernel in the neonatal rodent spinal cord and be involved in generation of locomotor rhythmic activity (Dougherty et al., 2013;Kiehn, 2016).
Regardless of the intrinsic bursting properties of single neurons involved in rhythmogenesis, mutual interactions between the neurons are critical to allow synchronization of neuronal activity leading to the generation of populational rhythmic activity. These excitatory interactions in the spinal cord are mediated by both chemical and electrical synapses (Tresch and Kiehn, 2000), the importance and role of which appear to be age-dependent (Walton and Navarrete, 1991;Chang et al., 1999). A recent study has shown that functionally related Shox2 interneurons in neonatal mice are locally connected with each other bidirectionally via electrical synapses or gap junctions (Ha and Dougherty, 2018). This bidirectional electrical coupling may represent a mechanism for neuronal synchronization and for populational oscillations in the neonatal spinal cord. Bidirectional gap junctional coupling within each population of genetically identified neurons, i.e., Shox2 non-V2a neurons, is found with high incidence in the neonatal spinal cord. With age, this electrical coupling decreases in both incidence and strength, which may occur in parallel (and possibly in connection) with an increase of the unilateral connections via excitatory chemical synapses (Ha and Dougherty, 2018). This change in the type, probability and strength of neuronal interactions may lead to specific changes in some characteristics of locomotor rhythm and pattern and the range of generated frequencies.
To theoretically investigate the potential roles of different (electrical and excitatory chemical) connections between Shox2 interneurons and their interaction with I NaP in generation of rhythmic populational activity, we developed a model of a population of neurons incorporating I NaP in a subset of cells that are sparsely connected by electrical and/or excitatory chemical synapses and investigated the dependence of frequency of populational bursting on the connection type and strength, the maximal conductance and distribution of I NaP , and other neuronal characteristics. We demonstrated that the frequency of populational bursting activity depends on the relative expression of particular types of connections between neurons. The model proposes a mechanistic explanation for emergence of a synchronized rhythmic activity in a population of Shox2 neurons and provides insights into the mechanisms of the locomotor rhythm generation.

Experimental Methods
All experimental procedures followed NIH guidelines and were reviewed and approved by the Institutional Animal Care and Use Committee at Drexel University.
Here, 25 bursts immediately before carbenoxolone application and 30 min after the start of carbenoxolone application were subjected to analysis previously described (Dougherty et al., 2013). Briefly, ventral root activity was rectified and smoothed (time constant 0.2 s) in Spike2 (Cambridge Electronic Design). Recordings included either 2 ventral roots or 3 ventral roots. The cycle and amplitude values were taken from a rostral lumbar root (L1, L2, or L3) ipsilateral to a recorded caudal lumbar root (L4 or L5). The frequency was calculated from the mean cycle period, each measured from burst onset to onset. Amplitude was the measured peak amplitude of each burst. The mean values were normalized to the no drug condition to control for interexperiment variation including root size, suction electrode size and seal.

Neuron Model
All neurons were simulated in the Hodgkin-Huxley style as single-compartment neuron models. The membrane potential, V, in neurons with I NaP was described by the following differential equation: where C is the membrane capacitance and t is time.
The membrane potential, V, in neurons without I NaP was described as follows: The ionic currents in Equations (1) and (2) were described as follows: where I Na is the fast sodium current with maximal conductancē g Na ; I NaP is the persistent (slowly inactivating) sodium current with maximal conductanceḡ NaP ; I K is the delayed-rectifier potassium current with maximal conductanceḡ K ; I L is the leakage current with constant conductance g L . E Na , E K , and E L are the reversal potentials for sodium, potassium and leakage currents, respectively; m and h with indexes indicating currents are, respectively, the activation and inactivation variables of the corresponding ionic channels. Dynamics of activation and inactivation variables of voltagedependent ionic channels (Na, NaP, and K) in Equation (3) was described by the following differential equations: where m ∞i (V) and h ∞i (V) define the voltage-dependent steadystate activation and inactivation of the channel i, respectively, and τ mi (V) and τ hi (V) define the corresponding time constants. Activation of the fast and persistent sodium channels was instantaneous. The expressions for channel kinetics in Equation (4) are described as follows: Parameters for channel kinetics were derived from Brocard et al. (2013). The maximal conductances for the fast sodium, delayedrectifier potassium, leakage currents were, respectively:ḡ Na = 80 nS,ḡ K = 100 nS, and g L = 1 nS. The maximal conductance and inactivation time constant for I NaP ,ḡ NaP and T h max , and the leakage reversal potential, E L , in different simulation were varied. The synaptic current (I Syn with conductance g Syn and reversal potential E Syn ) for the j-th neuron was described as follows: ..., N, (6) where g Syn = 1 nS, E Syn = 0, and w Syn,ij = w Syn if the i-th and jneurons were synaptically coupled and equal to 0 otherwise. The value of w Syn was varied.
The synaptic variable s was governed by first order activation scheme (Compte et al., 2003): where α s = 1; τ s = 15 ms; V s12 = −20; and k s = 2. The electrical coupling current I Gap for the j-th neuron was described as follows: where the coupling strength g Gap,ij = g Gap if the i-th and j-th neurons share a gap junction and is equal to 0 otherwise. The value of g Gap was varied.
The following general neuronal parameters were assigned: C = 40 pF; E Na = 55 mV; E K = -80 mV.

Two-Cell Model
Each neuron in the two-cell model represented either a bursting neuron incorporating the persistent sodium current, I NaP , or simple tonically spiking neuron without I NaP . Dynamics of the membrane potential of the neurons was described by Equations (1)-(5) whereḡ NaP = 4 nS, if not otherwise stated. These neurons were coupled either bidirectionally by electrical synapses as described by Equation (8) (6) and (7), w Syn = {0.5, 1, 2}]. The initial operational regime of each neuron was defined by the value of its leakage reversal potential, E L , indicated in Tables 1, 2.

Neuron Population
The neuron population contained 100 neurons and included two subpopulations: neurons incorporating I NaP (subpopulation 1) and neurons non-incorporating I NaP (subpopulation 2). In our simulations, N NaP , number of neurons incorporating I NaP, varied from 10 to 40 and the observed results were qualitatively similar.
Heterogeneity of neurons within the population was provided by random distributions of the baseline values of the mean conductance of the leakage channel, g L , and the leakage reversal potentials, E L0i (i = 1, 2) for two subpopulations, mean maximal conductance of the persistent sodium channel,ḡ NaP , in the subpopulation of neurons incorporating I NaP , and initial conditions for the values of membrane potential and channel kinetics variables. The baseline values of distributed parameters and initial conditions were assigned prior to simulations from their defined average values and variances using a random number generator (normal distribution), and a settling period of 10-40 s was allowed in each simulation. In our simulations we varied E L01 = {−73, −74, −75} for subpopulation of neurons with I NaP ; E L02 = {−68, −70, −72} for subpopulation of neurons without I NaP ; parameterḡ NaP was varied in the range of [2,5].  Figure 1. We also run simulation with different values of the maximal time constant of the persistent sodium current inactivation, T h max = {7, 9, 10} s. In simulations for the populational model shown in this paper: N NaP = 40; E L01 = −74 ± 14.8 mV; E L02 = −70 ± 14 mV; g L = 1 ± 0.2;ḡ NaP = 4 ± 0.8; T h max = 9 s, if not stated otherwise. Random connections between neurons in the population were assigned prior to each simulation based on probability of connections (p). A random number generator was used to define the existence of each synaptic connection. Typically, we used p Gap = 0.3 for gap junctional connections (Ha and Dougherty, 2018) and p Syn = 0.1 (Dougherty et al., 2013) for excitatory chemical connections. Similar probabilities of excitatory synaptic interconnections were estimated for neurons in the lumbar region of adult turtles (Radosevic et al., 2019) and in the pre-Bötzinger Complex involved in the generation of rhythmic inspiratory activity (Carroll and Ramirez, 2013). After the distribution of connections was set, the connection strengths of electrical and/or excitatory chemical synapses (g Gap and w Syn ) were altered. To estimate dependence of the model behavior on coupling strength and type when both types of connections were present on the model, we introduced two variables G and W (G = N × p Gap × g GaP and W = N × p Syn × w Syn ), characterizing the average total network electrical or chemical synaptic input to each neuron in the population. To define electrical and chemical synapses in the populational model we used standard definition of these synapses by Brain2 (Stimberg et al., 2014(Stimberg et al., , 2019. The coupling strengths for electrical and chemical connection (g Gap and w Syn or G and W) are specified for each simulation.

Computer Simulations
The simulations for the two-cell model were run by using a custom Matlab script (The Mathworks, Inc., Matlab 2020a). The simulations for the populational model were performed using Frontiers in Neural Circuits | www.frontiersin.org the custom script written in Python and implemented in Brian2 environment (Stimberg et al., 2014(Stimberg et al., , 2019. Differential equations were solved using the second order Runge-Kutta integration method. Simulation results were saved as the ASCI file containing the values of g Gap and w Syn and calculated mean values and standard deviation over simulation time for the period and amplitude of populational bursting. We also saved the figures showing populational activity, raster plot of neuron spiking, and voltage traces of sample neurons over simulation time to visually estimate the results of simulations. The scripts and model configuration files to create the simulations presented in the paper are available at https://github.com/RybakLab/Neuralinteractions-in-spinal-networks.

Data Analysis in Computer Simulations
The results of simulations were processed by custom Matlab scripts (The Mathworks, Inc., Matlab 2020a). For each neuron in the two-cell model, period of bursting was estimated as the difference between two consecutive burst onsets and averaged for the time course of simulation. Frequency was calculated as reciprocal to the period. To assess the populational model behavior, the averaged integrated activity of neurons in the population (average number of spikes per neuron per second) was used to calculate the oscillation period. The timing of burst onsets was determined at a threshold level equal to 30% of the average populational burst amplitude in the current simulation. For some simulation with less stable populational activity we determined the threshold level by visual inspection of the results. The bursting period was defined as the duration between two consecutive burst onsets. Duration of individual simulations depended on the values of parameters g Gap and w Syn , and to robustly estimate average value of the oscillation period, the first 10-20 transitional cycles were omitted to allow stabilization of model variables, and the values of the bursting period and amplitude were averaged for the next 10-20 consecutive cycles. To build 2D diagrams we estimated stability of the populational bursting activity in each simulation by calculating values of the period and amplitude of populational bursts and their standard deviations. All simulation results were divided in four groups: no populational activity, low amplitude (<10 spikes per neuron per second) or highly unstable populational activity (SD of calculated period ≥50%), stable populational bursting, and tonic populational activity. For the stable populational bursting, in each simulation we calculated amplitude of populational bursts as the mean of maximal values of the averaged integrated activity of neurons during the burst. The frequency of populational oscillations was determined as reciprocal to the mean bursting period. The calculated frequency and amplitude were colorcoded to build 2D diagrams in G x W parameter space. Reinitialization of the randomized parameters within the same ranges resulted in qualitatively similar simulation results.

RESULTS
The major focus of the present study was on the specific roles of, and possible cooperation between, neural interactions (through gap junction coupling and/or excitatory chemical synapses) and persistent (slowly inactivating) sodium current (I NaP ) in neuronal synchronization and generation of populational rhythmic activity related to the locomotor-like oscillations in the mammalian spinal cord. We started by simulating a pair of neurons operating in different regimes and connected by bidirectional gap junctions or unidirectional excitatory chemical synapses in order to investigate the effect of the type and strength of connections and the presence/absence of I NaP conductance on the activity of both neurons. Then, we considered a heterogeneous population of 100 neurons which incorporated a subpopulation of neurons containing I NaP and hence could intrinsically generate rhythmic bursting depending on neuronal excitability. The remaining neurons were simple spiking neurons. The neurons had parameters randomly distributed over the population and were sparsely connected by electrical and excitatory chemical synapses. The probabilities of connections between neurons were set based on the data from Shox2 neurons (Dougherty et al., 2013;Ha and Dougherty, 2018). Our modeling study focused on neuronal interactions, synchronization of neuronal activity, and control of frequency of populational bursting depending on the relative expression of different types of neuron connectivity and connection strength.

Gap Junctional Coupling in Two-Cell Model
We first simulated a pair of interneurons (In1 and In2) modeled in the Hodgkin-Huxley style and coupled bidirectionally by gap junctions with conductance g Gap. Each neuron represented either a conditionally bursting neuron with I NaP or a simple tonically spiking neuron without I NaP . The equations describing neuron dynamics and neuron parameters are specified in section Materials and Methods.
The baseline operating regime of each neuron (prior to coupling) was dependent on its resting membrane potential, which in turn was defined by the value of the leakage reversal potential, E L (see Table 1). With an increase in E L, the isolated conditionally bursting neuron changed its operating regime from silence to bursting and then to tonic spiking. The isolated "simple" neuron without I NaP could switch from silence to tonic spiking when E L increased. The strength of gap junctional coupling between the neurons (g Gap ) was varied. The time course of the neuron membrane potentials for various baseline operating regimes and the strength of coupling is illustrated in Figure 1. The frequencies of bursting, as well as neuron resting potentials at various values of g Gap are shown in Table 1.
In our simulations shown in Figures 1A1-A3, both neurons (In1 and In2) incorporate I NaP . In all three panels, In1, if uncoupled, operates in a bursting mode, whereas In2 is silent (at low level of excitation, panel A1), exhibits bursting (at intermediate level of excitation, panel A2), or is tonic spiking (at a highest level of excitation, panel A3). The effect of gap junctional coupling depends on the difference between the membrane potentials of coupled neurons and the strength of coupling, and hence can be depolarizing for one neuron and hyperpolarizing for the other one. In all three cases, gap junctional coupling leads to synchronization of neuronal bursting from the burst ratio M:1, M>1 at a lower g Gap (0.375 < g Gap < 0.66 nS) to the ratio of 1:1 at higher g Gap (g Gap ≥ 0.66 nS) values. Note that coupling may reduce bursting frequency relative to the uncoupled case (see panel A1) or increase it (panels A2 and A3, see also Table 1) depending on the relative values of the resting membrane potentials of coupled neurons. This is seen in Figure 1A1 when In1 depolarizes In2 and induces bursting in this neuron with common bursting frequency at higher g Gap . At the same time, In2 hyperpolarizes In1 resulting in a general decrease of common bursting frequency relative to the uncoupled case with increased g Gap (see Table 1). In Figures 1A2,A3, we see the opposite situation. In these cases, both neurons operate in bursting regimes. However, the level of excitation of In1 is lower than that of In2, and prior to coupling, In2 generates bursts with a higher frequency than In1 in Figure 1A2 and is tonically active in Figure 1A3. Electrical coupling of these neurons leads to depolarization of In1 and hyperpolarization of In2. As a result, burst frequency of In1 increases while burst frequency of In2 decreases. In Figure 1A3, at lower g Gap In2 switches from tonic to bursting mode, when coupled, with the burst ratio 1:2. Similar to the simulation shown in Figure 1A1, at higher g Gap the two-cell model shows a synchronized bursting activity; however, the resulting frequency is higher than the initial In1 frequency (see Table 1). Hence gap junctional coupling between neurons with I NaP leads to synchronized rhythmic bursting with bursting frequency dependent on neuronal excitability and the strength of coupling.
In simulations shown in Figures 1B1-B3, only one of the two electrically coupled neurons (In1) has I NaP . Prior the coupling, this neuron can be silent (panel B1) or operate in a bursting regime (panels B2 and B3), whereas the other, "simple" neuron (In2), depending on the level of its excitation, can be silent (panels B1 and B2) or exhibit tonic spiking (panel B3). In all these simulations, the In2 neuron has higher resting potential if uncoupled and depolarizes In1 if the two neurons are connected by gap junctions (see Table 1). This results in the emergence of bursting activity (as in panel B1) or in an increase of the bursting frequency with increasing g Gap (panels B2 and B3). In all three cases, In1 hyperpolarizes In2 during interburst intervals but induces high-frequency spiking in In2 during bursts, resulting in synchronized activity of the two-cell model replicating bursting activity of In1.
The two neurons in Figure 2 are similar to those shown in Figure 1B1 and only one of the two electrically coupled neurons (In1) has I NaP . Here, instead of varying the strength of gap junctional coupling, it was set at 0.1 nS and the effects of varying the strength of I NaP are illustrated. Atḡ NaP below the level shown in Figure 2, both neurons remain silent. However, with increasingḡ NaP , bursting emerges and increases in frequency. The results demonstrate that the presence of I NaP in one neuron can initiate synchronized bursting in electrically coupled neurons even if both neurons when uncoupled are silent.
Overall, gap junctional coupling leads to equilibration of membrane potentials between the coupled neurons, regardless of the presence/absence of I NaP . Depending on the direction of change (depolarization or hyperpolarization), this coupling will either increase or decrease bursting frequency of coupled neurons with I NaP and bring neurons without I NaP from silent or tonic to bursting mode in synchrony with neurons with I NaP . In turn, I NaP in some neurons can amplify the effect of electrical coupling of these neurons with other neurons. Finally, both the presence of I NaP in some neurons and electrical coupling between neurons promotes neuronal synchronization and both are critical for initiation and support of populational bursting.

Excitatory Chemical Synaptic Connections in Two-Cell Model
In the case of excitatory connections by chemical synapses, we considered only unidirectional connections because reported chemical synapses between spinal interneurons have been unidirectional (Dougherty et al., 2013;Chopek et al., 2018;Ha and Dougherty, 2018). Figure 3 and Table 2 illustrate effects of excitatory chemical synaptic connections in the two-cell model when both neurons (source neuron In1 and target neuron In2) are bursting neurons incorporating I NaP (Figures 3A1-A3), or when I NaP is incorporated in only In1 (Figures 3B1,B2) or only In2 neuron (Figures 3C1,C2).
In Figures 3A1-A3, the target neuron (In2) that incorporates I NaP operates in various regimes without connection: silent ( Figure 3A1), bursting (Figure 3A2), and tonic spiking ( Figure 3A3). The synaptic input from In1 to In2 leads to bursting activity in In2 with a burst ratio of 2:1 at a low w Syn (0.5) or synchronized bursting of both neurons (at higher w Syn values in Figures 3A1-A3, see also Table 2).
In Figures 3B1,B2, the target neuron In2 does not have I NaP and when uncoupled is either silent (Figure 3B1) or generates tonic spiking (Figure 3B2). In these cases, when In2 receives rhythmic excitatory synaptic input from the source neuron In1, it follows In1 activity and generates bursts synchronously with the source. In Figure 2B2, the bursting activity in In2 occurs on the background of baseline spiking activity.
In Figures 3C1,C2, the source neuron In1 is a "simple" neuron without I NaP that generates a sustained spiking activity providing constant excitatory synaptic drive to the target In2 neuron with I NaP . This drive depolarizes In2 and this depolarization increases with the strength of synaptic input (w Syn) . Therefore, depending on the baseline regime of the target In2 neuron, this input results either in the emergence of bursting activity (panel C1) or the increasing of bursting frequency (panel C2, see also Table 2) in the target In2 neuron.
The above simulations have shown that chemical connections in the two-cell model promote neuronal synchronization and have variable effects on the resultant common bursting frequency.

Populational Models
To study the synchronization and frequency control in a heterogeneous population of neurons connected by electrical and/or excitatory chemical synapses, we developed a model of a population of 100 neurons consisting of neurons with the intrinsic bursting properties (based on the incorporated I NaP ) and "simple" neurons without such properties (with base ratio of 40/60%). All neurons were described in the Hodgkin-Huxley style as in the two-cell models described above. The neurons in the population were sparsely connected by electrical and/or excitatory chemical synapses with probabilities of p Gap = 0.3 (Ha and Dougherty, 2018) and p Syn = 0.1 (Dougherty et al., 2013; see also Carroll and Ramirez, 2013;Radosevic et al., 2019). Following Butera et al. (1999b), we chose the leakage reversal potential, E L , the leakage conductance, g Leak , and the maximum conductance of the persistent sodium current,ḡ NaP , to be randomly distributed over population to provide natural heterogeneity between single neuron states. These parameters are the ones that define the basic level of neuronal excitation and their ability to intrinsically generate bursting activity (Butera et al., 1999a,b;Del Negro et al., 2002a;Koizumi and Smith, 2008). E L ,ḡ NaP , and g Leak for neurons were distributed normally with a priory assigned means and standard deviations (see section Materials and Methods). The major focus of our studies was on the generation of populational activity and frequency control of populational bursting for different types of connections and depending on the relative expression of different types of neuron connectivity and connection strength.

Performance of Populational Model With Electrical Coupling
To investigate the role of gap junction in populational bursting, we simulated a population of 100 neurons with bidirectional electrical coupling with coupling probability of 0.3 (Ha and Dougherty, 2018). Figure 4 shows the model performance for different values of coupling (g Gap ) in the population. In each panel, the upper row shows the averaged populational activity, the middle row shows raster plots of spikes elicited by all neurons in the population, and the bottom row shows traces of membrane potentials of several sample neurons with different excitabilities. Neurons with I NaP are shown by red and neurons without I NaP are shown by black. In Figure 4A, neurons in the population are uncoupled (g Gap = 0). Because of the random distribution of neuronal parameters, neurons operate in different basic regimes: some neurons with I NaP are silent, some generate bursting activity with different burst frequencies, and some are tonically active. Similarly, some neurons without I NaP are silent and some demonstrate tonic spiking activity with different frequencies. The integrated neuron activity is sustained (upper row in Figure 4A). With increased electrical coupling (g Gap = 0.033 nS, Figure 4B), most of neurons incorporating I NaP synchronize their bursting activity and some neurons without I NaP become involved in populational bursting. Further increase of g Gap (g Gap = 0.066 nS, Figure 4C) results in synchronization of bursting activity of all neurons incorporating I NaP and involvement of the majority of neurons without I NaP in populational bursting, which increases frequency and amplitude of populational bursting activity (see upper row in Figure 4C). However, after all neurons are involved in populational bursting, further increase of g GaP does not lead to increasing frequency of oscillations. Instead, burst duration starts to increase and interburst intervals to decrease and eventually all neurons switch to sustained tonic activity (not shown).
Interestingly, the frequency and amplitude of populational bursting depends on expression of I NaP in neurons of the population. Figure 5 shows the dependence of the frequency and amplitude of populational bursting on the average value of the I NaP conductance (ḡ NaP ) for three different values of N NaP, defining the number of neurons with I NaP in the population. As N NaP increases, the populational bursting starts at lower values ofḡ NaP . The burst frequency increases with increasing values of bothḡ NaP and N NaP (Figure 5A), and amplitude of populational bursts is higher at greater N NaP values but depends onḡ NaP only at lowerḡ NaP values ( Figure 5B).

Gap Junctional Coupling and Locomotor-Like Activity in the Spinal Cord
Previous experimental studies have shown that blocking gap junctions with carbenoxolone reduces the frequency of ventral root bursts of locomotor-like activity in the isolated rodent cord evoked by the mixture of NMDA and 5-HT (Ha and Dougherty, 2018). Figures 6A1,A2 show an example of ventral root recordings during drug-evoked locomotion prior to and 30 min after carbenoxolone application begins. Burst frequency and amplitude are quantified in Figures 6B1,B2 demonstrating that blocking gap junctions in these preparations reduces the frequency of locomotor oscillations and, in most cases, the amplitude of integrated locomotor activity. To simulate the effect of blocking electrical synapses, we reduced the strength of gap junctional coupling between neurons by 50% in our model. Figures 6C1,C2 demonstrate the average populational activity for a sample simulation when g Gap was reduced by 50%. To generally estimate the effect of this reduction on frequency and amplitude of synchronized bursting, we ran several simulations for populations with different randomization of neuron parameters. As can be seen in Figures 6D1,D2, regardless of parameter distribution, our model exhibits changes in the frequency and amplitude of oscillations with a decrease of g Gap similar to experimental results shown in Figures 6A1,A2,B1,B2.

Performance of Populational Model With Chemical Coupling
In the next stage of this study, we considered the same population of 100 neurons, but replaced the bidirectional gap junction FIGURE 3 | Effects of the excitatory chemical coupling on activity of two neurons (In1 and In2). (A1-A3) Both neurons incorporate the persistent sodium current (I NaP ). In1 is in a bursting mode and In2, if uncoupled, is silent (A1), or operates in a bursing regime (A2), or is tonically active (A3). (B1,B2) Only In1 incorporates I NaP and operates in a busting regime and In2, if uncoupled, is either silent (B1) or tonically active (B2). (C1,C2) In1 does not incorporate I NaP and and is tonically active. In2 incorporates I NaP and if uncoupled, is either silent (C1) or operates in a bursting regime (C2). The weight of synaptic connection increases from left to right.
connections within the population with unilateral excitatory chemical connections with probability of p Syn = 0.1 (Carroll and Ramirez, 2013;Dougherty et al., 2013;Radosevic et al., 2019).  number of tonically active neurons both among neurons with and without I NaP (Figures 7A,B). Simultaneously, the bursting neurons start to consolidate in clusters including synchronously active neurons with I NaP and partially involving neurons without I NaP in bursting activity. Further increase of w Syn increases frequency of bursting in neurons with I NaP and synchronizes all neurons, resulting in populational synchronous bursting of high frequency and amplitude. In contrast to the population of electrically coupled neurons, in the population with chemical coupling, an increase in w Syn is accompanied by a simultaneous increase of the background tonic activity.

Contribution of Neuronal Coupling by Electrical and Chemical Synapses to Populational Activity
Experimental studies have shown that electrical coupling between Shox2 neurons is age-dependent. It is prevalent between functionally related Shox2 neurons in neonatal spinal cords and decreases in incidence and strength with age (Ha and Dougherty, 2018). To investigate how the activity of the population changes when both electrical and excitatory chemical connections are varied, we simulated a population in which both electrical and excitatory chemical synapses are present and sparsely distributed among neurons with probabilities of p Gap = 0.3 and p Syn = 0.1, and introduced two variables G and W, characterizing in average the total network electrical or chemical synaptic input to each neuron in the population, respectively (see Materials and Methods).
For each pair (G,W), we estimated stability of the populational bursting activity by calculating values of the period and amplitude of populational bursts and their standard deviations. All simulation results were divided in four groups: no populational activity, low amplitude or highly unstable populational activity, stable populational bursting, and tonic populational activity (see Materials and Methods). Figures 8A,B show regimes of populational activity and 2D diagrams of colorcoded frequency and amplitude of stable populational bursting in [G, W] plane. As seen in Figure 8A, if W is fixed and G is increasing, for most values of W, frequency of populational bursting does not increase smoothly but there are distinct areas in which the frequency does not change much and transition points when frequency increases abruptly. For example, at lower values of W (W < 7) there are one or two distinct transitions from the dark blue to light blue to cyan or the dark blue to cyan color in the diagram. This frequency transition is even more sharp at intermediate values of W (i.e., W = ∼10): dark blue abruptly changes to yellow. Our simulations have shown that in the areas in which frequency does not change much, the increasing G leads predominantly to synchronization of bursting neurons with I NaP and involving in the populational bursting some neurons without I NaP . However, at some higher values of G, the remaining tonically active neurons depolarize neurons with I NaP that have not been initially involved in bursting, leading to abrupt transitions to a higher frequency. Then again increasing G results predominantly in synchronization of bursting neurons, this time at a higher frequency of populational bursting. Depending on particular distribution of neuronal parameters these transitions might happen up to three times when W is fixed and G is increasing. After all neurons are involved in populational bursting and synchronized, further increase of G does not produce increases in frequency of populational bursting. On the contrary, at higher G, the length of the populational bursts is increasing and the interburst interval is decreasing, leading to smooth transition to the slower populational bursting (shown by smooth transition from the cyan/yellow to blue in Figure 8A) and eventually the population switches to sustained tonic activity (shown for higher W values by black in Figures 8A,B).
With the increase of W, the unilateral excitatory chemical connections allow for a higher bursting frequency. If W is high enough, (if W > 12 in our simulations) the model is capable of generating synchronous populational bursting at low values of G or even if electrical coupling in the population is absent. With a fixed G, the frequency of populational bursting is increasing with an increase of the total strength of the excitatory chemical connections (W). However, the range of G in which the model generates populational bursting activity decreases with increasing W due to the earlier transition to sustained tonic activity of the majority neurons (black area in Figures 8A,B). Figure 8B shows that the average changes in burst amplitude depending on the values G and W are quite smooth because the amplitude depends in a greater extent on how wellneuron activity in the population is synchronized and neuronal synchronization directly depends on increasing total strengths of the electrical and/or excitatory chemical connections (see also Figures 4, 7). Figures 8C1-C3 demonstrate the results of simulations for three pairs of G and W values shown by small red circles in Figures 8A,B. In Figure 8C1, W = 0 and G = 6 nS. In this case, G is high enough to generate synchronized populational bursting activity in the absence of excitatory chemical connections. In Figure 8C2, G = 3 nS and at that strength of electrical coupling the model cannot generate synchronized populational bursting in the absence of excitatory chemical connections. However, populational bursting can be reestablished with increasing W, which is illustrated in Figure 8C2 for W = 7.5. On the other hand, if W is high enough, populational bursting activity emerges even without electrical coupling of neurons in the population as can be seen in Figure 8C3 at W = 15 and G = 0 nS. It should be noted that in the latter case, a higher level of background activity is present.
Altogether, our simulations indicate that both electrical and excitatory chemical connections between neurons in the heterogeneous population incorporating a subset of intrinsically bursting neurons support synchronization of activities of individual neurons resulting in the rhythmic populational bursting activity. However, the effects of increasing strength of electrical and excitatory chemical coupling of frequency of populational bursting are different.  1 in (B), and 1.6 in (C). The probability of connection is 0.1. Arrangement and notation are the same as in Figure 4.

DISCUSSION
This study is based on a previous, experimental investigation by Ha and Dougherty (2018) that revealed the presence of electrical coupling between Shox2 non-V2a interneurons, which have been suggested to be critically involved in neuronal synchronization and generation of populational bursting in the neonatal rodent spinal cords. Here, we used computational modeling to investigate the potential roles of different neural interactions (gap junctional coupling and/or excitatory chemical synapses) and their interactions with the cellular I NaP -dependent bursting properties in the neuronal synchronization and generation of locomotor-like oscillations. We simulated and analyzed the behavior of pair of neurons operating in different regimes as well as heterogeneous populations of 100 neurons in different states, sparsely connected by either gap junctions or excitatory chemical synapses, and studied the possible role of electrical and chemical synaptic connections in populational bursting and control of frequency and amplitude of populational oscillations. The different aspects of our simulation results in connection with previous experimental findings as well as limitations of current studies are discussed below. Simulation results are divided in four groups: no populational activity (white areas), low amplitude or highly unstable populational activity (gray), stable populational bursting (colors), and tonic populational activity (black). For the stable populational bursting, for each (G, W), the amplitude of populational bursting was determined as the mean of maximal values of the averaged integrated activity of neurons (average number of spikes per neuron per second) during the bursts. The frequency of stable populational oscillations was determined as reciprocal to the mean bursting period. (C1-C3) Examples of simulations for three pairs of G and W values indicated in (A,B) by small red circles. Arrangement and notation in (C1-C3) are the same as in Figure 4.

The Role of Electrical Coupling Between Spinal Interneurons in Rhythmic Activity in the Isolated Spinal Cord of Neonatal Rodents
Electrical coupling between neurons is known to be widely present in different regions of the central nervous system of mammals, including rodent spinal cord (Tresch and Kiehn, 2000), especially in the early stages of the development (Walton and Navarrete, 1991;Chang et al., 1999). The role of gap junctions in neuronal synchronizations and generation of populational oscillations has been previously studied using computational neuron models of different complexities (Kepler et al., 1990;Sherman and Rinzel, 1992;Pfeuty et al., 2003;Kopell and Ermentrout, 2004). Gap junctional coupling has been demonstrated to interconnect spinal Shox2 non-V2a neurons (Ha and Dougherty, 2018) as well as Hb9 neurons (via other unidentified neurons, Wilson et al., 2007). It has been suggested that electrical coupling within each of these genetically identified populations significantly contributes to neuronal synchronization and populational bursting during spinal circuit development. Although the role of this coupling in the generation of locomotor-like oscillations in neonatal mouse cords has not been directly demonstrated for the specific spinal neuronal populations, pharmacological suppression of gap junctions in the cord in general resulted in a decrease of frequency and amplitude of locomotor-like oscillations recorded from ventral roots of the neonatal spinal cord during drug-evoked fictive locomotion (Ha and Dougherty, 2018). properties in the generation of these oscillations (Butera et al., 1999a,b;Del Negro et al., 2001, 2002aSmith et al., 2007Smith et al., , 2013Koizumi and Smith, 2008). The critical role of I NaP -dependent pacemaker oscillations in the pre-BötC has been challenged by some studies (Del Negro et al., 2002b;Feldman et al., 2013) and currently remains debated (Del Negro et al., 2018, but Jasinski et al., 2013Bacak et al., 2016;Molkov et al., 2016;Ausborn et al., 2018;Phillips and Rubin, 2019). Based on the previous computational models on the generation of respiratory oscillations in the pre-BötC in vitro (Butera et al., 1999a,b;Rybak et al., 2004), a series of computational models of spinal locomotor circuits has been developed by suggesting the critical role of I NaP in the generation of locomotor rhythmic activity in the spinal cord (Rybak et al., 2006a,b;McCrea and Rybak, 2007). This suggestion has been partly supported by several independent studies confirming the presence of I NaP in the neonatal rodent spinal cords (Tazerart et al., 2007(Tazerart et al., , 2008Zhong et al., 2007;Ziskind-Conhaim et al., 2008;Brocard et al., 2010Brocard et al., , 2013Harris-Warrick, 2010;Tong and McDearmid, 2012). Consequently, our recent computational models of spinal locomotor circuits which incorporated I NaP considered this current to be critically involved in the spinal cord rhythmogenesis (Zhong et al., 2012;Brocard et al., 2013;Rybak et al., 2013Rybak et al., , 2015Shevtsova et al., 2015;Danner et al., 2016Danner et al., , 2017Danner et al., , 2019Shevtsova and Rybak, 2016;Ausborn et al., 2019). The persistent sodium channels, including their voltage-dependent activation, have been characterized in spinal Hb9 interneurons present in the spinal cord (Brocard et al., 2013) and presumably contribute to the generation and/or propagation of locomotor-like oscillations (Wilson et al., 2007;Ziskind-Conhaim et al., 2008;Brocard et al., 2013). Although I NaP has not yet been characterized in Shox2 neurons, there is evidence for the expression of persistent inward currents in these neurons which likely includes I NaP (Ha et al., 2019). In this study, we followed the previous computational models of the spinal locomotor circuits developed by Rybak et al. by suggesting that I NaP plays critical role in the generation of locomotor-like oscillations. We, however, suggested that this current is present in, and randomly distributed over, only a subset of all neurons involved in generation of populational oscillations (with basic ratio of 40/60%). Our particular interest was in considering the behavior of neurons with and without I NaP connected by gap junctions and the possible roles of this current and electrical coupling in initiation of synchronized bursting activity. Our modeling results show that electrical coupling can lead to equilibration of membrane potentials between the coupled neurons. Depending on the membrane potentials of neurons, the electrical coupling can increase and decrease bursting frequency and bring neurons without I NaP from silent or tonic to bursting mode, hence amplifying neuronal synchronization. We have also shown that the presence of I NaP in coupled neurons can amplify the effect of electrical coupling, leading to emergence of bursting activity in initially silent neurons (see Figure 2). This corresponds to conclusion made from other studies that I NaP is a strong modulator of electrical synapses (Haas and Landisman, 2012). Finally, we have shown that the presence of I NaP in some neurons in cooperation with electrical and/or chemical coupling between neurons promotes neuronal synchronization in heterogenous neuronal population and that this current can be critical for initiation and support of populational bursting. Since incorporation of I NaP into only a subset of neurons in the modeled population was enough to generate a synchronized populational bursting when bidirectional connections at incidences seen electrophysiologically were incorporated, we would predict that to generate rhythmic locomotor-like oscillations in the neonatal spinal cord, it is sufficient for I NaP to be present in a subset of rhythm generating neurons but not necessarily in all of them.

Gap Junctions vs. Chemical Synaptic Connections
In isolated cords of newborn mice, suppressing chemical synaptic interactions between neurons with lower extracellular Ca 2+ concentrations or blocking action potentials using TTX did not eliminate rhythmic motor output (Tresch and Kiehn, 2000), hence confirming the important role of electrical coupling in neonatal cords. Similarly, in our model, oscillations are generated in the population without chemical synaptic interactions (at W = 0). Lowering extracellular Ca 2+ concentrations not only suppresses chemical synaptic interactions but has also been shown to amplify I NaP -dependent oscillations in spinal neurons (Brocard et al., 2013). Although not simulated by the present simplified model, this may serve as additional indirect evidence for the role of I NaP in oscillations generated by networks of electrically coupled neurons in the neonatal spinal cord. This, of course, can only be the case under conditions that I NaP in the spinal cord is TTX-resistant (see Cummins et al., 1999;Herzog et al., 2001;Dai and Jordan, 2011). These issues remain for future study. Electrical transmission has been shown to be prevalent during early postnatal period but decreases as the animal matures (Walton and Navarrete, 1991;Chang et al., 1999). Although the decline is seen within the first postnatal week in motor neurons (Walton and Navarrete, 1991;Chang et al., 1999), gap junctional coupling within both Shox2 non-V2a neuron and Hb9 neuron populations remains relatively high throughout the first two postnatal weeks (Hinckley and Ziskind-Conhaim, 2006;Ha and Dougherty, 2018). The probability of coupling among neighboring neurons was found to be ∼0.3 (Ha and Dougherty, 2018) and this is the value that was chosen for the modeling study. The probability value of 0.1 was set for chemical synapses and falls more in line with a prior Shox2 connectivity study (Dougherty et al., 2013). There were a few key differences between the reports of Shox2 connectivity. First, the experiments showing connectivity by chemical synapses at a rate of ∼10% (Dougherty et al., 2013) was performed in mice where the Chx10 absence/presence was not distinguished, so it included the entire Shox2 population as a single sampling group. Connections were tested in a dorsal horn removed preparation and were typically within the field of view of a 60x objective but not necessarily neighboring neurons and there were no instances of electrical connections reported. In Ha and Dougherty (2018), pairs of neurons recorded were in close proximity (∼50 µm on average) and were selected based on the visualization of processes crossing and in close apposition. Further, these neurons were distinguished based on the presence/absence of Chx10 expression. Therefore, the 33% incidence of electrical synapses is only among very close neighboring neurons of the same group. In this study, chemical connections were only found at 2%, indicating that neighboring neurons are less likely to be connected by chemical synapses. Both values of chemical connections are likely to be underestimated, particularly the one in which only near neighboring neurons were recorded. Interestingly, the discrepancies in the results of these studies reinforces the idea that electrical connectivity is mechanism for local synchronization.
Gap junctions behave as low pass filters and efficiently convey membrane potential fluctuations. Additionally, action potentials in one Shox2 neuron reliably lead to spikelets in gap junctionally connected neurons. However, chemical synapses between neonatal Shox2 neurons had high failure rates (Ha and Dougherty, 2018). Therefore, it is possible that there is a higher incidence of chemical synaptic connections in older ages, not due to new synapses but increased efficacy at synapses established early in development. The experimental study considered single pairs using dual recordings and responses to current injection in the presynaptic neuron. It is unlikely that the postsynaptic cell receives input from a single presynaptic Shox2 neuron. Gap junctional coupling between multiple presynaptic neurons will synchronize their firing, as demonstrated in the current study, likely leading to a more reliable excitation of the postsynaptic neuron, even if excitatory chemical synapses are sparse and single synapses have high synaptic failure rates. This could provide a mechanism for increased efficacy of excitatory transmission not only within the population but in transmission to downstream targets.
Using our model, we investigated how the activity of the population might change with changing electrical and chemical synaptic connections. Our simulations demonstrated that both electrical and excitatory chemical connections between neurons in the heterogeneous population support synchronization of activities of individual neurons resulting in the rhythmic populational bursting activity. This implicitly confirms that the locomotor-like oscillations in spinal cord can be maintained over a continuum. Thus, the results suggest that gap junctional coupling is capable of strongly driving synchronization of the population. As bidirectional electrical coupling declines with age, other mechanisms can take over to maintain rhythmic bursting within the population. Our simulations demonstrate that the change to a dominance of unidirectional chemical synaptic interactions, particularly with I NaP present in a subset of these neurons, is sufficient to fulfill this role and to expand the range of attainable frequencies. Interestingly, our simulation results have demonstrated that a strong, synchronous population bursting can be produced in a population when as low as 10% of neurons are connected by excitatory chemical synapses.
Our simulations have also demonstrated that the effects of increasing the strength of electrical and excitatory chemical coupling on the frequency of populational bursting are different. In both cases, a subset of neurons in the population generates spiking activity and provides an excitatory tonic input to other neurons which is enhanced with increasing connection strength and results in an increased frequency of populational bursting. In addition, in the case of bidirectional electrical coupling, tonically active neurons directly depolarize other neurons, elevating neuronal resting potentials and leading to an increase in frequency of populational bursting. However, in the case of bidirectional electrical coupling, bursting neurons have a hyperpolarizing effect on tonically firing neurons and eventually an increase of g Gap turns spiking activity of majority of tonically firing neurons into bursting, as was shown in our simulations (see Figures 1A3,B3, 4) and in earlier studies (Sherman and Rinzel, 1992). This prevents further increase in bursting frequency. After all of the neurons become involved in populational bursting and their membrane potentials reach equilibrium, further increase of g Gap does not increase the frequency of populational oscillations (see also Kepler et al., 1990). Instead, the frequency starts to decrease. In contrast, in populations with chemical coupling, elevation of w Syn increases the activity of tonically active neurons which depolarize bursting neurons, hence leading to the increase the frequency of populational bursting.

Limitations of the Present Studies and Future Directions
In this modeling study we limited single neuron models by incorporating only minimal number of ionic channels necessary for generation of neuronal spiking activity, such as fast sodium (Na) and potassium rectifier (K), and persistent sodium (NaP) involved in bursting rhythmic activity. Therefore, these models remain largely generic. At the same time, we know that spinal neurons contain other ionic channels, including T-type calcium and h currents and different Ca 2+ and Ca 2+ -dependent potassium currents that can also support pacemaker properties and have been implicated in locomotor-like activity (Wilson et al., 2005;Anderson et al., 2012;Kiehn, 2016;Brocard, 2019). In addition, there is indirect evidence of an important role of Na + /K + pump in CPG operation (Kueh et al., 2016) and particularly in the neonatal spinal cords of mice (Picton et al., 2017). The Na + /K + pump current has been previously included in the computational models of Hb9 neurons (Brocard et al., 2013) and can be included in Shox2 neuron models in the future. This current is known to interact with h-current (Kueh et al., 2016) and can specifically affect neuronal coupling via electrical and chemical synapses. This awaits additional experimental and computational investigations. A more realistic model of different types of spinal interneurons that includes different ionic channels with experimentally measured characteristics should be developed and investigated (see for e.g., Sharples et al., 2020), which represents one of the directions of our future studies.
The other limitation of the present study is the consideration of only excitatory interactions within the spinal network. The potential role of electrical coupling between inhibitory interneurons and the role of inhibitory network interactions in cooperation with gap junction (Bennett, 1997;Bennett and Zukin, 2004;Kopell and Ermentrout, 2004;Connors, 2017) will also be a focus of our future investigations.
The contribution of gap junctional coupling within specific populations to locomotor rhythmicity remains to be tested experimentally. The pharmacological studies are limiting in that the gap junction blockers are acting at all electrical synapses, which have been shown between several classes of interneurons (Hinckley and Ziskind-Conhaim, 2006;Wilson et al., 2007;Zhong et al., 2010;Chopek et al., 2018) and motor neurons (i.e., Walton and Navarrete, 1991;Rash et al., 1996;Chang et al., 1999;Tresch and Kiehn, 2000). Additionally, carbenoxolone, the most common gap junction blocker, has been shown to have additional actions which may affect rhythmicity (Elsen et al., 2008;Tovar et al., 2009;Connors, 2012). Therefore, conditional knockout studies are required to determine implications more precisely.
Finally, the mutual interactions between different classes of genetically identified excitatory and inhibitory interneurons present in the mammalian spinal cord also awaits further experimental and computational investigations.

DATA AVAILABILITY STATEMENT
The modeling datasets for this study are publically available at https://github.com/RybakLab/Neural-interactionsin-spinal-networks.

ETHICS STATEMENT
This animal study was reviewed and approved by Drexel University Institutional Animal Care and Use Committee.

AUTHOR CONTRIBUTIONS
NS, NH, IR, and KD: conceptualization, investigation, and writing (review and editing). NS, IR, and KD: writing (original draft). NH and KD: experiments and data curation. NS: modeling, formal analysis, and software. NS and KD: visualization. IR and KD: supervision, project administration, and funding acquisition. All authors contributed to the article and approved the submitted version.

FUNDING
This study was supported by grants from the National Institutes of Health (R01 NS090919; R01 NS095366; R01 NS100928; and R01 NS110550).