The current study complies with ethical regulations set by The Independent Research Ethics Committee at the Medical Faculty of the Heinrich-Heine-University of Duesseldorf.Participants and study designOur analyses were conducted on the publicly available data of healthy young adults from the Human Connectome Project (HCP) S1200 release (http://www.humanconnectome.org/)33. We selected subjects with available functional, T1, and T2 neuroimaging data, resulting in a final sample of 1000 individuals (536 females) with a mean age of 28.73 ± 3.71 years, used for all analyses. The sample included 284 monozygotic twins (MZ), 184 dizygotic twins (DZ), 443 non-twin siblings, and 89 unrelated individuals, and the sociodemographic breakdown of the participants is additionally reported in Supplementary Table 3. Subjects were all born in Missouri but recruited in an attempt to broadly reflect the racial and ethnic composition of the United States population. Recruitment efforts aimed to yield a subject pool capturing a wide range of variability—in socioeconomic and behavioral terms—in order to be representative of the general healthy population. The term “healthy” was thus broadly defined. Individuals with documented neurodevelopmental and psychiatric disorders, or reporting physiological illnesses such as high blood pressure or diabetes were excluded from the HCP study recruitment protocol, but not individuals who reported smoking, being overweight, or a history of recreational drug use or heavy drinking (if they had not experienced severe symptoms). Informed consent was obtained for all study subjects and recruitment procedures were approved by the Washington University institutional review board. More detailed information about the HCP study design and recruitment procedure is available elsewhere33,72. Despite the use of the term gender in the HCP Data Dictionary, we use the term sex in this article given that the HCP study collected self-reported information on biological sex instead of gender identification, as reported elsewhere73. We have not used genetic information to verify the self-reported sex.Structural MRI acquisition and preprocessingThe HCP’s MRI data was acquired on a customized 3T Siemens Skyra ConnectomeScanner with a 32-channel head coil at Washington University across four scanning sessions held over two days. Structural MRI images were acquired on the same day via high resolution T1-weighted (T1w) and T2-weighted (T2w) sequences. Two separate T1w images were acquired and averaged, with identical scanning parameters using a 3D MPRAGE sequence (0.7 mm isovoxels, FOV = 224 mm, matrix = 320 × 320 mm, 256 sagittal slices; TR = 2400 ms, TE = 2.14 ms, TI = 1000 ms, flip angle = 8°, BW = 210 Hz per pixel, ES = 7.6 ms). Two separate T2w images were acquired and averaged, with identical scanning parameters using a variable flip angle turbo spin-echo (3D T2-SPACE) sequence, with the same isotropic resolution, matrix, FOV, and slices as for the T1w sequence (TR = 3200 ms, TE = 565 ms, BW = 744 Hz per pixel, total turbo factor = 314). The preprocessing steps included co-registering the T1w and T2w images, bias field (B1) correction, registration to MNI space, segmentation, and surface reconstruction. See refs. 33,72,74 for more detail on the HCP’s MRI protocols and the FreeSurfer segmentation pipeline.Functional MRI (fMRI) acquisition and preprocessingThe HCP’s fMRI data was collected after the structural sequences and following the HCP’s minimal processing pipeline, as described above. A total of 1 h of resting-state functional data was collected across four identical 15 min scanning sessions, equally split over two days (LR1, RL1, LR2, RL2), with a gradient echo EPI sequence at a resolution of 2 mm isotropic (FOV = 208 × 180 mm, matrix = 104 × 90 mm, 72 slices covering the whole brain, TR = 720 ms, TE = 33 ms, multiband factor of 8, FA = 52°). The multimodal surface matching algorithm (MSMAll) was used to co-register the data to the HCP template 32k_LR surface space, consisting of 32,492 nodes per hemisphere (59,412 nodes excluding the medial wall). A more detailed description of the resting state fMRI data acquisition and analysis protocol is available elsewhere74,75.Functional connectivity (FC) and the sensory-association (S-A) axis of functional organizationThroughout this work, we used the Schaefer 400 parcellation (clustered into 7 networks: visual, somatomotor, dorsal attention, ventral attention, limbic, frontoparietal, DMN40). This widely used functionally-derived parcellation scheme was originally obtained via a gradient-weighted Markov Random Field model integrating local gradient and global similarity approaches76. The vertex-wise functional timeseries were therefore averaged within the Schaefer 400 cortical regions. FC matrices (400 × 400) were then computed at the individual level—per scanning session—by correlating cortical timeseries in a pairwise manner using the Pearson product moment. We normalized the correlation coefficients using Fisher’s z-transformation. Final FC matrices were obtained by averaging each subject’s matrices across their four scanning sessions. From these FC matrices and for each subject, we computed the S-A axis of functional organization, as described below.We conducted data reduction on the FC matrices to yield macroscale gradients of functional organization10. For this, we used diffusion map embedding, a nonlinear manifold learning algorithm that reduces complex, high-dimensional structures of data (in our case affinity matrices) to low-dimensional representations combining geometry with the probability distribution of data points34. Thus, cortical regions that are strongly interconnected (i.e., whose timeseries show high correlations) are represented closer together in the resulting low dimensional manifold of FC data, whereas parcels with low covariance are represented farther apart, as indexed by the cortical regions’ gradient loadings. To this end, we used the BrainSpace Python toolbox35 to compute ten gradients with the following parameters: 90% threshold (i.e., only considering the top 10% row-wise z values of FC matrices, representing each seed region’s top 10% maximally functionally connected regions), α = 0.5 (α controls whether the geometry of the set is reflected in the low-dimensional embedding— i.e., the influence of the sampling points density on the manifold, where α = 0 (maximal influence) and α = 1 (no influence)), and t = 0 (t controls the scale of eigenvalues). These parameters were selected for consistency with previous studies and represent choices that are recommended to retain global relations between datapoints in the embedded space whilst being relatively robust to noise10,42. To confirm the robustness of the S-A axis computed at the 90% threshold, we further show with a sensitivity analysis that the mean S-A axis computed at the 90% threshold shows high correlations with mean S-A axes computed at different thresholds (from 10% to 90%, in steps of 10%), with r values ranging from 0.84 to 0.93 (see Supplementary Fig. 5). In order to increase comparability for further between-subject analyses, Procrustes alignment was used to align individual gradients to mean gradients. Mean gradients were computed by applying diffusion map embedding—with the same parameters listed above—to the mean FC matrix (i.e., FC matrices averaged across all subjects). The computation of these FC gradients was carried out independently per hemisphere (i.e., considering the top 10% row-wise z values of only half of the FC matrices, shaped 200 × 200) and the gradient loadings resulting from each hemisphere were subsequently concatenated. This decision was made for consistency and comparability reasons within our study, so that the top 10% functional connections selected for data reduction corresponded to those considered in the calculation of the mean geodesic distance of connectivity profiles—which were only computed per hemisphere– as described further below). We verified and confirmed the stability FC gradients when computing them per hemisphere versus at the whole brain level, as shown by the spatial correlation of mean gradient loadings (r = 0.98, pspin < 0.001). Finally, we took the well-replicated principal gradient explaining the most variance in the data and spanning from visual to DMN regions10, which we labeled the S-A axis and used to represent functional organization for subsequent analyses.We also computed, for each subject, mean FC strength at the parcel level in a seed-wise fashion, by averaging the row-wise z values of each seed region’s top 10% maximally functionally connected regions—again per hemisphere—and subsequently concatenated the hemispheric mean FC strength values to reconstruct whole brain data.Cortical microstructure and microstructural profile covariance (MPC)Microstructural properties—including myelin and cellular characteristics—show depth-dependent variation along cortical columns, as reported by histology42,77,78 as well as in vivo and post mortem neuroimaging37,41,42,78, which illustrate a cortical hierarchy11. Similar to previous work37, we quantified cortical microstructure, or “microstructural profile intensity” (MPI), using the myelin-sensitive MRI contrast obtained from the T1w/T2w ratio from the HCP minimal processing pipeline described above74 (a reliability check is reported in the Supplementary Fig. 6). The T1w/T2w ratio uses the T2w image to correct for inhomogeneities in the T1w image50. Then, we followed the previously described protocol37,41,42 to compute our measurement of MPC, which reflects the variation of MPI, across cortical depths. In short, we generated 14 equivolumetric surfaces within the inner and outer cortical surfaces, then excluded the inner- and outer-most surfaces, thus remaining with 12 surfaces representing cortical layers. Surface generation was based on a model compensating for cortical folding by altering the pairwise Euclidean distance (ρ) of intracortical surfaces throughout the cortex and thus preserving fractional volume between the surfaces. For each surface, ρ was calculated as defined in Eq. 1.$${{{\rm{\rho }}}}=\frac{1}{{A}_{{out}}-{A}_{{in}}}\cdot \left({-A}_{{in}}+\sqrt{\alpha {A}_{{out}}^{2}+\left(1-\alpha \right){A}_{{in}}^{2}}\right)$$
(1)
for which α denotes a fraction of the total volume of the segment that the surface accounts for, while Aout and Ain respectively denote the surface areas of the outer and inner cortical surfaces.Across the whole cortex and from the outer to the inner surfaces, we systematically sampled MPI values layer-wise for each of the 64,984 vertices of the HCP template 32k_LR surface space, which we then averaged within each of the 400 Schaefer cortical areas, per layer. Following a previously described protocol41, we constructed subject level 400 × 400 matrices using pairwise Pearson partial correlation on the MPI profiles of cortical parcels (i.e., correlating the MPI values across 12 layers between parcels), controlling for overall mean cortical MPI, followed by log transformation. We then used these matrices to compute MPC gradients by following the same procedure and using the same toolbox and parameters as for computing the FC gradients34,35, as described above and previously done37,41,42. For consistency and comparability with the computation of FC gradients and mean geodesic distance of connectivity profiles, we also computed MPC gradients independently per hemisphere and subsequently concatenated the gradient loadings resulting from each hemisphere. We verified and confirmed the stability MPC gradients when computing them per hemisphere versus at the whole brain level, as shown by the spatial correlation of mean gradient loadings (r = 0.99, pspin < 0.001). We also selected the principal gradient of MPC explaining the most variance in the data and spanning from sensory to paralimbic regions, which we labeled the MPC axis and used to represent microstructural organization in subsequent analyses.Measures of brain sizeIn our analyses we included different measures of brain size typically used in the literature, including intracranial volume (ICV), total cortical volume (TCV), and total surface area. For ICV, we used the FreeSurfer output measure IntraCranialVol, which is an estimate of ICV based on the Talairach transform. We computed our own measure of TCV by summing the volumes of the TotCort_GM_Vol and Tot_WM_Vol FreeSurfer output measures, which we considered relevant to our study’s focus on cortical functional organization (thus excluding the volumes of subcortical structures). We computed total surface area by using the FreeSurfer mri_surf2surf tool to resample cortical gray matter surface for each subject.Geodesic distance of connectivity profilesGeodesic distances, representing the shortest distance between two vertices along the folded cortical mantle’s curvature, were computed using the Micapipe toolbox79, and following the previously described protocol80. In short, geodesic distance matrices were computed for each subject along their native cortical midsurface. The first step consisted in defining a centroid vertex for each cortical parcel, identified as the vertex having the shortest summed Euclidean distance from all other vertices within the parcel. Then, Dijkstra’s algorithm81 was used to compute geodesic distances between the centroid vertices and all other vertices on the on the native midusrface mesh. The vertex-wise geodesic distance values were then averaged within each parcel to form the geodesic distance matrices. From these individual matrices, we finally averaged—parcel-wise—the geodesic distance values of each seed cortical region’s top 10% maximally functionally connected regions per hemisphere, thus obtaining for each subject the mean geodesic distance of functional connectivity profiles by region.Statistical analysisGiven that the HCP sample includes different levels of kinship, we used linear mixed effects models (LMMs) to account for sibling status (MZ, DZ, non-twin siblings) and family relatedness. In fact, all LMMs mentioned in this work consistently included sex, age, and total surface area as covariates (unless otherwise mentioned), and controlled for random nested effects of family relatedness and sibling status. In addition, effects on cortical data obtained via LMMs underwent false discovery rate (FDR) correction (q < 0.05), thus correcting for multiple comparisons across the 400 Schaefer cortical regions. Throughout this work, we also tested for associations in brain-wide patterns displayed in the form of cortical maps, for which we used Spearman-rank correlation followed by spin-permutation tests to control for spatial autocorrelation82. All statistical tests were two-sided.After computing the S-A axis of functional cortical organization, we tested for sex differences in the S-A axis loadings with an LMM. Then, we investigated which measure of brain size (out of ICV, TCV, and total surface area) had the largest effects on the S-A axis loadings using separate LMMs (respectively only including ICV, TCV, or total surface area as a covariate, in addition to sex, age and the random nested effect of family relatedness and sibling status). The reason underlying our decision to systematically include total surface area as a covariate in all our LMMs (as the measure of brain size) is that it showed the most widespread effects on the S-A axis loadings out of the three tested measures. Then, we investigated associations between the S-A axis and cortical morphometry, namely the MPC axis and the mean geodesic distance of connectivity profiles, using both LMMs and Spearman-rank spatial correlations of cortical maps.To probe whether sex differences in cortical morphometry may be associated with sex differences in the S-A axis, we tested whether sex differences in the S-A axis loadings were moderated by total surface area by modeling an additional interaction term of sex by total surface area on the S-A axis loadings within the original LMM. We also tested for sex differences in the MPC axis and in the mean geodesic distance of connectivity profiles, and conducted Spearman-rank correlations of cortical β-maps for the sex contrast in the S-A axis and in the morphometric measures. We also modeled two additional interaction terms within the original LMMs of sex by MPC axis loadings and sex by mean geodesic distance to show their effects on the S-A axis loadings. Finally, we conducted sensitivity analyses to test for sex effects on the S-A axis yielded by an LMM including all morphometric measures as covariates (i.e., including the MPC axis and the mean geodesic distance of connectivity profiles, in addition to total surface area), as well as an LMM not including any morphometric measures as covariates (i.e., also excluding total surface area). We then tested the similarity of both these sex effects with the original sex effects on the S-A axis with a Spearman-rank spatial correlation of the cortical β-maps.In order to probe the potential intrinsic functional underpinnings of sex differences in the S-A axis, we tested for sex differences in FC strength (also with an LMM), as well as sex differences in FC profiles, i.e., the presence of sex differences in the top 10% maximally functionally connected regions used to compute the S-A axis. To this end, we built 400 × 400 binary matrices at the subject level—based on the subjects’ individual FC matrix z values—in which we marked in a seed-wise fashion (along the matrix rows) whether the given cortical region (along the matrix column) belongs to the given seed’s top 10% maximally functionally connected regions, where 1 indicated that the parcel belongs to the seed’s top 10% maximally functionally connected regions and 0 indicated that the parcel does not belong to the seed’s top 10% maximally functionally connected regions. We then summed the binary matrices separately within sexes in order to fill 160,000 contingency matrices—one for each cell (i.e., functional connection) of the 400 × 400 FC matrix—denoting the number of males and females for which a given cortical region belongs or does not belong to the seed region’s top 10% maximally functionally connected regions (see Table 1 for a visual representation of the contingency matrix structure).Table 1 Contingency matrix structure for the computation of sex differences in functional connectivity profilesWe then conducted the Chi-square (χ2) test of independence (degrees of freedom = 1) on each contingency table to test for sex differences in the odds of each parcel of belonging to the top 10% maximally functionally connected regions of each seed region. Given the large number of tests conducted here (400 × 400 = 160,000), we controlled for multiple comparisons using FDR correction. We quantified the size of these sex effects with the odds ratio (OR), calculated as defined in Eq. 2.$${{{\rm{OR}}}}=\frac{{{{\rm{Cm}}}}/{{{\rm{NCm}}}}}{{{{\rm{Cf}}}}/{{{\rm{NCf}}}}}$$
(2)
where OR > 1 indicates greater male odds—and OR < 1 indicates greater female odds—of a given region of belonging to a given seed’s top 10% maximally functionally connected regions.We also tested for sex differences in network topology, i.e., how nodes are physically organized in networks and how networks are physically organized along the S-A axis. For this, we computed two measures of network dispersion: between- and within-network dispersion. Between-network dispersion is defined as the Euclidean distance between a pair of network centroids, where a higher value indicates that networks are more segregated from one another along the S-A axis. Within-network dispersion is defined as the sum squared Euclidean distance of network nodes (i.e., S-A axis regional loadings) to the network centroid, where a higher value indicates wider distribution and segregation of a given network’s nodes along the S-A axis. At the individual level, we thus computed between-network dispersion between all networks in a pairwise fashion (21 pairs), and within-network dispersion for all 7 networks, by defining network centroids as the median of the S-A axis loadings of all parcels belonging to a given network, following a previously described method43. Then, we computed sex differences in each of the 21 between-network dispersion metrics and 7 within-network dispersion metrics using LMMs. For each model, we computed a null distribution of β coefficients for sex differences using 1000 spherical rotations of the Schaefer parcellation scheme in order to shuffle the network labels82, against which we computed our p-value to determine statistical significance. We then assed pspinvalues against Bonferroni-corrected two-tailed α-levels of 0.001 (0.025/21) and 0.004 (0.025/7) for between-network and within-network dispersion sex contrasts, respectively.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Sex differences in functional cortical organization reflect differences in network topology rather than cortical morphometry
