Effective Stimulation Type and Waveform for Force Control of the Motor Unit System: Implications for Intraspinal Microstimulation

The input–output properties of spinal motoneurons and muscle fibers comprising motor units are highly non-linear. The goal of this study was to investigate the stimulation type (continuous versus discrete) and waveform (linear versus non-linear) controlling force production at the motor unit level under intraspinal microstimulation. We constructed a physiological model of the motor unit with computer software enabling virtual experiments on single motor units under a wide range of input conditions, including intracellular and synaptic stimulation of the motoneuron and variation in the muscle length under neuromodulatory inputs originating from the brainstem. Continuous current intensity and impulse current frequency waveforms were inversely estimated such that the motor unit could linearly develop and relax the muscle force within a broad range of contraction speeds and levels during isometric contraction at various muscle lengths. Under both continuous and discrete stimulation, the stimulation waveform non-linearity increased with increasing speed and level of force production and with decreasing muscle length. Only discrete stimulation could control force relaxation at all muscle lengths. In contrast, continuous stimulation could not control force relaxation at high contraction levels in shorter-than-optimal muscles due to persistent inward current saturation on the motoneuron dendrites. These results indicate that non-linear adjustment of the stimulation waveform is more effective in regard to varying the force profile and muscle length and that the discrete stimulation protocol is a more robust approach for designing stimulation patterns aimed at neural interfaces for precise movement control under pathological conditions.

Introduction

Recent advances in neural interface technology have allowed the direct modulation of nervous system functions by injecting currents into specific compartments of individual neurons (Holsheimer, 1998; Radivojevic et al., 2016). In regard to neuromuscular systems, individual motor units, consisting of a single motoneuron and its innervating muscle fibers (i.e., the muscle unit), could represent a target modulated at the spinal cord to accurately evoke proper movements (Wagner et al., 2018). However, the prediction of stimulation patterns and their effects on force control have remained difficult, mainly due to the non-linearities inherited in the fundamental elements comprising neuromuscular systems (Heckman and Enoka, 2012).

The input–output relationship of spinal motoneurons is highly non-linear in various species, including rats, mice, cats, and turtles (Heckman et al., 2003). The firing mode of the motoneuron may transition from quiescent to regular firing or from low- to high-frequency firing in response to brief excitatory current injection at the soma, thereby revealing the occurrence of bistability (Hounsgaard et al., 1988). While slowly increasing and decreasing current injection at the soma, low-threshold motoneurons (presumably slow-type motoneurons) tend to exhibit notable counterclockwise hysteresis and a self-sustaining firing behavior below the firing initiation threshold in the descending stimulation phase, whereas high-threshold motoneurons (presumably fast-type motoneurons) tend to exhibit slight clockwise hysteresis under the absence of self-sustaining firing behavior (Lee and Heckman, 1998b). The underlying mechanism of this non-linear input–output relationship has been suggested as the spatiotemporal interaction between action potential-producing membrane mechanisms at the soma and plateau potential-generating calcium channels in dendritic areas (i.e., 300–800 μm from the soma) (Kim et al., 2014). These dendritic calcium channels (presumably L-type Ca_v1.3 channels) are actively involved through monoaminergic neuromodulation due to the brainstem regarding normal motor behavior and through endogenous monoamines in regard to spinal cord injury (Heckmann et al., 2005).

The input–output relationship of muscle fibers has also been demonstrated to be highly non-linear in frogs, rats, and cats (Mrowczynski et al., 2006). Compared to its low frequency, a greater muscle force is produced in response to high-frequency current stimulation, resulting in a sigmoid curve of the relationship between the stimulation frequency and force output. In addition, the muscle force is maximized at the optimal muscle length over smaller and larger muscle lengths, thereby revealing a bell-shaped relationship between the muscle length and force output during isometric contraction (Winters et al., 2011). This non-linear relationship is attributable to the complex interactions between calcium dynamics, cross-bridge formation, and length variation in the sarcoplasm of muscle fibers (Kim et al., 2015). In summary, the dynamics of muscle activation are non-linearly related to both the stimulation frequency and muscle length and greatly decrease under shorter-than-optimal lengths at the physiological stimulation frequency (i.e., <20 Hz) (Perreault et al., 2003).

Intraspinal microstimulation exhibits the capability of activating a specific set of neurons within the spinal cord to modulate muscle activity in different body parts, including the hindlimb (Mushahwar and Horch, 1998), forelimb (Sunshine et al., 2013), and respiratory system (Sunshine et al., 2018). Recently, the relationship between current stimulation and motoneuron firing has been investigated via the direct application of direct or pulsed current injection to the spinal cord. The intradural region of the spinal cord in mice has been targeted for direct current stimulation purposes to elucidate the influence of the stimulation polarity on the firing outputs of hindlimb motoneurons (Ahmed, 2016). Computational and imaging studies have further demonstrated that the asymmetric waveform of a biphasic stimulus current pulse may enhance the target selectivity, and the stimulus frequency may alter the neuronal output (McIntyre and Grill, 2000; Wang et al., 2012). However, little is known regarding the stimulation pattern to produce desired force profiles via skeletal muscles at the motor unit level, which is the smallest element underlying all movements.

Here, we theoretically investigated the effective waveforms during intraspinal microstimulation to control the muscle force under two types of stimulation protocols, namely, continuous and discrete (or impulse) current stimulation protocols, during isometric muscle contraction within a full physiological range of the muscle length and force level. The present study focused on the direct activation of a spinal motoneuron with a microelectrode in close proximity to its initial segment and cell body. Extracellular microstimulation of other components, such as axons or dendrites, was not involved in the current analysis. We hypothesize (1) that non-linear stimulation waveforms are needed to elicit linear force generation at the motor unit level and (2) that the discrete current stimulation protocol more effectively prevents the non-linearities induced by the active channels in motoneuron dendrites, including muscle afferent inputs due to muscle length variation. Simulation analysis demonstrated the systematic stimulation waveform non-linearity depending on the stimulation type, muscle length, and force profile. These simulation results provide a basis for the design of stimulation patterns for neural interfaces to enhance motor precision control under pathological conditions.

Materials and Methods

Motor Unit Model

Stimulation waveforms were investigated in regard to extracellular microstimulation with a sharp-tipped electrode near the initial segment and cell body of a spinal motoneuron (Figure 1). A motor unit model comprising a two-compartment motoneuron model and a three-module muscle unit model was physiologically constructed and simulated in Python-based Motor Unit Simulator (PyMUS) software, which was developed for virtual experiments on single motor units under a wide range of physiological conditions (Kim and Kim, 2018).

FIGURE 1

Figure 1. Simulation setup. (A) Reduction in the anatomical motoneuron (top, dendrogram) into a two-compartment framework (bottom) as a function of the path length (D_path) from the center of the cell body (D_path = 0). The left arrow indicates the mapping of the soma/axonal hillock/initial segment of the anatomical model into the somatic compartment of the reduced model. The right arrow indicates the mapping of all points over the narrow band (a mean value of 0.6 mm) of the anatomical dendrites into the dendritic compartment of the reduced model. C_m, G_m, G_a G_C, V_i, and V_e indicate the membrane capacitance, leak conductance, voltage-gated conductance, coupling conductance, and intracellular and extracellular potentials, respectively. The subscripts S and D indicate the somatic and dendritic compartments, respectively. I_ext, the current passing through the microelectrode, is located close to the initial segment and cell body. In the current analysis, V_e,S is determined by I_ext, while the influence of I_ext on V_e,D is assumed to be negligible. (B) Model motor unit comprising the two-compartment motoneuron model and three-module muscle unit along with the muscle spindle afferent (I_a). I_eq is the equivalent current intracellularly injected at the somatic compartment of the motoneuron to evoke the same transmembrane potential induced by I_ext. [CaT], A, F_T, and X_m indicate the concentration of calcium bound to troponin, the muscle activation level, the force output, and the muscle length, respectively.

