Experimental protocol

Eighteen healthy and young adults were recruited (11 male, 7 female, height 176 ± 7 cm, body mass 71 ± 13 kg, age 24 ± 3 years, means ± s.d.). All the participants were regularly active and did not use orthotic insoles. None had any history of neuromuscular or musculoskeletal impairments, or any head or spine injury at the time of the measurements or in the previous six months. This study was reviewed and approved by the Ethics Committee of the University of Kassel (approval code E05201602). All the participants gave written informed consent for the experimental procedure, in accordance with the Declaration of Helsinki.

The recordings were conducted on two treadmills: a standard one, equipped with an even-surface (ES) belt (Laufergotest, Erich Jäger, Würzburg, Germany) and one equipped with a custom-made uneven-surface (US) belt (Woodway®, Weil am Rhein, Germany, Fig. 1, Supplementary video 1). The US treadmill’s belt consisted of terrasensa® classic modules (Sensa® by Huebner, Kassel, Germany) aiming to simulate uneven ground conditions. The kinematics data were recorded through a six-camera motion capture system (Oqus 3+, Qualisys AB, Gothenburg, Sweden) operating at 300 Hz. The muscle activity of 13 ipsilateral muscles was recorded using a 16-channel wireless EMG system (myon m320, myon AG, Schwarzenberg, Switzerland), with a frequency of 1000 Hz.

Figure 1 Sketch of the uneven-surface treadmill employed in this study. Full size image

The participants completed two different tasks on the ES and US treadmill, in random order: walking (shod, fixed velocity; 1.1 m/s female, 1.2 m/s male) and running (shod, fixed velocity; 2.0 m/s female, 2.2 m/s male). The difference in velocity between male and female have been chosen after a pilot study in which we estimated the average gender-specific comfortable running and walking velocity on the US treadmill. The participants were instructed to keep looking at a fixed spot in front of them and avoid looking at the treadmill’s belt.

Gait cycle assessment

Ten reflective markers were placed bilaterally on the leg, foot and spine. Namely, the greater trochanter, the Achilles tendon insertion on the calcaneus, the dorsal margin of the fifth metatarsal head, the second, seventh and tenth thoracic and the second lumbar vertebræ were marked. The gait cycle breakdown was obtained from kinematic data. We used the information coming from the calcaneus and fifth metatarsal markers. These data were low-pass filtered using a 4th order IIR Butterworth zero-phase filter with cut-off frequency of 50 Hz52. For estimating touchdown, we used the modified foot contact algorithm developed by Maiwald et al.52. For estimating lift-off, we developed the foot acceleration and jerk algorithm. The jerk algorithm searches for the global maximum of the fifth metatarsal vertical acceleration between two consecutive touchdown events to estimate the lift-off (LOe, where the “e” stays for “estimated”). This estimation, however, does not provide an accurate identification of the lift-off and needs some refinement. To get closer to the “real” lift-off timing, a characteristic minimum in the vertical acceleration (i.e. when the jerk equals zero) of the fifth metatarsal marker is identified in a reasonably small neighbourhood of the LOe. We found [LOe – 90 ms, LOe + 25 ms] for walking and [LOe – 50 ms, LOe + 100 ms] for running to be the sufficiently narrow intervals needed to make the initial lift-off estimation. Both approaches have been validated using force plate data (AMTI BP600, Advanced Mechanical Technology, Inc., Watertown, MA, USA) from 15 participants walking and running overground at six different velocities. True errors were of −8 ± 8 ms (−1.1% ± 1.0% of the stance phase) for touchdown and 12 ± 18 ms (1.6% ± 2.4%) for lift-off for walking and −1 ± 15 ms (−0.5% ± 5.4%) and −16 ± 23 ms (−5.7% ± 8.3%) respectively for running (means ± s.d.).

Local dynamic stability assessment

Participants were allowed for an accommodation period of maximum 60 s53. During walking, two trials of 180 s were recorded for all participants, while two trials of 120 s were recorded during the running task. A 4th order IIR Butterworth zero-phase filter with low-pass cut-off frequency of 20 Hz was applied to the computed coordinates of the four spine markers (i.e. the second, seventh and tenth thoracic and the second lumbar vertebræ). We adopted the maximum finite-time Lyapunov exponent (MLE) to assess the local dynamic stability of the human system during walking and running. The coordinates of the walking trials were downsampled to 100 Hz to improve computational performance. The vertical coordinates of these markers were then clustered to be used for further analysis and the calculation of the MLE. We calculated the MLE using the vertical coordinate data of the clustered spine markers, which we tested for stationarity54,55. Running on a treadmill restricts the movement at the anteroposterior direction due to the participants seeking to match the velocity of the treadmill and similarly, the treadmill width restricts the movement on the mediolateral direction. To avoid dependencies on step frequency, we identified the maximum number of shared steps (i.e. 0.5 of gait cycle) for all trials and participants56; 291 and 287 steps for walking and running, respectively. The coordinates of the data segments corresponding to the exact number of steps were then isolated for each trial. Following, the data segments were normalised to a uniform length based on the average number of data for each step, amounting ~16000 in walking and ~33000 data points in running. The high number of steps analysed ensured the reliability of our measurements57.

