Data Analysis

For each patient, age-standardized scores were calculated for letter-based fluency and for category-based fluency. To do so, we first calculated the total number of words produced for the three letter cues (letter-based sum) and for the three category cues (category-based sum). We then converted these sums into age-standardized Z-scores, using published norms for adult Hebrew speakers (Kavé, 2005).

Prior to this calculation, fluency data were preprocessed as follows: Errors and repetitions were excluded. When two homonyms were provided, the second mention was counted only if the patient pointed out the alternate meaning explicitly. Words inflected in both masculine and feminine forms were counted as one. Synonyms were counted as separate words. In category-based fluency, names of subcategories (e.g., a bird) were not counted if specific items within that subcategory (e.g., dove, eagle) were also provided. Slang terms, as well as foreign words, were generally considered acceptable (Kavé, 2005).

Data preprocessing was conducted using the open sourced “mrDiffusion” package¹ and MATLAB 2012b (The Mathworks, Natick, MA, United States). Preprocessing followed the same standard steps as in our previous publications (Dougherty et al., 2007; Yeatman et al., 2011; Blecher et al., 2016; Kronfeld-Duenias et al., 2016), as detailed below.

First, T1 images were aligned to an ac-pc orientation: the locations of the anterior and posterior commissures were identified manually on the T1 of each patient and these points were used to align the anatomical T1 volume to a canonical ac–pc orientation, using a rigid body transformation (no warping was applied). Second, distortions in the diffusion-weighted images due to eddy currents and subject motion were corrected by a 14-parameter constrained non-linear co-registration algorithm based on the expected pattern of eddy-current distortions (Rohde et al., 2004). Third, diffusion images were registered to the ac-pc aligned T1 anatomical images. Alignment was achieved by registering the b0 images to the resampled T1 image using a rigid-body, mutual information maximization algorithm (implemented in SPM5; Friston and Ashburner, 2004). At this final alignment stage, the combined transform resulting from motion correction, eddy current correction and anatomical alignment was applied to the raw diffusion data once, and the data were resampled at exactly 2.6 mm isotropic voxels. Next, the table of gradient directions was appropriately adjusted to fit the resampled diffusion data (Leemans and Jones, 2009). Finally, we fitted a tensor model to the diffusion data in each voxel using a standard least-squares algorithm, and extracted the eigenvectors and eigenvalues (λ₁, λ₂, λ₃) of the tensor. Given the single b-value used (b = 1000), the tensor model is the most appropriate for the analysis of our data. Importantly, at this b-value, the tensor model provides high accuracy, similar to more complicated shapes (Rokem et al., 2015). Using the eigenvalues extracted from each tensor, we calculated the FA in each voxel as the weighted standard deviation of the three eigenvalues (Basser and Pierpaoli, 1996). Additional complementary measures were calculated, including axial diffusivity (AD, λ1) and radial diffusivity [RD (λ2+λ3)/2]. AD is defined as the diffusivity along the principal axis of diffusion, and RD as the average diffusivity along the two remaining minor axes. Note that these calculations all took place at the individual patient level in the native space of each patient.

We focused on a small set of preselected tracts, defined individually in each patient’s native space. Tracts of interest included the fronto-temporal arcuate fasciculus (AF_ft), frontal aslant tract (FAT), IFOF and uncinate fasciculus (UF), bilaterally. These tracts were selected based on their known involvement in phonological processing (AF_ft, e.g., Yeatman et al., 2011), semantic processing (IFOF, UF, e.g., Duffau, 2013; Nugiel et al., 2016), or oral fluency (FAT, Catani et al., 2013).

In order to identify these tracts and quantify their diffusion parameters, we used the Automatic Fiber Quantification (“AFQ”) package, an automated segmentation and quantification tool (Yeatman et al., 2012). AFQ consists of the following steps: (1) Whole brain fiber tractography, (2) Tract segmentation based on region-of-interest (ROI) and automatic cleaning of fiber outliers, and (3) Quantification of diffusion properties along the tracts. For whole brain tracking (step 1), we used deterministic Streamlines Tractography (STT), with a 4th Runge–Kutta path integration method and 1 mm fixed step size (Mori et al., 1999; Basser et al., 2000; Press et al., 2002). Deterministic tractography was used in order to avoid issues pertaining to tract selection (Pestilli et al., 2014) that are not well addressed with a single, relatively low b-value scan. Deterministic methods proved to be reliable for the purpose of identifying such large, well-known tracts as the ones identified here (Yeatman et al., 2012). Tract segmentation (step 2) was done in the native space of each patient, using ROIs defined on a T1 template (ICBM, 2009a Non-linear Asymmetric template; Fonov et al., 2011), which were back-transformed into the patient’s native space (see Figure 1 and Supplementary Figure S1 for the definition of all ROIs). Whole brain fibers were restricted to only those that passed through both ROIs, for each tract (following Wakana et al., 2007 for the AF_ft, IFOF and UF; see Kronfeld-Duenias et al., 2016 for the procedure to segment the FAT using AFQ). After tract segmentation, an automatic cleaning strategy was applied, removing fibers longer than 4 standard deviations from the mean fiber length and those that spatially deviated more than 5 standard deviations from the core of the tract (see Yeatman et al., 2012 for details regarding the automatic segmentation and cleaning procedures).

