Cross-Frequency Coupling in Descending Motor Pathways: Theory and Simulation

Coupling of neural oscillations is essential for the transmission of cortical motor commands to motoneuron pools through direct and indirect descending motor pathways. Most studies focus on iso-frequency coupling between brain and muscle activities, i.e., cortico-muscular coherence, which is thought to reflect motor command transmission in the mono-synaptic corticospinal pathway. Compared to this direct pathway, indirect corticobulbospinal motor pathways involve multiple intermediate synaptic connections via spinal interneurons. Neuronal processing of synaptic inputs can lead to modulation of inter-spike intervals which produces cross-frequency coupling. This theoretical study aims to evaluate the effect of the number of synaptic layers in descending pathways on the expression of cross-frequency coupling between supraspinal input and the cumulative output of the motoneuron pool using a computer simulation. We simulated descending pathways as various layers of interneurons with a terminal motoneuron pool using Hogdkin–Huxley styled neuron models. Both cross- and iso-frequency coupling between the supraspinal input and the motorneuron pool output were computed using a novel generalized coherence measure, i.e., n:m coherence. We found that the iso-frequency coupling is only dominant in the mono-synaptic corticospinal tract, while the cross-frequency coupling is dominant in multi-synaptic indirect motor pathways. Furthermore, simulations incorporating both mono-synaptic direct and multi-synaptic indirect descending pathways showed that increased reliance on a multi-synaptic indirect pathway over a mono-synaptic direct pathway enhances the dominance of cross-frequency coupling between the supraspinal input and the motorneuron pool output. These results provide the theoretical basis for future human subject study quantitatively assessing motor command transmission in indirect vs. direct pathways and its changes after neurological disorders such as unilateral brain injury.

Introduction

The human motor system is a highly cooperative network comprised of different groups of neurons. Neural coupling, i.e., the synchronization of neural activity across these groups, is key to signal transmission among functionally related, though anatomically distant, neuronal groups (e.g., the motor cortices and spinal motoneuron pool) through direct and indirect descending pathways (van Wijk et al., 2012). Over decades, most researchers investigating neural coupling in the motor descending pathways have focused on the synchronization between cortical oscillations and muscle activities at the same frequency (i.e., iso-frequency coupling), known as the cortico-muscular coherence (Mima and Hallett, 1999). It is thought to reflect motor command transmission in the mono-synaptic corticospinal tract (Schoffelen et al., 2005). Previous simulation and in vivo studies demonstrated that, in this direct descending pathway, despite the non-linearity of individual neurons, neural oscillation of the supraspinal input could be linearly transmitted to the cumulative output of the motoneuron pool at the same frequency (Negro and Farina, 2011a,b). These previous studies explained the origin of iso-frequency coupling between the supraspinal input and the motoneuron pool output with respect to the use of the monosynaptic corticospinal tract as the fastest, direct descending pathways in healthy individuals.

However, the corticospinal tract is not the only motor pathway in humans. There are other indirect pathways (e.g., cortico-reticulospinal tract, rubrospinal tract) in parallel with the direct corticospinal tract (Dum and Strick, 1991; Jang and Seo, 2014). Although contributions from these indirect motor pathways are relatively small compared to the corticospinal tract in healthy individuals, they do still play important roles in various motor control tasks such as postural control during movement (Drew et al., 2004). Furthermore, in some neurological disorders, such as unilateral brain injury, the reliance on these indirect motor pathways may increase due to losses of corticospinal projections (Fries et al., 1993; Jang et al., 2013; Owen et al., 2017). The injury-induced increased reliance on these indirect motor pathways is likely associated with motor impairments (e.g., abnormal limb synergies and spasticity) post unilateral brain injury (Ellis et al., 2012, 2017; McPherson et al., 2018a,c; Li et al., 2019). Thus, investigating the neural coupling in these indirect motor pathways will allow for a more complete understanding of the transmission of motor commands from the brain to muscles, and may pave the way for quantitative assessments of the usage of indirect motor pathways in both normal and pathological motor control.

Compared to the direct corticospinal tract, these indirect motor pathways involve multiple synaptic connections via interneurons. Neuronal processing of synaptic inputs can lead to the modulation of inter-spike intervals which produces cross-frequency coupling, i.e., synchronization across different frequencies between input and output (Koch and Segev, 2000; Markram, 2003; Yang et al., 2018). Our previous work on multi-synaptic ascending sensory pathways (Yang et al., 2016b; Tian et al., 2018), as well as a recent opinion article (Yang et al., 2018), argued that multi-synaptic interaction in a neural pathway can lead to a substantial expression of cross-frequency coupling. However, insights into possible mechanisms underlying neural coupling in the multi-synaptic descending motor pathways are currently lacking. Focusing on the iso-frequency coupling (e.g., cortico-muscular coherence) only one previous study indicated that the input from the indirect motor pathways can reduce the iso-frequency coupling between the cortical input and motoneuron pool output (Negro and Farina, 2011a) while no insight has been provided into the neural mechanisms of cross-frequency coupling. This study aims to systematically evaluate the effect of the number of synaptic connections or interneuron layers on the expression of cross-frequency coupling between supraspinal input and output of the motoneuron pool using computer simulations. We hypothesize that multi-synaptic interaction in an indirect descending motor pathway increases the non-linear distortion of efferent motor signal transmission, resulting in enhanced cross-frequency coupling over iso-frequency coupling.

To test our hypothesis, we simulated descending pathways as various layers of interneurons in cascade with a terminal motoneuron pool, using Hodgkin-Huxley styled neuron models (Booth et al., 1997; Rybak et al., 2006). Both cross- and iso-frequency coupling between the input (which comprised of a supraspinal drive with an independent membrane noise) and the output of the motoneuron pool were computed using a recently developed generalized coherence method (Yang et al., 2016b). The ratio of cross- to iso-frequency coupling was calculated to determine which type (cross- or iso-frequency) of neural coupling is dominant and how it changes with an increasing number of synaptic connections or interneuron layers.

Methods

Motoneuron and Interneuron Models

We simulated descending pathways as various layers of interneurons in a cascade with a terminal motoneuron pool. All neurons were modeled in the Hodgkin–Huxley style. A two-compartment model comprising of a soma and a dendrite was used for simulating motoneurons (Booth et al., 1997). Because of the lack of adequate experimental data, a single compartment simplification of this model was used to simulate the interneurons (Rybak et al., 2006).

The motoneuron model incorporated the following ionic currents (with the corresponding channel conductances): fast sodium (I_Na with maximal conductance g_Na), persistent (slowly inactivating) sodium (I_NaP with maximal conductance g_NaP), delayed-rectifier potassium (I_K with maximal conductance g_K), calcium-N (I_CaN with maximal conductance g_CaN), calcium-L (I_CaL with maximal conductance g_CaL), calcium-dependent potassium (I_K,Ca with maximal conductance g_K,Ca), and leakage (I_L with constant conductance g_L) currents (Lee and Heckman, 2001; Darbon et al., 2004; Rybak et al., 2006; Streit et al., 2006):

