Biological oscillations for learning walking coordination: dynamic recurrent neural network functionally models physiological central pattern generator

Introduction

Neuronal modules profoundly influence many aspects of motor behavior. However, little is currently known about the control mechanisms that allow the coordination of these modular entities. In this theoretical context, human locomotion is a challenging movement because of the numerous neuroanatomical modules implicated in the different aspects of whole body movement, ranging from the cyclic propulsion of the limb to balance control. In spite of these intricate movement components and neuronal modules involved, it has been proposed that the overall control can be realized by reducing the number of degrees of freedom of the system into global variables (Bernstein, 1967; Lacquaniti et al., 1999, 2002; Flash and Hochner, 2005; Latash, 2008). The fact that the elevation angles of the three main lower limb segments are coordinated during gait within a covariation plane (Borghese et al., 1996), forming an elliptic loop corroborates the idea that control of locomotion is also submitted to the general law of reducing variables (Barliya et al., 2009). This general law is also applicable for different walking speeds (Bianchi et al., 1998), for forward and backward directions (Grasso et al., 1998), rectilinear or curvilinear walking paths (Courtine and Schieppati, 2004), walking with stilts (Dominici et al., 2009; Leurs et al., 2011), or with a transfemoral prosthesis walk (Leurs et al., 2012), with different levels of body weight unloading (Ivanenko et al., 2002) and for running (Ivanenko et al., 2007). Notably, this inter-segmental coordination is not present in newly walking toddlers (Cheron et al., 2001a,b; Ivanenko et al., 2005), but covariation planarity rapidly emerges over the first few days of independent walking experience, while the shape of the loop and plane orientation “mature” gradually over several years. This evolution indicates that the attractor plane and the shape of the loop are neurophysiologically defined, rather than being imposed by biomechanical constraints (see Hicheur et al., 2006; and Ivanenko et al., 2008 for discussion). More recently, the developmental study of this complex behavior in a new born baby (Dominici et al., 2011) has permitted revisiting the concept of locomotor modules coding for the control of movement primitives.

This modular approach raises the question of the dynamic coordination of modules in the context of oscillatory properties of neuronal ensembles. Indeed, the dynamical structure of these modules must logically obey a common principle for movement generation: the production of oscillatory activity. Although this principle is well accepted in case of the rhythmic nature of locomotion (Georgopoulos and Grillner, 1989; Grillner, 2006), the recent study of Churchland et al. (2012) surprisingly demonstrates that non-periodic movements such as reaching are also generated by neuronal oscillation. This means that there is a strong possibility that the spinal modules organized in a central pattern generator (CPG) could be dynamically controlled by cortical and/or supraspinal oscillations. Interestingly, Barliya et al. (2009) recently modeled the time course of elevation angles of the three lower limb segments in terms of simple oscillators coupled with appropriate time shifts reproducing the orientation of the plane and their elliptical shape. The oscillators were obtained by taking, after Fourier transform, the first three harmonics of the elevation angle kinematics. Each of these oscillators could be interpreted in term of oscillatory module acting in such a way that the orientation of the plane and the elliptic shape of the coordination are conserved. It could thus be possible that oscillatory signals derived from these harmonics are used for activating CPG modules.

The existence of CPG in the spinal cord in mammals has been proposed a century ago (Brown, 1911). In essence, it represents a group of neurons acting as a network to produce coordinated patterns of rhythmic activity. New evidence has shown the presence of CPG in the spinal cord in humans (Calancie et al., 1994; Bussel et al., 1996). The characteristics of such CPG modules are their adaptability and robustness that lead to the production of different gait patterns adapted to their current environmental context. For example, young infants (less than 1-year old) are already able to walk over a split-belt treadmill with different types of coupling (Yang et al., 2005). Some of them were even able to walk in opposite directions. Mimicking physiology, the robotic and neuroscientific community developed artificial CPGs that are commonly used to animate robots from multi-legged insect-like to humanoids (for a review see Ijspeert, 2008). In their pioneering work in the cat, Patla et al. (1985) proposed an analytic model of limb locomotor pattern generator based on recorded muscle activity induced by electrical stimulation over the mesencephalic locomotor region (MLR) in the decerebrate cat. In this model, locomotor like patterns of six muscles resulted from six independent oscillators with dedicated parameters. These oscillators were reduced to simple sine and cosine functions fed by a tonic input. Since then, different methods have been used from coupled non-linear oscillators (Ijspeert et al., 2007; Duvinage et al., 2011, 2012a), to highly detailed biophysics of small groups of interconnected neurons (Hellgren et al., 1992) and rhythm genesis of larger groups of neurons (Wadden et al., 1997) mainly in animal models. While human locomotion has often been reproduced computationally in the robotic field using equations that are numerically integrated (Righetti et al., 2005; Righetti and Ijspeert, 2006; Ceccato et al., 2009; Duvinage et al., 2011, 2012a), few methods involving neuron modeling for human gait generation have been studied so far. Among them Prentice and coauthors (1998, 2001) have successfully transformed fundamental timing signals (sine and cosine inputs) into individual muscles activation bursts related to gait locomotion at different speeds using a feedforward neural network. Our group used the electromyographic (EMG) signals of the lower limb muscles as input for a dynamic recurrent neural network (DRNN) producing as output the lower kinematics during locomotion (Cheron et al., 2003, 2012). However, the possibility to produce the motion of the lower limb segments by means of oscillations derived from the three harmonics of the Fourier transform of walking kinematics has not yet been assessed by means of the same DRNN. We show here that after learning based on different walking velocities, the DRNN is able to reproduce the lower limb kinematics of both legs. The DRNN can also generalize to the unlearned walking velocities. The analysis of the required network structure (e.g., number of units, distribution of time constant, and synaptic sign) provides a basis for the discussion about the constraints required for the elaboration of a CPG.

Methods

Experimental Setup

Seven healthy subjects aged from 25 to 35 years (mean age: 28 years) participated in this experiment. The protocol consisted of walking on a treadmill at 11 different velocities [from 1 (0.28 m/s) to 6 km/h (1.67 m/s) stepped by 0.5 km/h (0.14 m/s)] leading to 11 trials per subject (total of 77 trials over all subjects). Whole body kinematics was recorded using Vicon system with six infrared Bonita cameras recording at 100 Hz during 40 s for each trial. The tracking consisted of recording 29 markers placed over the whole body. The position of the markers in an orthogonal reference was computed using a custom biomechanical model. From the position of the markers of both legs, the kinematic (elevation) angles relative to the vertical plane of the laboratory have been calculated bilaterally for thighs, shanks, and feet. In this study we performed two experiments. The first one, “proof of concept,” was done to ascertain the feasibility of a DRNN to learn elevation angles from pure sine waves. This part only includes supervised learning on a single pattern (Figure 1A). The second experiment was performed in order to study the possibility to learn multiple patterns of walking (i.e., velocities) and to predict kinematics from unlearned patterns (Figure 1B). After that, the DRNN structures were analyzed.

FIGURE 1

Figure 1. Realization of a CPG module based on experimental recordings of human walking. (A) In the first step, the elevation angles of the thigh, shank, foot of the two legs are transformed in the time frequency domain using fast Fourier transform (FFT). Then the signal characteristics are computed for the first three harmonics and back-transformed into temporal space using sine waves formulation. These data are fed as input to the DRNN learning the transformation in kinematic signals (elevation angles). (B) A successful learning will permit the DRNN to predict kinematics based on sine wave inputs only. These output signals can produce biologically plausible walking patterns in a virtual reality avatar or an actual robotic exoskeleton.

