Structural and Functional Network-Level Reorganization in the Coding of Auditory Motion Directions and Sound Source Locations in the Absence of Vision

Auditory stimuli.

To infer the direction or location of sounds in the horizontal plane of the head, the auditory system extracts interaural differences in arrival time and sound level (Middlebrooks and Green, 1991; Blauert, 1997). Up-down localization notably relies on the interaction of sound waves within the pinnae, resulting in idiosyncratic direction-dependent spectral acoustic filters (Musicant and Butler, 1984; Wightman and Kistler, 1989; Middlebrooks and Green, 1992; Hofman et al., 1998). To create an externalized ecological sensation of sound source location and motion (Møller et al., 1996), we relied on individual in-ear stereo recordings that were recorded in a semianechoic room and from 30 loudspeakers on horizontal and vertical planes, mounted on two semicircular wooden structures with a radius of 1.1 m (Fig. 1A; Battal et al., 2020). Binaural in-ear recordings allow binaural properties, such as interaural time and intensity differences, as well as participant-specific monaural filtering cues, and serve to create reliable and ecological auditory space (Pavani et al., 2002). Participants were seated in the center of the apparatus with their head on a chin rest, such that the speakers on the horizontal plane were at the participant's ear level and those on the vertical plane were aligned with the participant's mid-sagittal plane.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Stimuli and experimental design. A, Acoustic apparatus used to present auditory moving and static sounds while binaural recordings were conducted from each participant's ear before the fMRI session. B, Auditory stimuli presented inside the fMRI consisted of the ollowing eight conditions: leftward (L), rightward (R), downward (D), and upward (U) moving sounds; and left, right, down, and up static sounds. Each condition was presented for 15 s (12 repetitions of 1250 ms sound, no ISI), and followed by a 7 s gap for indicating the corresponding direction/location in space and 8 s of silence (total interblock interval was 15 s). Sound presentation and response button press were pseudorandomized. Participants were asked to respond as accurately as possible during the gap period. C, The behavioral performance was recorded inside the scanner; error bars indicate the standard deviation, dashed line indicates the chance level.

Auditory stimuli were prepared using a custom-designed MATLAB script (R2013b, MathWorks). During the presentation of stimuli, the audio was recorded using binaural in-ear omnidirectional microphones (“flat” frequency range, 20–20,000 Hz; model SP-TFB-2, Sound Professionals) attached to a portable Zoom H4n digital wave recorder (16 bit, stereo, 44.1 kHz sampling rate). Microphones were positioned at the opening of a participant's left (L) and right (R) auditory ear canals. Then, these recordings were replayed to the participants when they were inside the functional MRI (fMRI) scanner. By using in-ear recordings, auditory stimuli automatically convolved with each individual's own pinna and head-related transfer function to produce a salient auditory perception in external space. The recorded auditory stimuli were used in both the main auditory experiment and the auditory motion localizer. All participants wore a blindfold throughout the auditory recordings and were instructed to keep their eyes closed. Before the recordings, the sound pressure level (SPL) was measured from the subject's head position and ensured that each speaker conveys 65 arbitrary dB SPL.

Stimuli recordings.