State space reconstruction was achieved through delay coordinate embedding58,59, for each point of the time series and its time-delayed copies as follows:

$$S(t)=[z(t),z(t+\tau ),\ldots ,z(t+(d-1)\tau )]$$ (1)

with S(t) being the d-dimensional reconstructed state vector, z(t) the input 1D coordinate series, τ the time delay and d the embedding dimension. Time delays were calculated for each time series from the first minimum of the mutual-information curve, based on the Average Mutual Information function60. The number of embedding dimensions was extracted through a Global False Nearest Neighbours analysis for each time series, with a threshold of one per thousand data points61.

Different values of τ and d can yield very different state-space reconstructions54,62,63. It is therefore suggested that optimised values of τ and d are necessary to best represent a dynamical system64. In the current study time delays were individually chosen for each participant64, by first calculating the optimal delay of the four time-series (two trials on the even and two on the uneven surface) and using the averaged value. Time delays were in the range of 17–20 data points (~0.34 of the average step) in the walking task and 36–40 data points (~0.33 of the average step) in the running task.

Following the reconstruction of the times series, the Rosenstein algorithm was used to compute the average exponential rate of divergence by calculating the linear distance of each point’s trajectory divergence to its closest trajectory65,66. The MLE were then calculated from the slope of the linear fit in the resulting divergence curves from 0 to 1 step. Analysis of the data was performed on MATLAB 2014b (Mathworks Inc., Natick, United States).

Spinal motor output assessment

For each condition, the muscle activity of the following 13 ipsilateral (right side) muscles was recorded: gluteus medius (ME), gluteus maximus (MA), tensor fasciæ latæ (FL), rectus femoris (RF), vastus medialis (VM), vastus lateralis (VL), semitendinosus (ST), biceps femoris (long head, BF), tibialis anterior (TA), peroneus longus (PL), gastrocnemius medialis (GM), gastrocnemius lateralis (GL) and soleus (SO). We recorded two trials of 60 s for each condition and analysed the first 50 gait cycles of each acquisition. The EMG signals were high-pass filtered and then full-wave rectified and low-pass filtered using a 4th order IIR Butterworth zero-phase filter with cut-off frequencies 50 Hz (high-pass) and 20 Hz (low-pass for the linear envelope) using R v3.4.1 (R Found. for Stat. Comp.). The amplitude was normalised to the maximum activation recorded for each participant across all conditions43,67,68. Each gait cycle was then time-normalised to 200 points69, assigning 100 points to the stance and 100 points to the swing phase.

For the spinal motor output characterisation, we mapped the EMG activity onto the estimated rostrocaudal location of alpha-motoneurons (MNs) pools in the segments from the second lumbar vertebra (L2) to the second sacral vertebra (S2) of the spinal cord47,70,71. The contribution of each muscle to the total estimated activity of the spinal segments was implemented using the myotomal charts developed by Kendall et al.72. This method shows the organisation of the efferent MNs network directed to the muscles, assuming a common spinal topography among the investigated participants. The motor output of each spinal segment S j was estimated using the Equation (2) introduced by La Scaleia et al.70:

$${S}_{j}=\frac{{\sum }_{i=1}^{{m}_{j}}(\frac{{k}_{ji}}{{n}_{i}}\times EM{G}_{i})}{{\sum }_{i=1}^{{m}_{j}}(\frac{{k}_{ji}}{{n}_{i}})}\times M{N}_{j}$$ (2)

where m j are the muscles innervated by each segment, n i is the number of spinal levels that innervate the ith muscle, k ij is a weighting coefficient specific to each muscle and spinal segment (e.g. k ij = 1 or k ij = 0.5 if S j is a major or minor MN source, respectively) and EMG i is the normalised recorded EMG, specific for each participant and trial70,72. This approach accounts for size differences at each spinal level in every MN j pool.

Modular organisation assessment