Dynamical Recurrent Neuronal Network

The DRNN is capable of modeling time-varying input–outputs and has varying weights as well as varying time constants for the artificial neurons (Pearlmutter, 1990). The adaptive time constants make it dynamic (Draye et al., 1996). The DRNN is governed by the following equations:

$\begin{array}{cc}{T}_i\frac{d{y}_i}{dt}=-y_i+F(x_i)+I_i& (1)\end{array}$

where F(a) is the squashing function F(a) = (1 + e^−a)⁻¹, y_i is the state or activation level of unit i, I_i, is an external input (or bias), and x_i is given by:

$\begin{array}{cc}{x}_i={\displaystyle \sum _j{w}_{ij}{y}_i}& (2)\end{array}$

which is the propagation equation of the network (x_i is called the total or effective input of the neuron i, w_ij is the synaptic weight between units i and j). The time constants T_i will act like a relaxation process. In order to make the temporal behavior of the network explicit, an error function is defined as:

$\begin{array}{cc}E={\displaystyle \underset{t_0}{\overset{t_1}{\int }}q(y(t),t)dt}& (3)\end{array}$

where t₀ and t₁ give the time interval during which the correction process occurs. The function q(y(t), t) is the cost function at time t which depends on the vector of the neuron activations y and on time. We then introduce new variables p_i (called adjoint variables) that will be determined by the following system of differential equations:

$\begin{array}{cc}\frac{d{p}_i}{dt}=\frac{1}{T_i}{p}_i-e_i-{\displaystyle \sum _j\frac{1}{T_i}}{w}_{ij}{F}^{\prime }(x_j)p_j& (4)\end{array}$

with boundary conditions p_i(t₁) = 0. After the introduction of these new variables, we can derive the learning equations:

$\begin{array}{l}\begin{array}{cc}\frac{\delta E}{\delta w_{ij}}=\frac{1}{T_i}\underset{t_0}{\overset{t_1}{{\displaystyle \int }}}{y}_i{F}^{\prime }(x_j)p_{j}dt& (5)\end{array}\hfill \\ \begin{array}{cc}\frac{\delta E}{\delta T_i}=\frac{1}{T_i}\underset{t_0}{\overset{t_1}{{\displaystyle \int }}}{p}_i\frac{d{y}_i}{dt}& (6)\end{array}\hfill \end{array}$

The sinusoid signals derived from the Fast Fourier transform (FFT) kinematic data are sent as input to a DRNN (Figures 1, 2 cf. Experiments 1, 2) that has to learn from these data to produce the kinematics specified as elevation angles (Figure 1A). Successful trainings were also used to produce kinematic patterns from unknown inputs (Figure 1B) aiming to produce walking for multiple purposes, such as virtual avatars or robotic exoskeletons. The training is supervised, involving learning rule adaptations of synaptic weights and time constant of each unit (Draye et al., 1995, 1996). A specific training procedure using Almeida algorithm was used to optimize learning performance (Cheron et al., 2011). The DRNN involves a looping mechanism (fully connected structure) which enables this network to learn and store information (memory). This equips the network with the ability to model complex situations with multiple influences. The DRNN was successfully used to map EMG signals into kinematics during walking (Cheron et al., 2003), for the identification of the triphasic EMG pattern in subjects performing fast elbow flexion movements (Cheron et al., 2007) or to predict specific muscular activity in elite fencers compare to amateurs (Cheron et al., 2011).

FIGURE 2

Figure 2. (A) Elevation angles (in degree) of the three segments (thigh, shank, and foot) for each leg as a function of time for a subject walking at 3 km/h (0.83 m/s). The whole pattern is presented for duration of 5 s. The black lines represent a pattern for two gait cycles used to determine the FFT characteristics. (B) The mean FFT transformation for six joints for 40 s (in gray) and the two gait cycles (in black). Note that the two gait cycle patterns are selected so as to preserve the FFT characteristics in terms of amplitude and frequency (stars) for the overall pattern of 40 s. These characteristics are used as parameters to generate sine waves as input of the DRNN.

Experiment 1: Proof of Concept

DRNN computation has been used to prove that simple sine waves with specified temporal characteristics can be used as input to an artificial neural network to be transformed into elevation angles obtained from kinematic recordings during locomotion (Figure 3).

FIGURE 3

Figure 3. DRNN learning method for Experiment 1. For each subject, we trained 100 DRNNs to learn a pattern of kinematic [corresponding to a velocity of 3 km/h (3.83 m/s)] from sine waves (Equations 7 and 8). Three sets of learning were defined as the input differed (SEA, SEB, SEC). The structures of DRNNs were modeled with 30 hidden units for each set of training. Each network iterated 10,000 times to change its synaptic weights and time constants. For each subject and structure we selected the network with the highest SI value. For each condition, the design of the network is then processed with 6 inputs, 30 hidden units, and 6 outputs.

Input: construction of sine waves

As the learning phase of the DRNN is a time-consuming process, we had to select appropriate sample from the whole available data as input. Moreover, even if the kinematics during walking is relatively stable, it may vary too much to feed the DRNN during the learning phase. For these reasons, we extracted two consecutive gait cycles from the 40 s of experimental data recorded in each trial (Figure 2A, black curves). They were chosen so as to be representative in terms of frequency of the whole recording set. Then, to determine the kinematic characteristics of gait, we transformed the data of the lower limb segments elevation angles into the time frequency domain using the Matlab fft function to get the FFT power amplitude and their related frequency values of the first three harmonic peaks (Figure 2B). It has been shown previously that the first two Fourier harmonics accounted for approximately 98% of the experimental variance of the thigh, shank, and foot angles (Bianchi et al., 1998). We decided to create sine waves based on the characteristics of the first three harmonics to ensure that the signal proposed as learning input to the DRNN contains enough information to match the desired output precisely.

We extracted the values of the three frequencies (f1, f2, f3) corresponding to the maximal amplitudes (a1, a2, a3). It was verified that f2 = 2f1 and f3 = 3f1. We then artificially produced six sinusoidal signals using frequency values as parameters (Equations 7 and 8).

$\begin{array}{l}\begin{array}{ll}{y}_{fi,1}=\sin (2\pi \times fi\times t)\hfill & (7)\hfill \\ y_{fi,2}=\sin (-(2\pi \times fi\times t))\hfill & (8)\hfill \end{array}\\ \text{ For}i=1,2,3\hfill \end{array}$

These six sine waves thus correspond to the frequency characteristics of the kinematic patterns that will be further used as pattern to be learned. For the sake of clarity we called the set of original inputs set SEA (for set of Equations A). Additionally, we produced six new sine waves using the preceding computations (Equations 7 and 8) where f1′ = f1 + 0.25 Hz and six other sine waves where f1″ = f1 − 0.25 Hz, respectively called SEB (for set of Equations B) and SEC (for set of Equations C). Please note that in the latter two cases the original relations f2 = 2f1 and f3 = 3f1 were respected and hence f2′ = 2f1 + 0.50 Hz; f3′ = 3f1 + 0.75 Hz in the set SEB and f2″ = 2f1 − 0.50 Hz; f3″ = 3f1 − 0.75 Hz in the set SEC. Three different input sets (SEA, SEB, SEC) were thus defined for learning.

Pattern to learn: kinematic data

The pattern to learned corresponds to the elevation angles of the right and left thigh, shank, foot segments in the two gait cycles of a 3 km/h (0.83 m/s) walk, normalized between −1 and 1. The outputs were the same for SEA, SEB, and SEC whereas inputs differed.

DRNN learning phase

