Building Functional Connectivity Neuromarkers of Behavioral Self-Regulation in Children with and without Autism Spectrum Disorder

Children with Autism Spectrum Disorder (ASD) are known to struggle with behavioral self-regulation, which associates with greater daily-life challenges and an increased risk for psychiatric comorbidities. Despite these negative outcomes, little is known about the neural expression of behavioral regulation in children with and without ASD. Here, we examined whole-brain linear associations between brain functional correlations (FC) and behavioral regulation through connectome predictive modelling (CPM), a data-driven protocol for developing predictive models of brain–behavior relationships from data, assessing ‘neuromarkers’ using cross-validation. Using data from two sites of the ABIDE II dataset comprising 276 children with and without ASD (8-13 years), we identified functional brain networks whose FC predicted individual differences in two, of three, behavioral regulation subdomains. These distributed network models predicted novel individuals’ inhibition and shifting from resting-state FC data both in a leave-one-out, as well as split halves, cross-validation. We observed commonalities and differences in the functional networks associating with these subdomains, with inhibition relying on more posterior networks, shifting relying on more anterior networks, and both involving regions of the DMN. Our findings present a substantial addition to our knowledge on the neural expressions of inhibition and shifting across the spectrum of children with and without ASD, demonstrating the utility of this trans-diagnostic modelling approach. Given the numerous cognitive and behavioral issues that can be quantified dimensionally in neurodevelopmental disorders, further refinement of whole-brain neuromarker techniques may thus pave a way for functional neuroimaging to meaningfully contribute to individualized medicine.