Segmented white matter tracts. Analyzed tracts are visualized in a single patient with MS (Male, 20 years old) overlaid on his T1 images. (A) Shown are three left hemisphere tracts: fronto-temporal arcuate fasciculus in blue, uncinate fasciculus in yellow, and inferior fronto-occipital fasciculus in red. (B) Left and right frontal aslant tracts are shown in green in the same participant. Dashed lines represent the location of the regions of interest (two for each tract) used to segment the fibers. AF_ft, fronto-temporal arcuate fasciculus; UF, uncinate fasciculus; IFOF, inferior fronto-occipital fasciculus; FAT, frontal aslant tract.

Quantification of diffusion properties (step 3) was applied as follows: For each participant, for each tract, a FA profile was calculated by sampling 100 equidistant nodes along the core of the tract, between the two ROIs used to segment it. Additionally, a mean tract-FA value was calculated by averaging over all the streamlines within each tract, end to end. The resulting FA-profile and mean tract-FA were subject to further statistical analyses (see section “Brain-Behavior Correlation Analysis” below). Tract profiles provide increased sensitivity to specific clusters of brain-behavior correlations (see, e.g., Travis et al., 2015; Blecher et al., 2016; Kronfeld-Duenias et al., 2016). For tract visualization, we used “Quench,” an interactive tract visualization tool (Akers, 2006).

MR images were analyzed by an experienced user (SM) to identify MS lesions using in-house developed lesion segmentation software (MS Analyze, MATLAB 7.5). The number and volume of brain lesions were calculated on axial T2 and T2-FLAIR images (slice thickness 3 mm; no gap). Lesions were identified on axial slices and assigned to specific lobes by comparison to a navigated anatomical MRI atlas².

We used the Kolmogorov–Smirnov test to assess the normality of the data (Corder and Foreman, 2009). Based on the results, Pearson correlation coefficients were calculated between mean-tract FA and each fluency measure, separately. We controlled the false discovery rate (FDR) across the 8 tracts of interest, at q = 0.05. Second, for each tract, Pearson correlation coefficients were calculated at each node along the trajectory of the tract. Significance was corrected for 800 correlations using a non-parametric permutation method, yielding a family-wise error (FWE) corrected alpha-value of 0.05 (Nichols and Holmes, 2002). This correction produced an FWE significant cluster size and a corrected alpha value for each tract of interest. We consider a segment significant if (1) the segment includes a cluster of adjacent nodes, each showing a correlation with p < 0.05 (uncorrected), and the cluster size is equal or larger than the critical cluster size determined by FWE; or (2) the segment includes any number of nodes that show a correlation with a P-value smaller than the corrected alpha determined by FWE. To visualize the pattern of co-variation between FA and fluency scores in significant segments, we extracted the mean FA value within the significant cluster for each patient [in case (2), this was achieved by defining a window of 17 nodes, centered on the most significant node]. Then, we plotted the data from the significant cluster or window against the relevant fluency score. A window size of 17 nodes was selected a priori as a reasonable size that balances generality and specificity, but very similar scatter plots were observed with cluster sizes of 13, 15, and 19 (not shown).

Significant correlations along the tract profiles were followed up with multiple regression models and partial correlations. Multiple regression models attempted to explain the variance in each fluency score using as predictors the mean FA in the significant cluster, together with age and education. Partial correlations considered the correlation between fluency and FA while controlling for (one at a time) a variety of demographic, cognitive and clinical factors, including age, gender, education, disease duration, executive function and attention. Executive function and attention were weakly correlated (r = 0.38), but we chose to enter them separately into this analysis because they stand for different cognitive components. We applied FDR correction at q < 0.05 for those 6 partial correlations. The FA-profile analysis was also followed up with partial correlations controlling for the number of lesions. This procedure resulted in three partial correlations for each significant segment (number of lesions in total, frontal and temporal areas). We applied FDR correction at q < 0.05 for those 3 × 2 partial correlations. In addition, significant correlations with one fluency score (letter-based or category-based) were followed up with partial correlations while controlling for the other fluency score.