It was expected that the DRNN would learn how to transform the input to output thanks to its dynamical and recurrent structure of 30 hidden neurons. For each of the seven subjects, the network iterated 10,000 times. This process was repeated 100 times, leading to 100 different networks. At the end of the learning phase, we selected for each subject the network for which the difference between predicted and real kinematics was minimal. This computation was performed by calculating a similarity index (SI) (Bengoetxea et al., 2010) defined by the following equation:

$\begin{array}{cc}SI=\frac{{\displaystyle {\int }_{-\frac{Tc}{2}}^{\frac{Tc}{2}}{p}_1}(t)p_2(t)dt}{{\left[{\displaystyle {\int }_{-\frac{Tc}{2}}^{\frac{Tc}{2}}{p}_1}{(t)}^{2}dt{\displaystyle {\int }_{-\frac{Tc}{2}}^{\frac{Tc}{2}}{p}_2}{(t)}^{2}dt\right]}^{\frac{1}{2}}}& (9)\end{array}$

where Tc is the period of the limit cycle, p₁ and p₂ being synchronized, i.e., the matching between both patterns is based on the maximum of each pattern. Note that if both functions are identical, SI = 1. SI was calculated with the recorded pattern of elevation angles and the output of the DRNN.

Experiment 2: Multiple Pattern Learning and Prediction

In this experiment sine waves were modulated in frequency and amplitude and transformed into kinematic data using multi-pattern training. The prediction of kinematic pattern from unknown sine wave pattern was also tested.

Input

As for the proof of concept methods, we extracted two gait cycles of each trial for multiple velocities (in km/h) (v = {1.5, 2.5, 3.5, 4.5, 5.5)}. We transformed the kinematic data into the time-frequency domain to get their frequency (f1, f2, f3) and amplitude (a1, a2, a3) (Figure 2B) parameters using the following set of Equations (10 and 11).

$\begin{array}{l}\begin{array}{ll}{y}_v,f{i}_v,a{i}_{v,1}=a{i}_v\times (sin(2\pi \times f{i}_v\times t))\hfill & (10)\hfill \\ y_v,f{i}_v,a{i}_{v,2}=a{i}_v\times (sin(-(2\pi \times f{i}_v\times t)))\hfill & (11)\hfill \end{array}\\ \text{ For}i=1,2,3\hfill \end{array}$

Patterns to learn: kinematic data

The patterns to be learned consisted of elevation angles of the right and left thigh, shank, and foot segments corresponding to the two gait cycles, normalized between –1 and 1. These patterns were obtained from recordings at multiple velocities (in km/h) (v = {1.5, 2.5, 3.5, 4.5, 5.5}).

DRNN training phase

We used the possibility to train the DRNN in a multi-pattern purpose (Bengoetxea et al., 2005). For each subject, the DRNN was trained to match the inputs/outputs patterns corresponding to five different velocities within a single DRNN structure. Hundred networks were made per subject, each of them iterating 50,000 times. This operation was assigned for 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, and 80 hidden units. When this was completed, we selected the best network for each subject and each number of hidden units using the maximal SI values.

Output prediction

We built sine waves from intermediary pattern velocities (km/h) (v = {2, 3, 4, 5}) as explained above. We fed the best DRNN structures previously obtained after the training phases with these unlearned inputs sine waves and analyzed DRNN behavior by calculating the predicted output with experimental data using SI values.

Neuronal Properties and Connectivity after Learning

Introduction of timing allows modeling of more complex frequency behavior, improves the non-linearity effect of the sigmoid function and the memory effect of time delays (Draye et al., 1995). The distribution of the time constant and the synaptic weights between units (Draye et al., 1996) after learning was analyzed after multiple pattern learning and prediction. In particular, we recorded the number of positive and negative connections inside DRNN structures of the best networks. Additionally we studied the distribution of neuronal time-constants.

Statistical Analysis

Statistical analysis was performed using Statistica software (Statsoft, www.statsoft.com). For each test performed and described in the result section, we firstly verify whether the data were distributed normally using Kolmogorov–Smirnov test. When the data were normally distributed we use ANOVA with repeated measures and post-hoc Fisher analyses. When it was not possible to use ANOVA we used non-parametric tests such as Friedman ANOVA or sign tests.

Computations