brackets) are given for the total samples comprised of both TD and ASD participants, which were used to build the models, as well as 4 motion artifacts. We chose a conservative threshold ('aggressive') in order to decrease the chance of 13 false positives. In other words, more components are removed as this threshold is more conservative 14 about what is retained. Noise components identified by AROMA were removed from the data. Third, 15 images were de-noised by regressing out the six motion parameters, as well as signal from white 16 matter, cerebral spinal fluid and the global signal, as well as their first-order derivatives [59]. While 17 there is currently no gold standard [60] regarding the removal of the global signal, we chose to 18 remove it based on recent evidence that it relates strongly to respiratory and other motion-induced 19 signals, which persist through common denoising approaches including ICA and models that attempt 20 to approximate respiratory variance [61]. Motion (defined as each participants's absolute maximum 21 displacement) was substantially reduced following this procedure (before: 1.28mm ± 0.85mm; after: 22 0.05mm ± 0.07mm). As a final step, as described in more detail below, head motion was incorporated 23 12 into models by removing connections that remained significantly (p<0.05) associated with z-scored 1 motion before cleaning in a Pearson's correlation [62]. 2 3 2.7 Connectome-Predictive Modelling 4 To elucidate how behavioral regulation skills are reflected in children's whole-brain FC profiles (or 5 'connectomes'), and how they vary across children with and without ASD, we utilized a protocol 6 termed Connectome Predictive Modelling (CPM). CPM is an algorithm for building predictive 7 models based on participants' FC matrices, and for testing these models using cross-validation of 8 novel data. Scripts are written in MATLAB and are freely available at 9 www.nitrc.org/projects/bioimagesuite. The CPM protocol is described in detail in [46] and has 10 previously been applied to pediatric data sets including data from the ABIDE sample [42, 62-64]. We 11 followed the CPM protocol [46], as well as recent recommendations for predictive modelling [65], in 12 calculating each participant's FC profile, building models of behavioral regulation, and in running the 13 following analyses: (1) a leave-one-out cross-validation to evaluate the potential of models to predict 14 an unseen participant's score, where N-1 participants are used to build the predictive model and the 15 model is subsequently tested on the left-out participant; (2) a split-halves prediction where all 16 available data was randomly split and models built in the first half were used to predict individuals in 17 the second half and vice versa; (3) a site-to-site prediction where models built in the GU site were 18 used to predict individuals in KK site and vice versa. We describe how we calculated FC matrices, 19 built the models and ran these analyses in the following sections. 20 21

Calculation of FC profiles 22
A functionally defined atlas, consisting of 268 cortical and subcortical regions-of-interest ('nodes') 23 13 that cover the whole brain [66], was used to define ROIs. For each child, we extracted the timecourse 1 of each node by taking the mean across voxels and a 268x268 connectivity matrix was calculated 2 between timecourses of node pairs using Pearson's correlation followed by Fisher's Z transformation. 3 Thus, each connection (or 'edge') in the matrix represents the strength of FC between two nodes, and 4 the matrix as a whole represents a child's FC profile or functional connectome. 5 6 2.7.2 Building FC Models of Behavioral Regulation 7 Models were built relating z-scored behavioral regulation subscales (emotional control, shift, inhibit) 8 to FC matrices across participants with and without ASD from both sites. Prior to modeling, effects of 9 motion and site were eliminated from participants' FC profiles by masking out connections that were 10 significantly (p<0.05) associated with motion in a Pearson correlation or different between sites in a t-11 test. 6915 out of a possible 35778 (=268x267, adjusted for the diagonal, divided by 2, because 12 matrices are symmetric) nodes were eliminated at this step due to motion, leaving 28863 valid nodes 13 in the matrix; accounting for site brought this number down to 23100. For model building, each edge 14 in the matrix was correlated with the behavioral regulation measures (again in a Pearson's 15 correlation), and only significantly correlated edges (p<0.01) were selected and retained. These 16 selected edges were first separated into positively and negatively associated edges based on the 17 direction of the correlation, as they may be interpreted differently in terms of their functional roles, 18

Cross-Validation: Leave-One-Out Prediction 3
To evaluate the potential of models to predict an unseen participant's score, one participant was 4 removed from the dataset and the remaining participants (N-1) are used to build the predictive model. 5 The left-out-participant's score was predicted based on the N-1 sample's fit of the linear regression 6 model, and this step is repeated in an iterative manner with a different participant left out in each 7 iteration. Spearman's r s was used to evaluate model performance i.e. comparing actual to predicted 8 scores because it is less sensitive to the effect of outliers than Pearson's r and because CPM 9 predictions are best considered relative rather than absolute [62,63]. Only models that showed a 10 significant correlation at p<0.05 between observed and predicted scores at this step were subjected to 11 follow-up testing. 12 13

Evaluation of the Predictive Model 14
The predictive potential is assessed by comparison of the predicted and observed scores in the full 15 model, and statistical significance is assessed using permutation testing (5000 iterations). Permutation 16 (i.e., randomization) testing was used to assess significance because the assumption underlying the 17 standard r-to-p conversion that was employed in the leave-one-out cross-validation (see above) is As a further test of the models, data from both sites were randomly split while retaining the same 5 number of TD and ASD participants and the same number of participants from each site. Models 6 were built in the first half and used to predict individuals in the second half and vice versa. Participant 7 characteristics for the two split halves are given in Table 2 deviations (in brackets) are given for the total samples comprised of both TD and ASD participants, 12 16 which were used to build the models, as well as for TD and ASD participants separately. Motion 1 (mm) refers to the absolute maximum displacement at any timepoint in the resting-state fMRI scan 2 prior to motion mitigation and denoising procedures. n=number of participants; m=male; f=female; 3 L=left-handed; A=ambidextrous; R=right-handed; IQ=Intelligence Quotient; FIQ=full scale IQ; 4 VIQ=verbal IQ; PIQ=performance IQ; SRS=Social Responsiveness Scale. SRS, inhibition, shifting 5 and emotion control are given as T scores. *denotes deviating numbers in the Inhibition models, for 6 which an additional six outliers (>3 SD in score; all ASD participants) were removed. **SRS scores 7 were not available for seven participants in Split Half 1; 1 participant with ASD was removed as an 8 outlier from Split Half 1. 9 10

Cross-Validation: From Site 1 to Site 2 Prediction 11
To evaluate the potential of models built in one site to predict an unseen participant's score from the 12 other site, models were built in the GU site and used to predict individuals in the KK site and vice 13 versa. 14 15

Model Specificity to Behavioral Regulation 16
To evaluate whether our models capture behavioral regulation dimensionally or are driven by the 17 categorical difference in scores due to ASD diagnoses, we examined the relationship between 18 predicted and observed scores for TD children and children with ASD separately. We further 19 examined (a) whether our models could be driven by motion or IQ, (b) how they relate to the core 20 symptoms of ASD, and (c) whether models built for one domain of behavioral regulation were 21 specific to that domain by computing cross-correlations between (c1) predicted scores and SRS 22 scores, and (c2) predicted and observed scores of different behavioral regulation domains. 23 17 1 2 3 RESULTS 3

Sample Characteristics 4
Characteristics for all samples (GU site, KK site, combined sites; split half 1, split half 2) are given in 5 Tables 1 and 2. Scores for all three subscales of behavioral regulation were significantly higher in 6 children with ASD than for TD children in all analyzed samples, reflecting greater challenges with 7 inhibition, shifting, and emotional control (Table 3). We further observed significantly higher scores 8 on social responsiveness, reflective of greater ASD symptoms, and several significant differences 9 within and across some of the samples in age, IQ, sex and head motion ( Table 2). As expected, the 10 three subscales of behavioral regulation exhibited correlations with each other as well as to social 11 responsiveness and, to a lesser degree, IQ and head motion (Table 4). There were no significant 12 differences in handedness between TD children and children with ASD in any of the samples.   Table 3. Significance of group differences between TD children and children with ASD. Scores 2 for all three subscales of behavioral regulation were significantly higher in children ASD than for TD 3 children in all analyzed samples, reflecting more issues with inhibition, shifting, and emotional 4 control. We observed that in some, but not all samples, ASD children moved more and had lower IQs. 5 In addition, akin to the ratio in the general population, ASD children tended to be predominantly 6 male. Results are given as p-values of independent sample t-tests, which are adjusted in cases of 7 unequal variance as assessed through Levene's test. Motion (in mm) refers to the absolute maximum 8 displacement at any timepoint in the resting-state fMRI scan before and after motion mitigation and 9 denoising procedures. SRS=Social Responsiveness Scale (total score); IQ=Intelligence Quotient; 10

FC Models of Behavioral Regulation 1
Significant models were built using negative edges for inhibition (r s =.23, p=0.0001) and shifting 2 (r s =.19, p=0.001), using leave-one-out prediction (N-1 at every iteration). Positive edge models were 3 not significant for these measures and neither positive nor negative edge models were significant for 4 emotion control (r s <.05, p>0.4). The FC model of inhibition was significant by permutation testing 5 (r s =.23, p=0.037), while the FC model of shifting fell just shy of significance (r s =.19, p=0.067). As 6 seen in Figure 1, inhibition revolved around edges in the somato-motor, visual and cerebellar 7 networks and was more posterior, while shifting appeared more focused on edges around the 8 frontoparietal and dorsal attention networks and was more anterior. Both inhibition and shifting 9 included a number of edges connecting with DMN regions as well as the temporal lobe. Note that in 10 the leftmost panels higher rank refers to a lower score, i.e. lower symptoms. In negative edge models, 11 lower FC associates with higher ranked scores. p=0.71). Therefore, no cross-prediction from GU to KK and vice versa was attempted. 18 19

Model Specificity to Behavioral Regulation 20
In the combined model, Spearman rank correlations between observed and predicted score ranks for 21 TD children (n=191) and children with ASD (n=85) separately were insignificant for the smaller ASD 22 group in both inhibition (r s =.14, p=0.22) and shifting (r s =.10, p=0.36), but near significant in the TD 23 23 group in inhibition (r s =.14, p=0.053) and significant in shifting (r s =.22, p=0.002) (Figure 2). 1 Spearman correlations between predicted shifting or inhibition and motion or IQ were insignificant 2 (r s 's<.10). Both predicted inhibition (r s =.16, p=0.01) and shifting (r s =.23, p<0.001) scores associated 3 significantly with total SRS scores; this association fell just shy of significance in a partial correlation 4 when controlling for diagnosis (inhibition and SRS r s =.11, p=0.069; shifting and SRS r s =.12, 5 p<0.051). In addition, predicted shifting scores significantly associated with observed emotional 6 control scores (r s =.14, p=0.018), but not inhibition (r s =.06, p=0.33). Similarly, predicted inhibition 7 scores weakly associated with observed shifting scores (r s =.12, p=0.04), but not emotional control 8 scores (r s =.1, p=0.09). Neither shifting nor inhibition significantly associated with motion, neither 9 before cleaning (r s <0.05) nor after (r s <.11). Note all associations reported in this section are 10 uncorrected for the number of comparisons performed. Children with ASD are known to struggle with behavioral self-regulation, which associates with 2 greater daily-life challenges and an increased risk for psychiatric comorbidities [11, 13, 67, 68]. In a 3 fully cross-validated, data-driven analysis in a large sample of children with and without ASD that 4 was compiled across two data collection sites, we identified distributed patterns of FC whose strength 5 predicted individual differences in two behavioral regulation subdomains. These whole-brain network 6 models predicted novel individuals' inhibition and shifting, but not emotional control scores from 7 resting-state FC data both in a leave-one-out, as well as in a split halves cross-validation, providing 8 evidence that meaningful correlates of behavioral regulation in intrinsic brain patterns exist. 9 Indicating the limitations of this approach, whole-brain network models could, however, not be built 10 within the smaller and less balanced samples collected at the respective sites. We further found that 11 although models captured within-group variation, models built on inhibition and shifting also 12 predicted ASD symptoms more generally, and this relationship persisted at a trend-level after taking 13 diagnosis into account. Overall our results, showing that complex brain network models predict 14 different measures of behavioral regulation across a sample of children with and without ASD, 15 demonstrate that whole-brain FC data can serve as a holistic neural index of inhibition and shifting. Our work further evidences two major challenges that remain towards achieving one of the primary 21 goals of human neuroimaging -to identify generalizable neuromarkers of clinical utility. First, whole-22 brain network models could not be built successfully within the smaller and less balanced samples 23 27 collected at the respective sites. It is a known problem in ASD research that there are site effects that 1 appear to influence generalizability, meaning that it may be unlikely to build generalizable predictive 2 markers using data from a single site. Second, models built on inhibition and shifting scores also 3 predicted ASD symptoms more generally. While this is to be expected given the known association 4 between ASD core social symptoms and behavioral regulation, it is worthwhile to note that since 5 these relationships persisted at trend-level when diagnosis was controlled for, they are unlikely to be 6 entirely driven by diagnostic category. In addition, relationships between observed and predicted 7 shifting scores remained significant when assessed only within the TD population, again indicating 8 that the predictions have some specificity to the constructs under investigation. Nonetheless, the 9 correlation between social and behavioral regulation traits makes it more difficult to tease apart the 10 unique and shared features between these constructs. The current study has several distinct strengths, which include the use of three measures of behavioral 8 regulation in two relatively large, "matched" groups of children with and without ASD and novel 9 preprocessing techniques. Both samples were acquired in close spatial and temporal proximity, that is, 10 around the same time and around the same geographic location (Washington, D.C. and Baltimore). 11 The assessment of behavioral regulation is well validated [16] and although parent-and self-reports 12 are subjective, they capture a measure of behavior integrated over a longer time frame than can be 13 The study also has several weaknesses, including differences within and between the GU and KK 21 sites. These are differences in (i) fMRI data acquisition procedures, (ii) in-/exclusion criteria, (iii) 22 numbers of participants, (iv) TD vs. ASD ratios, (v) sex ratios, (vi) IQs, (vii) BRIEF scores, (viii) 23 29 comorbidities, and (ix) medications. All of these serve as a reminder of just how divergent ASD 1 populations are, how difficult it is to assess whether a sample reflects the "true patient population" 2 and to what extent research findings are truly generalizable. Considering that these two sites contain a 3 notable number of participants also highlights how common datasets with unbalanced groups are -a 4 known problem when performing cross-validation and assessing prediction performance [65]. 5 Another weakness potentially rests in the removal of motion-associated edges where motion 6 correlates with the behavior that is investigated. While in the combined sample, motion did not 7 significantly associate with the behavioral regulation measures after multiple comparison correction, 8 there remains a possibility that edges that would have associated to, and strengthened the behavioral 9 regulation models, were removed. The same holds true for site-specific associations. Finally, 10 emotional control scores were not significantly predicted for novel individuals. It is perhaps 11 unsurprising that emotional control results did not follow the same pattern as the other measures of 12 behavioral regulation, as evidence suggests that definitions of emotional control abilities are broad, 13 difficult to dissociate from emotion generation processes, and fairly hard to capture reliably [25, 99-14 101]. Another possibility is that the neural mechanisms of emotional control are simply not reflected 15 in whole-brain functional connectivity patterns consistently across individuals with ASD. 16 17 It should be noted that like most methodological approaches, CPM has both advantages and 18 challenges. One major advantage is that because CPM models are defined and validated with 19 independent data, they promise to improve our ability to uncover generalizable brain-behavior 20 associations [65, 102]. A major challenge is that predictive models based on FC will only ever 21 account for a fraction of the variance, because they are limited by how much information the signal 22 can capture as well as the chosen phenotypic measure. They are also bound by linearity assumptions: 23 30 Linear models built across TD children and children with ASD can capture dimensional associations, 1 but may miss categorical differences, which are distinct from dimensional associations (please see 2 [94,95] for a discussion on this). Finally, one may consider the current CPM framework as perhaps a 3 bit simplistic in that it yields only one summary value for 'positive' and 'negative' networks, cannot 4 capture flexible brain network dynamics, and has no 'blueprint' for how to tie together predictions on 5 a multitude of behavioral measures. With this being said, as Scheinost  The characterization of differential presentation of ASD core symptoms and ASD associated 5 symptoms along the spectrum has been of immense interest to researchers, as heterogeneity is a 6 defining characteristic of ASD that makes it difficult to assess and treat within traditional diagnostic 7 frameworks. In this study, we utilized a data-driven approach to develop objective quantitative FC 8 models that elucidate and predict performance in behavioral regulation subdomains in ASD. We 9 observed both commonalities and differences in the functional organization of inhibition and shifting 10 across TD children and children with ASD, with inhibition relying on more posterior and shifting 11 relying on more anterior brain networks. We also demonstrate the generalizability and trans-12 diagnostic utility of this approach, as well as its clear limits to date. Given the numerous cognitive 13 and behavioral issues that can be quantified dimensionally in neurodevelopmental disorders -and 14 literally any disorder or disease -further refinement of whole-brain neuromarker techniques will 15 prove vital to neuroimaging's role in individualized medicine. CR and SB designed the study and analyzed and interpreted the data. SK made substantial 7 contributions to the analysis through adaptation of scripts and computational testing. CR wrote the 8 manuscript with input from SB. All authors read and approved the final manuscript. 9 10 Competing Interests 11 The authors have no conflict of interest to declare. 12 13