Briefly, in terms of the reduced motoneuron model, one compartment (referred to as the somatic compartment) represents the soma and axonal hillock/initial segment, and the other compartment (denoted as the dendritic compartment) represents the dendritic regions, including persistent inward current-generating voltage-gated channels (Booth et al., 1997). These two compartments are connected via the coupling conductance, representing the electrical distance between the soma and dendritic regions. Five passive parameters, including the specific membrane conductance and capacitance (i.e., G_m,S and C_m,S of the somatic compartment, G_m,D and C_m,D of the dendritic compartment, and G_C in regard to the coupling conductance between these two compartments), were analytically determined to capture five electrotonical properties, including the somatic input resistance (i.e., R_N,S), system time constant (i.e., τ_m), and three voltage attenuation factors between the soma and all dendritic sites separated by a similar path length (i.e., D_path) (i.e., ${\text{VA}}_{\text{SD}}^{\text{DC}}$ and ${\text{VA}}_{\text{SD}}^{\text{AC}}$ describing soma-to-dendrite propagation of direct and alternative currents and ${\text{VA}}_{\text{DS}}^{\text{DC}}$ describing dendrite-to-soma propagation of a direct current). Action and plateau potentials were generated in the somatic and dendritic compartments via the incorporation of Hodgkin–Huxley-style active channels [i.e., G_Naf, G_Nap, G_Kdr, G_K(Ca), and G_Can, with dynamic changes in the calcium reversal potential of the soma, and G_Cal, with a constant calcium reversal potential of the dendrites]. The passive parameter values were first determined to capture the passive properties [R_N,S = 1.29 MΩ, τ_m = 7.2 ms, ${\text{VA}}_{\text{SD}}^{\text{DC}}$(D_path = 0.6 mm) = 0.76, ${\text{VA}}_{\text{SD}}^{\text{AC}}$(D_path = 0.6 mm) = 0.27, and ${\text{VA}}_{\text{DS}}^{\text{DC}}$(D_path = 0.6 mm) = 0.75] measured at the soma and dendrites as previously reported in Zengel et al. (1985); Spruston et al. (1994), and Kim et al. (2014). Then, the active parameters were determined to replicate the active properties (spike height = 92.3 mV, rheobase current = 10.5 nA, afterhyperpolarization duration = 98.5 ms, and depth = 3.1 mV under passive dendrites, and effective persistent inward current = −22 nA under active dendrites) obtained from the cell body as previously mentioned in Zengel et al. (1985); Hochman and McCrea (1994), and Hornby et al. (2002). The voltage–current (i.e., VI) and current–voltage (i.e., IV) properties of the model motoneuron were validated via a comparison to those of an anatomically reconstructed motoneuron model that reproduces non-linear behaviors observed from in vivo cat α-motoneurons with low thresholds during voltage- and current-clamping at the cell body (Kim et al., 2014) (please refer to Figure 2A for non-linear behavior of the model motoneuron).

FIGURE 2

Figure 2. Non-linear behavior of the motoneuron (A) and motor unit (B). (A) Current stimulation applied to the somatic compartment of the motoneuron model (top), instantaneous firing rate (middle), and dynamics of the PIC activation level in the motoneuron dendrite (bottom). The arrows indicate the responses at the three muscle lengths (minimal length with X_m = −16 mm, optimal length with X_m = −8 mm, and maximal length with X_m = 0 mm). (B) Frequency of impulse current stimulation applied to the muscle unit model (top), force normalized with the peak force at the optimal length (middle), and dynamics of the muscle activation level (bottom). The arrows indicate the responses at the three muscle lengths (minimal length with X_m = −16 mm, optimal length with X_m = −8 mm, and maximal length with X_m = 0 mm).

Regarding the modular muscle unit model, the muscle force was produced through three physiological procedures in response to action potentials originating from the motoneuron. Module 1 transforms neural signals into dynamics of the calcium concentration (i.e., Ca_SP) in the sarcoplasm, including calcium release and uptake of the sarcoplasmic reticulum, calcium-buffering proteins, and calcium-bound troponin complex. Module 2 then transforms the concentration dynamics of calcium-bound troponin (i.e., Ca_SPT) into the degree of cross-bridge formation, representing muscle activation dynamics (i.e., A). Finally, module 3 transforms the above muscle activation dynamics into force based on Hill-type muscle–tendon mechanics, reflecting length–and velocity–tension relationships. Model parameter values were determined for each module to reproduce the force generation process. These values reflect sarcoplasmic calcium dynamics (Westerblad and Allen, 1994) in module 1, calcium–force relationship (Shames et al., 1996) in module 2, length–and velocity–tension properties (Scott et al., 1996) in module 3, and muscle activation dependence on muscle length variation (Sandercock and Heckman, 1997) in modules 1 and 3. The input–output properties of the model muscle unit have been validated through a comparison to those of adult cat soleus muscles within a full physiological range of stimulation rates (i.e., 1–100 Hz) and muscle lengths (i.e., approximately -16 to 0 mm) during isometric, isokinetic, and dynamic contractions (Kim et al., 2015) (please refer to Figure 2B for non-linear behavior of the model muscle unit).

The axonal nerve coupling the motoneuron with the muscle unit was modeled via a single parameter representing the delay time (∼10 ms) required for the transmission of action potentials from the motoneuron to the muscle unit, assuming perfect action potential transmission from the motoneuron to the muscle unit (Kim, 2017). The overall input–output properties of the model motor unit were validated by comparison to those of a physiologically reconstructed model for an adult cat slow-type motor unit in the previous study [please refer to Figures 8, 9 in Kim and Kim (2018) for non-linear behavior of the model motor unit]. The system equations and parameter values used for this study are presented as supplementary material (Supplementary File 1).

Stimulation Protocols

The non-linearity of the input–output relationship was investigated separately for the motoneuron and muscle unit. In regard to the motoneuron, a linearly ascending and descending current with a peak value of 20 nA over 10 s was injected at the soma of the motoneuron (as shown in Figure 2A). In regard to the muscle unit, a train of current impulses was applied such that the stimulation frequency increased and decreased according to a triangular shape with a peak value of 40 Hz over 10 s (see Figure 2B). The peak stimulation frequency was selected to ensure that the model muscle unit produced the maximal force during isometric contraction at the optimal muscle length.

In the present study, stimulation waveforms were investigated under extracellular motoneuron microstimulation with a sharp-tipped electrode near the initial segment and cell body (Figure 1). Under these conditions, the extracellular potential (i.e., V_e,S) at the soma is determined by the current (i.e., I_ext) passing through the microelectrode, whereas the extracellular potential (i.e., V_e,D) across the dendritic area (i.e., D_path > 0.6 mm) distal to the soma tends to not depend on I_ext (McIntyre and Grill, 2002). The variation in transmembrane potential (i.e., V_S) at the soma due to I_ext was simulated via intracellular injection of the equivalent current (i.e., I_eq) at the soma, as proposed in a previous study (Warman et al., 1992). With the application of Kirchhoff’s current law to the somatic compartment of the two-compartment model in addition to approximation of V_e,D as zero under passive membrane conditions, I_eq can be derived as the following analytical form:

$$\begin{array}{c}\hfill I_{\text{eq}}=C_{m,\text{S}}\cdot \text{d}{V}_{\text{S}}\mathrm{/}\text{d}t+G_{m,\text{S}}\cdot V_{\text{S}}+G_{\text{C}}\cdot V_{\text{S}}=-G_{\text{C}}\cdot V_{\mathrm{e},\mathrm{S}}\hfill \end{array}$$