$\begin{array}{c}\hfill I_{\text{Na}}=g_{Na}\times m_{Na}^3\times h_{Na}\times (V-E_{Na});\\ \hfill I_{\text{NaP}}=g_{NaP}\times m_{NaP}\times h_{NaP}\times (V-E_{Na});\\ \hfill I_{\text{K}}=g_K\times m_K^4\times (V-E_K);\end{array}$

$\begin{array}{l}\hfill I_{\text{CaN}}=g_{CaN}\times m_{CaN}^2\times h_{CaN}\times (V-E_{Ca});\\ \hfill I_{\text{CaL}}=g_{CaL}\times m_{CaL}\times (V-E_{Ca});\\ \hfill I_{\text{K},\text{Ca}}=g_{K,Ca}\times m_{K,Ca}\times (V-E_K);\\ \hfill I_{\text{L}}=g_L\times (V-E_L);\end{array}$

where V is the membrane potential of the corresponding neuron compartment [i.e., soma (V_(S)) or dendrite (V_(D))] in two-compartment models, or the neuron membrane potential V in the one-compartment interneuron model which is explained later). E_Na, E_K, E_Ca, and E_L are the reversal potentials for sodium, potassium, calcium and leakage currents, respectively. The variables m and h (with subscripts indicating ionic channels) represent the activation and inactivation variables of the corresponding ionic channels, as described by the following differential equations:

$\begin{array}{l}{\tau }_{mi}(V)\frac{d}{dt}{m}_i=m_{\infty i}(V)-m_i\\ {\tau }_{hi}(V)\frac{d}{dt}{h}_i=h_{\infty i}(V)-h_i\end{array}$

where i indicates the name of the channel, m_∞i(V) and h_∞i(V) represent the voltage-dependent steady-state activation and inactivation, and τ_mi(V) and τ_hi(V) are the corresponding time constants (see Booth et al., 1997; Rybak et al., 2006 for details of these parameters). The instantaneous value of m_K,Ca was calculated from the intracellular Ca²⁺ concentration of the corresponding compartment as (Booth et al., 1997):

$\begin{array}{l}{m}_{K,Ca}=\frac{Ca}{Ca+K_d}\end{array}$

where Ca is the Ca²⁺ concentration of the corresponding compartment of the neuron and K_d is the half-saturation level of this conductance. The kinetics of intracellular Ca²⁺ concentration (|Ca|) were computed separately for each compartment according to the following equation:

$\begin{array}{l}\frac{d}{dt}|Ca|=-f\times (\alpha I_{Ca}+k_{Ca}|Ca|)\end{array}$

where f defines the percentage of free to total Ca²⁺, α converts the total Ca²⁺ current, I_Ca, to Ca²⁺ concentration and k_Ca represents the Ca²⁺ removal rate.

The maximal channel conductances, equilibrium potentials and membrane capacitance of the neuron models were set with the same values as in Rybak et al. (2006). The details are specified in the Appendix. The equilibrium leakage potentials of the motoneurons and interneurons were set as described in section Simulations.

The dendrite–soma coupling currents (with conductance g_C) for soma (I_C(S)) and dendrite (I_C(D)) were calculated as (Booth et al., 1997):

$\begin{array}{l}{I}_{\text{C}(\text{S})}=\frac{g_C}{p}(V_{(D)}-V_{(S)})\\ I_{\text{C}(\text{D})}=\frac{g_C}{1-p}(V_{(S)}-V_{(D)})\end{array}$

where p is the parameter defining the ratio of somatic surface area to the total neuronal surface area.

We used conductance-based excitatory post-synaptic potentials (EPSPs) for simulating the synaptic inputs to each motoneuron. The synapses were modeled as exponentially decaying injected currents (I_SynE with peak conductance g_SynE and reversal potential E_SynE): I_SynE = g_SynE × (V – E_SynE) into the soma compartment (Negro and Farina, 2011a). The time constant τ_synE for the decay was 5 ms (Rybak et al., 2006). The peak conductance value for synapses on motoneurons was adjusted to produce an EPSP peak of 100 μV (Finkel and Redman, 1983).

With the inclusion of I_NaP to the motoneuron dendrite (Rybak et al., 2006), the membrane potentials of the motoneuron soma (V_(S)) and dendrite (V_(D)) were computed from the following equations:

$\begin{array}{l}C\frac{d{V}_{(S)}}{dt}=-I_{Na(S)}-I_{K(S)}-I_{CaN(S)}-I_{K,Ca(S)}-I_{L(S)}\\ -I_{C(S)}-I_{SynE}\\ C\frac{d{V}_{(D)}}{dt}=-I_{NaP(D)}-I_{CaN(D)}-I_{CaL(D)}-I_{K,Ca(D)}-I_{L(D)}\\ -I_{C(D)}\end{array}$

where C is the membrane capacitance and t is time.

The interneurons (single-compartment models) contain only a minimal set of ionic currents (Rybak et al., 2006):

$\begin{array}{l}C\frac{dV}{dt}=-I_{Na}-I_K-I_L-I_{SynE}\end{array}$

There is no existing literature reporting experimentally observed values of interneuron EPSPs of the descending pathways. Since interneurons are usually much smaller than motoneurons, they have higher input resistances and smaller somatic surface areas (Bui et al., 2003). Thus, we adjusted peak conductance of synapses on interneurons to produce an EPSP peak of 500 μV.

Input Signal and Connection Configuration

We used a probabilistic connection model (Ferrario et al., 2018) to simulate descending pathways with various layers and 100 neurons per layer (Lüscher et al., 1983) (see Figure 1). Each neuron in the first layer was fed by a time-varying injected current into the somatic compartment. The supraspinal input was designed as a Gaussian signal in beta band (15–35 Hz) [mimicking the cortical oscillations observed experimentally during motor tasks (Pfurtscheller and Da Silva, 1999)] with an added membrane noise (see Figure 2). The membrane noise was modeled as a bandlimited (1–100 Hz) Gaussian noise, which was independent for each neuron (Maltenfort et al., 1998). The total variance of this stochastic input was a percentage of the constant current injection to produce a mean ISI-CoV of (i) 0.55 for the 1st interneuron layer in case of multi-synaptic pathways (Prut and Perlmutter, 2003), and (ii) 0.2–0.3 for the motoneuron pool in case of the mono-synaptic pathway (Tanji and Kato, 1972; Sturm et al., 1997; Mattei and Schmied, 2002).

FIGURE 1

Figure 1. Simulation of descending pathways with various layers (N = 0, 1, 2, …) of interneurons in cascade with a terminal motoneuron pool. In, interneuron; Mn, motoneuron.

FIGURE 2

Figure 2. Simulated supraspinal input comprising of a Gaussian signal (15–35 Hz) with added Gaussian noise (1–100 Hz band-limited) with signal-to-noise ratio of ~−7.5 dB.

For the successive layers, the input to each neuron was the sum of output spike trains (convolved with the EPSP) of neurons randomly sampled from the previous layer. The number of neurons which contributed to the input of each interneuron was set to obtain a mean firing rate in the range of 19–24 spikes/s for the whole interneuron layer. This is in line with previous experimental observations in primate models during flexion/extension tasks (Prut and Perlmutter, 2003). The number of inputs to each motoneuron was set to 100 i.e., the sum of inputs from all interneurons of the terminal interneuron layer. The range of the firing rates was adjusted as explained in section Simulations. The mean firing rate of the active motoneurons (>8 spikes/s; Negro and Farina, 2011a) was thus obtained to be in the range of 16–19 spikes/s. Such a connection model resembles the anatomical course of various descending motor pathways which, via a varying number of interneuron layers, terminate on spinal motoneuron pools (Matsuyama et al., 2004).

