Converting DICOM files into NIfTI images. The input files of PANDA can be in either DICOM or NIfTI format. If the input files are in NIfTI format, this conversion step will be skipped. Otherwise, DICOM files will be converted into NIfTI format during this step. The dcm2nii tool embedded in MRIcron accomplished this task.

Estimating the brain mask. This step yields the brain mask by using the bet command of FSL (Smith, 2002). The brain mask is required for the subsequent processing steps. Here, the b0 image without diffusion weighting was used for the estimation.

Cropping the raw images. To reduce the memory cost and speed up the processing in subsequent steps we cut off the non-brain space in the raw images, leading to a reduced image size. The acquired brain mask was used to determine the borders of the brain along the three dimensions. The fslroi command of FSL was then applied to remove the non-brain spaces.

Correcting for the eddy-current effect. Eddy-current induced distortion of diffusion weighted images (DWI), as well as simple head-motion during scanning, can be corrected by registering the DW images to the b0 image with an affine transformation. To achieve this, the flirt command of FSL was used. Notably, this registering procedure was applied to all images, with the b0 image of first acquisition used as the target if multiple DWI acquisitions were scanned. It is worth mentioning that while the eddy_correct command of FSL is not called here, the result of this step is exactly the same as the output of eddy_correct. Basically, PANDA just splits the 4D file (the input file of eddy_correct) into a number of 3D files and then performed the affine-registration exactly like eddy_correct. The purpose of this implementation is to avoid the large memory demand when the 4D file size is huge. Finally, the gradient direction of each DWI volume was rotated according to the resultant affine transformations (Leemans and Jones, 2009).

Averaging multiple acquisitions. This step will be skipped if there is only one DWI acquisition. Otherwise, after eddy-current correction, the aligned multiple DWI was averaged by calling the fslmaths command of FSL.

Calculating diffusion tensor (DT) metrics. This step involves a voxel-wise calculation of the tensor matrix and the DT metrics, including FA, MD, axial diffusivity (AD), and radial diffusivity (RD) (Pierpaoli and Basser, 1996; Song et al., 2002). The dtifit command of FSL was applied.

Normalizing. To allow for comparison across subjects, location correspondence has to be established. To end this, registration of all the individual images to a standardized template is always applied. Here, PANDA non-linearly registered individual FA images of native space to the FA template in the MNI space by calling the fnirt command of FSL. The resultant warping transformations were then used to resample the images of the diffusion metrics (i.e., FA, MD, AD, and RD) into the MNI space with a customized spatial resolution (e.g., 1 × 1 × 1 mm or 2 × 2 × 2 mm). This resampling step was implemented by the applywarp command of FSL.