Sound stimuli consisted of 1250 ms pink noise (50 ms rise/fall time). In the motion condition, the same pink noise was presented moving in the following four directions: leftward, rightward, upward, and downward. Moving stimuli covered 120° of the space/visual field in horizontal and vertical axes. To create the perception of smooth motion, the 1250 ms of pink noise were fragmented into 15 equal length pieces with each 83.333 ms fragment being played every two speakers, and moved one speaker at a time, from the outer left to the outer right (rightward motion), or vice versa for the leftward motion. For example, for the rightward sweep, sound was played through speakers located at −60° and −52° consecutively, followed by −44°, and so on. A similar design was used for the vertical axis. This resulted in participants perceiving moving sweeps covering an arc of 120° in 1250 ms (speed, 96°/s; fade-in/fade-out, 50 ms) containing the same sounds for all four directions. The choice of the movement speed of the motion stimuli aimed to create a listening experience relevant to conditions of everyday life. Moreover, at such velocity it has been demonstrated that human listeners are not able to make the differences between concatenated static stimuli from motion stimuli elicited by a single moving object (Poirier et al., 2017), supporting the participant's report that our stimuli were perceived as smoothly moving (no perception of successive snapshots). In the static condition, the same pink noise was presented separately at one of the following four locations: left, right, up, and down. Static sounds were presented at the second-most outer speakers (−56° and +56° in the horizontal axis, and +56° and −56° in the vertical axis) to avoid possible reverberation differences at the outermost speakers. The static sounds were fixed at one location at a time instead of presented in multiple locations (Krumbholz et al., 2005; Smith et al., 2004, 2010; Poirier et al., 2017). This strategy was purposely adopted for two main reasons. First, randomly presented static sounds can evoke a robust sensation of auditory apparent motion (Strybel and Neale, 1994; Lakatos and Shepard, 1997; for review, see Carlile and Leung, 2016). Second, presenting static sounds located in a given position and moving sounds directed toward the same position allowed us to investigate whether moving and static sounds share a common representational format using multidimensional scaling (see below), which would have been impossible if the static sounds were randomly moving in space.

Participants were instructed to listen to the stimuli, without performing any task. Stimuli recordings were conducted in a session that lasted ∼10 min, requiring the participant to remain still during this period. All participants reported strong sensation of auditory motion and were able to detect directions and locations with high accuracy during the scanner session (Fig. 1C).

Auditory experiment.

To obtain blood oxygenation level-dependent activity related to auditory spatial processing, we computed single-subject statistical comparisons with fixed-effect general linear model (GLM). In the GLM, we used eight regressors from each condition (four motion directions, four sound source locations). The canonical double-gamma hemodynamic response function implemented in SPM8 was convolved with a boxcar function to model the above-mentioned regressors. Motion parameters derived from realignment of the functional volumes (three translational motion and three rotational motion parameters), button press, and four auditory response cue events were modeled as regressors of no interest. During the model estimation, the data were high-pass filtered with a cutoff at 128 s to remove the scanner drift and low-frequency fluctuations from the time series. To account for serial correlation because of noise in the fMRI signal, an autoregressive (AR 1) model was used.

At the fixed-effect individual-subject level, to obtain activity related to auditory processing in the whole brain, the contrasts tested the main effect of each condition, as follows: left motion, right motion, up motion, down motion, left static, right static, up static, and down static. Next, to identify regions responding preferentially to the auditory motion and static stimuli, we compared the response of all motion conditions to all static conditions (motion > static, and static > motion). These linear contrasts generated statistical parametric maps [SPM(t)] that were further spatially smoothed (Gaussian kernel, 8 mm FWHM) before being entered into a second-level group analysis, using a random effect model, accounting for intersubject variance.

At the group level, a series of one-sample t tests was implemented to examine the main effects of each condition (motion, static) and motion processing (motion > static) for each group. A conjunction analysis isolated brain areas jointly activated for the contrast motion > static in both groups (EB and SC). Two-sample t tests were then performed to compare these effects between groups (SC > EB, EB > SC).