All DRNN computations were performed in the Hydra computing center hosted in ULB (https://cc.ulb.ac.be/hydra/). We allocated 1 node and 10 Gb of memory per computation (i.e., per subject per condition in the Experiment 1, per subject per structure in the Experiment 2). The computations were run in parallel to optimize the learning duration. Afterwards we simulated 5% of the overall experiment in the same conditions to estimate the learning time. The overall duration of the process was obtained by linear interpolation (Figure 4).

FIGURE 4

Figure 4. Illustration of the DRNN structure and computational time. (A) Computational learning duration for one subject for experiment 2 (in hours). The gray circles represent the estimation of time required per number of hidden units from a recalculation of 5% of the procedure. Please note that the overall process for the seven subjects has taken approximately 200 days to be computed. (B), (C), and (D) represent the connectivity of the hidden layer with 5, 10, and 20 neurons, respectively. As the number of neurons increases, the number of connections increases by a factorial multiplication as well as the time required to adjust connection weight and time constant. Please note that only the units of the hidden layers are represented without input or output neurons.

Results

Experiment 1: Proof of Concept

A statistical test was designed to compare SI values from different input types (SEA, SEB, or SEC) (Figures 3,5). The Kolmogorov–Smirnov did not reject the hypothesis that SI values were normally distributed (D = 1.4414, p > 0.20) when analyzing together values of the different inputs (SEA, SEB, or SEC). We used an ANOVA with repeated measures where dependant variable was SI and the independent variable chosen was the type of input. The analysis showed an effect of the input frequency in the DRNN prediction (SI value) [H F_{(2, 12)} = 38.110, p = 0.00001]. Post-hoc analysis confirmed that SI values of original group of unchanged frequency input (SEA) were higher than the 2 modified groups where frequency inputs have changed (SEB and SEC).

FIGURE 5

Figure 5. Planar covariation between normalized thigh, shank, and foot for all participants. (A) Recording data—the patterns correspond to the real kinematics of 2 gait cycles for all participants. Best kinematic patterns predicted by the DRNN for each subject where SEA (B), SEB (C), and SEC (D) were fed as input. $\overline{SI}$ is the average of similarity index value (SI) in each condition for all subjects when compared with the recording data.

Experiment 2: Prediction of Intermediary Velocities

As we decided to use frequency (f1, f2, f3) and amplitude (a1, a2, a3) characteristics to modulate inputs for multiple learning procedures (Equations 10 and 11), we have verified that there was a statistical significance of these parameters for different velocities. The Kolmogorov–Smirnov test for f1 (D = 0.14370, p < 0.1), f2 (D = 0.14370, p < 0.1), and f3 (D = 0.14370, p < 0.1) was not clear enough to reject the fact their population may follow a normal law. When looking at the distribution, they tend to be normal and it is possible that the significance of the test is due to the weak number of values. According to similar test, the values for parameters a1 (D = 0.08476, p > 0.2), a2 (D = 0.08758, p > 0.2), and a3 (D = 0.09035, p > 0.2) were normally distributed. We then used two ANOVA with repeated measures where the dependent variables were, respectively, the amplitude and frequency values and the within-subject factors were velocity of walking and the specific harmonic (1, 2, or 3).

ANOVA shows significant changes in amplitude [F_{(10, 60)} = 31.351, p < 0.0001] and frequency [F_{(10, 60)} = 69.276, p < 0.0001] with an increase in velocity. Post-hoc analyses revealed an increase in f1, f2, a1, a2, a3 and a decrease of f3 with an increase in the velocity. These significant differences justified their use for building specific sine waves for different walking velocities (Figure 6).

FIGURE 6

Figure 6. Evolution of the frequency (A) and amplitude (B) of the first three harmonics when transforming elevation angles into the time frequency domain (FFT) at different walking velocities [from 1 to 6 km/h (0.28 to 1.67 m/s)]. The FFT peaks of each angle are averaged for each velocity for both frequency and amplitude. Squares, circles, and triangles represent the mean value for 7 subjects and whiskers represent the 95% confidence interval.

Concerning the performance of the DRNN outputs, we wanted to verify if the SI value applied for the best networks was different for learning and prediction. Additionally we wanted to statistically check if the number of hidden units of the networks increases the SI value (Figure 7). The Kolmogorov–Smirnov did not verify that the populations of SI values among the learning (D = 0.18467, p < 0.01) or the prediction (D = 0.20684, p < 0.01) were normal. Thus, to compare the SI values between learning and prediction process we chose to use a sign test as the structure of the network (weights and time constant) were the same. This test reveals no significant differences in SI values between the two populations except when the network contains 1, 2, 4, 5, 6, and 7 neurons. Moreover we use a Friedman test to analyze SI values (dependent variable) according to the number of hidden units (independent variables) of the network. There is an effect of the number of hidden units to the SI both for learning [χ²₍₁₆₎ = 111.5630, p < 0.0001] and for prediction [χ²₍₁₆₎ = 109.3721, p < 0.0001].

FIGURE 7

Figure 7. SI values of all multi-patterns learning and prediction for all participants with different number of hidden neurons. Whiskers correspond to 95% confidence interval.

We also analyze the output of the DRNN for specific velocities where the hidden layer was the biggest with 80 units. A Kolmogorov–Smirnov run together with learning and prediction values did not verify that the populations of SI was normal (D = 0.18342, p < 0.05).

We use a Friedman test to analyze SI values (dependent variable) according to the velocity (independent variables) of the network. There is an effect of the velocity to the SI [χ²₍₈₎ = 30.22, p = 0.00019] (Figure 8). A post-hoc analysis at 0.05 significance level reveals that SI values of 4.5 km/h were different from 1.5 and 2 km/h. SI values of 4 km/h were also different from 2 km/h. An example of prediction of intermediate velocities in one subject is illustrated in (Figure 9).

FIGURE 8

Figure 8. SI values for learning and prediction with a 80 hidden units network according to walking velocity. Whiskers correspond to 95% confidence interval.

FIGURE 9

Figure 9. Prediction of intermediate velocities in one subject. (A) Planar covariation of the normalized thigh, shank, and foot for one participant walking at 2 km/h. (B) Planar covariation of the normalized thigh, shank, and foot for the same participant walking at 5 km/h. (C) Kinematics of walking of the normalized thigh, shank, and foot for the same participant walking at 2 km/h. (D) Kinematics of walking of the normalized thigh, shank, and foot for the same participant walking at 5 km/h. In each subplot the red dotted line corresponds to the experimental data whereas black thick line corresponds to the output of the DRNN.

DRNN Structures from Experiment 2

Weight distribution analysis

The percentage of positive and negative weights was calculated for each best network per subject per condition. We wanted to verify if the distribution of the positive and negative weights were different (Figure 10A). The Kolmogorov–Smirnov test did not verify that the distribution of positive weights were normal (D = 0.15022, p < 0.01) nor the distribution of negative weights (D = 0.15022, p < 0.01). We then used a non-parametric sign-test to compare the two-distribution as the two samples were dependent. Regardless of the number of hidden units, the test shows a difference between the distributions of the two populations (Z = 9.482, p < 0.0001). When the number of units in the test were included, it appears that the population of negative weights is higher than population of positive weights when the network possesses more than 3 hidden units (for 4 hidden units, Z = 2.268, p = 0.023342). When the structure of the DRNN reaches 80 hidden neurons, 70.6 ± 0.84% of the synapses are negative.

FIGURE 10

Figure 10. Distribution of weights (A) and time constants (B) among the DRRN structures. The horizontal axis represents the number of hidden neurons and vertical the percentage of population. The bars represent the mean distribution among the best DRNN for each subject. Whiskers represent the 95% confidence intervals.

Time constant distribution analysis

The distribution of the time constant was represented by the median of neuronal time constants of the best networks per subject and number of hidden units (Figure 10B). The Kolmogorov–Smirnov displayed a normal distribution of the time constant median (D = 0.09110, p > 0.20). ANOVA with repeated measures was designed with the number of hidden units as independent variable and time constant median as dependent variable. ANOVA analysis shows an interesting effect of the number of hidden neurons [F_{(16, 96)} = 3.6245, p < 0.0001]. Overall, post-hoc analysis revealed that the distribution of the time constants were different when the network was small (less than 5 hidden units). It also reveals that the networks with 80 hidden neurons were differently distributed than medium-sized network (8–10 hidden neurons).

Discussion

Main Finding

We show here that a fully connected recurrent neural network is able to reproduce human walking pattern based on oscillatory properties of kinematics. Although, this network is a black box model without prewired structure mimicking a physiological CPG, its actual performance allows direct comparison with CPG dedicated structure and related algorithm (Duvinage et al., 2012a). Moreover, by the inherent input-output mapping, the DRNN models not only the CPG but also neural feedback pathways and the musculoskeletal system. For simplicity, we consider this neural network as “CPG-like structure” here. We proved that the DRNN is capable of generating the kinematics as elevation angles pattern of walking for both limbs (six degrees of freedom) from simple oscillations corresponding to the three main harmonics of the walking kinematics. Moreover, by modulating those frequencies and tuning them in amplitude as input, the DRNN was able to learn and reliably predict walking kinematics at different velocities (Figure 9). After this appropriate learning the DRNN can thus be considered as a CPG-like structure that would continuously receive oscillation inputs to produce the relevant elevation patterns of the six leg segments. Another interesting result is observed when looking at the structure of best CPGs obtained after learning. Hence, it appears that all of them contain a major part of negative connection weight between units.

Limitations of the Present Approach

Obviously, there are infinite different ways to train the DRNN and this is a strength of this approach. However, it implies a corollary limitation as not all possibilities could be tested in the present study. For example, in order to better document its generalization ability, the model could be trained for a defined low range of velocities, e.g., from 1 to 4 km/h, and then tested with unlearned oscillation input corresponding to higher velocities. A reverse procedure could also be made, i.e., from faster training to slower predictions. Furthermore, inter-subjects generalization has not been studied in the present investigation. The actual usefulness of performing this would largely depend on the basic or application purposes. Another limitation of the present work is the lack of feedback testing, which necessitates a priori identification of a reliable signal and a new operational strategy for learning. Future work will address these aspects specifically.

Neurophysiological Similarity Between Modeled CPG and CPG in Humans and Other Mammals

The understanding of CPG mechanisms remains central in locomotion study (Grillner, 2006; Kiehn, 2006; Rossignol and Frigon, 2011). The CPG is a spinal network of neurons capable of generating a rhythmic pattern of alternate activities between flexor and extensor motoneurons on the same side with reciprocal activation of homologous motoneurons in the contralateral limb. This intrinsic spinal circuitry has been well described in many invertebrate and vertebrate animals, and is highly conserved even in humans, where greater cortical control of spinal modules is required working in conjunction with sensory feedback (Calancie et al., 1994; Bussel et al., 1996; Duysens and Van de Crommert, 1998; Drew et al., 2002; Rossignol et al., 2009). The unique characteristics of human walking probably reflect a complex neural mechanism responsible for pattern production. It is therefore difficult to directly extend experimental findings obtained in quadruped animals to human walking (Barbeau et al., 1998; Capaday, 2002).

The fact that some patients with incomplete spinal injury can move their legs in a rhythmic fashion (Dietz et al., 1995) and that the primary sensorimotor cortex provides oscillatory commands toward the spine during walking (La Fougère et al., 2010) motivates new experiments in which different types of oscillatory signals could be used as input to the CPG-like DRNN. In this context, recent studies have showed EEG oscillations in relation to the gait cycle phase including event-related spectral perturbation in the alpha-beta and gamma bands (Gwin et al., 2011; Haefeli et al., 2011; Cheron et al., 2012; Wagner et al., 2012).

These results are consistent with a top down control of locomotion (Capaday, 2002) and demonstrate the feasibility of extracting EEG signals from the sensorimotor cortex controlling the contralateral foot placement during walking. Although the distinction between the brain signals directly linked to the motor commands and those related to the treatment of multiple sensory signals is a hard task. In this context, Petersen et al. (2012) have found evidence of synchrony in the frequency domain between the primary motor cortex and the tibialis anterior muscle prior to heel strike during the swing phase of walking signifying that rhythmical cortical activity is transmitted via the corticospinal tract to the active muscles. Additionally Wagner et al. (2012), showed a significant difference in the alpha (8–12 Hz) and beta (18–21 Hz) rhythm recorded over the central midline area between passive and active walking with exoskeleton. The role played by specific oscillations related to the initiation and control of human locomotion coming from supraspinal structure was recently demonstrated by local field potential recordings performed in the pedunculopontine nucleus in parkinsonian patients during rest and unconstrained walking (Thevathasan et al., 2012). Alpha oscillation recorded in the caudal part of this nucleus is correlated with gait speed and permits to suppress “task irrelevant” distraction for improving gait performance. Moreover, these authors showed that gait freezing of parkinsonian patient was associated with the attenuation of these alpha waves.

Consistent with this aspect of gait physiology, in our model, input sine waves are sufficient to predict successful output with the DRNN and offer the possibility to mimic such type of supraspinal oscillatory input. It could be advanced that nonlinear mapping of sinusoidal oscillations to kinematic pattern should be realized by other mathematical functions, such as by a Taylor series, but such multi-dimensional approximation seems to be highly difficult to obtain and does not permit testing different network configurations mimicking biological organization such as the CPG. In the present study, we focused on input oscillation derived from the first three harmonics of kinematic signals, which were slower than the alpha frequency range. However, it could be possible to extract slower oscillation from alpha or beta derived signals (envelope) in order to activate the DRNN in our future work.

Structural Similarities of Modeled CPG and Neurophysiological CPGs in Animals

It has been suggested that the biological CPG is considered to serve two basic functions: rhythm generation (RG) and pattern formation (PF).

Initially proposed by Perret and Cabelguen (1980), the idea that the biological CPG is composed of a rhythm and a pattern-amplitude generator is now widely accepted (Kriellaars et al., 1994; Guertin et al., 1995; Perreault et al., 1995; Grillner, 2006; Kiehn, 2006; Talpalar et al., 2011) and paved the way to more complex models of multi-level CPG (McCrea and Rybak, 2008, see below).

It is well-recognized that rhythm generating networks can be realized by means of (1) pacemaker neurons with intrinsic membrane properties such as those described in the stomatogastric ganglion of crustaceans or in the mammalian thalamus (Steriade and Llinás, 1988) or (2) most simplistic neurons without intrinsic pacemaker properties but interacting with inhibitory synapses for producing oscillation as emergent properties of this neuronal population (Geisler et al., 2005). Both neuronal systems thus present the fundamental ability to oscillate. Firstly described in the tadpole and lamprey CPGs, glutamatergic excitatory neurons distributed along the cord (Grillner, 2003) assume the function of rhythmic generator by driving motor neurons and other ipsilateral and commissural inhibitory neurons coordinating the different CPG modules. By blocking the inhibitory networking in the lamprey and also in rodent and cat, many authors (see Kiehn, 2006 for a review) have demonstrated that the glutamatergic burst neurons are the generators of the CPG rhythm.

In addition to intrinsic RG properties, the walking CPGs need to integrate the ipsilateral coordination of flexors and extensors across the same or different joints in a limb and perform inter-limbs coordination. It has been proposed (Zhong et al., 2012) that a subpopulation of neuronal CPG that drives extensor activity is tonically active and is regulated via inhibitory interactions with another CPG rhythmic structure responsible for flexors activity in the same hemicord. This assumption may explain why, during experimental recordings on the neonatal mouse isolated spinal cord, spontaneous deletions of extensor activity do not perturb rhythmic flexor activity. Thus, the inhibitory interneurons play a major role in the temporal sculpting and coordination of the CPG units. The interneurons and the Renshaw cells are involved in this function and in the regulation of walking speed. In addition, the left-right coordination is assumed by a complex network of excitatory and inhibitory commissural interneurons acting on both motor neurons and inhibitory interneurons of the contralateral side (Kiehn, 2006). Interestingly, we have shown that a great percentage of artificial neurons became inhibitory neurons (negative synaptic weight) when the number of neurons progressively increases in the DRNN structure. In this context, it was recently demonstrated in awake mice that the spiking activity of inhibitory neurons of the barrel cortex is organized in order to balance excitation and prevent explosive activity in the recurrently connected cortical microcircuit (Gentet et al., 2010). This physiological mechanism can also be proposed in the present case of the emergent structure of artificial DRNN circuit. Another, not exclusive explanation can reside in the prevalence of inhibitory recurrent connections for producing network oscillation (Geisler et al., 2005; Wildie and Shanahan, 2011). In their review, Nishimaru and Kakizaki (2009) have proposed that inhibitory interneurons play a major role in the CPG of rodent spinal cord. The interneurons are likely to control the bursting of motor neurons during locomotion and it appears that the synaptic transmission mediated by glycine and GABA shifts from excitatory to inhibitory during the prenatal period. It was recently demonstrated that in the absence of glutamatergic synaptic transmission, the flexor-extensor alternation appears to be generated by the inhibitory interneurons, mediating reciprocal inhibition from muscle proprioceptors to antagonist motor neurons (Talpalar et al., 2011).

The present artificial model does not pretend to mimic the complexity of the CPG structure. Instead, it presents a highly simplified recurrent organization from which CPG-like dynamic function emerges, following appropriate learning. Sinusoidal inputs serve as temporal referent to produce rhythmic angles patterns. This model could correspond to the RG structure described previously as a higher order structure that determines rhythmic output of the system (McCrea and Rybak, 2008; Zhong et al., 2012) since sine waves are transformed into kinematics. Another, lower order structure responsible for the phasing and intensity coordination (McCrea and Rybak, 2008; Zhong et al., 2012) could be assumed by another model of DRNN transforming theoretical kinematics into practical muscle orders. We have already studied such relation where EMG signals from walking where used to predict kinematics (Cheron et al., 2003). To conclude to this point, we propose that the two model driving two specific DRNNs (one for producing elevation angles from sine waves and one producing muscular patterns from elevation angles during walking) could act as a complementary top-down pathway to produce adequate coordinated patterns as it has been proposed to model the locomotion in spinal mouse (Zhong et al., 2012).

The present results can also be discussed in the light of the electrical stimulations performed in the MLR) inducing locomotor behavior in decerebrate cat (Shik et al., 1969) or in lamprey (McClellan and Grillner, 1984). In mice, prolonged rhythmic stimulation on the midline of the caudal hindbrain or the ventral spinal cord (C1–C4) induces a stable locomotor activity (Talpalar et al., 2011). Typically, low-frequency stimulation leads to slow-frequency movements and inversely fast-frequency stimulation leads to fast-frequency movements. Our model is in accordance with this physiological behavior, when amplitude of the artificial sine wave inputs increases, the amplitude of stepping increases as well, leading to a change in walking velocities (Figure 6B). In terms of neurological development, there is some evidence for the existence of CPG very early in CNS maturation (Yang et al., 1998; Dominici et al., 2011). Neonatal, so-called “infant” stepping has been ascribed to similar EMG patterns activity in different directions inducing stereotyped yet non-functional walking patterns (Lamb and Yang, 2000). This leads the authors to conclude that the same CPG controls different stepping in human infants in contrast with some studies in adults (Thorstensson, 1986; Grasso et al., 1998). Interestingly, we found that DRNN with only four hidden artificial neurons can generate walking pattern, whereas at least 50 hidden neurons are the prerequisite to generate accurate movements (Figure 7). Obviously, recruitment and training of such high numbers of neurons requires long computational time. For example, with a 4 hidden units DRNN structure, the learning process lasts about 5 min based on our computer performance while 160 min are necessary for a DRNN containing 80 units (Figure 4).

