Neural plasticity in amplitude of low frequency fluctuation, cortical hub construction, regional homogeneity resulting from working memory training

Working memory training (WMT) induces changes in cognitive function and various neurological systems. Here, we investigated changes in recently developed resting state functional magnetic resonance imaging measures of global information processing [degree of the cortical hub, which may have a central role in information integration in the brain, degree centrality (DC)], the magnitude of intrinsic brain activity [fractional amplitude of low frequency fluctuation (fALFF)], and local connectivity (regional homogeneity) in young adults, who either underwent WMT or received no intervention for 4 weeks. Compared with no intervention, WMT increased DC in the anatomical cluster, including anterior cingulate cortex (ACC), to the medial prefrontal cortex (mPFC). Furthermore, WMT increased fALFF in the anatomical cluster including the right dorsolateral prefrontal cortex (DLPFC), frontopolar area and mPFC. WMT increased regional homogeneity in the anatomical cluster that spread from the precuneus to posterior cingulate cortex and posterior parietal cortex. These results suggest WMT-induced plasticity in spontaneous brain activity and global and local information processing in areas of the major networks of the brain during rest.


Procedure
As described in our previous study 1  Hence, subjects used the program on their personal computers, and were encouraged to perform the WMT tasks every day, in addition to the two weekly training sessions in the laboratory. Details are described in our previous study 1 . Subjects in the no-intervention group participated in pre-and post-experiments based on (a) the recommendation of the previous review 2 , (b) the lack of effects in the active control group in our previous study 3 , in which sample characteristics and experimental settings were similar to the present study, and (c) the overwhelming number of studies that, in contrast with popular belief, failed to show effects in active control groups [4][5][6] . Moreover, active control groups are inappropriate controls due to specific effects of active control training on measures where they aren't supposed to have specific effects 2 . The lack of an active control or no control groups combined with within-training group analyses with training variables are gold standard procedures for imaging studies of cognitive training 2 . In particular, no previous studies show that these well-controlled active interventions affect resting state fMRI measures, and to the best of our knowledge, there is no persuasive scientific theory that suggests changes to resting state fMRI measurements. Hence, inclusion of both active control and no intervention groups is ideal 2 , but is usually suboptimal from a statistical perspective 2 .

Training tasks
As described in our previous study 1 , four WMT tasks were presented during each training session. In all training tasks, difficulties (number of items to be remembered) were modulated based on subjects' performance. The four WMT tasks were as follows: (a) a visuospatial WM task, (b) an auditory backward operation span task, (c) a dual WM task, (d) a dual N-back task.
(a) In the visuospatial WM task, circles are presented one at a time (1/s rate) in an interface of 10 irregularly distributed squares (circles are presented in one of these squares). After the presentation of stimuli, subjects indicate the location and order of the presented stimuli by clicking on a computer screen with a mouse. (b) In an auditory backward operation span task, pairs of single digits (0-9) are verbally presented at a rate of 1 pair/3 s. Within this 3-s period, one digit is presented per second, but in the final second, no stimuli are presented. Therefore, four pairs of digits are presented as follows: 1, 3, no stimulus; 4, 9, no stimulus; 3, 7, no stimulus; and 2, 5, no stimulus. Subjects were required to remember the value that was the sum of the presented pairs of digits, and the order in which they were presented (in the above example, they were required to remember 4, 3, 0, 7). After presentation, subjects were required to repeat the sequence by pressing numbered buttons on the screen in reverse order (7, 0, 3, 4 in the above example). (c) In the dual WM task, which is similar to that of a previous study 7 , subjects concurrently performed a visuospatial WM task and an auditory digit span task.
In that study, circles were presented one at a time at a rate of 1/3 s in the same interface used for the task (a). After the presentation of stimuli, subjects indicated the location and order of the presented stimuli by clicking on a computer screen with a mouse. One digit (0-9) was verbally and simultaneously presented with the circle presentation. After the presentation of stimuli, subjects indicated the digits and the order of stimuli by pressing numbered buttons on the screen in the order of presentation. Subjects could perform either task first. (d) In the dual WM task, which was similar to that in a previous study of WMT 8 , squares at eight different locations were sequentially presented on a computer screen at a rate of 3/1 s (stimulus length, 500 ms; inter-stimulus interval, 2,500 ms). Simultaneously, one of eight consonants was sequentially presented using headphones and a response was required whenever one of the presented stimuli matched the stimuli previously presented at the (n) position in the sequence. Additional details regarding the practical aspects of the task procedures and training are described in our previous study 1 . Training times for each task were not necessarily exactly equal among the four tasks, although tasks were performed in the fixed order [(a)-(d)]. The same task was not used for assessment prior to the intervention period.

