Abstract
The mechanisms subserving motor skill acquisition and learning in the intact human brain are not fully understood. Previous studies in animals have demonstrated a causal relationship between motor learning and structural rearrangements of synaptic connections, raising the question of whether neurite-specific changes are also observable in humans. Here, we use advanced diffusion magnetic resonance imaging (MRI), sensitive to dendritic and axonal processes, to investigate neuroplasticity in response to long-term motor learning. We recruited healthy male and female human participants (age range 19–29) who learned a challenging dynamic balancing task (DBT) over four consecutive weeks. Diffusion MRI signals were fitted using Neurite Orientation Dispersion and Density Imaging (NODDI), a theory-driven biophysical model of diffusion, yielding measures of tissue volume, neurite density and the organizational complexity of neurites. While NODDI indices were unchanged and reliable during the control period, neurite orientation dispersion increased significantly during the learning period mainly in primary sensorimotor, prefrontal, premotor, supplementary, and cingulate motor areas. Importantly, reorganization of cortical microstructure during the learning phase predicted concurrent behavioral changes, whereas there was no relationship between microstructural changes during the control phase and learning. Changes in neurite complexity were independent of alterations in tissue density, cortical thickness, and intracortical myelin. Our results are in line with the notion that structural modulation of neurites is a key mechanism supporting complex motor learning in humans.
SIGNIFICANCE STATEMENT The structural correlates of motor learning in the human brain are not fully understood. Results from animal studies suggest that synaptic remodeling (e.g., reorganization of dendritic spines) in sensorimotor-related brain areas is a crucial mechanism for the formation of motor memory. Using state-of-the-art diffusion magnetic resonance imaging (MRI), we found a behaviorally relevant increase in the organizational complexity of neocortical microstructure, mainly in primary sensorimotor, prefrontal, premotor, supplementary, and cingulate motor regions, following training of a challenging dynamic balancing task (DBT). Follow-up analyses suggested structural modulation of synapses as a plausible mechanism driving this increase, while colocalized changes in cortical thickness, tissue density, and intracortical myelin could not be detected. These results advance our knowledge about the neurobiological basis of motor learning in humans.
Introduction
The production of well coordinated movements is among the most remarkable abilities of humans, and understanding the neural basis underlying motor learning is a core topic of the neurosciences (Dayan and Cohen, 2011). Research in animals has shown that motor learning is accompanied by changes in the size, shape and volume of the brain (Rosenzweig and Bennett, 1996; Markham and Greenough, 2004; Zatorre et al., 2012). At a finer scale, it is assumed that several neuronal and non-neuronal processes contribute to these macrostructural alterations (Zatorre et al., 2012). In the gray matter (GM), synaptic remodeling is a well established mechanism underlying motor learning (Xu et al., 2009; Hayashi-Takagi et al., 2015), whereas myelin-related changes have been conclusively shown to occur in white matter (WM) tissue (Sampaio-Baptista et al., 2013; McKenzie et al., 2014; Bacmeister et al., 2020). Of note, neuroplastic processes in GM and WM were demonstrated to be highly dynamic with transient retraction and expansion tendencies (Xu et al., 2009; Wenger et al., 2017; Bacmeister et al., 2020), requiring highly specific methods to detect them.
Because of methodological limitations in in vivo sampling of brain tissue, it has long been impossible to investigate whether the learning-related plasticity processes demonstrated in animals, such as changes in neurite structures (i.e., dendrites and axons), also apply in humans. Over the past two decades, longitudinal studies using magnetic resonance imaging (MRI) with sensitized contrast to the mobility of water yielded promising achievements in this respect (for review, see Zatorre et al., 2012; Tardif et al., 2016; Lerch et al., 2017). Usually, these studies employed diffusion sensitization at low b-values (≤1000 s/mm−2) and fitted the measured signal to a diffusion tensor (DTI; Basser and Pierpaoli, 1996). Previous studies of motor learning based on this approach have shown changes in diffusivity and anisotropy in both GM (Sagi et al., 2012) and WM (Scholz et al., 2009).
Albeit sensitive to change, tensor-derived metrics are biologically ambiguous because they reflect a voxel average of the degree of hindrance or restriction experienced by the diffusing water (Jones et al., 2013; Lerch et al., 2017). Key neuroplastic processes related to the brain's neurite architecture or myelin remodeling can therefore not be directly inferred from DTI. Moreover, DTI is not well suited to deal with highly dispersed and complex tissue microstructure present in the vast majority of WM, but especially in GM voxels (Jones et al., 2013; Tardif et al., 2016; Lerch et al., 2017; Assaf, 2019).
Recent advances in diffusion MRI do now offer the possibility to characterize key neuroplastic processes and their dynamics in unprecedented specificity, while possibly also offering higher sensitivity to detect neuroplastic changes and brain-behavior-correlations (Kamiya et al., 2020). For instance, multicompartment models like Neurite Orientation Dispersion and Density Imaging (NODDI) can cope with complex fiber arrangements and can differentiate between intraneurite restricted (e.g., axons and dendrites) and extraneurite hindered (e.g., cell bodies, glial cells) diffusion. Extensive validation studies have shown a high correspondence between NODDI-based estimates of neurite density (Grussu et al., 2017; Wang et al., 2019; Gong et al., 2020) and orientation dispersion (Grussu et al., 2017; Mollink et al., 2017; Sato et al., 2017; Schilling et al., 2018) with their histologic counterparts.
Here, we investigate whether key neuroplastic mechanisms present in animals can also be observed in humans. Specifically, we aim to identify brain areas showing correlated changes (C. Thomas and Baker, 2013) between NODDI-derived microstructural indices and concurrent motor learning. To achieve this, we use an established dynamic balancing task (DBT; Wulf et al., 2001) as a learning paradigm, whose acquisition process is expected to result in structural reorganization, particularly in sensorimotor, supplementary motor and prefrontal areas and associated nerve fiber tracts, according to previous study results (Taubert et al., 2010, 2016; Lehmann et al., 2020, 2022). To narrow down possible neurobiological correlates underlying plasticity (Tardif et al., 2016), we use information from multiple state-of-the-art MRI modalities.
Materials and Methods
Participants and experimental design
We used a controlled within-subject design (see below) to investigate whether learning the DBT would change the microstructure of the brain's GM and WM. Twenty-four cognitively healthy adults [21 ♂, 3 ♀; age: (arithmetic) mean (M) = 22.21, SD = 3.05, range 19–29; body mass index (BMI): M = 23.55, SD = 2.52, range 18.99–29.40] with no history of neurologic, psychiatric or systemic diseases were included in this trial. Further exclusion criteria were contraindications to MRI, body mass index (BMI) > 30 kg/cm2, a history of neuropsychiatric diseases, left-handedness, self-reported physical activity of >2 h/week, prior experience with the DBT, and past or present performance-oriented participation in endurance and/or coordinative-demanding sports. The sample was a subset of the subjects that participated in our previous reliability studies (Lehmann et al., 2021; Aye et al., 2022). The study was performed in accordance with the ethical standards as laid down in the 1964 Declaration of Helsinki and its later amendments. Approval was granted by the Ethics Committee of Otto von Guericke University Magdeburg (approval number 106/98). Written informed consent was obtained from all individual participants included in the study.
The experimental design consisted of a control phase and a learning phase of equal duration. MRI measurements were conducted three times at four-week intervals. In the first four-week period, we assessed the reliability of our MR imaging protocol. After this control period subjects engaged in four consecutive weeks of learning the DBT. At the end of the learning period subjects received MRI scans identical to those in MRI1 and MRI2. To avoid acute influences of motor training on MRI (e.g., increased functional activity or blood flow), training sessions (TSs) and MRI measurements were performed at least 24 h apart. Note that all subjects participated in both phases of the experiment, such that each subject serves as its own control (cf. A.G. Thomas et al., 2009), reducing the influence of potential confounding variables. In other words, the MRI signal changes during the learning phase (MRI2–MRI3) can be related to the changes during the control phase (MRI1–MRI2), which functions as the baseline. The within-subjects design also ensures optimal statistical power because of the removal of variance associated with between-subject differences (Poldrack, 2000; C. Thomas and Baker, 2013; Szucs and Ioannidis, 2020).
Whole-body dynamic balancing task (DBT) and quantification of motor learning
After the four-week scan-rescan control interval was completed, subjects commenced with learning the DBT over four weeks, with two training sessions (TSs) per week at least 24 h apart. DBT training was performed on an unstable seesaw-like platform (stability platform, model 16030, Lafayette Instrument) moveable in a medio-lateral direction with a maximum deviation of 26° on either side of the pivot. Standing with both feet on the platform, subjects were instructed to keep the board in a horizontal position for as long as possible during each trial. Every TS consisted of 15 trials with a duration of 30 s each and an intertrial break of 90 s to avoid fatigue. The behavioral outcome measure was the time (millisecond timer) in which subjects kept the platform within a deviation range of ±3° from the horizontal (time balancing, BAL). After each trial, subjects received verbal feedback on how long the board was in the target zone (knowledge of results), while no feedback was provided on strategy or other aspects of the task (discovery learning approach; Orrell et al., 2006). During task execution, participants' attention was directed to a fixation cross affixed to the wall in front of them (external focus of attention).
To get a summary measure of learning performance over eight training sessions, we first averaged the 15 BAL values belonging to each training session. Second, we fitted a general power function
MR image acquisition
MRI data were acquired on a 3T MAGNETOM Prisma system (Siemens Healthcare) using a 64-channel head coil. The same measurement protocol was used for each volunteer and scanning session. Subjects were asked to relax, keep their mind free of any thoughts, and to move as little as possible. A pillow was placed surrounding the sides and the back of the head to minimize head motion and within-subject as well as between-subject differences in positioning.
Whole-brain diffusion MR images were obtained with a monopolar single-shot spin echo echo planar image (EPI) sequence: echo time (TE) = 74 ms; repetition time (TR) = 4970 ms; flip angle α = 90°; parallel GRAPPA acceleration factor = 2, matrix: 130 × 130; field-of-view (FOV) = 208 × 208 mm; nominal spatial resolution = 1.6 × 1.6 × 1.6 mm; multiband acceleration factor = 2; phase-encoding direction: anterior ≫ posterior). The sampling scheme was designed according to Caruyer et al. (2013; http://www.emmanuelcaruyer.com/q-space-sampling.php) and consisted of 228 isotropically distributed diffusion sensitization directions including three diffusion shells with b-values of 1000 s/mm2 (38 directions), 2000 s/mm2 (76 directions), and 3000 s/mm2 (114 directions). Note that diffusion sensitization in this protocol is considered strong enough to sensitize to pools of water with restricted diffusion (especially neuritic elements) and to reduce signal contributions of the faster moving extraneurite water (Jones et al., 2013; Weiskopf et al., 2021). Fourteen images without diffusion weighting (b = 0 s/mm2) were collected interleaved throughout the acquisition. To generate appropriate fieldmaps to correct for susceptibility-induced distortions, nine b = 0 s/mm2 images with reversed phase encoding (posterior ≫ anterior) were also acquired. The total scan duration was 22 min 31 s.
To quantify cortical macrostructure, we used a standard T1w three-dimensional magnetization-prepared rapid gradient echo sequence (Mugler and Brookeman, 1990) based on which cortical thickness was estimated. The imaging parameters used were as follows: 240 sagittal slices, inversion time, TI = 1100 ms; TR = 2600 ms; TE = 5.18 ms; readout pulse flip angle, α = 7°; parallel GRAPPA acceleration factor = 2; acquisition matrix = 320 × 320; FOV = 256 × 256 mm; nominal spatial resolution = 0.8 × 0.8 × 0.8 mm; scan duration = 7 min 25 s.
Not least, to test myelin-related hypotheses (Tardif et al., 2016; Lerch et al., 2017; Natu et al., 2019; Lazari and Lipp, 2021), we used magnetization transfer saturation (MTsat) maps (Helms et al., 2008) based on MT-weighted images. MT-weighted images were acquired (as part of a multiparameter mapping protocol) with a multiecho FLASH scan with TR = 37 ms/α = 7°. Multiple gradient echoes were acquired with alternating readout polarity at six equidistant TEs between 2.46 and 14.76 ms. Other acquisition parameters were: 224 sagittal slices, nominal spatial resolution = 0.8 × 0.8 × 0.8 mm, FOV = 230 × 230 mm. Transmit and receive field correction acquisition was done before with the following imaging parameters: 56 sagittal slices, FOV = 230 × 230 mm, TR = 4.1 ms, TE = 1.98 ms for B1− and TR = 2000 ms, TE1 = TE2 = 14 ms, 24 sagittal slices, slice thickness = 5 mm, α = 90°, 120°, 60°, 135°, 45° for the RF map which was used for the B1+ correction. MT-weighted images were further preprocessed as outlined in Aye et al. (2022).
Preprocessing of diffusion-weighted images
Preprocessing of diffusion-weighted images followed the standard FSL diffusion processing pipeline (Smith et al., 2004). This pipeline starts with identifying susceptibility-induced distortions based on reversed phase-encoding b = 0 s/mm2 volumes using the topup tool (Andersson et al., 2003). Afterwards, corrections for eddy currents and subject movement were conducted using eddy (Andersson and Sotiropoulos, 2016), and the fieldmap that resulted from topup was applied. Realignment of images in the course of motion correction was accompanied by appropriate correction of gradient directions (Leemans and Jones, 2009).
Fitting of microstructural maps
In the next step, diffusion imaging data were fitted with NODDI (Zhang et al., 2012) and with the conventional diffusion tensor (Basser and Pierpaoli, 1996). Note that both modeling approaches were used in WM, because it has been suggested that the tensor-model might be particularly sensitive to change in tissue environments with approximately parallel fiber bundles (Tardif et al., 2016; Alexander et al., 2019; Kamiya et al., 2020). Unlike WM, GM exhibits an even more heterogeneous neural environment and shows above that a problematic contamination with CSF (Tardif et al., 2016; Assaf, 2019). From this it follows that several assumptions of the tensor model, such as the presence of one principal diffusion direction per voxel or single-exponential signal decay, are unlikely to be fulfilled in GM (Jones et al., 2013; Lerch et al., 2017; Fukutomi et al., 2019). On the contrary, NODDI is particularly well suited to assess neuroplastic change in tissue environments with highly dispersed and complex microstructure (e.g., crossing, fanning, or kissing fibers) because it allows for distributed fiber orientations within voxels (Tardif et al., 2016; Lerch et al., 2017; Assaf, 2019; Wang et al., 2019). Furthermore, NODDI-derived tissue indices are free from CSF contamination, making it particularly well suited for the study of plasticity in GM (Assaf, 2019; Fukutomi et al., 2019). Since NODDI is also considered to be more sensitive to microstructural change in GM than the tensor (Tavor et al., 2013; Vukovic et al., 2021), only NODDI maps were used to assess plasticity in this particular tissue class.
NODDI is a multicompartment model of diffusion aiming at the estimation of biologically interpretable microstructural tissue properties in each image voxel (Zhang et al., 2012). It assumes that the diffusion signal arises from the sum of three distinguishable tissue environments, namely intraneurite space, extraneurite space, and free water (for details of the model, see Zhang et al., 2012; Kraguljac et al., 2023). There are three microstructural maps emerging from this model that are used for statistical analyses: the isotropic volume fraction (ISO), the neurite density index (NDI), and the orientation dispersion index (ODI). ISO represents the signal from freely diffusing water or CSF (as opposed to tissue) in a voxel. Note that ISO and the tissue volume fraction, which comprises the intraneurite and extraneurite compartments, sum up to 1. This means that a higher (lower) value of ISO must be accompanied by a lower (higher) value of the tissue volume fraction. NDI is the proportion of tissue that is intended to represent contributions from axons, dendrites or neural processes (Zhang et al., 2012; Assaf, 2019). The relative volume fraction of the extraneurite compartment (1-NDI) is assumed to mainly comprise contributions from cellular membranes of somas and glial cells (Zhang et al., 2012). Not least, ODI is a metric based on the Watson distribution concentration parameter κ (Zhang et al., 2012). It describes the spatial configuration of the intraneurite space and represents how widely the orientations of the neurites spread out in space (Kraguljac et al., 2023).
In the present study, NODDI parameter maps were estimated from corrected multishell diffusion images (b = 0 s/mm2, b = 1000 s/mm2, b = 2000 s/mm2, and b = 3000 s/mm2). To fit the model in WM voxels, we used the NODDI MATLAB Toolbox v1.0.1 (http://nitrc.org/projects/noddi_toolbox) with default settings, implementing the model formulation of Zhang et al. (2012). Fukutomi et al. (2018) and Guerrero et al. (2019) have recently suggested that the fixed NODDI model parameter intrinsic free diffusivity (din), which reflects the rate of diffusion along individual neurites, needs to be adjusted in GM voxels. To this end, GM-optimized NODDI maps with din set to 1.1 × 10−3 mm2/s (instead of the WM optimum of 1.7 × 10−3 mm2/s) were calculated using the Accelerated Microstructure Imaging via Convex Optimization (AMICO) v1.5.2 toolbox running in python (Daducci et al., 2015).
To assess microstructural plasticity in WM, we also fitted a diffusion tensor (Basser and Pierpaoli, 1996) at each voxel of the preprocessed images (b = 0 s/mm2 and b = 1000 s/mm2 shells only) using FSL's dtifit. The diffusion indices fractional anisotropy (FA), mean diffusivity (MD), and radial diffusivity (λ⊥) were computed from the eigenvalues of the diffusion tensor with the standard formulas as described by Pierpaoli and Basser (1996).
Image normalization using tract-based and gray matter-based spatial statistics
To ensure optimal sensitivity of longitudinal exploratory whole-brain analyses, it is paramount to have accurate alignment of the microstructural maps within each subject as well as cross-subject voxel-based correspondence in standard (MNI152) space. To this end, we used skeleton-based registration approaches optimized for the conditions of WM and GM, respectively. The underlying idea of these approaches is to project diffusion imaging-derived microstructural maps to an alignment-invariant sample-specific skeleton which either represents the center of the main WM tracts (Smith et al., 2006) or the maximally probable GM voxels (Nazeri et al., 2015), as outlined below. Besides improving registration, skeletonization is also hoped to minimize the influence of partial-volume contamination of microstructural maps (Smith et al., 2006; Nazeri et al., 2015; Lerch et al., 2017).
For WM microstructural maps, we used a reliable (Madhyastha et al., 2014; Lehmann et al., 2021) and sensitive (Lehmann et al., 2020) TBSS-based (Smith et al., 2006) processing routine, starting with determining an unbiased midspace between the three FA images of each participant (Smith et al., 2002). In the next step, the native space FA images of each participant were linearly registered to the individual midspace and subsequently averaged to generate a within-subject FA template (Madhyastha et al., 2014; Lehmann et al., 2020, 2021). Afterwards, each subject's FA template was nonlinearly aligned to every other one to identify the most representative template of the sample (target). After warping each subject's template to the target, images were linearly registered to MNI152 space (FMRIB58 1-mm template) using affine transformation, and a group-average FA image was thinned and binarized with an FA-value of >0.25 (skeletonization). Finally, the previously created warp fields were applied to all midpoint-space registered NODDI/DTI maps of the three measurement points, and the aligned NODDI/DTI maps in MNI space were projected onto the groupwise skeleton.
For GM microstructural maps, we used a modified gray matter-based spatial statistics (GBSS) pipeline inspired by TBSS (Nazeri et al., 2015; bash scripts available at https://github.com/arash-n/GBSS). GBSS starts with a segmentation of brain tissues using the NODDI and DTI data. Specifically, the CSF volume fraction (fCSF) corresponds NODDI's ISO map, while the WM volume fraction (fWM) was derived from a two-tissue classification (WM/non-WM) of FA images using the Atropos segmentation tool implemented in ANTs (Avants et al., 2011). GM fraction maps were then calculated according to the formula:
Note that we used midspace-registered average images of each subject for the processing steps mentioned above. Next, fGM and diffusion microstructural maps were warped to MNI152 space using the registrations estimated in the course of TBSS. Afterwards, fGM maps in MNI space were averaged resulting in a mean GM image which was subsequently thinned and binarized (skeletonization). Only voxels with fGM > 0.5 in >50% of the subjects were retained. Note that this threshold was set to a slightly more liberal level than the one originally proposed by the authors (cf. Nazeri et al., 2015) to ensure an adequate representation especially of the motor regions on the skeleton, which are known to have a relatively thin cortex. Next, each subject's maximally probable local GM voxel perpendicular to the skeleton was identified. Finally, applying the estimated projections from the previous step, GM-optimized NODDI and MTsat maps (after boundary-based intermodal registration to native diffusion space, Greve and Fischl, 2009) of each subject and measurement point were registered onto the skeleton.
Quality assurance and reliability assessment
We calculated a recently proposed index of diffusion MR image quality, the temporal signal-to-noise ratio (tSNR), from the preprocessed and brain-extracted DWI data. TSNR is a metric that captures both scanner-induced and subject-induced imaging artifacts and allows to reliably differentiate between poor and good image quality (Roalf et al., 2016). For each subject and imaging session, tSNR was estimated by first calculating the mean and standard deviation of each voxel's intensity over time, and then averaging the resulting values across all brain voxels to yield a single metric of image quality (Farrell et al., 2007; Roalf et al., 2016). Note that this calculation was performed separately for the b = 0 s/mm2 and b = 1000 s/mm2 shells, respectively (Farrell et al., 2007; Roalf et al., 2016). tSNR was not calculated for the shells with higher diffusion weightings because their signal intensity and noise profiles differ considerably from volumes without diffusion weighting (Farrell et al., 2007).
TSNR data were subsequently compared between measurement sessions (MRI1, MRI2, MRI3) by means of a one-way repeated-measurements ANOVAs (RM-ANOVAs). Furthermore, we checked for the presence of tSNR outliers using the absolute deviation around the median (MAD) framework (Leys et al., 2013). Specifically, an individual's tSNR value was categorized as an outlier if it was lower than 2.5 times the MAD below the group median (Leys et al., 2013).
Furthermore, we calculated the reliability metrics within-subject coefficient of variation (CVWS, Hopkins, 2000) and intraclass coefficient (two-way mixed model, single measure, absolute agreement; cf. McGraw and Wong, 1996) based on each individual's median voxel on the GM skeleton at time points MRI1 and MRI2 (i.e., in the absence of training). To do so, we employed the procedures and formulas that we have previously used to assess NODDI reliability in WM (Lehmann et al., 2021).
Not least, we tested for the presence of systematic error or drift in the global signal by comparing the median voxel of NODDI maps on the GM skeleton between measurement sessions (MRI1, MRI2, MRI3) with a one-way RM-ANOVA. Note in this respect that a global change in microstructural maps after training is not expected. Rather, we assume that plasticity effects occur in a regional fashion in networks that were previously identified to change with stabilometer training (Taubert et al., 2010, 2016; Lehmann et al., 2020, 2022).
Assessment of behaviorally relevant plasticity
Whole-brain voxel-wise statistical analyses on skeletonized microstructural maps were conducted using the Permutation Analysis of Linear Models (PALM) v. alpha115 toolbox (Winkler et al., 2014, 2016) running in a MATLAB R2017B environment. Specifically, we aimed to identify clusters of voxels simultaneously showing both (1) changes in microstructural indices over time and (2) a correlation between neuroplasticity and motor learning (C. Thomas and Baker, 2013). To this end, the global hypothesis of behaviorally relevant neuroplasticity was broken down into a set of four subhypotheses (general linear models), each sensitive to the empirical predictions of the theory (for a graphical overview, see Fig. 1). PALM's modified nonparametric combination (NPC) framework was then used to collate the information from these partial tests to produce a single measurement summarizing the statistical evidence for the underlying complex theory (Winkler et al., 2016).
Overview of the statistical analyses using the NPC framework. For each imaging modality, change images between the MRI measurements before (MRI1, MRI2) and after the learning phase (MRI3) were calculated. Afterwards, both resulting change images (ΔMRI1–MRI3 and ΔMRI2–MRI3) were subjected to a one sample t test against zero to address learning-related structural plasticity (top row). Likewise, two submodels addressing the correlation between brain changes and DBT learning rate were set up (bottom row). Union-intersection tests (UIT) were then conducted to identify clusters of voxels in which there was, on aggregate, evidence for learning-related plasticity in conjunction with behavioral relevant plasticity. To this end, the UIT outputs a single measurement that summarizes evidence across the four submodels in every voxel.
To practically implement the abovementioned complex hypothesis test, we first calculated percentage change images between MRI1 and MRI3 as well as between MRI2 and MRI3 and created design matrices as outlined in the following. The first set of designs addressed whether motor learning did induce changes in WM or GM tissue microstructure. The corresponding hypothesis was formalized by means of one-sample t tests comparing voxelwise percentage change maps against zero. To emphasize the aspect of reproducibility, the t test was not only conducted on percentage change images between MRI2–MRI3, but also for the time interval MRI1–MRI3 (Fig. 1). The idea behind this is that the neuroplastic effect should be present regardless of whether MRI1 or MRI2 is defined as the baseline, provided that the maps were unchanged and reliable between MRI1 and MRI2 (Lehmann et al., 2021). The other design matrices postulate that learning-related neuroplasticity correlates with DBT learning. To this end, we used regression-type models to test for a linear relationship between (percentage) changes in the microstructural maps and (residualized) DBT learning rate. Again, these tests were conducted for the time intervals MRI1–MRI3 and MRI2–MRI3. We used directional (one-sided) contrasts based on the anticipated pattern of results (for example, increase in microstructure and positive correlation with behavior, or vice versa) in all submodels.
We then used the modified NPC framework as implemented in PALM (Winkler et al., 2014, 2016) to jointly test and aggregate significance across the statistical submodels introduced above. In brief, NPC tests the null hypothesis that the null hypotheses of all partial tests are true (alternative hypothesis: any is false). NPC works by first testing each component hypothesis separately using synchronized permutations. Because of this random resampling, NPC relies on minimal statistical assumptions (e.g., with respect to the distribution of the data). Second, the resulting statistics for each permutation are recorded and based on them, the empirical null distribution for each of the partial tests is constructed. Finally, at each voxel, the test statistics resulting from the previous step are combined across subtests, in our case by using Fisher's (1932) combining function. The joint statistic is considered to be significant if an aggregate of the partial tests is significant, therefore forming a union-intersection test. Note also that the NPC procedure implicitly deals with nonindependence between partial tests under the null hypothesis owing to the use of synchronized permutations.
We ran the modified NPC procedure with 5000 permutations of the data to build up the empirical null distribution from which statistical inference was performed. Cluster-like structures in the images were formed using the threshold-free cluster enhancement (TFCE) approach (Smith and Nichols, 2009). The resulting statistical maps were then subjected to cluster-based family-wise error (FWE) correction by using the distribution of the maximum statistic (Smith and Nichols, 2009; Winkler et al., 2014). Results were considered to be significant at p < 0.05 (FWE-corrected), while potentially meaningful statistical trends were assessed with a more liberal FWE-threshold of p < 0.1.
Intermodal relationships of microstructural changes
For significant results emerging from whole-brain analyses, we assessed intermodal relationships between microstructural change maps to gain a more granular insight into the plasticity process. Specifically, we were interested in processes that might drive the plasticity results. In addition to looking at the correlations between learning-related changes in the NODDI maps (ODI, NDI, ISO), here we also consider myelin-sensitive MTsat images (Lazari and Lipp, 2021) to examine in detail whether changes in the NODDI maps covary with myelin-related processes or are relatively independent of them. To this end, within significant clusters, we extracted the median voxel intensities of the microstructural change maps (ΔMRI2–MRI3) from the skeleton. Afterwards, we fitted multiple linear regression models with the modality of the significant cluster as dependent variable and the remaining modalities as independent variables.
It has been suggested that MRI-derived microstructural maps sample different, yet correlated, biophysical properties of the local tissue environment, depending on the MR sequence parameters (Weiskopf et al., 2021). Therefore, it must be assumed that the amount of shared variance (collinearity) between the predictor variables is high, which complicates inferences about the relative importance of predictors. To circumvent this problem, we have calculated relative weights (RWs) instead of “standard” regression coefficients. Briefly, from the original regressors, the RW analysis approach creates a new set of regressors that is as highly correlated as possible with the original set of variables but orthogonal (uncorrelated) to each other (Johnson, 2000). Therefore, an RW reflects the contribution of a predictor to variance of the dependent variable, considering the predictor's unique contribution as well as its common contribution with the other predictors in the model. RW analyses decomposes the coefficient of determination (multiple R2) of the model into contributions from the individual regressors (Johnson, 2000). Regression models and RWs were calculated using R's standard library (v 4.1.1) and the tool RWA Web with bootstrap confidence interval (CI) estimation (Tonidandel and LeBreton, 2015). For statistical significance testing, 10,000 samples were drawn using the bias-corrected and accelerated (BCa) method of bootstrapping.
Presence of macrostructural and tissue volume changes
To be able to interpret microstructural changes especially in GM, it is paramount to also consider colocalized changes in cortical macrostructure and tissue volume. For example, a loss of neurite structures in combination with cortical thinning might manifest itself in an increase of the NODDI parameter neurite density (Colgan et al., 2016). Furthermore, as the NODDI metrics ODI and NDI are corrected for the contribution of free water (ISO), interpreting their changes during learning should consider potential changes in tissue volume (1-ISO).
With respect to cortical macrostructure, we estimated cortical thickness based on T1-weighted images using the default settings of the CAT12 toolbox v12.7 r1738 (Christian Gaser, Structural Brain Mapping Group, Jena University Hospital; http://www.neuro.uni-jena.de/cat12/) within SPM12 v7771 (Statistical Parametric Mapping, Wellcome Trust Centre for Neuroimaging; http://www.fil.ion.ucl.ac.uk/spm/software/spm12/) for MATLAB R2017b (MathWorks). After preprocessing, the Human Connectome Project Multi-Modal Parcellation (HCP-MMP; Glasser et al., 2016) atlas was mapped to the estimated surfaces of each individual and measurement point by using the spherical registration parameters obtained from the surface-based processing (Gaser et al., 2022). Afterwards, cortical thickness was calculated for selected task-relevant regions of interest (ROIs) in native space. A particular brain region qualified as an ROI if both of the following criteria were given: (1) presence of a significant neuroplastic effect showing up in both hemispheres based on the exploratory whole-brain (NPC) analyses in the present study, (2) evidence of neuroplastic change after DBT training found in previous studies (Taubert et al., 2010, 2016; Lehmann et al., 2020). For each ROI, we calculated percentage change in cortical thickness between measurement points MRI2 and MRI3 and compared the obtained values against zero by means of one-sample t tests. Finally, the resulting p-values from the previous step were combined using Brown's (1975) method, which takes into consideration the internal correlation (dependency) of the data (Cinar and Viechtbauer, 2022). If the aggregated p-value is significant, this would lend support to the notion of regional specificity of the neuroplastic effect. Similarly, we calculated percentage change scores between MRI2 and MRI3 in NODDI-derived free water maps (ISO) based on the median voxel within (volumetric) ROIs, conducted one-sample t tests against zero and pooled the results with Brown's (1975) method; t tests and p-value combinations were done using R v4.1.1 and the R package poolr (Cinar and Viechtbauer, 2022).
Results
Motor skill practice improved performance in the DBT
We started our analyses by analyzing motor performance of the DBT. We found that the general power function yielded an adequate fit to individual's learning data (Fig. 2) with a median coefficient of determination of R2 = 0.93 (25th percentile: 0.89, 75th percentile: 0.95). Participants' learning rates (exponent of the power function) were significantly different from zero, t(23) = 13.84, p < 0.001, Hedges' g = 2.73, 95% CI [1.59,3.86], suggesting that the DBT performance of the sample improved during training (large effect).
DBT performance significantly improved over eight sessions of learning. Colored dots represent the averaged within-session performance of the n = 24 participants. To assess learning, each individual's data were fitted with a general power function (colored lines).
Quality and reliability of MR metrics
Before commencing with the analysis of training-induced brain plasticity, we compared the MRI measurement sessions for differences in tSNR (Farrell et al., 2007; Roalf et al., 2016) and checked the data for the presence of outliers. By means of One-Way RM-ANOVAs (including MRI1, MRI2, MRI3), we found that tSNR was stable over time (b = 0 s/mm2 data: F(1.43,32.78) = 0.62, p = 0.49, Greenhouse-Geisser corrected; b = 1000 s/mm2 data: F(2,46) = 0.25, p = 0.78). Likewise, individual tSNR values did not fall outside the predefined rejection criterion of 2.5*MAD below the median (Leys et al., 2013), thus indicating the absence of issues with respect to MR image quality.
Scan-rescan reliability was assessed between measurements MRI1 and MRI2, a time interval in which no change was expected because of the absence of training. With respect to microstructural maps in WM, our analyses revealed high measurement precision and remarkable properties to discriminate between individuals. Because the reliability of NODDI metrics was restricted to WM in our previous paper (cf. Lehmann et al., 2021), we calculated absolute (CVWS) and relative (ICC) reliabilities only for NODDI maps in the GM. To this end, we extracted the median voxel of GM-optimized NODDI maps (Fukutomi et al., 2018; Guerrero et al., 2019) from the GBSS-skeleton of the measurements MRI1 and MRI2. Two-way mixed model ICC (McGraw and Wong, 1996) and CVWS (Hopkins, 2000) revealed good-to-excellent reliability of the NODDI metrics in GM tissue (Fig. 3). We also subjected the NODDI maps to a one-way RM-ANOVA to test whether median voxels extracted from the GBSS-derived skeleton means would differ depending on measurement point (MRI1, MRI2, MRI3). Again, no differences in NDI, F(2,46) = 0.24, p = 0.79, ODI, F(2,46) = 1.38, p = 0.26 and FISO, F(2,46) = 0.02, p = 0.98 were found, therefore indicating the absence of systematic error (Fig. 3).
Assessment of reliability and systematic error of NODDI metrics in the gray matter. All statistics and graphs presented are based on the median of voxels on the GBSS whole-brain skeleton (50th percentile) extracted separately for each of the n = 24 participants at the different time points of the experiment (MRI1, MRI2, MRI3). A, Examplary scatter plot of the test-retest correlation of GM-ODI during the control phase of the experiment (MRI1–MRI2). B, Relative (ICC) and absolute (CVWS) test-retest-reliability during the control phase of the experiment (MRI1–MRI2). Error bars represent 95% asymmetric CIs/confidence limits around the point estimates. All ICCs were significant at p < 0.001. C, Global values of NODDI metrics over the time course of the experiment. No significant differences between measurement points were found. ns = not significant.
Behaviorally relevant changes in cortical microstructure induced by motor learning
Using the modified NPC framework (Winkler et al., 2014, 2016), we continued by performing whole-brain analyses aimed at identifying clusters of voxels whose changes were both significantly different from zero (plasticity) and correlated with the DBT learning rate (behavioral relevance).
With respect to GM microstructure, we found a pattern of results consistent with the assumption that behaviorally relevant structural changes occurred in the brain during DBT learning. Specifically, a learning-associated increase in the NODDI metric ODI and a colocalized positive correlation between changes in ODI and DBT learning rate were found (Fig. 4; Table 1). The cluster extended mainly over the primary motor cortex, the anterior cingulate cortex, and the medial prefrontal cortex, as well as several subdivisions of the dorsolateral prefrontal cortex. Because GM-ODI within this cluster was unchanged and reliable in the time interval without training (Fig. 4), it is unlikely that the neuroplastic effect was because of random measurement noise. Likewise, quality of diffusion MR images (averaged tSNR across MRI2 and MRI3 of the b = 0 s/mm2 data) was unrelated to percentage change in ODI during the learning period (MRI2–MRI3), r = −0.11, p = 0.61. Also note that within-cluster ODI changes between MRI2 and MRI3 (as well as MRI1 and MRI3) were correlated with DBT learning, whereas changes between MRI1 and MRI2 did not predict DBT learning (Fig. 4).
Peak coordinates and localization of significant clusters emerging from the voxel-based NPC analysis
Training-related, behaviorally relevant changes in the cortical orientation dispersion index (ODI) derived from the NODDI model. The figure demonstrates that cortical ODI changed noticeably as a result of DBT learning; furthermore, the degree of change in cortical ODI correlated with improved performance (i.e., faster learning). A, Results from the whole-brain union-intersection test (UIT) on Δ_ODI maps based on the NPC methodology. Significant clusters depict voxels, in which the UIT revealed evidence for a training-induced increase in ODI as well as a correlation between Δ_ODI and DBT learning rate. For visualization, the clusters derived from PALM were fattened with the “tbss_fill” script and subsequently projected from MNI152 space onto the surface of the FSaverage brain using the CAT12 toolbox (Gaser et al., 2022). Clusters are displayed at p < 0.05, FWE-corrected (TFCE). Panels B and C show descriptive data illustrating the underlying idea of the NPC analysis. B, The data underlying the significant NPC cluster are presented as connected dots (within-cluster median of ODI of the n = 24 subjects at each measurement point, extracted after back-projection to native space). Boxplots are overlaid for each time point. Note that ODI in the time interval without learning (MRI1–MRI2) was highly reliable, both in terms of within-subject variation, CVWS = 1.83%, 95% CI [1.48,2.43] and intraclass correlation, ICC(C,1) = 0.86, 95% CI [0.71,0.94], p < 0.001. Also note the increase in ODI visible at MRI3 (postlearning). C, Scatterplots with lines of best fit show the relation between within-cluster Δ_ODI (calculated based on the data shown in panel B in different time intervals) and DBT learning rate. Positive correlations are apparent for the time intervals involving the learning period (i.e., MRI1–MRI3 and MRI2–MRI3), while data were randomly scattered for the correlation in the time interval without training (MRI1–MRI2; r = −0.13, p = 0.55). The correlations between ODI changes and DBT learning rate were significantly stronger in the time intervals with training than in the time intervals without training (MRI1–MRI3 vs MRI1–MRI2: z = 4.60, p < 0.001; MRI2–MRI3 vs MRI1–MRI2: z = 3.29, p = 0.001), as determined by Steiger's (1980) z-test.
Cortical macrostructure and tissue volume fraction unaltered by motor learning
Next, we looked for the presence of changes in cortical thickness and free water (ISO), the latter as an indicator of tissue volume, during motor learning. Based on the results of the NPC whole-brain analysis in the present study and the evidence from previous studies (Taubert et al., 2010, 2016; Lehmann et al., 2020, 2022), we chose the bilateral primary motor cortex (area four from the HCP-MMP atlas; Glasser et al., 2016) and the bilateral supplementary motor areas (6mp, 6ma, and SCEF according to Glasser et al., 2016) as ROIs. Separate for each modality, we then compared percentage change during learning (MRI2–MRI3) against zero by means of one-sample t tests. The eight p-values retrieved were then pooled within-modality using Brown's (1975) method on the basis of the calculated correlation structure. No significant training-induced changes in cortical thickness, χ2(7.98) = 8.91, p = 0.35 and free water (ISO), χ2(9.23) = 8.12, p = 0.54 were found.
Changes in orientational coherence are correlated with changes in neurite density, but not with intracortical myelin
While ODI is a measure representing geometric complexity (angular variation) of neurite orientation, several other microstructural processes have been proposed to contribute to structural plasticity in GM (Blumenfeld-Katzir et al., 2011; Keifer et al., 2015; Asan et al., 2021; Schmidt et al., 2021; Matsuda et al., 2022; Mediavilla et al., 2022). To gain more granular insight into the neuroplastic process during motor learning, we sought to predict changes in ODI by a linear combination of NDI, ISO, and MTsat change (Fig. 5). The multiple linear regression model significantly explained variance in the criterion Δ_ODI, multiple R2 = 0.40, adjusted R2 = 0.30, F(3,20) = 4.36, p = 0.02. Further analyses revealed that Δ_NDI was positively correlated with Δ_ODI and accounted for 77.61% of explainable variance in the criterion, which is a statistically significant contribution, RW = 0.31, 95% BCa CI [0.03,0.52]. On the contrary, the RWs for Δ_ISO and Δ_MTsat were not significant. Therefore, Δ_NDI was the most important predictor, differing significantly from the RWs obtained for Δ_ISO, RW_diff = −0.22, 95% BCa CI [−0.50,−0.09] and for Δ_MTsat, RW_diff = −0.31, 95% BCa CI [−0.59,−0.06].
Microstructural correlates of Δ_ODI in the cluster from NPC whole-brain analysis. A, Results of a multiple regression model with changes in Δ_NDI, Δ_ISO, and Δ_MT as predictors of Δ_ODI. Relative weights analysis was conducted to partition multiple R2 among the predictor variables, i.e., RWs of all predictors (depicted as points) sum to the multiple R2 of the model. Error bars represent 95% BCa CIs used to test the statistical significance of the RWs based on 10,000 bootstrap samples. B, Scatterplot with line of best fit shows the positive correlation between Δ_NDI and Δ_ODI.
Results of NODDI and DTI metrics in white matter
Whole-brain NPC analyses did not show nominally significant results in the WM as assessed with TBSS. However, we observed statistical trends for radial (λ⊥; pFWE < 0.06) and mean diffusivities (pFWE < 0.09) in prefrontal fiber tracts. NODDI results showed no trend regarding possible behaviorally relevant neuroplasticity (all ps > 0.1).
Discussion
Establishing behaviorally relevant structural changes in the human brain is a challenging endeavor. In addition to calls for rigorous study designs and appropriate analytical approaches to address this topic (C. Thomas and Baker, 2013), most studies to date have used MRI pulse sequences that do not allow for a biologically specific interpretation of neuroplastic changes (Tardif et al., 2016; Lerch et al., 2017). Using a powerful within-subject design in combination with multishell diffusion MRI, we show that training a complex DBT over four weeks increases the organizational complexity of cortical microstructure. Importantly, the microstructural changes during training correlated with the change in DBT performance, and structural remodeling of the brain was specific, meaning that it selectively affected brain areas associated with motor processes, in the absence of systematic global changes of tissue microstructure over time.
To the best of our knowledge, only two previous studies have investigated neuroplasticity in GM using multishell diffusion MRI combined with advanced biophysical modeling accounting for non-Gaussian diffusion. In response to two short learning sessions of a car-racing game, Tavor et al. (2013) report an increase in restricted diffusion in prefrontal, temporal, and parietal areas based on the CHARMED model (Assaf and Basser, 2005). Furthermore, Vukovic et al. (2021) detected widespread diffusion kurtosis (Jensen et al., 2005) changes following a brief language learning session. However, to our knowledge, the present study is the first to find cortical microstructural plasticity after a lengthy period of motor learning, as well as the first to report a correlation between microstructural change in GM and concomitant behavioral change.
Learning-induced changes in GM microstructure were predominantly present in the primary sensorimotor cortex, the supplementary and cingulate motor areas, the premotor cortex, but also extended to prefrontal areas. These regions are generally known to be involved in movement planning, motor control, and learning (Strick et al., 2021), but they have also been identified in previous work as important nodes in the structural brain network underlying DBT learning (Taubert et al., 2010, 2016; Lehmann et al., 2020, 2022).
Based on findings in animal models, it has been hypothesized that structural changes in neurites might also underlie learning in humans (Zatorre et al., 2012). Previous studies have indeed shown that diffusion MRI in general and the orientation distribution function derived from it in particular are sensitive to neuritic elements and their organizational complexity (Colgan et al., 2016; Assaf, 2019; Wang et al., 2019; Mak et al., 2021; Reveley et al., 2022). In light of such findings, a previous study in humans interpreted a long-lasting learning-induced reduction in neocortical mean diffusivity as an indication that a relatively stable memory engram had formed (Brodt et al., 2018). Although interpretation of NODDI indices in terms of the underlying cellular correlates remains challenging (Kamiya et al., 2020), our finding of a behaviorally relevant increase in the variability of neurite orientations in task-related brain areas seems also consistent with the notion that a motor engram has been formed. In this sense, on the one hand, it has been shown in animal models that motor learning is closely linked to structural changes in dendrites (Xu et al., 2009; Hayashi-Takagi et al., 2015; Hwang et al., 2022), and on the other hand, NODDI's ODI has commonly been interpreted as an index of dendritic complexity (Colgan et al., 2016; Assaf, 2019; Wang et al., 2019; Nazeri et al., 2020; Mak et al., 2021). Further supporting this interpretation is the fact that several studies of neuroplasticity in animals combining histology and MRI found that dendritic remodeling was accompanied by MRI-detectable signal changes (Blumenfeld-Katzir et al., 2011; Keifer et al., 2015; Schmidt et al., 2021; Matsuda et al., 2022). Another candidate mechanism that could give rise to a reduced orientational coherence of water diffusion is axonal growth (Lerch et al., 2011), which, however, has been observed less frequently compared with dendritic changes.
The interpretation that neuritic processes were affected by training is additionally supported by the positive correlation between changes in ODI and colocalized changes in NDI, an index of the relative volume fraction of the intraneurite compartment in the tissue (Zhang et al., 2012). However, in contrast to previous findings in humans (Tavor et al., 2013; Vukovic et al., 2021), no evidence for a significant increase in restricted diffusion was detected. We speculate that the timescales at which neuroplastic processes operate might have played a role here. For example, mice learning a forelimb reaching task show an increase in dendritic spine density within an hour, which returns to pretraining baseline within two weeks (Xu et al., 2009). Similar rapid plasticity processes might have taken place in the studies by Tavor et al. (2013) and Vukovic et al. (2021), who found changes in GM diffusion within just a few hours or days of training. In contrast, in our study a possible increase in neurite volume early in the learning process could have already disappeared after four weeks of training (Wenger et al., 2017).
In humans, morphometric changes of the brain's GM because of training are thought to result from changes in cortical thickness, surface area, tissue density, cellular composition of tissue microstructure or some permutation of these factors (Wenger et al., 2017; Asan et al., 2021). In the present study, there was no evidence that cortical thickness was changed after training or that cortical tissue was more densely packed. This is consistent with the findings of several studies on aging and neurodegenerative diseases, which collectively suggest that there is limited covariation between tissue properties derived from T1w imaging (GM volume/cortical thickness) and those from NODDI (Nazeri et al., 2015; Mak et al., 2021; Bai et al., 2022), underlining the position that diffusion imaging-derived indices may be sensitive for detecting aspects of neuroplasticity in GM that go beyond “conventional” methods (Tavor et al., 2013; Assaf, 2019; Kamiya et al., 2020; Vukovic et al., 2021).
Besides learning-induced changes in neuritic elements and neuroglia, research conducted in animals suggests that increased myelin density may be another cellular mechanism that significantly contributes to plasticity in the neocortex (Mediavilla et al., 2022). Although it has been suggested that apical dendrites are the main driver of restricted diffusion in the neocortex (Reveley et al., 2022), NODDI metrics have shown to be associated with myelin staining in postmortem samples of rodent brains (Sepehrband et al., 2015; Wang et al., 2019) and human spinal cord (Grussu et al., 2017). Likewise, a colocalization of cortical neurite density as assessed with NODDI and myelin content as estimated based on the T1w/T2w ratio has been reported (Fukutomi et al., 2018). In the present study, however, there was no significant association between training-induced changes in myelin-sensitive magnetization transfer maps and concurrent changes in ODI, suggesting that changes in intracortical myelination did not contribute to alterations in neurite complexity.
It is plausible that some limitations might have influenced the results obtained. First, in the controlled within-subjects design we used, as opposed to a parallel randomized controlled trial, it cannot be ruled out that confounders influenced the control and learning phases of the experiment differently. However, explanations for our results other than training appear unlikely because (1) direct interaction with the balancing task, not just life-as-usual, was required to induce microstructural plasticity, (2) NODDI metrics were reliable with no signs of undue influence caused by systematic and unsystematic error, (3) NODDI changes during the training but not the control period correlated with DBT learning. Second, it has been argued that some of the modeling assumptions that NODDI makes may be oversimplified (Jones et al., 2013; Kamiya et al., 2020), especially in complex tissues such as GM. One point concerns the fact that NODDI does not explicitly model glial cells (Kamiya et al., 2020), although this would be desirable because of the fact that neuroglial changes are closely coupled to neuronal changes (Perez-Alvarez et al., 2014). Such justified criticisms must be weighed against the facts that NODDI is highly reproducible (Lehmann et al., 2021; Coelho et al., 2022) and extensively validated by histologic examinations (cf. Kraguljac et al., 2023). Third, we unexpectedly did not find statistically significant evidence for neuroplastic changes in WM. Instead, we observed a trend toward decreased diffusivity in prefrontal fiber tracts, consistent with previous findings (Taubert et al., 2010; Lehmann et al., 2020). This might suggest that our study has been underpowered to detect medium-size effects in exploratory whole-brain analysis of the WM.
In summary, this study shows that advanced diffusion MRI combined with biophysical modeling is a promising approach to gain highly specific insights into the underlying mechanisms of neuroplasticity in humans. From an applied point of view, it is interesting to note that neurite orientation dispersion, which we found to increase after motor training, is often lowered in aging and neurodegenerative conditions (Nazeri et al., 2015; Nazeri et al., 2020; Mak et al., 2021; Bai et al., 2022). In perspective, it is therefore hoped that NODDI will increasingly be used as an MRI index for evaluating treatments to maintain brain health or to ameliorate the effects of disease.
Footnotes
This work was supported by the Deutsche Forschungsgemeinschaft (DFG) Grant CRC 1436 TPC01. N.A. was supported by the innovation fund of the Otto von Guericke University Magdeburg (2039236003).
The authors declare no competing financial interests.
- Correspondence should be addressed to Nico Lehmann at nico1.lehmann{at}ovgu.de