Where V_S is the potential difference between the intracellular (V_i,S) and extracellular sides (V_e,S) of the somatic compartment in the two-compartment neuron model.

Linear force production of the motor unit was induced considering two types (continuous and discrete) of current stimulation (i.e., I_eq) in the somatic compartment of the motoneuron model over 10 s. The shape and amplitude of continuous current stimulation were first determined with a linear function in a piecewise manner so that the force linearly increased and decreased at the various speeds over 10 s. Then, piecewise linear functions were fitted with a continuous function ensuring linear force production by the motor unit (please refer to Table 1). In the case of discrete current stimulation, the timing and amplitude of the current pulses were adjusted such that the discrete current simulation process reproduced the temporal evolution of force production induced under continuous current stimulation condition. In this study, the amplitude of the current pulses indicated the minimal current intensity yielding a desired force profile. The width of the current pulse was set to 0.5 ms, ensuring a one-to-one stimulus-firing reaction at the stimulus amplitude threshold. To compare the different states between the various muscle lengths, all muscle forces induced under current stimulation in the somatic compartment of the motoneuron model were normalized based on the maximum force (i.e., P₀) actively produced by the model motor unit during isometric contraction at the optimal muscle length. In the present study, the simulation data are not shown when the model motor unit could not produce the desired peak force (i.e., 20–100% of P₀).

TABLE 1

Table 1. Piecewise equations for the time courses of continuous current stimulation (i.e., I_eq) presented in Figures 3, 5, 7.

Muscle Length Variation

The stimulation patterns required for triangular force production by the motor unit during isometric contraction were predicted under three different modes of the muscle length: physiologically minimal (X_m = −16 mm), optimal (X_m = −8 mm), and maximal (X_m = 0 mm) muscle lengths. The physiological range of the variation in muscle length was determined based on the locomotor data of adult cats reported in a previous study (Goslow et al., 1973). The variation in muscle length influences both the motoneuron and muscle unit by transmitting muscle spindle signals to the motoneuron and adjusting a certain kinetic parameter [i.e., K5 in Kim et al. (2015)] of the formation of calcium-binding troponin in the sarcoplasm of the muscle unit model as a function of the muscle length.

Muscle Afferent Signals

The feedback signals originating from the muscle spindle were simulated by summing the excitatory synaptic inputs (i.e., G_esyn) over both the soma and dendrite (Segev et al., 1990). Under passive membrane conditions, the peak conductance for muscle afferent inputs was determined to match the current (i.e., I_N,PASS) experimentally measured during voltage clamping over the soma under a 10-mV hyperpolarized potential below the resting potential while slowly varying the muscle length from the physiologically minimum length to the maximum length (Lee et al., 2003). The peak conductance for muscle afferent inputs was set to 0 mS/cm² to yield an I_N,PASS value of 0 nA at −16 mm, 0.0064 mS/cm² to yield an I_N,PASS value of 2.5 nA at −8 mm, and 0.0128 mS/cm² to yield an I_N,PASS value of 5 nA at 0 mm under isometric conditions.

Simulations

All simulations were performed in the PyMUS software environment (version 2.0.1) via the Python integration method (a Variable-coefficient Ordinary Differential Equation (VODE) was applied to the motoneuron, and the Livermore Solver for Ordinary Differential Equations (LSODE) was applied to the muscle unit) at a fixed time step (0.1 ms) on a desktop computer operated by 64-bit Microsoft Windows 10. Default values of the model parameters and simulation conditions were implemented unless noted otherwise in the text. The computer codes of the models and simulations in this study are presented in the supplementary information (Supplementary File 2) and are publicly available from the public repository of GitHub¹. The data employed for the estimation of the stimulation waveforms required for force production by the motor unit at the various speeds and levels in the present study are presented in the supplementary material (Supplementary File 3).

Results

Non-linear Input–Output Relationship of the Motoneuron and Muscle Unit

We first evaluated whether the model motor unit captures the non-linearity of the input–output function of the motoneuron and muscle unit. Figure 2 shows the non-linear response of the motoneuron and muscle unit under linearly increasing and decreasing excitation effects. The model motoneuron suitably replicated experimental observations of non-linear firing of spinal motoneurons during triangular current injection at the cell body [please refer to Figure 3 in Hounsgaard et al. (1984)]. In particular, the firing rate was enhanced at a current intensity higher than the recruitment threshold in the ascending stimulation phase, and firing was sustained below the recruitment threshold in the descending stimulation phase (as shown in the top and middle panels of Figure 2A). This result was mainly attributed to the activation of persistent inward current (PIC)-generating calcium channels at the dendrite underlying the non-linear input–output relationship of the motoneuron (as shown in the bottom panel of Figure 2A). In terms of the model muscle unit, a non-linear relationship of the stimulation frequency and force genesis has been experimentally verified during electrical stimulation of soleus muscles in adult cat preparations [please refer to Figure 8 in Rack and Westbury (1969)]. The muscle unit model reproduced the sigmoidal shape of force production well when slowly increasing and decreasing the stimulation rate (as shown in the top and middle panels, respectively, of Figure 2B). This was mainly attributed to the non-linear dynamics of muscle activation (as shown in the bottom panel of Figure 2B). Furthermore, the non-linear input–output relationship of the motoneuron and muscle unit was modulated by the variation in muscle length, indicating the influence of the muscle afferent inputs and the length dependence of muscle activation. These results indicate that the stimulation waveforms required for the motoneuron to control the muscle force may strongly depend on the muscle length.

Continuous Stimulation Waveforms for Force Control at the Various Muscle Lengths

Figure 3 shows the waveforms under continuous current stimulation in regard to linear force generation at the various speeds and levels during isometric contraction at the minimal, optimal, and maximal muscle lengths. Overall, the current intensity exponentially increased and decreased in regard to linear force development and relaxation, respectively, at all muscle lengths. The current stimulation amplitude increased with increasing force speed and magnitude and with decreasing muscle length. The length dependence of the stimulation amplitude was attributed to the reduction in muscle afferent inputs and muscle activation level with decreasing muscle length from the physiologically maximal length. In accordance, the peak firing rate of the motoneuron was the highest at the highest speed and largest magnitude of force production in the shortened muscle state.

FIGURE 3

Figure 3. Continuous stimulation waveforms for linear force production. Intensities of continuous current stimulation (I_eq) applied to the somatic compartment (top), instantaneous firing rate (upper middle) in the motoneuron, and normalized muscle force (lower middle) produced by the muscle unit at the three muscle lengths (bottom): physiologically minimum (A), optimal (B), and maximum muscle lengths (C). The level of force production produced by the motor unit is denoted as a percentage of the maximal isometric force (P₀) at the optimal muscle length and indicated with different colors.

