The Influence of Footwear on the Modular Organization of Running

For most of our history, we predominantly ran barefoot or in minimalist shoes. The advent of modern footwear, however, might have introduced alterations in the motor control of running. The present study investigated shod and barefoot running under the perspective of the modular organization of muscle activation, in order to help addressing the neurophysiological factors underlying human locomotion. On a treadmill, 20 young and healthy inexperienced barefoot runners ran shod and barefoot at preferred speed (2.8 ± 0.4 m/s). Fundamental synergies, containing the time-dependent activation coefficients (motor primitives) and the time-invariant muscle weightings (motor modules), were extracted from 24 ipsilateral electromyographic activities using non-negative matrix factorization. In shod running, the average foot strike pattern was a rearfoot strike, while in barefoot running it was a mid-forefoot strike. In both conditions, five fundamental synergies were enough to describe as many gait cycle phases: weight acceptance, propulsion, arm swing, early swing and late swing. We found the motor primitives to be generally shifted earlier in time during the stance-related phases and later in the swing-related ones in barefoot running. The motor primitive describing the propulsion phase was significantly of shorter duration (peculiarity confirmed by the analysis of the spinal motor output). The arm swing primitive, instead, was significantly wider in the barefoot condition. The motor modules demonstrated analogous organization with some significant differences in the propulsion, arm swing and late swing synergies. Other than to the trivial absence of shoes, the differences might be deputed to the lower ankle gear ratio (and the consequent increased system instability) and to the higher recoil capabilities of the longitudinal foot arch during barefoot compared to shod running.

Introduction

In the last decade, the study of locomotion in evolutionary anthropology has been increasingly focusing on endurance running. Humans, compared to non-human primates, show exceptional endurance running speeds (Bramble and Lieberman, 2004). However, the advent of modern running shoes is contemporary history compared to the two million-years-old fossil evidence of running as a derived capability of the genus Homo (Bramble and Lieberman, 2004; Lieberman et al., 2010). Running barefoot or in minimal footwear has been the predominant condition for most of the human history (Lieberman et al., 2010). Hence, it can be expected that the strategies adopted to run barefoot might differ from those employed to run shod.

During running, the foot can strike the ground in multiple ways, called foot strike patterns (FSPs). Rearfoot (RS), midfoot (MS), and forefoot (FS) strike are the common classifications, depending on the location of the first contact area with the ground (Hasegawa et al., 2007). We recently found that if almost 90% of the population adopts a RS pattern when running shod, only half maintain it when switching to barefoot (Santuz et al., 2016), changing to either MS or FS (joined in an unique pattern and indicated as mid-forefoot strike, MFS for brevity). Adopting MFS patterns can increase the plantarflexors activity, reduce the ground contact times and affect the kinetics and kinematics of the whole gait cycle (Komi, 1984, 1992; von Tscharner et al., 2003; Hasegawa et al., 2007; Lieberman et al., 2010). Therefore, we can argue that switching between the two conditions of running shod and barefoot does not only imply kinematic and kinetic changes, but might involve a different organization of movement. From a motor control perspective, this assumption can be investigated by analyzing the modular organization of muscle activity before and after altering the running condition.

Since the late 1960s (Bernstein, 1967) it has been accepted that the central nervous system can simplify the production of movements by avoiding the activation of each muscle separately (Mussa-Ivaldi et al., 1994; Bizzi et al., 2008). This important feature might be implemented by reducing the degrees of freedom through a linear combination (Mussa-Ivaldi et al., 1994) of specific muscle activation patterns, called synergies (Bizzi et al., 2008). A number of studies were able to show that synergies reside in the brain stem or spinal cord and follow a modular organization (Tresch et al., 2002; Hart and Giszter, 2004; Bizzi et al., 2008; Roh et al., 2011; Bizzi and Cheung, 2013). Recently, a study in mice using optogenetics to isolate the excitatory and inhibitory neuronal populations could show a strong specificity in the spinal cord topography (Hägglund et al., 2013). The synergies as low dimensional units, via descending or afferent pathways, produce a complex electromyographic (EMG) pattern in muscles (Tresch et al., 2002; Bizzi and Cheung, 2013), creating a locomotor drive mediated by a certain amount of supraspinal control (Roh et al., 2011). During walking, the same amount of basic activation patterns could be found in patients with spinal cord injury and in healthy participants at different speeds and gravitational loads (Ivanenko et al., 2006). Synergies similar to those found in humans at a spinal (Ivanenko et al., 2006) or muscular level can be observed also in the motor cortex of the primate and cat (Yakovenko et al., 2011; Overduin et al., 2015). Moreover, studies on the excitability of the corticospinal system showed that training can improve task-specific brain organizations (Moscatelli et al., 2016a,b; Monda et al., 2017). This suggests a high degree of cooperation within the central nervous system's structure at all levels. In this study, we used an unsupervised learning method called non-negative matrix factorization (NMF) (Lee and Seung, 1999) for reducing the high dimensional EMG input into a small number of synergies. We focus on the comparison between the modular organization of shod and barefoot running. Compared to the analysis of direct EMG signals, the muscle synergies concept has the clear advantage of being a high-throughput approach for analyzing muscle activities. In fact, it does not only provide the researcher with an automatic, low-dimensional, clustering of the activations during the gait cycle, but it also identifies the weighted contribution of each muscle for producing a certain movement.

