OUP user menu

Human spinal locomotor control is based on flexibly organized burst generators

Simon M. Danner, Ursula S. Hofstoetter, Brigitta Freundl, Heinrich Binder, Winfried Mayr, Frank Rattay, Karen Minassian
DOI: http://dx.doi.org/10.1093/brain/awu372 awu372 First published online: 12 January 2015


Constant drive provided to the human lumbar spinal cord by epidural electrical stimulation can cause local neural circuits to generate rhythmic motor outputs to lower limb muscles in people paralysed by spinal cord injury. Epidural spinal cord stimulation thus allows the study of spinal rhythm and pattern generating circuits without their configuration by volitional motor tasks or task-specific peripheral feedback. To reveal spinal locomotor control principles, we studied the repertoire of rhythmic patterns that can be generated by the functionally isolated human lumbar spinal cord, detected as electromyographic activity from the legs, and investigated basic temporal components shared across these patterns. Ten subjects with chronic, motor-complete spinal cord injury were studied. Surface electromyographic responses to lumbar spinal cord stimulation were collected from quadriceps, hamstrings, tibialis anterior, and triceps surae in the supine position. From these data, 10-s segments of rhythmic activity present in the four muscle groups of one limb were extracted. Such samples were found in seven subjects. Physiologically adequate cycle durations and relative extension- and flexion-phase durations similar to those needed for locomotion were generated. The multi-muscle activation patterns exhibited a variety of coactivation, mixed-synergy and locomotor-like configurations. Statistical decomposition of the electromyographic data across subjects, muscles and samples of rhythmic patterns identified three common temporal components, i.e. basic or shared activation patterns. Two of these basic patterns controlled muscles to contract either synchronously or alternatingly during extension- and flexion-like phases. The third basic pattern contributed to the observed muscle activities independently from these extensor- and flexor-related basic patterns. Each bifunctional muscle group was able to express both extensor- and flexor-patterns, with variable ratios across the samples of rhythmic patterns. The basic activation patterns can be interpreted as central drives implemented by spinal burst generators that impose specific spatiotemporally organized activation on the lumbosacral motor neuron pools. Our data thus imply that the human lumbar spinal cord circuits can form burst-generating elements that flexibly combine to obtain a wide range of locomotor outputs from a constant, repetitive input. It may be possible to use this flexibility to incorporate specific adaptations to gait and stance to improve locomotor control, even after severe central nervous system damage.

  • central pattern generation
  • epidural spinal cord stimulation
  • human
  • modular organization
  • spinal cord injury


In quadrupedal mammals, neural networks in the lumbar spinal cord, termed central pattern generators, provide the basic sequential neural activity to hindlimb motor pools as required for stepping, even when experimentally isolated from supraspinal and sensory modulation (Brown, 1911; Grillner, 1981, 2006). Evidence for spinal locomotor pattern generation in humans is inherently indirect. Automatic rhythmic activation of lower-extremity muscles, i.e. spinal myoclonus, has been observed in subjects with chronic spinal cord injury (Bussel et al., 1988; Calancie et al., 1994; Calancie, 2006). Further, the lumbar spinal cord circuitry can be activated by epidural stimulation to generate rhythmic motor outputs in motor-complete spinal cord injured individuals lying supine (Dimitrijevic et al., 1998; Shapkova and Schomburg, 2001; Gerasimenko et al., 2002; Minassian et al., 2004). The current view is that some elements of the human lumbar locomotor circuitry are similar to those of other mammals (Duysens and Van de Crommert, 1998; Grillner, 2011; Guertin, 2013). Yet, no conclusions have been drawn on its functional organization.

There are two leading conceptual models for the organization of the mammalian locomotor pattern generator (Grillner and Zangger, 1975; McCrea and Rybak, 2008; Guertin, 2009). One model is the two-level hierarchical central pattern generator with a half-centre-like rhythm generator acting as a ‘master-clock’ for pattern formation networks that, in turn, distribute excitation to the different motor neuron pools (Burke et al., 2001; Rybak et al., 2006; Zhong et al., 2012). The second model is a network of coupled unit burst generators, each forming a separate rhythmogenic module (Grillner, 1981, 2003; Hägglund et al., 2013). The consensus is that the central pattern generating networks, one per hemicord and limb, show a modular organization where burst-generating elements are combined to produce integrated motor outputs (McCrea and Rybak, 2008; Hägglund et al., 2013).

The modular organization of human motor control has been studied by statistical analysis of motor patterns recorded during specific locomotor tasks (Giszter and Hart, 2013). A motor module generates basic control signals to produce functionally relevant patterns of muscle activation (Dominici et al., 2011; Fox et al., 2013). The use of a small set of modules by the nervous system is thought to simplify the control of stereotyped motor behaviour (Clark et al., 2010; Fox et al., 2013; Zelik et al., 2014). These motor modules agree with the concepts of functional reflexes, motor primitives, and muscle synergies while some elements of modularity conform to central pattern generator models (Dominici et al., 2011; Duysens et al., 2013; Giszter and Hart, 2013). Decomposition of EMG activity recorded during walking identified two basic activation patterns across the developmental period of human gait, an ‘extension’ and a ‘swing’ pattern, that were retained and tuned in adult locomotion (Dominici et al., 2011) and were suggested to be controlled at spinal level (Grillner, 2011).

