Participants
Neuroimaging data was acquired as part of the Baby Connectome Project53. A total of 319 participants spanning 694 sessions were available for analysis. Some infants were scanned multiple times across sessions with age. In order to adhere to assumptions regarding statistical independence of data in our models we included data from only one session per participant, choosing the session with the highest whole-brain temporal signal to noise ratio (tSNR) of the resting-state data. Fifty-four participants were excluded due to poor data quality quantified by excessive motion (AFNI’s @1dDiffmag values > 0.3 mm/TR) or low whole-brain tSNR (< 30). An additional 10 participants were excluded due to failed co-registration during data preprocessing that could not be remedied without extensive manual alteration. Finally, as there were few participants older than 25 months (n = 38 spanning 26–69 months of age), these participants were also excluded. The final sample consisted of 212 infants between the ages of < 1 month and 25 months. We divided participants into four age bins spanning 6–7 months in age: 0–6 months (n = 50), 7–12 months (n = 59), 13–18 months (n = 58), and 19–25 months (n = 45). The decision to bin infants into discrete categorical age ranges was made to allow the exploration of the relationship between age and connectivity without making assumptions about the shape of the association (i.e. linear or nonlinear, as required when modelling continuous data). The age bin cutoffs were constructed to allow a comparable number of participants in each age bin, with roughly equal age ranges. A histogram of the specific ages of infants in each age bin can be viewed in Supplementary Fig. 1.
Procedure and data acquisition
Detailed information regarding the imaging protocol is described by Howell and colleagues (2019). All imaging was acquired on 3T Siemens Prisma MRI scanners at the Center for Magnetic Resonance Research (CMRR) at the University of Minnesota and the Biomedical Research Imaging Center (BRIC) at the University of North Carolina at Chapel Hill, using a Siemens 32 channel head coil. We here utilize the T1w (MPRAGE, 0.8mm isotropic voxels, TE = 2.24, TR = 2400/1060), T2w (3D variable flip angle turbo spin-echo sequence, 0.8mm isotropic voxels, TE = 564, TR = 3200), and two resting-state T2*-weighted functional scans (gradient-echo EPI acquisition: Multiband 8, TR = 800ms, TE = 37ms, flip angle = 52 degrees, FOV = 208mm, 104x91 matrix, 2mm isotropic voxels) with anterior-to-posterior and posterior-to-anterior encoding (420 TRs, 5 min 36 sec for each encoding; total EPI duration: 840 TRs, 11 min 12 sec) per participant.
Preprocessing
Data processing was conducted with Freesurfer (RRID: SCR_001847), infant Freesurfer54, and AFNI (RRID:SCR_005927). The T1w anatomical data were processed using either Freesurfer (for participants > 10 months of age) or infant Freesurfer (for participants < = 10 months of age) for tissue segmentation. Different versions of Freesurfer were required for different age ranges, as the grey matter white matter contrast is poor in infants less than a year old, often resulting in failed convergence for adult Freesurfer. Anatomical scans and segmented grey matter, white matter, and CSF volumes were linearly aligned to the EPI data using AFNI’s 3dAllineate. In infants > 10 months, warping the T1w data to the EPI data produced a good coregistration. For infants < = 10 months, the T1w scan was first aligned to the T2w scan (which had better grey matter/white matter contrast). The T2w scan was then co-registered with the EPI data, and the warp was applied to the T1w data to align it with the EPI data. The functional data were despiked (3dDespike), adjusted for slice-time acquisition (3dTshift), warped to the mid-point of the P-A and A-P encodings and volume registered using AFNI’s unWarpEPI, and frame-to-frame motion was estimated. The data were then smoothed with a 4mm FWHM kernel (3dmerge) and rescaled to % signal change. We then applied a modified form of the ANATICOR (Jo et al., 2010) denoising approach with a nuisance regression (3dTfitter) containing: 6 motion parameters, local average of WM signal (averaged within a 15 mm radius sphere), average ventricle signal, as well as the top 3 principal components of the time series from all white matter and ventricle voxels (modified aCompCor, Behzadi et al., 2007; Stoddard et al., 2016). All regressors were detrended with a fourth-order polynomial before denoising (3dDetrend), and the same detrending was applied during nuisance regression to the voxel time series. The residual time series were concatenated across the A-P and P-A resting-state runs. Subject-specific anatomical and functional data were first registered to an age-specific MNI template, and subsequently, to the median age T1w MNI template (11–14 months) using nonlinear warps (3dQwarp).
Region of interest definition
The hippocampus was isolated from the Freesurfer or infantFreesurfer subcortical segmentations for each participant and manually corrected where necessary. The hippocampus was then manually segmented into anterior and posterior segments at the uncal notch55 to create anterior and posterior masks for each hemisphere for every participant. The uncal notch was identifiable in even the youngest infants. The resulting masks were co-registered to the group analysis space (median age template) using the warps created during preprocessing.
Functional connectivity computation
We conducted a seed-to-voxel functional connectivity analysis using the left and right anterior and posterior hippocampal masks as seeds. Specifically, the average time-course across voxels within each hippocampal region was calculated using AFNI’s 3dmaskave. Next, the Pearson’s correlation between each hippocampal ROI time-course and the time-course of each voxel in the brain was calculated using 3dTcorr1D. The resulting correlation maps (one for each ROI, reflecting whole-brain functional connectivity), were Fisher-z transformed to normalize the data for group analyses.
Group-level contrasts
We ran a linear model using AFNI’s 3dMVM to test for group differences in connectivity across the age bins. Specifically, connectivity was predicted as a function of long-axis (anterior/posterior), age bin (0-6mo/7-12mo/13-18mo/19-25mo), and hemisphere (left/right). We also included whole brain mean tSNR and mean motion (estimated with AFNI’s 3dDiffmag) as regressors of no interest to adjust effects for variability in tSNR and motion across participants. tSNR is useful for inclusion in these analyses as it is a composite of all sources of nuisance variability that impact the BOLD signal and functional connectivity, including unmodeled effects of motion and respiration. We were particularly interested in the interaction between long-axis compartments and age bin. To directly explore the degree of differential anterior-posterior connectivity at each age, we additionally specified a-priori pairwise contrasts comparing anterior > posterior hippocampal connectivity within each age bin. Finally, we contrasted hippocampal connectivity across successive age bins (7-12mo > 0-6mo, 13-18mo > 7-12mo, 19-25mo > 13-18mo) to better understand the main effect of age bin on hippocampal-cortical connectivity. Voxels surviving a threshold FDR-corrected to q < 0.01 corrected across all voxels in the brain were considered statistically significant. We chose this method of correction because it is robust to inflated false positives due to multiple comparisons. To further mitigate false positives, we additionally conservatively implemented a cluster-forming threshold to remove very small clusters of less than 10 contiguous voxels.
Hierarchical clustering analysis
Having identified several neocortical clusters of voxels where connectivity was changing along the long-axis as a function of age bin (quantified by a significant long-axis x age bin interaction), we next sought to describe the patterns driving the interaction for each region identified. To facilitate interpretation, we clustered together neocortical regional clusters evincing a similar direction of interaction to form what we term here, “superclusters”. As we did not know the number of superclusters a-priori, we used a hierarchical clustering approach. We extracted the average connectivity across participants and hemispheres between the anterior and posterior hippocampi and each neocortical cluster where there was a significant interaction for each age bin. We subtracted posterior from anterior hippocampal connectivity to get a measure of the difference in connectivity along the long axis for each neocortical cluster for each age bin. We calculated the difference in connectivity along the long-axis because we were interested specifically in the interaction between long-axis and age bin and did not want the clustering solution to be driven by the main effects of long-axis (as in the case of including anterior and posterior connectivity values separately). The result was a 44 neocortical clusters x 4 age bins matrix of anterior-posterior differential connectivity. We then transformed this matrix into a distance matrix by quantifying the Euclidean distance between each of these observations. Finally, we clustered the distance matrix into superclusters using Ward’s method.
The agglomerative clustering process starts with one cluster for each observation (in this case, the connectivity values between the hippocampus and each neocortical cluster) and combines observations with the smallest distance between them to form new clusters (i.e. superclusters). Ward’s method sums the squared Euclidean distances between each member of a cluster and the cluster centroid, joining clusters which increase the sum of within-cluster square errors the least56. A new distance is calculated between new clusters, and the process is repeated until there is only one cluster with all observations in it. We used the agglomeration schedule to determine the optimal number of superclusters, choosing the point in the agglomerative process where the distance between conglomerated clusters jumps, indicating the combination of distant clusters that are different enough that they likely shouldn’t be combined.
Statistical Analyses
Statistical analyses of connectivity or each supercluster were computed in R studio version 2023.12.1.40257( http://www.posit.co/). We used a linear mixed effects model to test connectivity as a function of supercluster, long-axis, and age bin with a random intercept for each participant using the nlme package in R (version 3.1.158; https://CRAN.R-project.org/package = nlme)58. As model residuals were originally heteroskedastic, we specified the variance structure to allow the variance to vary across levels of the heteroskedastic fixed effects, as identified using Levene’s test for unequal variance as part of the car package in R59 (version 3.1.0; https://cran.r-project.org/web/packages/car/index.html). Denominator degrees of freedom and p-values were estimated using the containment method as implemented in nlme. Post-hoc comparisons were interrogated using pairwise t-tests on the estimated marginal means from the omnibus model, using the emmeans package (version 1.7.5, https://cran.r-project.org/web/packages/emmeans/index.html), and were FDR corrected to control for multiple comparisons. To confirm that results from these statistical analyses were not driven by subject-level bias, we repeated our statistical tests on weighted averages computed in a leave-one-out validation approach, as described in the Supplementary Methods.
Cluster validation
We took several steps to validate the hierarchical clustering solution. First, we visually verified that the superclusters were distinct by plotting the first two dimensions of a principal components analysis. Second, we plotted differential anterior-posterior hippocampal connectivity over age bins for each neocortical cluster comprising the superclusters, to visually check that neocortical clusters were being grouped together sensibly based on the interaction between long-axis and age bin. Third, we calculated silhouette statistics to assess overall supercluster quality, and to quantify the degree to which neocortical clusters within each supercluster fit the supercluster they were assigned to versus other superclusters. In order to perform statistics to compare the patterns of each supercluster, it was important to verify that the comparisons would not be biased by subject-level noise as a result of performing the selection and comparisons on the same data (i.e. double-dipping; e.g. Kriegeskorte et al., 2009). To evaluate this, the hierarchical clustering was repeated while leaving each subject out once (212 times). We then tabulated which clusters were grouped together into superclusters across all 212 iterations, which by definition does not depend on subject-level noise. We then constructed “bias-free” estimates of the supercluster patterns using a weighted average of only those pairs of clusters that were grouped together 100% of the time. The bias-free, cross-validated patterns could then be compared to the original patterns based on the entire group of subjects.