The objective of the current study was to investigate the modular organization of shod and barefoot running using muscle synergies in order to gain new knowledge about the neurophysiological factors underlying human locomotion. Based on reported changes in the kinematic, kinetic (Lieberman et al., 2010) and EMG (Komi, 1984, 1992; von Tscharner et al., 2003; Cappellini et al., 2006, 2010; Ivanenko et al., 2008) characteristics of shod and barefoot running, we hypothesized that there is an alteration in the modular organization between the two conditions that might be associated with the specificity of the respective foot strike with the ground. In a first step we calculated the similarities between trials of the same condition using the coefficient of determination (R²) and set their repeatability (intraday) thresholds (Santuz et al., 2017). Afterwards, we investigated the similarities between the two conditions and compared with the intraday thresholds. This two-step process allowed an improved objective and quantitative interpretation of the conventionally qualitative output given by methods like the NMF.

Materials and Methods

Experimental Protocol

Twenty healthy and young adults were recruited (10 male, height 180 ± 5 cm, body mass 77 ± 8 kg, age 31 ± 7 years, 10 female, height 169 ± 8 cm, body mass 60 ± 8 kg, age 28 ± 5 years). All the participants were regularly active and did not use orthotic insoles and did not have any previous experience of barefoot running. None showed or reported any history of neuromuscular or musculoskeletal impairments, or any head or spine injury at the time of the measurements or in the previous 6 months. This study was reviewed and approved by the Ethics Committee of the Humboldt-Universität zu Berlin. All the participants gave written informed consent for the experimental procedure, in accordance with the Declaration of Helsinki.

The muscle activity of 24 ipsilateral muscles was recorded using one 16-channel (myon m320, myon AG, Schwarzenberg, Switzerland) and one 8-channel (myon RFTD E08, myon AG, Schwarzenberg, Switzerland) wireless surface-EMG systems. The acquisition frequency was set to 1,000 Hz (16 ms latency, constant). Vertical ground reaction forces (VGRFs) were recorded at 120 Hz through a pressure plate (FDM-THM-S, zebris Medical GmbH, Isny im Allgäu, Germany) integrated with a treadmill (mercury, H-p-cosmos Sports and Medical GmbH, Nussdorf, Germany). The pressure plate data were acquired using the proprietary software (WinFDM-T v2.5.1, zebris Medical GmbH, Isny im Allgäu, Germany) and then extracted in a raw format for autonomous post-processing using a validated custom algorithm (Santuz et al., 2016) written in R version 3.4.1 (R Foundation for Statistical Computing, R Core Team, Vienna, Austria). The EMG devices and the plate were synchronized using an analog signal.

The participants completed a self-selected warm-up on the treadmill, in order to choose their comfortable shod-running speed. The procedure to find the comfortable speed was implemented using the method of limits (Treutwein, 1995). The speed was randomly increased with steps of 0.02 to 0.05 m/s at varying time intervals (around 5 to 10 s) until the participant was comfortable with a specific pace. The operation was then repeated starting from a faster speed and randomly decreasing it as previously done. If the comfortable value did not differ more than 10% from the previous, the average of the two values was taken as the preferred. Otherwise, the whole procedure was iterated. The warm-up protocol typically lasted between 5 and 10 min. After being instructed about the protocol, the participants completed two different tasks, in random order: shod running at the preferred running speed (2.9 ± 0.4 m/s for male, 2.6 ± 0.2 m/s for female) and barefoot running at the same speed.

For each condition, the muscle activity of the 24 ipsilateral (right side) muscles was recorded: splenius capitis (SP), trapezius (descending, TR), latissimus dorsi (LD), deltoid (anterior, DA), deltoid (posterior, DP), biceps brachii (BB), triceps brachii (long head, TB), erector spinæ (longissimus, L1 vertebra, ES), rectus abdominis (RA), abdominal external oblique (AE), gluteus medius (ME), gluteus maximus (MA), adductor longus (AL), 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). Around 50 gait cycles (49 ± 4) (Oliveira et al., 2014) were recorded after an accommodation period of maximum 60 s (White et al., 2002). Between the trials there was a break necessary to change shoes before and after running barefoot. The same randomized protocol was repeated after 15 min of rest for use in the (intraday) repeatability analysis, without removing the electrodes.

Foot Strike Patterns Assessment

For every trial, the FSP and the strike index were calculated using a validated algorithm based on the numerical analysis of foot pressure distribution (Santuz et al., 2016). As we previously suggested (Santuz et al., 2016), the FSPs have been grouped into two categories rather than three: RS and MFS (including MS and FS patterns). The strike index, as originally defined by Cavanagh and Lafortune (Cavanagh and Lafortune, 1980), was calculated as the distance from the heel to the center of pressure at impact relative to total foot length.

Spinal Motor Output Assessment

For the spinal motor output characterization, we mapped the 24 measured EMG activities onto the estimated rostrocaudal location of alpha-motoneurons (MNs) pools in the segments from the second cervical vertebra (C2) to the second sacral vertebra (S2) of the spinal cord (Ivanenko et al., 2006; La Scaleia et al., 2014). The segments T2, T3, and T4 have been excluded from the analysis since they do not innervate any of the considered muscles. The wireless EMG systems had a built-in band-pass filter (5–500 Hz, 3 dB/oct, 4th order). 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 normalized to the maximum activation recorded for each participant across all conditions (Karamanidis et al., 2004; Bizzi et al., 2008; Devarajan and Cheung, 2014). Each gait cycle was then time-normalized to 200 points (Cappellini et al., 2016), assigning 100 points to the stance and 100 points to the swing phase. The cervical segments (C2 to C8) mainly innervate upper limb and neck muscles. The thoracic segments (T1 to T12) connect to the trunk muscles, while the lumbar (L1 to L5) and sacral (S1 and S2) segments innervate the lower limb muscles. 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. (2005). This method shows the organization of the efferent MNs network directed to the muscles, assuming a common spinal topography among the investigated participants. Without accounting for size differences in MN pools at each spinal level, the motor output of each spinal segment S_j was estimated using the Equation (1.1) (La Scaleia et al., 2014):

