Resting functional connectivity of language networks: characterization and reproducibility.

The neural basis of language comprehension and production has been associated with superior temporal (Wernicke's) and inferior frontal (Broca's) cortical areas, respectively. However, recent resting-state functional connectivity (RSFC) and lesion studies have implicated a more extended network in language processing. Using a large RSFC data set from 970 healthy subjects and seed regions in Broca's and Wernicke's, we recapitulate this extended network that includes not only adjoining prefrontal, temporal and parietal regions but also bilateral caudate and left putamen/globus pallidus and subthalamic nucleus. We also show that the language network has predominance of short-range functional connectivity (except posterior Wernicke's area that exhibited predominant long-range connectivity), which is consistent with reliance on local processing. Predominantly, long-range connectivity was left lateralized (except anterior Wernicke's area that exhibited rightward lateralization). The language network also exhibited anti-correlated activity with auditory (only for Wernicke's area) and visual cortices that suggests integrated sequential activity with regions involved with listening or reading words. Assessment of the intra-subject's reproducibility of this network and its characterization in individuals with language dysfunction is required to determine its potential as a biomarker for language disorders.

INTRODUCTION

Superior temporal (Wernicke’s area) and inferior frontal (Broca’s area) cortices have been classically associated with language comprehension and production. However, lesion[1] and functional magnetic resonance imaging (fMRI) studies[2] have identified additional temporal, parietal and prefrontal regions, supporting the involvement of a more extended language network.[3,4] This network seems to be organized around a central axis of at least two interconnected heteromodal epicenters (Wernicke’s and Broca’s areas),[5] and abnormalities in its flexible parallel architecture might help explain a variety of clinical manifestations in language disorders (aphasia).[6]

Wernickes’ area (Brodmann areas, BAs, 22, 39 and 40) is traditionally associated with language comprehension and its damage results in Wernicke’s aphasia (receptive or fluent aphasia). Broca’s area (posterior inferior frontal gyrus; BA 45 and 44) is traditionally associated with language production, and its damage results in Broca’s aphasia (expressive or non-fluent or agrammatic aphasia).

Lesion and fMRI studies in healthy subjects indicate that speech comprehension and production is lateralized to the left brain hemisphere,[4] and departures from this leftward asymmetry have been implicated in developmental dyslexia,[7] schizophrenia[8] and autism spectrum disorder (ASD).[9] There is also preliminary evidence that heteromodal association cortices in the language network[10] might have a lower proportion of long-range connectivity axons.[11] Thus, a better characterization of the functional connectivity of Broca’s and Wernicke’s areas and their laterality patterns might help in the development of biomarkers for language disorders.

The low-frequency fluctuations of spontaneous brain activity might play an important role in the language network.[12] Brief (5–10 minutes) MRI scanning captures this intrinsic low-frequency activity and is now extensively used to evaluate the “resting-state” functional connectivity (RSFC) among brain regions.[13] The increasing availability of large public RSFC databases allow multicenter studies evaluating brain network properties in large samples with increased statistical power.[14,15] This method was recently used to evaluate learning effects[16] and the functional organization[17,18] of the language network. However, the reproducibility of RSFC from Broca’s and Wernicke’s areas, and the relative contribution of short- and long-range connections in the language network are still largely unknown.

Recently, we proposed functional connectivity density mapping (FCDM), a voxelwise data-driven technique to measure the short-range (local) functional connectivity density (FCD).[19] Here we used a standard seed-voxel correlation approach to map the strength of the RSFC with Wernicke’s and Broca’s areas, as well as FCDM to map short- and long-range FCD in 970 subjects from a large public MRI database.[14] We assessed the variability of the RSFC patterns across subjects and institutions and mapped the laterality of short- and long-range FCD patterns from Broca’s and Wernicke’s areas. We hypothesized that the functional connectivity of the Wernicke’s area would include a bilateral temporal parietal network (BAs 20, 22, 39, and 40), and that of Broca’s area would be left-lateralized and include inferior frontal cortical regions (BAs 44, 45 and 47). We further hypothesized that this pattern would be highly reproducible across research institutions and would exhibit increased contribution of short-range FCD.

METHODS

Subjects

Functional scans collected in “resting-state” conditions from 970 healthy subjects (Table 1) obtained from 22 research sites of the “1000 Functional Connectomes” Project (http://www.nitrc.org/projects/fcon_1000/) were included in the study. Datasets from other research sites that were not available at the time of the study (pending verification of IRB status) did not report demographic variables (gender and age), exhibited image artifacts that prevented the study of short-and long-range FCD, or did not meet the imaging acquisition criteria (3s ≥ TR, full brain coverage, time points > 100, spatial resolution better than 4-mm) were not included in the study.

Table 1

Demographic data, acquisition parameters of the resting-state functional connectivity (RSFC) datasets.

Dataset	Subjects	Age [years]	B [T]	Tp	TR [s]
Ann Arbor	20M/3F	15–41	3	295	1.0
Baltimore	8M/15F	20–40	3	123	2.5
Bangor	20M/0F	19–38	3	265	2.0
Beijing	76M/122F	18–26	3	225	2.0
Berlin	13M/13F	23–44	3	195	2.3
Cambridge	75M/123F	18–30	3	119	3.0
Cleveland	11M/20F	24–60	3	127	2.8
Dallas	12M/12F	20–71	3	115	2.0
ICBM	17M/25F	19–85	3	128	2.0
Leiden	23M/8F	20–27	3	215	2.2
Leipzig	16M/21F	20–42	3	195	2.3
MIT	17M/18F	20–32	3	145	2.0
Newark	9M/10F	21–39	3	135	2.0
New York	8M/12F	18–46	3	175	2.0
Ontario	11 subj.	N/A	4	105	3.0
Orangeburg	15M/5F	20–55	1.5	165	2.0
Oulu	37M/65F	20–23	1.5	245	1.8
Oxford	12M/10F	20–35	3	175	2.0
Palo Alto	2M/15F	22–46	3	235	2.0
Queensland	11M/7F	21–34	3	190	2.1
Saint Louis	14M/17F	21–29	3	127	2.5
TaipeiA	13 subj.	N/A	3	295	2.0
TaipeiB	8 subj.	N/A	3	175	2.0

M: males; F: females; B: magnetic field strength; T: tesla; Tp: time points; TR: repetition time.

Image preprocessing

The statistical parametric mapping package SPM2 (Wellcome Trust Centre for Neuroimaging, London, UK) was used for image realignment and spatial normalization to the stereotactic space of the Montreal Neurological Institute (MNI). For this purpose a 12-parameters affine transformation with medium regularization, 16-nonlinear iterations, voxel size of 3 × 3 × 3 mm3 and the EPI.mnc template included in the SPM2 package were used. Other preprocessing steps were carried out using IDL (ITT Visual Information Solutions, Boulder, CO). A multilinear regression approach that used the time-varying realignment parameters (3 translations and 3 rotations) was applied to minimize motion related fluctuations in the MRI signals,[19] and the global signal intensity was normalized across time points. Band-pass temporal filtering (0.01–0.10 Hz) was used to remove magnetic field drifts of the scanner and minimize physiologic noise of high frequency components.[19] Voxels with signal-to-noise (as a function of time) < 50 were eliminated to minimize unwanted effects from susceptibility-related signal-loss artifacts. These MRI time series, MRI(t), reflecting the preprocessing steps were saved in hard drive for subsequent analyses.

RSFC analysis

We aimed to map networks functionally connected to Broca’s and Wernicke’s areas from MRI(t) using a standard seed-voxel correlation approach. Two cortical voxels were selected using probabilistic atlases of human brain anatomy[20,21] as seeds for subsequent correlation analyses (Fig 1): a voxel near the center of mass of the left pars triangularis (BA 45), xyz = (−51, 27, 18) mm, was selected to represent Broca’s area; a voxel in the left supramarginal gyrus (at the boundaries of BAs 39, 40 and 20), xyz = (−51, −51, 30) mm, was selected to represent Wernicke’s area. Since the boundaries of Broca’s and Wernicke’s areas vary across individuals we also defined cubic regions-of-interest (ROI) with 125 voxels (volume = 3.375 cc) centered at the coordinates of the abovementioned voxels for complementary correlation analyses addressing the robustness of the results as a function of the location of the seeds (Fig 1). The strength of the functional connectivity was computed through Pearson correlations between time-varying signals at the seed location and those at other brain voxels; the average signal in the ROI was used for the complementary ROI-voxel correlation analyses. The Fisher transform was used to convert the step distributed correlation coefficient, R, into normally.distributed correlation coefficients (Broca: RB; Wernicke: RW). These rescaled correlation maps were spatially smoothed using an 8-mm Gaussian kernel to minimize the differences in the functional anatomy of the brain across subjects. Group analyses of functional connectivity included all 970 subjects and were carried independently for each seed region in SPM2 using t-test. Clusters with pcorr < 0.05, corrected for multiple comparisons using a family-wise error (FWE) threshold, were considered significant for group analysis in SPM.

Fig 1

Location of seeds (intersecting blue lines) and ROIs (red cubes) used to represent Broca’s and Wernicke’s areas superimposed on three orthogonal views of the human brain (B and E) and on standard probabilistic atlases depicting the anatomical location of Brodmann areas (A and D) and of the gyri implicated in the language network (C and F).

Combined network

Brain voxels significantly correlated with both Broca’s and Wernicke’s seeds were used in subsequent analyses of the variability of the connectivity measures across research sites. A simple t-test, including all rescaled RB and RW maps, was used to identify voxels significantly correlated with both seeds and a binary mask of the combined network was created using a t-score threshold t > 14 (p< 10−31).

Network modularity

We used a weighted-connectivity null model[22] to further explore the optimal partitions of the language network. Specifically, Pearson correlation coefficient were computed among clusters identified by the RSFC statistical analysis, independently for each subject, from average time-varying signals, averaged in 9-mm isotropic ROIs centered at the cluster coordinates in Table 2; then the average connectivity matrix, M, was computed across all subjects. The Brain Connectivity Toolbox (https://sites.google.com/a/brain-connectivity-toolbox.net/bct/Home) was used in order to quantify the degree to which the language network can be subdivided into clearly delineated groups that reflect the optimal network structure. A Louvain iterative algorithm[23] was used to find the initial partition vector corresponding to the undirected and weighted correlation matrix M. Thus, the modularity of the language network was optimized iteratively in a way that the number of optimal partitions was found automatically.[22]

Table 2

MNI coordinates and statistical significance (t-scores) for relevant clusters showing significant correlated and anticorrelated activity with Broca’s, RB, and Wernicke’s, RW, seeds as well as the statistical significance for the differential RSFC coupling with the seed regions, RB–Rw, the differential functional connectivity density (FCD) between short-range (SR) and long-range (LR) connections, and the lateralization index (LI) of long-range FCD.

Brain region	BA	X [mm]	y [mm]	z [mm]	RB	RW	RB-RW	FCDSR > LR	LI (LR)	LI (SR)
CORRELATED
Wernicke’s area	39/40	−51	−51	30	15.8	63.2	−41.4	8.4	−9.0	−5.9
Inferior parietal	40	57	−51	36	NS	41.9	−30.2	9.4	7.5	11.9
Broca’s area	45	−51	27	18	69.7	16.5	55.9	12.6	18.1	NS
Pars triangularis	45	51	30	18	27.8	6.3	20.0	20.9	−15.8	NS
Middle frontal	46	−39	18	45	18.6	34.5	−12.2	5.9	8.3	NS
Pars opercularis	44	42	21	42	7.1	27.2	−15.6	17.8	−8.8	NS
Pars orbitalis	47	−45	39	−12	29.7	29.0	NS	6.6	10.5	NS
Pars orbitalis	47	45	39	−15	15.0	20.3	NS	14.4	−7.6	5.6
Inferior temporal	21/20	−57	−30	−15	10.6	26.6	−12.8	10.5	NS	6.1
Inferior temporal	21/20	63	−30	−12	5.7	22.2	−12.9	12.3	NS	NS
Superior frontal	8	−3	36	45	22.5	31.3	−5.7	7.9	−6.6	−7.9
Caudate		−12	9	15	17.5	19.2	NS	24.7	−9.5	NS
Caudate		12	12	12	13.3	16.8	NS	21.7	8.6	NS
Putamen/globus pallidus		−18	0	9	13.9	12.1	NS	22.8	NS	NS
Ventral Thalamus		−9	−9	0	9.7	12.8	NS	17.4	−5.7	NS
Cerebellum	crus	15	−81	−30	18.9	22.6	NS	15.5	NS	NS
ANTICORRELATED
Striate	17	6	−75	−6	−14.9	−21.4	5.8	−9.0	NS	−8.8
Extrastriate	18	21	−69	−15	−14.1	−20.2	5.4	−8.4	NS	NS
Posterior parietal	7	6	−81	45	−16.5	−11.3	NS	−18.9	−8.9	−17.7
Superior Parietal	5	3	−51	57	−15.5	−9.4	NS	NS	NS	NS
Superior temporal	42	−63	−18	9	NS	−18.1	11.7	−5.4	−9.9	−17.4
Superior temporal	42	60	−21	12	NS	−15.0	11.7	NS	NS	8.1
Cingulate	24	0	0	48	−6.5	−15.0	6.4	NS	NS	NS

BA: Brodmann area. Statistical analyses: one-way ANOVA model in SPM2.

Functional connectivity density (FCD)

The preprocessed image time series MRI(t) underwent FCDM[19,24] to compute the strength of the lFCD and the global functional connectivity density (gFCD). The number of functional connections, k(x0), was determined through Pearson correlations between time-varying signals at x0 and those in other voxels using an arbitrary threshold R > 0.6; this correlation threshold was selected in our previous work because R < 0.4 increased false positive rate and CPU time, and that R > 0.7 leaded to lFCD maps with reduced dynamic range and lower sensitivity; thus we fixed R = 0.6 for all calculations.[19]

The lFCD at a given voxel x0 was computed as the local k(x0) between x0 and its neighbor voxels using a “growing” algorithm developed in IDL. Specifically, A voxel (xj) was added to the list of voxels functionally connected with x0 (xN; N = {i}) only if it was adjacent to a voxel that was linked to x0 by a continuous path of functionally connected voxels and R0j > 0.6. This calculation was repeated for all voxels that were adjacent to voxels that belong to the list of voxel functionally connected to x0 in an iterative manner until no new voxels could be added to the list. The lFCD at x0 was computed as the number of elements in the local functional connectivity cluster, k(x0). Then the calculation is initiated for a different x0. Whereas this calculation is performed for all N voxels in the brain, the number of necessary correlations to compute a map of the lFCD is reduced by a large factor (~ 1000).

The gFCD at x0 was computed as the global k(x0) between x0 and all other N = 57,712 voxels in the brain. Two voxels were considered functionally connected if the correlation factor R > 0.6, consistently with the computation of the lFCD (note that we verified that the gFCD was stable under threshold variations in the range 0.5 < R < 0.7). This calculation was repeated for all x0 voxels in the brain and involved the computation of a N(N−1)/2 correlations.

Parallelization of gFCD

A simple approach based on data parallelism was implemented to speed up the calculation of the gFCD by taking advantage of multiprocessor computer architectures. The data was distributed across processing nodes in a way that all nodes executed the same 32-bit code on different imaging data (two subjects per core). The 188kB standalone code that efficiently allocated 128MB of dynamic memory via pointer reference was developed and compiled using Visual C++ 6.0 (Microsoft Corp, Redmond, WA). The computation of the gFCD maps required only 120 minutes per subject in a single core PC, which was significantly faster than what could be achieved using user-friendly high-level programming languages such as IDL or Matlab (~ few days per subject).

A workstation with two Intel® Xeon® X5680 processors (six cores per processor, 12MB L3 Cache, 64-bit, 3.33 GHz), which accounted a total of 12-cores and allows 24 processing threads with Hyper-Threading Technology (HTT), running Windows 7 (64 bit) was used to compute the gFCD maps for each subject. Note that HTT enabled the operating system to address two virtual processors for each of the 12 cores that were physically present. Twenty-four batches of jobs were submitted simultaneously by terminal commands; Windows 7 distributed the workflow among the 24 virtual processors. This simple but efficient parallelization approach allowed computation of 24 subjects (one subject per virtual processor) at once and required in average only five minutes/subject to complete.

Short- and long-range FCD

Since lFCD predominantly reflects the regional short-range functional connectivity we defined: short-range FCD = lFCD. Since the gFCD included both local and distal connections we defined long-range FCD = gFCD – lFCD.[25] These maps inherited the radiological convention (left is right) of the native images. Short- and long-range FCD maps were spatially smoothed (8-mm) in SPM to minimize the differences in the functional anatomy of the brain across subjects. The short- and long-range FCD distribution were rescaled to their average values in the whole-brain in order to minimize variability across institutions. A one-way analysis of variance (ANOVA) SPM model with two conditions (short-and long-range FCD) was used for statistical analyses of FCD and a conservative PFWE < 0.05 was used for statistical significance.

Laterality index

Short- and long-range FCD maps with neurological convention (R: right is right) were additionally created by flipping (as mirror-reversal of) the radiological FCD maps (L) across the x-axis. Then, the strength and direction of short- and long-range FCD asymmetries were computed voxel-by-voxel using a laterality index map[26,27]: which partially accounts for inter-subject variability and differential acquisition protocols among research sites. Negative LI values indicate leftward asymmetry and positive LI values indicate rightward asymmetry. The Fisher transform was used to normalize the step distributed LI values. The short- and long-range LI maps were spatially smoothed (8-mm). Two-way ANOVA with two conditions (short- and long-range FCD) and two groups (males and females) was used for voxelwise SPM analyses of FCD laterality. Voxels with a conservative pFWE < 0.05 were considered significantly lateralized across subjects.

Anatomical labeling was based on the Automated Anatomical Labeling (AAL) atlas and the population-average landmark- and surface-based atlas of the cerebral cortex provided in the installation package of the MRIcron image viewer (http://www.nitrc.org/projects/mricron).

RESULTS

Seed-voxel correlation patterns

The correlation patterns across 970 subjects involved prefrontal, temporal and parietal cortices and were similar for all 22 research sites, independently for the Broca and Wenicke seeds (Fig 2 and Table 2). Independent t-tests for each seed demonstrated the high significance of the correlations between time varying signals at the seeds and those in brain regions classically associated with language processing (Figs 2A and 2B). Simple pairwise correlation analyses demonstrated a striking reproducibility of the patterns across all research institutions (R > 0.64; Fig 3). Specifically, the time varying signals in the Broca’s seed were positively correlated with those in a bilateral network that comprised Wernicke’s area, as well as prefrontal regions (pars triangularis, opercularis and orbitalis, and inferior, middle and superior frontal gyri) and subcortical brain regions (bilateral caudate, left putamen, ventral thalamus and right cerebellum crus) and were negatively correlated with signals in visual, somatosensory and anterior cingulate cortices (Table 2, Fig 4A). The MRI signals in the Wernicke’s seed were correlated with the same regions and, additionally, negatively correlated with those in the primary auditory cortex (Table 2, Fig 4B). Figures 4C and 4D show that the normalized correlation coefficients (Fisher’s Z) in the positively correlated network had approximately normal distribution across subjects; the normalized correlation coefficients in the negatively correlated network also exhibited approximate normal distributions (not shown). A t-test analysis of all correlation maps identified a left-lateralized network of regions that were positively correlated with both seeds. A binary mask of the combined network (Fig 5A, solid line) was created for subsequent analyses that evaluated the variability of this network pattern across research sites (SI Fig 1A; T-score > 14, P < 10−31; t-test, df = 1939). A direct comparison of RSFC patterns revealed that the correlations were stronger in inferior parietal cortex, middle frontal gyrus, pars opercularis, and inferior temporal gyrus and weaker in pars triangularis for the Wernicke seed than for the Broca seed (PFWE < 0.05; Fig 5B and Table 2 and SI Fig 1B; t-test).

Fig 2

RSFC maps reflecting the average spatial distribution of temporal correlations with time-varying signals in (A) left Broca’s area (RB) and (B) left Wernicke’s area across subjects for each research site in Table 1. Population-average landmark- and surface-based (PALS) atlas of the cerebral cortex; left brain hemisphere. Sample: 970 healthy subjects.

Fig 3

Reproducibility of RSFC patterns from Broca (top panel) and Wernicke (bottom panel). Scatter plots across brain voxels showing the proportionality between the strength of the connectivity for each institution (horizontal axis) and the average connectivity across institutions (averaged across all 970 healthy subjects; vertical axis).

Fig 4

Statistical significance of correlations with time-varying signals in Broca’s (A) and Wernicke’s (B) areas, rendered on lateral and medial cortical views of the brain hemispheres across 970 healthy subjects. Line plots reflect the distribution across subjects of correlations with time-varying signals in Broca’s (RB; C) and Wernicke’s (RW; D) at specific regions-of interest (labels identify Brodman areas in the left, L, and right, R, hemispheres).

Fig 5

Statistical significance (color-coded t-score windows) for: RSFC (RB and RW) (A), differential (short-range – long-range) FCD (B), differential (Broca vs. Wernicke) RSFC (C), and LI for long-range FCD (left hemisphere versus right) (D) superimposed on lateral views of the cortical brain surface. Sample size: 970 (A and B) and 947 (C and D) healthy subjects. Contour lines are the boundaries of the network positively correlated with the seeds (solid line) and potential subdivisions of the network (dashed lines).

Anticorrelated activity

Several regions showed anticorrelated activity with Broca’s and Wernicke’s areas (Fig 5A, green-blue; Table 2, RB and RW), and this pattern was significantly stronger in visual (BAs 17–19), auditory (BAs 41, 42 and 22), somatosensory (BAs 1–3 and 5) and motor (BAs 4 and 6) cortices for the Wernicke’s seed than for the Broca’s seed (Table 2 and Fig 6).

Fig 6

Brain surface rendering of statistical significance maps showing the brain regions where the negative correlations were stronger for Wernicke’s seed than for Broca’s seed (t-test; 970 subjects). B: Brain surface rendering of Brodmann Areas (BA) classically associated with primary/association visual (BAs 17–19), auditory (BAs 41, 42 and 22), somatosensory (BAs 1–3 and 5) and motor (BAs 4 and 6) cortices.

ROI-voxel correlation patterns

The RSFC patterns from ROIs in Broca’s and Wernicke’s areas were similar to those from the corresponding seeds (SI Fig 2A). However, the average strength of the positive functional connectivity across voxels in the combined network was stronger for ROI-voxel correlation patterns than for seed-voxel correlation patterns (Broca: RROI = 0.145 ± 0.073, Rseed = 0.082 ± 0.039; Wernicke: RROI = 0.133 ± 0.074, Rseed = 0.082 ± 0.038). The average strength of the negative functional connectivity across voxels in the anticorrelated network was also stronger for ROI-voxel correlation patterns than for seed-voxel correlation patterns (Broca: RROI = – 0.050 ± 0.014, Rseed = – 0.032 ± 0.008; Wernicke: RROI = – 0.088 ± 0.013, Rseed = – 0.059 ± 0.007). Functional connectivity differences between ROI-and seed-voxel methods were statistically significant (SI Figs 2C and 2F) and reflect the higher signal-to-noise ratio due to noise averaging across the 125 voxels in the ROIs.

Modularity

The modularity-weighted algorithm independently identified four optimal partitions (modules or communities) for the language network (Fig 7A and 7B). Module #1 grouped Wernicke’s area with parietal, temporal and frontal cortex (bilateral) and the cerebellum; Module #2 grouped bilateral ROIs in visual, auditory, somatosensory and default-mode network (precuneus and cingulum); Module #3 grouped Broca’s area with pars traingularis (right hemisphere); and Module #4 grouped the thalamic and striatal ROIs. This modularity of the language network was highly consistent across institutions (Fig 7C).

Fig 7

Fig 7A: The averaged modular correlation matrix across 970 subjects consists of Pearson correlation coefficients, R, among 23 ROIs (Table 2) and identified four non-overlapping modules: 1-Wernicke’s module, yellow; 2-Anticorrelated module; blue; 3-Broca’s module, red; and 4-thalamic-striatal module. B: Potential spatial partitions of the language network. C: The modular correlation matrices were highly reproducible across the 22 research institutions.

Functional connectivity density patterns

The average distribution of the short-range FCD was maximal in posterior cingulate/ventral precuneus (BA 23/31; SI Fig 3). Other brain regions that included high short-range FCD included cuneus and middle occipital gyrus (BA 18), cingulate gyrus (BA 24), inferior parietal cortex (BA 40), as well as precentral (BA 6), inferior and middle frontal gyri (BAs 6 and 9), thalamus and cerebellum. The average distribution of the long-range FCD was maximal in the visual cortex, posterior cingulate and ventral precuneus (SI Fig 4). Angular gyrus, superior and inferior parietal, medial prefrontal (BA 10 and 11) and temporal cortices also had significant long-range FCD. A direct comparison of rescaled distributions revealed that the long-range FCD was higher in posterior parietal (BAs 3, 5, 7, 37, and 40), occipital (BAs 17, 18, 19), temporal (BAs 21 and 22) and frontal (BA 10) cortices and lower in dorsolateral prefrontal (BAs 6, 8, 9, 11, 44–47), inferior parietal (BA 40), temporal (BA 20), and limbic (BA, 23, 24, 32) cortices, as well as subcortical regions (putamen, caudate, thalamus, midbrain, and cerebellum) (PFWE < 0.001; t-test; Fig 5C and SI Fig 5) than the short-range FCD. Note that areas that showed higher long-range than short-range FCD also showed anticorrelated activity with that of the combined network (visual, auditory and posterior parietal cortices; Table 2).

Lateralization Index

The LI patterns in the combined network revealed that the long-range FCD was left-lateralized in Broca’s area, middle frontal gyrus, pars opercularis and orbitalis, posterior parietal cortex (angular gyrus), caudate and ventral thalamus, but contrary to our hypothesis it was right-lateralized in Wernicke’s area (BA 40) and in superior frontal (BA 8) and temporal (BA 22) cortices (PFWE < 0.05; Fig 5D and Table 2 and SI Fig 6).

Gender effects

Broca’s area exhibited stronger RSFC with left pars triangularis (P = 0.007; t-test; positive connectivity), cingulate and left superior temporal cortices (P < 0.05; negative connectivity) for females than for males. Other gender effects on RSFC strength in the language network were not statistically significant.

DISCUSSION

Here we document a striking reproducibility of RSFC patterns across 970 healthy subjects from 22 research institutions worldwide, despite demographic differences, variable MRI technology and acquisition parameters (Fig 5A). We also uncovered a pattern of anticorrelated activity in primary sensory and motor cortices and in precuneus with that in the combined network. The pars triangularis showed stronger (gender-related) coupling with the Broca’s seed than with the Wernicke’s seed (Fig 5B) and four distinct functional modules were identified within the language network (Fig 7). Functional connectivity density mapping revealed an overall predominance of short-range over long-range FCD in extended regions of the combined network (Fig 5C); only the posterior region of Wernicke’s area exhibited predominance of long-range connectivity. Furthermore, pronounced lateralization emerged in prefrontal and posterior parietal language regions (leftward) as well as auditory and medial and dorsolateral prefrontal cortices (rightward) for long-range (Fig 5D) but less so for short-range FCD.

RSFC patterns

The MRI signals in Broca’s and Wernicke’s areas were positively correlated with signals in an extended network of regions that included the inferior frontal network (pars triangularis, pars orbitalis and pars opercularis), middle frontal gyrus, inferior temporal and temporo-parietal areas (cited as supramarginal gyrus, planum temporale, Sylvian parieto-temporal, superior temporal gyrus, and inferior parietal cortex) as well as subcortical regions that included bilateral caudate and left putamen/globus pallidus, as well as ventral thalamic and subthalamic nuclei. Most of the cortical areas identified in this network are consistent with those reported on a recent review on language processing[2] that support their association with speech comprehension and articulation.

The RSFC analysis also revealed that the lateral orbitofrontal cortex (BAs 11 and 47) was functionally connected to the Broca’s and Wernicke’s seeds. This finding is consistent with the hypothesis that the orbitofrontal cortex is crucial for semantic processing[28], and patients with language disorders (semantic dementia) have shown cortical thinning in orbitofrontal cortex.[29] In addition, the RSFC analysis revealed several subcortical (bilateral caudate, left putamen/globus pallidus and left ventral thalamus/subthalamic nuclei) and right cerebellar regions that were functionally connected to Broca’s and Wernicke’s areas for which there is also evidence for their involvement in language processing. The caudate, for example, was implicated in monitoring and controlling lexical and language alternatives in language production tasks,[30] and its lesion in a polyglot patient was associated with involuntarily switching from one language to another.[31] Moreover recent fMRI studies have suggested that the caudate might predict how well a second language could be learned.[32] The putamen/globus pallidus were implicated in the temporal synchronization of cortically generated contextual and linguistic modules[33] and lesions in caudate and left putamen have been associated with aphasia.[34] The subthalamic nucleus was implicated in the establishment of the firing mode of thalamic nuclei underpinning the engagement of language-dedicated cortical components.[35] Finally there is growing evidence supporting the involvement of the right cerebellum in language processing[36] and language impairment in ASD.[37]

The variability in the strength of the RSFC in the language network allowed the differentiation between regions that were preferentially coupled to Wernicke’s area (inferior, lateral and ventral parietal, middle frontal and inferior temporal cortices, planum temporale and prefrontal regions BA 9, 10 and 46) from those preferentially coupled to Broca’s area (pars triangularis, orbitalis and opercularis, BA 6 and BA 37). These patterns are consistent with findings in non-fluent aphasia patients with degeneration in frontal or occipital cortices.[38]

The MRI signals in Broca’s and Wernicke’s areas were also negatively correlated with signals in auditory (BA 41 and BA 42), visual (BA 17 and 18) and precuneus (BA 7) cortices. These findings are consistent with the central role of anticorrelated functional networks in the functional organization of the brain[39] which cannot be attributed to signal processing artifacts.[40] This suggests a functional temporal segregation of auditory and visual perception (frequencies and loudness of speech sounds, seeing word forms) and that of speech comprehension and production. Note that superior temporal cortex damage in stroke patients reduced the ability to comprehend spoken sentences[41] and that the visual cortex has been its implicated in the identification of visual letter strings.[42]

The stronger anticorrelated activity in visual and auditory cortices for the Wernicke’s seed than for the Broca’s seed is consistent with language processing, which relies on visual and auditory cortices much more than the motor component of speech. Indeed anticorrelated activity with the auditory cortex (BA 42) was only significant for Wernicke’s seed suggesting that it may uniquely relate to the role of the auditory integration cortex (BA 42) in listening to spoken words.

Functional connectivity density

Broca’s network demonstrated higher proportion of short-range than long-range FCD. The predominant reliance of Broca’s network on short-range connectivity is consistent with the expansion in size and complexity of cortical association areas during the evolution of the human brain[43], and with the novelty-routinization theory that postulates that routine information processing in modality-specific association cortex favors local connectivity between adjacent cortical regions.[11] Since the length of the axons has been associated with their vulnerability to degenerative processes, as is the case for amyotrophic lateral sclerosis and spastic paraplegias[44], and postmortem studies have reported decreases in myelinated fiber length with age in the human brain[45] (10% decrease per decade), the low proportion of long-range FCD in Broca’s network might increase its resilience to neurodegenerative processes. On the other hand the high proportion of long-range FCD hubs in posterior regions of Wernicke’s area (angular gyrus) might provide critical gateways for semantic word processing[46] and could underlie the decline in semantic abilities with age and its disruption in early Alzheimer’s Disease.[47] The predominant long-range FCD in angular gyrus is consistent with its role in modulating information flow from the left dorsal occipital gyrus that processes words recognition into the left supramarginal gyrus[48] which could explain why disconnection of the angular gyrus is associated with dyslexia[49] and decreased connectivity with poor reading abilities.[50] The angular gyrus has also been implicated in modulating information from prefrontal cortices to facilitate speech comprehension as a function of semantic context.[51]

Laterality

The analysis of the laterality patterns revealed a leftward lateralization for the long-range connectivity in Broca’s area and in posterior Wernicke’s (angular gyrus), which is consistent with previous RSFC studies[52] and supports the lateralization of language to the left hemisphere. This lateralization might reflect the optimization of neural processing involved in the control of complex movements for speech production that might not be efficiently performed simultaneously by both hemispheres. Abnormalities in the leftward asymmetry of the language areas have been implicated in developmental dyslexia,[7] schizophrenia[8] and ASD.[9] However, we also document an unexpected rightward lateralization of the anterior Wernicke’s region for long-range connectivity that suggests a predominant involvement of the right hemisphere in language comprehension processed through the supramarginal gyrus. Indeed there is increasing evidence that the right hemisphere is involved with processing figurative language (ie metaphors, idiomatic expressions) in contrast to processing literal interpretation[53] and of the right supramarginal gyrus involvement in phonological decisions.[54]

Network Modularity

Recent studies have demonstrated the test-retest reliability of the modularity of the resting-state networks.[55] The functional parcellation analysis in the present study identified four interacting modules that reflect the fundamental organization of the language network: Module #1 grouped regions strongly coupled to Wernicke’s area and are associated with speech comprehension;[2] Module #2 grouped regions anticorrelated with Wernicke’s area; Module #3 only included pars triangularis, which is traditionally associated with speech production;[2] and Module #4 contained the striatum and anterior thalamus. The high reproducibility of these partitions across research institutions further supports the segregation of the language network.

Females had stronger Broca’s area connectivity with pars triangularis (module 3), cingulum and left auditory cortex (module 2) than males. This is also consistent with the females’ higher short-range connectivity in heteromodal cortices[56] and higher long-range connectivity in the default-mode network,[24] as well as with previously reported gender effects on white matter[57,58] and with the gender differences in the lateralization of the functional connectivity in inferior frontal and superior temporal cortices.[27] These differences might reflect genetic factors[59] as well as testosterone levels during fetal development (higher for males than for females) that could influence neural connectivity by averting programmed cell death during neural development.[60] However, men and women have different communication skills and styles.[61-64] Therefore, the higher connectivity of the speech production network for females than for males could also reflect differences in socialization between genders.

Study limitations

Some of the variability of the RSFC patterns could reflect differences in respiratory volume per unit time (RVT), which causes BOLD modulation patterns.[65] We were unable to control for this effect because RVT data was not available in the “1000 Functional Connectomes” dataset. Future resting-state studies on gFCD should control for RVT. Since a large fraction of the signal fluctuations are shared throughout the brain,[40] the global signal is another potential source of variability for studies on functional connectivity. Global signal fluctuations were removed during preprocessing steps in the present study in order to minimize their contribution to the variability of the gFCD patterns. However, recent studies suggests that much of the global signal reflects neuronal activity.[66] Thus the removal of global signal fluctuation could have decreased the amplitude of the long-range FCD in this study. Single-voxel RSFC is prone to arbitrary seed location, there is no clear delineation of the Broca’s and Wernicke’s areas and their precise localization might vary across individuals. The complementary ROI-voxel correlation analyses assessed the reliability of the results in 125 voxels that surround the seeds. The striking similarity of the ROI-voxel and seed-voxel correlation patterns together with the increased strength of the RSFC for ROI-voxel method compared to seed-voxel method demonstrates that the results of this study were robust under variable seed location.

Thus, here we document an extended network with positive correlated activity that includes Broca’s and Wernicke’s areas as well as frontal, parietal and temporal cortices, caudate, putamen and subthalamic nuclei that is consistent with reports of language impairment in patients with lesions in these brain regions. We also show anticorrelated activity in temporal and occipital cortices that have been implicated in auditory and visual language processing, respectively, and that suggests temporal segregation between these two components.