## 2 Materials and Methods
### Data Analysis  
**Gait event detection and walking speed**. The markers trajectories were sampled at 100 Hz. The x, y, and z components were smoothed with a fourth-order Savitzky-Golay filter with a cut-off frequency of 9 Hz [27], which was used to obtain also the subsequent derivative terms. The gait cycle events, i.e., the heel strike (HS) for initial stance and the toe off (TO) for initial swing, were defined based on the position and speed of the foot. Specifically, we used the algorithm proposed in [28], considering the markers located in correspondence of the fifth metatarsal phalangeal joint and of the lateral malleolus instead of the toe and heel markers. For both normal and exoskeleton-assisted gait, we validated the algorithm by comparing it with the timing of the same gait events, as determined with the use of the force plates, in a preliminary testing and in all the testing trials (15% for Ekso strides and 40% normal strides) where the kinetic data were available. We found a difference <4% for all gait conditions, as reported by [29; 20]. These gait events were used to compute the stride duration or gait cycle duration i.e., the time elapsed between two consecutive HS of the same leg (left and right strides). The swing durations were measured as the time between the TO and HS when the considered foot is touching the ground. Here, we reported the percentage of the swing phase with respect to the stride duration. The mean speed of the subjects during assisted and unassisted gait was computed from the speed of the marker on C7 during each stride. **Single muscle activations. ** EMG data were recorded at 1 kHz, processed using a ninth-order band pass Butterworth filter between 20-450 Hz [30; 31]. The filtered data were rectified, and we applied a low-pass fourth-order Butterworth filter with 4-Hz cut-off frequency to obtain their envelope. The EMG envelopes were segmented in correspondence of each gait cycle (from each HS to the following one). Then, we time-interpolated each gait cycle over a time base with 101 points [20]. The position of the electrodes did not change during the entire acquisition (see Figure 1b,c), allowing for direct comparisons of each muscle among all the conditions. To directly compare and average the EMG data from different subjects and legs, we normalized each muscle signal for its maximum value computed over all the acquisitions under the different conditions that we considered in the analyses (e.g., for each subject/leg/muscle = 6 strides \(\times\) 2 aids \(\times\) 5 Ekso - modalities + 3 strides \(\times\) 2 legs for the normal walking). We also verified that a different normalization (e.g., for the median value [32; 33]) did not change the main results we obtained. **Muscle synergy analysis. ** For each participant and walking modality, we extracted muscle synergies from a matrix obtained by the normalized EMG envelopes, by using the non-negative matrix factorization (NNMF) algorithm [34; 35]. The NNMF algorithm decomposes the normalized EMG envelopes in a defined number of positive components, or muscle synergies, each composed of weight coefficients (W) and activation profiles (H): the first (W) represents the participation of each muscle in each synergy [36], and the second (H) represents the timing of activity of each muscle synergy. The NNMF was implemented using the multiplicative update rules [37] with W and H initialized with random non-negative values. The implementation of the NNMF algorithm is based on the minimization of the \"Euclidian distance\" between the muscle synergies and the combination of W and H.  
In order to minimize the possibility that the iterative algorithm converges to a local and not a global minimum, the extraction was repeated 50 times, and we selected the solution explaining the highest overall amount of variance. We obtained eight sets of muscle synergies for each subject in each walking modality.