2.4. Data analysis

Data were analyzed using Matlab (version 2017b, The Mathworks, Natwick, MA). For the (f)MRI processing we employed FreeSurfer (version 6.0.0, Laboratory for Computational Neuroimaging, Boston, United States, see Fischl (2012)) the FMRIB (Functional MRI of the Brain) Software Library [version 5.0, Analysis Group, Oxford, UK, see Jenkinson, Beckmann, Behrens, Woolrich, and Smith (2012)], connectome workbench (version 1.2.3, WU‐Minn HCP Consortium, USA, see Marcus et al. (2011)) and the preprocessing scripts of the human connectome project (version 3.24.0, WU‐Minn HCP Consortium, USA, Glasser et al. (2016, 2013)).

We quantified performance via the strength of frequency locking between the visual cue and the force produced. As corresponding measure we used the normalized spectral overlap (Daffertshofer, Peper, & Beek, 2000), that is, the similarity ψxy between the power spectrum of the cueing signal, Px, and that of the produced force, Py, after rescaling the frequency axis by the factor ρ=p:q. This measure reads

It is bounded to the interval [0,1] with 1 indicating maximum similarity. We Fisher‐transformed this value prior to statistical evaluation to stabilized normality.

Statistical analysis of the behavioral data was performed with SPSS (IBM Corp. Released 2015. IBM SPSS Statistics for Macintosh, Version 23.0 Armonk, NY, United States: IBM Corp.). Frequency locking values per trial for both groups were subjected to an ANOVA with repeated measures. Within‐subject factors were the 10 task trials and the retention test. Age group was included as between‐subject factor. The significance threshold was set to α = 0.05. Sphericity was tested using Mauchly's test and, when violated, we applied a Greenhouse–Geisser's correction. Post hoc t‐tests were performed whenever a main effect of Age or Trial was detected to evaluate effects of motor learning. In the first case we applied an independent samples t‐test between age groups and in the latter case a dependent t‐test on trials.

EEG data were segmented with BrainVision Analyzer software (Brain Products) into three parts: the experimental learning task, pre‐ and post‐motor task resting state. We removed MR‐scanner artifacts as detailed in Supplementary Material S1, down‐sampled to 512 Hz after anti‐aliasing filtering, and finally band‐pass filtered between 0.1 and 100 Hz. Subsequent EEG preprocessing consisted of removing and interpolating bad channels, removing the ballisto‐cardiogram (pulsatile blood movement causing body and electrode movement inside the scanner) and excessive eye blinks as well as movement artifacts using independent component analysis (ICA) (Hyvarinen, 1999). Since ballistocardiograph artifacts can generally be expected to covary with the ECG, we combined these two signals with the EEG channels, conducted principal component analysis (PCA) rather than ICA and removed all components that were dominated by the ECG.

MRI data were preprocessed following the human connectome project preprocessing pipeline (Glasser et al., 2013, 2016). Structural T1 and T2 weighted images were aligned, bias field corrected, skull stripped and nonlinearly registered to MNI space. We applied FreeSurfer's recon‐all to reconstruct cortical gray matter surfaces and a subcortical gray matter volume segmentation. Myelin maps across the surface were created by taking the ratio of T1w/T2w from voxel intensities inside the cortical gray matter (Glasser & Van Essen, 2011). The first five fMRI scans were removed, and motion correction was performed by aligning the first image to the series. In this stage we also checked for potential scanner artifacts due to heating. Motion parameters served as regressors when cleaning the fMRI data from remaining motion‐related artifacts via ICA. The fMRI data were corrected for echo‐planar imaging distortion using a gradient echo field map, aligned to the T1w image as well as MNI space and bias field corrected. The fMRI data were converted into the CIFTI file format, with time series of voxels inside the cortical gray matter ribbon being mapped onto cortical surface vertices and the subcortical gray matter being resampled onto a standard volume mesh. Next, fMRI data were cleaned using FSL's FIX tool (Griffanti et al., 2014; Salimi‐Khorshidi et al., 2014). In brief, using a pre‐trained classifier automatically we labeled independent components as either neural activity or artifacts. The artifact components as well as motion parameters served as regressors to remove the corresponding co‐variate from the fMRI. Finally, the FIX classifier was trained on hand‐labeled data from all the 23 subjects of this study. Subjects’ individual cortical surfaces were registered to a parcellation template (Glasser et al., 2016) of the human connectome project following a multimodal surface registration approach (Robinson et al., 2014, 2018). The features used for the surface registration were cortical folding, myelin maps, resting state network locations and visuotopic maps. The pre‐ and post‐task fMRI data were used to identify subject individual resting state network locations and visuotopic maps. Aligned fMRI data entered statistical analysis.

Source localization of the β‐rhythm was realized using dynamical imaging of coherent sources beamformers (Gross et al., 2005) with tissue segments (skin, skull, cortical spinal fluid, gray and white matter) determined through finite‐element‐modeling. For both we employed the open‐source Fieldtrip toolbox and used their template MRI (Oostenveld, Fries, Maris, & Schoffelen, 2011).

