Participants
A total of 43 healthy, right-handed adult volunteers participated in the study (7 male and 36 female: mean age (± SD) was 22.9 ± 4.4 years). Handedness was assessed using the Edinburgh Handedness Inventory62. None of the participants had a history of neurological or psychiatric diseases. All participants provided written informed consent for participation in the experiment. The study was conducted according to the Declaration of Helsinki and was approved by the Ethical Committee of the National Institute for Physiological Sciences, Japan.
Experimental design
We carried out MRS-fMRI experiments using a 7T MRI scanner (MAGNETOM 7T, Siemens Healthineers, Erlangen, Germany) with a 32-channel receiving head coil and a single-channel transmitting coil (Nova Medical Inc., MA, USA). All participants underwent resting-state fMRI and MRS scans before and after the motor sequence learning tasks, as well as one MRS and four fMRI scans during motor sequence learning tasks in the task session (Fig. 1A). Dielectric pads (CaTiO3)63 were placed around each participant’s head while scanning at 7T to improve the B1 transmit field inhomogeneity. All scans were performed within the Specific Absorption Rate (SAR) limit of the normal operation mode.
Motor sequence learning task
Thirty participants were asked to perform pre-determined five-digit sequences “4-1-3-2-4” (n = 17) or “2-3-1-4-2” (n = 13) as quickly and accurately as possible in the MRI scanner (Fig. 1B)9–11. Additionally, 13 participants were asked to perform 120 different sequences to assess the non-specific learning as the control condition. The sequence “4-1-3-2-4” corresponds to “index-little-middle-ring-index.” The motor sequence task consisted of six 30-s tapping epochs followed by 30-s rest epochs that were repeated five times (Fig. 1B). The visual feedback signals were displayed using a projector (Optoma EH503; Optoma Inc., Fremont, CA, USA) with a lens (APO 50-500 mm F4.5-6.3 DG OS HSM; SIGMA, Kanagawa, Japan) on a screen viewed by the participants via a mirror mounted to the receiving head coil. Response time was measured using Presentation software version 16.4 (Neurobehavioral Systems, NY, USA; RRID: SCR_002521). The rest epoch started with the appearance of the instruction “Rest” on the screen for 500 ms, followed by a 500-ms presentation of four blue circles aligned within an equally spaced horizontal array. The instruction “Task” appeared for 2 s at the end of the rest epoch as a signal to the participants to retrieve motor sequences and prepare for their execution (Fig. 1B). The task epoch began with four closed white circles presented for 500 ms, which changed into open circles. During the task epoch, participants tapped the button box (Current Design, Philadelphia, USA) according to the sequence shown at the top of the screen (i.e., “4-1-3-2-4”). Visual feedback of correct tapping was provided by filling the white circle corresponding to the tapped finger. When the participant provided an incorrect response, the visual feedback signal remained at the previous position until the correct button was tapped. Task performance was measured using transition time (TT), defined as the average time between two correct button responses per epoch. The performance improvement was calculated using the following equation:
where TT1 indicates the transition time at block 1 and TT5 indicates the transition time at block 5. The task performance data were analyzed using repeated-measures analysis of variances (ANOVA), with block as a factor, performed using the Statistical Package for the Social Sciences software version 25.0.0 (SPSS, IBM Corp., NY, USA; RRID: SCR_002865). Two participants were excluded due to a statistical outlier in the TT values (> 2 SD).
Structural data acquisition
Three dimensional T1-weighted (T1w) images were acquired for anatomical reference (Magnetization Prepared Rapid Acquisition Gradient Echo [MPRAGE]64, TR/TE = 3,000/3.08 ms; inversion time [TI] = 1,200 ms; field of view = 240×225 mm2; matrix size = 320×320; slice thickness = 0.75 mm; 224 slices; generalized auto-calibrating partially parallel acquisitions [GRAPPA]65 acceleration factor = 3; bandwidth = 230 Hz/Px; flip angle = 14°; acquisition time = 4 min 50 s).
fMRI data acquisition
fMRI images were acquired before, during, and after the motor sequence learning tasks using a multiband gradient-echo echo-planar imaging sequence66. The scan parameters were set as per the human connectome project (HCP) 7T protocol67 (TR/TE = 1,000/22.2 ms; field of view = 208×208 mm2; matrix size = 130×130; slice thickness = 1.6 mm; 85 slices; multi-band/GRAPPA acceleration factor = 5/2; bandwidth = 1,924 Hz/Px; flip angle = 45°). The spin echo field map was acquired68 (TR/TE = 3,000/60 ms; field of view = 208×208 mm2; matrix size = 130×130; slice thickness = 1.6 mm; 85 slices; multi-band/GRAPPA acceleration factor = 5/2; bandwidth = 1,924 Hz/Px; flip angle = 180°; acquisition time = 1 min 26 s). A B1 transmit field map in the center of the brain, around the slice of the M1 hand knob area, was acquired for each participant to optimize the input power for accurately producing a 90° pulse for all fMRI scans. In particular, participants were instructed to keep their eyes open while viewing a fixation cross and to avoid having any specific thoughts or falling asleep during resting-state fMRI scans.
MRS data acquisition
A 2×2×2 cm3 volume of interest was centered over the right M1 hand knob area (Fig. 2A), without dura, on T1w MPRAGE images. The hand knob area was identified using fMRI during a sequential finger opposition task with the left hand (TR/TE = 1,000/24 ms; field of view = 192×192 mm2; matrix = 96×96; slice thickness = 2 mm; 20 slices; GRAPPA acceleration factor = 2; bandwidth = 2,170 Hz/Px; flip angle = 45°; acquisition time = 3 min 30 s). Ultra-short TE MRS data were acquired before, during, and after the motor sequence learning task using the STimulated Echo Acquisition Mode (STEAM) sequence (TR/TE = 5,000/5.68 ms; mixing time = 40 ms; vector size = 2,048; bandwidth = 4,000 Hz/Px; average = 64) with VAriable Power RF pulses with Optimized Relaxation delays (VAPOR) water suppression25,69. The STEAM sequence was combined with outer volume suppression to improve localization performance. A 4-average water reference signal was acquired for eddy current correction and absolute quantification of the metabolites. Before data acquisition, all first- and second-order shim terms were automatically adjusted with the fast automatic shim technique using echo-planar signal readout for mapping along with projections (FASTMAP)70,71. In addition, B1 transmit field strength for localization pulses and VAPOR water suppression was adjusted for individual participants.
HCP-style structural data acquisition with 3T MRI and preprocessing
In addition to the MRS-fMRI data acquisition using 7T MRI, the HCP-style structural data of all participants were obtained using a 3T MRI scanner (Magnetom Verio, Siemens Healthcare, Erlangen, Germany) with a 32-channel receiving head coil (Siemens Healthcare, Erlangen, Germany). The obtained 3T MRI data were utilized to correct the geometric distortion of the 7T MR data72 (see below, fMRI preprocessing). Scan parameters were as per the HCP 3T protocol with minor modifications73. Three-dimensional T1w images were acquired (MPRAGE64, TR/TI/TE = 2,400/1,060/2.24 ms; field of view = 256×240 mm2; matrix size = 320×320; slice thickness = 0.8 mm; 224 slices; GRAPPA acceleration factor = 2; bandwidth = 210 Hz/Px; flip angle = 8°; acquisition time = 6 min 38 s; measurement = 2). Three-dimensional T2 weighted (T2w) images were acquired (Sampling Perfection with Application optimized Contrast using different angle Evolutions [SPACE]74, TR/TE = 3,200/560 ms; field of view = 256×240 mm2; matrix size = 320×320; slice thickness = 0.8 mm; 224 slices; GRAPPA acceleration factor = 2; bandwidth = 744 Hz/Px; turbo factor = 167; acquisition time = 6 min; measurement = 2). All data were processed using the structural pipeline (PreFreeSurfer, FreeSurfer, and PostFreeSurfer) of the minimal HCP preprocessing pipeline version 4.0.0-alpha.5, including the following steps: gradient magnetic field nonlinearity distortion correction, T2w images to T1w image registration, and Montreal Neurologic Institute (MNI) volume registration73.
MRS data analysis
Raw MRS data were post-processed using MATLAB R2018a (The MathWorks, Inc., MA, USA; RRID: SCR_001622). Motion-corrupted data were removed to improve the spectral quality. To quantify the proportion of gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF) fractions in the volume of interest, segmentation in SPM was applied to the T1w MPRAGE images. Eddy current correction and frequency correction were performed using a water reference scan, and the zero and first-order phases of the array coil were aligned using the cross-correlation method of MRspa (RRID: SCR_017292). Subsequently, LCModel version 6.3-1N75,76 (Stephen Provencher, Inc., ON, Canada; RRID: SCR_014455) analysis was used to quantify the concentration of neurochemicals within the chemical shift range of 0.5 to 4.1 ppm76. Other parameters in the LCModel were as reported previously77. The concentrations of GABA and Glu were normalized to that of total creatine (tCr). The change in glutamate to GABA ratio (Glu/GABA) after the motor sequence learning task was calculated using the following equation:
where Glu/GABApre and Glu/GABApost indicate the Glu/GABA ratio at pre-task and post-task, respectively. The distribution of GABA and Glu concentrations was visualized using the RainCloudPlots Python-script78 (https://github.com/RainCloudPlots/RainCloudPlots). A repeated-measures ANOVA was carried out in SPSS, with the concentrations of GABA and Glu at different time points (pre-task, during-task, and post-task) as a factor. The Cramer–Rao lower bounds (CRLB), water linewidth at FWHM, and signal-to-noise ratio (SNR) were used for the quality control of spectra76. The CRLB and SNR were calculated using LCModel, and water linewidth was obtained by fitting to the additional water spectrum using MATLAB. Data were excluded when CRLB > 15 % (n = 1), linewidth > 19 Hz (n = 1), or SNR < 30. A repeated-measures ANOVA was performed on the CRLB, water linewidth, and SNR time points (pre-task, during-task, and post-task) with a within-subjects factor using SPSS.
fMRI preprocessing
All fMRI data were processed using the functional pipeline (fMRIVolume) of the minimal HCP preprocessing pipeline72. This pipeline included the following steps: motion correction, gradient magnetic field nonlinearity distortion correction, field map-based distortion correction (Topup)68, nonlinear registration into 3T MNI structure data, and grand-mean intensity normalization. Finally, volume-based smoothing with a 5-mm full width at half maximum (FWHM) Gaussian kernel was applied.
Task fMRI data analysis
Task fMRI data analysis was performed using Statistical Parametric Mapping (SPM12; RRID: SCR_007037) in MATLAB R2018a. A general linear model (GLM) was fitted to the fMRI data for each participant79,80. The fMRI time series for preparation phases 2 s before task execution and execution phases were modeled with boxcar functions convolved with the canonical hemodynamic response function. Each block consisted of six execution-related and preparation-related regressors. The design orthogonality between the execution and preparation phases was −0.0137 ± 0.054 for block 1, −0.0141 ± 0.054 for block 2, −0.0137 ± 0.054 for block 3, and −0.0139 ± 0.054 for block 4 (mean ± SD). Temporal high-pass filtering with a cutoff frequency of 1/128 Hz was applied. Using a first-order autoregressive model, the serial autocorrelation was estimated from the pooled active voxels with the restricted maximum likelihood procedure and subsequently used to whiten the data81. Several nuisance covariates, including six head motion parameters and CSF time-series, were incorporated into the model. The parameter estimates for each execution-related and preparation-related regressors were assessed using constant and predefined linear contrasts. Increasing contrast vectors were defined numerically as an increment of one per block, maintaining the mean equal to zero.
For group-level analysis of task fMRI data, one-sample t-tests of participants’ contrast images were performed82. The resulting set of voxel values for each contrast constituted the SPM{t}. We calculated the T-score of linear increment in preparation-related activity in right M1 in non-specific learning. The statistical threshold was set at p < 0.05, FWE-corrected at the voxel-level83, unless otherwise specified.
Anatomical labeling and visualization
MRIcron (RRID: SCR_008264) was used to display fMRI activation maps on a standard brain image. The Automated Anatomical Labeling atlas was used for anatomical labeling84.
Resting-state fMRI data analysis
Resting-state functional connectivity analysis was conducted using the CONN toolbox version 17 in SPM1285 (RRID: SCR_009550). An anatomical component-based noise correction method (aCompCor)86 was applied to remove the five components of signals from WM, CSF, and residual head motion-related signals through linear regression. A temporal bandpass filtering of 0.008–0.090 Hz was applied.
Seed-to-voxel correlation analysis was performed at the individual level. We selected the preparation-related increased voxels in M1 (MNI: x = 36, y = −25, z = 51), determined in the second-level analysis of task fMRI (FWE voxel-level corrected p < 0.05), as a seed region of interest (ROI) (Fig. 6A). An individual seed-based functional connectivity map was obtained by computing Pearson’s correlation coefficients between the time-series from the M1 seed ROI and the time-series of all other voxels across the whole brain. Fisher’s r-to-z transformation was used to convert the correlation coefficients into z-scores. M1-seeded functional connectivity changes were integrated using the following equation using AFNI version 18.1.32. (RRID: SCR_005927):
where Connectivitypre and Connectivitypost are the pre-task and post-task functional connectivity values, respectively.
We calculated the changes in functional connectivity within ROIs of the SMN and FPN defined from CONN’s ICA analyses of the HCP dataset of 497 individuals. The SMN includes the supplementary motor cortex and bilateral sensorimotor cortex, whereas the FPN consists of the bilateral LPFC and PPC. The correlations between Glu/GABA changes within M1 and M1 seed-based functional connectivity changes were analyzed using linear regression analysis.