Simulations

Simulations were run at a sampling rate of 1 kHz using 200 epochs with a 1-s duration per epoch. The resulting data were sufficient for a robust neural coupling analysis (Hagihira et al., 2001). Our simulated multi-synaptic pathways represented the part of descending pathway involving only spinal interneurons and motoneurons, since no reticular neurons were simulated due to the complete lack necessary parameters in the existing literature (McDougal et al., 2017). A previous study reported heterogeneity of excitatory spinal interneuron populations based on their firing rates. It found a non-monotonous decline in the mean firing rate histogram with a local peak at ~50 spikes/s (Prut and Perlmutter, 2003). In the simulations, we mimicked this histogram (see Figure 3A) by combining two random exponential distributions of leakage potential (E_L) values for each interneuron layer. The range of E_L was adjusted so that its mean was around −64 mV (Rybak et al., 2006). Motoneuron firing rates have been experimentally reported to be predominantly in the range of 5–30 spikes/s during isometric contractions of limb muscles (for contraction levels ≤60% of maximum voluntary torque) (De Luca and Hostage, 2010). We adjusted the leakage potentials of the motoneurons in the pool to generate a distribution (Rybak et al., 2006) of firing rates in a similar range (see Figure 3B).

FIGURE 3

Figure 3. Firing rate of simulated neurons. (A) Histogram of mean firing rates of the simulated interneuron layer with a secondary peak at ~ 50 spikes/s. (B) Scatter plot showing the distribution of mean firing rates of the simulated pool of motoneuron.

Neural Coupling Analysis

We used our recently developed generalized coherence measure, i.e., n:m coherence (NMC) (Yang et al., 2016b), to assess cross- and iso-frequency coupling between the simulated input and output signals. The n:m coherence is a straightforward extension of the linear coherence used in corticomuscular coherence (Mima and Hallett, 1999) based on high-order statistics (Nikias and Mendel, 1993) for distinguishably determining cross- and iso-frequency coupling between signals. Thus, the iso-frequency coupling of our results obtained by this method would be comparable to previous corticomuscular coherence studies (Mima and Hallett, 1999; Mima et al., 2001; Yang et al., 2016a, 2018).

Let X(f), Y(f) be the Fourier Transform of two time series (e.g., the input and output signals). The NMC between them is defined as:

$\begin{array}{l}NMC(f_X,f_Y)=\frac{\left|S_{XY}(f_X,f_Y)\right|}{\sqrt{S_X^n\left(f_X\right)S_Y^m\left(f_Y\right)}}\end{array}$

for assessing cross-frequency (f_X ≠ f_Y) and iso-frequency (f_X = f_Y) coupling between signals, where m/n is the simple whole number ratio of f_X/f_Y (e.g., if f_X = 8, f_Y = 16 then m = 1, n = 2) and

$\begin{array}{l}{S}_{XY}(f_X,f_Y)=<X^n(f_X){\left(Y^m(f_Y)\right)}^{*}>,\hfill \\ S_X^n(f_X)=<X^n(f_X){\left(X^n(f_X)\right)}^{*}>\hfill \end{array}$

where < · > represents the averaging over epochs and ${\text{X}}^{\text{n}}=\underset{n}{\underbrace{X\left(f_x\right)\cdot X\left(f_x\right)\cdot \dots \cdot X\left(f_x\right)}}.$

The NMC reflects the strength of iso- or cross-frequency coupling between signals. When f_X = f_Y, we have m = n = 1, then the NMC is equivalent to the classical (linear) coherence for iso-frequency coupling (Yang et al., 2016a). When f_X ≠ f_Y, then the NMC indicates the non-linear coupling between signals across different frequency components (i.e., cross-frequency coupling) (Yang et al., 2015). Thus, the n:m mapping can generate harmonic (m = 1) and subharmonic coupling (m > 1) between the input and the output in the frequency domain (Yang et al., 2016b). As a generalized coherence method, the NMC is a metric indicating cross-frequency coherence between signals, which is different from other cross-frequency coupling methods such as the phase-amplitude coupling (De Hemptinne et al., 2013) reflecting how a low-frequency phase modulates a high-frequency amplitude.

According to Cauchy-Schwarz-inequality, we have:

$\begin{array}{l}\left|\langle X^n(f_X){(Y^m(f_Y))}^{*}\rangle \right|\le {\left(\langle {\left|X^n(f_X)\right|}^2\rangle \right)}^{1/2}{\left(\langle {\left|Y^m(f_Y)\right|}^2\rangle \right)}^{1/2}\end{array}$

Thus, the NMC is bounded by 0 and 1, where 1 indicates that two signals are perfectly coupled at the tested frequency pair (f_X, f_Y). As the NMC values are computed by comparing different frequency pairs between signals, the significant threshold was adapted with a Bonferroni correction to control the type I error (family-wise error rate: 0.05) (Yang et al., 2016b). There are 2,100 frequency pairs that were included for Bonferroni corrections, i.e., 21 frequencies in the input (from 15 to 35 Hz at 1 Hz resolution) × 100 frequencies in the output (from 1 to 100 Hz at 1 Hz resolution). More details of the NMC method is available in Yang et al. (2016b).

Since the supraspinal input had added independent noise for each 1st layer neuron, each coupling analysis was repeated 100 times, each time with a different realization of the independent noise (as described in section Input Signal and Connection Configuration) added to the supraspinal input in the same signal-to-noise ratio as the original input (i.e., ~-7.5 dB). To compare the dominance of cross- vs. iso-frequency coupling, we defined the cross-frequency coupling over iso-frequency coupling index as COI = (CFC–IFC)/(CFC+IFC), where CFC is the sum of all significant cross-frequency coupling values and IFC the sum of all significant iso-frequency coupling. We included only “significant” CFC and IFC values to exclude false positives in the coherence analysis. The range of COI is [−1, 1], where a larger COI indicates a more dominant cross-frequency coupling.

To examine the effect of the number of synaptic/interneuron layers on neural coupling of descending pathways, we computed NMC, IFC, CFC, and COI between the given supraspinal input and the cumulative spike train (CST) output of the simulated motoneuron pool (derived as the sum of individual motoneuron spike trains following Negro and Farina, 2011a). In addition, the IFC, CFC, and COI between the supraspinal input and the CST of each successive interneuron layer was also computed to evaluate how they change across layers. Furthermore, we also examined the combined effect of both mono- and multi-synaptic pathways on neural coupling by varying the weight of the input from either pathway (with the same supraspinal input) to the terminal motoneuron pool.

Results

Neural Coupling Between the Supraspinal Input and the Cumulative Output From Motoneuron Pool

