DPARSF: A MATLAB Toolbox for “Pipeline” Data Analysis of Resting-State fMRI

Resting-state functional magnetic resonance imaging (fMRI) has attracted more and more attention because of its effectiveness, simplicity and non-invasiveness in exploration of the intrinsic functional architecture of the human brain. However, user-friendly toolbox for “pipeline” data analysis of resting-state fMRI is still lacking. Based on some functions in Statistical Parametric Mapping (SPM) and Resting-State fMRI Data Analysis Toolkit (REST), we have developed a MATLAB toolbox called Data Processing Assistant for Resting-State fMRI (DPARSF) for “pipeline” data analysis of resting-state fMRI. After the user arranges the Digital Imaging and Communications in Medicine (DICOM) files and click a few buttons to set parameters, DPARSF will then give all the preprocessed (slice timing, realign, normalize, smooth) data and results for functional connectivity, regional homogeneity, amplitude of low-frequency fluctuation (ALFF), and fractional ALFF. DPARSF can also create a report for excluding subjects with excessive head motion and generate a set of pictures for easily checking the effect of normalization. In addition, users can also use DPARSF to extract time courses from regions of interest.


INTRODUCTION
Resting-state functional magnetic resonance imaging (fMRI) has been more and more widely used since Biswal et al. (1995) fi rstly reported the presence of spatially coherent activity in the resting-state blood oxygen level-dependent (BOLD) fMRI signal. Resting-state fMRI is considered as a powerful tool for investigating the spontaneous neuronal activity which consumes most of the brain's energy (Fox and Raichle, 2007). In addition, resting-state fMRI is also a convenient way for clinical studies since it has advantages of reasonable spatial and temporal resolution and non-invasiveness, as well as its simplicity that does not need to set complicated cognitive tasks.
Functional connectivity (FC) is widely used in resting-state fMRI studies (Biswal et al., 1995;Lowe et al., 1998;Xiong et al., 1999;Cordes et al., 2000;Greicius et al., 2003;Fox et al., 2005Fox et al., , 2006Fransson, 2005;Vincent et al., 2006). While FC measures the signal synchrony among remote brain areas, the regional spontaneous activity could be examined by several metrics, such as the regional homogeneity (ReHo, Zang et al., 2004), the amplitude of low-frequency fl uctuation (ALFF, Zang et al., 2007) and the fractional ALFF (fALFF, Zou et al., 2008). All the aforementioned methods could be calculated by a toolbox Resting-State fMRI Data Analysis Toolkit 1 (REST). As an easy-to-use MATLAB toolbox, REST is compatible with Statistical Parametric Mapping 2 (SPM). The data could be preprocessed by SPM and then entered into REST's analysis. Although SPM is a powerful tool, lots of complicated and timeconsuming operations are needed when analyzing large sample data set. In SPM, the parameters need to be set step-by-step and subject-by-subject. These manual procedures may be time-con-DPARSF: a MATLAB toolbox for "pipeline" data analysis of resting-state fMRI important information of left hemisphere and right hemisphere. If this converting option is checked, DPARSF will convert the data by calling dcm2nii in MRIcroN software 3 . If users have converted the data previously, this option is not needed.

REMOVE FIRST 10 (MORE OR LESS) TIME POINTS
The fi rst few volumes of the functional images are often discarded for signal equilibrium and to allow the participants' adaptation to the scanning noise. If this option is checked, DPARSF will delete the specifi ed number of time points for each participant.

SLICE TIMING
Most fMRI data are acquired using two-dimensional pulse sequences that acquire images one slice at a time, thus all slices are acquired at different time within a repeat time (TR). Timing differences are especially problematic for longer TR. Hence the differences in image acquisition time between slices need to be corrected. The number of slices, slice order and reference slice need to be specifi ed, then DPARSF will do slice timing by calling functions in SPM.

HEAD MOTION CORRECTION
The goal of motion correction is to adjust the time series of images so that the brain is in the same position in every image (Huettel et al., 2004). If this option is checked, the time-series of images will be motion-corrected by calling functions in SPM. Since excessive head motion may induce large artifact in fMRI time-series, participants with excessive head motion need to be excluded from further analysis. DPARSF will create a report of head motion based on the realign parameters estimated by SPM (as shown in "ExcludeSubjects.txt" in the "RealignParameter" directory).