Perspectives

Such a tool can be used to produce gait kinematics in numerous and various applications. For rehabilitation it can be used to train people for recovering a walking pattern corresponding to their physical characteristics by training with an appropriate feedback. Specifically dedicated DRNN based on the proper dynamics of participants could be used for medical applications such as in prosthesis and exoskeleton control (Cheron et al., 2012). It can also be integrated in BCI applications in which higher order commands can be used, e.g., from steady state visual or somatosensory evoked potentials (Cheron et al., 2012) or P300 (Castermans et al., 2011; Duvinage et al., 2012b). This neuronal avenue might lead to the decoding of higher neuronal commands that govern CPG mechanisms. Since these CPG can be trained using specific sinusoidal frequency signals, it might be possible to extract this type of signals from specific EEG rhythms since the brain itself is an effective machine for producing oscillations (Buzsáki and Draguhn, 2004). One of the strengths of this approach is that it is not necessary to determine in advance the topology and the timing sequences between the artificial neurons. This contrasts with other CPGs, such as a recently developed ones (Duvinage et al., 2011) based on coupled oscillators (Righetti and Ijspeert, 2006), where adjustment of intrinsic parameters by optimization techniques was necessary.

In future studies, by introducing an informational delay (Draye et al., 1997) or an artificial distance based on a Gaussian factor modulating the weights between the different neurons (Draye et al., 2002), it will be possible to analyze the self-tailored organization of the links between neurons and the possible emergence of specific topologies. In this case it will also be possible to select different modular architectures of the DRNN.