Both iso-frequency coupling and cross-frequency coupling were detected in the simulated motor pathways (see Figure 4). The detected cross-frequency coupling includes harmonic coupling (i.e., output frequency over input frequency ratio n/m is an integer) and non-integer n:m coupling. This result is in line with previous experimental studies reporting both harmonic and non-integer coupling in the human sensorimotor system (Daffertshofer et al., 2000; Yang et al., 2016a). A higher amount of cross-frequency coupling was observed in the multi-synaptic pathways where there are one or more interneuron layers.

FIGURE 4

Figure 4. Neural coupling between the supraspinal input and cumulative spike train output of the terminal motoneuron pool for the descending motor pathways with 0, 1, and 3 interneuron layers. The iso-frequency (1:1) coupling is indicated by the green dashed line. The detected cross-frequency coupling including harmonic coupling (i.e., output frequency over input frequency ratio n/m is an integer, e.g., 2:1, 3:1, 4:1, indicated by the red dashed lines) and non-integer n:m coupling (other points in the map).

To examine how iso-frequency coupling and cross-frequency coupling evolved with increasing interneuron layers, we computed the IFC, CFC, and COI between the given supraspinal input and the motoneuron pool output for simulated pathways with various layers (N = 0, 1, 2, 3, …) of interneurons in cascade with a terminal motoneuron pool (see Figure 5). Using one-way ANOVA we found that the number of interneuron layers had significant effect on IFC [F_{(10, 1089)} = 3613.36, p < 0.001], CFC [F_{(10, 1089)} = 2934.50, p < 0.001] and COI [F_{(10, 1089)} = 7108.90, p < 0.001]. We used Tukey's honest significant difference (HSD) criterion for post-hoc comparisons, with Bonferroni correction to control the type I error. Hence, we adjusted the threshold p-value as 0.05/k to control the family-wise error rate to be <0.05, where k is the number of post hoc comparisons (k = 10). We found that the IFC decreased with increasing number of interneuron layers in the pathways (p < 0.05/10), while the CFC increased in the pathway with up to three interneuron layers (p < 0.05/10). Their combined effect resulted in an initial increase in COI with interneuron layers (p < 0.05/10) and a saturation for the pathways with more than three interneuron layers. It was also observed from the COI values that the iso-frequency coupling is only dominant (COI < 0) in the mono-synaptic descending pathway where the supraspinal input directly drives the motoneuron pool without passing through any interneuron layer. Consequently, the cross-frequency coupling became dominant (COI > 0) when there were interneuron layers in the descending pathway.

FIGURE 5

Figure 5. IFC, CFC, and COI between the supraspinal input and the motoneuron pool output in simulated pathways with various layers of interneurons in cascade with a terminal motoneuron pool. The depicted values represent the mean (with ±2 standard deviation as indicated by error bars) calculated from 100 repetitions of the coupling analysis, each being run with a different realization of additive noise to the supraspinal input (as described in section Neural Coupling Analysis). Tukey's test was performed to test for significant decrease in IFC and increase in CFC and COI between the pathways containing n − 1 and n (n = 1, 2, 3…) interneuron layers. To control the family-wise error rate, we set the threshold p = 0.05/k (k = 10, i.e., number of comparisons). Asterisks in superscript of the n-th layer number indicate a significant change of the results in the n-interneuron layer motor pathway in comparison to (n − 1)-interneuron layer motor pathway (**p < 0.05/10).

Neural Coupling Between the Supraspinal Input and the Output From Successive Neuron Layers

Using one-way ANOVA, we also found that the number of interneuron layers had a significant effect at p < 0.05 level on the IFC [F_{(9, 990)} = 3235.46, p < 0.001], CFC [F_{(9, 990)} = 2113.54, p < 0.001] and COI [F_{(9, 990)} = 5597.45, p < 0.001] between the supraspinal input and the output from successive interneuron layers. Using Tukey's HSD criterion with Bonferroni correction (k = 9), the IFC was observed to decrease across successive interneuron layers (p < 0.05/9) while the CFC and COI increased (p < 0.05/9) up to the fifth and sixth layer, respectively (see Figure 6). Additionally, in multi-synaptic pathways, IFC dropped more at the terminal motoneuron layer of the n-layer pathway in comparison to that of the terminal interneuron layer of the n + 1 layer pathway (for n = 1–11, unpaired t-test, p < 0.001 in all cases).

FIGURE 6

Figure 6. IFC, CFC, and COI between the supraspinal input and cumulative spike train output of successive neuron layers in simulated pathways with various layers of interneurons in cascade with a terminal motoneuron pool. The 1-layer pathway represents the monosynaptic tract. The depicted values represent the mean (with ±2 standard deviation as indicated by error bars) calculated from 100 repetitions of the coupling analysis, each being run with a different realization of additive noise to the supraspinal input (as described in section Neural Coupling Analysis). Tukey's test was performed to test for significant decrease in IFC and increase in CFC and COI between the n − 1th and nth (n = 2, 3, 4…) interneuron layer. To control the family-wise error rate, we set the threshold p = 0.05/k (k = 9, i.e., number of comparisons) Asterisks in superscript of the n-th layer number indicate a significant change of the result in that layer in comparison to that of the previous layer (**p < 0.05/9).

Combined Effect of Mono-Synaptic and Multi-Synaptic Pathways

In reality, the motor system contains both the mono-synaptic corticospinal tract and multi-synaptic indirect motor pathways. Hence, to examine the combined effect of both types of descending pathways in a motor system, we performed simulations in a system having dual input to the terminal motoneuron pool. The dual input is comprised of (i) a direct supraspinal drive (resembling the monosynaptic corticospinal tract) and (ii) an indirect drive from the same supraspinal input after being passed through two layers of interneurons (resembling a multi-synaptic descending pathway). The relative weights of these two drives were systematically varied to examine their effects on the neural coupling between the supraspinal input and the motoneuron pool output. Using one-way ANOVA, we found that the proportion of direct vs. indirect drive has a significant effect on the IFC [F_{(5, 594)} = 640.63, p < 0.001], CFC [F_{(5, 594)} = 5552.81, p < 0.001] and COI [F_{(5, 594)} = 7748.80, p < 0.001]. Using Tukey's HSD criterion with Bonferroni correction (k = 5), IFC was found to reduce with increased indirect drive (p < 0.05/5), while the CFC increased (p < 0.05/5). Their combined effect resulted in the progressive increase of COI (p < 0.05/5) (see Figure 7).

FIGURE 7

Figure 7. IFC, CFC, and COI between the supraspinal input and cumulative spike train output of the terminal motoneuron pool for the simulated dual input system. The dual input is comprised of (i) a direct supraspinal drive (resembling the monosynaptic corticospinal tract) and (ii) an indirect drive from the same supraspinal input after being passed through two layers of interneurons (resembling a multi-synaptic descending pathway). The relative proportion (in terms of signal power) of the indirect drive in the composite input was systematically varied from 0 to 100% in steps of 20% increments. The depicted values represent the mean (with ±2 standard deviation as indicated by error bars) calculated from 100 repetitions of the coupling analysis, each being run with a different realization of additive noise to the supraspinal input (as described in section Neural Coupling Analysis). Tukey's test was performed to test for significant decrease in IFC and increase in CFC and COI between the successive steps of increase in indirect drive. To control the family-wise error rate, we set the threshold p = 0.05/k (k = 5, i.e., number of comparisons). Asterisks in superscript of the percentage values denote significant change (**p < 0.05/5).