Here, we used epidural spinal cord stimulation (SCS) to generate rhythmic motor patterns in the absence of a volitional locomotor task or task-specific peripheral feedback. We analysed the characteristics and modularity of these SCS-induced rhythmic motor outputs in individuals with motor complete spinal cord injury to uncover the functional organization of spinal locomotor control.

Materials and methods

Under approval of the Ethics Committee of the City of Vienna, Austria, in accordance with the Declaration of Helsinki, and with signed informed consent, 10 otherwise healthy adults aged 28.2 ± 11.8 (mean ± SD) years with post-traumatic and clinically motor-complete spinal cord injury in a chronic condition (≥1 year post-injury) were studied. Stretch and cutaneo-muscular reflexes of the legs were preserved in all. The neurological status was evaluated according to the International Standards for Neurological Classification of Spinal Cord Injury (Table 1).

View this table:
Table 1

Demographic data and neurological classification of spinal cord injury

SubjectGenderAge (years)Years post injuryVertebral level(s) of fractureAISNeurological level of injuryTotal motor score*Pin prick scoreLight touch score
1F20.13.3T4, T5AT4504844
2M24.31.0T3, T4AT5505050
3M58.43.3T7, T8BT9508886
4M25.31.6C5, C6BC8276868
5F33.02.7T4–T6, T10AT5505054
9M21.22.6T7, T8AT7506161
10M18.03.0C4, C5, T7BC6166464
  • AIS = American Spinal Injury Association Impairment Scale; *Lower-limb motor scores were 0 in all.

Additionally, a neurophysiological evaluation with surface-EMG recordings from the adductors, quadriceps, hamstrings, tibialis anterior and triceps surae muscles bilaterally was applied (Sherwood et al., 1996; McKay et al., 2012). None of the subjects had the ability to induce EMG activity by the attempt of isolated unilateral ankle dorsi- and plantar-flexion. Non-specific residual excitatory influence was tested by the activation of muscles below the lesion in response to a reinforcement manoeuvre, a forceful neck flexion against resistance (Table 2). Residual inhibitory supraspinal influence was assessed by evoking foot withdrawal by non-noxious mechanical plantar stimulation without and then with the volitional attempt to suppress the response (Table 2).

View this table:
Table 2

Neurophysiological evaluation of translesional influence

Reinforcement (neck flexion)Volitional suppression of foot withdrawal reflex10-s sample of rhythmic EMG patterns generated by SCS
SubjectNumber of activated muscles (max = 5 + 5)% of reduction*Locomotor-likeRhythmic co-activation across all musclesOther mixed-synergy patterns
  • *Reduction of root mean squares of EMG activity; **Subject 10 did not respond to plantar stimulation.

  • x = pattern present.

Spinal cord stimulation system

All subjects had SCS systems implanted for the control of spinal spasticity affecting the legs (Pinter et al., 2000). The system included a linear electrode array (Pisces-Quad electrode, Model 3487A, Medtronic), connected to either an external test stimulator (Model 3625 Test Stimulator, Medtronic), or an implanted pulse generator (Itrel 3, Model 7425, Medtronic). The electrode array consisted of four evenly spaced ring electrodes, each 1.27 mm in diameter and 3 mm long, with an interelectrode spacing of 6 mm. For their identification, these contacts will be referred to as 0 to 3, 0 being the most rostral. The electrode array was placed in the dorsal epidural space. Target position was over the midline of the lumbar spinal cord and, across the subjects, the vertebral levels of the electrode positions ranged from T11 to L1. The stimulators delivered monophasic rectangular constant-voltage pulses of 210-µs width, each followed by a longer second phase of small amplitude to adjust charge balance for each stimulus and avoid delivery of direct current (Rattay et al., 2000). Available stimulation frequencies were 2–130 Hz, and intensities up to 10.5 V. Electrode impedance of the SCS system is ∼1000 Ω (Minassian et al., 2004). Each electrode of the array could be set as cathode (‘–’), anode (‘+’), or to be inactive, allowing for various bipolar electrode combinations. Monopolar stimulation was carried out with one of the electrodes selected as cathode and the active area on the implanted pulse generator case (‘c’) as anode. The rostro-caudal positions of the active cathode with respect to the spinal cord segments were estimated based on the relative thresholds of posterior root-muscle reflexes electromyographically recorded as stimulus time-locked compound muscle action potentials from quadriceps and triceps surae elicited by 2-Hz stimulation. Lower thresholds of posterior root-muscle reflexes in quadriceps than in triceps surae indicate cathode positions over the L2–L4 lumbar spinal cord segments, equal thresholds a cathode over the L5–S2 lumbosacral segments, and lower thresholds in triceps surae a cathode caudal to the S2 segment (Minassian et al., 2007).

Electromyographic recordings