In terms of the shortened muscle, the motor unit produced various force levels ranging from 20 to 80% of the maximal force (i.e., P₀) generated at the optimal muscle length (Figure 3A). The 100% P₀ level was not achievable due to the bell-shaped length-tension muscle properties, indicating a force decline with decreasing or increasing muscle length from the optimal length. The stimulation waveforms were almost symmetric up to 60% of P₀. At 80% of P₀, however, the stimulation intensity rapidly decreased near the recruitment threshold to achieve linear force relaxation in the descending stimulation phase. This was required to prevent the full activation of the PIC channels located across the motoneuron dendrite in the descending stimulation phase. At the optimal muscle length, the current threshold for firing initiation was lower than that in the shortened muscle case due to the excitatory muscle afferent inputs to the motoneuron (Figure 3B). Brief injection of an excitatory step current at the beginning of current stimulation was required to slowly produce force up to 20% of P₀. This result indicates that the PIC may partially be activated and slowly oscillate, leading to slow firing at the motoneuron. Notably, the model motor unit could not fully reach the maximal force level due to the limitation of the stimulation intensity to prevent full PIC activation in the motoneuron. Once the PIC was fully activated in the ascending stimulation phase, the deactivation of PIC channels became uncontrollable, resulting in difficulty in achieving linear force relaxation in the descending stimulation phase. This result could be attributed to the severe attenuation of the electrical signals generated near the cell body when transmitted to distal dendritic areas in the motoneuron. In terms of the lengthened muscle, the current threshold for firing initiation further decreased due to the increase in excitatory muscle afferent inputs to the motoneuron (Figure 3C). Force production was realized from the peak level of 40% of P₀ to 80% of P₀ owing to the twitch force greater than 20% of P₀ and the bell-shaped length-tension muscle properties.

All these results suggest that the exponential waveform under a continuous current intensity effectively controlled linear force production and that the continuous type of current stimulation could limit the range of force generation in shorter-than-optimal muscles.

Discrete Stimulation Waveforms for Force Control at the Various Muscle Lengths

After stimulation waveform estimation under continuous stimulation conditions, the stimulation waveforms required for force production were assessed under discrete stimulation conditions (Figure 4). In general, the rate of current impulse stimulation should exponentially increase and then decrease to produce forces according to a triangular shape at all muscle lengths. Similar to the continuous stimulation case, the amplitude and peak rate of the current impulses increased with increasing speed and magnitude of force production and with decreasing muscle length. The length dependence of the stimulation amplitude and rate was attributed to the increase in muscle afferent inputs and muscle activation level with increasing muscle length from the shortened state. In accordance, the peak firing rate of the motoneuron was the highest at the minimal muscle length and the highest speed and level of force production. However, in contrast to the continuous stimulation case, the stimulation waveforms were almost symmetric at all force production levels during isometric contraction at the minimal muscle length (see Figure 4A). Furthermore, the model motor unit fully produced the maximum force at the optimal muscle length without full activation of the PIC channels across the motoneuron (Supplementary File 4).

FIGURE 4

Figure 4. Discrete stimulation waveforms for linear force production. Frequency and amplitude (two top panels) of impulse current stimulation (I_eq) applied to the somatic compartment, instantaneous firing rate (upper middle) in the motoneuron, and normalized muscle force (lower middle) produced by the muscle unit at the three muscle lengths (bottom): physiologically minimum (A), optimal (B), and maximum muscle lengths (C). The level of force production produced by the motor unit is denoted as a percentage of the maximal isometric force (P₀) at the optimal muscle length and indicated with different colors.

All these results indicate that the exponential waveform due to the current impulse frequency is effective in regard to linear production of the muscle force and that the discrete type of current stimulation is suitable for the full range of force generation and muscle length levels.

Influence of the PICs Over the Motoneuron Dendrites on the Stimulation Waveforms

We further investigated the influence of the PIC channels responsible for the plateau potentials over the motoneuron dendrites on the stimulation waveforms required for force control. To this end, we compared the stimulation waveforms obtained with and without PIC channels over the dendritic compartment of the motoneuron model, as shown in Figures 3, 4, respectively. Figure 5 shows the stimulation waveforms under continuous stimulation conditions without PIC channels over the motoneuron. Overall, the current stimulation amplitude increased to compensate for the lack of an intrinsic PIC current in the motoneuron when the motoneuron PICs were excluded at all muscle lengths. In contrast to the case involving PIC channels over the motoneuron dendrite (please refer to the dashed-line boxes in Figure 5), the stimulation waveforms were symmetric between the ascending and descending force production phases at all muscle lengths (Figure 5A). Notably, the model motor unit without PIC channels over the motoneuron fully reached the maximal force level at the optimal muscle length during linear force development (Figure 5B). However, the force level of 20% of P₀ was unattainable because of the increase in the recruitment threshold and firing rate in the absence of motoneuron PIC channels (Figure 5B).

FIGURE 5

Figure 5. Influence of the PIC on the stimulation waveforms under the continuous conditions. Intensities of continuous current stimulation (I_eq) applied to the somatic compartment (top), instantaneous firing rate (upper middle) in the motoneuron, and normalized muscle force (lower middle) generated by the muscle unit at the three muscle lengths (bottom): physiologically minimum (A), optimal (B), and maximum muscle lengths (C). The level of force production produced by the motor unit is denoted as a percentage of the maximal isometric force (P₀) at the optimal muscle length and indicated with different colors. The dashed-line rectangles indicate the simulation results with the motoneuron PIC channels, as shown in Figure 3, for the purpose of comparison.

The comparison analysis results revealed the causal relationship of PIC saturation over the motoneuron dendrites and the current intensity limitation and stimulation waveform disruption observed under continuous stimulation conditions.

Figure 6 shows the stimulation waveforms under discrete stimulation conditions without PIC channels over the motoneuron dendrites. Overall, the stimulation waveforms were comparable to those in the case with PIC channels over the motoneuron dendrites. The stimulation rate first exponentially increased and then decreased during linear force development and relaxation within the full range (i.e., from 20 to 100% of P₀). However, the amplitude of the current impulses increased to compensate for the absence of an intrinsic PIC current in the motoneuron. The amplitude increment increased with increasing force level at all muscle lengths (please refer to the dashed-line boxes in Figure 6 representing the PIC activation case). It should be noted that under discrete current stimulation, the model motor unit produced 20% of P₀ (Supplementary File 5), in addition to the maximal force at the optimal muscle length in the absence of PIC channels over the motoneuron dendrites.

FIGURE 6

Figure 6. Influence of the PIC on the stimulation waveforms under discrete conditions. Frequency and amplitude (two top panels) of impulse current stimulation (I_eq) applied to the somatic compartment, instantaneous firing rate (upper middle) in the motoneuron, and normalized muscle force (lower middle) produced by the muscle unit at the three muscle lengths (bottom): physiologically minimum (A), optimal (B), and maximum muscle lengths (C). The level of force production by the motor unit is denoted as a percentage of the maximal isometric force (P₀) at the optimal muscle length and indicated with different colors. The dashed-line rectangles indicate the simulation results with the motoneuron PIC channels, as shown in Figure 4, for the purpose of comparison.

These results reinforce the robustness of discrete stimulation under short-width current pulses for muscle force control purposes regardless of the presence of PICs over the motoneuron dendrites.

Influence of the Muscle Afferent Inputs on the Stimulation Waveforms

We also assessed the influence of muscle afferent feedback on the stimulation waveforms for force speed and magnitude control purposes. Figure 7 shows the stimulation waveforms estimated with and without muscle afferent inputs over the motoneuron model, as shown in Figures 3, 4, respectively. Under both continuous and discrete stimulation conditions, the current stimulation amplitude increased to compensate for the lack of muscle afferent inputs to the motoneuron at the optimal and lengthened muscle lengths. In addition, the peak current intensity under continuous current stimulation conditions was relatively high to prevent the full activation of the PIC channels over the motoneuron in the ascending force production phase. This result could be explained by the absence of muscle afferent-mediated facilitation of PIC activation over the motoneuron dendrites (Figure 2A). However, the model motor unit without muscle afferent feedback did not fully reach the maximal force level at the optimal muscle length under continuous motoneuron current stimulation conditions (Figure 7A). Importantly, in the case of discrete current stimulation, the stimulus rate increased until the maximal force level was reached without the full activation of motoneuron PIC channels (Supplementary File 6).