$\begin{array}{ll}{S}_j=\frac{{\displaystyle {\sum }_{i=1}^{m_j}\left(\frac{k_{ji}}{n_i}\times EM{G}_i\right)}}{{\displaystyle {\sum }_{i=1}^{m_j}\left(\frac{k_{ji}}{n_i}\right)}}& (1.1)\end{array}$

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 normalized recorded EMG, specific for each participant and trial (Kendall et al., 2005; La Scaleia et al., 2014).

Modular Organization Assessment

The gait cycle breakdown was obtained from the pressure plate's raw data. Using a custom algorithm (Santuz et al., 2016), the touchdown was identified as the first non-zero pressure matrix after the last toe-off. The EMG signals were pre-processed using the filtering and normalization conditions reported above.

Muscle synergies data were extracted through a custom script (Santuz et al., 2017) (R v3.4.1, R Found. for Stat. Comp.) using the classical Gaussian NMF algorithm (Lee and Seung, 1999) from the first circa 50 gait cycles of each acquisition (Oliveira et al., 2014). EMG data were pre-processed using the same filtering conditions reported in the previous paragraph. The m = 24 time-dependent muscle activity vectors were grouped in an m × n matrix V, factorized such that V ≈ V_R = WH. V_R represents the new reconstructed matrix, which approximates the original matrix. The motor primitives matrix H (Dominici et al., 2011; Santuz et al., 2017) contained the time-dependent coefficients of the factorization with dimensions r × n, where r represents the number of synergies necessary to reconstruct the signal and n the number of data points (n = 200 number of cycles). The motor modules matrix W (Gizzi et al., 2011; Santuz et al., 2017) 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 (2.1, 2.2). The limit of convergence was reached when a change in the calculated R² between V and V_R was smaller than the 0.01% in the last 20 iterations (Cheung et al., 2005; Santuz et al., 2017). To choose the minimum number of synergies required to represent the original signals, the curve of R² values vs. synergies was fitted using a simple linear regression model, using all 10 synergies. The mean squared error 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⁻⁵ (Santuz et al., 2017).

