We recruited 13 right-handed volunteer participants (Table 1); participants did not present any pathology that affected the upper limb, spine or posture. Volunteers were excluded if they reported neck, shoulder, or arm pain (> 2 on a 1–10 verbal scale) within the last 3 months. The University of Auckland Human Participants Ethics Committee approved the research protocol and methods of the study (reference number 022246), and informed consent was gained before participation in any procedure.

The participants attended a single session where they performed three upper limb isometric contractions tasks: maximal voluntary force (MVF), multidirectional trials exploring a high spatial volume, and synergy-tuned trials aimed at individual synergies’ preferred spatial directions. A similar protocol has been described previously (Ortega-Auriol et al. 2018).

The MVF tasks consisted of a maximum average force from three trials of external shoulder rotation. Shoulder’s external rotation is the weakest degree of freedom (DoF) for force development of the shoulder; consequently, forces in other DoFs will require less force relative to the MVF. Tasks involving low-to-moderate forces do not modulate coherence (Poston et al. 2010; Mima et al. 1999), but high force levels can shift the observed CMC from β to the lower γ band (Brown et al. 1998; Roh et al. 2012). A standardisation by shoulder external rotation allowed for low-to-moderate force levels across all directions.

During the multidirectional and synergy-tuned tasks, the participants were seated and exerted an isometric force with their dominant arm. The task trials were directed in specific spatial directions at 40% of the MVF using a handle instrumented with a force transducer (Fig. 1A). This force level was selected to ensure that participants could sustain the required isometric contraction in all directions without substantial effort, minimising the presence of fatigue and postural compensations. Trials were accepted if a 3D spatial target was achieved during four consecutive seconds within a resultant force range of ± 5 N. The handle was located at a position calculated as 40% of the arm length (acromion—3rd metacarpal head) in front of each participant’s shoulder, providing a comfortable position for participants to accomplish the desired protocol. The tasks consisted of matching a movable sphere to a target in a 3D visual representation, facilitated by real-time virtual reality feedback on a screen in front of them. This virtual reality feedback was complemented by shadows underneath the spheres, enhancing depth perception and guiding participants in reaching targets located in front or behind the initial sphere. Participants were given a training period prior to the tasks, allowing them to familiarise themselves with the task requirements and environment, ensuring a high level of performance and comfort during the actual trials. The distance vector between the spheres was 40% of the MVF, creating a clear and consistent objective for the participants during the task trials.

Experimental setup. A EMG sensor placement (black dots, grey dots are located ventrally), virtual reality feedback (VRF), and instrumented handle. B 10–20 EEG setup schematic. C Representation of target directions of the multidirectional task. D Screenshot of the VR feedback displayed on the screen. Each VRF wall is located 100 N away from the centre

The multidirectional task consisted of isometric contractions to match targets in 26 different directions evenly distributed in a sphere (Fig. 1C). A significant number of muscle synergies were extracted from the multidirectional trials’ concatenated electromyography (EMG). From the activation coefficients of the extracted synergies, we determined the spatial tuning, known as the preferred direction (PD), of each of the extracted synergies. The synergy’s PD was determined as the direction in which the magnitude of activation coefficient of the synergy was maximal.

The synergy-tuned task consisted of several target matching trials towards the specific PD of a muscle synergy. The participants performed 50 trials in each synergy PD. The order of the trials was randomised and self-paced. The participants could rest between trials to avoid fatigue effects. A new synergy extraction was made from concatenated trials of each synergy’s PD. To corroborate the equivalency and conservations of synergies from the synergy-tuned trials, synergies were extracted again from the synergy-tuned task’s concatenated trials. Finally, CMC, CSC and IMC were calculated from the synergy-tuned trials.

Force was recorded at the handle with a 6-axis force transducer (Omega160, ATI Industrial Automation, USA) at 120 Hz. A Python-based custom software interface recorded signals. EEG signals were recorded from a 32 Ag/AgCl electrode EEG system (EasyCap; Brain Products GmbH, Germany). The electrodes were positioned according to the 10–20 system, referenced to the FCz channel, and offline to a common reference. Signals were acquired with BrainVision Recorder software (Brain Products GmbH, Germany) at 5 kHz.

Surface EMG signals were recorded from 16 single differential channels and sampled at 2 kHz using a Trigno device (Delsys Inc., USA). EMG activity was recorded from muscles of the participant’s dominant upper limb: superior (ST) and middle trapezius (MT), infraspinatus (Inf), teres minor (TM), serratus anterior (SA), anterior (AD), middle (MD), and posterior deltoid (PDel), pectoralis major (PM, clavicular fibres), short (BS) and long (BL) heads of biceps brachii, long (TL) and lateral (Tlat) heads of triceps brachii, brachioradialis (Braq), extensor carpi radialis (ECR), and flexor carpi radialis (FCR). These muscles were chosen based on their force capability and likely contribution to the required task, as essential considerations for accurate reconstruction of synergies (Steele et al. 2013). Electrodes were positioned according to SENIAM and Cram’s guidelines (Hermens et al. 1999; Criswell 2010). The participant’s skin was prepared with a gentle abrasive gel to clean and improve transmission before placing the electrodes.