FIGURE 7

Figure 7. Influence of muscle spindle feedback on the stimulation waveforms. (A,B) Intensities of continuous current stimulation (I_eq) applied to the somatic compartment (top), instantaneous firing rate (upper middle) in the motoneuron, and normalized muscle force (lower middle) produced by the muscle unit at the two muscle lengths (bottom): physiologically optimal and maximum muscle lengths. (C,D) Frequency and amplitude (two top panels) of impulse current stimulation (I_eq) applied to the somatic compartment, instantaneous firing rate (upper middle) in the motoneuron, and normalized muscle force (lower middle) generated by the muscle unit at the two muscle lengths (bottom): physiologically optimal and maximum muscle lengths. The level of force production produced by the motor unit is denoted as a percentage of the maximal isometric force (P₀) at the optimal muscle length and indicated with different colors. The dashed-line rectangles indicate the simulation results with the motoneuron PIC channels, as shown in Figures 3, 4 for the purpose of comparison.

The comparison analysis results indicate that muscle afferent feedback may exacerbate the limitation of the current intensity to prevent full PIC activation during continuous stimulation, and this limitation may also be avoided under discrete stimulation conditions with short-width current pulses.

Discussion

We estimated the waveforms under extracellular microstimulation near the initial segment and cell body of a spinal motoneuron that control the speed and level of force production by the motor unit during isometric contraction at the various muscle lengths. In terms of both the continuous and discrete stimulation types, the non-linearity of the stimulation waveform systematically increased with increasing force speed and level and with decreasing muscle length. Furthermore, the presence of neuromodulatory inputs from the brainstem and afferent feedback originating from the muscle spindle reduced the current stimulation intensity for linear force production. These results may provide a template for the design of stimulation waveforms under varying force profile, muscle length, and neuromodulation conditions.

Continuous Versus Discrete Stimulation Protocols

The non-linearity of the stimulation waveform was assessed by evaluating the difference between the areas under the linear and stimulation curves, i.e., between the lowest and highest points of the stimulation curve (Emancipator and Kroll, 1993). In the present study, the variation in the non-linearity of the stimulation waveform was approximated as the variation in the highest point of the stimulation curve because the lowest point of the stimulation curve remained similar at all produced force levels during isometric contraction at a specific muscle length. Under both the continuous and discrete stimulation protocols, the non-linearity of the current stimulation waveform substantially increased during linear force generation, particularly higher than 60% of the maximal isometric force (i.e., P₀) at shorter-than-optimal muscle lengths (Figure 8). In addition, the current amplitude increased in the absence of motoneuron PIC channels (please refer to Figures 5, 6) and muscle afferent inputs (Figure 7), indicating the influence of neuromodulation from the brainstem and sensory feedback originating from the muscle on the current stimulation pattern.

FIGURE 8

Figure 8. Relationship between current stimulation and force genesis at the various muscle lengths. (A) Peak current amplitude eliciting the peak force during triangular force production under continuous stimulation conditions. (B) Peak impulse frequency eliciting the peak force during triangular force production under discrete stimulation conditions. The peak force is normalized based on the maximal force isometrically produced at the optimal muscle length. The three symbols (circles, squares, and triangles) indicate the physiological minimal (X_m = −16 mm), optimal (X_m = −8 mm), and maximal (X_m = 0 mm) muscle lengths, respectively.

The non-linear current stimulation waveforms at the soma of the motoneuron were captured with an exponential function within a broad range of speeds and force production levels at muscle lengths larger than the optimal length (Table 1). In terms of the shortened muscles, however, the stimulation waveforms were considerably distorted to prevent the saturation of motoneuron PICs to achieve linear force relaxation in the descending stimulation phase under continuous stimulation conditions (see Figure 3A). Furthermore, the model motor unit did not fully reach the maximum force at the optimal muscle length without full PIC activation over the motoneuron dendrites in the ascending phase under continuous stimulation (Figure 8A). In contrast, the model motor unit reached all force levels within the full physiological range of the muscle length under discrete stimulation conditions (Figure 8B). This result indicates that at a high level of force generation, continuous current stimulation may sufficiently depolarize the membrane potential of the motoneuron dendrites to fully activate the dendritic PIC channels. However, the activation kinetics of dendritic PIC channels have been reported to be slow (i.e., ∼60 ms) (Carlin et al., 2000). Furthermore, attenuation of the alternating electrical signals tends to be much more severe than that of the direct electrical signals when propagating to the distal dendritic areas in the motoneuron (Kim and Jones, 2012). Thus, discrete current stimulation with a small pulse width (i.e., 0.5 ms) may effectively control the firing output of the motoneuron, thus avoiding the full activation of PIC channels over the motoneuron dendrites.

To this end, under the various conditions of brainstem neuromodulation and muscle length, the discrete stimulation protocol is more suitable than the continuous stimulation protocol for the precise modulation of the muscle force at the motor unit level.

Comparison to Previous Studies

The motor unit model considered in the present study accurately replicated the non-linear input–output relationship of the motoneuron and muscle unit, as separately investigated in previous experimental studies. The motoneuron model captured the acceleration in firing rate above the recruitment threshold in the ascending phase and the persistence of firing below the recruitment threshold in the descending phase during triangular current injection at the soma in the presence of monoamines (Lee and Heckman, 1998a, 1999). The muscle unit model also reproduced the sigmoidal force production in response to a linear increase and decrease in the stimulation frequency, as experimentally reported in previous studies (Rack and Westbury, 1969; Mrowczynski et al., 2006). The muscle length dependence of the non-linear behavior of the motoneuron and muscle fibers is consistent with the results obtained in a previous computational study involving an anatomically realistic model of the spinal motoneuron (Kim, 2017). Thus, the model motor unit constructed in the PyMUS software environment in this study may provide an efficient computational platform for the estimation of the stimulation waveforms required for force control at the motor unit level under a wide range of physiological conditions, including neuromodulation from the brainstem and sensory feedback originating from the muscle spindle.

To our knowledge, the waveform under intraspinal microstimulation has not yet been systematically investigated for muscle force control under various force profile and muscle length conditions. Seminal work has demonstrated, based on both extracellular and intracellular microelectrodes, that spinal motoneurons may be directly activated at the lowest threshold with a stimulating electrode near the initial segment (Gustafsson and Jankowska, 1976). Recently, a field model of the spinal cord coupled with realistic motoneuron models has been analyzed, suggesting that local cells may be selectively activated without activation of bypassing nerve fibers through asymmetric modulation of the duration and amplitude of the cathodic and anodic phases of the biphasic stimulus current pulse (McIntyre and Grill, 2002). This form of the stimulus current pulse may be applied to the discrete stimulation waveforms estimated in this study to improve the selectivity of the target motoneuron. In addition, the previous study has revealed that the stimulus amplitude must be increased to elicit firing in a one-to-one manner at high stimulation frequencies (McIntyre and Grill, 2002). This phenomenon was also observed in the current study with increasing force level (Figure 4). Furthermore, the motoneuronal output during extracellular microstimulation near the cell body has been evaluated with respect to several factors, including the dendritic active conductance and synaptic contact. The previous analysis has suggested that dendritic active conductance and synaptic contact activation may reduce the current threshold for firing initiation with little effect on the motoneuronal output (McIntyre and Grill, 2002). These results are consistent with those obtained in the present study. The amplitude of the current pulses was lowered with PIC activation over the motoneuron dendrites (Figure 6). The influence of synaptic activation was evaluated by increasing the excitatory synaptic conductance over the somatic compartment of the reduced motoneuron. This increase in the excitatory synaptic conductance reduced the amplitude of the current impulses with little difference observed between the stimulation waveforms (not shown).