NORMALIZATION
The brain size, shape, orientation, and gyral anatomy vary largely across participants. For inter-subject comparison to be feasible, the individual brain is usually transformed or spatially normalized into a standardized template. SPM provides two optional ways to normalize the functional images into the Montreal Neurological Institute (MNI) space: (1) using echo-planar imaging (EPI) template (Ashburner and Friston, 1999) and (2) using unifi ed segmentation on T1 image (Ashburner and Friston, 2005). The latter way could improve the accuracy of spatial normalization (Ashburner and Friston, 2005) but is a little complicated in SPM. It contains three steps -coregistration, segmentation and writing normalization parameters. DPARSF has integrated these three steps into one. It's important to check the effect of normalization of each individual since some data may meet problem in normalization. DPARSF can generate a set of pictures (Figure 2) for easily checking the effect of normalization. It should be noted that DPARSF provides a very simple way for visual inspection. Users should check the effect of spatial normalization carefully.

SMOOTHING
Smoothing is used as a preprocessing step to suppress noise and effects due to residual differences in functional and gyral anatomy during inter-subject averaging. The most common smoothing technique is the Gaussian fi lter which has the shape of a normal distribution. DPARSF will smooth the data with the specifi ed width at half of the maximum value (full-width-half-maximum, or FWHM) by calling functions in SPM.

REMOVE LINEAR TREAD
Long-term physiological shifts, movement related noise remaining after realignment or instrumental instability may contribute to a systematic increase or decrease in the signal with time (Turner et al., 1997;Lowe and Russell, 1999). The exact cause for the drift of the baseline signal is not completely understood (Smith et al., 1999), and how this structured trend affect further analysis is an interesting issue. If this option is checked, DPARSF will remove the systematic drift or trend using linear model as does in Analysis of Functional Neuroimage (AFNI) (Cox, 1996) by calling functions in REST.

FILTERING
Low frequency (0.01-0.08 Hz) fl uctuations (LFFs) of the restingstate fMRI signal were reported to be of physiological importance (Biswal et al., 1995) and also were suggested to refl ect spontaneous neuronal activity (Lu et al., 2007). Zuo et al. (2010) also reported low frequency oscillations (0.01-0.073 Hz) were primarily detected within gray matter. In contrast, relatively high frequency oscillations (0.073-0.25) were primarily restricted to white matter. It's reported that respiratory and aliased cardiac signals fall in the range of relatively high frequency band (Cordes et al., 2001). Thus, the data is usually bandpass (e.g., 0.01-0.08 Hz) fi ltered to reduce the effect of very low frequency and high frequency physiological noise. Of note, the data should not to be fi ltered when calculating fALFF, because fALFF is a ratio of low frequency amplitude to full band amplitude. If this option is checked, DPARSF will fi lter the data by calling ideal fi lter functions in REST (which is in accordance with AFNI).

REGIONAL HOMOGENEITY
While FC analysis measures the signal synchrony of LFF activity among different brain areas, it does not provide information of regional spontaneous activity. The ReHo method (Zang et al., 2004), unlike the connectivity-based methods typically used in most resting-state fMRI studies, is suitable for exploring regional brain activity by examining the degree of regional synchronization of fMRI time courses. This is accomplished on a voxel-by-voxel basis by calculating Kendall's coeffi cient of concordance (KCC, Kendall and Gibbons, 1990) of time series of a given voxel with those of its nearest neighbors. A larger value of ReHo indicates a higher regional synchronization. In order to reduce the global effects of variability across participants, as did in PET studies (Raichle et al., 2001), the ReHo of each voxel was divided by the global mean ReHo value within the whole-brain mask (default brain mask provided in REST which was thresholded at 50% on the SPM5's a priori brain mask). It needs to be noted that, smoothing before ReHo calculation will largely increase the regional similarity. We recommend that the smoothing procedure is performed after ReHo calculation. However, it is still an open issue. If this option is checked, DPARSF will calculate ReHo and then smooth the ReHo results by calling functions in REST and SPM.

