A multimodal approach to capture post-stroke temporal dynamics of recovery

Objective. Several training programs have been developed in the past to restore motor functions after stroke. Their efficacy strongly relies on the possibility to assess individual levels of impairment and recovery rate. However, commonly used clinical scales rely mainly on subjective functional assessments and are not able to provide a complete description of patients’ neuro-biomechanical status. Therefore, current clinical tests should be integrated with specific physiological measurements, i.e. kinematic, muscular, and brain activities, to obtain a deep understanding of patients’ condition and of its evolution through time and rehabilitative intervention. Approach. We proposed a multivariate approach for motor control assessment that simultaneously measures kinematic, muscle and brain activity and combines the main physiological variables extracted from these signals using principal component analysis (PCA). We tested it in a group of six sub-acute stroke subjects evaluated extensively before and after a four-week training, using an upper-limb exoskeleton while performing a reaching task, along with brain and muscle measurements. Main results. After training, all subjects exhibited clinical improvements correlating with changes in kinematics, muscle synergies, and spinal maps. Movements were smoother and faster, while muscle synergies increased in numbers and became more similar to those of the healthy controls. These findings were coupled with changes in cortical oscillations depicted by EEG-topographies. When combining these physiological variables using PCA, we found that (i) patients’ kinematic and spinal maps parameters improved continuously during the four assessments; (ii) muscle coordination augmented mainly during treatment, and (iii) brain oscillations recovered mostly pre-treatment as a consequence of short-term subacute changes. Significance. Although these are preliminary results, the proposed approach has the potential of identifying significant biomarkers for patient stratification as well as for the design of more effective rehabilitation protocols.


Introduction
Stroke is the leading cause of adult long-term disability in Western societies. More than 1.5 million people are affected every year in Europe [1]. Even though acute stroke care and intensive rehabilitation are improving, two-thirds of chronic stroke survivors have to cope with persisting neurologic deficits, and only 20% of them are able to go back to their normal professional and private life [2].
The most common impairments in the acute and chronic stages are cognitive conditions and motor deficits contralateral to the affected brain hemisphere [3]. A profound neuromuscular reorganization occurs after stroke [4,5]. The affected limb is typically characterized by spasticity [6], stereotyped movement patterns, mainly caused by abnormal muscle co-activation and an enlarged activity of the antagonist muscles [7,8], which result in a reduced range of motion against gravity [5], and, thus, to a limited workspace in three-dimensional reaching movements [9,10].
A key factor in stroke recovery is the intensity of training, especially in the acute phase, to enhance functional restoration and prevent inactivity-related complications [11]. Yet because of the patient-specific clinical picture, treatment programs might vary in duration, intensity, and frequency [12]. Therefore, the success of the rehabilitation processes depends on the ability of the clinician to discern the individual levels of impairment and responses to treatment with simple, robust, and effective methods.
Currently, patients are evaluated mainly using clinical scales, with Fugl-Meyer being one of the most adopted measures of motor impairment after stroke. Yet, the precision of these clinical tests are limited by inter-rater and intra-rater reliability [13][14][15], as well as by floor and ceiling effects [16,17]. Moreover, some of them require a consider amount of time to be administered. Clinical scales should, therefore, be integrated with targeted neuro-biomechanical assessments, in order to provide a more detailed description of the patients' clinical status.
Many instrumental approaches that investigate different domains of the hierarchical organization of the neuromusculoskeletal system can be employed for this purpose, including measures of kinematics, muscular, and brain activities.
In a recent work, Thrane and colleagues coupled standard clinical tests with upper extremity kinematics to reduce the ceiling effects of the Fugl-Meyer Assessment [18]. They showed that post-stroke participants with near or fully recovered sensorimotor function still showed deficits in movement kinematics that were not captured by the clinical assessment. These kinematic measures provide a detailed and quantitative description of motor behaviors.
However, the neural deficits may be masked at the kinematic level by compensatory strategies and similar movements may be produced through different neuromuscular mechanisms. Therefore, muscular control strategies should be also considered in the patients' clinical picture. Indeed, several behavioral studies on animals and humans [19][20][21][22] have showed that muscle synergies extracted from the factorization analysis of electromyographic signals (EMG) can reveal underlying patterns in muscle activity that may reflect different levels of neural functions and their integrity has been proposed as a physiological marker of cortical damage. For instance, Cheung et al observed that in severe stroke patients there is a lack of preservation of muscle synergies in the affected side [23], and their level of preservation correlates with the level of motor impairment [24]. Finally, recordings of neural activity can also inform about the central nervous system reorganization after brain damage. Indeed, previous electrophysiological studies in post-stroke patients showed that compromised cortical areas are characterized by augmented slow potentials, which are predictive of long-term poststroke outcome [25,26].
Despite these interesting and promising results, each of the proposed approach investigates a specific domain of the neuro-musculoskeletal system. Therefore, each methodology provides a detailed but sectorial assessment. Conversely, merging all the domains may provide a comprehensive framework for a more complete and quantitative patient profiling [27,28].
To this end, the main goal of this study was to develop a multivariate analysis method to couple clinical evaluations with multimodal instrumental evaluations in order to provide a deeper characterization of the neuro-biomechanical status of stroke patients undergoing different rehabilitation protocols. We first introduced an exhaustive and diverse set of measures extracted from different sources, including motor performance as well as muscle activity during a motor task and brain activity at rest. Then, we presented a methodological approach combining and integrating these parameters, which allowed to identify neuro-biomechanical features modifications at different time points during intervention and recovery. We tested this methodology on a small cohort of stroke subjects that went through a period of intense motor training using a 3D reaching task. We believe that the features extracted with our comprehensive multivariate analysis can increase the understanding of the mechanisms underlying motor impairments and recovery, and, additionally, they can potentially help design more effective rehabilitative interventions and monitor the progress of the disease as well as the effects of rehabilitative treatments [29,30].