SCS-induced surface EMG activity was acquired from the quadriceps, hamstrings, tibialis anterior, and triceps surae muscles bilaterally using pairs of silver-silver chloride recording electrodes (Intec Medizintechnik), each placed centrally over the muscle bellies and oriented along the long axis of the muscles with an interelectrode distance of 3 cm. Additional EMG recordings were taken from the lumbar paraspinal trunk muscles. Abrasive paste (Nuprep, Weaver and Company) was used for skin preparation to reduce contact resistance. EMG signals were amplified (Grass Instruments) with a gain of 2000 and filtered to a bandwidth of 30–700 Hz. Data were digitized at 2002 samples per second and channel using a Codas ADC system (Dataq Instruments).

Study protocol

Subjects were lying supine on an examination bed. For a given combination of active electrodes, stimulation was initially applied at the lowest frequency available, and the stimulus intensity was increased in 1 V increments until posterior root-muscle reflexes were detected electromyographically in quadriceps, hamstrings, tibialis anterior, and triceps surae bilaterally. At this intensity, the stimulation frequency was increased step-wise, including frequencies around 5, 10, 16, 21, 25, 31 and 40 Hz. The exact stimulation frequencies were monitored by recording stimulation artefacts via the paraspinal EMG electrodes. The frequency variation was repeated for higher intensities up to a maximum of 10 V. Subsequently, stimulus intensity and frequency variation were repeated for different active electrode combinations. For each subject, data from two separate recording sessions were analysed.

Data selection

Ten-second segments of unilateral rhythmic EMG activity elicited by SCS in quadriceps, hamstrings, tibialis anterior, and triceps surae with established and consistent rhythmicity were selected from the data pool on the basis of the quality of rhythmicity across the four unilateral muscles. Right and left legs were analysed independently. From each data section with unchanged stimulation parameters, all 10-s segments, each shifted by one interstimulus interval, were extracted. The SCS-induced EMG activities were series of stimulus time-related posterior root-muscle reflexes (Minassian et al., 2004) that can be normally identified as separate responses for SCS frequencies of up to 40 Hz (Figs. 1A and B). Envelopes of the EMG data were thus calculated from the peak-to-peak amplitudes of the series of posterior root-muscle reflexes (Fig. 1A). The ranking parameter r of each 10-s segment of rhythmic EMG activity was calculated as r = max(fQ)·max(fHam)·max(fTA)·max(fTS), where f is the frequency spectrum between 0.2 Hz and 3 Hz of the EMG envelope of the respective muscles (Q, quadripceps; Ham, hamstrings; TA, tibialis anterior; TS, triceps surae). All 10-s segments of a recording session were sorted according to their r-values. Samples with the highest r-values of each recording session were visually inspected and finally selected for further analysis so that only one example of a given EMG pattern was chosen from a recording segment with unchanged stimulation parameters.

Figure 1

Non-negative matrix factorization (NMF): data preparation and interpretation of the basic patterns and weights. Construction of EMG profiles that served as input data for the NMF, shown for tibialis anterior (AE). (A) Envelopes of the selected rhythmic EMG segments were calculated from peak-to-peak amplitudes of series of responses. (B) Section of same EMG trace with enlarged time-scale shows that bursts are comprised of modulating reflex-compound muscle action potentials time-related to the stimulus pulses (stim.). (C) Low-pass filtered envelopes of the same tibialis anterior bursts were used to define the cycles and their separation into extension- (Ext) and flexion-like (Fl) phases. (D) Filtered envelopes of all complete cycles within the 10-s sample. These subdivisions were extrapolated to 100 data points for the Ext and Fl phase each, amplitude normalized, and averaged to construct the EMG profile (E). Calculated for all muscles and samples, these EMG profiles (X) served as the input to the NMF (F). The NMF identifies a preselected number of k basic patterns (P) and their weights (W), the latter specifying the contribution of each basic pattern to a given EMG profile. The various EMG profiles (black solid lines) can be approximated (dotted lines) by the weighted sum of the small number of basic patterns, which represent common underlying temporal components. pk,j are individual basic patterns, and wk,j their weights, consecutively numbered by index j.

Decomposition of rhythmic EMG activities

A dimensionality reduction technique was applied to the EMG profile of each muscle of the selected 10-s segments to identify basic patterns underlying the muscle activation. To construct the EMG profiles, the EMG envelopes were low-pass filtered (zero-phase Equiripple filter: Fpass = 0.125, Fstop = 0.167 normalized frequency units, Apass = 2, Astop = 5 dB; Fig. 1C) and separated into subdivisions, each covering a complete cycle of rhythmic activity. Periods of tibialis anterior activity were used to define the flexion-like phases based on the hypothesis that the uniarticular tibialis anterior would be explicitly controlled by a flexor-related drive. Indeed, previous findings have shown that phases of tibialis anterior activity separated SCS-induced rhythmic activity into two different functionally organized phases (Minassian et al., 2004, 2007; Danner et al., 2013). The cycles of rhythmic activity were defined as the intervals between the ends of successive tibialis anterior bursts. The flexion phase of rhythmic activity was defined as the time between the onset and end of each tibialis anterior burst, and the remaining phase of the cycle as the extension phase. The onset of a tibialis anterior burst was the time when the data of the low-pass filtered EMG envelope exceeded the local minimum preceding the burst by 50% of the difference between the peak value of the burst and this local minimum (Fig. 1C). The end of a tibialis anterior burst was defined accordingly. Each subdivision of the EMG envelopes was extrapolated to 200 data points, 100 values for the extension and the flexion phase, each. The phase normalization was required because of the considerable variation of the relative extension- and flexion-phase durations observed across samples. These data were normalized to the peak-amplitude. The EMG profiles resulted from averaging the individual traces of all complete cycles detected within the 10-s segments (Fig. 1D and E).