Output for voxel-based analysis. The resultant images of the diffusion metrics in the standard space are ready for voxel-based statistical analysis. However, in the framework of voxel-based analysis, these images are typically smoothed to some degree, which can reduce the effect of image noise and misalignment between subjects. Accordingly, PANDA smoothed the images with a given Gaussian kernel, which was realized by calling the fslmaths command of FSL. The smoothed diffusion metric images can then be directly used for voxel-based statistical analysis with any preferred tools, e.g., FSL (http://www.fmrib.ox.ac.uk/fsl/), SPM (http://www.fil.ion.ucl.ac.uk/spm/), or AFNI (http://afni.nimh.nih.gov/afni/).

Output for atlas-based analysis. In addition to the popular voxel-based method of analysis, diffusion metrics can be analyzed at the level of region of interests (ROI), which may provide better statistical sensitivity in some cases (Faria et al., 2010). Recently, a few WM atlases (e.g., the ICBM-DTI-81 WM labels atlas and the JHU WM tractography atlas) have been proposed (Mori et al., 2008). These WM atlases in the standard space allow for parcellation of the entire WM into multiple ROIs, each representing a labeled region in the atlas. To support ROI-based analysis, PANDA calculates the regional diffusion metrics (i.e., FA, MD, AD, and RD) by averaging the values within each region of the WM atlases. These resultant ROI-based data (saved as text files) can be statistically analyzed with SPSS (http://www-01.ibm.com/software/analytics/spss/) and other statistical packages.

Output for TBSS-based analysis. The TBSS framework avoids the necessity of choosing a spatial smoothing procedure during voxel-based analysis and may provide better sensitivity and interpretability when it is applied to multi-subjects dMRI datasets (Smith et al., 2006). To support this type of analysis, PANDA follows the standard TBSS framework. Firstly, the mean of all the aligned FA images was created and skeletonized, resulting in a mean FA skeleton. Secondly, the diffusion metric data from individual subjects were projected onto the skeleton. Finally, individual images with data on the skeleton were created. The resultant images can be directly used for voxel-wise statistical analysis on the skeleton. Here, the fslmaths and tbss_skeleton commands of FSL were employed.

Two basic elements are required for a network: a node and a connection. Thus, the central tasks for constructing brain networks are: (1) defining network nodes and (2) defining connections between nodes. The schematic flowchart of network construction is demonstrated in Figure 2.

Defining network nodes. Typically, the entire brain is divided into multiple regions using a prior gray matter (GM) atlas, where each region represents a network node (Bullmore and Sporns, 2009). However, the prior atlases are generally defined in the standard space and need to be transformed to the native dMRI space of each individual. To address this, PANDA uses the framework proposed by Gong et al. (2009a). Specifically, the individual FA image in native space was co-registered to its corresponding structural image (i.e., T1-weighted) using an affine transformation. The individual structural image was then non-linearly registered to the ICBM152 template. Based on the resultant transformations in these two steps, an inverse warping transformation from the standard space to the native dMRI space can be obtained. Prior atlases in the standard space were then inversely warped back to individual native space by applying this inverse transformation. Currently, PANDA provides two well-defined atlases: the Automated Anatomical Labeling (AAL) (Tzourio-Mazoyer et al., 2002) atlas and the Harvard-Oxford atlas (HOA) (http://www.cma.mgh.harvard.edu/fsl_atlas.html). Notably, users can import customized atlases into PANDA to define the network nodes. During this step, the flirt, fnirt, inwarp, and applywarp commands of FSL were used.

Constructing networks using deterministic tractography. In general, deterministic tractography assumes a deterministic fiber orientation at every location during tracking, typically ending up with 3D trajectories for reconstructed WM tracts. Here, the dti_recon and dti_tracker commands of the Diffusion Toolkit (http://trackvis.org/dtk/) were called to reconstruct all possible fibers within the brain by seeding from all the WM voxels. For every pair of brain nodes/regions defined above, fibers with two end-points located in their respective masks were considered to link the two nodes. Based on the linking fibers, PANDA calculated three basic weighted matrices: number-weighted matrix (M^N), FA-weighted matrix (M^FA), and length-weighted matrix (M^L). In the matrices, each row or column represents a brain region/node. The values of the elements M(i, j)^N, M(i, j)^FA, and M(i, j)^L represent the number, averaged FA and averaged length of linking fibers between node i and node j, respectively. The resultant matrices were saved as a MATLAB data file and can be directly used for topological analysis with graph theoretic approaches (Bullmore and Sporns, 2009; Bullmore and Bassett, 2011).

Constructing networks using probabilistic tractography. In contrast, probabilistic tractography typically runs the tracking procedure many times, and fiber orientation is determined probabilistically. This type of tractography may improve tracking sensitivity, particularly for non-dominant fibers. The probabilistic tractography proposed by Behrens et al. (2003, 2007) has been implemented in FSL and is called by PANDA for network construction. This process involves two steps as follows:

BedpostX. Using the Markov Chain Monte Carlo sampling technique, this module estimated the local probability distribution of fiber direction at each voxel, a prerequisite for running subsequent probabilistic tractography (Behrens et al., 2003). In PANDA, bedpostX was realized by calling the xfibres command of FSL.

Probabilistic Tractography and Network Construction. Network construction using FSL-based probabilistic tractography has been previously described (Gong et al., 2009b). Briefly, for each defined brain region/node, probabilistic tractography was performed by seeding from all voxels of this region. For each voxel, 5000 fibers were sampled. To achieve this, the probtrackx command of FSL was called. The connectivity probability from the seed region i to another region j was defined by the number of fibers passing through region j divided by the total number of fibers sampled from region i. The connectivity probability of each node to the other nodes within the brain network can be calculated by repeating the tractography procedure for all nodes. This leads to an individual-specific weighted matrix, whose rows and columns represent the brain nodes and whose elements represent the connectivity probability between nodes. This matrix can also be directly used for various network analyses.

The test included data from 23 young adults (males, 11; females, 12; age, 17–24 years) and 17 elderly individuals (males, 8; females, 9; age, 54–77 years). All subjects were recruited from the campus and the local community. Subjects with a history of neurological or psychiatric disorders were excluded from this study. Written informed consent was obtained from each subject, and the protocol was approved by the Ethics Committee of the State Key Laboratory of Cognitive Neuroscience and Learning, Beijing Normal University.

All scans were performed using the 3-T Siemens Tim Trio MRI scanner in the Imaging Center for Brain Research, Beijing Normal University. Diffusion MRI was acquired using a single-shot echo-planar imaging-based sequence with following parameters: coverage of the whole brain; slice thickness, 2 mm; no gap; 68 axial slices; repetition time (TR), 9000 ms; echo time (TE), 92 ms; flip angle, 90°; 66 non-linear diffusion weighting directions with b = 1000 s/mm² and one image without diffusion weighting (i.e., b = 0 s/mm²); 4 repetitive acquisitions; acquisition matrix, 128 × 124; field of view (FOV), 256 × 248 mm²; resolution, 2 × 2 × 2 mm. Three-dimensional T1-weighted images with high resolution were obtained using a three-dimensional magnetization prepared rapid gradient echo (MP-RAGE) sequence with the following parameters: 1 mm slice thickness without gap; 176 sagittal slices; TR, 1900 ms; TE, 3.44 ms; flip angle, 9°; acquisition matrix, 256 × 256; FOV, 256 × 256 mm²; resolution, 1 × 1 × 1 mm.

The whole pipeline procedure of PANDA was run on all dMRI datasets with an in-house computing cluster of 6 nodes, each with 30GB of memory and 12 Intel Xeon E5649 2.53 GHz cores. For each pipeline step, default parameters were chosen.

Graph theoretical approaches have been applied to characterize the topology of brain networks that are derived from neuroimaging data (Bullmore and Sporns, 2009). Here, we focus on two topological network parameters: global efficiency and local efficiency. Global efficiency was defined as the average of the inverse of the “harmonic mean” of the characteristic path length, which represents global information transferring ability within the network (Latora and Marchiori, 2001). Local efficiency quantifies the ability of the network fault tolerant, corresponding to the efficiency of the information flow between nodal neighbors. Specifically, local efficiency was defined as the average of nodal local efficiency that is computed as the global efficiency of the sub-graph composed by its nearest neighbors (Latora and Marchiori, 2001).

For diffusion metric, we tested the group difference on FA across the entire WM. Specifically, normalized and smoothed (6 mm Gaussian kernel) FA images produced by PANDA were employed for this voxel-based analysis. A general linear model (GLM) with gender being taken as a covariate was applied to each WM voxel. For multiple comparison correction, false discovery ratio (FDR) was applied, and p < 0.01 was considered as significant.

For each subject, the FA-weighted matrix generated from PANDA was selected for topological analysis. Each matrix is 78 × 78 and represents the WM network of cerebral cortex. Each row or column of the matrix represents a cortical region of the AAL template (Gong et al., 2009a,b). The global efficiency and local efficiency were then calculated. To test the group effect on the global and local efficiency, a GLM with gender and brain size as covariates was applied, and p < 0.05 was chosen as the significant level.