Discussion

This study investigated neural coupling in descending motor pathways using computer simulations. We simulated the pathways as various layers of interneurons in a cascade with a terminal motoneuron pool, using Hodgkin–Huxley neuron models, to examine the effect of the number of synapses or interneuron layers on the expression of cross-frequency coupling, as well as its ratio over iso-frequency coupling.

Most studies investigating neural coupling in the descending pathways mainly focus on the mono-synaptic corticospinal tract using iso-frequency coupling measures such as cortico-muscular coherence (Mima and Hallett, 1999; Salenius and Hari, 2003; Negro and Farina, 2011b; van Wijk et al., 2012). In this simulation study, we examined both iso- and cross-frequency coupling in the mono-synaptic descending pathway. Our results confirmed the dominance of iso-frequency coupling (as indicated by COI < 0) in the mono-synaptic pathway, though cross-frequency coupling is also present (see Figure 4). This result is in line with our previous experimental work using the NMC to assess cross- and iso-frequency coupling between brain and muscle signals during a low-effort (1 Nm) isotonic wrist flexion in healthy young participants, showing that the motor task using the mono-synaptic corticospinal tract mainly generates iso-frequency coupling (Yang et al., 2016a).

Using linear coherence measure alone, a previous modeling study indicated that the recruitment of multi-synaptic indirect motor pathways can reduce the iso-frequency coupling between the supraspinal input and motoneuron pool output (Negro and Farina, 2011a). Consistent with the previous study, we found that the IFC decreases in the multi-synaptic pathways: the more interneuron layers in the pathway, the smaller IFC between the supraspinal input and the motoneuron pool output. However, we also found that the CFC initially increases in the multi-synaptic pathways and is then followed by a saturation after passing a few neuron layers. The combined effect of changes in IFC and CFC leads to the dominance of cross-frequency coupling (as shown by COI > 0) in the multi-synaptic pathways.

The mechanism underlying the changes in IFC and CFC could be associated with the information distortion that occurs across neuron layers leading to decorrelation of the supraspinal input and the motoneuron pool output (Negro and Farina, 2011a). Such distortion is likely caused by the modulation of inter-spike intervals when the motor command is passing through multiple synaptic layers (Koch and Segev, 2000; Markram, 2003; Yang et al., 2018). This modulation could be attributed to (1) heterogeneous recruitment thresholds and spike after-hyperpolarizations of the individual neurons which results in different firing rates for the same steady-state drive (Powers and Binder, 2001; Heckman and Enoka, 2012; Yang et al., 2018), and (2) the consequent interplay between the time-varying input (as shown Figure 2, which is added over the steady-state drive) and different firing rates (as shown in Figure 3) of the neurons over the entire pool (Thompson et al., 2018). Thus, besides the input frequecies, the neurons also generate responses at other frequencies, which contain the components that are cross-frequency coupled with the input signal, as well as a certain amount of noise that is not phase-locked to the suprapinal input. Not only the cross-frequency coupled components but also the noise can be cumulatively enhanced when the signal is passing from one layer to the next. After passing a few neuron layers, the reduced signal-to-noise ratio then leads to the saturation of CFC.

In the multi-synaptic pathways, a sharp decrease of IFC was found at the terminal motoneuron pool in comparison to the last interneuron layer. This is likely caused by different neuronal processing properties of the simulated interneurons and motoneurons. The motoneurons modeled in this study had an active dendrite with a persistent inward current as well as calcium dependent potassium currents in both the soma and the dendrite compartments (Heckman et al., 2008). In contrast, the interneurons were modeled without such conductances and had a single compartment only. These differences may have given rise to lower IFC (and higher COI) in motoneuron outputs due to their effects on the neurocomputational properties. Indeed, this study opens a broad new area for exploring the origin of different types of neural coupling at the single neuron level and a detailed analysis of the role of individual ionic conductances on IFC and CFC can be the scope of future studies.

The proposed COI measure reflects the dominance of iso-frequency coupling vs. cross-frequency coupling. Interestingly, the iso-frequency coupling is only dominant in the mono-synaptic pathway, while the cross-frequency coupling is dominant in multi-synaptic pathways. After a unilateral brain injury, damage to the mono-synaptic corticospinal tract can increase the reliance on multi-synaptic indirect motor descending pathways (e.g., cortico-reticulospinal tracts for upper limbs) (Owen et al., 2017; McPherson et al., 2018a; Karbasforoushan et al., 2019). The simulated “dual-drive” model mimicks this pathological condition by varing the ratio of multi-synaptic drive vs. mono-synpatic drive. Our results show that the increased input from the indirect drive leads to a more dominant cross-frequency coupling as reflected by an increased value of COI over the increased percentage of the indirect motor drive. This result is in line with our pilot work on eight participants with hemiparetic stroke. The COI between the brain and muscle signal increases when participants with a unilateral stroke progressively lift the weight of their paretic arm (Yang et al., 2019), thereby enhancing the recruitment of indirect motor pathways to compensate for the loss of corticofugal (i.e., corticospinal and corticobulbar) projections from the lesioned hemisphere (McPherson et al., 2018a).

Thus, the COI can be used as a quantitative measure to indicate the relative usage of multi-synaptic indirect motor pathways vs. mono-synaptic direct corticospinal tract. This measure could have a significant impact on future neuro-pathophysiological studies on individuals with an unilateral brain injury, since recent studies have indicated that motor impairments after a unilateral brain injury could be associated with an increased reliance on multi-synaptic indirect motor pathways following a lesion-induced loss of direct corticospinal projections (Owen et al., 2017; McPherson et al., 2018a; Karbasforoushan et al., 2019). Therefore, a measure that quantitatively determines the usage of indirect motor pathways over direct corticospinal drive could be crucial (1) for evaluating motor recovery following unilateral brain injuries, and (2) for determining the effect of targeted therapeutic interventions (Ellis et al., 2018; McPherson et al., 2018b) that aim to reduce the maladaptive reliance on indirect motor pathways after a hemiparetic stroke. In the future, we will examine both cross-frequency and iso-frequency coupling, as well as the COI, between the brain and muscle signals to characterize the relative ratio of the recruitment of indirect vs. direct motor pathways following unilateral brain injuries, such as hemiparetic stroke and unilateral celebral palsy.

Limitations