Each of the selected 10-s segments of rhythmic activity provided four EMG profiles. The EMG profiles of all selected data segments were assembled into a matrix X (with each column being an EMG profile) that served as input for the non-negative matrix factorization (NMF; cf. supporting online material for Dominici et al., 2011). The aim was to test whether the EMG profiles of the rhythmic lower-limb activities could be closely reconstructed by linear combinations of a small number of non-negative basic activation patterns (Fig. 1F).

NMF is a method that factorizes a matrix X into two matrices, P and W, NMF(X, k) → P·W, so that X = P·W + e, and |e| → min, where P consists of a predefined number of k basic activation patterns (activation timing profiles, primitives, shared activation pattern) and W is the matrix of weights (loadings). X is an n × m, P an n × k and W a k × m matrix, with n = 200, i.e. the count of data points (time samples) of each EMG profile, and m the number (index) of the EMG profiles. All matrices are non-negative. Models using 1–8 basic activation patterns were analysed here. The coefficient of determination (R2) and the Akaike Information Criterion with correction for finite sample sizes (AICC; Akaike, 1974; Burnham and Anderson, 2002) were calculated to assess the model’s goodness of fit. The AICC is rooted in information theory and represents an estimate of the relative information loss by representing the process underlying the data with a given model (Burnham and Anderson, 2002).

Statistical analysis

Mean values were compared using Student’s t-tests or analyses of variance. Categorical data were analysed using Pearson’s χ2-test. Pearson’s product moment correlations were calculated to investigate the relation of two scalar values and linear regressions were calculated to investigate whether one scalar value predicts the outcome of another. Similarly, logistic regressions were used to analyse scalar predictors on categorical values. If model assumptions were not met, their respective non-parametric equivalents were used. All post hoc tests were Bonferroni corrected. An α-error of P < 0.05 was regarded as significant.


Sustained epidural stimulation induced rhythmic EMG activity in quadriceps, hamstrings, tibialis anterior, and triceps surae with various synchronous and reciprocal relationships between the muscles. We identified 39 samples of 10 s duration with rhythmic activity generated in all recorded muscle groups of one limb. Figure 2 gives four examples with different patterns and cycle durations. The rhythmic patterns were found in Subjects 1–7 (Tables 1 and 2). Some rhythmic activity was observed in Subjects 8–10 as well, but failed to persist for the required 10 s or was not present in all four muscles. Remaining subclinical supraspinal influence over lumbosacral spinal motor excitability was not related to the occurrence of rhythmic muscle activation (cf. Table 2). In 23 of 39 samples, rhythmic activity was expressed unilaterally only, with the muscles of the contralateral leg showing either unmodulated or irregularly modulated motor output. The unilateral samples of rhythmic activity were found on the side with lower posterior root-muscle reflex thresholds in 16 cases and on the side with the higher thresholds only once.

Figure 2

Example 10-s segments of rhythmic EMG activity evoked by epidural spinal cord stimulation in quadriceps, hamstrings, tibialis anterior, and triceps surae. (A) Subject 2, right side, electrode pairing and polarity of 0+ and 1−, 7 V, 30.1 Hz; (B) Subject 3, left side, electrodes 0+ and 3−, 10 V, 22.5 Hz; (C) Subject 4, left side, electrodes 0+ and 3−, 6 V, 31.5 Hz; (D) Subject 2, right side, electrodes 1+ and 3−, 10 V, 27.7 Hz.

Rhythmic patterns could be elicited by stimulation using either bipolar or monopolar configurations. The effective cathode location was over the L2–L4 spinal cord segments in 34 cases, over L5–S2 in four cases, and caudal to S2 in one case. The mean relative stimulus intensity evoking the rhythmic patterns was 2.8 ± 1.1 times the posterior root-muscle reflex threshold of quadriceps and the mean SCS frequency was 29.5 ± 4.85 Hz, with a minimum of 22.5 Hz.

The burst-like activity had a constant phase relation between the muscles within each sample. The burst frequency was 0.71 ± 0.41 Hz, and varied between samples from 0.27 Hz to 1.84 Hz. An analysis of variance revealed significant interindividual differences of the burst frequency [F(6,32) = 6.551, P < 0.001]. Subject 5 had significantly higher frequencies than Subjects 1–4, whereas all other post hoc tests did not yield significant results. The burst frequency was not significantly correlated with the applied SCS frequency (r = 0.125, P = 0.450), but was predicted by the relative stimulation intensity [constant: B = 0.414, SE = 0.188; intensity: B = 0.105, SE = 0.062, β = 0.269, t(37) = 1.697, P = 0.098, R2 = 0.072], although not significantly, with a tendency of increased burst frequencies with increasing stimulation intensity.