The present study was conducted in regard to the direct activation of a spinal motoneuron via a microelectrode placed near the initial segment and cell body (Gustafsson and Jankowska, 1976). In contrast, indirect transsynaptic activation of the motoneuron has been verified to mostly occur when the microelectrode is located in dendritic regions likely including PIC channels (Gustafsson and Jankowska, 1976). The aforementioned dendritic location of the excitatory synaptic contacts may greatly facilitate the saturation of dendritic PIC channels, resulting in non-responsiveness of the dendrites to synaptic activation at a low force production level (Heckman et al., 2008). Thus, the effective type and waveform required for extracellular microstimulation might differ in the case of transsynaptic activation of the motoneuron and should be further investigated considering the addition of an additional excitatory synaptic conductance to the dendritic compartment of the reduced motoneuron model constructed in this study.

With the above considerations, the present study extends previous investigations by proposing the effective type (i.e., discrete) and waveform (i.e., exponential) of current stimulation for spinal motoneurons that may lead to the precise control of force production at the motor unit level.

Potential Contributions to Spinal Cord Modulation

Spinal cord stimulation has been applied to effectively restore complex motor functions through the modulation of neuronal networks in the spinal cord (Harkema et al., 2011; Megia Garcia et al., 2020). Intraspinal microstimulation with a multielectrode array may represent an efficient technique to stimulate and modulate the motor output of targeted motor units for the precise control of muscle forces and voluntary movements. From this perspective, this study may provide a basis to design the type and waveform of current stimulation for the precise control of voluntary muscle contractions. The predictions obtained in this study may be testable with recently developed experimental techniques. The force output of a single motor unit in response to current stimulation intracellularly injected at the soma of a spinal motoneuron has been characterized via in vivo mouse preparation (Manuel and Heckman, 2011). In addition, the multielectrode array interface developed for spinal cord stimulation (Meacham et al., 2011) and nerve cuff and flexible split ring electrode developed for selective stimulation of the motor nerve (Deurloo et al., 2003; Lee et al., 2017) could also be employed to test the stimulation waveforms required for force control predicted at the motor unit level in the present study. It should be noted that further analysis is needed for the indirect control of spinal motoneurons through current stimulation via the skin on the back (Gogeascoechea et al., 2020) or the epidural portion (Ahmed, 2016) of the spinal cord, which typically involves a variety of interneurons and afferents.

Modeling Considerations

The fundamental limitations of the reduced modeling approach of spinal motoneurons and muscle fibers have been fully addressed in previous studies (Kim et al., 2014, 2015). To estimate the stimulation pattern more physiologically, several issues should be resolved by improving the current motor unit model in future studies. First, the influence of the dendritic structure and PIC distribution on the motoneuron output should be reflected in the prediction of the stimulation pattern under pathological conditions (Elbasiouny et al., 2010). Second, the active membrane mechanisms over the dendrites of spinal motoneurons should be considered in detail, including calcium-activated potassium currents (Li and Bennett, 2007) and hyperpolarization-activated cation channels (Manuel et al., 2007). Third, the present study was conducted particularly considering slow fatigue-resistible motor units. To evaluate other types of motor units, such as fast fatigable motor units, the stimulation waveforms estimated in this study should be adjusted considering the discrepancy in electrical and mechanical properties, particularly in regard to the somatic input resistance, system time constant, rheobase current, afterhyperpolarization potential, and PIC dynamics in spinal motoneurons (Lee and Heckman, 1999; Heckman and Enoka, 2012), twitch rate and amplitude, progressive force reduction phenomenon (i.e., the sag observed over short contraction time <2 s and fatigue over a long contraction time >2 min), and length– and velocity–tension properties of muscle fibers (Brown et al., 1999; Burke, 2011). Fourth, the present study did not consider the shape of the current pulse (Anderson et al., 2019). Thus, the waveforms predicted in the current study might vary according to the shape of the current pulse. Finally, further investigation is required in terms of the stimulation waveforms to control the force production with the whole muscle, considering the organization of spinal motoneurons exhibiting different anatomical and electrical properties (Zengel et al., 1985; Capogrosso et al., 2013). In addition, the influence of spinal interneurons and peripheral afferents should be considered, particularly in the case of skin or epidural stimulation of the spinal cord (Harkema et al., 2011; Takei and Seki, 2013; Gogeascoechea et al., 2020).

Concluding Remarks

In principle, the waveform obtained with a microelectrode near the motoneuron cell body to realize isometric force control of the motor unit strongly depends on the stimulation type, neuromodulatory input, and muscle length. Compared to continuous current stimulation, discrete current stimulation more suitably prevents the non-linearities induced by the full activation of PIC-generating channels over the motoneuron dendrites, including muscle spindle feedback. Thus, motor unit function modulation may be more effectively achieved by stimulating spinal motoneurons with short-width current pulses. The model-based approach and waveform equations developed in the present study may provide a basis for the design of stimulation patterns applicable to spinal motoneurons to better understand neuromuscular physiology and neuromodulation of movement disorders via a neural interface.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

HK conceived and designed the study and wrote the manuscript. YJ and HK performed the simulations and analyzed the data. Both authors contributed to the article and approved the submitted version.

Funding

DGIST R&D Program of the Ministry of Science and ICT, South Korea (20-BT-06 and 21-BT-06).

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.

Acknowledgments

The authors thank Minjung Kim for her help with updating the PyMUS software for this study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2021.645984/full#supplementary-material

Supplementary File 1 | Equations and descriptions of the model motor unit.

Supplementary File 2 | Computer codes of the models and simulations.

Supplementary File 3 | Dataset of the current stimulation levels applied to the motoneuron for force production of the muscle unit at the various speeds and levels, as shown in Figures 3–7.

Supplementary File 4 | Production of 100% of the maximal force by the model motor unit, as shown in Figure 4, at the optimal muscle length under discrete current stimulation.

Supplementary File 5 | Production of 20% of the maximal force by the model motor unit, as shown in Figure 6, at the optimal muscle length under discrete current stimulation conditions.

Supplementary File 6 | Production of 100% of the maximal force by the model motor unit, as shown in Figure 7, at the optimal muscle length under discrete current stimulation conditions.

Abbreviations

VI, voltage–current; IV, current–voltage; VODE, variable-coefficient ordinary differential equation; LSODE, livermore solver for ordinary differential equations; PIC, persistent inward current; PyMUS, python-based motor unit simulator.

Footnotes

1. ^ https://github.com/NMSL-DGIST/PyMUS

References

Ahmed, Z. (2016). Modulation of gamma and alpha spinal motor neurons activity by trans-spinal direct current stimulation: effects on reflexive actions and locomotor activity. Physiol. Rep. 4:e12696. doi: 10.14814/phy2.12696

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, D. J., Kipke, D. R., Nagel, S. J., Lempka, S. F., Machado, A. G., Holland, M. T., et al. (2019). Intradural spinal cord stimulation: performance modeling of a new modality. Front. Neurosci. 13:253. doi: 10.3389/fnins.2019.00253

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Brown, I. E., Cheng, E. J., and Loeb, G. E. (1999). Measured and modeled properties of mammalian skeletal muscle. II. The effects of stimulus frequency on force-length and force-velocity relationships. J. Muscle Res. Cell Motil. 20, 627–643.

Google Scholar

Burke, R. E. (2011). Motor Units: Anatomy, Physiology, and Functional Organization. Comprehensive Physiology. Hoboken, NJ: John Wiley & Sons, Inc.

Google Scholar