We acknowledged that there are a few limitations of the current study. First, the interneuron model has only a basic set of ionic conductance since the details of the ionic conductances of spinal and reticular interneurons are yet to be explored. However, such reductionist interneuron models have been used in other simulation studies as well, and this simplification is not expected to change the overall results of this study (Maltenfort et al., 1998; Cisi and Kohn, 2008; Williams and Baker, 2009; Negro and Farina, 2011a). Second, there is no existing literature detailing the connection pattern of individual interneuron layers in multi-synaptic descending pathways. However, the probabilistic model used in this simulation has been previously demonstrated to capture global connectivity properties in motor descending pathways well (Humplik and Tkačik, 2017). Thirdly, the number of neurons in each layer of multi-synaptic descending tracts have not yet been experimentally determined. However, compound EPSP recording on motoneurons from spinal interneurons has shown the number of inputs to be ~ 100 (Lüscher et al., 1983). Hence, it is reasonable to assume that the number of descending inputs on the motoneurons from each tract should be in this order. The motoneuron pool size can vary from muscle to muscle (Karpati et al., 2001), in a range from around 100 (e.g., first dorsal interosseous in humans) (Buchthal and Schmalbruch, 1980) to around 800 (e.g., biceps brachii in humans) (Feinstein et al., 1955). Thus, a size of 100 units for a motoneuron pool falls in the lowest part of the range. In this study, the size of 100 units per neuron layer was adopted also for uniformity and computational convenience. Finally, we did not consider synaptic and transmission delays in this work. The overall delay in a motor pathway may be determined based on the onset latency of Transcranial Magnetic Stimulation (TMS) induced Motor Evoked Potential (MEP) in the targeted muscle (Schwerin et al., 2011). However, the latency of MEP can be affected by coil orientation: difference may exist between direct vs. trans-synaptic activation of the pyramidal cells and the measurement of MEP responses in proximal vs. distal muscles. In short, it is still hard to get a “precise” assessment of the delay in a motor pathway. Meanwhile, there is no available experimental data for the delays in each neuron layer that we can included in this simulation work. Moreover, the time delay only has the effect on the relative phase between signals. This will not result in additional resonance components with new frequencies. Thus, the time delay issue will not affect our current results and overall conclusion of this paper.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

NS conducted the simulation and data analysis under the supervision of YY, JD, and CH. NS and YY drafted and wrote the manuscript. JD and CH commented and revised the manuscript.

Funding

The research leading to these results has received funding support from the support of the Dixon Translational Research Grants Initiative (PI: YY) at Northwestern Medicine and Northwestern University Clinical and Translational Sciences Institute (UL1TR001422), NIH 1R21HD099710 (PIs: JD, YY), R01HD039343 and R01NS105759 (PI: JD).

Conflict of Interest

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

References

Booth, V., Rinzel, J., and Kiehn, O. (1997). Compartmental model of vertebrate motoneurons for Ca2+-dependent spiking and plateau potentials under pharmacological treatment. J. Neurophysiol. 78, 3371–3385. doi: 10.1152/jn.1997.78.6.3371

PubMed Abstract | CrossRef Full Text | Google Scholar

Buchthal, F., and Schmalbruch, H. (1980). Motor unit of mammalian muscle. Physiol. Rev. 60, 90–142. doi: 10.1152/physrev.1980.60.1.90

PubMed Abstract | CrossRef Full Text | Google Scholar

Bui, T. V., Cushing, S., Dewey, D. E., Fyffe, R. E., and Rose, P. (2003). A comparison of the electrotonic and morphological properties of Renshaw cells, Ia inhibitory interneurons and motoneurons in the cat. J. Neurophysiol. 90, 2900–2918. doi: 10.1152/jn.00533.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Cisi, R. R., and Kohn, A. F. (2008). Simulation system of spinal cord motor nuclei and associated nerves and muscles, in a web-based architecture. J. Comput. Neurosci. 25, 520–542. doi: 10.1007/s10827-008-0092-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Daffertshofer, A., Peper, C., Frank, T., and Beek, P. (2000). Spatio-temporal patterns of encephalographic signals during polyrhythmic tapping. Hum. Mov. Sci. 19, 475–498. doi: 10.1016/S0167-9457(00)00032-4

CrossRef Full Text | Google Scholar

Darbon, P., Yvon, C., Legrand, J. C., and Streit, J. (2004). INaP underlies intrinsic spiking and rhythm generation in networks of cultured rat spinal cord neurons. Euro. J. Neurosci. 20, 976–988. doi: 10.1111/j.1460-9568.2004.03565.x

PubMed Abstract | CrossRef Full Text | Google Scholar

De Hemptinne, C., Ryapolova-Webb, E. S., Air, E. L., Garcia, P. A., Miller, K. J., Ojemann, J. G., et al. (2013). Exaggerated phase–amplitude coupling in the primary motor cortex in Parkinson disease. Proc. Natl. Acad. Sci. U.S.A. 110, 4780–4785. doi: 10.1073/pnas.1214546110

PubMed Abstract | CrossRef Full Text | Google Scholar

De Luca, C. J., and Hostage, E. C. (2010). Relationship between firing rate and recruitment threshold of motoneurons in voluntary isometric contractions. J. Neurophysiol. 104, 1034–1046. doi: 10.1152/jn.01018.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Drew, T., Prentice, S., and Schepens, B. (2004). Cortical and brainstem control of locomotion. Prog. Brain Res. 143, 251–261. doi: 10.1016/S0079-6123(03)43025-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Dum, R. P., and Strick, P. L. (1991). The origin of corticospinal projections from the premotor areas in the frontal lobe. J. Neurosci. 11, 667–689. doi: 10.1523/JNEUROSCI.11-03-00667.1991

PubMed Abstract | CrossRef Full Text | Google Scholar

Ellis, M. D., Carmona, C., Drogos, J., and Dewald, J. (2018). Progressive abduction loading therapy with horizontal-plane viscous resistance targeting weakness and flexion synergy to treat upper limb function in chronic hemiparetic stroke: a randomized clinical trial. Front. Neurol. 9:71. doi: 10.3389/fneur.2018.00071

PubMed Abstract | CrossRef Full Text | Google Scholar

Ellis, M. D., Drogos, J., Carmona, C., Keller, T., and Dewald, J. P. (2012). Neck rotation modulates flexion synergy torques, indicating an ipsilateral reticulospinal source for impairment in stroke. J. Neurophysiol. 108, 3096–3104. doi: 10.1152/jn.01030.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Ellis, M. D., Schut, I., and Dewald, J. P. (2017). Flexion synergy overshadows flexor spasticity during reaching in chronic moderate to severe hemiparetic stroke. Clin. Neurophysiol. 128, 1308–1314. doi: 10.1016/j.clinph.2017.04.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Feinstein, B., Lindegård, B., Nyman, E., and Wohlfart, G. (1955). Morphologic studies of motor units in normal human muscles. Cells Tissues Organs 23, 127–142. doi: 10.1159/000140989

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferrario, A., Merrison-Hort, R., Soffe, S. R., and Borisyuk, R. (2018). Structural and functional properties of a probabilistic model of neuronal connectivity in a simple locomotor network. Elife 7:e33281. doi: 10.7554/eLife.33281.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Finkel, A. S., and Redman, S. J. (1983). The synaptic current evoked in cat spinal motoneurones by impulses in single group 1a axons. J. Physiol. 342, 615–632. doi: 10.1113/jphysiol.1983.sp014872

PubMed Abstract | CrossRef Full Text | Google Scholar

Fries, W., Danek, A., Scheidtmann, K., and Hamburger, C. (1993). Motor recovery following capsular stroke: role of descending pathways from multiple motor areas. Brain 116, 369–382. doi: 10.1093/brain/116.2.369

PubMed Abstract | CrossRef Full Text | Google Scholar