Extension- and flexion-phase durations amounted to 1.29 ± 0.75 s and 0.59 ± 0.27 s, respectively, with relative extension-phase durations of 65.6 ± 10.4% of the full cycle durations. A linear regression showed that the cycle duration predicted the absolute duration of the extension phase [constant: B = −0.163, SE = 0.062; cycle duration: B = 0.775, SE = 0.030, β = 0.974, t(37) = 26.243, P < 0.001; R2 = 0.949] and flexion phase [constant: B = 0.163, SE = 0.062; cycle duration: B = 0.225, SE = 0.030, β = 0.782, t(37) = 7.633, P < 0.001; R2 = 0.612] as well as the relative extension-phase duration [constant: B = 0.527, SE = 0.030; cycle duration: B = 0.069, SE = 0.014; β = 0.624, t(37) = 4.855, P < 0.001; R2 = 0.389]. Thus, with longer cycle durations the absolute extension- and flexion-phase durations increased with different slopes (Fig. 3A), and the relative extension-phase duration increased while the relative flexion-phase duration decreased (Fig. 3B).

Figure 3

Relation between cycle duration and (A) absolute extension- and flexion-phase duration as well as (B) relative extension-phase duration of spinal cord stimulation-induced rhythmic activity. Solid lines are regression lines fitted to the observed data. Dashed lines indicate changes of cycle duration with equal extension- and flexion-phase durations.

A variety of rhythmic EMG patterns was intra- and interindividually generated. Each muscle normally displayed one major burst per cycle of rhythmic activity (e.g. Fig. 2A), with few examples of more complex patterns (e.g. Fig. 2D). Correlation of the EMG envelopes showed that, with respect to the tibialis anterior activity, the burst-like activities of the other ipsilateral muscles were largely coordinated either synchronously or alternatingly. Quadriceps were coactive with tibialis anterior in 69% of all 10-s segments of rhythmic activity, hamstrings with tibialis anterior in 36% and triceps surae with tibialis anterior in 38%. A logistic regression model showed that SCS frequency did not predict the phase-relation of the activity in either muscle group with the tibialis anterior bursts, nor between quadriceps and hamstrings bursts. On the other hand, stimulation intensity predicted the coactivation of triceps surae with tibialis anterior [constant: B = 3.424, SE = 1.250; intensity: B = −1.140, SE = 0.416, Cox and Snell R2 = 0.227, odds ratio = 0.320, χ2(1) = 10.019, P = 0.002], i.e. increasing the relative stimulation intensities increased the probability of generating alternating rhythmic activation of tibialis anterior and triceps surae.

All possible permutations of alternating or synchronous relations between rhythmic activities of quadriceps, hamstrings and triceps surae with respect to tibialis anterior were found, except for the pattern with quadriceps in reciprocity with the other muscles. The various patterns were not consistently generated across subjects (Table 2). Rhythmic coactivation of all unilateral muscle groups (Fig. 2A) was the most common pattern (n = 15). The remaining patterns occurred two to five times each. ‘Locomotor-like’ patterns with alternating activation of antagonists (Fig. 2B) were found four times. They were generated in Subjects 2, 3 and 5, i.e. those subjects with the most widespread activation of lower-limb muscles by a reinforcement manoeuvre. All locomotor-like patterns were generated with the cathode over the L2–L4 spinal cord segments, at frequencies of 25.5 Hz (n = 3) and 30.1 Hz, and relative stimulus intensities of 3.3 and 4.5 times the posterior root-muscle reflex threshold of quadriceps (n = 2, each).

Decomposition of the EMG patterns

NMF of the pooled data including all EMG profiles (n = 156) demonstrated that the variety of EMG patterns could be closely reproduced by linear combinations of a small number of basic activation patterns with muscle- and sample-specific weights (cf. Fig. 1F). According to the AICC and the Akaike weights AICw (Table 3), the model based on three basic activation patterns (k = 3) had the highest relative probability of all calculated models to best represent the process that had generated the data. The second best solution was the k = 4 model. All other solutions were judged as unlikely.

View this table:
Table 3

Goodness of fit measures of the NMF models

  • k = number of basic activation patterns of the non-negative matrix factorization (NMF) model; R2 = coefficient of determination; AICC = Akaike Information Criterion with correction; AICw = Akaike weights.

Nonetheless, the model with only two basic activation patterns explained 83% of the total variance across the EMG profiles of all rhythmic samples (R2 in Table 3). The two basic patterns were of sinusoidal-like shapes and had an alternating temporal structure, with the first basic pattern, p2,1, peaking at 60% of the extension phase, and the second one, p2,2, at 55% of the flexion phase, respectively (Fig. 4A, Model k = 2; pk,j is the j-th basic pattern of model k). The weights of these two basic activation patterns, required to reconstruct the various EMG profiles, had a highly significant negative correlation (Fig. 4B, Model k = 2).