Upon visual inspection of motor‐related potentials, we selected the maximum increase in hand force as central events for the subsequent EEG analysis because the surround epochs reveal the maximum difference between β‐synchronization and de‐synchronization (Figure 3). Epochs of ±400 ms before and after these events served as contrast for subsequent statistics. To identify significance of beamformer power, we followed a Monte Carlo approach with cluster‐based test statistics for both groups (with a significance threshold of α _cluster = 0.01, Nichols and Holmes (2002)). Cluster‐level statistics were determined for the separate groups as the sum of t‐values per cluster. Probabilities were determined by collecting trials of pre‐ and post‐event intervals and test statistics were computed on randomly chosen partitions. These steps were repeated 8,192 times to construct the permutation density of the test statistics, which allowed for using a dependent samples t‐test between pre‐ and post‐event epochs (significance threshold α = 0.05). Following the same approach but using an independent t‐test, we tested for significant differences between groups (α < .05). While we focused on β‐band activity, we would like to note that we provide all the corresponding findings for the α‐frequency band (8–14 Hz) as Supplementary Material S2. Beamformer outcomes were parcellated using the SPM anatomical atlas (Eickhoff et al., 2005). ROIs were defined via significant outcomes of the aforementioned tests and corresponding virtual sensors were defined via the mean spatial filter over each ROIs onto which the EEG data were projected after band‐pass filtering in the frequency band under study (β‐band in the main text, α‐band in Supplementary Material S2). The motor‐event‐related β‐amplitude modulation per trial and for the two resting state blocks (based on randomly placed virtual events) was evaluated as described in Houweling, Beek, and Daffertshofer (2010); see also Neuper, Wörtz, and Pfurtscheller (2006). In a nutshell, we normalized the β‐band time‐series of the virtual sensor to baseline, that is, the first resting state recording, computed the mean Hilbert amplitude and performed a time‐locked averaging (David, Kilner, & Friston, 2006) over the aforementioned ±400 ms epochs. We further evaluated with a generalized estimating equation whether the motor event‐related β‐modulation correlated with the behavioral performance per trial.

Illustration event definition. Left panel: Pressure inside the air‐filled rubber bulb; maximum changes of hand force were defined as central events in the motor event‐related beamformer approach. Epochs of ±400 ms before and after these events served as contrast. Right panel: Corresponding event‐related synchronization (ERS) and desynchronization (ERD) in the β‐band of the EEG signal (electrode C3)

For both groups we also computed grand‐average source‐level β‐power across learning vs. pre‐learning rest, also referred to as task‐related power (Andres et al., 1999; Gerloff et al., 1998; Serrien, Cassidy, & Brown, 2003; Toro et al., 1994). With this we tested whether β‐power during motor task execution differed between groups, again using a voxel‐wise permutation tests with an independent samples t‐test.

We analyzed the fMRI data using general linear modeling (GLM) as outlined below. When combining fMRI with EEG, the source localized EEG served as an additional regressor. For this, we estimated the instantaneous Hilbert amplitude in the frequency‐band of interest. For every scan epoch (700 in total) we utilized the mean β‐amplitude 800 ms around every motor event and averaged over the events per scan. To accommodate for delays and dispersions in the BOLD‐responses, the regressors were convolved with the double gamma hemodynamic response function (HRF) (Grinband, Steffener, Razlighi, & Stern, 2017) using the FMRIB Software Library; see Figure 4.

Computation of the regressors. Illustration of the computation of the beta amplitude and the resulting regressors. From the source‐localized EEG (top panel), the Hilbert amplitude was determined and the mean β‐amplitude 800 ms around the motor events was averaged for every scan (middle panel, blue) and convolved with the hemodynamic response function (HRF) to determine the HRF‐convolved regressor (bottom panel, blue). The red curves show the regressor of the task design (in the bottom panel after HRF‐convolution)

Prior to statistical analysis, the fMRI data were spatially smoothed both on the surface and in the subcortical volume using a 4 mm FWHM Gaussian Kernel to suppress spatial noise and to increase the signal‐to‐noise ratio. To remove low‐frequency noise, the data were also filtered with a Gaussian‐weighted linear high‐pass filter with a cutoff of 135 s. On single‐subject level, we fitted the GLM: Y=βkXk+ε, using the grayordinates‐wise BOLD data Y in CIFTI format, that is, the time series from cortical vertices and subcortical voxels (Friston et al., 1995). Here, βkXk denotes the design matrix and ε the residual error. To build the design matrix we used the task paradigm and combined it with the EEG source‐localized β‐amplitude (Figure 4, bottom panel). The latter stemmed either from contralateral M1 (left area 4a) or from frontal cortex (left area 6 ~ premotor area), both according to the SPM atlas (Eickhoff et al., 2005). We contrasted each regressor against baseline and both regressors against each other to estimate statistical effects of interest. When both regressors—task and EEG—are placed in a single GLM, the parameter estimates will reflect the activation in BOLD of one regressor adjusted for the effect of the other. By consequence, the variance explained by both regressors will be removed, isolating the effect of the EEG regressor. Note that by conducting orthogonalization, the shared information of the two regressors is attributed to the regressor that is not orthogonalized (Mumford, Poline, & Poldrack, 2015); hence, orthogonalization was not applied. Based on these single‐subject contrasts, a mixed‐effects group‐level analysis was performed using a paired or unpaired two group t‐test, to contrast within‐ and between‐group activation, respectively. The latter contrasts were masked to reveal only those areas that deemed significant using the single‐group activation contrast. Finally, we thresholded the statistical maps using a false discovery rate of q = 0.01 (Genovese, Lazar, & Nichols, 2002).