Statistical inferences were performed using familywise error (FWE) correction for multiple comparisons using p < 0.05 over the entire brain volume or over small spherical volumes (radius, 15 mm) located around regions of interest (ROIs; Table 1) using a minimal cluster size threshold of 20 contiguous voxels (Worsley et al., 1996). Significant clusters were anatomically labeled using the xjView MATLAB toolbox (http://www.alivelearn.net/xjview) or structural neuroanatomy information provided in the Anatomy Toolbox (Eickhoff et al., 2007).

Independent visual and auditory motion localizer: region of interest definition.

We used independent auditory and visual motion localizers to functionally define hPT and hMT⁺/V5 regions. Preprocessing steps were similar to whole-brain univariate analysis (see section Univariate fMRI analysis). Single-subject statistical comparisons were made using a fixed-effect GLM for each participant with two regressors (motion, static), and motion parameters (six regressors of no interest). The canonical double-gamma hemodynamic response function implemented in SPM8 was convolved with a boxcar function for each regressor. Motion parameters derived from realignment of the functional volumes (three translational motion and three rotational motion parameters); button press was modeled as a regressor of no interest. During the model estimation, the data were high-pass filtered with cutoff of 128 s to remove the scanner drift and low-frequency fluctuations from the time series. To account for serial correlation because of noise in fMRI signal, autoregressive (AR 1) was used.

One-sample t tests were conducted to characterize the main effect of motion processing (motion > static). This linear contrast generated statistical parametric maps that were further spatially smoothed (Gaussian kernel 8 mm FWHM) and entered a second-level group analysis using a random-effects GLM. Group-level peak coordinates of bilateral hMT⁺/V5 and hPT were defined by contrasting the main effects of localizer scan (motion vs static), surviving a whole-brain familywise error correction (p < 0.05). Peak coordinates from the auditory and visual motion localizers were used to create a sphere of 6 mm radius (110 voxels) around the following four ROIs: left hMT⁺/V5, right hMT⁺/V5, left hPT, and right hPT. The four ROIs were defined functionally but constrained by anatomic landmarks of the regions. hPT was selected within the triangular region lying caudal to the Heschl's gyrus on the supratemporal plane, while hMT⁺/V5 was constrained with the ascending limb of the inferior temporal sulcus (Zeki et al., 1991; Watson et al., 1993).

To prevent the possible spurious overlap between the visual and auditory responsive regions arising from smoothing (Jiang et al., 2014), we performed the statistical inferences on the β-parameter estimates that were extracted from the unsmoothed data.

Because defining the visual motion area hMT⁺/V5 is obviously impossible in EB, we opted for the peak coordinates from the visual motion localizers in sighted subjects in normalized in MNI to be applied to blind people similarly. We, then, extracted data from the same spatially aligned ROIs in both groups to perform functional univariate and multivariate analyses.

Multiclass direction and location decoding.

To investigate motion direction and static location information in areas hMT⁺/V5 and hPT in sighted and blind participants, four classes of classifiers were trained and tested to discriminate between the response patterns of the four auditory motion directions and four sound source locations, respectively.

Preprocessing steps were identical to the steps performed for univariate analyses, except for functional volumes that were smoothed with a Gaussian kernel of 2 mm (FWHM). MVPA were performed in CoSMoMVPA [http://www.cosmomvpa.org/ (Oosterhof et al., 2016), which implements LIBSVM software (http://www.csie.ntu.edu.tw/∼cjlin/libsvm)]. A general linear model was implemented in SPM8, where each block was defined as a regressor of interest. A t-map was calculated for each block separately. Two multiclass linear support vector machine classifiers with a linear kernel with a fixed regularization parameter of C = 1 were trained and tested for each participant separately within each group. The two multiclass classifiers were trained and tested to discriminate between the response patterns of the four auditory motion directions and locations, respectively.

For each participant, the classifier was trained using a cross-validation leave-one-run-out procedure where training was performed with n – 1 runs and testing was then applied to the remaining run. In each cross-validation fold, an ANOVA-based feature selection was applied on the training data to select a subset of voxels (n = 110) that showed the most significant variation between the categories of stimuli (in our study, between orientations). The selected features were used for training and testing. This feature selection process not only ensures a similar number of voxels within a given region across participants, but, more importantly, identifies and selects voxels that are carrying the most relevant information across categories of stimuli (Cox and Savoy, 2003; De Martino et al., 2008); therefore, minimizing the chance to include voxels that carry noises unrelated to our categories of stimuli. The t-maps in the training set were normalized (z scored) across conditions, and the estimated parameters were applied to the test set. To evaluate the performance of the classifier and its generalization across all the data, the previous step was repeated 12 times where in each fold a different run was used as the testing data and the classifier was trained on the other 11 runs. For each region per subject, a single classification accuracy was obtained by averaging the accuracies of all cross-validation folds.

Binary direction and location decoding.

To investigate the preference of “axis of motion/space” in both hMT⁺/V5 and hPT, binary classifiers were used to discriminate brain activity patterns for motion direction within and across axes. Four binary classifiers were used to discriminate brain activity patterns for left versus right motion, left versus right static, up versus down motion, and up versus down static (hereafter called within-axis classification). We used eight additional classifiers to discriminate across axes [left vs up, left vs down, right vs up, and right vs down motion directions; left vs up, left vs down, right vs up, and right vs down sound source locations (hereafter called across-axes classification)].

To assess whether hMT⁺/V5 and hPT regions demonstrate axis-of-motion characteristic tuning for auditory motion, we averaged the two within-axes motion accuracies from leftward versus rightward, and upward versus downward binary classifications. We also averaged the motion accuracies from four across-axes motion binary classification. To assess the existence of axis-of-space, we performed averaging on the sound source location classification accuracies with two within-axes and four across-axes binary classifications to obtain one accuracy value per subject per axis.

Multidimensional scaling.

To visualize the pairwise decoding accuracies of motion directions and sound source locations, the 28 decoding accuracy values were represented in the form of a dissimilarity matrix. Each column and each row of the matrix represented one motion direction or sound location. The matrix elements were the pairwise decoding accuracies. The accuracy values were used as a dissimilarity measure. We used multidimensional scaling (MDS) to project the high-dimensional dissimilarity matrix space onto two dimensions with the pairwise decodings that were obtained from hMT⁺/V5 and hPT in both sighted and blind individuals. The 32 neural dissimilarity matrices (1 per participant) for each of the four ROIs were used as neural input for MDS visualization. Note that MDSs are for visualization purposes only and were not used for statistical inferences.

Statistical analysis.

Statistical analyses were performed in MATLAB (for nonparametric tests) and R [for repeated-measures ANOVAs, and linear mixed-effects models (LMMs)]. Tests for normality were conducted using a Shapiro–Wilk test. All the within-group comparisons were conducted using paired t tests.

Statistical significance in the multivariate classification analyses was assessed using nonparametric tests permuting condition labels and bootstrapping (Stelzer et al., 2013). Each permutation step included shuffling of condition labels and rerunning the classification, which was repeated 100 times on the single-subject level. Next, we applied a bootstrapping procedure to obtain a group-level null distribution that is representative of each group. For each group, from each subject's null distribution one value was randomly chosen and averaged across all the subjects. This step was repeated 100,000 times, resulting in a group-level null distribution of 100,000 values. The classification accuracies across subjects were considered as significant if the p value was <0.05 after corrections for multiple comparisons using the falser discovery rate (FDR) method (Benjamini and Yekutieli, 2001). The group comparison was also tested for significance by using permutation (100,000 iterations).

Classification accuracies were entered into a LMM, as computed through the lmer function in R (afex package; Singmann et al., 2021). We performed an LMM with group (EB, SC), condition (motion, static), region (hMT⁺/V5, hPT), and hemisphere (left, right) as fixed effects, and subject as a random effect.

To assess axis-of-motion preference across groups and across ROIs, the pairwise decoding accuracies were submitted to LMMs. For each ROI, we performed an LMM with group (EB, SC), condition (within axis, across axis), and hemisphere (left, right) as fixed effects, and subject as a random effect.

Preprocessing of diffusion data.

Data preprocessing was implemented in MRtrix 3.0 (Tournier et al., 2012; www.mrtrix.org), and in FSL version 5.0.9 (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSL). Data were denoised (Veraart et al., 2016), corrected for Gibbs ringing, Eddy current distortions, and head motion (Andersson and Sotiropoulos, 2016), and for low-frequency B1 field inhomogeneities (Tustison et al., 2010). Spatial resolution was upsampled by a factor of 2 in all three dimensions (1.15 mm isotropic) using cubic b-spline interpolation, and intensity normalization across subjects was performed by first deriving scale factors from the median intensity in select voxels of white matter, gray matter, and CSF in b = 0 s/mm² images, and then applying these across each subject image (Raffelt et al., 2012a). This step normalizes the median white matter b = 0 intensity (i.e., non-diffusion-weighted image) across participants so that the proportion of one tissue type within a voxel does not influence the diffusion-weighted signal in another. The T1-weighted structural images were nonlinearly registered to the diffusion data in ANTs (Avants et al., 2008) using an upsampled FA map (1 × 1×1 mm³) and segmented in maps for white matter, gray matter, CSF, and subcortical nuclei using the FAST algorithm in FSL (Zhang et al., 2001). This information was combined to form a five-tissue type image to be used for anatomically constrained tractography in MRtrix3 (Smith et al., 2012). These maps were used to estimate tissue-specific response functions (i.e., the signal expected for a voxel containing a single, coherently oriented fiber bundle) for gray matter, white matter, and CSF using multishell multitissue (MSMT)-constrained spherical deconvolution (CSD; Jeurissen et al., 2014). Fiber orientation distribution functions (FODFs) were then estimated using the obtained response function coefficients averaged across subjects to ensure subsequent differences in the amplitude of fODFs will only reflect differences in the diffusion-weighted signal. Note that by using MSMT-CSD in our single shell, data benefitted from the hard non-negativity constraint, which has been observed to lead to more robust outcomes (Jeurissen et al., 2014). Spatial correspondence between participants was achieved by generating a group-specific population template with an iterative registration and averaging approach using FOD images from all the participants (Raffelt et al., 2011). Each subject's FOD image was registered to the template using FOD-guided nonlinear registrations available in MRtrix (Raffelt et al., 2012b). These registration matrices were used to transform the seed and target regions from native diffusion space to template diffusion space, where tractography was conducted. We chose to conduct tractography in the template diffusion space, as FOD-derived metrics of microstructural diffusivity can only be computed in that space. Subsequently, we extracted the following quantitative measures of microstructural diffusivity for all the fixels (fiber populations within a voxel) in the brain: fiber density (FD), fiber-bundle cross-section (FC) and a combined measure of FD and cross-section (FDC; Raffelt et al., 2017). For further details on these metrics, see the FOD-derived microstructural metrics section.

Preparation of hMT⁺/V5 and hPT for tractography.

As for the fMRI analyses, hMT⁺/V5 was functionally defined using the visual motion localizers in the 13 SCs included in the analysis of the diffusion data. To functionally define hPT, we used the peak coordinate of the motion versus static contrast in the auditory experiment. To avoid the location of hPT being biased toward SC, we used the conjunction of the SC and EB participants who underwent the fMRI and dMRI acquisitions.

We computed the warping images between the standard MNI space and the native structural space of each participant by conducting a nonlinear registration in ANTs (Gurtubay-Antolin et al., 2021). Using these registration matrices, we transformed the group peak coordinates from the standard MNI space to the native structural space of each participant. We then transformed these coordinates from the native structural space to the native diffusion space using the same transformation used to register the T1-weighted structural images to the diffusion data (described in the diffusion data preprocessing section). Once in native diffusion space, the peak coordinates were moved to the closest white matter voxel (FA, >0.25; Blank et al., 2011; Benetti et al., 2018; Gurtubay-Antolin et al., 2021), and a sphere of 5 mm radius was centered there. To ensure that tracking was done from white matter voxels only, we masked the sphere with individual white matter masks. Last, ROIs were transformed from native diffusion space to template diffusion space, where tractography was conducted.

Tractography: hMT⁺/V5–hPT connections.

Probabilistic tractography between hMT⁺/V5 and hPT was conducted for both hemispheres. In addition to the hPT, a region just anterior to hMT⁺/V5 (called hMTa; Rezk et al., 2020) is also selectively recruited for the processing of moving sounds (Poirier et al., 2005; Saenz et al., 2008; Battal et al., 2019; Rezk et al., 2020). To avoid the notion that hMT⁺/V5–hPT connections could include fibers connecting hPT with hMTa (which lies near hMT⁺/V5 and responds to auditory motion), individually identified hMTa (as localized by auditory motion localizer task) was used as an exclusion mask.

We computed tractography in symmetric mode (i.e., seeding from one ROI and targeting the other, and conversely). We then merged the tractography results, pooling together the reconstructed streamlines. We used two tracking algorithms in MRtrix (“iFOD2” and “Null Distribution2”). The former is a conventional tracking algorithm, whereas the latter reconstructs streamlines by random tracking. The iFOD2 algorithm (second-order integration over FODs) uses a Bayesian approach to account for more than one fiber orientation within each voxel and takes as input an FOD image. Candidate streamline paths are drawn based on short curved “arcs” of a circle of fixed length (the step size), at a tangent to the current direction of tracking at the current points (Tournier et al., 2010). The Null Distribution2 algorithm reconstructs connections based on random orientation samples, identifying voxels where the diffusion data are providing more evidence of connection than that expected from random tracking (Morris et al., 2008). We used the following parameters for fiber tracking: randomly placed 5000 seeds for each voxel in the ROI, a step length of 0.6 mm, FOD amplitude cutoff of 0.05, and a maximum angle of 45° between successive steps (Tournier et al., 2010, 2012). We applied the anatomically constrained variation of this algorithm, whereby each participant's five-tissue type, segmented T1 image provided biologically realistic priors for streamline generation, reducing the likelihood of false positive results (Smith et al., 2012). The set of reconstructed connections was refined by removing streamlines whose length was 2.5 SDs longer than the mean streamline length or whose position was >2.5 SDs away from the mean position, as in previous studies (Gurtubay-Antolin et al., 2021; Yeatman et al., 2012; Takemura et al., 2016). To calculate the distance of a streamline from the core of the tract, we resampled each streamline to 100 equidistant nodes and treat the spread of coordinates at each node as a multivariate Gaussian. The core of the tract was calculated as the mean of the x, y, and z coordinates of each fiber at each node.

Statistical analysis.

Statistical significance in the presence of hMT⁺/V5–hPT connections was assessed comparing the number of streamlines reconstructed by random tracking (Null Distribution2 algorithm), with those generated by conventional tracking (iFOD2 algorithm; Morris et al., 2008; McFadyen et al., 2019; Gurtubay-Antolin et al., 2021). We calculated the logarithm of the number of streamlines [log(streamlines)] to increase the likelihood of obtaining a normal distribution (Müller-Axt et al., 2017; Tschentscher et al., 2019), which was tested before application of parametric statistics using the Shapiro–Wilk test in RStudio (Allaire, 2015). The log-transformed number of streamlines were compared using two-sided paired t tests. To control for unreliable connections, we calculated the group mean and SD of the log-transformed number of streamlines for each connection, and we discarded participants whose values were >2.5 SDs away from the group mean for the respective connection. Connections were only considered reliable when the number of streamlines reconstructed with the iFOD2 algorithm were higher than the ones obtained with the null distribution algorithm. Significance was thresholded at p = 0.05, Bonferroni corrected for multiple comparisons (p = 0.025, two hemispheres).

Whole-tract analysis.

Differences between SCs and EBs in hMT⁺/V5–hPT connections were assessed using macrostructural and microstructural measures. Macrostructural characterization included the trajectory and the connectivity index, whereas we relied on diffusivity measures derived from fiber orientation distribution for the microstructure assessment.

To assess whether hMT⁺/V5–hPT connections followed the same trajectory in the SC and EB groups, we computed their spatial overlap by means of the Dice similarity coefficient (DSC) (Dice, 1945). The DSC measures the spatial overlap between regions A and B, and it is defined as DSC (A, B) = 2(A ∩ B)/(A + B) where ∩ is the intersection. We calculated the DSC of hMT⁺/V5–hPT connections, using as region A the sum of binarized tract density images of hMT⁺/V5–hPT connections (thresholded at six subjects) in standard space and in the SC group. Region B was the analogous image in the EB group.

The connectivity index was determined by the log-transformed number of streamlines from the seed that reached the target divided by the log-transformed product of the generated sample streamlines in each seed/target voxel (5000) and the number of voxels in the respective seed/target mask (Müller-Axt et al., 2017; Tschentscher et al., 2019; Gurtubay-Antolin et al., 2021), as follows: Connectivity index (CI) =log(streamlines)log (5000 * (Vseed + Vtarget)).

Microstructure assessment relied on the FOD-derived quantitative metrics FD, FC, and FDC (fixel-based analysis pipeline in www.mrtrix.org). These estimates characterize a microstructure isolating the contribution of crossing fibers within the same voxels of interest, providing more accurate values than tensor-derived metrics such as fractional anisotropy (Behrens et al., 2007; Tournier et al., 2007; Raffelt et al., 2017). Briefly, FD is a quantitative measure of fiber density within a voxel, given that the integral of the FOD along a particular direction is proportional to the intra-axonal volume of axons aligned in that direction. It was calculated as the integral of the FOD along a particular direction, and it is sensitive to alterations at the fixel level (Raffelt et al., 2017). FC reflects changes in the intra-axonal volume of a fiber bundle that are manifested as a difference in the number of fixels in the spatial extent that the fiber bundle occupies (cross-sectional area). FCs for each fixel were calculated in the plane perpendicular to the fixel direction and were estimated by using the nonlinear warps used to register each subject's FOD to the template. Last, multiplying FD and FC, we computed the metric FDC, which combines both sources of information (within-voxel FD and changes in FC). These estimates were computed in all the voxels that were crossed by fibers connecting hMT⁺/V5 and hPT and then averaged for each bundle.

All measures (i.e., connectivity index, FD, FC, FDC) were compared between SCs and EBs using two-sided paired t tests, after testing for normality using the Shapiro–Wilk test. Significance threshold was set at p = 0.006 (α = 0.05 Bonferroni corrected for 2 hemispheres and four metrics) and subjects' values were >2.5 SDs away from the group mean were considered outliers. A possible lateralization of group differences in the microstructural metrics was tested by a 2 × 2 ANOVA with Group as a between-subject factor (2 levels: SC, EB) and Hemisphere as a within-subject factor (2 levels: R, L)

Along-tract analysis.

An along-tract analysis was performed on the reconstructed tracts using the Along Tract Statistics toolbox (Colby et al., 2012) in MATLAB (R2016b). A mean tract geometry (MTG) was estimated for each tract after resampling the streamlines to 50 points, so that all streamlines were composed by the same number of points. For each tract and participant, the mean of diffusion scalar values of corresponding points (i.e., streamline points assigned to the same MTG point) was computed per each MTG point, and individual along-tract profiles were estimated per each tract. To control for the possibility that group differences could be driven by scalar outliers, we excluded tracts whose mean scalar values deviated >2.5 SDs from the mean of each subject or deviated for >10 consecutive points from the group mean. To further exclude aberrant values, we rejected participants presenting profiles in which at five consecutive positions along the tract the profile values deviated ±2.5 SDs for the respective scalar analysis (in a group-wise fashion; Novello et al., 2018).

Per each of the considered diffusion scalars, a linear mixed-effects model was then adopted with fixed effects of group, position (i.e., MTG point), and a group and position interaction to assess possible diffusion scalar changes affecting the two groups differently. A group and position interaction was considered significant for p < 0.008 (α = 0.05. Bonferroni-corrected two hemispheres and three diffusion scalars). We adopted a permutation (n = 5000) approach to control the type 1 error and adjust p-values accordingly: a null distribution of maximum T-values across all MTG points under the null hypothesis of no group differences was empirically estimated by permuting group labels and recording the maximum observed T-statistics. For each MTG point, the group analysis T-statistics were then compared against the null distribution to get an adjusted p-value (Colby et al., 2012). Results were then considered significant for an FEW-corrected p value <0.05.