Conflict of Interest Statement

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

This work was funded by the Belgian Federal Science Policy Office, the European Space Agency, (AO-2004, 118), the Belgian National Fund for Scientific Research (FNRS), the research funds of the Université Libre de Bruxelles and of the Université de Mons (Belgium), the FEDER support (BIOFACT), and the MINDWALKER project (FP7 – 2007–2013) supported by the European Commission. M. Duvinage is a FNRS Research Fellow. This paper presents research results of the Belgian Network DYSCO (Dynamical Systems, Control, and Optimization), funded by the Interuniversity Attraction Poles Program, initiated by the Belgian State, Science Policy Office. The scientific responsibility rests with its author(s).

References

Barbeau, H., Pépin, A., Norman, K. E., Ladouceur, M., and Leroux, A. (1998). Review?: walking after spinal cord injury: control and recovery. Neuroscientist 4, 14–24.

Barliya, A., Omlor, L., Giese, M. A., and Flash, T. (2009). An analytical formulation of the law of intersegmental coordination during human locomotion. Exp. Brain Res. 193, 371–385.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bengoetxea, A., Dan, B., Leurs, F., Cebolla, A. M., De Saedeleer, C., Gillis, P., et al. (2010). Rhythmic muscular activation pattern for fast figure-eight movement. Clin. Neurophysiol. 121, 754–765.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bengoetxea, A., Leurs, F., Cebolla, A., Wellens, S., Draye, J. P., and Cheron, G. (2005). A dynamic recurrent neural network for drawing multi-directional trajectories. Comput. Methods Biomech. Biomed. Engin. 8, 29–30.

Bernstein, N. A. (1967). The Co-Ordination and Regulation of Movements. Oxford; New York: Pergamon Press.

Bianchi, L., Angelini, D., Orani, G. P., and Lacquaniti, F. (1998). Kinematic coordination in human gait: relation to mechanical energy cost. J. Neurophysiol. 79, 2155–2170.

Pubmed Abstract | Pubmed Full Text

Borghese, N. A., Bianchi, L., and Lacquaniti, F. (1996). Kinematic determinants of human locomotion. J. Physiol. 494, 863–879.

Pubmed Abstract | Pubmed Full Text

Brown, T. G. (1911). The intrinsic factors in the act of progression in the mammal. R. Soc. Lond. Proc. Ser. B 84, 308–319.

Bussel, B., Roby-Brami, A., Néris, O. R., and Yakovleff, A. (1996). Evidence for a spinal stepping generator in man. Electrophysiological study. Acta Neurobiol. Exp. (Wars) 56, 465–468.

Pubmed Abstract | Pubmed Full Text

Buzsáki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks. Science 304, 1926–1929.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Calancie, B., Needham-Shropshire, B., Jacobs, P., Willer, K., Zych, G., and Green, B. A. (1994). Involuntary stepping after chronic spinal cord injury evidence for a central rhythm generator for locomotion in man. Brain 117, 1143–1159.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Capaday, C. (2002). The special nature of human walking and its neural control. Trends Neurosci. 25, 370–376.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Castermans, T., Duvinage, M., Petieau, M., Hoellinger, T., Saedeleer, C. D., Seetharaman, K., et al. (2011). Optimizing the performances of a P300-based Brain-Computer Interface in ambulatory conditions. IEEE J. Emerg. Sel. Topics Circuits Syst. 1, 566–577.

Ceccato, J.-C., Azevedo-Coste, C., and Cazalets, J.-R. (2009). Real time control of a CPG-based model of the human trunk in different walking conditions. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2009, 1388–1391.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cheron, G., Bengoetxea, A., Bouillot, E., Lacquaniti, F., and Dan, B. (2001a). Early emergence of temporal co-ordination of lower limb segments elevation angles in human locomotion. Neurosci. Lett. 308, 123–127.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cheron, G., Bouillot, E., Dan, B., Bengoetxea, A., Draye, J. P., and Lacquaniti, F. (2001b). Development of a kinematic coordination pattern in toddler locomotion: planar covariation. Exp. Brain Res. 137, 455–466.

Pubmed Abstract | Pubmed Full Text

Cheron, G., Cebolla, A., Bengoetxea, A., Leurs, F., and Dan, B. (2007). Recognition of the physiological actions of the triphasic EMG pattern by a dynamic recurrent neural network. Neurosci. Lett. 414, 192–196.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cheron, G., Duvinage, M., Castermans, T., Leurs, F., Cebolla, A., Bengoetxea, A., et al. (2011). “Toward an integrative dynamic recurrent neural network for sensorimotor coordination dynamics,” in Recurrent Neural Networks for Temporal Data Processing, Chapter 5, ed H. Cardot (INTECH), 67–80.

