## Materials and Methods
### Synergy Extraction and Data Analyses  
All EMG and kinematic data were processed using custom MATLAB scripts (MATLAB R2014b, MathWorks(r), Natick, MA, USA). The EMG signals were high-pass filtered with a 1 st order Butterworth filter with a cut-off at 20 Hz, full-wave rectified and smoothed with a 0.1-s moving average window [10; 38]. EMG signals from forward and lateral stepping were normalized to the respective maximum activation obtained over all trials performed in the congruent direction so that activation could not exceed 100% [39; 40]. The EMG signals were time-synchronized with the kinematically defined start and end points and subsequently normalized over time to define 0%-100% of the step cycle. Finally, because EMG data were only collected during step ascent and therefore represented intervals, rather than continuous activation patterns (as would be the case during gait trials), signals were averaged over the three repetitions performed in each condition. The choice to average signals rather than concatenating them was made in order to obtain the best reconstruction quality for our relatively short intervals, at the risk of losing information on step-to-step variability [40].  
Muscle synergies were extracted from the individual average EMG data matrix using non-negative matrix factorization (NNMF). NNMF calculates muscle synergies (W) and their relative temporal activation patterns (C), resulting in muscle activations being represented as W \(\times\) C + e. W represents the relative muscle co-activation, defined as the relative weight of each muscle per synergy, and is constructed as an \(m\times n\) matrix where \(m\) is the total number of muscles and \(n\) is the selected number of synergies. C represents the temporal activation patterns and is constructed as an \(n\times t\) matrix where \(t\) represents the number of data points over normalized time (100 per individual trial) and e is the residual error matrix [22; 41]. The algorithm was repeated 1000 times for each subject to avoid local minima. The appropriate number of synergies was defined using two criteria. First, using an iterative process where the number of synergies varied between 1 and 10, the minimum number of synergies was selected based on the number required to reach \(\geq\) 85% of group-averaged variance accounted for (VAF). As an additional local criterion, synergies had to account for \(\geq\)75% VAF for each individual muscle [42]. This double criterion approach was selected in order to adequately reproduce relevant features of the synergy compositions. VAF was defined as the uncentered Pearson correlation coefficient between W \(\times\) C and the EMG amplitude time series. To compare spatio-temporal characteristics, individual synergies obtained from different subjects were pooled and matched based on the correlation of their structure (muscle weightings in each synergy of W) using a custom cluster analysis algorithm [43]. If a synergy showed equal correlation to more clusters, that synergy remained in the pool it was initially assigned to. Each synergy of W and C was subsequently averaged over all participants in that age group. For comparisons between age groups, the group-averaged synergies were also matched based on their structure using cluster analysis. Finally, we computed time-averaged standard deviations of the synergy activation patterns, with a fixed number of synergies, to assess if age-related differences in group-averaged VAF could be affected by between-subject variability of temporal activation patterns.