The positive–negative mode link between brain connectivity, demographics and behaviour: a pre-registered replication of Smith et al. (2015)

In mental health research, it has proven difficult to find measures of brain function that provide reliable indicators of mental health and well-being, including susceptibility to mental health disorders. Recently, a family of data-driven analyses have provided such reliable measures when applied to large, population-level datasets. In the current pre-registered replication study, we show that the canonical correlation analysis (CCA) methods previously developed using resting-state magnetic resonance imaging functional connectivity and subject measures (SMs) of cognition and behaviour from healthy adults are also effective in measuring well-being (a ‘positive–negative axis') in an independent developmental dataset. Our replication was successful in two out of three of our pre-registered criteria, such that a primary CCA mode's weights displayed a significant positive relationship and explained a significant amount of variance in both functional connectivity and SMs. The only criterion that was not successful was that compared to other modes the magnitude of variance explained by the primary CCA mode was smaller than predicted, a result that could indicate a developmental trajectory of a primary mode. This replication establishes a signature neurotypical relationship between connectivity and phenotype, opening new avenues of research in neuroscience with clear clinical applications.


Introduction 1.Background
Understanding how brain functional connectivity is linked to human behavior is a major goal of neuroscience for its potential to elucidate the mechanisms of the brain that lead to human cognition [1].By modelling functional magnetic resonance imaging (fMRI) data as a region-to-region functional network, researchers can leverage mathematical tools from the study of complex systems (e.g., graph theory) to probe these high-dimensional relationships [2][3][4][5].Smith et al. (2015) In their landmark 2015 study, Smith et al. investigated the high-dimensional relationship between functional connectivity and behavioral and phenotypic measures using data on young adults aged 22-35 from the Human Connectome Project (HCP) [6].To study the relationship between these subject measures (SMs) and functional connectomes, the authors calculated a 200-dimension group-ICA functional parcellation for the 461 HCP subjects in their study.From this parcellation, they derived subject-level functional connectomes, which are symmetric matrices of the edge weights between all 200 nodes.Of the SM data available in the HCP dataset, Smith et al. chose 158 measures that were quantitatively and qualitatively appropriate for analysis.

Original Study by
Using canonical correlation analysis (CCA), which can simultaneously consider two variable sets from different modalities to uncover essential hidden and high-dimensional associations between the sets [7], the authors discovered that one statistically significant CCA mode (the primary mode) explained more of the observed covariance between functional connectivity and subject measures (SMs) compared to the other CCA modes.The authors found a strong positive correlation (r=0.87,p<10 -5 ) between the primary CCA mode connectome weights and SM weights.From the primary CCA mode, they derived a 'positive-negative axis' that quantified the correlation of various SMs to the CCA mode, linking lifestyle, demographic, and psychometric SMs to each other and to a specific pattern of brain connectivity [8].SMs commonly considered to be positive qualities (e.g., higher income, greater education level, and high performance on cognitive tests) were positively correlated with the primary CCA mode, and SMs commonly considered to be negative qualities (e.g., substance use, rule-breaking behavior, and anger) were negatively correlated with the CCA mode.
A hierarchical clustering analysis of the 200 connectome nodes revealed clusters of nodes in four regions: one cluster in the sensory, motor, insula, and dorsal attention regions; and three clusters covering the default mode network, subcortical, and cerebellar regions [8].The brain areas that contributed most strongly to variations in connectivity were similar to those associated with the default mode network (DMN), suggesting that functional connectivity within (and, to some degree, between) the DMN may be important for higher-level cognition and behavior [7].
Finally, to test the predictive value of their model, the authors conducted a train-test split in which 80% of subject data (randomly selected, respecting family relationships) was used to train a CCA model, which was then tested on the remaining 20% of the data.This process was repeated 10 times, with statistical significance estimated via permutation testing (1,000 permutations) of family structures in each iteration.This train-test process yielded a statistically significant mean correlation of r=0.25 between the primary CCA mode connectome and SM weights in the testing dataset, strongly supporting the primary CCA mode finding [8].
The findings from Smith et al., 2015 suggest that high-dimensional patterns of positive and/or negative SMs can predict high-dimensional patterns of brain connectivity and vice versa.If one (SM or brain connectivity) can be used to predict the other, researchers move one step closer toward understanding the complex relationship between brain, lifestyle, and behavior [9].Establishing a baseline neurotypical signature for this relationship will have important implications for understanding sub-clinical, clinical populations, and other populations (e.g., developmental).
As a precursor to the current study, using code published by the original authors, we computationally replicated their findings in both the HCP 500 and 1200 subject releases to validate our own analysis pipeline [10].

Conceptual replication in an independent dataset
In the current study, we aim to conceptually replicate Smith et al. 2015's positive-negative axis and primary CCA mode findings in novel data from the Adolescent Brain Cognitive Development (ABCD) study.This is not a replication in the sense that we seek to find results identical to those of the original authors (an exact computational replication of the original study is documented in our preprint [10]).Rather, we apply the methodology of Smith et al. 2015 to the ABCD dataset to determine if there exists a strong correlation between connectomes and SMs and dominant CCA mode that explains a significant amount of covariance within a dataset with very different subject and scanner characteristics.The longitudinal ABCD study is an ongoing effort to study the environmental influences on behavioral and brain development in pre-adolescents, recruiting and following over 11,000 children aged 9-10 over a period of 10 years [11].This multi-site study leverages multiple imaging modalities, including resting-state fMRI, to characterize brain development, providing a large neuroimaging dataset for exploration using the tools of network neuroscience.
We selected the ABCD dataset to replicate Smith et al. 2015 for multiple reasons.The 22-35 year-old subjects in the HCP dataset, which was fairly homogenous in terms of subject demographics and scan quality [6], are vastly different from the pre-adolescents in the ABCD dataset in terms of brain development and demographics.Thus, this dataset presents a tremendous opportunity to search for the presence of a positive-negative axis in a heterogenous set of subjects who are still in neurodevelopmental stages.Additionally, the ABCD dataset is one of the largest available that has comparable subject measures to the HCP dataset; since the original study was only conducted using 461 subjects, this increases the statistical power of the analysis.A study by Xia et al. used neuroimaging data of youth (8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22) years old) and found positively correlated patterns of functional connectivity and psychiatric symptoms across 4 dimensions of psychopathology.Each dimension was associated with a distinct pattern of abnormal connectivity while still showing common features (e.g., loss of network segregation) distinguishing them from normal brain function.To account for potential confounds, the authors controlled for age, sex, race, and motion, and they also found that these findings replicated in an independent dataset [12].The implications of their study are threefold: 1) it is possible to perform a CCA-type analysis on highly heterogeneous youth data and obtain interpretable results; 2) the analysis is sensitive enough to distinguish relationships among groups of related SMs and functional connectomes; and 3) there appears to be a common neurobiological mechanism underlying vulnerability to a wide range of psychiatric symptoms, echoing the positive-negative axis finding.Given their findings, we felt that the ABCD dataset was appropriate for the current study despite the age difference and heterogeneity compared to HCP data.
Finding a primary CCA mode and positive-negative axis in the ABCD data would have implications beyond speaking to the robustness of the original study and the presence of the phenomena in a dataset other than HCP.It would suggest that the high-dimensional relationship between brain connectivity and phenotype may be independent of age, demographics, scanner type, pre-processing methods, and even the subject measures, which can vary in how they index cognitive constructs.Such observations would be pivotal for the establishment of a signature neurotypical relationship between connectivity and phenotype, which could then be applied to a variety of clinical populations, potentially increasing its practical application in medicine.Most importantly, finding the positive-negative axis in the baseline ABCD dataset would enable researchers to study it longitudinally over subsequent ABCD releases for the next 10 years.In doing so, researchers could explore how the positive-negative axis changes with regard to age, behavioral characteristics, mental and physical health outcomes, education, demographics, and environment, and possibly even determine which factors form the basis for the positive-negative axis.We aim to design the analysis pipeline such that it can easily be re-run by other researchers who wish to explore how these factors impact the CCA results over time as subsequent ABCD datasets are released.Given these potential implications, we feel that replicating the positive-negative in a heterogeneous dataset is important not only to lend validity to the original study, but also to open new avenues of research in neuroscience.