Figure 4

Basic activation patterns pk,j and relations between their weights wk,j identified by non-negative matrix factorization from the pooled EMG-profile data across all subjects, muscles and samples of rhythmic EMG patterns. (A) Amplitude-normalized basic activation patterns shown for the models with k = 2, 3, and 4, ordered according to the relative timing of their peaks. (B) Negative (lines ending with small filled circles) and positive (lines ending with arrowheads) correlations between the weights of different basic activation patterns within each model indicate the tendency of reciprocal or simultaneous loading of two patterns. (C) Positive correlations in-between models, informing on preservation or splitting of basic activation patterns with increasing number of k. In B and C, only (Bonferroni corrected) correlations with P < 0.05 are illustrated. Index j in pk,j and wk,j enumerates the basic activation patterns based on their peak-timing; Ext = extension and Fl = flexion phase.

The k = 3 model identified a basic activation pattern that peaked at 56% of the extension phase, and a second and third basic pattern that peaked at 15% and 74% of the flexion phase, respectively (Fig. 4A, Model k = 3). The first and third basic patterns closely resembled p2,1and p2,2, and the highly significant negative correlation between their weights was retained (Fig. 4B, Model k = 3). The weights of the second basic pattern p3,2 that peaked during early flexion were not significantly correlated to those of either p3,1 or p3,3.

Comparison between the k = 3 and k = 2 models demonstrated a highly significant positive correlation between the weights of p3,1 and p3,3 with those of the respective basic extensor and flexor activation patterns of the lower-order model, respectively (Fig. 4C). Further, there was a near-even relationship between the weights of the early flexor component p3,2 with those of both, the extensor and flexor components of the k = 2 model. Hence, muscle activations explained by the additional early flexor component of the k = 3 model were otherwise approximated by the basic extensor or flexor component in the lower-order model.

The four-basic-pattern solution was characterized by a new component, p4,1, with its maximum at 15% of the extension phase (Fig. 4A, Model k = 4). The second basic pattern peaked at 61% of the extension phase, and the third and fourth basic activation patterns were almost identical to the two flexor components of the k = 3 model and peaked at 14% and 72% of the flexion phase, respectively. There was a highly significant positive correlation between the weights of the two extensor components and a highly significant negative correlation between the weights of these two patterns with the weights of the late flexor pattern p4,4 (Fig. 4B, Model k = 4). As in case of the k = 3 model, the weights of the early flexor pattern p4,3 were not significantly correlated to those of either of the other basic patterns. Correlations of weights between models suggested that the basic activation patterns of the lower-order models were retained in the k = 4 model, and the new pattern p4,1 was related to the basic extensor pattern of the k = 3 and k = 2 models (Fig. 4C).

Across the variety of rhythmic EMG patterns, the basic extensor pattern of the k = 2 model preferentially loaded on the hamstrings and triceps surae muscle groups and least on tibialis anterior, whereas there was no muscle-specificity of the basic flexor pattern apart from the bias to contribute to the tibialis anterior activity (Fig. 5). These relations of weights were retained across the k = 3 and k = 4 models. The early flexor pattern was more preferentially distributed to hamstrings than to quadriceps and triceps surae (Fig. 5, w3,2 and w4,3). The relative distributions to the muscles were similar for the first two basic patterns of the k = 4 model that peaked during the extension phase.

Figure 5

Correlations between weights wk,j of a given basic activation pattern across muscles and all selected EMG samples. Relations are given for the models with k = 2, 3, and 4 basic activation patterns for quadriceps (Q), hamstrings (Ham), tibialis anterior (TA), and triceps surae (TS). Bars indicate group mean and standard error, normalized to the individual maximum weight of the respective basic activation pattern and model. Index j in wk,j enumerates the basic activation patterns based on their peak-timing. All group tests P < 0.001, asterisks indicate significant post hoc tests: *P < 0.05, **P < 0.01, ***P < 0.001.


We studied rhythmic EMG patterns produced by the functionally isolated human lumbar spinal cord under repetitive, constant electrical stimulation. They included a variety of coactivation, mixed synergy and locomotor-like patterns with a range of burst frequencies. Statistical decomposition of the EMG profiles across subjects, muscles, and the various rhythmic patterns revealed a common, flexible functional organization of their underlying spinal interneuronal circuitry.

Central and peripheral influences upon the lumbar spinal cord circuitry

Some subclinical residual excitatory and/or inhibitory translesional influence on the excitability of the lumbosacral spinal cord circuitry was detected in all subjects studied here (Table 2). This observation suggests that some descending long-tract fibres or propriospinal connections transiting the lesion may have survived their spinal cord injury clinically classified as motor complete, as shown possible in earlier studies neuropathologically (Kakulas, 2004) and neurophysiologically (Sherwood et al., 1992; McKay et al., 2004). Such residual influence was not a decisive component for the generation of the rhythmic patterns in response to SCS. The absence of consistent rhythmic EMG patterns in Subjects 8–10 might have been due to individual differences in the state of the spinal networks (cf. Harkema, 2008). Nevertheless, the three individuals in whom locomotor-like EMG patterns were generated (Subjects 2, 3 and 5) most probably had the highest level of base excitatory supraspinal support of the lumbar circuitry (Table 2, reinforcement).