Cheron, G., Duvinage, M., De Saedeleer, C., Castermans, T., Bengoetxea, A., Petieau, M., et al. (2012). From spinal central pattern generators to cortical network: integrated BCI for walking rehabilitation. Neural Plast. 2012:375148. doi: 10.1155/2012/375148

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cheron, G., Leurs, F., Bengoetxea, A., Draye, J. P., Destrée, M., and Dan, B. (2003). A dynamic recurrent neural network for multiple muscles electromyographic mapping to elevation angles of the lower limb in human locomotion. J. Neurosci. Methods 129, 95–104.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Churchland, M. M., Cunningham, J. P., Kaufman, M. T., Foster, J. D., Nuyujukian, P., Ryu, S. I., et al. (2012). Neural population dynamics during reaching. Nature 487, 51–56.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Courtine, G., and Schieppati, M. (2004). Tuning of a basic coordination pattern constructs straight-ahead and curved walking in humans. J. Neurophysiol. 91, 1524–1535.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Dietz, V., Colombo, G., Jensen, L., and Baumgartner, L. (1995). Locomotor capacity of spinal cord in paraplegic patients. Ann. Neurol. 37, 574–582.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Dominici, N., Daprati, E., Nico, D., Cappellini, G., Ivanenko, Y. P., and Lacquaniti, F. (2009). Changes in the limb kinematics and walking-distance estimation after shank elongation: evidence for a locomotor body schema? J. Neurophysiol. 101, 1419–1429.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Dominici, N., Ivanenko, Y. P., Cappellini, G., D' Avella, A., Mondì, V., Cicchese, M., et al. (2011). Locomotor primitives in newborn babies and their development. Science 334, 997–999.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Draye, J., Pavisic, D. A., Cheron, G. A., and Libert, G. A. (1996). Dynamic recurrent neural networks: a dynamical analysis. IEEE Trans. Syst. Man Cybern. B Cybern. 26, 692–706.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Draye, J. P., Cheron, G., Libert, G., and Godaux, E. (1997). Emergence of clusters in the hidden layer of a dynamic recurrent neural network. Biol. Cybern. 76, 365–374.

Pubmed Abstract | Pubmed Full Text

Draye, J.-P., Pavisic, D., Cheron, G., and Libert, G. (1995). Adaptative time constants improve the prediction capability of recurrent neural networks. Neural Process. Lett. 2, 12–16.

Draye, J.-P., Winters, J. M., and Cheron, G. (2002). Self-selected modular recurrent neural networks with postural and inertial subnetworks applied to complex movements. Biol. Cybern. 87, 27–39.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Drew, T., Jiang, W., and Widajewicz, W. (2002). Contributions of the motor cortex to the control of the hindlimbs during locomotion in the cat. Brain Res. Rev. 40, 178–191.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Duvinage, M., Castermans, T., Hoellinger, T., and Reumaux, J. (2012a). Human walk modeled by PCPG to control a lower limb neuroprosthesis by high-level commands. J. Syst. Cybern. Inform. 10, 70–80.

Duvinage, M., Castermans, T., Jimenez-Fabian, R., Petieau, M., Seetharaman, K., Hoellinger, T., et al. (2012b). “A five-state P300-based foot lifter orthosis: proof of concept,” in The 3rd IEEE Biosignals and Biorobotics conference (ISSNIP) (Manaus), 1–6.

Duvinage, M., Castermans, T., Hoellinger, T., Cheron, G., and Dutoit, T. (2011). “Modeling human walk by PCPG for lower limb neuroprosthesis control,” in 5th International IEEE/EMBS Conference on Neural Engineering (NER) (Cancun), 317–321.

Duysens, J., and Van de Crommert, H. (1998). Neural control of locomotion; the central pattern generator from cats to humans. Gait Posture 7, 131–141.

Pubmed Abstract | Pubmed Full Text

Flash, T., and Hochner, B. (2005). Motor primitives in vertebrates and invertebrates. Curr. Opin. Neurobiol. 15, 660–666.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Geisler, C., Brunel, N., and Wang, X.-J. (2005). Contributions of intrinsic membrane dynamics to fast network oscillations with irregular neuronal discharges. J. Neurophysiol. 94, 4344–4361.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gentet, L. J., Avermann, M., Matyas, F., Staiger, J. F., and Petersen, C. C. H. (2010). Membrane potential dynamics of GABAergic neurons in the barrel cortex of behaving mice. Neuron 65, 422–435.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Georgopoulos, A. P., and Grillner, S. (1989). Visuomotor coordination in reaching and locomotion. Science 245, 1209–1210.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Grasso, R., Bianchi, L., and Lacquaniti, F. (1998). Motor patterns for human gait: backward versus forward locomotion. J. Neurophysiol. 80, 1868–1885.

Pubmed Abstract | Pubmed Full Text

Grillner, S. (2003). The motor infrastructure: from ion channels to neuronal networks. Nat. Rev. Neurosci. 4, 573–586.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Grillner, S. (2006). Biological pattern generation: the cellular and computational logic of networks in motion. Neuron 52, 751–766.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Guertin, P., Angel, M. J., Perreault, M. C., and McCrea, D. A. (1995). Ankle extensor group I afferents excite extensors throughout the hindlimb during fictive locomotion in the cat. J. Physiol. (Lond.) 487(Pt 1), 197–209.

Pubmed Abstract | Pubmed Full Text

Gwin, J. T., Gramann, K., Makeig, S., and Ferris, D. P. (2011). Electrocortical activity is coupled to gait cycle phase during treadmill walking. Neuroimage 54, 1289–1296.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Haefeli, J., Vögeli, S., Michel, J., and Dietz, V. (2011). Preparation and performance of obstacle steps: interaction between brain and spinal neuronal activity. Eur. J. Neurosci. 33, 338–348.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hellgren, J., Grillner, S., and Lansner, A. (1992). Computer simulation of the segmental neural network generating locomotion in lamprey by using populations of network interneurons. Biol. Cybern. 68, 1–13.

Pubmed Abstract | Pubmed Full Text

Hicheur, H., Terekhov, A. V., and Berthoz, A. (2006). Intersegmental coordination during human locomotion: does planar covariation of elevation angles reflect central constraints? J. Neurophysiol. 96, 1406–1419.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ijspeert, A. J. (2008). Central pattern generators for locomotion control in animals and robots: a review. Neural Netw. 21, 642–653.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ijspeert, A. J., Crespi, A., Ryczko, D., and Cabelguen, J.-M. (2007). From swimming to walking with a salamander robot driven by a spinal cord model. Science 315, 1416–1420.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ivanenko, Y. P., Cappellini, G., Dominici, N., Poppele, R. E., and Lacquaniti, F. (2007). Modular control of limb movements during human locomotion. J. Neurosci. 27, 11149–11161.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ivanenko, Y. P., D'Avella, A., Poppele, R. E., and Lacquaniti, F. (2008). On the origin of planar covariation of elevation angles during human locomotion. J. Neurophysiol. 99, 1890–1898.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ivanenko, Y. P., Dominici, N., Cappellini, G., and Lacquaniti, F. (2005). Kinematics in newly walking toddlers does not depend upon postural stability. J. Neurophysiol. 94, 754–763.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ivanenko, Y. P., Grasso, R., Macellari, V., and Lacquaniti, F. (2002). Control of foot trajectory in human locomotion: role of ground contact forces in simulated reduced gravity. J. Neurophysiol. 87, 3070–3089.

Pubmed Abstract | Pubmed Full Text