ALFF AND fALFF
The regional spontaneous activities can be examined by the ALFF. Biswal et al. (1995) found that the root mean square of the LFF in the white matter was reduced by about 60% relative to the gray matter. The power spectrum of the LFF (equivalent to the square of the ALFF) has been used to indicate the magnitude of neural activity (Kiviniemi et al., 2000;He et al., 2007;Hoptman et al., 2010;Zhang et al., in press). However, it has been shown that ALFF is signifi cantly higher than the global mean ALFF in cisterns and vicinity of large blood vessels . That was apparently induced by the large fl uctuations of high frequency physiological noise. Thus an improved measure fALFF  is defi ned as the ratio of total amplitude within the low-frequency range (0.01-0.08 Hz) to the total amplitude of the entire detectable frequency range. It was found that fALFF can better reveal the default mode network (DMN) within groups (i.e., by one-sample t-test). However, which of the two measures (ALFF vs. fALFF) is better for between-groups studies is still unknown. Of note, ALFF measures have higher test-retest reliability in gray matter regions than fALFF, while more susceptible to possible artifactual fi ndings in the vicinity of blood vessels and the cerebral ventricles . Similar to standardization procedure of ReHo analysis, the ALFF or fALFF of each voxel was divided by the global mean ALFF or fALFF value within the whole-brain mask.
In addition, different frequency bands are considered to be generated by distinct oscillators, each with specifi c properties and physiological functions, as the neuronal oscillation classes are arrayed linearly when plotted on the natural logarithmic scale (Penttonen and Buzsáki, 2003;Buzsáki and Draguhn, 2004). Thus, ALFF or fALFF of different frequency bands could also be investigated. For example, ALFF and fALFF of four frequency bands, namely slow-5 (0.01-0.027 Hz), slow-4 (0.027-0.073 Hz), slow-3 (0.073-0.198 Hz), and slow-2 (0.198-0.25 Hz), were examined, and the results showed that fALFF in the slow-4 (0.027-0.073 Hz) band is relatively specifi c to the basal ganglia . If these options are checked, DPARSF will calculate ALFF and/or fALFF of the specifi ed frequency band by calling functions in REST.