The lumbosacral neural circuits received sensory input from the legs, yet adapted to the supine position, where hip-joint extension is limited and axial leg loads are absent. Thus, sensory information otherwise signalling phase transitions in the locomotor cycle and entraining the step frequency both in the chronic spinal cat (Grillner and Rossignol, 1978; Hultborn and Nielsen, 2007) and in people with motor-complete spinal cord injury (Dobkin et al., 1995; Harkema et al., 1997; Dietz et al., 2002) were not provided. Additionally, antidromic volleys generated by SCS block or mitigate peripheral feedback. Thus, we obtained a model to study central spinal locomotor mechanisms, i.e. a model approximating immobilized, ‘fictive’ locomotion (cf. Pearson and Gordon, 2013) with the constraints inherent to human studies. An additional yet conceptually different step would be to investigate rhythmogenesis and its interaction with sensory feedback by allowing movement to occur with minimal mechanical constraints (Gurfinkel et al., 1998; Selionov et al., 2009). Physiologically appropriate phasic afferent feedback could have coordinated our motor outputs into more locomotor-like patterns (cf. Fedirchuk et al., 1998).

Neural processes initiated by spinal cord stimulation

In individuals with motor-complete spinal cord injury, SCS-effects upon lower-limb muscle activation are initiated by bilateral stimulation of afferent fibres within posterior rootlets/roots of several lumbar spinal cord segments (Rattay et al., 2000; Minassian et al., 2007; Angeli et al., 2014) and the subsequent trans-synaptic activation of motor neurons and spinal circuitry (Minassian et al., 2004; Capogrosso et al., 2013). The rhythmic EMG patterns reported here were strictly composed of series of stimulus-triggered, modulating responses. We have previously proposed that repetitive SCS activates spinal interneuronal circuits that in turn modulate the concomitantly elicited posterior root-muscle reflexes (Minassian et al., 2004, 2007; cf. Degtyarenko et al., 1998; Kinoshita and Yamaguchi, 2001).

Effects of spinal cord stimulation parameters

Eighty-seven per cent of the rhythmic EMG segments were generated with the cathode over the L2–L4 spinal cord segments. The role of these segments in the generation of rhythmic motor patterns was previously recognized (Dimitrijevic et al., 1998; Shapkova, 2004). Either the L2–L4 posterior-root fibres have specifically effective projections onto the circuitry that generated the rhythmic outputs, or key elements of these functional neural networks are located within these segments (Grillner, 2006). Moreover, 59% of the rhythmic EMG patterns described here were only expressed unilaterally, supporting earlier indications that side-specific rhythm-generating networks exist bilaterally in the human spinal cord (Yang et al., 2005).

Increasing stimulus intensities significantly increased the odds of changing the EMG pattern from synchronous bursts of activity across all ipsilateral muscles to any other type of pattern. The elicitation of synchronous activity in antagonists may suggest that elements of the spinal locomotor network were not evenly activated (Fedirchuk et al., 1998). This view is also supported by the observation that fictive locomotion often starts with synchronous activity in flexor and extensor nerves and then with time, the interneuronal network organizes itself to generate flexor-extensor alternation (Meehan et al., 2012).

The minimum stimulation frequency needed to elicit rhythmic EMG patterns was 22.5 Hz, suggesting that certain repetition rates of inputs and subsequent temporal and spatial synaptic summation processes were required to activate the rhythm-generating networks. Yet, the stimulation frequency alone neither determined the burst frequency nor the specific phase coordination of the generated motor outputs to the different muscles. A similar effective SCS-frequency range of 25–50 Hz was reported previously (Dimitrijevic et al., 1998; Gerasimenko et al., 2002). Such frequencies also facilitate rhythmic activity generated by step-related feedback during assisted treadmill stepping (Minassian et al., 2005; Harkema et al., 2011) and may enable synergistic lower limb movements initiated by residual supraspinal influence (Angeli et al., 2014). Shapkova (2004) found two ranges of stimulation frequencies, 3–12 Hz and 20–30 Hz, to be effective in inducing fast air-stepping in paraplegic children lying supine with their legs elastically suspended in a semi-flexed position. The latter is close to the frequency range reported here and probably induced similar rhythm generating mechanisms. The lower frequencies would have evoked strong twitch-contractions (cf. Murg et al., 2000; Jilge et al., 2004; Sayenko et al., 2014) that together with the suspension may have resulted in the fast oscillatory movements described by Shapkova and Schomburg (2001).

Rhythm generation

