## Materials and Methods
### Patient Multi-Domain Assessment  
#### 2.3.1 Kinematic and EMG Assessments  
At T0 and T1, an instrumental data recording was performed with 6-infra-red cameras and a BTS Smart-D system, while patients performed HTMM movements.  
Upper-limb kinematics was obtained by tracking the position of passive markers placed on D5 and C7 vertebas, acromion, lateral elbow epicondyle and styloid process of the ulna. Electromyography was recorded from eight muscles, namely, upper trapezius, pectoralis major, anterior deltoid, middle deltoid, posterior deltoid, triceps brachii caput medialis, biceps brachii caput longus, and brachioradialis. The sampling rate was 140 Hz and 1000 Hz for kinematics and EMG, respectively. Non-assisted movements were chosen in order to evaluate the range of motion, to characterize the physical status of the patient. Moreover, this choice allowed to assess the voluntary strategy (spatial and temporal) of muscle recruiting in a realistic task.  
Kinematics Analysis  
The 3D coordinates of the markers were filtered with a low-pass, 3 rd order Butterworth filter, with a cut-off frequency set at 6 Hz. They were used to compute movement phases, and only the forward phase, when the subject drew the hand toward the mouth, were analyzed. The following parameters were chosen as a selection of relevant measures based on kinematics:  
* mean duration of the forward phase;
* shoulder flexion angle, projected into the sagittal plane;
* elbow flexion angle;
* normalized jerk, as a measure of smoothness, computed according to the formula presented in [PERSON] et al. [35].  
Details of the kinematic protocol are available in a previous paper [21].  
EMG Analysis and Synergies Extraction  
Consistently with kinematic analysis, only forward phases were considered for EMG processing. Each signal from the 8 channels was high-pass filtered with a cut off frequency of 50 Hz, rectified and low-pass filtered at 10 Hz by means of a Butterworth filter (5 th order). Removing low frequencies avoids motion artifacts, while the low pass filter allowed to remove noise and obtain EMG envelopes. Afterwards, all envelopes were aligned and resampled (1000 samples), so that they could be directly compared.  
Defining the number of patients in the study as \(nPt=5\), the number of experimental conditions (pre and post therapy) as \(c=2\), the resampling size of each phase as \(rate=1000\), and the number of considered muscles as \(nm=8\). The EMG envelopes from each experimental condition were pooled together in a group of \(nPt\times c\) aggregated matrices (each one having size \(rate\times nm\)) and synergies were extracted using the non-negative matrix factorization (NMF) algorithm for the decomposition of each of the aggregated matrices [36]. The NMF outputs two matrices: one represents time-invariant, spatial synergies, while the other represents time-variant activations signals for each synergy. The factorization order \(r\), used as input for the NMF, was selected in the range from 1 to 8. For each \(r\), the NMF algorithm was applied 1000 times in order to avoid local minima and the repetition accounting for the higher variance of the signal was chosen as the representative of order \(r\). The order representing spatial synergies and temporal components was chosen as the minimum explaining at least 0.80 of the variance of the signal. Further synergies were added to the dataset only if increasing \(r\) increased the overall explained variance of at least 0.05.