Hagihira, S., Takashina, M., Mori, T., Mashimo, T., and Yoshiya, I. (2001). Practical issues in bispectral analysis of electroencephalographic signals. Anesth. Analg. 93, 966–970. doi: 10.1097/00000539-200110000-00032

PubMed Abstract | CrossRef Full Text | Google Scholar

Heckman, C., Johnson, M., Mottram, C., and Schuster, J. (2008). Persistent inward currents in spinal motoneurons and their influence on human motoneuron firing patterns. Neuroscientist 14, 264–275. doi: 10.1177/1073858408314986

PubMed Abstract | CrossRef Full Text | Google Scholar

Heckman, C. J., and Enoka, R. M. (2012). Motor unit. Compr. Physiol. 2, 2629–2682. doi: 10.1002/cphy.c100087

PubMed Abstract | CrossRef Full Text | Google Scholar

Humplik, J., and Tkačik, G. (2017). Probabilistic models for neural populations that naturally capture global coupling and criticality. PLoS Comput. Biol. 13:e1005763. doi: 10.1371/journal.pcbi.1005763

PubMed Abstract | CrossRef Full Text | Google Scholar

Jang, S. H., Chang, C. H., Lee, J., Kim, C. S., Seo, J. P., and Yeo, S. S. (2013). Functional role of the corticoreticular pathway in chronic stroke patients. Stroke 44, 1099–1104. doi: 10.1161/STROKEAHA.111.000269

PubMed Abstract | CrossRef Full Text | Google Scholar

Jang, S. H., and Seo, J. P. (2014). The distribution of the cortical origin of the corticoreticular pathway in the human brain: a diffusion tensor imaging study. Somatosens. Mot. Res. 31, 204–208. doi: 10.3109/08990220.2014.917292

PubMed Abstract | CrossRef Full Text | Google Scholar

Karbasforoushan, H., Cohen-Adad, J., and Dewald, J. (2019). Brainstem and spinal cord MRI identifies altered sensorimotor pathways post-stroke. Nat. Commun. 10:3524. doi: 10.1038/s41467-019-11244-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Karpati, G., Hilton-Jones, D., and Griggs, R. C., (eds). (2001). Disorders of Voluntary Muscle., 7th edition Cambridge: Cambridge University Press.

Koch, C., and Segev, I. (2000). The role of single neurons in information processing. Nat. Neurosci. 3:1171. doi: 10.1038/81444

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, R., and Heckman, C. (2001). Essential role of a fast persistent inward current in action potential initiation and control of rhythmic firing. J. Neurophysiol. 85, 472–475. doi: 10.1152/jn.2001.85.1.472

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, S., Chen, Y.-T., Francisco, G., Zhou, P., and Rymer, W. Z. (2019). A unifying pathophysiological account for post-stroke spasticity and disordered motor control. Front. Neurol. 10:468. doi: 10.3389/fneur.2019.00468

PubMed Abstract | CrossRef Full Text | Google Scholar

Lüscher, H., Ruenzel, P., and Henneman, E. (1983). Composite EPSPs in motoneurons of different sizes before and during PTP: implications for transmission failure and its relief in Ia projections. J. Neurophysiol. 49, 269–289. doi: 10.1152/jn.1983.49.1.269

PubMed Abstract | CrossRef Full Text | Google Scholar

Maltenfort, M. G., Heckman, C., and Rymer, W. Z. (1998). Decorrelating actions of renshaw interneurons on the firing of spinal motoneurons within a motor nucleus: a simulation study. J. Neurophysiol. 80, 309–323. doi: 10.1152/jn.1998.80.1.309

PubMed Abstract | CrossRef Full Text | Google Scholar

Markram, H. (2003). Elementary principles of nonlinear synaptic transmission. Computational Models for Neuroscience (Springer), 125–169. doi: 10.1007/978-1-4471-0085-0_5

CrossRef Full Text | Google Scholar

Matsuyama, K., Mori, F., Nakajima, K., Drew, T., Aoki, M., and Mori, S. (2004). Locomotor role of the corticoreticular–reticulospinal–spinal interneuronal system. Prog. Brain Res. 143, 239–249. doi: 10.1016/S0079-6123(03)43024-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Mattei, B., and Schmied, A. (2002). Delayed and prolonged effects of a near threshold EPSP on the firing time of human α-motoneurones. J. Physiol. 538, 849–865. doi: 10.1113/jphysiol.2001.012701

PubMed Abstract | CrossRef Full Text | Google Scholar

McDougal, R. A., Morse, T. M., Carnevale, T., Marenco, L., Wang, R., Migliore, M., et al. (2017). Twenty years of ModelDB and beyond: building essential modeling tools for the future of neuroscience. J. Comput. Neurosci. 42, 1–10. doi: 10.1007/s10827-016-0623-7

PubMed Abstract | CrossRef Full Text | Google Scholar

McPherson, J. G., Chen, A., Ellis, M. D., Yao, J., Heckman, C., and Dewald, J. P. (2018a). Progressive recruitment of contralesional cortico-reticulospinal pathways drives motor impairment post stroke. J. Physiol. 596, 1211–1225. doi: 10.1113/JP274968

PubMed Abstract | CrossRef Full Text | Google Scholar

McPherson, J. G., Ellis, M. D., Harden, R. N., Carmona, C., Drogos, J. M., Heckman, C. J., et al. (2018b). Neuromodulatory inputs to motoneurons contribute to the loss of independent joint control in chronic moderate to severe hemiparetic stroke. Front. Neurol. 9:470. doi: 10.3389/fneur.2018.00470

PubMed Abstract | CrossRef Full Text | Google Scholar

McPherson, J. G., Stienen, A. H., Drogos, J. M., and Dewald, J. P. (2018c). Modification of spastic stretch reflexes at the elbow by flexion synergy expression in individuals with chronic hemiparetic stroke. Arch. Phys. Med. Rehabil. 99, 491–500. doi: 10.1016/j.apmr.2017.06.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Mima, T., and Hallett, M. (1999). Corticomuscular coherence: a review. J. Clin. Neurophysiol. 16:501. doi: 10.1097/00004691-199911000-00002

PubMed Abstract | CrossRef Full Text | Google Scholar

Mima, T., Toma, K., Koshy, B., and Hallett, M. (2001). Coherence between cortical and muscular activities after subcortical stroke. Stroke 32, 2597–2601. doi: 10.1161/hs1101.098764

PubMed Abstract | CrossRef Full Text | Google Scholar

Negro, F., and Farina, D. (2011a). Decorrelation of cortical inputs and motoneuron output. J. Neurophysiol. 106, 2688–2697. doi: 10.1152/jn.00336.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Negro, F., and Farina, D. (2011b). Linear transmission of cortical oscillations to the neural drive to muscles is mediated by common projections to populations of motoneurons in humans. J. Physiol. 589, 629–637. doi: 10.1113/jphysiol.2010.202473

PubMed Abstract | CrossRef Full Text | Google Scholar

Nikias, C. L., and Mendel, J. M. (1993). Signal processing with higher-order spectra. IEEE Signal Process. Mag. 10, 10–37. doi: 10.1109/79.221324

CrossRef Full Text | Google Scholar