REMOVE EFFECT OF NUISANCE COVARIATES
In the past few years, there has been increased attention to the anti-correlation phenomenon of resting-state fMRI. A typical case is that, while the global (whole-brain) signal was removed, many researchers consistently observed that there were signifi cant anti-correlations between the components of the default-mode and attention networks (Greicius et al., 2003;Fox et al., 2005;Fransson, 2005). Recently, the global signal has been found to be associated with respiration-induced fMRI signal (Birn et al., 2006). To reduce the effect of the physiological artifacts, the whole-brain signal would be removed by a regression analysis before FC analysis (Greicius et al., 2003;Fox et al., 2005;Fransson, 2005). Of note, it FIGURE 2 | Pictures for checking normalization. The normalized functional image was overlaid on a high resolution 3D anatomical image (the opaque one with skull. From "Colin Holmes, " http://imaging.mrc-cbu.cam.ac.uk/ downloads/Colin/, also distributed with MRIcroN as ch2) in the MNI space. Users can easily check the accuracy of spatial normalization by visual inspection.
is still an ongoing controversy since removal of the global brain signal causes the re-distribution of correlation coeffi cients and the interpretation of biological mechanisms of negative correlations is ambiguous (Murphy et al., 2009). In addition to the global mean signal, six motion parameters, the cerebrospinal fl uid (CSF), and the white matter signals would also be removed as nuisance variables to reduce the effects of head motion and non-neuronal BOLD fl uctuations (Fox et al., 2005;Kelly et al., 2008). It is still an open issue that where the ROIs should be located to represent the white matter and CSF. REST provided a few default masks made from SPM5's a priori masks, i.e., the whole brain mask (brainmask.nii) thresholded at 50%, the white matter mask (white.nii) thresholded at 90%, and the CSF mask (csf.nii) thresholded at 70%. It should be noted that the removal of nuisance covariates is not for ReHo and ALFF analysis. It is not clear yet how the nuisance covariates affect the ReHo or ALFF results.

FUNCTIONAL CONNECTIVITY
FC is widely used in resting-state fMRI (Biswal et al., 1995;Lowe et al., 1998;Xiong et al., 1999;Cordes et al., 2000;Greicius et al., 2003;Fox et al., 2005Fox et al., , 2006Fransson, 2005;Vincent et al., 2006). The correlations in spontaneous BOLD fl uctuations may refl ect the inter-regional correlations in neuronal variability (Friston et al., 1993;Horwitz, 2003). If this option is checked, the averaged time course will be obtained from a specifi ed seed region and the correlation analysis will be performed in a voxel-wise way to generate the FC map. The correlation coeffi cient map will be converted into z map by Fisher's r-to-z transform to improve the normality (Rosner, 2006) by calling functions in REST.

ILLUSTRATIONS
To validate and illustrate the usage of DPARSF, we performed the ReHo, ALFF, fALFF and FC analyses.

DATA
Data were selected from a large sample resting-state fMRI dataset of our group, which has been publicly released in the "1000 Functional Connectomes" Project 4 . We selected 86 young healthy volunteers (48 females: 20.8 ± 1.6 years old, range 18-25; and 38 males: 20.7 ± 1.7 years old, range 17-25) with head motion less than 2.0 mm displacement in any of the x, y, or z directions or 2.0° of any angular motion throughout the resting-state scan. All are righthanded and had no history of neurological and psychiatric disorders. Written informed consent was obtained from each participant, and the study was approved by the Institutional Review Board of Beijing Normal University Imaging Center for Brain Research.
MRI data were acquired using a SIEMENS TRIO 3-Tesla scanner in the Beijing Normal University Imaging Center for Brain Research. The participants lay supine with the head snugly fi xed by straps and foam pads to minimize head movement. During the resting-state session, the participants were instructed to keep as motionless as possible and not to think systematically. The functional images were obtained using an EPI sequence with the following parameters: 33 axial slices, thickness/gap = 3/0.6 mm, in-plane resolution = 64 × 64, TR = 2000 ms, TE = 30 ms, fl ip angle = 90°, FOV = 200 × 200 mm. In addition, a T1-weighted sagittal three-dimensional magnetization-prepared rapid gradient echo (MPRAGE) sequence was acquired, covering the entire brain: 128 slices, TR = 2530 ms, TE = 3.39 ms, slice thickness = 1.33 mm, fl ip angle = 7°, inversion time = 1100 ms, FOV = 256 mm × 256 mm, and in-plane resolution = 256 × 192.

PREPROCESSING
Data were processed by using DPARSF pipeline analysis as introduced in the last section. Briefl y, after converting DICOM fi les to NIFTI images, the fi rst 10 time points were discarded. Then slice timing and head motion correction were performed. The data were then normalized to MNI space by using unifi ed segmentation of T1 image and re-sampled to 3-mm isotropic voxels. After smoothing with a 4 mm FWHM Gaussian kernel (for ALFF, fALFF, FC except for ReHo), the linear trend of time courses were removed and then temporally band-pass fi ltering (0.01-0.08 Hz) (with an exception of fALFF) was performed.

ReHo, ALFF, fALFF AND FC
As indicated in the corresponding sections, spatial smoothing was performed after ReHo calculation, but for the other three (ALFF, fALFF and FC) methods, spatial smoothing was performed before their calculation. fALFF was calculated based on preprocessed data without fi ltering since fALFF is the value divided by the total power in the entire detectable frequency range. The ReHo, ALFF or fALFF of each voxel was divided by the global mean value within the whole-brain mask.
Before FC calculation, nine nuisance covariates including six head motion parameters, the global signal, the white matter signal and the CSF signal were removed from the preprocessed data. A sphere (radius = 6 mm) in the posterior cingulate cortex (PCC) (−5, −49, 40) were defi ned as the seed region for each participant in line with a previous study (Fox et al., 2005). The averaged time course was then obtained from the sphere ROI and the correlation analysis was performed in a voxel-wise way to generate the FC of the PCC, called the PCC-FC map. Finally, the correlation coefficient map was converted into z maps by Fisher's r-to-z transform to improve the normality (Rosner, 2006).

STATISTICAL ANALYSIS
To explore the within-group patterns, one-sample t-tests were performed on the ReHo, ALFF, fALFF, and FC maps, respectively, in a voxel-wise way by using SPM. For ReHo, ALFF and fALFF maps, the one-sample t-tests were to fi nd regions showing signifi cantly higher ReHo, ALFF and fALFF, respectively, than the global mean value. After being divided by the global mean value, each individual ReHo, ALFF and fALFF map has a "new" global mean of "1". Thus the onesample t-tests were performed against "1" other than "0". Since the module for one-sample t test in SPM just can compare values with base "0", we subtracted "1" from the ReHo, ALFF and fALFF maps which had been divided by global mean value of the whole brain, and then perform one-tailed one-sample t-tests on the subtracted maps in SPM. For FC maps, two-tailed one sample t-test was performed on the z maps to show both the DMN and its anti-correlated network patterns. The within-condition statistical threshold was set at t > 3.89 (P < 0.0001) for one-tailed t-tests (for ReHo, ALFF and fALFF) and |t| > 4.08 (P < 0.0001) for two-tailed t-test (for FC) and cluster size >135 mm 3 , which corresponds to a corrected P < 0.0001. This correction was confi ned within the whole-brain mask (size: 1912437 mm 3 ) and was determined by Monte Carlo simulations (Ledberg et al., 1998) that were performed by the program REST AlphaSim (which is based on AlphaSim in AFNI 5 ).

RESULTS AND DISCUSSION
One-sample t-tests showed that the DMN which included PCC/precuneus, medial prefrontal cortex and bilateral inferior parietal lobule, exhibited signifi cantly higher ReHo, ALFF and fALFF than the global mean (Figures 3A-C). This pattern was consistent with that in previous studies which showed high ReHo , ALFF and fALFF Zou et al., 2008;Zuo et al., 2010) within the DMN. Previous studies have consistently demonstrated that the DMN regions show task-independent deactivation across a wide range of cognitive tasks compared with the resting-state (Shulman et al., 1997;Binder et al., 1999;Mazoyer et al., 2001), and these areas had signifi cantly higher blood fl ow and oxygen consumption than the global mean value (Raichle et al., 2001). The current results that DMN showed high activity is consistent with the conclusion that these regions represent the functional core underlying resting brain dynamics (Ghosh et al., 2008;Honey et al., 2009;Zuo et al., 2010). In line with previous studies Zuo et al., 2010), we also found ALFF measure may be more susceptible to possible artifactual effect in the vicinity of blood vessels and the cerebral ventricles than fALFF. However, which of the two measures, ALFF vs. fALFF, is better for between-groups studies or sensitive to abnormal spontaneous brain activity needs to be further investigated.
FC analysis showed that the medial prefrontal cortex and bilateral inferior parietal lobule had signifi cantly positive correlation with PCC ( Figure 3D). The dorsal anterior cingulate cortex, bilateral insula, bilateral middle temporal cortex and bilateral dorsolateral FIGURE 3 | Within-condition patterns of ReHo (A), ALFF (B), fALFF (C) and PCC-FC (D). All these methods revealed the pattern of the default mode network. The numbers below the images refer to the MNI z coordinates. The statistical threshold was set at t > 3.89 (P < 0.0001) for one-tailed t-tests (for ReHo, ALFF and fALFF) and |t| > 4.08 (P < 0.0001) for two-tailed t-test (for FC) and cluster size >135 mm 3 , which corresponds to a corrected P < 0.0001. LH, left hemisphere; RH, right hemisphere. prefrontal cortex showed negative correlation with PCC. These results are consistent with previous studies that suggested a competitive relationship between the DMN and the anti-correlated network (Fox et al., 2005;Fransson, 2005;Long et al., 2008;Yan et al., 2009).
The illustrations of ReHo, ALFF, fALFF and FC analyses by using DPARSF pipeline analysis validated its correctness and demonstrated its effectiveness.