Study data and participants

All data were acquired from the Human Connectome Project (HCP; http://www.humanconnectome.org/) open-access data initiative offering high-quality anatomical and functional MRI of the human brain. We used the HCP Young Adults (HCP-YA 1200 Subjects) data release as it contains a large sample of healthy adults for whom both language task fMRI and diffusion MRI sequences were acquired. The dataset comprised 1,200 healthy adults, aged 22–35 years. Each participant underwent an identical imaging protocol acquired on the same MRI scanner. Individuals with neuropsychiatric or neurologic disorders, diabetes, high blood pressure, premature birth, and severe symptoms associated with substance use were excluded from data collection (Van Essen et al., 2013). The present study focused on language fMRI and diffusion MRI data only. Individuals were only selected for inclusion if they had fMRI data available for the language story task (see below) and had corresponding 3T diffusion MRI data. This resulted in a sample size of 1,040 participants (562 females), with a mean age of 28.74 (SD = 3.69) years. According to the Edinburgh Handedness Inventory (Oldfield, 1971), 962 (92%) participants preferred their right hand, scoring at least 10 on a scale of −100 (left) to 100 (right). Eighty-five participants preferred left, scoring below −10, and two were ambidextrous, scoring zero.

Data acquisition

HCP data were acquired on a Siemens 3T Skyra system, with a 32-channel (SC72) head coil. Task fMRI data were collected using gradient-echo echo-planar imaging (EPI) with an isotropic resolution of 2.0 mm (TR, 720 ms; TE, 33.1 ms; matrix, 104 × 90, 72 slices; flip angle, 52°; BW, 2,290 Hz/Px; FOV, 208 × 180 mm, 72 slices; multiband accelerator factor, 8; Marcus et al., 2013). The HCP dMRI data were acquired using three shells (b = 1,000, 2,000, and 3,000 s/mm2) with 90 diffusion gradient directions and five b0 volumes with RL phase encoding direction (TE, 89.5 ms; TR, 5,520 ms; flip angles, 78/160°; isotropic voxel size, 1.25 mm3, multiband factor, 3; Sotiropoulos et al., 2013). A list of technical abbreviations is provided in Table 1.

Language paradigm

The language comprehension task used in the Human Connectome Project was designed by Binder et al. (2011). The task consists of two 3.8 min runs. Each run has four blocks of story tasks alternating with four blocks of math tasks. The story and math tasks are matched in terms of length, word and phoneme rate, speaking style, and prosodic features. The story blocks present subjects with 5–9 auditory sentences, followed by questions about the content of the story. The math task requires participants to perform arithmetic operations followed by equals and two choices. Since arithmetic tasks do not engage temporal lobe activity (Baldo and Dronkers, 2007), we decided to use a STORY–MATH contrast, as it effectively isolates regions responsible for language comprehension without “masking” temporal lobe activity. Additionally, the temporal lobe is involved in high-level processes of normal consciousness (Spitsyna et al., 2006); thus, we avoided passive tasks as a baseline to reduce the risk of masking activities in this region (which is essential for language comprehension). This contrast allowed us to cancel out the regions that are jointly activated in both tasks (such as low-level auditory and phonological input), isolating the regions involved in narrative processing including semantic and nonspeech-related aspects of language, theory of mind, and inference processing.

fMRI preprocessing and analysis

The preprocessed task fMRI data were retrieved from the HCP database (https://db.humanconnectome.org). The HCP preprocessing included fMRIVolume and fMRISurface pipelines, which were primarily built using tools from FSL (Jenkinson et al., 2012; http://www.fmrib.ox.ac.uk/fsl), FreeSurfer (Fischl, 2012), and the HCP Workbench (Marcus et al., 2013). Details of the preprocessing steps have been described previously (Glasser et al., 2013). The goal of the first fMRIVolume pipeline was to generate a 4D whole-brain time series. This was accomplished by (1) removing spatial distortions by gradient nonlinearity distortion correction, (2) realigning volumes using rigid-body motion correction using a single-band reference image as the target, and (3) estimating (using FSL toolbox “topup”) and correcting field map-based EPI distortions. The resulting EPI data were (4) registered to a T1-weighted scan and then (5) nonlinearly (FNIRT) to Montreal Neurological Institute (MNI) space, and (6) blood oxygenation level−dependent (BOLD) signal intensity was normalized by the average. This process resulted in individual subjects being mapped with a notable degree of left–right symmetry (Elam et al., 2021), which aligns with laterality research recommendations (Vingerhoets et al., 2023).

The goal of fMRISurface pipeline was to transform the resulting 4D time series into Connectivity Informatics Technology Initiative (CIFTI) grayordinate space, encompassing cortical, subcortical, and cerebellar gray matter collectively (Pham et al., 2022). This was accomplished by mapping fMRI data within cortical gray matter ribbon onto the native cortical surface, registering it to CIFTI grayordinate space (surface representation with 32,492 vertices on each hemisphere), and mapping the set of subcortical gray matter voxels from each subcortical parcel in each individual to a standard set of voxels in each atlas parcel, resulting in 2 mm average surface vertex and subcortical volume voxel spacing. Finally, grayordinate space data were smoothed using the Gaussian kernel.

We used a fully processed task-based STORY–MATH fMRI activation contrast of parameter estimates (COPE) map, which was generated by FSL FEAT and is readily available on https://db.humanconnectome.org as part of the “S1200 Subjects” dataset. Considering the spatial heterogeneity of the individual brain scans, the MSM-ALL (Multimodal Surface Matching) registered dataset was used, which uses information on areal features derived from the resting state network, myelin maps, and alignment of folding. The motivation for using MSM-ALL over MSM-SULC (cortical folding-based registration) came from previous studies that demonstrated a weaker correlation between sulcal depth and local curvature with regions responsible for higher cognitive functions, including Broca's area (Van Essen, 2005; Fischl et al., 2008), compared with the MSM-ALL registration, which showed improved cross-subject alignment of independent task fMRI datasets (Robinson et al., 2018).

Language comprehension laterality quotient

Grayordinates localized regions of interest (ROIs) on the “inflated” brain surface (Glasser et al., 2016; Van Essen and Glasser, 2016). A laterality quotient (LQ) was calculated to assess HLD for each participant's task fMRI activation using the CIFTI toolbox in MATLAB in ROIs associated with language comprehension. The analyses were conducted separately for frontal and temporal regions, given the well-documented phenomenon of crossed language dominance, where a participant may exhibit dominance in one hemisphere for frontal regions and the opposite hemisphere for temporal regions (Seghier, 2008). For the frontal ROIs, Brodmann areas 44 and 45 were selected due to their established high reliability in determining language dominance during semantic tasks (Sabbah et al., 2003; Seghier et al., 2008). In our temporal lobe laterality analyses, ROIs were chosen within the anterior temporal lobe (TGd, TGv, TE1a, TE2a, STGa, STSva, STSda) because these areas have been shown to be heavily involved in language comprehension (Binder et al., 2011; Fig. 1).

Figure 1.

Regions of interest (ROIs) were selected to calculate LQ based on the Jaccard–Tanimoto index (more details provided in text; Seghier, 2019). The LQ value is expressed as a percentage, ranging from −1 to 1. Values greater than 1/3 indicate left lateralization, values less than −1/3 indicate right lateralization, and values between −1/3 and 1/3 indicate bilateral orientation (Seghier, 2019). TE1, temporal area 1, found in middle temporal gyrus; STS, superior temporal sulcus; STG, superior temporal gyrus; TG, temporal gyrus; TE2, temporal area, including ventral and dorsal parts of inferior temporal gyrus.

To work with CIFTI files, we generated dscalar files for each ROI using wb_command, imported them into MATLAB, and extracted z-values from ROIs using the CIFTI toolbox. The z-values were thresholded for each participant by including only grayordinates with values greater than the median in each ROI (Dietz et al., 2016). To account for the unequal number of grayordinates between the left and right hemispheres (approximately 100 more grayordinates on the left than the right), we corrected these regional differences to ensure that comparisons between hemispheres were not skewed by differences in their sizes. This adjustment involved dividing the total sum of thresholded z-values by the number of grayordinates in each hemisphere for both frontal and temporal ROIs separately. The laterality quotient (LQ) was then computed for each participant's normalized z-values using the equations below:LQ=(L−R)max(L,R) where L and R represent the normalized z-values in the left and right ROIs, respectively. We chose to employ an innovative LQ formula based on the Jaccard–Tanimoto index to provide a more sophisticated approach in evaluating and classifying language lateralization (Seghier, 2019). This revised formula defines LQ as a metric of distance that adheres to a consistent distribution pattern, thus enhancing its sensitivity toward hemisphere activity differences, accentuating the distinctions between the two hemispheres. The values above +1/3 indicate left language dominance (LLD), values below −1/3 indicate right language dominance (RLD), and values between −1/3 and +1/3 indicate bilateral language representation (BLR), ensuring an equal cumulative probability in each dominance category.

Diffusion processing

Diffusion data were downloaded from the HCP S1200 Young Adult Data Release and preprocessed using the HCP Diffusion preprocessing pipeline using the FMRIB diffusion toolbox in FSL. Briefly, the pipeline included b0 image intensity normalization, removing EPI susceptibility-induced field distortions with FSL's “topup” algorithm (Andersson and Sotiropoulos, 2016), correcting for eddy current distortions, head movements, and gradient nonlinearities (Glasser et al., 2013). Quality control of the preprocessed diffusion MRI data was performed using DSI Studio software (http://dsi-studio.labsolver.org). An automatic quality control routine then checked the b-table to ensure its accuracy (Schilling et al., 2019). The diffusion data were coregistered in MNI space using q-space diffeomorphic reconstruction (Yeh and Tseng, 2011) to obtain the spin distribution function (SDF) with a recommended length ratio of 1.25, as specified in the original study (Yeh et al., 2010; Fig. 2).

Figure 2.

Flowchart of the methods pipeline. The preprocessed diffusion MRI data were reconstructed in an MNI space. The outputs of the reconstruction and SDFs were calculated to obtain the fiber orientations using DSI Studio. Then, two different approaches were used to examine the white matter tracts associated with language laterality. The connectometry approach involved obtaining a local connectome matrix and finding out its association with LQ. Shape analysis involved the recognition of the WM tracts using the HCP atlas and mapping eleven WM fiber bundles important for language function. The measures of key shape features, such as curl and volume, were extracted and linear regression analyses were used to look at the associations between shape metrics and LQ.

Connectometry analysis

We applied a whole-brain group connectometry analysis using DSI Studio as described in previous applications (Rahmani et al., 2017; Barnes-Davis et al., 2020, 2022; Dresang et al., 2021) to study the relationship between regional white matter quantitative anisotropy (QA) and language lateralization measures derived from LQs (Fig. 2). The connectometry approach derives the QA measure from the SDF in each fiber orientation, which defines the number of anisotropic spins along that direction in each streamline (Yeh et al., 2010, 2013). The anisotropy in each section of a white matter tract is then correlated with the study variable (Yeh et al., 2016). Unlike a voxel-based FA metric, which attributes identical anisotropy values to all fiber orientations within a voxel, QA demonstrates a discerning capability by identifying specific axonal orientations in each peak orientation of the SDF (Yeh et al., 2013).

Our connectometry analyses were conducted in two phases: the initial phase focused on examining the lateralization of frontal regions during language comprehension and the subsequent phase investigated temporal regions. Initially, connectometry analyses were performed on all participant groups concurrently, followed by post hoc analyses on three distinct groups separately to aid in interpretation and capture varying effects related to different degrees of laterality. Specifically, the first post hoc analysis included participants with LLD and BLR, the second consisted of participants with RLD and BLR, and the third included individuals with both LLD and RLD. The linear effect of handedness, sex, and age was mitigated using a partial linear correlation. A nonparametric Spearman partial correlation was used to derive the continuous segments correlating with an LQ (Yeh et al., 2016). Each reconstructed white matter tract within a voxel was tracked to extract a QA map for each participant (Yeh et al., 2013). A T-score threshold was assigned to the highest level of three to reduce the possibility of false positive results (Ashraf-Ganjouei et al., 2019). The tracks were filtered by topology-informed pruning with 16 iterations to remove implausible spurious connections (Yeh et al., 2019). Given the large sample size in our study, and to prevent false positives, a conservative false discovery rate (FDR) correction for multiple comparisons was employed with a threshold of 0.01 to select tracks showing significant associations between LQ and QA. To estimate the false discovery rate, 5,000 randomized permutations were applied to the group label to obtain the null distribution of the track length. After the correlational results were obtained, additional categorical analyses were performed at the group level (LLD/RLD, LLD/BLR, RLD/BLR). Short tracts (<20 mm) were removed for easier interpretation of our results.

Shape analysis

The SDF maps generated from the connectometry analysis were used for tract shape analysis, and automatic fiber tractography was performed using a deterministic fiber tracking algorithm utilizing DSI Studio software (Yeh, 2020). Eleven white matter tract bundles that are part of language comprehension networks (Friederici et al., 2007; Harvey et al., 2013; Rollans and Cummine, 2018; Shin et al., 2019; Ivanova et al., 2021; Forkel et al., 2022; Zhong et al., 2022) were then automatically tracked and recognized based on the HCP-842 tractography atlas (Yeh et al., 2018; Fig. 3). These include the arcuate fasciculus (AF), corpus callosum body, corpus callosum forceps major (splenium), corpus callosum forceps minor (genu), inferior fronto-occipital fasciculus (IFOF), frontal aslant tract (FAT), inferior longitudinal fasciculus (ILF), the three branches of the superior longitudinal fasciculus (dorsal SLF1, middle SLF2, and ventral SLF3), and the uncinate fasciculus. All white matter bundles were independently tracked within the left and right hemispheres, while the corpus callosum bundles were tracked as a whole. The diffusion sampling length ratio was set at 1.25, and the output resolution was resampled to 2 mm isotropic. To remove false connections, topology-informed pruning was applied with 32 iterations (Yeh et al., 2019). We decided to exclude participants for whom we could not reconstruct at least one of their ROI bundles. As a result, 290 participants were excluded, leaving us with a final sample size of 750 participants. Finally, after identifying all white matter tracts of interest, the following shape metrics were extracted: tract length, span, curl, elongation, diameter, volume, and surface area were extracted (Fig. 4).

Figure 3.

Eleven white matter tracts were reconstructed for shape analysis based on the HCP-842 atlas computed on 1,065 healthy people (Yeh et al., 2018).

Figure 4.

Schematic illustration of the shape analysis of the white matter tracts. a, The area metrics used in the included surface area (mm). b, The length metrics used in the study included mean tract length (mm) as well as span bundle (mm) and diameter (mm) of the bundle. c, The volume metrics used in the study included branch volume (mm3; blue dotted line), trunk volume (mm3; right dotted line), and total bundle volume (mm3; black dotted line). d, The shape metrics used in the study included curl and elongation.

To evaluate the statistical significance of differences among various laterality groups (LLD, RLD, BLR), we conducted an analysis of variance (ANOVA). This analysis utilized the same laterality groupings based on frontal and temporal ROIs and included covariates consistent with those used in the connectometry analyses. All computations were performed using R (version 4.4.1).