Synchronisation across devices, EEG, EMG, and force acquisition were performed with Python-based custom software (https://dragonflymessaging.org/applications.html, U. of Pittsburgh). Data analysis was performed in MATLAB 9.3 (MathWorks, USA) using custom-made scripts and the FieldTrip toolbox (Oostenveld et al. 2011). A schematic representation of the workflow to process the data is shown in Fig. 2.

Data processing pipeline from MVF to statistical analysis

To process the EMG and EEG signals, all trials were split into a ramp and a hold phase based on the force profile. The ramp phase of a trial was defined as the time window from the initial movement of the feedback sphere until the force trace’s inflexion point (knee). A trial’s hold phase was defined as the intermediate 2 s of the required 4 s of target matching of a multidirectional or synergy-tuned trial. To trim the trials, force data (Fig. 3) were low-pass filtered (Butterworth, second-order, 5 Hz), and the inflexion point of the force traces was calculated by a custom algorithm and corrected through visual inspection if necessary.

Multidirectional task results. A Frequency of occurrence of extracted synergies per participant. Since five was the highest number of synergies across participants, the same number of clusters were extracted for subsequent analyses. B Mean VAF across participants for the extraction of N synergies, the synergies’ mode is noted with a red circle. C Average force traces across all participants and trials

The pre-processing of EMG signals for synergy extraction is described in detail in a previous report (Ortega-Auriol et al. 2018). Briefly, EMG signals were: band-pass filtered (bidirectional Butterworth, 2nd order, 5–400 Hz), demeaned, full-wave rectified, normalised to maximum activation across trials for each muscle to preserve relative contribution, converted to unit variance, low-pass filtered again to obtain an envelope (Butterworth, 2nd order, 5 Hz), and, only for synergy extraction, rebinned into 100 data points.

Non-negative matrix factorisation (Lee and Seung 1999) was applied individually to the processed concatenated EMG signals from the multidirectional and synergy-tuned tasks. NMF can be modelled as D = W × C + ϵ, where D is the original data set, W the synergy structure and C the activation coefficients, and ϵ is the variance not explained by the synergies. NMF was implemented using the multiplicative rule (Berry et al. 2007). The final solution was the resultant of 20 consecutive iterations with a difference of EMG reconstruction error smaller than 0.01% among them. First, to determine a significant number of synergies, the algorithm iterated from one until the number of muscles minus one. Second, we used the VAF metric (Cheung et al. 2005) to determine the number of synergies that achieved the best reconstruction of the original data. VAF was applied as a global (quality of original dataset reconstruction) and local criterion (quality of individual muscles’ signal reconstruction). A significant number of synergies was determined when global VAF ≥ 90% and local VAF ≥ 80% (Cheung et al. 2005; Chvatal and Ting 2012; Kim et al. 2016; Ortega-Auriol et al. 2018).

Once a significant number of synergies were determined from the multidirectional task, the PD of each synergy was computed as the average of each trial’s target direction scaled by the activation coefficient of that synergy during that trial (Eq. 1):

$$\overset$}}}_ }} = \frac \times C_ )}},$$

(1)

where Qi is the direction unit vector of the ith trial, Cri is the activation coefficient of the rth synergy of the ith trial, and T is the total number of trials. In a previous publication (Ortega-Auriol et al. 2018), our results showed that synergies’ activation coefficients are directionally scaled in the explored spatial volume, allowing the identification of synergy PDs.

EEG and EMG data from the synergy-tuned task were band-pass filtered (bidirectional Butterworth, second-order, 5–400 Hz), demeaned, EMG data were rectified via Hilbert Transform, and each trial was split into a ramp and a hold phase. Signals were rectified to emphasise the grouping and timing of action potentials within the EMG signals, as recommended for correlation and coherence analysis (Boonstra and Breakspear 2011; Farina et al. 2013; Ward et al. 2013). EEG signal impedance was checked at two instances during preparation: before starting data collection and in the time window between the multidirectional and synergy-tuned tasks. Electrodes with impedance levels over 15 Ohms were adjusted. EEG was downsampled to 2 kHz to match the EMG sampling frequency and optimise processing times. Independent Component Analysis was applied to EEG to remove electro-oculographic artefacts. To remove electro-oculographic artefacts while conserving the integrity of the data, a single component which visually presented artefacts and was located in the frontal region, was subtracted from the EEG data reconstruction.