$\{\begin{array}{c}H=H\frac{\left(W^{T}V\right)}{\left(W^{T}WH\right)}(2.1)\\ W=W\frac{\left(V{H}^T\right)}{\left(WH{W}^T\right)}(2.2)\end{array}$

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 activation (Santuz et al., 2017). When two or more fundamental synergies are blended into one, a combined synergy appears. Due to the lack of consent in the literature on how to interpret combined synergies, they were excluded from the analysis. An example of combined synergies is reported in Figure 1.

FIGURE 1

Figure 1. Example of two fundamental synergies combined into one. The histograms in the panels A,B represent the two fundamental sets of motor modules and the curves in the panels D,E the two respective primitives, with arbitrary x- and y-axis units. The combined motor modules and primitives are shown in panels C,F, respectively.

The fundamental synergies recognition was implemented using a custom learning algorithm based on a curve-fitting model (Santuz et al., 2017). The first implementation step consists in choosing some examples of single-peaked activation patterns, which might represent a fundamental primitive. The code is then provided with this training set and a search of similar shapes is done across the whole dataset of factorized curves. With a first iteration, the primitives that have a high similarity (R² > 0.95) with the ones present in the training set are added to it. The number of fundamental primitives is then selected by clustering similar motor modules. After updating the training set, the code starts the recognition across the entire dataset searching, synergy-by-synergy, for similar primitives (we found R² > 0.50 to be a good threshold). Non-recognized curves can then be visually inspected with an interactive routine or automatically identified as new fundamental or combined primitives. This approach, validated in a pilot study, can reproduce the results of a completely manual selection of the curves with a margin of error of ± 5%.

Metrics for Comparison of Curves

We evaluated the center of activation (CoA) and full width 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 center of mass of that circular distribution (Cappellini et al., 2016). The polar direction represented the gait cycle's phase, with angle 0 ≤ θ_t ≤ 2π. The following equations define the CoA:

$\{\begin{array}{c}A={\displaystyle \sum _{t=1}^p\left(\cos {\theta }_t\times P_t\right)}(3.1)\\ B={\displaystyle \sum _{t=1}^p\left(\sin {\theta }_t\times P_t\right)}(3.2)\end{array}$

$\begin{array}{ll}CoA=\text{arctan}(B/A)& (3.3)\end{array}$

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 minimum (Cappellini et al., 2016). For every trial, both parameters were calculated at each gait cycle and then averaged to proceed with the statistical analysis. A maximum of 50 gait cycles for each acquisition were selected for analysis. The CoA and FWHM were analyzed for stance and swing distinctively for spinal maps and over the whole gait cycle for the motor primitives.

Statistics

A two-way ANOVA with repeated measures, followed by a Tukey post-hoc analysis with false discovery rate p-value adjustment, was used to investigate CoA and FWHM between conditions. The same statistics was used for the motor modules, using the muscles and the conditions (shod vs. barefoot) as independent variables. To assess the similarities between the fundamental motor primitives of shod and barefoot running, we used the coefficient of determination R². We calculated the similarity values between the pairs of trials recorded during the same day (intraday repeatability) in shod and barefoot running. Then, we compared them with the similarity values between the two conditions (shod and barefoot running). Type A uncertainty was expressed as u_A = s/_√n. All the significance levels were set to α = 0.050 and the statistical analyses were conducted using R v3.4.1 (R Found. for Stat. Comp.).

Results

Foot Strike Patterns and Gait Parameters

Out of 20 participants, 14 (7 male, 7 female) transitioned from RS (shod) to MFS (barefoot). Three kept the MFS pattern in both conditions and three retained a RS pattern in both shod and barefoot running. The participants demonstrated significant (p < 0.001) differences in the average strike index, presenting values of 0.15 ± 0.17 in shod running and 0.53 ± 0.18 in barefoot running (with 0 denoting the most posterior and 1 the most anterior point of the shoe, see Figure 2). Also the average contact times of 301 ± 36 ms and 274 ± 32 ms as well as the average cadence (step frequency) of 162 ± 10 and 166 ± 11 steps/min were significantly different (p < 0.001) between shod and barefoot running, respectively. The mean values of the left and right VGRFs normalized to body weight were significantly lower in the barefoot condition (1.82 ± 0.20 for the shod and 1.75 ± 0.16 for the barefoot condition, p < 0.001). The impulse (mean values of left and right sides) was significantly lower in the barefoot condition (201 ± 39 N·s vs. 186 ± 35 N·s, p < 0.001), but the flight time was contrarily higher (70 ± 24 ms for shod and 89 ± 22 ms for barefoot running, p < 0.001).

FIGURE 2

Figure 2. Sagittal view of a typical rearfoot (A) and forefoot (B) strike patterns during shod and barefoot running, respectively. The strike index values extracted from plantar pressure distribution for these two representative cases are presented as well.

Spinal Motor Output

Figure 3 depicts the average spatiotemporal spinal motor output for shod and barefoot running. The two-way ANOVA identified statistically significant differences in the FWHM of the mapped EMG activities when comparing shod and barefoot running for both the stance (p = 0.018) and swing (p = 0.019) phase of the gait cycle (Table 1). The post-hoc analysis showed significantly lower FWHM in the barefoot condition of segment L4's spinal motor output, innervating the muscles ME, AL, FL, RF, VM, VL, ST, TA, and PL. The CoA was not significantly different between conditions in neither the stance (p = 0.107) or the swing (p = 0.091) phase (Table 1).

FIGURE 3

Figure 3. The average spatiotemporal spinal motor output is presented for shod and barefoot running, normalized in amplitude to the maximum of each segment. These curves have been obtained by mapping each of the 24 muscle activations onto the relative spinal segment (cervical from C2 to C8, thoracic from T1 to T12, lumbar from L1 to L5 and sacral from S1 and S2). The two level plots show the average alpha-motoneurons activity for each condition, giving additional information about the absolute activation level (normalization to the maximum of each condition). The stance and swing phases have been temporally normalized to the same amount of data points (100 each). Values are the means across all subjects and all trials.

TABLE 1

Table 1. Differences between shod and barefoot running in the center of activity (CoA) and full width at half maximum (FWHM) of the electromyographic activities mapped onto the estimated rostrocaudal location of the spinal cord (segments C2 to S2).

Modular Organization

The average number of recognized fundamental synergies during running was significantly different between the two conditions (3.9 ± 0.6 for shod and 3.6 ± 0.6 for barefoot running, p < 0.001). However, in both conditions, five fundamental activation patterns could be identified (Figures 4, 5). The five fundamental synergies extracted during both shod and barefoot running, were associated with temporally different phases of the gait cycle and ordered according to the timing of each motor primitive's global maximum (Figures 4, 5). The first synergy (peak at ~8% of the stance phase) functionally referred to the body weight acceptance, with a major involvement of knee extensors and plantarflexors. The second synergy (peak at ~27% of the stance phase) described the propulsion phase, to which the plantarflexors mainly contributed. The third synergy (peak at ~90% of the stance phase) was associated with the arm swing, when the upper body muscles played an important role. The fourth synergy (peak at ~22% of the swing phase) identified the early swing, showing contributions from upper body muscles, stabilizing muscles of the lower limb and the start of foot dorsiflexors activation. The fifth and last synergy (peak at ~72% of the swing phase) reflected the late swing and the landing preparation, highlighting the relevant contribution of knee flexors, foot dorsiflexors (in the shod condition) and plantarflexors (in the barefoot condition).

FIGURE 4

Figure 4. Motor modules and motor primitives for shod and barefoot running at the comfort speed. The motor modules are presented on a normalized y-axis base. For the motor primitives, the x-axis full scale represents one gait cycle (stance and swing normalized to the same amount of points and divided by a vertical line) and the y-axis the normalized amplitude. The motor modules represent the average normalized values across all the participants. The mean motor primitives are represented with a thick black line, while all the trials are denoted by thin gray lines.

FIGURE 5

Figure 5. Average motor modules and motor primitives of the five fundamental synergies for shod and barefoot running at the comfort speed. The motor modules are presented on a normalized y-axis base. For the motor primitives, the x-axis full scale represents one gait cycle (stance and swing normalized to the same amount of points and divided by a vertical line) and the y-axis the normalized amplitude. Asterisks denote significant differences between shod and barefoot running.

The motor primitives of the weight acceptance, propulsion and early swing synergies were significantly dissimilar (p = 0.023, 0.002, and <0.001, respectively; Figure 5, Table 2). The motor modules exhibited significant differences in the propulsion (p < 0.001), arm swing (p = 0.023) and late swing (p < 0.001) synergies (Figure 5). The muscles responsible for said changes where mainly the upper and lower leg muscles in the propulsion (higher contribution in the shod condition), the trunk muscles in the arm swing (higher contribution in the barefoot condition), the knee flexors and foot plantarflexors in the late swing synergy (higher contribution in the barefoot condition, Figure 5).

TABLE 2

Table 2. Similarities, indicated as $R_{\text{S},\text{B}}^2$, between the motor primitives of shod and barefoot running as mean of intraday repetitions.

The CoA of the motor primitives for all the synergies, except from the early swing one, moved significantly in time. The CoA values were lower in barefoot running (anticipated activation) for those synergies related completely or partially to the stance phase. For those synergies describing the only swing phase, the CoA values were instead bigger in the barefoot compared to the shod condition (Table 3). Further, we found a significant (p < 0.001) decrease in the FWHM of the propulsion primitives and an increase (p < 0.001) of the arm swing primitives in barefoot compared to shod running (Table 3).

TABLE 3

Table 3. Differences between shod and barefoot running in the center of activity (CoA) as well as in the relative full width at half maximum (FWHM) of motor primitives.

Discussion

In this study, we analyzed the modularity of the neuromuscular control of shod and barefoot running. We hypothesized a different modular organization of motion mainly due to the presence or absence of shoes in the two conditions. We found that the motor primitives (or fundamental activation patterns) were generally shifted earlier in time during the stance-related phases and later in the swing-related ones. The motor primitives were found to be significantly wider in the arm swing phase but not in the propulsion, where the basic activation was significantly of shorter duration (peculiarity confirmed by the analysis of the spinal motor output). Moreover, the motor modules (or muscle weightings) demonstrated analogous organization with some significant differences in the propulsion, arm swing and late swing synergies.

The cadence and the strike index significantly increased when changing from shod to barefoot running. Contact times and VGRFs decreased accordingly in the barefoot compared to the shod condition. These results agree with previous studies (De Wit et al., 2000; Lieberman et al., 2010) on the comparison of shod and barefoot running. It is well known that the gear ratios of the ankle joint muscles [i.e., the ratio between the ground reaction force and the muscle force moment arms Carrier et al., 1994] do not only vary through the running stance phase (Carrier et al., 1994), but also when switching from the shod to the barefoot condition (Braunstein et al., 2010). In the last 20% of the stance phase the gear ratio at the ankle joint is lower during barefoot compared to shod running (Braunstein et al., 2010). Lower gear ratios at the ankle joint decrease the contact time while running (Lee and Piazza, 2009) and provide an explanation for the shorter contact times found during barefoot running. Further, a lower gear ratio at the ankle joint induces a reduction in the potential of the plantarflexors to generate efficient muscle force due to the force-velocity relationship (Carrier et al., 1994). In inexperienced runners, this may initiate a dynamic instability in the whole system (including the upper body), requiring stabilization achieved through feedback- as well as predictive-based motor control. We recently found a significant decrease in the dynamic stability of running by switching from shod to barefoot (Ekizos et al., 2017). Moreover, it has been reported that the intrinsic foot muscles show higher absolute activation levels during stance in shod compared to barefoot running (Kelly et al., 2016). This difference produces an alteration in the longitudinal arch compression during the stance phase, leading to higher recoil capabilities in barefoot running (Kelly et al., 2016). This increase in the capacity of the foot to store and return energy is likely an odd feature for the unexperienced barefoot runner and might be another mechanism driving the system to an increased instability.

These very same factors (i.e. different gear ratios, dynamic stability and foot's recoil capabilities) could as well partly explain the differences we found in the duration of the motor primitives. First of all, the reduction in duration of the propulsion-related primitive might be a direct consequence of the lower gear ratios and, possibly, of the increased energy storage and return capabilities when running barefoot. However, this does not explain the increase in the duration of the motor primitives in the arm swing synergy. It has been recently shown that the FWHM of EMG activity undergoes, during gait, a systematic decrease with age in typically developing children (Cappellini et al., 2016). Conversely, very limited age-related changes appear in children affected by cerebral palsy. Moreover, cerebral palsy and typically developing children show a comparable structure of motor modules (Cappellini et al., 2016). Analogously, a widening of the motor primitives can be found in adult patients with cerebellar ataxia and in healthy adults walking on a narrow beam and on slippery ground (Martino et al., 2015). This consolidation of the motor output, promoted by learning and impaired by pathology, might reflect the system's need of adding fail-safe robustness to cope with previously unexperienced running conditions (e.g., the absence of footwear).

Concerning motor modules, significant differences were found in the propulsion, arms swing and late swing synergies. The modules of the propulsion phase indicated that upper leg muscles and, most importantly, foot plantarflexors mainly contributed to the inequality. The relative contributions of these muscles were lower in the barefoot condition, indicating a higher specificity of the muscles more important for the propulsion. During arm swing, the TR, RA, and PL muscles were found to be significantly responsible for the identified changes. The relative contribution of TR and PL was higher in barefoot compared to shod running, while the contrary emerged for the RA. However, the intrinsic variability of this synergy's patterns is high and the EMG activities low compared to other gait cycle phases. Therefore, small adaptations in the strategy might translate in statistical differences. The changes in FSP are the cause for the alteration of the motor modules of foot dorsiflexors and plantarflexors in the late swing synergy. In agreement with the prediction based on one of our earlier studies (Santuz et al., 2016), 70% of the participants changed FSP without undergoing a specific training intervention when switching from shod to barefoot running. Most of the times participants automatically switched from RS (shod) to MFS (barefoot). In some cases, participants changed FSP after a few steps, reportedly due to the discomfort of striking the ground with the bare rearfoot. Specifically, 14 out of 20 participants transitioned from RS (shod) to MFS (barefoot). It is well known that the muscles TA and GM and GL play an important role in the final part of the swing phase, just before touchdown (Komi, 1984). In RS patterns, the TA has the twofold task of dorsiflexing the foot to prepare it for the strike and to control the plantarflexion immediately after the touchdown (von Tscharner et al., 2003). In MFS patterns, given the substantial impact loads at contact during running, a preactivation happens right before the strike and the subsequent activation in the early stages of the stance phase (Komi, 1992). Looking at the late swing synergy, it is evident that the TA contributed more in shod running, an activity that mostly involves a RS pattern. In contrast, the average pattern in barefoot running was a MFS, where the preactivation of GM and GL is predominant. These considerations might as well be extended to overground running, since it has been recently shown that treadmill and overground running share similar motor modules with minimal temporal shifts in the motor primitives (Oliveira et al., 2016).

We cannot exclude that habitual barefoot runners might be able to compensate for the differences in the modular organization of muscle activation found in our participants (which were all inexperienced barefoot runners). Although some effects of barefoot running habituation on FSP can be expected, we argue that the main alterations in the motor modules would remain visible in habitual barefoot runners. This mainly because of the predictable changes in the EMG activity (Komi, 1984, 1992; von Tscharner et al., 2003) and, consequently, in the motor modules associated to the kinematic and kinetic alterations induced by a MFS compared to a RS. Concerning motor primitives, however, we suggest that a training intervention focused on the practice of barefoot running might lead to an improvement in the accuracy of motor commands' timing, thus reducing the FWHM of those primitives that here appear wider. Given the characteristics of barefoot running that we discussed above, however, we do not expect that a retraining program would be able to affect consistently the propulsion motor primitive.

With this study using the muscle synergies concept, we bring new insights in the modular organization of shod and barefoot running. Investigating the differences between the synergies, we could confirm that although in both shod and barefoot running five fundamental synergies are enough to describe the running task, a dissimilarity exists in the modular organization of movement. Moreover, we found an increase in the FWHM of the motor primitives of the arm swing synergy: a possible indication of weak motor learning (Cappellini et al., 2016). These findings suggest a reorganization of the motor output possibly due to the nervous system's effort to cope with the biomechanical specificity of barefoot running. This specificity might be explained by a lower ankle gear ratio (Carrier et al., 1994), different FSP (Komi, 1984, 1992; von Tscharner et al., 2003) and increased instability (Ekizos et al., 2017) created by the absence of shoes. The results indicate a possible reorganization of movement when task's complexity either increases or is not well mastered. The required adjustments seem to go in the direction of an improved robustness of motor output guaranteed by longer activation patterns applied on similar muscle modules, showing some adaptability in such a task-specific structure as the muscle synergies.

Author Contributions

Conceptualization: AS, LJ, AE, and AA; Methodology: AS, LJ, and AA; Investigation: AS, LJ, and AE; Formal Analysis: AS and LJ; Writing – Original Draft: AS; Writing – Review and Editing: AS, LJ, AE, VB, and AA; Visualization: AS; Supervision: VB and AA.

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

The authors are grateful to the participants that always showed great commitment and interest during the measurements and to DAAD for the financial support during the publication process.

References

Bernstein, N. A. (1967). The Co-ordination and Regulation of Movements. Oxford: Pergamon Press Ltd.

Google Scholar

Bizzi, E., and Cheung, V. C.-K. (2013). The neural origin of muscle synergies. Front. Comput. Neurosci. 7:51. doi: 10.3389/fncom.2013.00051

PubMed Abstract | CrossRef Full Text | Google Scholar

Bizzi, E., Cheung, V. C.-K., D'Avella, A., Saltiel, P., and Tresch, M. C. (2008). Combining modules for movement. Brain Res. Rev. 57, 125–133. doi: 10.1016/j.brainresrev.2007.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Bramble, D. M., and Lieberman, D. E. (2004). Endurance running and the evolution of Homo. Nature 432, 345–352. doi: 10.1038/nature03052

PubMed Abstract | CrossRef Full Text | Google Scholar

Braunstein, B., Arampatzis, A., Eysel, P., and Brüggemann, G.-P. (2010). Footwear affects the gearing at the ankle and knee joints during running. J. Biomech. 43, 2120–2125. doi: 10.1016/j.jbiomech.2010.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Cappellini, G., Ivanenko, Y. P., Dominici, N., Poppele, R. E., and Lacquaniti, F. (2010). Migration of motor pool activity in the spinal cord reflects body mechanics in human locomotion. J. Neurophysiol. 104, 3064–3073. doi: 10.1152/jn.00318.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Cappellini, G., Ivanenko, Y. P., Martino, G., MacLellan, M. J., Sacco, A., Morelli, D., et al. (2016). Immature spinal locomotor output in children with cerebral palsy. Front. Physiol. 7:478. doi: 10.3389/fphys.2016.00478

PubMed Abstract | CrossRef Full Text | Google Scholar

Cappellini, G., Ivanenko, Y. P., Poppele, R. E., and Lacquaniti, F. (2006). Motor patterns in human walking and running. J. Neurophysiol. 95, 3426–3437. doi: 10.1152/jn.00081.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Carrier, D. R., Heglund, N. C., and Earls, K. D. (1994). Variable gearing during locomotion in the human musculoskeletal system. Science 265, 651–653. doi: 10.1126/science.8036513

PubMed Abstract | CrossRef Full Text | Google Scholar

Cavanagh, P. R., and Lafortune, M. A. (1980). Ground reaction forces in distance running. J. Biomech. 13, 397–406.

PubMed Abstract | Google Scholar

Cheung, V. C.-K., D'Avella, A., Tresch, M. C., and Bizzi, E. (2005). Central and sensory contributions to the activation and organization of muscle synergies during natural motor behaviors. J. Neurosci. 25, 6419–6434. doi: 10.1523/JNEUROSCI.4904-04.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Devarajan, K., and Cheung, V. C.-K. (2014). On nonnegative matrix factorization algorithms for signal-dependent noise with application to electromyography data. Neural Comput. 26, 1128–1168. doi: 10.1162/NECO_a_00576

PubMed Abstract | CrossRef Full Text | Google Scholar

De Wit, B., De Clercq, D., and Aerts, P. (2000). Biomechanical analysis of the stance phase during barefoot and shod running. J. Biomech. 33, 269–278. doi: 10.1016/S0021-9290(99)00192-X

PubMed Abstract | CrossRef Full Text | Google Scholar

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. doi: 10.1126/science.1210617

PubMed Abstract | CrossRef Full Text | Google Scholar

Ekizos, A., Santuz, A., and Arampatzis, A. (2017). Transition from shod to barefoot alters dynamic stability during running. Gait Posture 56, 31–36. doi: 10.1016/j.gaitpost.2017.04.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Gizzi, L., Nielsen, J. F., Felici, F., Ivanenko, Y. P., and Farina, D. (2011). Impulses of activation but not motor modules are preserved in the locomotion of subacute stroke patients. J. Neurophysiol. 106, 202–210. doi: 10.1152/jn.00727.2010

CrossRef Full Text | Google Scholar

Hägglund, M., Dougherty, K. J., Borgius, L., Itohara, S., Iwasato, T., and Kiehn, O. (2013). Optogenetic dissection reveals multiple rhythmogenic modules underlying locomotion. Proc. Natl. Acad. Sci. U.S.A. 110, 11589–11594. doi: 10.1073/pnas.1304365110

PubMed Abstract | CrossRef Full Text | Google Scholar

Hart, C. B., and Giszter, S. F. (2004). Modular premotor drives and unit bursts as primitives for frog motor behaviors. J. Neurosci. 24, 5269–5282. doi: 10.1523/JNEUROSCI.5626-03.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Hasegawa, H., Yamauchi, T., and Kraemer, W. J. (2007). Foot strike patterns of runners at the 15-km point during an elite-level half marathon. J. Strength Cond. Res. 21, 888–893. doi: 10.1519/R-22096.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanenko, Y. P., Cappellini, G., Poppele, R. E., and Lacquaniti, F. (2008). Spatiotemporal organization of α-motoneuron activity in the human spinal cord during different gaits and gait transitions. Eur. J. Neurosci. 27, 3351–3368. doi: 10.1111/j.1460-9568.2008.06289.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanenko, Y. P., Poppele, R. E., and Lacquaniti, F. (2006). Spinal cord maps of spatiotemporal alpha-motoneuron activation in humans walking at different speeds. J. Neurophysiol. 95, 602–618. doi: 10.1152/jn.00767.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Karamanidis, K., Arampatzis, A., and Brüggemann, G.-P. (2004). Reproducibility of electromyography and ground reaction force during various running techniques. Gait Posture 19, 115–123. doi: 10.1016/S0966-6362(03)00040-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Kelly, L. A., Lichtwark, G. A., Farris, D. J., and Cresswell, A. (2016). Shoes alter the spring-like function of the human foot during running. J. R. Soc. Interface 13:20160174. doi: 10.1098/rsif.2016.0174

PubMed Abstract | CrossRef Full Text | Google Scholar

Kendall, F. P., McCreary, E. K., Provance, P. G., Rodgers, M. M., and Romani, W. A. (2005). Muscles: Testing and Function, with Posture and Pain, 5th Edn. Philadelphia, PA: Lippincott Williams & Wilkins.

Komi, P. V. (1984). Physiological and biomechanical correlates of muscle function: effects of muscle structure and stretch-shortening cycle on force and speed. Exerc. Sport Sci. Rev. 12, 81–121.

PubMed Abstract | Google Scholar

Komi, P. V. (1992). “Strength and power in sport - Stretch–Shortening Cycle,” in The Encyclopaedia of Sports Medicine, ed P. V. Komi (Oxford: Blackwell Science Ltd.), 184–202.

La Scaleia, V., Ivanenko, Y. P., Zelik, K. E., and Lacquaniti, F. (2014). Spinal motor outputs during step-to-step transitions of diverse human gaits. Front. Hum. Neurosci. 8:305. doi: 10.3389/fnhum.2014.00305

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, D. D., and Seung, H. S. (1999). Learning the parts of objects by non-negative matrix factorization. Nature 401, 788–791. doi: 10.1038/44565

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S. S. M., and Piazza, S. J. (2009). Built for speed: musculoskeletal structure and sprinting ability. J. Exp. Biol. 212, 3700–3707. doi: 10.1242/jeb.031096

PubMed Abstract | CrossRef Full Text | Google Scholar

Lieberman, D. E., Venkadesan, M., Werbel, W. A., Daoud, A. I., D'Andrea, S., Davis, I. S., et al. (2010). Foot strike patterns and collision forces in habitually barefoot vs. shod runners. Nature 463, 531–535. doi: 10.1038/nature08723

CrossRef Full Text | Google Scholar

Martino, G., Ivanenko, Y. P., D'Avella, A., Serrao, M., Ranavolo, A., Draicchio, F., et al. (2015). Neuromuscular adjustments of gait associated with unstable conditions. J. Neurophysiol. 114, 2867–2882. doi: 10.1152/jn.00029.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Monda, V., Valenzano, A., Moscatelli, F., Salerno, M., Sessa, F., Triggiani, A. I., et al. (2017). Primary motor cortex excitability in karate athletes: a transcranial magnetic stimulation study. Front. Physiol. 8:695. doi: 10.3389/fphys.2017.00695

PubMed Abstract | CrossRef Full Text | Google Scholar

Moscatelli, F., Messina, G., Valenzano, A., Monda, V., Viggiano, A., Messina, A., et al. (2016a). Functional assessment of corticospinal system excitability in karate athletes. PLoS ONE 11:e0159846. doi: 10.1371/journal.pone.0159846

CrossRef Full Text | Google Scholar

Moscatelli, F., Messina, G., Valenzano, A., Petito, A., Triggiani, A. I., Messina, A., et al. (2016b). Differences in corticospinal system activity and reaction response between karate athletes and non-athletes. Neurol. Sci. 37, 1947–1953. doi: 10.1007/s10072-016-2693-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Mussa-Ivaldi, F. A., Giszter, S. F., and Bizzi, E. (1994). Linear combinations of primitives in vertebrate motor control. Proc. Natl. Acad. Sci. U.S.A. 91, 7534–7538. doi: 10.1073/pnas.91.16.7534

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliveira, A. S. C., Gizzi, L., Farina, D., and Kersting, U. G. (2014). Motor modules of human locomotion: influence of EMG averaging, concatenation, and number of step cycles. Front. Hum. Neurosci. 8:335. doi: 10.3389/fnhum.2014.00335

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliveira, A. S. C., Gizzi, L., Ketabi, S., Farina, D., and Kersting, U. G. (2016). Modular control of treadmill vs. overground running. PLoS ONE 11:e0153307. doi: 10.1371/journal.pone.0153307

PubMed Abstract | CrossRef Full Text | Google Scholar

Overduin, S. A., D'Avella, A., Roh, J., Carmena, J. M., and Bizzi, E. (2015). Representation of muscle synergies in the primate brain. J. Neurosci. 35, 12615–12624. doi: 10.1523/JNEUROSCI.4302-14.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Roh, J., Cheung, V. C.-K., and Bizzi, E. (2011). Modules in the brain stem and spinal cord underlying motor behaviors. J. Neurophysiol. 106, 1363–1378. doi: 10.1152/jn.00842.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Santuz, A., Ekizos, A., and Arampatzis, A. (2016). A pressure plate-based method for the automatic assessment of foot strike patterns during running. Ann. Biomed. Eng. 44, 1646–1655. doi: 10.1007/s10439-015-1484-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Santuz, A., Ekizos, A., Janshen, L., Baltzopoulos, V., and Arampatzis, A. (2017). On the methodological implications of extracting muscle synergies from human locomotion. Int. J. Neural Syst. 27:1750007. doi: 10.1142/S0129065717500071

PubMed Abstract | CrossRef Full Text | Google Scholar

Tresch, M. C., Saltiel, P., D'Avella, A., and Bizzi, E. (2002). Coordination and localization in spinal motor systems. Brain Res. Rev. 40, 66–79. doi: 10.1016/S0165-0173(02)00189-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Treutwein, B. (1995). Adaptive psychophysical procedures. Vision Res. 35, 2503–2522. doi: 10.1016/0042-6989(95)00016-S

PubMed Abstract | CrossRef Full Text | Google Scholar

von Tscharner, V., Goepfert, B., and Nigg, B. M. (2003). Changes in EMG signals for the muscle tibialis anterior while running barefoot or with shoes resolved by non-linearly scaled wavelets. J. Biomech. 36, 1169–1176. doi: 10.1016/S0021-9290(03)00086-1

PubMed Abstract | CrossRef Full Text | Google Scholar

White, S. C., Gilchrist, L. A., and Christina, K. A. (2002). Within-day accommodation effects on vertical reaction forces for treadmill running. J. Appl. Biomech. 18, 74–82. doi: 10.1123/jab.18.1.74

CrossRef Full Text | Google Scholar

Yakovenko, S., Krouchev, N., and Drew, T. (2011). Sequential activation of motor cortical neurons contributes to intralimb coordination during reaching in the cat by modulating muscle synergies. J. Neurophysiol. 105, 388–409. doi: 10.1152/jn.00469.2010

PubMed Abstract | CrossRef Full Text | Google Scholar