Subjects
Six stroke subjects (4 females, age 68 ± 18 yo) between 2 and 6 weeks from the occurrence of the stroke lesion, all with right hemiplegia and at least 10 • of residual motion in shoulder and elbow joints (details in table 1), and six healthy subjects (4 females, age 58 ± 16 yo) were enrolled in the study. All subjects were right-handed. The healthy subjects did not present any evidence or known history of skeletal or neurological diseases, and they exhibited intact joint range of motion and muscle strength.
The study was carried out in the Neurorehabilitation Unit of the University Hospital of Geneva (HUG), Switzerland and of the University Hospital of Pisa (Cisanello hospital), Italy. It was approved by the Commission Cantonale d'Ethique de la Recherche (CCER) de Genève, Switzerland, and by the Comitato Etico Area Vasta Nord Ovet (CEAVNO) in Pisa. The recordings were carried out in agreement with the Declaration of Helsinki and Good Clinical Practice norms. The study is registered with the number NCT02770300 in ClinicalTrials.gov. The participants were informed of the procedure and they signed an informed consent, which included the consent for the use of all data collected during the experiment in scientific publications.

Experimental set-up and procedure
The experimental protocol for the stroke patients consisted of four sessions of clinical and robotic assessment interleaved by four weeks of experimental training-three sessions of 30 min per weekproposed in addition to the habitual physical rehabilitative treatment of the patient. The stroke population completed two assessment sessions before (A1 and A2) and two after (A3 and A4) the experimental training. The initial assessment sessions were completed two weeks (A1, baseline) and one week (A2) before the beginning of the training. A2 was done in order to estimate the rate of the changes due to spontaneous recovery and inpatient therapy, since the patients were in the subacute phase of the stroke. The final two assessment sessions were completed one week after the end of the training (A3), in order to evaluate the effects of the rehabilitation protocol, and one month after the end of the training (A4) to evaluate the retention of the changes induced by the experimental training (figure 1(a)).
The enrolled stroke participants were receiving different therapies during the period of experimental training, including extra sessions of conventional therapy without the use of robotic devices; or standard robot-assisted rehabilitation therapy with an upper-limb exoskeleton; or automatic personalized robot-assisted rehabilitation therapy with an upper-limb exoskeleton [31] (see supplementary material for details on patients' division during the experimental training (available online at stacks.iop.org/JNE/17/045002/mmedia)).
During the robotic assessment sessions and the treatment sessions involving a robot, the subjects were assisted by an exoskeleton for the upper limb (Arm Light Exoskeleton Rehab Station, ALEx RS) developed by Wearable Robotics srl [32,33]. The task proposed while working with ALEx RS was a 3D point-to-point reaching task. Specifically, the participants were instructed to start from the center of the workspace, reach one of the eighteen outer targets, and then move back to the starting position at a comfortable speed. Visual feedback was displayed on a monitor placed in front of the subjects ( figure 1(b)). A yellow sphere corresponded to the position of the exoskeleton's end-effector, while a red sphere indicated the target to be reached. The outer targets were in total eighteen, equally distributed along a spherical workspace of 19 cm radius ( figure 1(c)). The selected radius of the sphere allows for a maximum exploration of the workspace, while maintaining the reaching movement executable for people of most body sizes. This design of the motor task allowed exploiting an extensive three-dimensional workspace. The spherical workspace was positioned so that its center was aligned with the acromion of the right arm (i.e. arm trained with the exoskeleton) mid-way between the center of the target panel and the subject's acromion. In order to preserve the depth perception, the dimensions of the target spheres were modified in accordance with their position in the 3D space. If a subject was unable to reach a target (i.e. the subject did not move for more than 3 s), ALEx RS activated its assistance mode to guide the subject towards the target according to a minimum jerk speed profile [33,34]. During the assessment sessions, the subjects were asked to reach all the eighteen targets as many times as possible within 30 min, while during the robotic training the subjects practiced only a subsection (8 targets) of the available workspace. During each training session, the choice of which target to reach among the available eighteen was done either by the physical therapist (standard robot-assisted rehabilitation therapy) or by an automatic-personalized algorithm. The algorithm was embedded in the control scheme of ALEx RS and selected the targets at the beginning of each session and changed them within the same session accordingly to the subject's performance (automatic personalized robot-assisted rehabilitation therapy) [31]. Patients that performed standard rehabilitation without ALEx RS were trained with a similar amount of upper limb movements.
In addition to the experimental training sessions described above, all the patients received habitual physical therapy at the stroke unit during the acute and sub-acute phase, while hospitalized. Specifically, all patients received two times 30 min of physical therapy per day on five days per week and five times 30 min of occupational therapy per week on an inpatient basis for 8 to 16 weeks. Following the end of the hospitalization, patients continued receiving outpatient treatments consisting of 1-4 h of physical and occupational therapy per week. Therapy was adapted by the therapists to the current capacities of each patient by choosing from a list of appropriate exercises comprising upper-extremity relaxation techniques, unilateral task-specific mobilizations, bilateral upper limb exercises with a wand, ball exercises, active ante/retropulsion exercises, active pronation/supination exercises and grasping exercises. All subjects were hospitalized until the assessment performed at the end of the experimental treatment (A3). At the time of A4 all subjects were no longer hospitalized and participated at the assessment as outpatients. Therefore, we can considerer the dosage of inpatient and outpatient therapies globally equivalent across subjects.
The healthy subjects went through a single robotic assessment session where they were asked to reach the eighteen targets 5 times.
During each robotic assessment session, both for stroke patients and healthy controls, muscular (EMG) and brain (EEG) activity were recorded. Resting state eye-closed EEG was recorded for 5 min before the beginning of each robotic assessment session. Kinematic data acquisition was synchronized with EMG signals acquisition by using trigger signals sent from ALEx RS at the following events: movement start, movement end, and at the occurrence of roboticassistance.
For the goal of this study we considered all 6 stroke subjects together independently from the therapy received during the experimental training period. Therefore, we are going to compare motor performance of two groups: healthy and stroke subjects.

Clinical evaluation
Along with the robotic assessments, and at the same time points (A1, A2, A3 and A4), the stroke subjects were evaluated with clinical tests by a therapist not directly involved in the study and blind to the group allocation. The sensorimotor status of the patient was evaluated using the upper limb section of the Fugl-Meyer Assessment (FMA) scale [35]; and the grip force measured using a Jamar dynamometer [36].

Data analysis
For the data analysis of all subjects, stroke and healthy, we considered only the six targets (figure 1(c), blue spheres) that all the stroke patients were able to actively reach, which were, in addition, representative of the 3 main movement directions: up and down (target 1 and target 5), right and left (target 3 and target 7), near to and far from the body (target 10 and target 13). We only considered the movements going from the center of the workspace to the outer targets. The selected set of movements differed in terms of elbow and shoulder joints combination in order to fully capture subjects' impairment. In particular, target 13 elicited elbow extension, which is generally more difficult than elbow flexion (i.e. target 10). Target 7 required extension of both elbow and flexion of the shoulder as opposite to target 3, which involved a coordination of elbow extension and shoulder flexion. Finally, movements toward target 5, the one placed at the bottom of the spherical workspace, were facilitated by gravity. All reaching movements performed by the subjects were requiring active gravity compensation of the subject's arm, while the exoskeleton was compensating only for the weight of its components. In all further analyses, we did not consider each selected target separately but we pulled all movements together.
For the healthy subjects the total number of movements was 30 (i.e. 5 repetitions of 6 targets), while for the stroke subjects it depended on the level of residual mobility at each assessment. The session lasted 30 min, during which the stroke subjects could perform a minimum number of one repetition of each target and a maximum of 5 repetition of 6 targets.

Kinematic analysis and measures
Kinematic parameters were computed from the handle (i.e. exoskeleton's end-effector) positions recorded by ALEx RS during each reaching movement. The start and end of the reaching movement were defined as the time points when the speed profile of the EE of the exoskeleton respectively exceeded or dropped below 2% of the local maximum value [33]. We adopted the following parameters [31,37]: MV, the mean tangential velocity of the handle; nMD, the mean absolute value of the distance between the actual trajectory and the straight line connecting the starting position with the target (theoretical path) normalized by the length of the theoretical path, which is a measure of movement accuracy; nPK, the number of peaks in the speed profile, a well-known parameter quantifying movement smoothness; the spectral arc-length metric (SAL), (expressed as a negative value), that uses a movement speed profile's Fourier magnitude spectrum to quantify movement smoothness [38]; the distance from target when the exoskeleton started to assist the movement normalized by the distance between the target and the center of the workspace (Dtrgt); the robot assistance frequency, i.e. the number of assisted movements (RAF, Figure 1. (a) experimental protocol: the protocol was organized in four assessments (A1 to A4) and four weeks of training. During the assessments, clinical scales were measured along with kinematic, muscular activity and brain activity before and while training with the exoskeleton. (b) experimental set up: the subject was sitting in front of a computer screen while wearing the exoskeleton and was given a visual feedback of the end effector of the robot in the shape of a yellow sphere. A red sphere indicated the position of a target that the subject had to reach for. (c) reaching task: the reaching task consisted in center-out reaching movements towards 18 targets equally distributed along a sphere. The blue targets were the ones considered in the multimodal performance analysis, while the set of blue and red targets together constituted the full training target set. in percentage); the percentage of workspace explored without the help of the robotic-assistance normalized by the ideal workspace corresponding to the volume of a sphere of radius 19 cm (WS); and the time to complete the task (t task , i.e. reaching of the outer target). MV, nMD, nPK have been computed only over the path's trajectory in which the patient moved actively.
EMG data were preprocessed offline using MAT-LAB (MathWorks, Natick MA). The raw EMG signals were detrended, band-pass filtered 50-500 Hz (Butterworth filter, 7th order), rectified, low-pass filtered with a cut-off frequency of 10 Hz (Butterworth filter, 7th order) to obtain the envelopes [41,42]. To correct the EMG-amplitude differences due to electrode placement and to ensure that the extraction of the synergies would not be biased against the low-amplitude muscles, the envelope of each muscle signal was normalized by the median computed for each individual across each session. The normalization based on the median value instead of the maximum is more robust to outliers [41]. Then, for each subject and session, EMG data were epoched considering the six different directions and were concatenated for muscle synergies analysis.

Muscle synergies extraction and measures
For each subject, muscle synergies were extracted by using the non-negative matrix factorization algorithm (NNMF) [43]. The NNMF algorithm decomposes the EMG envelope in a defined number of positive components or muscle synergies. The organization of a synergy is determined by the contribution (weight coefficient) of each muscle, as specified by the weight matrix W. Its activation profile is defined by the activation coefficients, specified by the matrix H [43]. Since, the iterative algorithm can find a solution as a local and not global minimum, each extraction was repeated 50 times, and the repetition with the solution explaining the highest overall amount of EMG variance was selected. For each subject, to objectively determine the minimum number of muscle synergies required to reconstruct the data set, we used the variance accounted for by the synergies model and chose as threshold the number of synergies at which the VAF was higher than 95% [41]. The same number of muscle synergies was retained for all healthy subjects to allow an easy intra-task comparison (i.e. mean number of muscle synergies across healthy subjects). For the patients, instead, the individual number of muscle synergies was considered. Indeed, the number of retained synergies has been proposed as a biomarker of cortical damage and recovery [23].
Muscle synergies were matched among subjects and sessions according to their similarity (determined by using normalized scalar products) with a set of reference synergies. The set of reference synergies was obtained by grouping the muscle synergies of the healthy subjects with a hierarchical clustering procedure based on minimization of the Minkowski distance between weighting coefficient vectors [41].
Muscle synergies structure was compared within the healthy population and between the stroke patients and the healthy subjects by using the scalar products (dot HEALTHY and dot STROKE respectively). When a stroke subject presented a lower number of synergies, the value of the dot product between the healthy synergy and the correspondent missing synergy of the stroke was set to zero.

Spinal maps extraction and measures
The same preprocessed EMG data were resampled on the minimum number of time points (2 s) across directions, sessions and subjects and were used to estimate the spinal maps along the rostro-caudal direction related to C2-T1. The spinal maps describe the spatiotemporal organization of the peripheral EMG signals at the level of the spinal cord, by estimating the motoneuronal (MN) activity for each spinal segment as previously described in literature to investigate the muscle activity in the lower and upper extremities [44][45][46]. The weight coefficients approximating the rostro-caudal distribution of the MN pools innervating the upper limb muscles included in the study were located in the segments from C2 to T1, as reported by Kandel [47] (see table 1 supplementary material).
In order to assess the similarity between two different spinal maps, we used the 2D correlation coefficient between two maps [46] (R MAP,H and R MAP,S when comparing respectively the spinal maps within the healthy population and between the healthy population and the stroke patients) and the root mean square error (RMSE MAP,H and RMSE MAP,S ), calculated as where MAP 1 and MAP 2 are the two spinal maps under comparison. N (i.e. number of spinal segments) is equal 8 and M (i.e. time samples) is equal 3000 samples.
We first computed these measures comparing the maps within the healthy population. Then, we compared the maps of each stroke subject and for each session with those of the healthy subjects.

EEG analysis 2.7.1. Data acquisition and pre-processing
EEG data were continuously acquired at 500 Hz using an Active II EEG system (Biosemi, Amsterdam) with 64 pre-amplified (active) EEG channels with standard 10-20 configuration. EEG data were preprocessed offline using MATLAB (MathWorks, Natick MA) and EEGLAB toolbox [48]. The raw EEG data were filtered (1 Hz to 40 Hz, Butterworth zero-phase 8th order IIR filter [49][50][51]) and down-sampled to 128 Hz. EEG electrodes with prolonged prominent artifacts (assessed by visual inspection) were removed and interpolated using spherical interpolation [52,53]. The data were, then, re-referenced to a common average. Finally, the data were visually inspected to remove periods contaminated by artifacts (i.e. amplitude >80 µV) and the remaining data were concatenated.

Resting-state EEG measures
Resting-state EEG measures were obtained by singular value decomposition (SVD) of the EEG signals [54,55]. We followed the same approach of our previous work [53,55]. The SVD of a real matrix is a factorization of the form E = MSN ′ . In this case, E ∈ R C × R T was the matrix of the pre-processed EEG signals concatenated across sessions and participants with C equals to 64 EEG channels and T equal to the total time of recordings summed over sessions and participants. The left-singular vectors of M ∈ R C × R C are a set of orthonormal eigenvectors that in our case represented the group-level EEG-SVD topographical maps. The right-singular values vectors of N ∈ R T × R C represented the group-level temporal courses. They are ranked according to their non-zero singular values (i.e. diagonal values of S). In order to obtain the subject-specific temporal courses of each EEG-SVD topography for each participant and session, the EEG activity of each individual for each session was projected on the group-level EEG-SVD topographical maps corresponding to 75% of the explained variance [53,55], which was calculated from the singular values as: For instance, the time courses for the first {1, · · · , p} SVD components for one subject was obtained as E 1 = M T 1,...,p E 1 . We assessed reproducibility of the group-level topographies by a split-half reproducibility analysis. We randomly split the original set of subjects and sessions (i.e. the six healthy controls and the four assessments of the patients) into 2 groups, each with 14 sessions. We generated ten different random splits. Subsequently, for each split, EEG-SVD topographies were computed for each group concatenating the data of all subjects and sessions within the group. EEG-SVD topographies obtained from the two groups were then matched using Hungarian algorithm, and their similarity was assessed by Pearson correlation [56].
For each topographical map and each subject/session, a time-frequency representation of the corresponding subject-specific temporal course was calculated from 1 to 64 Hz by convolving the signals with a complex-valued Morlet wavelet with 3 cycles. Time-frequency power was calculated as the squared magnitude of the complex wavelet-transformed data. We then computed coefficient of variation (CVs-i.e. ratio between variance and mean spectral power over time) for four typical frequency bands (i.e. δ: 1-4 Hz; θ: 4-8 Hz; α: 8-12 Hz; β: 15-30 Hz).
In order to identify whether these resting-state EEG measures represented reliable biomarkers of motor recovery, we deployed a multivariate analysis of correlation, i.e. canonical correlation (CCA), between the FMA score and the CVs of each EEG-SVD topography and frequency band. If we consider the CVs of the EEG-SVD and the Fugl-Meyer scores as vectors of random variables (X = (x 1 , . . . , x n ) with n = 1 and Y = (y 1 , . . . , y m ) with m = 12-i.e. 4 frequency bands per 3 EEG-SVD components), and there are correlations among these variables, then the CCA would find linear combinations (i.e. canonical components) of the CVs that have maximum correlation with the Fugl-Meyer. Specifically, the CCA seeks vectors a ∈ R n and b ∈ R m such that the random variables a T X and b T Y maximize the correlation p = corr(a T X, b T Y). We used a permutation test to very the significance of the found canonical components. In details, the CVs of each EEG-SVD topography and frequency band were permuted over subjects and we then used CCA to find correlation between these permuted data and the Fugl-Meyer. We repeated this procedure 1000 times and considered the 99th percentile of the resultant distribution of correlation as a significance threshold. Once identified the significant canonical components, the brain canonical scores (V) were obtained by projections over the canonical correlation vectors (V = b T Y). We then evaluated the changes of the brain canonical scores over assessment sessions for each individual patient.

Multimodal analysis
In order to quantitatively assess the presence of a relationship between the various parameters extracted from the different sources (kinematic, muscle and brain activity) we ran a correlation analysis among all of them. In this way, links between specific pairs of parameters extracted not only from the same type of signal but also from different ones (i.e. kinematic and brain, or muscle and brain) can be highlighted. In particular, we extracted Pearson correlation coefficients from a matrix Q (28 × 18), where each column was a parameter and each row an observation. Data of stroke and healthy subjects were pulled together. We obtained a matrix of correlation coefficients resulting from pairwise comparisons between the columns of Q and we retained only the correlation coefficients that resulted significant (p-value < 0.05).
To identify the neuro-biomechanical parameters that were most important for post-stroke motor recovery, we applied a multistep statistical procedure based on principal component analysis (PCA) [57]. Specifically, PCA was applied to all the parameters extracted from the kinematic, EMG and EEG analysis for all healthy and stroke subjects together. To avoid introducing bias due to the different scales of the various parameters, we normalized measures before running PCA, so as to have zero mean and a standard deviation equal to 1 (z-score). We retained the first three principal components (PCs), which explained more than 50% of the total variance, and we projected the original dataset in the 3D space defined by the constructed PC1-3. For each distribution, the coordinates of the centroid were computed by averaging all the coordinates of the points included in that distribution. We then computed the distance between the centroid of the distribution of the healthy controls and the projection of the data of each stroke subject in the identified PCs space, for each session. The distance was computed in the 3PCs space (dist ALL ) but also along each individual PC (dist PC1 , dist PC2 , dist PC3 ). Finally, to determine the relation between the clinical outcome and the identified neuro-biomechanical measures, we correlated the PC scores with the clinical scores computing Pearson correlations coefficients.

Statistical procedures
In order to summarize the information, the results reported were averaged across subjects. All data are reported as mean values +/−standard error

Clinical outcome improvement
Patients showed a reduction of impairment along with an increase of grip forces during the evaluation time. Significant improvements in the FMA (χ 2 = 16.68, p < 0.001) and in the grip strength (χ 2 = 9, p = 0.020) occurred throughout the four assessments for all six patients. Changes in the clinical scores ensued between the first (A1) and the second assessment (A2) (on average 6 ± 2.5 points, p = 0.031, in the FMA and 1.12 ± 0.64 Kg, p = 0.250, in the grip strength) and likely reflect the shortterm changes typical of the subacute phase (see table  2). Yet, stronger improvements occurred between A2 and the assessment after the training (A3). On average all subjects increased of 7.8 points (±2.1 SE, p = 0.031) in the FMA and 1.08 Kg (±0.49 SE, p = 0.250) in the grip strength. Finally, these measures did not change significantly at the assessment performed one month following the end of the training (A4; on average 2.4 ± 1.2 points, p = 0.250 in the FMA and 1.14 ± 0.49 Kg, p = 0.250, in the grip strength).

Motor performances improvement is captured by the kinematic measures
Motor performances significantly improved in all patients, with the exception of the dimension of the workspace ( figure 2, table 3). Between A1 and A2 patients' movement increased in speed (on average from 0.074 ± 0.012 m s −1 in A1 to 0.094 ± 0.012 m s −1 in A2) with a concurrent reduction, yet not significant, of the jerkiness and of the robotic-assistance (respectively from −5.44 ± 0.60 to −5.11 ± 0.36 and from 57.95 ± 18.84% to 47.15 ± 20.40%). After the month of training these improvements were stronger and statistically significant. Mean velocity and SAL reached level comparable with those of the controls (respectively 0.151 ± 0.013 m s −1 and −3.67 ± 0.21 for stroke and 0.153 ± 0.011 m s −1 and −4.18 ± 0.37 for healthy subjects, differences between the two groups were not significant: for the velocity p = 0.589 and SAL p = 0.240) and subjects S3, S5, and S6 did not require anymore the robotic assistance. Noticeably, subjects S2 and S4, even if they still needed robotic-assistance in all the assessments, they were able to complete a larger portion of the workspace without the robot, as highlighted by a decrease in D trgt , (overall along the four assessment from 1.39 ± 0.15 to 0.42 ± 0.07 for S2 and from 1.61 ± 0.75 to 0.25 ± 0.11 for S4). Finally, along the four assessments, patients performed straighter movements (i.e. reduction of nMD: from 0.191 ± 0.029 at A1 to 0.094 ± 0.011 at A4, χ 2 = 9.24 p = 0.026) in a shorter amount of time (from 0.191 ± 0.029 at A1 to 0.094 ± 0.011 at A4, χ 2 = 13.77 p = 0.003).

Muscle coordination improves with training as highlighted by muscle synergies
Six (5.8 ± 0.5) muscle synergies were found for each healthy participant, when considering the movements from the center of the workspace to the six main outer targets ( figure 3(a)). The muscle synergies structure was similar to that already reported in literature for analogous tasks [33].
Specifically, synergy 1 involved mainly the pectoralis, which was responsible for the arm flexion. Synergies 2 and 4 were dedicated to the extension of the arm and included respectively the anterior and medial part of the deltoid (Syn 2), and the two triceps and the posterior part of the deltoid (Syn 4). Synergy 5 was mostly composed of the activation of the forearm muscles (i.e. BRAD and PRO). Finally, synergies 3 and 6 accounted for the activity of the postural muscles (i.e. both trapezii, INFRA and RHO for Syn 3, and LAT, INFRA and RHO for Syn 6) involved in the elevation of the arm against gravity. The stroke subjects presented a significant increase in the number of synergies from A1 to A4, χ 2 = 8.54 and p = 0.036. Specifically, the number of synergies augmented between A2 (4.1 ± 0.4) and A3 (5.0 ± 0.3) and it was maintained at A4 (5.4 ± 0.4 synergies modules) ( figure 3(b)). This change highlights an increased level of complexity of muscle activity as expected during motor recovery [24]. Additionally, through time and training the structure of the synergies of the stroke subjects became more similar to the  structure of the controls. Indeed, the dot STROKE augmented significantly for each stroke subject, χ 2 = 9.49 and p = 0.023. Specifically, it was stable during the first two assessments (dot STROKE = 0.48 ± 0.07 and dot STROKE = 0.47 ± 0.04 at A1 and A2, respectively). At A3, instead, there was a significant improvement of this indicator with respect to A2 (dot STROKE = 0.58 ± 0.03, p = 0.031), and it continued to significantly evolve also at follow-up (dot STROKE = 0.67 ± 0.03, p = 0.031).

Motoneuronal activity in the spinal circuits increases with training
We further investigated the global muscular strategies adopted by the participants by computing spinal maps, which estimate spatiotemporal motoneuronal activation. For the controls, spinal maps were characterized by a main period of activation towards the end of the reaching task, localized around C7-T1 ( figure  4(a)). For the stroke subjects, instead, at A1 and A2 the motoneuronal activity was either almost absent because they were not able to accomplish the task without the help of the exoskeleton or it was diffused across all the spinal segments and not localized, in contrast to the controls ( figure 4(b)). Along the four assessments, the spinal maps of the stroke subjects became significantly more similar to those of the controls, both in term of correlation R MAP,S , χ 2 = 9.56 and p = 0.022, and of root mean square error RMSE MAP,S , χ 2 = 10.62 and p = 0.014. Yet, the greatest and the only significant changes for both metrics occurred between A2 and A3 (i.e. with training), where R MAP,S augmented from 0.18 ± 0.06 to 0.42 ± 0.08 (p = 0.031) and RMSE MAP,S diminished from 60.50 ± 5.17 to 46.33 ± 3.23 (p = 0.031). Importantly, the improvement remained stable at follow-up as indicated by a non-significant difference of the values of R MAP,S and RMSE MAP,S between A3 and A4.

Cortical activity gets closer to healthy level with training as highlighted by EEG topographies
We utilized SVD decomposition of the EEG signals concatenated in time across participants and sessions to identify reliable and reproducible topographical maps [53,55]. We selected the first three EEG-SVD components, which accounted for 75% of the variance ( figure 5(a)). We tested the consistency of the maps by split-half reproducibility analysis. All the three components were highly reproducible (average correlation ± STD across splits and components: 0.98 ± 0.01).
We computed canonical correlation between the FMA scores and the coefficients of variation of each EEG-SVD topography and frequency band. We found one significant canonical correlation component (p < 0.01 permutation test; correlation: 0.96, figure 5(b)) with highest coefficients for CV of delta and alpha bands for the 1st EEG-SVD component, theta and alpha bands for the 2nd component, and delta and beta bands for the 3rd component ( figure  5(a)). Through time and training the CVs got closer to those of the healthy controls, in particular for S3, S5, and S6. Yet, no significant difference was found between the three assessments when considering all subjects together.

Multimodal analysis reveals different temporal dynamics of recovery
Before running the multimodal analysis, we looked for correlations between specific pairs of parameters obtained both from the same domain and from different domains. Results are summarized in figure  6(a). We found that there were some significant correlations between couples of parameters extracted from the same domain, like between the nPK and time to reach the target or the linearity index (nMD) and the robotic assistance (RAF) obtained from the kinematic domain. Similarly, we found significant correlations between dot syn1 and dot syn3 or between the number of synergies and dot syn6 . Interestingly, we noticed significant correlations also between parameters of different domains, for example the percentage of workspace explored without robotic assistance and the number of muscle synergies or the time to reach the target and the RMSE MAP extracted from the spinal maps. Nevertheless, despite these significant correlations, the value of the correlation coefficient was, on average, relatively modest: r = mean 0.53 ± STD 0.12.
We then combined together all the measures extracted from the robotic assessment (kinematic, EMG and EEG metrics) in a multimodal analysis using PCA. The explained variance for the first three PCs was 59.75%. The metrics extracted from the kinematic and spinal maps were those mostly contributing to PC1. The second PC, instead, was mainly associated with the workspace coverage, the number of synergies, and the similarity of synergy 2 and 6. Finally, the factor loadings associated with  (a) Canonical correlation coefficients for the coefficients of variation of the first three EEG-SVD topographies and the four frequency bands (delta, theta, alpha, and beta). (b) 2D correlation plot of canonical correlation scores for FMA (y-axis) and resting-state EEG measures (y-axis). Black dots represent healthy subjects. Blue, cyano, and light and dark green dots represent stroke patients at A1, A2, A3, and A4, respectively. (c) Brain canonical scores for each subject for A1 (blue bars), A2 (cyano bars), A3 (light green bars), and A4 (dark green bars) and for controls (grey bars). the third PC were the number of peaks in the speed profile, movement duration, the workspace coverage, the structure of synergy 1, and the EEG scores ( figure 6(c)). When projecting the individual data points into the 3D space defined by the newly constructed variables PC1-3, clear differences emerged between the healthy subjects and the stroke patients ( figure 6(b)). Importantly, these differences were, also, affected by time and training (figures 6(b) and (d)). Indeed, the distance between the patients' and controls' data points in the 3D space significantly decreased along the four assessments, χ2 = 10.22 and p = 0.016. In particular, the distance decreased significantly (p = 0.031) from A2 (dist ALL = 4.07 ± 0.32) to A3 (dist ALL = 2.52 ± 0.37), and it remained stable at A4 (dist ALL = 2.17 ± 0.36). Interestingly, we found different temporal dynamics of recovery for the three principal components, summarized in figure 6(e). Indeed, the distance along PC1 changed along the four assessments, χ 2 = 9.61 and p = 0.022, with a significant decrease in particular between A1 and A2 (dist PC1 = 4.73 ± 0.98 and dist PC1 = 2.43 ± 0.56 at A1 and A2, respectively −p = 0.031) and between A2 and A3 (dist PC1 = 0.80 ± 0.17 −p = 0.031). PC2 significantly decrease between A2 and A3 (dist PC2 = 2.76 ± 0.29 and dist PC2 = 1.51 ± 0.23 at A2 and A3, respectively-p = 0.031). Instead, the distance along PC3 decreased throughout the four assessments, χ 2 = 8.39 and p = 0.03, yet mainly between A1 and A2 (dist PC3 = 1.91 ± 0.53 Projection of the dataset in the 3D space identified by the first three principal components (PCs). The data relative to the healthy subjects are in black, the data of the stroke subjects at A1 in blue, A2 light blue, A3 light green and A4 dark green. The centroid of each group of data is highlighted by a marker of bigger size and higher intensity with respect of the majority. (c) Coefficients of the three PCs. In the shades of red data relative to the kinematic, in the shades of brown data relative to muscle synergies, shades of purple data relative to the spinal maps and in orange data relative to EEG. (d) Top: Distance between the healthy and the stroke group in the space of PCs1-2-3 (distALL) and along each PC separately (distPC1, distPC2, distPC3). Bottom: Distance between the healthy group and each stroke subject in the space of PCs1-2-3 (distALL) and along each PC separately (distPC1, distPC2, distPC3). Highlighted in grey the regions where there was a significant change. (e) Summary image illustrating the different temporal dynamics of recovery identified by the proposed methodology. dist PC3 = 0.95 ± 0.14 at A1 and A2, respectivelyp = 0.031). When correlating the 3 PCs with the clinical scores we found a significant correlation between PC1 and the FMA, r = 0.74 p < 0.001, and between PC1 and grip force, r = 0.54 p = 0.003.

Discussion
The current study aimed to introduce a comprehensive set of movements, muscle, and brain measures as well as to present a multivariate methodological approach combining them in a unique framework. This approach can be applied to any dataset containing multidomain measures. Here we applied this method to a multidomain dataset, including a small group of subacute stroke subjects that went through a month of intense motor rehabilitation.

Multimodal assessment
Motor rehabilitation aims at maximizing the recovery and independence in daily living by discouraging dysfunctional compensatory behaviors and promoting the re-learning of appropriate motor control strategies [58]. For this, rehabilitation protocols have to adapt to the individual dynamics of recovery [59]. However, because of the multifaceted nature of stroke, recovery assessments should explore different aspects of the nervous system. For this reason, in this study, we simultaneously recorded kinematic, EMG, and EEG signals to gather a general understanding of the neuro-biomechanical state of the subjects and of its evolution over time after a rehabilitative treatment, independently from the nature of the latter (e.g. robotic or conventional therapy). Importantly, all domains of the neuromuscular system-kinematic, muscle, and neural activity-were affected by the motor training. Besides, significant correlations could also be observed across them, although they could not fully capture the complex variations pertaining to the different domains. As such, these observations support the feasibility and the need for a multimodal approach that could efficiently summarize these factors in a single metric. For this purpose, we deployed principal component analysis. We projected the individual data points into the 3D space defined by the first three principal components explaining a high percentage of variance, and we computed the distance between the healthy controls and the stroke subjects at each assessment. During the training, the distance between the patients' and controls' data points in the 3D space significantly decreased, yet with different dynamics. Previous studies already reported different evolutions of the recovery factors for movement speed, efficiency, and accuracy [60][61][62]. We here extended these findings to new physiological measures, hence providing a more complete characterization of the neuro-biomechanical status of the patients and of its evolution. We found that patients' kinematic and spinal maps parameters were mainly clustered over the first principal component (PC1) and improved continuously during the four assessments. These quantitative measures were also those exhibiting a higher correlation with the post-stroke impairments captured by the FMA assessment and by the grip force. A strong correlation between kinematics and the patient status as described by clinical scales is in general well accepted and expected since they measure similar parameters [63,64]. A more striking observation is, instead, the strong correlation with the spinal maps, as they have been so far rarely used to describe upper limb movements [33,65] and levels of impairment after stroke [46].
Together with spinal maps, we also estimated muscle synergies as metrics of motor coordination extracted from electromyographic signals. Interestingly, muscle synergies features clustered in a separated principal component (PC2), as compared to spinal maps, and they normalized already in the shorter period of the training, highlighting that motoneuronal activity and muscle coordination evolve differently over time.
Finally, the third component mainly included metrics related to some kinematics measures (speed profile, movement duration, and workspace coverage), the structure of synergy 1, and the brain features that normalized in the shorter period post-lesion (i.e. between A1 and A2).
Overall these results highlighted that post-stroke recovery develops at different stages for different aspects of the nervous systems. A short-term recovery (between A1-A2) likely due to a combination of spontaneous recovery and inpatient therapy; a medium-term recovery (between A2-A3) probably as result of our intervention and the inpatient rehabilitative therapy and finally, a long-term recovery (between A3-A4) that showed a maintenance of the functional improvements achieved at A3 probably supported as well by the outpatient therapy. While further studies should confirm these preliminary results in larger datasets, these different dynamics were previously not captured when considering standard clinical approaches or mono-parameter analysis and should be taken in consideration when designing new rehabilitation protocols.

Specific unimodal analysis
Stroke subjects performed the 3D point-to-point reaching task worst as compared to healthy controls. In accordance with previous studies [4,5,28,66] patients' movements were less smooth and slower, and often required assistance from the exoskeleton, especially during the A1 assessment, a few weeks after the injury. Already at the second assessment, which was performed before the extra dose of rehabilitative intervention (either robotics or standard therapy), performances improved. This progress immediately after the injury is typical of the recovery that usually arises in the acute phase [67]. Moreover, in this period patients were already receiving inpatient physical and occupational therapy. Yet, the strongest improvement occurred between A2 and A3 (i.e. after the treatment). At A3, kinematic performances became comparable to those of the controls and maintained an analogous pattern at follow-up (i.e. assessment A4).
The multimodal set-up deployed in this study linked this trend of kinematic improvements to the motoneuronal changes estimated using spinal maps. The use of spinal maps to characterize upper limb movements had so far been limited [33,65,68] and they were never exploited in stroke patients to explore upper limb impairments. Moreover, we recently demonstrated a strong correlation between spinal maps activation and blood-oxygen-level-dependent signal of the spinal cord captured by functional magnetic resonance imaging, which further supports the use of this analysis technique to infer motoneuronal activity [69]. In our stroke population, the motoneuronal activity described by the spinal maps was reduced after the stroke lesion with respect to the one of healthy participants, both in terms of temporal activation (R STROKE ) and magnitude of activation (RMSE STROKE ), but it increased during the training, getting closer to the one of the healthy participants at A3 and A4.
The number and structure of the patients' muscle synergies also became more similar to those of the healthy controls, in particular after the month of intense training. Similarly to observations reported in previous studies, the number of synergies for the stroke subjects increased through time with the recovery of motor functions [23,24,41,46] and this increase was also accompanied by changes in their structure. The synergies structure improved also for the subjects that did not recover all the six synergies characteristic of the healthy population (see also figure 1 in supplementary material for more details). These results are in accordance with studies showing that the conjunction between the spontaneous recovery and the intensive treatment improved motor performance, both in terms of kinematics and muscle activation patterns [24,70,71]. Interestingly, muscle synergies related to the control of the shoulder (i.e. Syn 2 and 4) were those with the strongest changes particularly after training with the exoskeleton (see figure 1 supplementary material) paralleling our previous findings [24].
The kinematic and muscle factors explored in this study provide a detailed description of the biomechanical status of the patients. Yet, they do not tap into the complex neural reorganization processes that occur after brain insult. To capture these changes, we deployed topography-based analysis of the EEG signals acquired at rest. We opted for a topographybased approach, in contrast to traditional EEG waveform analysis, as we believe this approach can better capture post-stroke large-scale neural processes without any a priori hypothesis on the spatial location of abnormal brain activity [52,72]. Indeed, stroke has been nowadays reconsidered as a distributed network disease with structural and functional changes occurring between brain areas distant to the lesion [29,73,74]. We deployed singular-value decomposition to extract EEG topographies that resemble previously identified microstates, whose preserved occurrence and duration has been shown to correlate with a better effective motor recovery [72]. Specifically, in our approach we concatenated in time the data recorded for healthy subjects and patients and subsequently derived group-level spatial maps that come with subject-specific time courses. Indeed, we were interested on preserving the subject-specificity only in the time courses, while considering the same spatial subspace across subjects. We further supported this assumption by performing SVD decomposition of the EEG signal of each participant and every session (see supplementary figure 2), which showed a high correlation across all subjects for all top-three components. We recently deployed similar analysis to discriminate patients with spatial neglect of different severity levels [53]. Here, we paralleled these previous results showing that spectral power of the SVD topographies, which has a typical 1/f fall-off typical of the EEG spectral power for all components, and particularly delta and alpha rhythms, correlated with level of motor impairments. Interestingly, through time and training the aberrant brain oscillation patterns were restored, in particular for the subjects that showed a stronger motor improvement. Indeed, the three patients (S1, S2, and S4) for which brain rhythms remained abnormal not only had lower FMA scores, but also aberrant motoneuronal activation. Furthermore, they still needed robotic assistance after the intensive training and at follow-up.
Overall kinematic, muscle, and cortical activity showed an improvement particularly evident during the training and correlated with the clinical status of the patients. Yet, when summarizing these affinities across the different physiological measures using a multivariate approach, the latter highlighted that post-stroke recovery develops at different stages for different aspects of the nervous systems, suggesting that kinematic, muscle, and brain features should be taken into consideration when designing new rehabilitation protocols. For this, our approach could be an effective method to summarize the different psychological measures in few metrics that could be monitored over time to personalize the treatment [31,37].

Limitations of the study
Although the findings presented in this study suggest that the proposed multimodal approach has the capacity to provide additional information regarding the evolution of the recovery after stroke, the limited sample size limits the generalizability of the results. Nevertheless, the main goal of this study was to propose a potential comprehensive method that included various types of information rather than probing its clinical efficacy. Moreover, the heterogeneity of the stroke subjects included in the study does not allow for an ideal homogeneous stratification of the population. However, this could not be precisely controlled, as we were working with subacute stroke and the recruitment was done on a continuous basis one patient at the time. Therefore, further studies including larger cohorts of participants would be necessary to draw meaningful conclusions about the clinical efficacy of the presented approach. One further limitation regards the number of movements performed by the patients recruited in the three different groups. To ensure comparability among the groups, we proposed the same duration of the session and we asked the therapist to provide a training with a similar intensity as the one provided by the robot. Yet, the number of moments performed highly depended on the ability level of each stroke subject.

Conclusion
The use of an instrumental assessment followed by a multimodal approach identified quantitative neurophysiological metrics correlating with clinical measures such as the FMA and grip force, as well as clustered metrics that evolve distinctly during recovery, underlining their functional and clinical relevance. A combined analysis of kinematic, muscular, and brain activity seems to be able to provide a good and accurate patient characterization in line with the outcome of the clinical scales. In the future, similar methods should be implemented in order to track the evolution of the neuro-biomechanical state of the patients after brain damage, to define suitable personalized rehabilitative intervention strategies and to provide a deeper insight into the recovery process after stroke.