Kiehn, O. (2006). Locomotor circuits in the mammalian spinal cord. Annu. Rev. Neurosci. 29, 279–306.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kriellaars, D. J., Brownstone, R. M., Noga, B. R., and Jordan, L. M. (1994). Mechanical entrainment of fictive locomotion in the decerebrate cat. J. Neurophysiol. 71, 2074–2086.

Pubmed Abstract | Pubmed Full Text

Lacquaniti, F., Grasso, R., and Zago, M. (1999). Motor patterns in walking. Physiology 14, 168–174.

Pubmed Abstract | Pubmed Full Text

Lacquaniti, F., Ivanenko, Y. P., and Zago, M. (2002). Kinematic control of walking. Arch. Ital. Biol. 140, 263–272.

Pubmed Abstract | Pubmed Full Text

La Fougère, C., Zwergal, A., Rominger, A., Förster, S., Fesl, G., Dieterich, M., et al. (2010). Real versus imagined locomotion: a [18F]-FDG PET-fMRI comparison. Neuroimage 50, 1589–1598.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lamb, T., and Yang, J. F. (2000). Could different directions of infant stepping be controlled by the same locomotor central pattern generator? J. Neurophysiol. 83, 2814–2824.

Pubmed Abstract | Pubmed Full Text

Latash, M. L. (2008). Synergy, Edn 1. Oxford: Oxford University Press.

Leurs, F., Bengoetxea, A., Cebolla, A. M., De Saedeleer, C., Dan, B., and Cheron, G. (2012). Planar covariation of elevation angles in prosthetic gait. Gait Posture 35, 647–652.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Leurs, F., Ivanenko, Y. P., Bengoetxea, A., Cebolla, A.-M., Dan, B., Lacquaniti, F., et al. (2011). Optimal walking speed following changes in limb geometry. J. Exp. Biol. 214, 2276–2282.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

McClellan, A. D., and Grillner, S. (1984). Activation of “fictive swimming” by electrical microstimulation of brainstem locomotor regions in an in vitro preparation of the lamprey central nervous system. Brain Res. 300, 357–361.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

McCrea, D. A., and Rybak, I. A. (2008). Organization of mammalian locomotor rhythm and pattern generation. Brain Res. Rev. 57, 134–146.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Nishimaru, H., and Kakizaki, M. (2009). The role of inhibitory neurotransmission in locomotor circuits of the developing mammalian spinal cord. Acta Physiol. (Oxf.) 197, 83–97.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Patla, A. E., Calvert, T. W., and Stein, R. B. (1985). Model of a pattern generator for locomotion in mammals. Am. J. Physiol. 248, R484–R494.

Pubmed Abstract | Pubmed Full Text

Pearlmutter, B. (1990). “Dynamic recurrent neural networks,” in Carnegie Mellon University Repository. Available online at: http://repository.cmu.edu/about.html

Perreault, M. C., Angel, M. J., Guertin, P., and McCrea, D. A. (1995). Effects of stimulation of hindlimb flexor group II afferents during fictive locomotion in the cat. J. Physiol. (Lond.) 487(Pt 1), 211–220.

Pubmed Abstract | Pubmed Full Text

Perret, C., and Cabelguen, J. M. (1980). Main characteristics of the hindlimb locomotor cycle in the decorticate cat with special reference to bifunctional muscles. Brain Res. 187, 333–352.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Petersen, T. H., Willerslev-Olsen, M., Conway, B. A., and Nielsen, J. B. (2012). The motor cortex drives the muscles during walking in human subjects. J. Physiol. (Lond.) 590, 2443–2452.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Prentice, S. D., Patla, A. E., and Stacey, D. A. (1998). Simple artificial neural network models can generate basic muscle activity patterns for human locomotion at different speeds. Exp. Brain Res. 123, 474–480.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Prentice, S. D., Patla, A. E., and Stacey, D. A. (2001). Artificial neural network model for the generation of muscle activation patterns for human locomotion. J. Electromyogr. Kinesiol. 11, 19–30.

Pubmed Abstract | Pubmed Full Text

Righetti, L., Buchli, J., and Ijspeert, A. J. (2005). “From dynamic hebbian learning for oscillators to adaptive central pattern generators,” in Proceedings of 3rd International Symposium on Adaptive Motion in Animals and Machines – AMAM 2005 (Ilmenau: Verlag ISLE), 45.

Righetti, L., and Ijspeert, A. J. (2006). “Programmable central pattern generators: an application to biped locomotion control,” in Proceedings of the 2006 IEEE International Conference on Robotics and Automation, BIOROB-CONF-2006-011 (Pisa), 1585–1590.

Rossignol, S., Barrière, G., Alluin, O., and Frigon, A. (2009). Re-expression of locomotor function after partial spinal cord injury. Physiology (Bethesda) 24, 127–139.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rossignol, S., and Frigon, A. (2011). Recovery of locomotion after spinal cord injury: some facts and mechanisms. Annu. Rev. Neurosci. 34, 413–440.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Shik, M. L., Severin, F. V., and Orlovsky, G. N. (1969). Control of walking and running by means of electrical stimulation of the mesencephalon. Electroencephalogr. Clin. Neurophysiol. 26, 549.

Pubmed Abstract | Pubmed Full Text

Steriade, M., and Llinás, R. R. (1988). The functional states of the thalamus and the associated neuronal interplay. Physiol. Rev. 68, 649–742.

Pubmed Abstract | Pubmed Full Text

Talpalar, A. E., Endo, T., Löw, P., Borgius, L., Hägglund, M., Dougherty, K. J., et al. (2011). Identification of minimal neuronal networks involved in flexor-extensor alternation in the mammalian spinal cord. Neuron 71, 1071–1084.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Thevathasan, W., Pogosyan, A., Hyam, J. A., Jenkinson, N., Foltynie, T., Limousin, P., et al. (2012). Alpha oscillations in the pedunculopontine nucleus correlate with gait performance in parkinsonism. Brain 135, 148–160.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Thorstensson, A. (1986). How is the normal locomotor program modified to produce backward walking? Exp. Brain Res. 61, 664–668.

Pubmed Abstract | Pubmed Full Text

Wadden, T., Hellgren, J., Lansner, A., and Grillner, S. (1997). Intersegmental coordination in the lamprey: simulations using a network model without segmental boundaries. Biol. Cybern. 76, 1–9.

Wagner, J., Solis-Escalante, T., Grieshofer, P., Neuper, C., Müller-Putz, G., and Scherer, R. (2012). Level of participation in robotic-assisted treadmill walking modulates midline sensorimotor EEG rhythms in able-bodied subjects. Neuroimage 63, 1203–1211.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wildie, M., and Shanahan, M. (2011). Establishing communication between neuronal populations through competitive entrainment. Front. Comput. Neurosci. 5:62. doi: 10.3389/fncom.2011.00062

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yang, J. F., Lamont, E. V., and Pang, M. Y. C. (2005). Split-belt treadmill stepping in infants suggests autonomous pattern generators for the left and right leg in humans. J. Neurosci. 25, 6869–6876.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yang, J. F., Stephens, M. J., and Vishram, R. (1998). Infant stepping: a method to study the sensory control of human walking. J. Physiol. (Lond.) 507(Pt 3), 927–937.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Zhong, G., Shevtsova, N. A., Rybak, I. A., and Harris-Warrick, R. M. (2012). Neuronal activity in the isolated mouse spinal cord during spontaneous deletions in fictive locomotion: insights into locomotor central pattern generator organization. J. Physiol. (Lond.) 590, 4735–4759.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text