Abstract
Oscillations in the beta and gamma bands (13–30 Hz; 35–70 Hz) have often been observed in motor cortical outputs that reach the spinal cord, acting on motoneurons and interneurons. However, the frequencies of these oscillations are above the muscle force frequency range. A current view is that the transformation of the motoneuron pool inputs into force is linear. For this reason possible roles for these oscillations are unclear, since if this transformation is linear, the high frequencies in the motoneuron inputs (e.g., 20 Hz from pyramidal tract neurons) would be filtered out by the muscle and have no effect on force control. A biologically inspired mathematical model of the neuromuscular system was used to investigate the impact of high-frequency cortical oscillatory activity on force control. The model simulation results evidenced that a typical motoneuron pool has a nonlinear behavior that enables the decoding of a high-frequency oscillatory input. An input at a single frequency (e.g., beta band) leads to an increase in the steady-state force generated by the muscle. When the input oscillation was amplitude modulated at a given low frequency, the force oscillated at this frequency. In both cases, the mechanism relies on the recruitment and derecruitment of motor units in response to the oscillatory descending drive. Therefore, the results from this study suggest a potential role in force control for cortical oscillations at frequencies at or above the beta band, despite the low-pass behavior of the muscles.
SIGNIFICANCE STATEMENT The role of cortical oscillations in motor control has been a long-standing question, one view being that they are an epiphenomenon. Fast oscillations are known to reach the spinal cord, and hence they have been thought to affect muscle behavior. However, experimental limitations have hampered further advances to explain how they could influence muscle force. An approach for such a challenge was adopted in the present research: to study the problem through computer simulations of an advanced biologically compatible mathematical model. Using such a model, we found that the well-known mechanism of recruitment and derecruitment of the spinal cord motoneurons can allow the muscle to respond to cortical oscillations, suggesting that these oscillations are not epiphenomena in motor control.
Introduction
Force generation by a muscle is dependent on the patterns of action potentials discharged by the respective α-motoneurons (MNs) as well as the biomechanical features of the muscle and the surrounding structures. The firing patterns of the MNs, on the other hand, are determined by several inputs, originating from different sources (Watanabe et al., 2013; Baudry and Duchateau, 2014), such as afferent signals and cortical descending commands, which include oscillatory activity. However, the most important input sources for an MN pool are the common synaptic inputs (CSIs), since they define the time-varying force generated by the muscle (De Luca and Erim, 1994; Farina et al., 2014).
Cortical oscillations in the beta (13–30 Hz) and gamma bands (35–70 Hz), have been observed in monkeys (Murthy and Fetz, 1992, 1996; Sanes and Donoghue, 1993; Baker et al., 1997, 2001, 2003) and humans (Conway et al., 1995; Schoffelen et al., 2005; Kristeva et al., 2007; Engel and Fries, 2010; Gwin and Ferris, 2012). Cortical oscillations are coupled in the beta and gamma bands to the MN pool activity through pyramidal tract axons (Baker et al., 2003). This coupling has been studied by indirect means in different motor tasks by the coherence between the EEG (or MEG or local field potentials) and the EMG (Conway et al., 1995; Baker et al., 1997; Salenius et al., 1997; Halliday et al., 1998; Mima et al., 1999; Kilner et al., 2000; Gilbertson et al., 2005; Schoffelen et al., 2005; Andrykiewicz et al., 2007; Kristeva et al., 2007; Chakarov et al., 2009; Naranjo et al., 2010; Yang et al., 2010; Gwin and Ferris, 2012; Mendez-Balbuena et al., 2012; Ushiyama et al., 2012). The significant coherence reported at beta and/or gamma frequencies indicates that muscular activity at these frequency bands is dependent, at least in part, on corticospinal activity. Despite the large number of experimental observations regarding cortical oscillations in different motor tasks, their functional roles are not clear (Fetz, 2013; Aumann and Prut, 2015).
Due to unsurmountable experimental difficulties to measure the inputs to MN pools in humans (Heckman and Enoka, 2012), the study of how the CSIs influence the MN spike trains was approached previously by means of computer model simulations (Williams and Baker, 2009a,b; Stegeman et al., 2010; Negro and Farina, 2011a,b, 2012; Watanabe et al., 2013; Farina et al., 2014; Gallego et al., 2015). This approach was also used to study oscillatory CSIs in general neuron populations (Knight, 1972; Fourcaud-trocme et al., 2003), and specifically MN populations (Lowery et al., 2007; Williams and Baker, 2009b; Stegeman et al., 2010; Negro and Farina, 2011a; Farina et al., 2014; Dideriksen et al., 2015; Gallego et al., 2015).
Previous research (Lowery et al., 2007; Stegeman et al., 2010; Negro and Farina, 2011a; Farina et al., 2014) proposed that a pool of MNs would act as a linear system, and hence the force output would mirror the time-varying CSIs of the MN pool for frequencies up to ∼6–10 Hz.
The above-mentioned experimentally obtained corticomuscular coherence in the beta (and/or gamma) bands and the postulate that the MN pool behaves linearly are apparently inconsistent. This is so because, according to the latter conjecture, the beta (or higher) band cortical output oscillations that are part of the CSIs to an MN pool would be filtered out by the muscle and hence contribute nothing to muscle force control (Farina et al., 2014).
To explore the issue of how beta (or higher) band CSI oscillations can influence force production, the present study relies on a biologically compatible neuromuscular computational model (Cisi and Kohn, 2008; Elias et al., 2012, 2014; Watanabe et al., 2013). The results from the simulation of this computational model indicate that the recruitment properties of a pool of MNs enable cortical oscillations at beta or higher frequencies to have a role in force control.
Materials and Methods
Biologically compatible model.
The effects of oscillatory CSIs at the beta or higher frequency ranges on the spiking activity of an MN pool were analyzed in the present study by means of simulations of a previously reported and validated biologically compatible neuromuscular model (Cisi and Kohn, 2008; Elias et al., 2012, 2014; Watanabe et al., 2013). Figure 1 shows block diagrams of the model at different levels of detail. The simulations of this model provide membrane potentials and spike times of all MNs, MN synaptic conductances, and muscle force. All parameter values were the same as those used by Watanabe et al. (2013).
Structure of the neuromuscular model used in this study. A, Schematic view of the overall model. The motor unit pool encompasses MNs and the respective MU. B, Schematic representation of muscle-specific motor units (MNs and MUs). Each of the 400 independent renewal point processes representing the descending axons activate synaptically ∼30% of the MNs. Besides the descending axons, each MN receives an IN input in the form of a Poisson point process. Each MU receives a train of action potentials from the respective MN and generates twitches that pass through a saturation nonlinearity and are summed with the contributions from the other MUs to produce the muscle force, F(t). C, Train of renewal point processes representing a train of spikes from a descending axon. D, Equivalent circuit used to represent each MN model. gsyn1 to gsynM, synaptic conductances of a dendrite activated by synapses 1 to M, respectively; gc, coupling conductance; gLd and gLs, dendritic and somatic leakage conductances, respectively; gNa, gKf, and gKs, conductances of Na+, fast K+, and slow K+, respectively; EL, leakage potential; ENa and EK, Na+ and K+ equilibrium potentials, respectively; Esyn1 to EsynM, reversal potentials for synapses 1 to M, respectively; Cs and CD, somatic and dendritic capacitances, respectively; VS and VD, somatic and dendritic membrane potentials, respectively. E, Spike trains of the MNs of the pool. F, force produced by a single MU.
The simulation results shown here are for a single MN pool parameterized for the soleus motor unit pool. However, it should be added that the results do not depend on the specifically chosen MN pool. Each of 400 descending axons connect to ∼30% of the MNs, randomly chosen for each simulation run. (Each MN will roughly receive a similar number of synaptic inputs.) Each motor unit encompasses an MN and a muscle unit (MU; Fig. 1B). The spike trains of descending axons were modeled as independent renewal point processes with interspike intervals (ISIs) following a Poisson distribution (Fig. 1C). Besides the descending axons, each MN receives a different independent noise (IN) source input in the form of a Poisson point process, with the ISIs following a Poisson distribution. The MN models are conductance based, with two compartments to represent the MN soma and the dendritic tree (Fig. 1D). The MN pool used in this study is composed of 800 slow MNs, 50 fast-fatigue-resistant MNs, and 50 fast-fatiguing MNs (Watanabe et al., 2013). All the synapses (from the CSIs and the IN inputs) on a given MN act through changes in the dendritic conductances (Watanabe et al., 2013). The net dendritic synaptic conductance of a randomly chosen MN provides an estimate of the overall CSIs to the MNs since, probabilistically, the MNs have similar net synaptic conductance signals. The force produced by a given MU is obtained by the convolution of the respective MN spikes (Fig. 1E) with the impulse response of a second-order critically damped system (Fig. 1F). This force signal is then passed through a nonlinear saturation (Watanabe et al., 2013) and summed to similarly obtained force signals from the other MUs, resulting in the muscle force F(t). All the differential equations were solved using a fourth-order Runge–Kutta method with fixed step of 0.05 ms. More details can be found in the studies by Cisi and Kohn (2008), Watanabe et al. (2013), and Elias et al. (2014).
Simulation protocols.
In the simulations, the mean basal firing rate of the descending axons was set to 65 spikes/s, and independent stochastic synaptic stimulation was added to each MN in the form of Poisson point processes with a mean ISI of 8 ms to produce a simulated contraction equivalent to ∼10% of the maximal voluntary contraction (MVC). These values of mean basal firing rate of the descending axons and ISI of the INs lead to an average conductance magnitude caused by the CSIs of one-third of the average conductance magnitude caused by the IN, similar to that used by Williams and Baker (2009a,b). All the inputs to the MNs were excitatory, although the presence of inhibitory inputs did not change the general conclusions of this study.
Two simulation protocols were used. The first aimed to simulate steady force conditions, and the second was used to simulate periodic contractions. For both protocols, in the first 60 s, the descending axons fired stochastically at 65 spikes/s. After 60 s, both protocols differed in the following manner.
The first protocol used a pure 20 Hz sinusoidal signal, with two different mean firing rates:
where FR(t) is the mean firing rate of each descending axon. The resulting signal FR(t) has a duration equal to 180 s and can be seen in Figure 2A.
Mean firing rate of each descending axon as used in the simulation protocols. The mean basal firing rate for both cases was 65 spikes/s. A, The firing rate modulation began after 60 s with a 20 Hz sinusoidal signal of amplitude 20 spikes/s. After 120 s, the mean firing rate decresased to 58 spikes/s. B, The firing rate modulation began after 60 s with a 20 Hz sinusoidal signal, which was amplitude modulated by a sinusoidal signal of 1 Hz.
In the second case,
in which the second term corresponds to what is known in communication theory as amplitude-modulated signal (AM; having a carrier frequency equal to 20 Hz; Carlson, 1975). The resulting signal has a duration of 120 s and is shown in Figure 2B. For both the pure 20 Hz oscillation and the AM oscillation, a populational synchronization measure proposed by Baker et al. (1999) was computed for several simulation runs, resulting in a mode value of 0.3, similar to the experimental data from the study by Baker et al. (1999).
The frequency of 20 Hz was chosen to simulate beta-band oscillations coming from the motor cortex and reaching the spinal cord (Conway et al., 1995; Kilner et al., 2000; Baker et al., 2003), but it could be any other frequency, e.g., in the gamma band or higher (Buzsáki and Draguhn, 2004). The amplitude-modulated signal was motivated by our observation that an increase in the beta oscillation amplitude led to an increase in force level.
Cumulative spike train.
The cumulative spike train (CST) is the summation (superposition) of the MN spike trains (Farina et al., 2014).
After the summation of the MN spike trains, each spike is transformed into a pulse of 0.05 ms duration and amplitude 20000 (i.e., the area of each rectangular pulse was equal to 1), for the purpose of computing the CST spectrum.
One hundred randomly chosen combinations of five MNs of a given simulation run were used to generate 100 different CSTs, which were used to compute their spectra and the coherence between MN inputs and the CST (see subsections Spectral estimation and Coherence, below).
Activation ratio.
The concept of activation ratio (AR; Héroux et al., 2014) was used to quantify differences of the fraction of time a given MN was active. This ratio was computed as the percentage of time the ISIs of a given MN were lower than 100 ms. This means that if an MN has ISIs >100 ms (<10 spikes/s) during the whole simulation, its AR is equal to 0. At the other end, if an MN fires >10 spikes/s over the simulation period (ISIs lower than 100 ms), its AR is equal to 1. Bursty activity of an MN corresponds to an AR value between 0 and 1.
Spectral estimation.
Spectra and cross-spectra of the synaptic conductance, muscle force, and the CST were estimated by using Welch's method (Bendat and Piersol, 2010). The first 3 s of each signal after time 0 and after changes of modulation type were removed to avoid transient effects. Therefore, 57 s of each signal were used for the spectral and cross-spectral computations. The number of FFT points was set to 60,000, a Hamming window was used with no overlap, and, as a consequence, the spectral resolution was 0.33 Hz. Before spectral and cross-spectral computations, all signals had their linear trends removed.
Coherence.
The corticomuscular coherence was computed by using the net dendritic synaptic conductance as an estimate of the MN inputs, as in the study by Williams and Baker, (2009a,b), and the CST of the MN spikes as the motor output. The coherence between the MN inputs (signal x) and the CST (signal y) was computed using the magnitude squared coherence (Bendat and Piersol, 2010):
where Sxy(f) is the cross-spectrum between signals x and y, and Sx(f) and Sy(f) are the spectra of individual signals x and y, respectively.
The 95% confidence level (CL) for a zero coherence hypothesis was computed by using the following statistics (Brockwell and Davis, 1991):
where K is the number of segments used in Welch's method, α = 0.05, and Fn1,n2 is the F-distribution with n1 and n2 degrees of freedom.
Further simulations.
To test the robustness and wider generality of the results obtained in this study, additional simulations were executed using a few different parameter values, two different type of descending drives, and two model modifications.
Instead of the Poisson point processes, gamma point processes were used to represent the firing behavior of the descending axons. This different point process statistic was used to analyze the robustness of the influence of the beta-band CSIs in muscle force control to the change of premotoneuronal command statistics. The shape factor of the gamma point process was 7, as in the study by Watanabe et al. (2013).
The first model modification was the addition of Renshaw cells to analyze whether the recurrent inhibition of the MNs changes how beta-band CSIs affect muscle force production. The Renshaw cell model used was identical to that described by Cisi and Kohn (2008), and the connectivity among Renshaw cell and MNs followed that used by Williams and Baker (2009b).
The second model modification was the use of active dendrites instead of passive dendrites, this being implemented by adding a representation of L-type Ca2+ channels into the dendritic compartment. This modification was tested to investigate whether the presence of persistent inward currents affects the influence of beta-band oscillations of the CSIs on muscle force control. The model used was that described in the studies by Elias et al. (2012) and Elias and Kohn (2013).
As an alternative to the AM descending drive to achieve a time-varying force, a 20 Hz oscillation in the mean rate of the input driving the MN pool was added to a low-frequency sinusoid of 1 Hz to achieve a 1 Hz force variation [FR(t) = 0.2 sin(2πt) + 20 sin(2π20t)]. This low-frequency input is the most commonly used in previous simulation studies (Farina et al., 2014; Gallego et al., 2015) of muscle force control.
Implementation details.
The biologically compatible neuromuscular computational model described above was written in the Java programming language (Oracle) and simulated in an open-source Web-based application developed in the Eclipse integrated development environment (The Eclipse Foundation). It may be downloaded from a public repository (http://dx.doi.org/10.6084/m9.figshare.1486375). All the data processing was performed in the MATLAB environment (MathWorks).
Results
Steady force protocol
Figure 3A–C shows the spectra of the input to the MN pool for the three situations in Figure 2A. As expected (Sanderson, 1980), a peak at 20 Hz appeared in Figure 3, B and C, due to the oscillatory CSIs at this frequency. Figure 3D–F shows the coherence between the synaptic inputs and the CST (Farina et al., 2014) from five randomly chosen MNs with the same MN pool inputs as in Figure 3A–C, respectively. The coherence values for the constant mean firing rate input (dark blue line) were not statistically different from zero, whereas the coherence values of ∼20 Hz were significantly above zero for the oscillatory input cases (Fig. 3D–F, dark brown and green lines). Figure 3G shows the force signal during the 180 s of simulation. When the input mean intensity to the MN pool was 65 spikes/s, the muscle produced a force of ∼375 N (in blue). Starting at 60 s, the mean input firing rate oscillated at 20 Hz around 65 spikes/s, and the resulting muscle force increased to ∼475 N (in brown). Starting at 120 s, the mean input firing rate oscillated at 20 Hz, but around 58 spikes/s, and the muscle force returned approximately to the initial force of 375 N (in green). Figure 3H shows the corresponding raster plots (aligned with the plots in G) of the MN spike trains during the different inputs; MNs are numbered according to the size principle (Henneman, 1957). It is noteworthy that there was a recruitment of higher twitch amplitude motor units between 60 and 120 s, caused by the 20 Hz oscillation of the input firing rate to the MN pool around 65 spikes/s (brown dots). This recruitment of higher twitch amplitude motor units occurred during the positive part of the 20 Hz cycle and caused the increase of the muscle force between 60 and 120 s. This can be better seen in the insets of Figure 3H, which also show that the lower-threshold MNs had their firing rates modulated by the 20 Hz input oscillation. An interesting observation is that the mean force levels were similar in the first and third time intervals (Fig. 3G, blue and green), even though the mean firing rate at the inputs was smaller in the third segment (58 spikes/s).
Results from the simulation using the firing rate of the descending command shown in Figure 2A. A–C, Net dendritic synaptic conductance power spectrum with the mean basal firing rate of 65 spikes/s and no oscillation in the input (A), with the mean basal firing rate 65 spikes/s and 20 Hz oscillation in the input (B), and with 20 Hz oscillation and the mean firing rate decreased to 58 spikes/s (C). D–F, Coherence between the synaptic conductance and the CST of five randomly chosen MNs with the same MN pool inputs as in A–C, respectively. Individual data from the 100 randomly chosen CSTs are in light color, and the overall mean is in dark color. The 95% confidence level is shown as a dashed line. G, Force signal during the simulation. H, Corresponding raster plot of the MN spike trains during the simulation. The MNs are ordered according to the size principle (Henneman, 1957). The insets show 100 ms windows of the raster plot.
Periodic force protocol
Figure 4, A and B, shows the input spectra for constant and oscillatory firing rate inputs, respectively. In the latter case (red line), three peaks appeared in the input spectrum, at 19, 20, and 21 Hz (Fig. 4B, inset). Figure 4, C and D, shows the coherence between the synaptic inputs and the CST computed from five randomly chosen MNs with the same MN pool inputs of Figure 4A and B, respectively. The coherence values for the constant firing rate mean input (dark blue line) were not statistically different from zero, whereas the coherence values around 20 Hz were statistically significant for the amplitude-modulated input (dark red line). The force spectra resulting from each of the inputs are shown in Figure 4, E and F. For the amplitude-modulated input, the force spectrum shows a component at 1 Hz, which was not present in the CSI (red line).
Spectra and coherences of inputs and outputs when the firing rate of the descending command is as shown in Figure 2B. A, B, Net synaptic conductance spectrum with no modulation in the input (A) and with the modulated input (B). C, D, Coherence between synaptic conductance and the CST of five randomly chosen MNs of the pool with no oscillation in the input (C) and with the modulated input (D). Individual data from 100 randomly chosen CSTs are in light color, and the mean is in dark color. The 95% confidence level is shown as a dashed line. E, F, Force spectrum with no modulation in the input (E) and with the modulated input (F).
The component at 1 Hz corresponds in the time domain to the 1 Hz muscle force oscillation starting at 60 s, shown in Figure 5A (red). In Figure 5B, the spike times of all MNs of the pool are shown, ordered according to the size principle (Henneman, 1957). The colors of the dots change according to the AR values of each MN (see color bars). When the input signal was amplitude-modulated, the MNs with the lowest AR (the darker red points) were recruited and derecruited at exactly 1 Hz, thereby generating the muscle force oscillation. The darker thick lines in Figure 5, C and D (averages of the corresponding individual spectra shown in thin lines) are the mean spectra of the CST of five MNs randomly chosen from 5% of those having the lowest AR (firing during a lower fraction of time). The lighter thick lines in Figure 5, C and D, are the mean spectra of the CST of five randomly chosen MNs from 5% of those with the highest AR (firing during a higher fraction of time). The spectral peak for the high-AR MNs (Fig. 5C) when the input mean rate was constant (i.e., the MN pool was driven by time-invariant point processes) indicates that the mean firing rate of the first recruited MNs is ∼16 spikes/s. For the amplitude-modulated input, the CST originating from the low-AR MNs (Fig. 5D, red line) has a significantly higher power spectrum peak at 1 Hz than that of the CST formed by high AR MNs (yellow line). On the other hand, the 20 Hz spectral component appears for both low- and high-AR MNs.
Results from the simulation using the firing rate of the descending command shown in Figure 2B. A, Muscle force signal during a 15 s window. When the synaptic inputs were modulated (starting at 60 s), the muscle force began to oscillate at 1 Hz. B, Raster plot of the MN spike trains during the 15 s window. The dot colors specify the values of the MN ARs, which are indicated by the color bars. C, D, Darker line, Mean spectrum of the CST of five randomly chosen MNs from the 5% of the MNs from the pool with the lowest AR; lighter line, mean spectrum of the CST of five randomly chosen MNs from the 5% of the MNs from the pool with the highest AR. The spectra were normalized by the maximal power of the each spectrum (C) and by the power at 20 Hz (D). Individual data from the 100 randomly chosen CSTs are in light color.
Further simulations
The choices regarding oscillation amplitude and frequency were in part arbitrary. For this reason, additional simulations were performed to evaluate the robustness of the extraction of low-frequency information from higher-frequency oscillation.
Figure 6 shows the results of simulations for different amplitudes of the firing rate sinusoidal oscillation at the MN pool input, evidencing that the effect of the oscillations is robust, i.e., they increase the steady force level. Figure 6A–C correspond, respectively, to point process intensities equal to FR(t) = 65 + 10 sin(2π20t), FR(t) = 65 + 20 sin(2π20t), and FR(t) = 65 + 30 sin(2π20t), respectively. The simulations had a 15 s duration, the rate oscillations of the MN inputs began at 3 s. The muscle mean force levels increased when the input oscillation amplitude was higher.
The dashed vertical line indicates the beginning of the 20 Hz input oscillation. The amplitude of the 20 Hz oscillation is different for each case. A, Amplitude of 10 spikes/s; FR(t) = 65 + 10 sin(2π20t). B, Amplitude of 20 spikes/s; FR(t) = 65 + 20 sin(2π20t). C, Amplitude of 30 spikes/s; FR(t) = 65 + 30 sin(2π20t). With the increase in oscillation amplitude, there is an increase in the steady-state mean force value.
Figure 7 shows simulation results for different carrier oscillation frequencies of the amplitude-modulated firing rate sinusoidal oscillation at the MN pool input. The frequencies were 20 Hz (beta band), 40 Hz (low gamma band), and 80 Hz (high gamma band), respectively, for Figure 7A–C. The amplitude-modulated oscillation of the MN inputs began after 3 s of simulation and lasted until 15 s. The muscle force oscillated at 1 Hz for all carrier frequencies, evidencing the robustness of the demodulation, albeit the fidelity of the force oscillation decreased for higher carrier frequencies.
The dashed vertical line indicates the beginning of the input amplitude-modulated oscillations. The carrier frequency is different in each case. A, Carrier frequency at 20 Hz; FR(t) = 65 + 20 sin(2π20t) · [1 + 0.5 sin(2πt)]. B, Carrier frequency at 40 Hz; FR(t) = 65 + 20 sin(2π40t) · [1 + 0.5 sin(2πt)]. C, Carrier frequency at 80 Hz; FR(t) = 65 + 20 sin(2π80t) · [1 + 0.5 sin(2πt)]. The 1 Hz oscillation in muscle force appears for all the carrier oscillation frequencies, although with a lower amplitude for higher frequencies.
In the simulations in Figure 8, the 20 Hz oscillation of the firing rate of the MN inputs began after 3 s. Figure 8A shows the muscle force for the descending point processes following a gamma distribution with shape factor 7 instead of Poisson point processes (Watanabe et al., 2013). Figure 8B shows the muscle force with Renshaw cells incorporated to the neuromuscular model (Cisi and Kohn, 2008). Figure 8C shows the muscle force for MNs with active dendrites (Elias et al., 2012; Elias and Kohn, 2013), instead of passive dendrites as used in the previous simulations. The general effect of a steady force increase and maintenance during the beta-band input oscillation was found in all three cases, albeit the force variability was considerably increased when there were Renshaw cells in the feedback. A more detailed analysis of the effects of Renshaw cell feedback or active dendrites on force variability and system behavior could be pursued in future studies using appropriate frequency domain analysis tools.
The dashed vertical line indicates the beginning of the input amplitude-modulated oscillations. A different component was added to the neuromuscular model in each case. The basal firing rate was adjusted for each case so the mean contraction level was 10% MVC. A, Descending commands with a gamma point process with shape factor 7 [basal firing rate, 65 spikes/s; FR(t) = 65 + 20 sin(2π20t); IN ISI, 8 ms]. B, Renshaw cells were added to the model [basal firing rate, 108 spikes/s; FR(t) = 108 + 40 sin(2π20t); IN ISI, 3 ms]. C, Active dendrites were added to the model [basal firing rate, 33 spikes/s; FR(t) = 33 + 20 sin(2π20t); IN ISI, 12 ms].
Figure 9 shows the results from a simulation using an alternative to the AM descending drive to achieve a time-varying force, by adding a 1 Hz sinusoid to a 20 Hz sinusoid. Figure 9A shows a 15 s window from the 60 s of this simulation. The force oscillates at 1 Hz, as in Figure 5A. However, the coherence between the MN synaptic inputs and the CST computed from five randomly chosen MNs, shown in Figure 9B, presents a peak at 1 Hz, in addition to the peak at 20 Hz. Figure 9C presents the mean spectra of the CST of five MNs randomly chosen from 5% of those having the lowest AR (black thick line) and of the CST of five randomly chosen MNs from 5% of those with the highest AR (gray thick line). Differently from Figure 5D, the power of both spectra are approximately equal. The thin lines indicate the individual CST spectra.
Simulation using a CSI formed by a 1 Hz sinusoid added to a 20 Hz sinusoid. A, Force signal during a 15 s window. The force oscillated at 1 Hz. B, Coherence between synaptic conductance and the CST of five randomly chosen MNs of the. Individual data from 100 randomly chosen CSTs are in gray, and the mean is in black. The 95% confidence level is shown as a dashed line. C, Darker line, Mean spectrum of the CST of five randomly chosen MNs from the 5% of the MNs from the pool with the lowest AR; lighter line, mean spectrum of the CST of five randomly chosen MNs from the 5% of the MNs from the pool with the highest AR. Individual data from the 100 randomly chosen CSTs are lighter gray.
Discussion
This study analyzed how beta (or higher frequency) band CSIs to an MN pool can influence the generation of static and oscillatory force. Due to the relatively high number of nonlinear elements (the MNs) and the stochastic nature of the spike trains, it is not easy to predict many behaviors that the population of such elements may have. Computer simulations of a biologically compatible large-scale model can provide quantitative and qualitative predictions and interpretations (mechanisms) that would be almost impossible to reach by intuition alone. The present modeling study provides new insights on how cortical output may influence force control.
The coherence results, relating the inputs (CSI) and outputs (CST) of the MN pool (Figs. 3D,E, 4D) are qualitatively remindful of those observed experimentally (Conway et al., 1995; Schoffelen et al., 2005; Andrykiewicz et al., 2007; Mendez-Balbuena et al., 2012; Perez et al., 2012; Ushiyama et al., 2012), indicating an effective transmission of the high-frequency (20 Hz) CSI to the MN pool (Conway et al., 1995). However, coherence measures are only able to capture linear relations. For example, the absence of corticomuscular coherence at low frequencies does not exclude a nonlinear influence of beta (or higher) band oscillatory CSIs in force control. Indeed, the mean force increase (60–120 s; Fig. 3G) and the appearance of a 1 Hz muscle force oscillation (starting at 60 s; Fig. 5A) correspond to power increases in the force spectrum at 0 and 1 Hz, respectively, indicating a nonlinear behavior, since in both cases no component at these frequencies was present in the CSIs (Figs. 3B, 4B).
The analysis of the firing behavior of the MNs (Figs. 3H, 5B) helps clarify what the main mechanism is behind the extraction of low-frequency information from the beta-band oscillations in the CSIs. It can be noted that the mean force increase (60 to 120 s; Fig. 3G, brown) is associated with the recruitment of additional motor units (having a higher twitch amplitude) in response to the 20 Hz oscillation (Fig. 3H, inset, brown dots). The MN firings, approximately synchronized to the 20 Hz oscillations, generate only a steady increase in force level because of the low-pass filtering by the muscle. A similar behavior can be observed for the 1 Hz amplitude modulation of the 20 Hz oscillation (Fig. 5). When the input signal was amplitude modulated, the MNs with the lowest AR (Fig. 5B, darker red points) were recruited and derecruited at exactly 1 Hz, generating the muscle force oscillation (Fig. 5A,B, vertical lines). This effect can be better analyzed with the CST spectra (Fig. 5D). Whereas the 20 Hz component appears in the spectra of all MNs, as shown by Farmer et al. (1993) and Negro and Farina (2012), the 1 Hz power depends on the MN thresholds. The CST power spectrum at 1 Hz from the low-AR (higher recruitment threshold) MNs (Fig. 5D, red curve) is significantly higher than that corresponding to high-AR MNs (yellow curve). This means that the last recruited MNs are the main contributors to the low-frequency muscle force oscillation, since the force frequency components are a low-pass-filtered version of the CST power spectrum (Negro et al., 2009).
A consequence of the nonlinear behavior of the MN pool is that a similar steady force level (Fig. 3G, blue, green) could be obtained from a lower mean firing rate input (58 spikes/s instead of 65 spikes/s) if a 20 Hz oscillation exists (Fig. 2A, green). As hypothesized by Baker et al. (1997) and Salenius et al. (1997), the present study shows that a 20 Hz CSI to an MN pool can produce an isometric force with minimum pyramidal tract discharge. This means that a lower metabolic cost can be achieved with the presence of beta oscillations in the CSI.
The nonlinear behavior of the MN pool contradicts the conclusions from the studies by Lowery et al. (2007), Stegeman et al. (2010), and Farina et al. (2014), which analyzed an MN pool operating at ∼5% MVC with most MNs recruited. These studies dismissed any role for the high-frequency oscillations, since the muscles cannot respond to frequencies in the MN CSTs higher than 6–10 Hz (Negro et al., 2009). The results presented here lead to the conclusion that the extraction of low-frequency information from higher-frequency oscillations is possible mainly due to the different recruitment threshold of each MN. Experimental data indicate that, in most of the situations, the MN pool is not fully recruited. For example, the soleus muscle has its MNs fully recruited at 100% MVC (Oya et al., 2009), the adductor pollicis at 50% MVC (Kukulka and Clamann, 1981), the first dorsal interosseus between 50 and 75% MVC (De Luca et al., 1982; Moritz et al., 2005), the biceps brachii at 88% MVC (Kukulka and Clamann, 1981), and the abductor hallucis at 58% MVC (Kelly et al., 2013).
The results presented in this study were from the neuromuscular model associated with the soleus muscle (Watanabe et al., 2013). However, the results are very robust, and are similar for differences in number and type of MNs, parameter value range, MN dendritic model, MN pool structure (e.g., with inhibitory feedback from Renshaw cells), IN ISI, input point process statistics, and high-frequency oscillatory frequency. Actually, we conjecture that the behavior described here should occur in any neuronal network having a gradation of firing thresholds.
A caveat is that the low-frequency amplitude modulation of the 20 Hz (or higher) oscillation is just one possibility to convey the commands for time-varying force production. Whether it is really used by the motor cortex remains to be verified. An alternative (among others) to the amplitude modulation would be low-frequency modulations of the CSI intensity, as used in the studies by Lowery et al. (2007), Stegeman et al. (2010), and Farina et al. (2014). We simulated this alternative using a 1 Hz sinusoid for the low-frequency modulation coexisting with a 20 Hz oscillation, the result being that the beta-band oscillation provided a steady force around which the 1 Hz force developed (Fig. 9A). Differently from the AM case (Fig. 4D), the coherence between the CSI and the MN CSTs had an additional peak at 1 Hz (Fig. 9B), and the spectra of the MN CSTs had a similar power at 1 Hz for both low- and high-AR MNs (compare with Figs. 5D, 9C). Although most of the existing experimental data on corticomuscular coherence for subjects producing an oscillatory force (Andrykiewicz et al., 2007; Chakarov et al., 2009; Naranjo et al., 2010; Mendez-Balbuena et al., 2012; Bayraktaroglu et al., 2013; Jacobs et al., 2015) do not show significant values at the force oscillation frequency, it could be that the EEG is not able to mirror, with a reasonable signal-to-noise ratio, what is actually coming out of the motor cortex at this frequency range. Even though the spiking activity of the pyramidal tract neurons contributes linearly to the spectra of local field potentials (Baker et al., 2003; Manning et al., 2009) and the EEG (Ball et al., 2008; Nelson and Pouget, 2012; Heitmann et al., 2013), these extracellular field potentials are influenced by many different sources (Buzsáki et al., 2012). If the idealized situation simulated here were changed to one with many sources of variability (e.g., phase noise in oscillations, additive and multiplicative random sources), the resultant coherence peak values would be more similar to those obtained experimentally. For this reason, comparisons of the coherence functions obtained in the simulations with those obtained experimentally (based on EEG, local field potential, or MEG) should be done only qualitatively (Williams and Baker, 2009a,b).
It is desirable that predictions of theoretical or computer-simulated results be experimentally testable. The demodulation of the CSI by the MNs can be tested in at least two ways during the generation of an oscillatory force. The first would be to compare the CST spectra of a set of the first recruited MNs with a set of the last recruited MNs. If the CSI is equivalent to an amplitude-modulated 20 Hz, the spectra should be similar to Figure 5D and not Figure 9C. The second manner would be to verify whether the firing rates of the cortical output neurons are being amplitude modulated (Eq. 2). Both suggestions require sophisticated recording techniques of the activities of large numbers of neurons simultaneously. Particularly for the MNs, it would be necessary to differentiate the set of the first recruited from the set of the last recruited.
In summary, our results suggest a putative role, at least for certain time windows, for high-frequency cortical output: it contributes to regulating forces minimizing energy expenditure by the motor cortex. The mechanism behind the capability of an MN pool to decode high-frequency cortical oscillations is the gradation of recruitment thresholds found in MN pools of humans and other species.
Footnotes
This work was funded by the Brazilian National Council for Scientific and Technological Development (Grant 303313/2011) and the São Paulo Research Foundation (Grant 2011/21103-7, R.N.W.). We thank Dr. Leonardo Abdala Elias for helpful discussions.
The authors declare no competing financial interests.
- Correspondence should be addressed to either Renato N. Watanabe or Andre F. Kohn, Laboratorio de Engenharia Biomedica, Escola Politécnica, Universidade de São Paulo, CEP 05508-010, São Paulo, Brazil. renato.watanabe{at}usp.br or andfkohn{at}leb.usp.br