The burst frequency of the SCS-induced output patterns appeared centrally determined as no timing information through the stimulation was provided. The observed interindividual differences suggest that the operation of the rhythm generator also depends on its central state of excitability. The wide range of burst frequencies of 0.27–1.84 Hz covered different rhythmic behaviours. The lower range was similar to the burst frequencies of spinal myoclonus of 0.3–0.6 Hz (Bussel et al., 1988; Calancie, 2006). The mean burst frequency of 0.71 ± 0.41 Hz was close to that of slow gait of adults (0.77 Hz, 1.54 step/s, 2.6 km/h) and the upper range of frequencies included fast gait (1.31 Hz, 2.61 steps/s, 6.0 km/h; Oberg et al., 1993).

The longer mean extension- than flexion-phase durations and the adaptation to decreasing cycle durations by decreasing relative extension phase durations suggest that the ‘asymmetric alternating activity’ (Grillner, 2006) characteristic of terrestrial locomotion, including human gait (Murray, 1967), may be at least partially controlled at the spinal level in humans.

Pattern generation

In studies of the modular organization of locomotor control, a motor module is considered a functional unit of the CNS producing specific motor outputs and, computationally seen, consists of a basic activation pattern along with its variable weights of distribution to different muscles (Clark et al., 2010; Dominici et al., 2011; Fox et al., 2013). It is assumed that motor modules are largely controlled by brainstem and spinal networks (Fox et al., 2013; Giszter and Hart, 2013). Yet, in studies of volitional locomotor tasks, cortical modulation (Dominici et al., 2011) and diverse locomotor behaviours (Zelik et al., 2014) contribute to the underlying motor modules, with further influence from step-related sensory feedback and biomechanical events of the gait cycle (Ivanenko et al., 2003). Here, the absence of such influences allows the basic activation patterns to be interpreted as intrinsic spinal drives (Dominici et al., 2011; Giszter and Hart, 2013), generated by spinal networks in response to the ongoing epidural stimulation.

The sinusoidal-like basic activation patterns identified by the two-pattern NMF model represent alternating bursts of activity. When reciprocally distributed to muscles, they describe patterns similar to alternating activity in extensor and flexor nerves recorded during fictive locomotion in experimentally reduced animal models (Hägglund et al., 2013). This similarity may suggest that the two basic activation patterns were generated by analogous spinal circuits as in other mammals (Grillner, 2011) and relate to the drives of extensor- and flexor-burst generators. Yet, each of the basic patterns contributed to the activation of several muscles in agonistic relationships as well. The two elementary extensor and flexor patterns were retained and further individualized across the NMF models with three and four basic activation patterns. The additional basic activation pattern that peaked in the early flexion phase explained muscle activations that, in the two-pattern model, were shaped by either extensor or flexor drives. Yet, within the higher-order models, it loaded independently from the basic extensor and flexor patterns and was characteristically associated with hamstrings activation. Together, these results are reminiscent of the complex central activation patterns of posterior biceps and semitendinosus muscles during fictive locomotion in reduced animal preparations (McCrea and Rybak, 2008; Guertin, 2009). The additional activation pattern identified by the four-pattern model resulted from the basic extensor pattern of the lower order models and mainly added details to the shapes of rhythmic EMG activities peaking during the extension-like phases.

Flexible organization of burst-generating elements

The basic patterns identified here resemble those underlying volitional locomotion in toddlers as well as adult rhesus monkeys, cats, rats, and guinea fowls (Dominici et al., 2011). In volitional locomotion, the distribution of the basic activation patterns to the muscles is invariantly specified by the motor task. Here, all bifunctional muscle groups could be driven by both, extensor- and flexor-related basic patterns—simultaneously or reciprocally—with their respective contributions varying between samples. This implies that the burst-generating elements driving the motor neuron pools are modularly organized and are flexibly combined to produce a variety of rhythmic patterns. Similarly, various types of fictive motor patterns were observed in the acutely spinalized marmoset monkey including alternating as well as synchronous activity in flexor and extensor nerves (Fedirchuk et al., 1998). These activities were interpreted as ‘component fictive patterns’, which if combined would resemble a true locomotor pattern.

In conclusion, the functionally isolated human lumbosacral spinal cord can organize a rich repertoire of rhythmic output patterns to leg muscles in response to repetitive, constant stimulation. This range of rhythmic patterns can be explained by variable combinations of a small number of common spinal neural drives, each with a distinct temporal structure, and may be realized by multiple burst generators with modifiable connection strengths, determining reciprocities and phase lags of the rhythmic drives (Grillner, 2006; Hägglund et al., 2013). To use this remarkable capability and flexibility of the spared, sublesional lumbosacral spinal networks, neurorehabilitation approaches should consequently target the configuration of the burst-generating elements in order to incorporate specific adaptations to gait and stance.


Vienna Science and Technology Fund (WWTF), grant number LS11-057, Wings for Life Spinal Cord Research Foundation, grant number WFL-AT-007/11, and Austrian Science Fund (FWF), grant number L512-N13.


We would like to express our appreciation to Milan R. Dimitrijevic and W. Barry McKay for critical reading of the manuscript and Renate Preinfalk and Maria Auer for data collection.

non-negative matrix factorization
spinal cord stimulation


View Abstract