Replication overview
In this conceptual replication, we followed the same methodology as Smith et al., 2015 to perform a CCA on ABCD subject connectomes and subject measure data (for full details, please refer to section 2. Methods).After applying similar preprocessing steps to the resting-state scan data, we derived a 200-dimension group-ICA functional parcellation.This parcellation was used to calculate subject-level connectomes, which were aggregated to form the connectome matrix used in the CCA.The choice of parcellation can have significant effects on connectivity results [13].For replication purposes we chose to identify functional networks using the same "soft" (e.g., data-derived spatial ICA) parcellation scheme as Smith et al. 2015, which has been shown to have greater predictive power compared to a priori functional and anatomical parcellations [5,14].
The ABCD study has over 60,000 SMs (this includes items from individual questions in an assessment to scan quality control measures), which we quantitatively filtered as per the original study.After removing quantitatively inappropriate SMs, we performed a one-to-one matching of ABCD SMs to their counterparts in the HCP dataset.This process identified 89 ABCD SMs, of which 74 were deemed appropriate for the CCA.We regressed out confounds similar to those in the original study prior to conducting the CCA.
After conducting the CCA, we performed multiple post-hoc analyses from the original study, including generating a positive-negative axis that correlates the SMs (all of those included in the CCA, and 4 which were not) with the primary CCA mode; a hierarchical analysis of the connectomes, and determining which functional connections correlate most strongly with the primary CCA mode; a clustering analysis to identify the major brain regions into which the 200 nodes fall; and a train-test split to evaluate the predictive performance of the CCA model on unseen data.
Our criteria for a successful replication of the primary CCA mode (i.e., the strongest mode) were as follows: 1.The primary CCA mode explains a statistically significant amount of variance in the connectomes and subject measures relative to the null distribution generated via permutation testing 2. The primary CCA mode z-scores for connectomes and SMs are at least a factor of 2 and 3 greater, respectively, than the next largest z-scores for connectomes and SMs (from any of the other modes) 3. A statistically significant correlation exists between the primary CCA mode SM and connectome weights (p<0.001) These criteria attempt to account for the large differences in the ABCD and HCP datasets by focusing on conceptually replicating the primary finding of the original study in the novel dataset, rather than replicating specific numerical results exactly.Although we are particularly interested in finding a single dominant CCA mode, given the large sample size of the ABCD dataset, it is possible that multiple significant modes exist.The appearance of these modes would not necessarily indicate a failure to replicate but would warrant a careful analysis of the results, since in PCA/CCA there can be an arbitrary rotation amongst components that can spread the original result validly across the multiple significant modes [7].Since the positive-negative axis, hierarchical analysis, and train-test split were conducted post-hoc in the original study, there was no definition for what constituted a successful replication for these aspects of the study.