Preprocessing and individual-level analysis of imaging data
Preprocessing of imaging data was performed using SPM8 implemented in Matlab and SPM8's extension software DPARSF (Data Processing Assistant for Resting-state fMRI). In the following procedures, we avoided co-registration and co-normalization procedures (which use the same normalization parameters for both pre-and post-images) for pre-and post-images, including registration of mean images to preand post-images, because of concerns of possible bias or problems occurring when image properties were substantially different between pre-and post-images 9 .
For each session and each subject, the first image of a series of BOLD images was skull stripped by masking the images using the threshold of a given signal intensity from spatially smoothed (using 8 mm FWHM) BOLD images. This skin-skull-stripping procedure was performed so that these parts were not treated as the outer edge of the brain parenchyma in the preprocessing procedures. Furthermore, the skull-stripped BOLD image was co-registered to a previously created custom made skull-stripped EPI template. The series of BOLD images for each session for each subject was slice-timing corrected and realigned using DPARSF. The series of BOLD images for each session for each subject were segmented and independently normalized on the basis of the modified diffeomorphic anatomical registration through exponentiated lie algebra (DARTEL)-based methods 10 to give images with 3.75 × 3.75 × 3.75 mm 3 voxels. In this process, a custom template was also created 10 , and the whole brain mask, which comprises voxels that show gray matter tissue probability + white matter tissue probability + cerebrospinal fluid (CSF) tissue probability >0.1, the mask of the areas that are likely to be white matter (white matter tissue probability >0.99, to avoid contamination of signals from white matter), and the mask of the areas that are strongly likely to be CSF (CSF tissue probability >0.99, to avoid the contamination of signals from other tissues) were created from the custom template.
The normalized series of BOLD images were processed by DPARSF for individual level analysis. Initially, 26 nuisance covariates (including mean signals from the voxels within the white matter mask, mean signals from the voxels within the CSF mask, and Friston 24 motion parameters were analyzed. The Friston 24-parameter model (six head motion parameters, six head motion parameters one time point before, and the 12 corresponding squared items) 11 was used to regress out head motion effects.
Recent work indicates that regressing out Friston 24-parameters is more effective than other movement correction methods, such as corrections for rigid-body using six parameters, derivative 12 parameters, and 12 voxel-specific regressors 12 . To further eliminate residual effects of motion on rsfMRI measurements, volume-level mean framewise displacements were computed and used as covariates in second-level analyses 13 .
Recent controversy regarding whether or not the whole brain signal should be regressed out in the rsfMRI analysis continues 14 , and regressing out the whole brain signal may have merits 12 and weaknesses 15 . In this study, we did not regress out the whole brain signal because it partly reflects global brain activity 16 , and it has been suggested that WMT may affect global brain activity 1 . Furthermore, although regressing out the whole brain signal facilitates the distinction between network specific brain activity and global brain activity, our main focuses fALFF and DC do not involve network specific activity and regressing out whole brain activity is particularly conceptually problematic for calculating DC, which is the correlation between one voxel and the whole brain.
Processed images were spatially smoothed with 8-mm FWHM, and the resulting images were masked with the whole brain mask that was created as described above.
Analyses of fALFF were performed using the DPARSF software as previously described 17,18 . The time series of each voxel was transformed into the frequency domain, and the power spectrum was obtained. Because the power of a given frequency is proportional to the square of the amplitude of that frequency component, the square root was calculated at each frequency of the power spectrum, and the average square root was then obtained across 0.01-0.08 Hz at each voxel. This average square root was taken as ALFF. For analyses of fALFF, the ratio of the power of each frequency at a low frequency range to that of the entire frequency range (fALFF) was computed as described previously 19,20 . Specifically, after preprocessing, time series for each voxel were transformed into the frequency domain without band-pass filtering. The square root was calculated at each frequency of the power spectrum, and the sum of the amplitudes across 0.01-0.08 Hz was divided by that of the entire frequency range.
After preprocessing, fMRI data were temporally band-pass filtered (0.01 < f < 0.08 Hz) to reduce low frequency drift and high frequency. Furthermore, weighted DC measures were calculated using DPARSF, as previously described 21 . Briefly, Pearson correlation coefficients were initially computed between the time series of all pairs of gray matter voxels, leading to whole-brain functional connectivity matrices for each individual. We calculated this across the whole brain regardless of area because previous studies using large sample sizes have shown the strength of functional connectivity with the periphery of functional networks and gray matter (rather than deep areas within functional areas or gray matter) are important for cognition 10 . Individual correlation matrices were transformed into a Z-score matrix using Fisher's r-z transformation to improve normality. We further computed weighted DC strengths of voxels as the sum of the connections (Z-values) with all other voxels. As described previously 22 , we conservatively restricted our analysis to positive correlations above a threshold of r = 0.25. A relatively higher threshold was chosen to avoid counting voxels with weak correlations that reflect signal noise. Finally, standardized weighted DC maps were acquired by subtracting mean values from within the abovementioned gray matter mask, and then dividing by the standard deviation of the whole gray matter mask 23 .
Voxels with higher DC values reflect the sum of the connections for a given voxel and all other voxels, and are thus indicative of their central roles in transferring information across brain regions.
Subsequently, we calculated regional homogeneity images according to previously described procedures 24 . Initially, spatially normalized rsfMRI images were band-pass filtered (0.01 < f < 0.08 Hz) and masked using the whole brain mask. A regional homogeneity value was then calculated to measure similarity of the time series of a given voxel with its nearest 26 voxels. After normalization, the resulting images were spatially smoothed using 8-mm FWHM. Normalization of regional homogeneity maps was performed by dividing the regional homogeneity among each voxel by the average regional homogeneity of the whole brain. This regional homogeneity value is known as Kendall's coefficient concordance (KCC) 25 , and is used to measure the similarity of a time series of a given voxel with that of its nearest voxels in a voxel-wise way, based on the assumption that voxels are temporally similar to their neighbors 26 .
Using this method, we measured temporal homogeneity of regional BOLD signals, which reflected the temporal homogeneity of neural activity. The formula to calculate KCC has been expounded in previous studies 27 . Data obtained from 1 subject whose final performance data were missing were removed from the calculation of the average in this task.