Owen, M., Ingo, C., and Dewald, J. (2017). Upper extremity motor impairments and microstructural changes in bulbospinal pathways in chronic hemiparetic stroke. Front. Neurol. 8:257. doi: 10.3389/fneur.2017.00257

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfurtscheller, G., and Da Silva, F. L. (1999). Event-related EEG/MEG synchronization and desynchronization: basic principles. Clin. Neurophysiol. 110, 1842–1857. doi: 10.1016/S1388-2457(99)00141-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Powers, R. K., and Binder, M. D. (2001). Input-output functions of mammalian motoneurons. Rev. Physiol. Biochem. Pharmacol. 143, 137–263. doi: 10.1007/BFb0115594

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Rybak, I. A., Shevtsova, N. A., Lafreniere-Roula, M., and Mccrea, D. A. (2006). Modelling spinal circuitry involved in locomotor pattern generation: insights from deletions during fictive locomotion. J. Physiol. 577, 617–639. doi: 10.1113/jphysiol.2006.118703

PubMed Abstract | CrossRef Full Text | Google Scholar

Salenius, S., and Hari, R. (2003). Synchronous cortical oscillatory activity during motor action. Curr. Opin. Neurobiol. 13, 678–684. doi: 10.1016/j.conb.2003.10.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Schoffelen, J.-M., Oostenveld, R., and Fries, P. (2005). Neuronal coherence as a mechanism of effective corticospinal interaction. Science 308, 111–113. doi: 10.1126/science.1107027

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwerin, S. C., Yao, J., and Dewald, J. P. (2011). Using paired pulse TMS to facilitate contralateral and ipsilateral MEPs in upper extremity muscles of chronic hemiparetic stroke patients. J. Neurosci. Methods 195, 151–160. doi: 10.1016/j.jneumeth.2010.11.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Streit, J., Tscherter, A., and Darbon, P. (2006). Rhythm generation in spinal cultures: Is it the neuron or the network? Adv. Network Electrophysiol. 377–408. doi: 10.1007/0-387-25858-2_16

CrossRef Full Text | Google Scholar

Sturm, H., Schmied, A., Vedel, J. P., and Pagni, S. (1997). Firing pattern of type-identified wrist extensor motor units during wrist extension and hand clenching in humans. J. Physiol. 504, 735–745. doi: 10.1111/j.1469-7793.1997.735bd.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanji, J., and Kato, M. (1972). Discharges of single motor units at voluntary contraction of abductor digiti minimi muscle in man. Brain Res. 45, 590–593. doi: 10.1016/0006-8993(72)90488-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Thompson, C. K., Negro, F., Johnson, M. D., Holmes, M. R., McPherson, L. M., Powers, R. K., et al. (2018). Robust and accurate decoding of motoneuron behaviour and prediction of the resulting force output. J. Physiol. 596, 2643–2659. doi: 10.1113/JP276153

PubMed Abstract | CrossRef Full Text | Google Scholar

Tian, R., Yang, Y., Van Der Helm, F. C., and Dewald, J. P. A. (2018). A novel approach for modeling neural responses to joint perturbations using the NARMAX method and a hierarchical neural network. Front. Comput. Neurosci. 12:96. doi: 10.3389/fncom.2018.00096

PubMed Abstract | CrossRef Full Text | Google Scholar

van Wijk, B. C., Beek, P. J., and Daffertshofer, A. (2012). Neural synchrony within the motor system: what have we learned so far? Front. Hum. Neurosci. 6:252. doi: 10.3389/fnhum.2012.00252

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, E. R., and Baker, S. N. (2009). Renshaw cell recurrent inhibition improves physiological tremor by reducing corticomuscular coupling at 10 Hz. J. Neurosci. 29, 6616–6624. doi: 10.1523/JNEUROSCI.0272-09.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Dewald, J., Van Der Helm, F. C., and Schouten, A. C. (2018). Unveiling neural coupling within the sensorimotor system: directionality and nonlinearity. Euro. J. Neurosci. 48, 2407–2415. doi: 10.1111/ejn.13692

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Schouten, A. C., Solis-Escalante, T., and Van Der Helm, F. C. (2015). Probing the nonlinearity in neural systems using cross-frequency coherence framework. IFAC-papersonline 48, 1386–1390. doi: 10.1016/j.ifacol.2015.12.326

CrossRef Full Text | Google Scholar

Yang, Y., Solis-Escalante, T., Van De Ruit, M., Van Der Helm, F. C., and Schouten, A. C. (2016a). Nonlinear coupling between cortical oscillations and muscle activity during isotonic wrist flexion. Front. Comput. Neurosci. 10:126. doi: 10.3389/fncom.2016.00126

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Solis-Escalante, T., Van Der Helm, F. C., and Schouten, A. C. (2016b). A generalized coherence framework for detecting and characterizing nonlinear interactions in the nervous system. IEEE Trans. Biomed. Eng. 63, 2629–2637. doi: 10.1109/TBME.2016.2585097

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Tian, R., Cummings, M., Hogenhuis, C. H. E. Y., Deol, J., Drogos, J., et al. (2019). Determine the Usage of Indirect Motor Pathways Following a Hemiparetic Stroke: A Pilot Study. Philadelphia, PA: Neuroscience 2019.

Appendix

Equilibrium Potentials

E_Na = 55mV; E_K = −80mV; E_Ca = 80mV;

E_L(S) = −65mV ±0.3; E_L(D) = −65 ±0.15mV (motoneurons*)

*Motoneuron equilibrium potentials were assigned from uniform random distributions with mean ± S.D as given above to obtain firing rates in the range of 5–30 spikes/s as described in section Simulations.

**Interneuron equilibrium potentials were set as described in section Simulations.

Neuron Paramteres

Motoneurons:

$g_{Na(S)}=120mSc{m}^{-2};g_{NaP(S)}=0.1mSc{m}^{-2};g_{K(S)}=100mSc{m}^{-2};$

$g_{CaN(S)}=14mSc{m}^{-2};g_{K,Ca(S)}=5mSc{m}^{-2};g_{L(S)}=0.51mSc{m}^{-2};$

$g_{CaN(D)}=0.3mSc{m}^{-2};g_{CaL(D)}=0.33mSc{m}^{-2};g_{NaP(D)}=0.1mSc{m}^{-2};$

$g_{K,Ca(D)}=1.1mSc{m}^{-2};g_{L(D)}=0.51mSc{m}^{-2};g_c=0.1mSc{m}^{-2};$

$p=0.1;f=0.01;\alpha =0.0009mol{C}^{-1}\mu m^{-1};k_{Ca}=2m{s}^{-1};K_d=0.2\mu M$

Interneurons:

$g_{Na}=120mSc{m}^{-2};g_K=100mSc{m}^{-2};g_L=0.51mSc{m}^{-2};$

Synapses:

$E_{synE}=-10mV;g_{synE}=0.01mSc{m}^{-2};{\tau }_{synE}=5ms;$

$E_{synE}=-10mV;g_{synE}=0.0075mSc{m}^{-2};{\tau }_{synE}=5ms;$

TABLE A1

TABLE A1. Steady-state activation and inactivation variables and time constants for voltage-dependent ionic channels (Rybak et al., 2006).