Methods
This article received in-principle acceptance (IPA) at Royal Society Open Science.Following IPA, the accepted Stage 1 version of the manuscript, not including results and discussion, was pre-registered on the OSF (http://doi.org/10.17605/OSF.IO/HMC35).Minor deviations from the protocol are identified in footnotes, and detailed explanations of deviations are provided in section 3.3.Deviations From Protocol.

Participants
The sample composition of the ABCD Release 2.0.1 before and after filtering out subjects is shown in Table 1.The filtering process for excluding subjects is discussed in section 2.3.Data Preprocessing.
Prior to excluding subjects, the sample contained a total of 11,875 pre-adolescents between ages 9 and 10.
After filtering, the sample contained 7,810 subjects . 1 The ABCD dataset contained multiple types of sibships (children with the same parents): monozygotic (MZ) and dizygotic (DZ) twins, as well as siblings and single children.As in the original study, the sibship data was used to generate permutations of the ABCD SM data for statistical testing.Twins were designated as MZ or DZ based on genomic data, and, when zygosity was unavailable, twins were simply treated as siblings.There were 315 subjects that were designated as twins (n=285) or triplets (n=30) but were missing zygosity data; these subjects were considered non-twin siblings for our study to match the procedure followed by  The ABCD functional MRI scans were acquired at 21 sites across the United States, using 3T scanners from three manufacturers: Siemens, Philips, and GE.The fMRI acquisition parameters were designed to be as similar as possible across all scanners: 90x90 matrix, 60 slices, spatial resolution 2.4mm isotropic, TR = 0.8s, multiband acceleration = 6, and multiband echo planar imaging (EPI).Subjects were scanned four times in 5-minute scan lengths for a total of 20 minutes of resting-state scan time.T1w and T2w scans with 1.0mm isotropic resolution were also acquired and used for co-registration purposes.See [15] for more information on the ABCD fMRI scanning protocol.

Subject measure data
The ABCD subject measure battery includes a wide range of assessments of imaging data, biomarkers, cognitive function, substance use and abuse, mental and physical health, and the youth's family and environment [16][17][18][19].Given the longitudinal nature of the ABCD study, the battery was designed such that researchers can characterize subjects' changes in key domains over time.Measures were selected if they could be used through early adulthood or had parallel versions which were appropriate for older ages without their results becoming invalid over time due to repeated assessments or practice effects.A comprehensive computerized battery was administered at baseline during an in-person meeting; shorter follow-up assessments were administered over the phone at 6-month intervals, and more comprehensive follow-up assessments at annual and biennial in-person meetings.Assessments included youth self-assessments, parent assessments of youth, parent self-assessments, and teacher assessments of youth.
For our analysis, only baseline subject measure data were used.

Data Preprocessing
The ABCD fMRI data were obtained from the National Institute of Mental Health Data Archive (NDA).We used Collection 3165 (https://nda.nih.gov/edit_collection.html?id=3165), processed by the Oregon Health & Science University Developmental Cognition And Neuroimaging Lab (DCAN, DOI: 10.5281/zenodo.2587210),which included preprocessed and concatenated resting-state fMRI data . 3 Preprocessing steps in the DCAN pipeline were in accordance with the HCP recommended minimal preprocessing [20].Briefly, the T1-and T2-weighted anatomical images were used to generate subject-specific surface meshes.The functional data were corrected for motion and gradient distortions, and then co-registered to the T1-weighted image and projected into CIFTI grayordinate space, which combines data from cortical gray matter projected onto surface mesh and subcortical gray matter in volumetric space [20].All subsequent preprocessing and analysis steps were conducted within CIFTI space.Finally, a nuisance regression and bandpass filter (0.009hz -0.08hz) were applied to remove artifacts related to head motion and respiration [21,22].Global signal regression was not performed.In addition to the fMRI preprocessing steps applied in Collection 3165, we further removed structured artifacts using an ICA-FIX procedure [23,24] .Collection 3165 also provides censor files for each 4 subject, which denote time points that exceed 0.3 mm in head motion and only segments of 5 or more time points are preserved.To generate the final time series files, we first removed time points flagged by the censor files and then truncated the data to include only the first 10 minutes to ensure that all subjects had the same number of time points.After truncating, we excluded subjects from analysis for having less than 10 minutes of data.This was performed in order to properly process the scans in the group-ICA pipeline, and to mimic the HCP dataset in which all scans were of the same length (albeit a much longer 60 minutes).
It is important to note that the ABCD imaging data used in the current study were processed with the DCAN pipeline which is nearly identical to the HCP minimal processing pipeline but has an additional nuisance regression for respiratory artifacts.In addition, we elected to perform motion censoring since children generally move more than adults during resting-state fMRI scans [25].Aside from these differences, the ABCD resting-state data was prepared in an analogous manner to the HCP resting-state data in the original study.
The subject measure data were obtained from the NDA ABCD RDS release (http://doi.org/10.15154/1504431).The subject measures were provided in a subjects-by-SMs matrix (.Rds file), containing all subjects and their available baseline and follow-up SM data up to the ABCD 2.0.1 release.The data were already processed with qualitative variables factored, summary variables generated, and missing data removed (for complete details, see https://github.com/ABCD-STUDY/analysis-nda).

Subject filtering based on imaging data
To determine which subjects were appropriate for our analysis, we applied the following data inclusion criteria, and included only subjects that met all of the criteria : 5 1.Availability of resting-state scan data and motion summary data in the NDA ABCD collection 3165 a. Subject must have preprocessed and concatenated resting-state fMRI data b.Subject must have a motion summary file associated with the resting-state scan data, and motion data and censoring data must be present in the file 2. Subject must have at least 10 minutes of "good" resting-state scan time (i.e., 750 remaining time points after frame removal/motion correction with a maximum frame displacement (FD) threshold of 0.3 mm) 3. Subject's mean FD must not be anomalous in the overall distribution (i.e., it does not fall within the top 0.25% or bottom 0.25% of the motion distribution) 4. Subject's scans must meet quality control (QC) requirements, based on QC data included in the SM data matrix a.At least one T1 anatomical scan which passed both quality control and protocol compliance checks, for registration b.At least two resting-state scans that passed both quality control and protocol compliance checks Subjects that were missing any of the aforementioned data (scan, motion, or QC metrics data), had less than 10 minutes of "good" scan time after motion censoring, had anomalous FD motion values, and/or did not pass the QC thresholds were removed from our sample.The cutoffs for anomalous mean FD motion values were 0.1663 mm for the top 0.25% and 0.0212 mm for the bottom 0.25% After this initial filtering, 434 subjects were dropped due to criteria one, 1,743 due to criteria two, 9 due to criteria three, and 40 due to criteria four, for a total of 2,226 subjects flagged for removal.As a result, 7,812 out of the original 10,038 subjects remained as possible candidates for our analysis after this stage 6 7 of filtering.

Subject measure filtering
Following the same procedure as Smith et al. 2015, a quantitative filtering was performed in order to remove SM data unsuitable for CCA due to insufficient variance, too much missing data, or potential outliers that could skew the CCA [7,8].Only the baseline measurements for subjects were kept (based on the earliest dated record for each subject), and all follow-up assessments were removed from consideration.
A SM was kept only if all of the following criteria, derived from Smith et al. 2015, were met: 1.There was enough data available a. Defined as at least 50% of subjects having data for a given SM 2. There was sufficient variation in the SM a. Defined as less than 95% of subjects having the same SM value 3. The SM did not contain an extreme outlier value based on the most extreme value from the median a. Specifically, a subject measure contained an extreme outlier if: max(Y s ) > 100*mean(Y s ), where X s is a vector of all subjects' values for an SM s, and vector Y s = (X s -median(X s )) 2   The ABCD 2.0.1 release included a total of 64,148 subject measures.After applying these criteria, 5,699 SMs were dropped due to criteria one, 2,765 due to criteria two, and 41,383 due to criteria three, for a total of 49,847 SMs dropped, resulting in 14,301 out of the original 64,148 SMs passing as quantitatively appropriate for our analysis.For a list of all SMs and why they were dropped according to the quantitative criteria, see Supplementary File 1.
From the remaining 14,301 SMs, we then identified ABCD SMs to include in our study via a one-to-one matching of the remaining ABCD SMs with the 461 HCP SMs published by the original authors (https://www.fmrib.ox.ac.uk/datasets/HCP-CCA/).When possible, exact SM matches were selected (for example, both studies utilized the NIH Toolbox, and thus have identical SMs for this construct); otherwise, SMs were selected to be as similar as possible.Where multiple candidate matching SMs were present, we included all reasonable matches.This one-to-one matching identified 89 ABCD SMs which were deemed to be exact or approximate matches to HCP SMs in the original study.
In the original study, a total of 9 confound variables were identified and regressed out of the SM matrix (reconstruction software version, head motion, weight, height, systolic and diastolic blood pressures, hemoglobin A1C level, brain volume, and intracranial volume) [8].However, some of these SMs were not available in the ABCD dataset or were deemed unsuitable by the quantitative exclusion; specifically, the ABCD dataset did not contain the systolic or diastolic blood pressures, hemoglobin A1C, or quarter/release (i.e., the type of reconstruction software used) data for subjects.Subject height was not included as a confound because it did not satisfy SM quantitative inclusion criteria three.In addition, we included SMs for scanner manufacturers and ABCD scan sites to account for any site-or scanner-dependent variation.Thus, the following ABCD SMs were identified as confounds to be regressed out of the SM data prior to the analysis (the ABCD SM name, as listed on the NDA, is given in parentheses): 1. Scanner site (abcd_site) 2. MRI scanner manufacturer (mri_info_manufacturer) 3. Mean frame displacement (remaining_frame_mean_FD) 4. Weight (antho_weight_calc) 5. BMI (anthro_bmi_calc) 6. Cube-root of total brain volume (smri_vol_subcort.aseg_wholebrain) 7. Cube-root of total intracranial volume (smri_vol_subcort.aseg_intracranialcolume)All 7 confound SMs were demeaned and any missing data were imputed as zero.Additional confounds were generated by demeaning and squaring confound SMs 3-7, for a total of 12 confound variables.This was done to account for any potentially nonlinear effects of these SMs [8].These 7 confound SMs were excluded from the CCA; however, they were regressed out of the SMs that were ultimately input to the CCA.Smith et al. identified and excluded an additional 45 variables that were considered "undesirable" as they were "not sufficiently likely to be measures relating to brain function" or "minor measures highly correlated with more major related measures".We matched eight of our SMs with these variables (sex, age, race, handedness, family income level, parent employment status, parent education level, parent marital status).To be consistent with the original study, they were excluded from the analysis.
After excluding the 7 confounds and 8 undesirable SMs, 74 SMs remained for use as the input to our analysis .At this point we removed any subjects that were missing more than 50% of the final 74 SMs 8 (n=2) .Thus, the final SM matrix had dimensions 7,810 by 74 (subjects-by-SMs) .Although only this 9 10 selected subset of 74 SMs was used as inputs in the CCA, we examined the following 4 SMs in our post-hoc positive-negative axis: BMI, age, household income level, and parent education level.These SMs were selected because they were analyzed post-hoc in the original study, and because we deemed the ABCD measures quantitatively appropriate for a correlation analysis.BMI and age are continuous variables whose magnitude has interpretable meaning.Household income level and education level were factored such that higher scores correspond to higher income or higher education level.For a list of the final 74 SMs, see Appendix A; for the one-to-one ABCD/HCP SM matching see Supplementary File 2.

Group-ICA of ABCD resting-state fMRI data
Using the preprocessed resting-state scan data, we calculated a 200-dimension group-ICA functional parcellation using FSL's MELODIC tool [26].We used MELODIC's temporal concatenation option for multi-session fMRI data (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/MELODIC).After calculating the group-ICA, we used FSL's dual regression utility to derive the individual subject parcellated time series.Specifically, dual-regression stage-1 was used to estimate the node-time series, in which the 200-dimension group-ICA map was used as a spatial regressor against the full time series data for each subject, estimating one parcellated time series per subject.The result is 200 nodes' time series of 750 time points (for 600 seconds of scan time) for each subject.Finally, we derived the subject-level functional connectivity matrices ("netmats") via network modeling using the FSLNets toolbox (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSLNets).Netmats were estimated for each of the 7,810 subjects 11 using partial correlation with L2 regularization (rho = 0.01) via FSLNets' Ridge Regression netmat option (ridgep).The symmetric netmats had dimensions 200x200, where each entry of the matrix was a number representing the strength of the edge between two nodes (indicated by the row and column indices).The netmat values were converted from Pearson's r values to z-scores via Fisher's transformation, and a group average partial correlation network matrix was estimated by averaging the z-scored netmats across all 7,810 subjects.Full correlation netmats were also calculated using FSLNets' full correlation option (corr), and a group average full correlation network was estimated by averaging the z-scored full correlation netmats across all 7,810 subjects.
As in Smith et al., the subject-level functional connectivity matrices were combined to create a single subjects-by-edges matrix.This was accomplished by half-vectorizing each subject's symmetric functional connectivity matrix (200 dimensions) to obtain a vector of length n(n-1)/2 (equal to 19,900 for n=200 dimensions) per subject, and concatenating these subject-level vectors to produce the final 11 The final sample was N=5,013, see section 3.3.Deviations From Protocol for details. 10 Due to the change in sample size (see section 3.3.Deviations From Protocol for details) and SM exclusion the final dimensions of the SM matrix was 5,013 x 73.
subjects-by-edges matrix with dimensions 7,810 x 19,900 .Preparation of this matrix can be 12 accomplished with the code provided in our project repository.

Data preparation
The CCA pipeline was based on the code provided by the original authors at their website (https://www.fmrib.ox.ac.uk/datasets/HCP-CCA/) and the steps described in the original study [8].The CCA script was validated through an exact computational replication of the original Smith et al., 2015 study [10].Aside from differing input connectome and SM data, sizes of matrices, and confounds used, our analysis is identical to that of the original study.We outline the analysis steps here for clarity (shown graphically in Figure 1).In order to run the pre-CCA PCA reduction [7,8] on the SM data, it was first necessary to account for missing data in matrix S 3 (1.73% of data was missing).To do so, we estimated a subjects-by-subjects covariance matrix one SM at a time for matrix S 3 , where, for any pair of two subjects, SMs missing in either subject resulted in the comparison being ignored.The approximated covariance matrix was then projected onto the nearest valid positive-definite covariance matrix using the nearestSPD toolbox [27] in Matlab, resulting in matrix S 4 (dimensions 7,810 x 74), which had no missing data and thus accounted for missing data without needing to impute missing SM values.There was a strong correlation (r=0.9999) between the before and after covariance matrices.
S 4 was input to a 70-dimension PCA, generating matrix S 5 which was the SM data matrix input to the CCA.Although a 100-dimension PCA (as in the original study) was not possible due to having fewer valid SMs than in the original study, Smith et al. noted that there were no statistically significant differences in the final CCA model when the pre-CCA reduction of SMs and netmats was run with a much smaller number of PCA components (specifically, 30 instead of 100).
The functional connectomes were processed in a similar manner.First, the 7,810 x 19,900 subjects-by-edges data matrix was imported into Matlab and stored in N 0 .From N 0 , two matrices (N 1 and 14 N 2 ) were generated, and then horizontally concatenated to form matrix N 3 .N 1 was formed by first demeaning N 0 column-wise, then globally variance normalizing the matrix.N 2 was formed by normalizing each column of N 0 relative to its mean, followed by dropping any columns whose mean values were too 14 Due to a change in sample size, the dimensions of N 0 were 5,013 x 19900. 13Due to a change in sample size and SM exclusion, the dimensions of S 1 were 5,013 x 73. 12 Due to change in sample size, the final subjects-by-edges matrix had the dimensions 5,01 3x1 9,900.

Canonical Correlation Analysis
The CCA was run using the Matlab canoncorr function (from the Statistics and Machine Learning Toolbox), using the 70-eigenvector matrices for the SMs (S 5 ) and connectomes (N 5 ) as inputs.The CCA estimated 70 modes, optimizing de-mixing matrices A and B such that the resulting subject-weight matrices U=N 5 *A and V=S 5 *B (both had dimensions 7,810 x 70 ) were maximally similar to each other.15 We then determined the correlation between the subject connectome and SM weights for each of the 70 CCA modes by correlating the corresponding column pairs in matrices U and V.The Pearson's r value for each CCA mode indicated the strength for which a mode of covariation exists in both the subject connectomes and subject measures.
To estimate the significance of the correlation between the weights of each CCA mode (i.e., correlation between the corresponding columns of matrices U and V), we performed a 100,000 permutation-based significance test in which the CCA was re-run on each permutation of the SM matrix to generate a null distribution of CCA results.As in the original study, the permutations (preserving family structure) were generated using the package hcp2blocks (https://github.com/andersonwinkler/HCP/blob/master/share/hcp2blocks.m)[28] .Sibships were 16 determined based on subject zygosity and their parent/family IDs; all related subjects were matched according to their family IDs.Within these family groupings, MZ twins were permuted, DZ twins were permuted, and non-twin siblings were permuted.After the subject-level permutations were created, entire families that contained the same types of sibships and number of children were permuted.As discussed in section 2.1.Participants, to match the procedures of the original study, any twins/triplets whose zygosity data were unavailable were treated as non-twin siblings when generating permutations.Since there are many related subjects in the ABCD dataset (1,339 subjects were twins, 1,245 were non-twin siblings) the family structure was maintained when permuting so that a valid null distribution could be built.
The null distributions were used to calculate thresholds for significance in the percent variance explained in the connectomes and SMs by each CCA mode.Specifically, the 5 th and 95 th percentiles of variance in the null distributions were used to determine which of the first 20 CCA modes explained a significant percentage of the total variance in the original SM (S 2 ) and connectome (N 0 ) data matrices.Connectomes and subject measures were analyzed separately.

Post-hoc analyses
To understand which SMs were most strongly correlated with the primary mode and determine the directionality of the relationship (a positive or negative correlation), we generated a positive-negative axis that quantified the level of variance explained by each subject measure.This analysis included the final 74 SMs that were input to the pre-CCA reduction, and 4 additional SMs (BMI, age, household income level, parent education level; see Appendix A) for a total of 78 SMs whose raw data was stored in a new matrix, S 7 (dimensions 7,810 x 78 subjects-by-SMs) .Matrix S 8 was then calculated by regressing the 12 17 confound variables (identified in section 2.4.1.Data Preparation) out of matrix S 7 (similar to how matrix S 3 was originally calculated).We correlated each column of matrix S 8 with the primary mode SM weights (i.e. the first column of matrix V) to find Pearson's r values for each of the 78 SMs against the primary mode, and then plotted these correlations on an axis.In addition to the Pearson's r correlation values, we calculated z-score for each SM and the percentage of variance in the primary CCA mode that was explained by each SM.The percent variance was calculated by solving the following expression for each SM: ) 17 Due to a change in sample size and SM exclusion, the dimensions of S 7 were 5,013 x 77. 16Instead of using hcp2blocks, we used a modified version of this script to generate our permutation blocks.The code was modified to handle data from multiple study sites, and can be found at https://github.com/andersonwinkler/ABCD/blob/master/share/abcd2blocks.m.
Where pinv() is Matlab's pseudo-inverse function, vector V 1 is a demeaned copy of the primary CCA mode weights for SMs (a 7,810 x 1 vector) , and vector (also a 7,810 x 1 vector) is a demeaned copy To calculate the correlation between the CCA connectome-modulation weights (from the initial CCA) and the original population mean connectome, we first obtained the relative weights of involvement of the original set of connectome edges by correlating the primary CCA mode connectome weight vector U 1 (i.e., the first column of matrix U) against N 0 , resulting in the "full length" edge weight vector A F1 (a 19,900 element vector), effectively mapping the primary CCA mode onto the original connectome data matrix.A F1 was then correlated against the original population mean connectome.
In addition to the 100,000 permutation-based significance test of the main CCA result, a 80-20 split train-test cross-validation was performed in which approximately 80% (~6,249) of subjects were used in a "training" CCA, and the remaining 20% (~1,562) of subjects were left out as a test validation set .Since We then evaluated the primary mode of the "training" CCA to determine how similar it was to the initial CCA result (the first CCA calculated with all 7,810 subjects, see section 2.4.2.Canonical Correlation 21 Analysis).The primary mode weight vectors for connectomes and SMs (U train,1 and V train,1 , respectively) were correlated against the primary mode weight vectors for connectomes and SMs from the original analysis (vectors U 1 and V 1 ).
To test the performance of the "training" CCA, the CCA connectome and SM de-mixing matrices, A train and B train , from the train dataset were multiplied into the left-out SM and connectome matrices in order to estimate subject weight matrices U test and V test for the test data set.The subject weight matrices were estimated as follows: 21 Due to a change in sample size, the first CCA was calculated with 5,013 subjects. 20Due to a change in sample size, 80% of the sample was N=~4010 and 20% of the sample was N=~1003. 19Due to a change in sample size, the dimensions of the vector were 5,013 x 1.
18 Due to a change in sample size, the dimensions of the vector were 5,013 x 1.
The matrices N 5,test and S 5,test were calculated from the testing subset of the connectome and SM data via steps similar to those outlined in 2.4.1.Data preparation.The connectome and SM weight vectors for the primary mode of the "testing" CCA (i.e. the first column of U test and V test , denoted U test,1 and V test,1 ) were then correlated.In order to measure the significance of correlation between the "testing" CCA primary mode connectome and SM weight vectors U test,1 and V test,1 , we performed a 1,000-permutation test (respecting family structure).This test-train and significance testing process was repeated 10 times, each time with a randomly selected subset of test and train subjects, and the average correlation of U test,1 and V test,1 from all 10 runs was calculated.

Code Availability
The code to run this entire pipeline is available on this project's Github repository (https://github.com/nih-fmrif/abcd_cca_replication).We have also uploaded a compressed archive of the code in the Github repository to the OSF project (http://doi.org/10.17605/OSF.IO/HMC35) which also contains the full text of the pre-registration.Questions about the code may be posted as issues on the github project or emailed to the corresponding author directly.

Canonical Correlation Analysis
The CCA yielded four modes whose correlation between the SMs and connectomes weights were statistically significant by permutation test (Mode 1: r = 0.53, P<10 -5 ; Mode 2: r = 0.44, P<10 -5 ; Mode 3: r = 0.26, P<10 -4 ; Mode 4: r = 0.25, P<0.01, corrected for multiple comparisons of all modes).Scatterplots in Figure 2 illustrate the relationship between SM and connectome weights in CCA Modes 1 and 2 (see Supplementary Figure S1 for Modes 3 and 4).Following the approach of Smith et al., 2015, we color the points by an example SM; here we chose fluid cognition as measured by the NIH Toolbox as it was most similar to the fluid intelligence measure used in Figure 1b from Smith et al., 2015.While fluid cognition does not appear to be associated with the relationship between SM and connectome weights in Mode 1, we observed a positive relationship in Mode 2 such that subjects with greater mode 2 weights for both SMs and connectomes have higher fluid cognition scores.However, the relationship between CCA weights and fluid cognition in Figure 2 is meant to be qualitative.In our post-hoc positive-negative axis analysis we quantitatively relate all SMs to the primary CCA mode.From these results, we determine that Mode 2 is the primary, "strongest" mode according to our replication criteria (see 1.3.1 Replication Overview), since it maximally explains variance in both SM and connectome matrices (criteria 1) and exhibits a statistically significant correlation between SM and connectome weights (criteria 3).The Z-score for SM variance explained by Mode 2 was greater than the next highest Z-score by more than a factor of 2 (criteria 2).
The only replication criteria not met is criteria 2, which addressed the variance explained by the connectome and SM weights.We expected that the connectome variance explained by the primary CCA mode would be a factor of 2 greater than the next highest mode.Instead, we observed that Modes 1 and 2 explained similar amounts of variance in the connectomes (Mode 1: Z=2.99, Mode 2: Z=3.16).Although not part of our criteria, we did note that both of these Z-scores were a factor of 2 larger than Mode 3 (Z=-0.34).With regard to the SM matrix, we expected the primary mode's variance explained Z-score to be a factor of 3 larger than the next highest mode.Instead, we observed that the Z-score of Mode 2 (8.56) was a factor of 2.35 greater than the next highest mode (Mode 3: Z=3.64).Thus while our results are numerically similar to our expectations, the magnitude of variance explained by the primary CCA mode did not reach our a priori hypotheses.
To validate the CCA results, which were calculated on the full SM and connectome matrices, we then submitted the data to a 80-20 train-test split cross-validation procedure.Here we implemented a CCA on 80% of the sample and then estimated the subject weight vectors for Mode 2 in the left-out 20%.The weight vectors were then correlated and compared against 1,000 permutations where we shuffled the SM matrix and performed the same analysis.This procedure was repeated 10 times.The mean correlation between U test and V test was 0.22 ± 0.03, whereas the mean correlation for the permutation was 0.002.In all 10 iterations, the cross-validated correlations were significantly different from the null-permuted correlations (P<0.001),thus confirming our main CCA results.Further, our 80-20 train-test split results are numerically similar to the results from Smith et al., 2015 (i.e., mean correlation between U test and V test in the original paper was 0.25).Notably, while still statistically significant, the mean correlation between the U test and V test (r = 0.22) is smaller than the corresponding correlation from the full CCA (r = 0.44).We attribute this decrease in correlation to the possibility that the CCA model defined on the full dataset may be overfit and a generalizable correlation between canonical variates would likely be smaller in magnitude.
The results presented thus far used data from a preprocessing pipeline that removed structured artifacts using an ABCD-trained ICA-FIX classifier.However, to remain as close to the original study as possible, we also repeated the main CCA and variance analysis on data preprocessed with an ICA-FIX classifier trained using data from the HCP.While there were small numeric differences in the results from these two pipelines, the overall findings are consistent, suggesting that our results are not dependent on which ICA-FIX classifier we use.See Supplementary S2.3 HCP-trained ICA-FIX Results.

Post-hoc Analyses
We next sought to determine the relationship between the CCA modes and specific SMs and connectome edges.Here we focus on Mode 2 as this is the primary mode according to our replication criteria (see

Replication Overview).
As in Smith et al., 2015, we found a positive-negative axis where positively-valenced SMs (e.g., measures of neurocognition, memory, executive function, parental education) correlated positively with the CCA mode, whereas negatively-valenced SMs (e.g., Child Behavior Checklist ADHD, rule breaking behavior, conduct disorder) correlated negatively with the CCA mode.Interestingly, no negatively-valenced SMs were on the positive end of the spectrum and no positively-valenced SMs were found on the negative end.
See Figure 4 for a depiction of the ABCD positive-negative axis for Mode 2 and Supplementary Figure S2 for an unthresholded positive-negative axis that includes all SMs.Further, we also find consistent positive-negative axis results in the HCP-trained ICA-FIX pipeline.See Supplementary Figure S5 for the axis thresholded at r = ±0.2 and Supplementary Figure S6 for the unthresholded axis.Finally, we examined the relationship between specific connectome edges and the primary CCA mode (Mode 2).Here we calculated correlations between subject-specific connectomes and the connectome weights from the primary CCA mode. Figure 5 shows the top 30 edges most strongly related to the primary mode, with node clusters derived from clustering the full set of connectomes using Ward's clustering algorithm.In Figure 5 we observe that a variety of nodes are represented within the strongest edges, including sensory regions as well as higher-order association cortices topographically similar to the posterior midline and temporo-parietal regions of the default mode network.Interestingly, compared to Smith et al., 2015 we do not observe a clear organization in which all nodes topographically similar to the default mode network (e.g., nodes 2, 5, 8, 10, 11, 14, 20, 23, 26 in Figure 5) are clustered together.Instead, these nodes are spread out among different clusters; however, as in Smith et al., 2015, the connectome edges between these nodes are all positively related to the primary CCA mode.See

Deviations From Protocol
On closer inspection of the ABCD dataset, we realized that deviations from our pre-registration were necessary to more closely adhere to the analysis methods used in Smith et al., 2015.All of the deviations detailed below were made while preprocessing the data.The final CCA mode calculation, on which successful replication was defined, was not conducted prior to the deviations and did not influence our decisions.The accepted Stage 1 protocol is available on the project's Open Science Foundation page (http://doi.org/10.17605/OSF.IO/HMC35).
In our pre-registration, we planned to use the fully preprocessed and concatenated ABCD resting-state data distributed by DCAN.The DCAN pipeline removes motion and respiratory artifacts through nuisance regression and filtration [21], whereas Smith et al., 2015 removed structured artifacts using ICA+FIX [23,24].If we were to run ICA+FIX on the already "cleaned" preprocessed DCAN pipeline data, this would be a significant departure from the methods of Smith et al., 2015.We therefore decided to use the DCAN's BIDS-formatted input data and their abcd-hcp-pipeline (https://github.com/DCAN-Labs/abcd-hcp-pipeline,version 0.0.1) to do preprocessing in accordance with the HCP recommended minimal processing (the first stages of the DCAN pipeline are a wrapper for the HCP minimal pipeline, which was used in Smith et al., 2015).
In our stage 1 pre-registered manuscript, we failed to specify which ICA-FIX classifier we would use to classify and remove structured artifacts.Therefore we elected to run two versions of our pipeline: one using an ICA-FIX classifier that we trained using ABCD data and another using the ICA-FIX classifier distributed with the ICA-FIX software which was trained using HCP data.We present the results from the ABCD-trained ICA-FIX pipeline in the main manuscript and the HCP-trained ICA-FIX pipeline in section S2.Since we elected to not use the pre-processed, concatenated fMRI time-series included in the DCAN ABCD release, it was necessary to examine the quality of the four individual fMRI runs for each subject before passing them to ICA-FIX.Here we noticed many runs had low numbers of time points (i.e., less than the full acquisition protocol) and/or high mean framewise displacement (FD) values, which could lead to spurious results.Thus, we applied the following revised data inclusion criteria.The addition of the following inclusion criterion 1 is the only deviation in proctol.Inclusion criteria 2 and 3 remained the same as the pre-registered protocol.Runs that passed the above criteria were then submitted to ICA-FIX to remove structured artifacts.Censor files were then created using an FD threshold of 0.3 mm.To create the censor files, we used head motion parameters generated during pre-processing to flag time points that exceeded an FD of 0.3 mm for removal.In addition, segments of the time series were also flagged for removal if fewer than 5 consecutive time points had an FD displacement below the threshold.After ICA-FIX, runs were censored (i.e., time points flagged for removal were removed from the time series) and concatenated.Finally, we truncated the timeseries to 10 minutes for all subjects.
After the 10,038 subjects in Collection 3165 were filtered according to the new criteria and only those with "good" time points were retained, 5,013 subjects remained as possible candidates for our analysis.This sample is smaller than the 7,810 described in our pre-registered protocol because our new analysis excluded individual runs that lacked sufficient time points or exceeded a mean FD of 0.3 mm.
The subject measure filtering protocol was unchanged.However, given the revised sample of 5,013 subjects, subject measure "resiliency_6b" was removed for having too much missing data (i.e., failed to meet SM inclusion criteria

Discussion
In this study we have shown that a single axis of co-variation spanning 'negative' and 'positive' attributes links diverse participant characteristics with specific patterns of brain connectivity in 9-and 10-year old children.We set out to replicate the landmark paper from Smith et al., 2015, which used the HCP dataset to examine the multidimensional relationship between functional brain network organization at rest and a variety of subject measures, including metrics of cognitive function and lifestyle.Using canonical correlation analysis, the authors found a primary mode of correlation that explained a statistically significant amount of variance in both the functional connectomes and SMs.Further, when correlating the primary mode with specific SMs, the authors found that these correlations were organized along a positive-negative axis, such that positively-valenced SMs (e.g., fluid intelligence) exhibited a positive correlation with the primary CCA mode, and the negatively-valenced SMs (e.g., anger) exhibited a negative correlation.The authors also found that connectivity within the default mode network was positively correlated with the primary CCA mode.
Our analysis met 2 of the 3 strict numerical replication criteria described in our pre-registration (see 1.3.1 Replication Overview).For our first criteria, we found that Mode 2 of our CCA explained a significant amount of variance in both the connectome and subject measure matrices, according to a permutation test.
For our third criteria, the connectome and SM CCA weights for Mode 2 exhibited a statistically significant positive correlation.Only our second criteria was unsuccessful.We expected that the primary CCA mode would explain a majority of the variance relative to the other modes, as indicated by Z-scores of a factor of 2 for the connectome matrix and a factor of 3 for the SM matrix greater than the next largest Z-score.In hindsight, we would judge the specific numeric magnitude factor criteria as overly conservative and not necessarily a critical finding of the original study.Further, one noteworthy difference between our findings and Smith et al., 2015 is that the primary mode in the original paper was Mode 1, whereas we observed Mode 2 as primary.CCA modes are ordered according to the magnitude of correlations between the corresponding linear low-rank projections of the left and right input matrices -that is, mode order relies on the strength of the relationship between the canonical variates [7].While our criteria 3 for a successful replication dealt with the statistical significance of the correlation between these canonical variates, our criteria 1 and 2 focused on the variance that the CCA modes explained in the original SM and connectome matrices separately.Thus, while our observed Mode 1, by definition, explained the most variance in the relationship between SM and connectome matrices, Mode 2 was more successful at explaining the variance in each matrix separately.As a result, Mode 2 met our criteria for the primary mode.Further, while Smith et al., 2015 found a single statistically significant mode, we found four, though only one met our criteria as primary.Subsequent work examining the multidimensional relationship between brain imaging and subject measures suggests that multiple interpretable modes may exist and account for a variety of brain-behavior dimensions [32,33].Despite the difference in which mode was primary, in our view our analysis shows that 9-to 10-year old children exhibit the same basic phenomenon that was reported by Smith et al. in young adults: a surprising amount of variation in brain connectivity can be explained by a single axis of 'positive' and 'negative' attributes.
Further, our analysis has shown that this brain-behavior relationship is robust to significant heterogeneity in acquisition and analyses.First, the HCP dataset used in the original study was collected at a single site on a single scanner [6], whereas the ABCD dataset was collected across 21 sites throughout the US and multiple scanner manufacturers and models [15].Second, we have shown that our findings remain consistent through two preprocessing pipelines: one which removed structured artifacts using an ICA-FIX classifier trained on ABCD data and another which leveraged the ICA-FIX classifier trained on the HCP dataset, which was used in the original study.With such heterogeneity in the ABCD dataset relative to the HCP dataset, it is surprising that the results from Smith et al., 2015 replicate so well, since previous work has shown poor reliability when combining data from multiple sites using methods such as intraclass correlation coefficient [34].While we did observe a lower correlation between SM and connectome weights compared to Smith et al., 2015 (original study: r = 0.87; current study: r = 0.44, for the primary mode), it is possible that the current analysis is not as susceptible to heterogeneity compared to other methods [34] since the data undergo both deconfounding and dimensionality reduction prior to the CCA (which is itself another form of dimensionality reduction between two sets of data [7]).The conceptual replication of a finding that originated using a relatively homogenous dataset such as HCP within the more heterogeneous ABCD dataset offers hope in the face of the so-called replication crisis in psychological science more broadly [35].
The findings from Smith et al., 2015 are also robust to the many differences in demographics between the subjects of HCP versus those of ABCD.The ABCD dataset attempts to reflect the race and ethnicity breakdown of the United States in general [36], whereas the HCP dataset is comparatively more homogeneous [6].Moreover, the current study extends the results from the young adults (ages 22 to 35 years) of HCP into the 9-and 10-year old baseline sample of ABCD.Thus, we find that a primary mode of correlation between functional connectomes and SMs generalizes across multiple sources of variability in data acquisition and subject demographics, suggesting that this high-dimensional relationship may in fact be characteristic of general human cognition.
As in the original study, we also found evidence for a positive-negative axis linking specific subject measures to the primary CCA mode.Here SMs typically regarded as positive qualities (e.g., working memory, executive function, parental income) were positively correlated with the primary mode, whereas SMs typically regarded as negative qualities (e.g., conduct disorder, ADHD symptoms, rule breaking behavior) were negatively correlated with the primary mode.Smith et al., 2015 noted that the primary mode could be evidence for a neural correlate of the general intelligence g factor [37].Interestingly, SMs on the positive end of the axis included both measures of neurocognition as well as two parental measures (highest level of education and income) that have been identified by the World Health Organization as key factors in the Social Determinants of Health [38].They also serve as proxy measures of socioeconomic status (SES) which has been previously shown to have a strong relationship with brain function [39,40].SES is a multidimensional construct [40] and many of the relationships between SES and brain function are nonlinear and/or moderated by other factors [41,42], which warrants some caution in the interpretation of the current study.However the fact that these measures were among the strongest correlations with the dominant mode demonstrates a need for expanding variables related to the social determinants of health in future studies to complement and contextualize standard neurocognitive and clinical measures.
On the negative-correlation side of the positive-negative axis, the appearance of clinical measures such as subscales of the Child Behavior Checklist [43] for ADHD, conduct disorder, rule breaking behaviors, and externalizing symptoms are of particular interest since many symptoms of behavior disorders emerge during adolescence [44].Future work should determine the predictive accuracy of this axis.For example, the ability to predict behavior disorder symptoms from functional connectivity (and vice-versa) could be of great utility to clinicians.In addition, future work should extend this analysis into younger children to examine a possible dynamic emergence and clinically predictive accuracy of the positive-negative axis.
Finally, we also found a specific pattern of functional connectivity edges that was related to the primary CCA mode.The nodes whose edges were most strongly correlated with the CCA included nodes topologically similar to sensory areas and the default mode network.As in Smith et al., 2015, we found that connectivity within the default mode network was positively correlated with the primary CCA mode.The default mode network is a network of distributed brain regions implicated in introspective and abstract thought, social cognition, and autobiographical memory [45].However, while Smith et al., 2015 found nodes representing a broad swath of the default mode network within the edges most strongly related to the CCA mode, we only found that edges between nodes within the posterior midline and temporo-parietal areas were most strongly related to the CCA mode.Further, while the default mode nodes in Smith et al., 2015 all clustered together, we observed that these nodes were represented in multiple clusters.This is consistent with the findings on the development of the default mode network and higher-order association cortices in general, such that these regions undergo a protracted development into early adulthood [46].That is, with age the DMN shifts from a weakly-connected set of nodes to a cohesive network, which could relate to the development of increasingly complex and abstract cognitive abilities [47].We speculate that the aforementioned SES measures are relevant to the instantiation and subsequent development of the brain's DMN and that the positive-negative axis strengthens throughout development, which could serve as a biomarker for at-risk children and adolescents.However, since we only observe that nodes within a subset of the default mode network are most strongly related to the primary mode and that the baseline ABCD sample is a relatively narrow age range (i.e., 9-10 year olds), it is possible that we are only tapping into one slice of a dynamic developmental process.Future work should leverage the upcoming longitudinal releases within ABCD to examine the relationship between functional connectivity edges and a primary CCA mode across adolescence, as well as the relationship between subject measures and the development of the DMN.
In conclusion, the current study sought to replicate Smith et al., 2015 and extends their findings into the more heterogeneous ABCD dataset.We found a primary mode of correlation between brain functional connectivity and subject measures meeting two of three pre-registered numerical criteria.We also found evidence of the positive-negative axis first reported by Smith et al., 2015 where positive subject measures were positive correlated with the primary mode and negative measures were negatively correlated.Finally, like Smith et al., 2015, we also found that connectivity within regions of the default mode network were positively correlated with the primary mode, although this pattern of results in the ABCD dataset was not as dominated by default mode, and connections involving several other brain regions emerged as significantly linked to the primary mode; this may reflect developmental differences.The current study is situated within more recent efforts to examine the multidimensional relationship between multiple imaging derived phenotypes (IDP), physiological measures, genomics, and behavior showing interpretable modes related to phenotypes such as fluid intelligence, handedness and cardiovascular disease [30,31].The ABCD dataset contains a similar diverse array of measures spanning these different modalities.Here we have demonstrated the feasibility of these ICA-based methods with this unique developmental dataset.Given the robustness of the current replication, we expect future studies (using more modalities and more recent approaches to data fusion) will show both replication of previous results (e.g.handedness) and discover new, interpretable modes of variation in brain structure and function related to child and adolescent development.
low (< 0.1), then demeaning the matrix column-wise and finally globally variance normalizing it.The 12 confound variables were then regressed out of N 3 , forming matrix N 4 .N 4 served as the input to the 70-dimension PCA to reduce dimensionality, estimating the top 70 eigenvectors of the connectomes to form matrix N 5 , which was the connectome data matrix input to the CCA.

FIGURE 1 :
FIGURE 1: Outline of data analysis, highlighting the pre-processing stages for subject measures (LEFT SIDE) and connectomes (RIGHT SIDE) immediately prior to the CCA.N 0 -N 5 and S 1 -S 5 are the names of variables in our MATLAB code.

20 the
ABCD dataset contains sibships, the training and test subsets were generated without splitting any families across the subsets.To calculate the "training" CCA, we applied the same steps outlined in 2.4.1.Data preparation to the training subset of the connectome and SM data to calculate the de-mixing matrices A train and B train such that U train =N 5,train *A train and V train =S 5,train *B train .

Figure 3 :
Figure 3: Variance explained in the original (a) connectome and (b) SM matrices.Blue dots and lines correspond to the observed variance explained.The mean variance explained of the null distribution is in black and the 5th and 95th percentiles of the null distribution are in gray.

Figure 4 .
Figure 4.The SMs most strongly related to the CCA Mode 2 are organized by correlation value.SMs were thresholded at r=+/-0.2.SM label font size corresponds to the amount of SM variance explained by the CCA mode (smallest font= 4%, largest font = 38%).Blue font color indicates that the SM was included in the CCA; gray indicates that the SM was not included.

TABLE 1 :
A The 7,810 x7 4 Subjects-by-SM data matrix was imported into Matlab and stored in matrix S 1 .Matrix S 2 13 was formed from S 1 by using a rank-based inverse Gaussian transformation, ensuring Gaussianity for each SM.Matrix S 3 was then generated by regressing out the 12 confound variables (see 2.3.2.Subject measure filtering) out of matrix S 2 .
values for the i-th SM in matrix S 8 .To facilitate interpretation of our ICA and connectivity results, we submitted the nodes derived from the group-ICA to a hierarchical clustering algorithm.As in Smith et al., 2015, we entered the group-averaged full correlation network (all 200 nodes) into Ward's clustering algorithm implemented in Matlab.The low-dimensional results were then examined and assessed for correspondence with large-scale clusters observed in Smith et al., 2015 (e.g., sensory, motor, dorsal attention, and default mode networks).
[29,30,31]st have at least one T1w image that passes quality control and protocol compliance metrics provided by ABCD Using published studies that use the ABCD fMRI data as a guide[29,30,31], a resting-state run was deemed "good" if it met all of the following criteria:1.Mean framewise displacement (FD) for the run is less than 0.3 mm 2. The pre-censoring run length (i.e., length of the raw scan) is at least 50% of the expected length of 380 time points (i.e., it must be at least 190 time points, or approximately 150 seconds long) 1, see section 2.3.2.), leaving 73 final SMs.No subjects in the new 5,013 subject sample were dropped due to missing 50% or more SM data.