An open-access dataset was used in this study¹, which contains a large sample of children (n = 122, 3.5–12 years old, 64 females, nine left-handed), and adults (n = 33, 18–39 years old, 20 females). No participants had any known cognitive or neural disorders. All adult participants and the parents/guardians of the child participants gave written consent. All protocols were approved by the Committee on the Use of Humans as Experimental Subjects at the Massachusetts Institute of Technology.

Participants underwent fMRI scans while watching a short silent version of “Partly Cloudy,” a 5.6-min animated movie with plots eliciting frequent inferences of characters’ mental states (beliefs, desires, emotions) and bodily sensations (particularly pain). The movie began after 10 s of a black screen and 10 s of opening credits. After fMRI scans, all children completed a socio-cognitive behavioral battery measuring theory of mind (ToM) abilities (Gweon et al., 2012). The custom-made 24-item ToM battery (available at https://osf.io/G5ZPV/) involved listening to an experimenter tell a story and answering prediction and explanation questions that required reasoning about characters’ mental states. The questions included varied aspects of ToM abilities including diverse desires and beliefs, false belief, emotion, moral judgment, and interpretation (Gweon et al., 2012). The ToM score used in the analysis corresponded to the proportion of questions answered correctly.

Whole-brain structural and fMRI data were acquired on a 3-Tesla Siemens Tim Trio scanner at the Massachusetts Institute of Technology. All participants were scanned using the standard Siemens 32-channel head coil except for those under age five, who used custom 32-channel phased-array head coils made for younger children. T1-weighted structural images were collected in 176 interleaved sagittal slices with 1 mm isotropic voxels (Adults: FOV = 256 mm; children: FOV = 192 mm). Functional data were collected with a gradient-echo EPI sequence (#slices = 32; TR = 2 s, TE = 30 ms, flip angle = 90°). The details of scanning protocols can be found in Richardson et al. (2018).

We used the data preprocessed by the original authors (Richardson et al., 2018). Specifically, all functional images were first registered to the Montreal Neurological Institute (MNI) template, and registration of each individual’s brain to the MNI template was visually inspected. The registered data were then smoothed using a 5 mm Gaussian kernel. The Artifact Detection Tools² were used to detect timepoints with more than 2 mm framewise displacement (FD) to the previous time point or with a fluctuation in the global signal that exceeded a threshold of three standard deviations from the mean global signal. Additionally, we applied temporal interpolation on artifactual timepoints and regressed out the first five principal components using CompCor (Behzadi et al., 2007). The residual time courses were then simultaneously band-pass filtered between 0.008–0.15 Hz and detrended with the first- and second-order polynomials using 3dTproject. We also excluded the volumes corresponding to the opening credits (first 10 TRs). We adopted a stringent exclusion criterion and removed participants with a mean FD greater than 0.5 mm, leaving 90 children and 32 adults for further analysis. We performed scrubbing by removing TRs with FD greater than 1 mm.

Following preprocessing, we extracted the ROI timeseries from the denoised data using the social brain atlas (Alcalá-López et al., 2018) containing 36 regions of interest (ROIs). The social brain atlas, developed using meta-analyses of social and affective abilities, consists of key regions responsible for social cognition, such as the amygdala, precuneus, medial prefrontal cortex (MPFC), and temporoparietal junction (TPJ), among many other regions responsible for social and affective information processing. The ROI masks were dilated by one voxel along each direction, and then masked by a binary group mask, resulting in 4 ROIs (bilateral temporal pole and cerebellum) being dropped due to little spatial coverage (fewer than five voxels). Next, we extracted the mean timeseries from each ROI mask. As shown in Figure 1, we computed the intra-subject FC by correlating the timeseries using Pearson correlation. For the inter-subject connectivity, Pearson correlation was computed across different individuals (i.e., child-to-adult) within a given ROI (ISC) and across two ROIs (ISFC). More specifically, we used the adults’ timeseries as a reference, and we correlated the children’s timeseries with those of the adults. We chose adults as the reference group, given the previous evidence showing stronger and more coherent neural coupling and activation patterns in adults than children (Cantlon and Li, 2013; Moraczewski et al., 2018). We averaged upper- and lower-diagonal ISFC values to obtain symmetric ISFC matrices. Both connectivity measures were then Fisher-transformed and normalized to have zero mean and unit variance. Since ISC corresponds to the diagonal terms of the ISFC matrix, hereinafter, we will use ISFC to refer to both ISC and ISFC.

PLS regression is a statistical method that identifies a linear relationship between explanatory variables (X, connectivity) and response (Y, ToM scores). Briefly, PLS regression (Wold et al., 1983) seeks to optimize the covariance between the response (Y ∈ R^n×q) and explanatory variables (X ∈ R^n×p) by finding latent components T ∈ R^n×k for ill-posed problems (n < p, fewer samples than explanatory variables). PLS regression aims to find some lower-dimensional latent components T ∈ R^n×k, to be regressed on Y, which can be formulated as follows:

where P ∈ R^p×k and Q ∈ R^q×k are loading matrices (in PLS regression, the loadings are not orthogonal), E and F are random errors following independent and identical normal distribution.

PLS regression is closely related to principal components regression, for which the latent components are simply derived by performing principal component analysis on X itself. By contrast, PLS regression uses both X and Y to construct latent components T as a linear transformation of X:

where W ∈ R^p×q is a matrix of weights.

In our case of a univariate response (q=1), PLS regression iteratively maximizes the squared sample covariance between Y and the latent components under the condition that the latent components are mutually empirically uncorrelated and the columns of W are of unit length. Formally, to find the i^thcolumn of W, PLS regression solves the following optimization question:

where i = 1, …, k and j = 1, …, i − 1. Here, X and Y have mean-centered columns, and the maximum number of latent components is min(n − 1, p).

Once obtaining W, we can then plug it into the Equation 1 and solve Q using ordinary least squares:

Finally, we can find the solution to the regression model Y^=T⁢QT=X⁢W⁢QT=X⁢βP⁢L⁢S by combining Equations 1, 4:

The beta weights β^P⁢L⁢S∈Rn×K are equivalent to the regression weights of multiple regression (although the problem at hand cannot be readily solved using OLS given n < p).

PLS regression is well-suited for questions where the dimensionality of independent variables (#connectivity) is much higher than the number of observations (#participants) and has been previously used in many associating brain patterns with behavior (Yoo et al., 2017; Ching Fong et al., 2018).

Here, we built three models, namely FC, ISFC, and a joint model (ISFC + FC), to predict the ToM scores after controlling for age (linear and quadratic effects), gender, and handedness. Ten-fold cross-validation was used to determine the optimal dimensionality of the latent components, and the cross-validation was repeated 100 times to minimize the influence of sampling variability. The model performance was evaluated by the mean absolute error (MAE) and Spearman correlation between the true and predicted ToM scores. To evaluate prediction performance for the chosen analyses, we ran 50,000 permutation tests to derive empirical null distributions of ToM prediction. We also conducted a partial correlation analysis to ensure that age and head motion did not drive our prediction.

We compared the predictive performance of three models across a wide range of hyperparameters (#components = 1–10). As shown in Figures 2A,B, the optimal predictive performance, as reflected by the lowest mean absolute error (MAE) and highest correlation between the true and predicted ToM scores, was achieved with the joint model (ISFC + FC) with the component number equal to 2. To rule out the possibility that the increased feature dimensionality drove the improved predictive performance, we randomly sampled half of the FC and ISFC edges and re-evaluated the model performance 1,000 times. Despite the slight drop in predictive performance, the joint model using only half of the FC and ISFC features outperformed the models using ISFC and FC separately (paired t-test ps < 0.001). Furthermore, permutation results revealed that the optimal joint model significantly predicted ToM scores (mean rho = 0.34, p < 0.001), as shown in Figure 2C.

We averaged 100 joint models with two components and visualized associated beta coefficients as shown in Figure 3A and showed identified latent components in Supplementary Figure 1. As is typical of a data-driven approach, regions with high beta values were diffusive across the social brain, with no single dominant anatomical pattern. To better delineate the key ROIs that contributed to the prediction accuracy, we conducted a leave-one-out (LOO) analysis by removing each ROI at a time and compared the predictive performance of the reduced model against the full joint model. As shown in Figure 3B, we found that excluding the anterior midcingulate cortex (aMCC), right middle temporal gyrus (rMTG), right TPJ (rTPJ), right anterior insula (rAI) and right supplementary motor area (rSMA), and left inferior frontal gyrus (lIFG) significantly lowered the predictive performance (Bonferroni corrected p < 0.05), suggesting the significance of these regions in our predictive model. We also conducted a follow-up LOO analysis excluding ISFC and FC edges separately. We found that excluding FC profiles of rMTG and lIFG as well as the ISFC profiles of rMTG, rSMA, and ventromedial prefrontal cortex (vMPFC) significantly lowered the model’s prediction performance (Bonferroni corrected p < 0.05).

Given the previous evidence showing high similarity between ISFC (excluding diagonal ISC) and intra-subject movie-viewing FC patterns (Kim et al., 2018; Demirtaş et al., 2019), we conducted a similar confirmation analysis by comparing the two connectivity patterns, which revealed a moderate positive correlation (mean rho = 0.34). Moreover, it should be noted that the FC-ISFC similarity is modulated by the diagonal ISC terms (for a detailed explanation, see Supplementary material), and higher ISC translates to more adult-like brain responses in children. We were interested in examining whether ISFC-FC similarity had any neurobiological significance, given the previous studies showing children with better task performance had more adult-like brain responses (Cantlon and Li, 2013; Cai et al., 2019). Thus, we carried out a follow-up exploratory analysis by correlating individuals’ ISFC-FC similarity with age and ToM abilities. We observed the similarity between ISFC and FC to be positively correlated with age (rho = 0.550, p < 0.001), while not correlated with ToM scores (rho = 0.12, p > 0.05).

By comparing the predictive efficacy of FC and ISFC, our results suggested that there may indeed be synergy between ISFC and FC to predict individual differences. As shown in Figure 2, the joint model using ISFC and FC not only significantly predicted participants’ ToM scores, but far exceeded individual models’ predictive accuracy. More importantly, the reduced joint model, which roughly matched the dimensionality of ISFC and FC, still outperformed individual models by a wider margin (ps < 0.001). We speculate that the superior performance of the model combining the connectivity measures may indicate that distinct cognitive factors may drive these two neural couplings. ISFC theoretically more precisely captures the shared connectivity patterns driven by the dynamic movie stimuli than standard task-FC (Simony et al., 2016), since the physiological and motion artifacts were uncorrelated across subjects. Moreover, empirical evidence has shown that the time-varying ISFC better tracked the movie cues than the time-varying FC (Bolton et al., 2020), and ISFC within the DMN exhibited reliable and distinct patterns during narrative processing (Simony et al., 2016), suggesting that ISFC may be a more cognitively relevant representation of movie-evoked connectivity patterns than the direct measurement of movie-FC. On the other hand, the intra-subject FC measured during movie viewing is a combination of task-evoked FC and individual FC fingerprint, with the FC fingerprint being the dominant one (Gratton et al., 2018). Thus, including FC measures in our model provides additional information about individuals’ FC fingerprints to facilitate the prediction. As a previous study has shown the improved performance of combining intra-subject FCs from multiple conditions for predictive modeling of phenotypic measures (Gao et al., 2019), we further demonstrated the synergy of intra- and inter-subject FC for predicting socio-cognitive abilities using a movie-viewing naturalistic paradigm.

To further elucidate the relationship between the two neural coupling measures, we examined the similarity of individuals’ ISFC (i.e., excluding the diagonal ISC terms) and FC patterns, which previous studies have focused on (Kim et al., 2018; Lynch et al., 2018; Demirtaş et al., 2019). Consistent with the prior work, we found that the two connectivity patterns were positively correlated (Kim et al., 2018; Demirtaş et al., 2019; c.f., Lynch et al., 2018). Both ISFC and FC patterns contained some degree of movie-evoked connectivity, which could have given rise to the overall positive ISFC-FC similarity. We also found that FC-ISFC similarity was significantly positively correlated with children’s age, suggesting more adult-like social brain responses in older children, consistent with the earlier observations (Richardson et al., 2018). Moreover, younger children could have less stable movie-evoked FC patterns as the underlying cognitive processes may be more variable and individualized, while stronger and more coherent movie-evoked FC patterns emerged in the older children, leading to a higher FC-ISFC similarity (Cantlon and Li, 2013; Moraczewski et al., 2018). Nevertheless, we noted that although FC-ISFC similarity was significantly correlated with age, there was no significant correlation between FC-ISFC and ToM scores, justifying the necessity of a more complex predictive model.

Our LOO analysis identified a few critical ROIs that significantly lowered the predictive accuracy in our final joint model. The majority of these key ROIs were found to be part of intermediate (i.e., rAI, aMCC, rSMA, and lIFG) and high-level processing subnetworks (i.e., rTPJ and rMTG), as defined by Alcalá-López et al. (2018). The high-level processing subnetwork consists of key regions closely associated with ToM, such as the TPJ, MTG, posterior cingulate cortex (PCC), and precuneus (Schurz et al., 2014). Per Alcalá-López and colleagues, this high-level subnetwork is more strongly connected within itself than any other subnetworks, and is most likely to be associated with the cognitive categories of social cognition in all subnetworks (Alcalá-López et al., 2018). The significance of the high-level processing network was confirmed when examining the latent components in the final model, which associates the connectivity profiles to ToM scores, as the second component was weighted heavily toward ROIs within the high-level processing subnetwork (Supplementary Figure 1). Indeed, the behavioral task included more higher-cognitive elements of ToM processing, such as false belief. Our LOO analysis also identified a few nodes from an intermediate-level subnetwork. The intermediate-level subnetwork, including the AI, aMCC, and IFG, among others, intertwines with the lower-level visual-sensory subnetwork that handles the preprocessed social-affective environmental inputs (Alcalá-López et al., 2018). AI is known for its role as a bridge between large-scale brain networks. Together with the aMCC, these two ROIs within the intermediate networks were involved in empathy and pain-related processing (Kurth et al., 2010; Lamm et al., 2011) and were also found among the ROIs that were highly predictive of ToM scores. Finally, the absence of lower-level visual-sensory ROIs among the most predictive regions suggests that although movie stimuli induce wide-spread connectivity changes in visual areas within and across individuals (Hasson et al., 2004; Lynch et al., 2018), such connectivity changes may be less socio-behaviorally relevant as compared to those in the high association cortices.

Several limitations require further consideration. Firstly, we used adults as the reference group when computing ISFC. Our choice was justified by previous studies showing children with more variable inter-subject neural coupling and weaker activation patterns than adults across many cortical regions (Cantlon and Li, 2013; Moraczewski et al., 2018). Therefore, we believe the adults may be well-suited as a reference group given their higher homogeneity. Future investigations could focus on the impact of choosing different reference groups and better delineating developmental effects. A second limitation is that we chose a set of a priori ROIs using the social brain atlas since we are interested in predicting individuals’ social abilities. Future studies could examine a different parcellation, such as whole-brain parcellations (Power et al., 2011; Craddock et al., 2012; Shen et al., 2013), given the success of whole-brain FC predictive models (Shen et al., 2017; Beaty et al., 2018; Lake et al., 2019; Xiao et al., 2021). Alternatively, future studies could use atlases developed specifically for children since the brain parcelation used in the current study was derived from adult studies or conduct hyperalignment for better functional alignment (Haxby et al., 2020). Thirdly, given the limited sample size for studying individual differences in behavior (Marek et al., 2022), our results are subject to further validation. Future work should use consortium-level datasets and include different movie clips. Another potential remedy is using more advanced learning models to improve predictive performance, such as sparse group LASSO (Bai et al., 2020) and sparse tensor decomposition (Zhang et al., 2021). Fourth, we only studied static FC and ISFC, while future work could investigate whether dynamic FC (Xie et al., 2018b) and dynamic ISFC (Bolton et al., 2018) can be jointly studied to better characterize neural dynamics and individual differences. Finally, the current dataset utilized a ToM measure that was a composite of varied ToM abilities, which may account for the varied networks predictive of performance. While ToM is often discussed as a unitary construct, much theoretical and empirical work has demonstrated its multifaceted nature, including weak or no associations between measures (Schaafsma et al., 2015; Gernsbacher and Yergeau, 2019; Warnell and Redcay, 2019). Thus, the regions predictive of ToM should be interpreted within the context of this measure. Different measures may identify distinct patterns of predictive performance.

The current study investigated the potential benefits of jointly studying intra- and inter-subject connectivity, namely, FC and ISFC, respectively, and whether the two connectivity measures combined could enhance the prediction of individual differences in social cognition measured by a ToM questionnaire. Using PLS regression, we showed that the joint connectivity (ISFC + FC) model outperformed individual models even after matching the feature dimensionality. Our results suggest that intra- and inter-subject connectivity may encode unique and complementary information about the individuals’ social abilities, and we shall make full use of both connectivity measures to gain enriched insight into neural processes underlying naturalistic fMRI paradigms.