Muscle synergies data were extracted through a custom script73 (R v3.4.1, R Found. for Stat. Comp.) using the classical Gaussian NMF algorithm50,73,74 from the first 50 gait cycles of each acquisition51. EMG data were pre-processed using the same filtering conditions reported in the previous paragraph and each gait cycle was time-normalised to 200 points69, assigning 100 points to the stance and 100 points to the swing phase74. The m = 13 time-dependent muscle activity vectors were grouped in an m × n matrix V (n equal to 10000), factorised such that V≈V R = WH. V R represents the new reconstructed matrix, which approximates the original matrix V. The motor primitives73,75 matrix H contained the time-dependent coefficients of the factorisation with dimensions r × n, where r represents the minimum number of synergies necessary to reconstruct the original signals (V). The motor modules73,76 matrix W with dimensions m × r, contained the time-invariant muscle weightings. H and W described the synergies necessary to accomplish a movement. The update rules for H and W are presented in Equations (3) and (4).

$${H}_{i+1}={H}_{i}\frac{{{W}_{i}}^{T}V}{{{W}_{i}}^{T}{W}_{i}{H}_{i}}$$ (3)

$${W}_{i+1}={W}_{i}\frac{V{({H}_{i+1})}^{T}}{{W}_{i}{H}_{i+1}{({H}_{i+1})}^{T}}$$ (4)

The limit of convergence was reached when a change in the calculated R2 between V and V R was smaller than the 0.01% in the last 20 iterations73,77. This was done for a number of synergies successively increased from 1 to 10. The computation was repeated 10 times for each synergy, each time creating new randomised initial matrices H and W, in order to avoid local minima73,78. The solution with the highest R2 was then selected for each of the 10 synergies.

To choose the minimum number of synergies required to represent the original signals, the curve of R2 values versus synergies was fitted using a simple linear regression model, using all ten synergies. The mean squared error73 was then repeatedly calculated, each time removing the lower synergy point, until only two points were left or until the mean squared error fell below 10−5. The aforementioned procedure allowed us to extract fundamental and combined synergies from the raw EMG data. A fundamental synergy can be defined as an activation pattern whose motor primitive shows a single peak of activation73,74. When two or more fundamental synergies are blended into one, a combined synergy appears. Combined synergies were excluded from the analysis. An example of combined synergies is reported in Fig. 2.

Figure 2 Example of two fundamental synergies combined into one. The histograms (a and b) represent the two fundamental sets of motor modules for seven muscles. The curves (d and e) show the two respective primitives, with arbitrary x- and y-axis units. The combined motor modules and primitives are presented in panels (c and f), respectively. Full size image

Metrics for comparison of curves

We evaluated the centre of activity (CoA) and full width at half maximum (FWHM) for the resulting curves of the extracted spinal maps and motor primitives (matrix H) in both conditions and types of locomotion. The CoA was defined as the angle of the vector (in polar coordinates) that points to the centre of mass of that circular distribution69,74. The polar direction represented the gait cycle’s phase, with angle 0 ≤ θ t ≤ 2π. The following equations define the CoA:

$$A=\sum _{t=1}^{p}(\cos \,{\theta }_{t}\times {P}_{t})$$ (5)

$$B=\sum _{t=1}^{p}(\sin \,{\theta }_{t}\times {P}_{t})\,$$ (6)

$$CoA=\arctan (B/A)$$ (7)

where p is the number of points of each gait cycle (p = 200) and P is the activation vector. The FWHM was calculated as the number of points exceeding each gait cycle’s half maximum, after subtracting the gait cycle’s minimum69,74. The first 50 gait cycles of each acquisition were selected for analysis. The CoA and FWHM were analysed step-by-step and then averaged for stance and swing distinctively for spinal maps and over the whole gait cycle for the motor primitives.

Statistics

To evaluate the differences in MLE and gait parameters between ES and US locomotion, we used a one-way ANOVA with repeated measures or non-parametric Friedman test in case the normality assumptions on the residuals were not satisfied. To investigate CoA and FWHM of the spinal motor output, we employed a two-way ANOVA with repeated measures, followed by a Tukey post-hoc analysis with false discovery rate p-value adjustment (the independent variables being the condition, i.e. ES or US, and the spinal segment). The same statistics was used to assess differences between the motor modules, using the muscles and the conditions (ES or US) as independent variables. To evaluate the similarities between the fundamental motor primitives of the ES and US conditions, we used the coefficient of determination R2. The analysis was conducted as follows: first, we calculated the similarity between the pairs of trials recorded during ES and US locomotion (i.e. the repeatability level when comparing two trials of the same condition)73,74. Then, we statistically compared the similarity between trials of different conditions (ES and US locomotion) with the repeatability values. To do so, we used a one-way ANOVA with repeated measures or non-parametric Friedman test in case the normality assumptions on the residuals were not satisfied. Type A uncertainty was expressed as \({u}_{A}=s/\surd n\). All the significance levels were set to α = 0.05 and the statistical analyses were conducted using R v3.4.1 (R Found. for Stat. Comp.).

Data availability

The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.