To group similar synergies across participants, a cluster analysis was applied (García-Cossio et al. 2014; Roh et al. 2015) to the pooled synergy structures of all participants from the multidirectional and synergy-tuned tasks. Cluster analysis was applied using a k-medoids algorithm (Park and Jun 2009) with a cosine function as the cluster distance metric between members and the cluster’s centroids. The number of clusters was fixed to the maximum number of extracted synergies across participants. Membership assignment within a cluster was constrained to prevent the inclusion of two or more synergies from a single participant within a cluster. In these cases, the closest synergy to the centroid of the next nearest available cluster was reassigned. The reassignment process was iterated until no further repetitions were found. The synergy grouping results across participants is displayed as a mean synergy for each cluster.

IMC, CMC and CSC were calculated from the concatenated trials from each set of synergy-tuned trials in a single direction. To confirm consistency between the extracted synergies from the synergy-tuned trials and those from the multidirectional task, we compared them using cross-correlation and dot product. The synergies’ structures from the synergy-tuned task were not different from those extracted from the multidirectional task. IMC was calculated between three different muscle groups: (A) all–all, representing the average IMC of all muscles pairs within a single synergy, (B) high–high, being the highest IMC across the three muscles with the highest weights within a synergy, and (C) high–low displaying the highest IMC across three highest and the lowest weight muscles within each synergy. Muscle selection for group comparison was based on the muscles’ weights within a synergy. Muscle weights were determined in two different ways, first based on the mean synergies across participants and secondly on the individual synergy structure of each participant. CMC was calculated between each EMG channel and EEG channel in the motor area. CSC was calculated between the activation coefficients of the single synergy extracted per set of synergy-tuned trials and the respective EEG data.

After pre-processing, EMG, synergies’ activation coefficients, and EEG signals were transformed into the frequency domain to calculate the coherence measures. A fast Fourier Transformation (FFT) was applied to the bandwidth between 3 and 50 Hz. FFT results consisted of 24 frequency bins between 3 and 50 Hz with steps of 2 Hz. To narrow the scope of CMC and CSC calculations, we only considered channels on and near the motor cortex area: FC5, FC1, Fz, Cz, C3, T7, CP5, and CP1. CMC and IMC were calculated for the available combinations between EEG–EMG and EMG–EMG channels. CSC was calculated between a single activation coefficient and a set of synergy-tuned trials of EEG data. From the subset of analysed EEG channels, the one with the highest average CMC or CSC value was used for further analysis. Raw coherence (Rosenberg et al. 1989) was calculated using the FieldTrip toolbox (Oostenveld et al. 2011) applying Eq. (2):

$$C_ (f) = \frac (f)|^ }} (f)P_ (f)||}},$$

(2)

and values were normalised by applying a z-transformation (Baker et al. 2003; Reyes et al. 2017) using Eq. (3):

$$Z = \frac \right)}} }},$$

(3)

where c is the raw coherence value, and N is the number of tapers, which was four, used in the coherence calculation. From the individual estimates of coherence, pooled CMC, CSC and IMC were calculated to produce a single global estimate of correlated CMC or IMC (Amjad et al. 1997).

A significance threshold based on a surrogate coherence analysis derived from the original EMG and EEG data was calculated to define a significance level for the coherence calculations. Once the original data were transformed into the frequency domain, the argument of the complex quantity (angle of the polar form) was independently shuffled. The shuffling was iterated 50 times across all trials, channels, and participants. Then, coherence was calculated as described previously. This procedure allows the conservation of the power spectrum original amplitude structure of the signal while only shifting the signal phase, uncorrelating the signals in the time and frequency domain (Faes et al. 2004; Marchis et al. 2015). Coherence threshold significance was established as above the 90th percentile of the resultant by-chance coherence distribution.

Data distributions were checked for normality using a Kolmogorov–Smirnov Test (k-test), assessing for skewness and visually inspecting a normality plot. Then we compared the average IMC between different muscle groups within each synergy (high–all, high–high, and high–low, Fig. 5) using Friedman’s ANOVA; this analysis was constrained within the most relevant frequency bands (7–20 Hz, Fig. 4) around the observed peak.

Individual (white overlaid bars) and cluster means of normalised synergies (greyscale bars) for ramp and hold phases. Five different clusters (S1–S5) of synergies were calculated according to the maximum number of identified synergies across participants. Letters above the white overlaid bars correspond to each participant, as in Table 1. Muscles are labelled in an abbreviated form: superior (ST) and middle trapezius (MT), infraspinatus (Inf), teres minor (TM), serratus anterior (SA), anterior (AD), middle (MD), and posterior deltoid (PDel), pectoralis major (PM), short (BS) and long (BL) heads of biceps brachii, long (TL) and lateral (Tlat) heads of triceps brachii, brachioradialis (Braq), extensor carpi radialis (ECR), and flexor carpi radialis (FCR)