Capogrosso, M., Wenger, N., Raspopovic, S., Musienko, P., Beauparlant, J., Bassi Luciani, L., et al. (2013). A computational model for epidural electrical stimulation of spinal sensorimotor circuits. J. Neurosci. 33, 19326–19340. doi: 10.1523/jneurosci.1688-13.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlin, K. P., Jones, K. E., Jiang, Z., Jordan, L. M., and Brownstone, R. M. (2000). Dendritic L-type calcium currents in mouse spinal motoneurons: implications for bistability. Eur. J. Neurosci. 12, 1635–1646. doi: 10.1046/j.1460-9568.2000.00055.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Deurloo, K. E., Holsheimer, J., and Bergveld, P. (2003). Fascicular selectivity in transverse stimulation with a nerve cuff electrode: a theoretical approach. Neuromodulation 6, 258–269. doi: 10.1046/j.1525-1403.2003.03034.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Elbasiouny, S. M., Amendola, J., Durand, J., and Heckman, C. J. (2010). Evidence from computer simulations for alterations in the membrane biophysical properties and dendritic processing of synaptic inputs in mutant superoxide dismutase-1 motoneurons. J. Neurosci. 30, 5544–5558. doi: 10.1523/jneurosci.0434-10.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Emancipator, K., and Kroll, M. H. (1993). A quantitative measure of nonlinearity. Clin. Chem. 39, 766–772. doi: 10.1093/clinchem/39.5.766

CrossRef Full Text | Google Scholar

Gogeascoechea, A., Kuck, A., van Asseldonk, E., Negro, F., Buitenweg, J. R., Yavuz, U. S., et al. (2020). Interfacing with alpha motor neurons in spinal cord injury patients receiving trans-spinal electrical stimulation. Front. Neurol. 11:493. doi: 10.3389/fneur.2020.00493

PubMed Abstract | CrossRef Full Text | Google Scholar

Goslow, G. E. Jr., Reinking, R. M., and Stuart, D. G. (1973). The cat step cycle: hind limb joint angles and muscle lengths during unrestrained locomotion. J. Morphol. 141, 1–41. doi: 10.1002/jmor.1051410102

PubMed Abstract | CrossRef Full Text | Google Scholar

Gustafsson, B., and Jankowska, E. (1976). Direct and indirect activation of nerve cells by electrical pulses applied extracellularly. J. Physiol. 258, 33–61. doi: 10.1113/jphysiol.1976.sp011405

PubMed Abstract | CrossRef Full Text | Google Scholar

Harkema, S., Gerasimenko, Y., Hodes, J., Burdick, J., Angeli, C., Chen, Y., et al. (2011). Effect of epidural stimulation of the lumbosacral spinal cord on voluntary movement, standing, and assisted stepping after motor complete paraplegia: a case study. Lancet 377, 1938–1947. doi: 10.1016/s0140-6736(11)60547-3

CrossRef Full Text | Google Scholar

Heckman, C. J., and Enoka, R. M. (2012). Motor unit. Compr. Physiol. 2, 2629–2682.

Google Scholar

Heckman, C. J., Hyngstrom, A. S., and Johnson, M. D. (2008). Active properties of motoneurone dendrites: diffuse descending neuromodulation, focused local inhibition. J. Physiol. 586, 1225–1231. doi: 10.1113/jphysiol.2007.145078

PubMed Abstract | CrossRef Full Text | Google Scholar

Heckman, C. J., Lee, R. H., and Brownstone, R. M. (2003). Hyperexcitable dendrites in motoneurons and their neuromodulatory control during motor behavior. Trends Neurosci. 26, 688–695. doi: 10.1016/j.tins.2003.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Heckmann, C. J., Gorassini, M. A., and Bennett, D. J. (2005). Persistent inward currents in motoneuron dendrites: implications for motor output. Muscle Nerve 31, 135–156. doi: 10.1002/mus.20261

PubMed Abstract | CrossRef Full Text | Google Scholar

Hochman, S., and McCrea, D. A. (1994). Effects of chronic spinalization on ankle extensor motoneurons. II. Motoneuron electrical properties. J. Neurophysiol. 71, 1468–1479. doi: 10.1152/jn.1994.71.4.1468

PubMed Abstract | CrossRef Full Text | Google Scholar

Holsheimer, J. (1998). Concepts and methods in neuromodulation and functional electrical stimulation: an introduction. Neuromodulation 1, 57–61. doi: 10.1111/j.1525-1403.1998.tb00018.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hornby, T. G., McDonagh, J. C., Reinking, R. M., and Stuart, D. G. (2002). Motoneurons: a preferred firing range across vertebrate species? Muscle Nerve 25, 632–648. doi: 10.1002/mus.10105

PubMed Abstract | CrossRef Full Text | Google Scholar

Hounsgaard, J., Hultborn, H., Jespersen, B., and Kiehn, O. (1984). Intrinsic membrane properties causing a bistable behaviour of alpha-motoneurones. Exp. Brain Res. 55, 391–394.

Google Scholar

Hounsgaard, J., Hultborn, H., Jespersen, B., and Kiehn, O. (1988). Bistability of alpha-motoneurones in the decerebrate cat and in the acute spinal cat after intravenous 5-hydroxytryptophan. J. Physiol. 405, 345–367. doi: 10.1113/jphysiol.1988.sp017336

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, H. (2017). Muscle length-dependent contribution of motoneuron Cav1.3 channels to force production in model slow motor unit. J. Appl. Physiol. 123, 88–105. doi: 10.1152/japplphysiol.00491.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, H., and Jones, K. E. (2012). The retrograde frequency response of passive dendritic trees constrains the nonlinear firing behaviour of a reduced neuron model. PLoS One 7:e43654. doi: 10.1371/journal.pone.0043654

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, H., Jones, K. E., and Heckman, C. J. (2014). Asymmetry in signal propagation between the soma and dendrites plays a key role in determining dendritic excitability in motoneurons. PLoS One 9:e95454. doi: 10.1371/journal.pone.0095454

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, H., and Kim, M. (2018). PyMUS: python-based simulation software for virtual experiments on motor unit system. Front. Neuroinform. 12:15. doi: 10.3389/fninf.2018.00015

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, H., Sandercock, T. G., and Heckman, C. J. (2015). An action potential-driven model of soleus muscle activation dynamics for locomotor-like movements. J. Neural. Eng. 12:046025. doi: 10.1088/1741-2560/12/4/046025

CrossRef Full Text | Google Scholar

Lee, R. H., and Heckman, C. J. (1998b). Bistability in spinal motoneurons in vivo: systematic variations in rhythmic firing patterns. J. Neurophysiol. 80, 572–582. doi: 10.1152/jn.1998.80.2.572

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, R. H., and Heckman, C. J. (1998a). Bistability in spinal motoneurons in vivo: systematic variations in persistent inward currents. J. Neurophysiol. 80, 583–593. doi: 10.1152/jn.1998.80.2.583

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, R. H., and Heckman, C. J. (1999). Paradoxical effect of QX-314 on persistent inward currents and bistable behavior in spinal motoneurons in vivo. J. Neurophysiol. 82, 2518–2527. doi: 10.1152/jn.1999.82.5.2518

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, R. H., Kuo, J. J., Jiang, M. C., and Heckman, C. J. (2003). Influence of active dendritic currents on input-output processing in spinal motoneurons in vivo. J. Neurophysiol. 89, 27–39. doi: 10.1152/jn.00137.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S., Sheshadri, S., Xiang, Z. L., Delgado-Martinez, I., Xue, N., Sun, T., et al. (2017). Selective stimulation and neural recording on peripheral nerves using flexible split ring electrodes. Sens. Actuators B Chem. 242, 1165–1170. doi: 10.1016/j.snb.2016.09.127

CrossRef Full Text | Google Scholar

Li, X., and Bennett, D. J. (2007). Apamin-sensitive calcium-activated potassium currents (SK) are activated by persistent calcium currents in rat motoneurons. J. Neurophysiol. 97, 3314–3330. doi: 10.1152/jn.01068.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Manuel, M., and Heckman, C. J. (2011). Adult mouse motor units develop almost all of their force in the subprimary range: a new all-or-none strategy for force recruitment? J. Neurosci. 31, 15188–15194.

Google Scholar

Manuel, M., Meunier, C., Donnet, M., and Zytnicki, D. (2007). Resonant or not, two amplification modes of proprioceptive inputs by persistent inward currents in spinal motoneurons. J. Neurosci. 27, 12977–12988. doi: 10.1523/jneurosci.3299-07.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

McIntyre, C. C., and Grill, W. M. (2000). Selective microstimulation of central nervous system neurons. Ann. Biomed. Eng. 28, 219–233. doi: 10.1114/1.262

CrossRef Full Text | Google Scholar

McIntyre, C. C., and Grill, W. M. (2002). Extracellular stimulation of central neurons: influence of stimulus waveform and frequency on neuronal output. J. Neurophysiol. 88, 1592–1604. doi: 10.1152/jn.2002.88.4.1592

PubMed Abstract | CrossRef Full Text | Google Scholar

Meacham, K. W., Guo, L., Deweerth, S. P., and Hochman, S. (2011). Selective stimulation of the spinal cord surface using a stretchable microelectrode array. Front. Neuroeng. 4:5. doi: 10.3389/fneng.2011.00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Megia Garcia, A., Serrano-Munoz, D., Taylor, J., Avendano-Coy, J., and Gomez-Soriano, J. (2020). Transcutaneous spinal cord stimulation and motor rehabilitation in spinal cord injury: a systematic review. Neurorehabil. Neural Repair 34, 3–12.

Google Scholar

Mrowczynski, W., Celichowski, J., and Krutki, P. (2006). Interspecies differences in the force-frequency relationship of the medial gastrocnemius motor units. J. Physiol. Pharmacol. 57, 491–501.

Google Scholar

Mushahwar, V. K., and Horch, K. W. (1998). Selective activation and graded recruitment of functional muscle groups through spinal cord stimulation. Ann. N. Y. Acad. Sci. 860, 531–535. doi: 10.1111/j.1749-6632.1998.tb09096.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Perreault, E. J., Heckman, C. J., and Sandercock, T. G. (2003). Hill muscle model errors during movement are greatest within the physiologically relevant range of motor unit firing rates. J. Biomech. 36, 211–218. doi: 10.1016/s0021-9290(02)00332-9

CrossRef Full Text | Google Scholar

Rack, P. M., and Westbury, D. R. (1969). The effects of length and stimulus rate on tension in the isometric cat soleus muscle. J. Physiol. 204, 443–460. doi: 10.1113/jphysiol.1969.sp008923

PubMed Abstract | CrossRef Full Text | Google Scholar

Radivojevic, M., Jackel, D., Altermatt, M., Muller, J., Viswam, V., Hierlemann, A., et al. (2016). Electrical identification and selective microstimulation of neuronal compartments based on features of extracellular action potentials. Sci. Rep. 6:31332.

Google Scholar

Sandercock, T. G., and Heckman, C. J. (1997). Force from cat soleus muscle during imposed locomotor-like movements: experimental data versus Hill-type model predictions. J. Neurophysiol. 77, 1538–1552. doi: 10.1152/jn.1997.77.3.1538

PubMed Abstract | CrossRef Full Text | Google Scholar

Scott, S. H., Brown, I. E., and Loeb, G. E. (1996). Mechanics of feline soleus: I. Effect of fascicle length and velocity on force output. J. Muscle Res. Cell Motil. 17, 207–219. doi: 10.1007/bf00124243

PubMed Abstract | CrossRef Full Text | Google Scholar

Segev, I., Fleshman, J. W. Jr., and Burke, R. E. (1990). Computer simulation of group Ia EPSPs using morphologically realistic models of cat alpha-motoneurons. J. Neurophysiol. 64, 648–660. doi: 10.1152/jn.1990.64.2.648

PubMed Abstract | CrossRef Full Text | Google Scholar

Shames, D. M., Baker, A. J., Weiner, M. W., and Camacho, S. A. (1996). Ca(2+)-force relationship of frog skeletal muscle: a dynamic model for parameter estimation. Am. J. Physiol. 271(6 Pt 1), C2062–C2071.

Google Scholar

Spruston, N., Jaffe, D. B., and Johnston, D. (1994). Dendritic attenuation of synaptic potentials and currents: the role of passive membrane properties. Trends Neurosci. 17, 161–166. doi: 10.1016/0166-2236(94)90094-9

CrossRef Full Text | Google Scholar

Sunshine, M. D., Cho, F. S., Lockwood, D. R., Fechko, A. S., Kasten, M. R., and Moritz, C. T. (2013). Cervical intraspinal microstimulation evokes robust forelimb movements before and after injury. J. Neural Eng. 10:036001. doi: 10.1088/1741-2560/10/3/036001

CrossRef Full Text | Google Scholar

Sunshine, M. D., Ganji, C. N., Reier, P. J., Fuller, D. D., and Moritz, C. T. (2018). Intraspinal microstimulation for respiratory muscle activation. Exp. Neurol. 302, 93–103. doi: 10.1016/j.expneurol.2017.12.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Takei, T., and Seki, K. (2013). Spinal premotor interneurons mediate dynamic and static motor commands for precision grip in monkeys. J. Neurosci. 33, 8850–8860. doi: 10.1523/jneurosci.4032-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner, F. B., Mignardot, J. B., Le Goff-Mignardot, C. G., Demesmaeker, R., Komi, S., Capogrosso, M., et al. (2018). Targeted neurotechnology restores walking in humans with spinal cord injury. Nature 563, 65–71.

Google Scholar

Wang, Q., Millard, D. C., Zheng, H. J., and Stanley, G. B. (2012). Voltage-sensitive dye imaging reveals improved topographic activation of cortex in response to manipulation of thalamic microstimulation parameters. J. Neural Eng. 9:026008. doi: 10.1088/1741-2560/9/2/026008

CrossRef Full Text | Google Scholar

Warman, E. N., Grill, W. M., and Durand, D. (1992). Modeling the effects of electric fields on nerve fibers: determination of excitation thresholds. IEEE Trans. Biomed. Eng. 39, 1244–1254. doi: 10.1109/10.184700

CrossRef Full Text | Google Scholar

Westerblad, H., and Allen, D. G. (1994). The role of sarcoplasmic reticulum in relaxation of mouse muscle; effects of 2,5-di(tert-butyl)-1,4-benzohydroquinone. J. Physiol. 474, 291–301. doi: 10.1113/jphysiol.1994.sp020022

PubMed Abstract | CrossRef Full Text | Google Scholar

Winters, T. M., Takahashi, M., Lieber, R. L., and Ward, S. R. (2011). Whole muscle length-tension relationships are accurately modeled as scaled sarcomeres in rabbit hindlimb muscles. J. Biomech. 44, 109–115. doi: 10.1016/j.jbiomech.2010.08.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Zengel, J. E., Reid, S. A., Sypert, G. W., and Munson, J. B. (1985). Membrane electrical properties and prediction of motor-unit type of medial gastrocnemius motoneurons in the cat. J. Neurophysiol. 53, 1323–1344. doi: 10.1152/jn.1985.53.5.1323

PubMed Abstract | CrossRef Full Text | Google Scholar