Stability and variation of brain-behavior correlation patterns across measures of social support

TL;DR

This study explores how brain-behavior patterns change with different types of social support, revealing flexible brain systems influenced by social relationships.

Contribution

The paper introduces a novel application of multilayer modeling to analyze brain-behavior correlations across social support measures.

Findings

Subcortical and control brain systems show high flexibility in brain-behavior associations based on social support types.

Community affiliations of brain network nodes vary across different social support measures.

The multilayer approach offers a generalizable method for studying brain-behavior correlation stability.

Abstract

The social environment has a critical influence on human development, cognition, and health. Research in health psychology and social neuroscience indicate an urgent need to understand how social relationships are associated with brain function and organization. To address this, we apply multilayer modeling and modularity maximization—both established tools in network neuroscience—to jointly cluster patterns of brain-behavior associations for seven social support measures. By using network approaches to map and analyze the connectivity between all pairs of brain regions simultaneously, we can clarify how relationships between brain regions (e.g. connectivity) change as a function of social relationships. This multilayer approach enables direct comparison of brain-behavior associations across social contexts for all brain regions and builds on both ecological and developmental neuroscientific findings and network neuroscientific approaches. In particular, we find that subcortical and control systems are especially sensitive to different constructs of perceived social support. Network nodes in these systems are highly flexible; their community affiliations, which reflect groups of nodes with similar patterns of brain-behavior associations, differ across social support measures. Additionally, our application of multilayer modeling to patterns of brain-behavior correlations, as opposed to just functional connectivity, represents an innovation in how multilayer models are used in human neuroscience. More than that, it offers a generalizable technique for studying the stability and variation of brain-behavior associations.

Figures

A schematic of the method for computing the correlations between functional connectivity and social support. (a) shows the 92 (one per subject) functional connectivity matrices, where elementijindicates the Pearson correlation between nodesiandj. (b) represents all 92 subjects scores on seven multi-item measures of perceived social support:Friendship, Stress, Rejection, Hostility, Emotional Support, andInstrumental Support. The matrices in (a) are vectorized and stacked on top of each other in (b) to yield a 92 (number of subjects) by 115921 (number of edges on each functional connectivity matrix) matrix. We compute the correlation between social support measures (b) and edge weights (c), resulting in seven (one for each measure) correlation matrices (d), where elementijin matrixkgives the correlation between social support measurekand the edge weight between nodesiandj. Because modularity maximization treats negative weights as belonging to just one community, we transform each correlation matrix to a similarity matrix (e). We do this by computing the similarity between every pair of rows in the correlation matrix. Each of the seven similarity matrices becomes itself a layer in the multilayer model (f). The multilayer model is a 2D tensor (f) whose on-diagonal elements are the seven correlation matrices and whose off-diagonal elements are the coupling between layers. We perform community detection on this multilayer model, which yields a community assignment for every node for each layer (g). We compute for each node its flexibility, or how much it changes its community affiliation between layers (h).

System-level correlations between edge weights and social support measures. We shuffled social support scores across subjects and recalculated the correlations with edge weights 1000 times. For each measure, we summarize the results of this permutation testing by showing the log proportions of edges of each canonical subsystem (e.g., ContA, ContB) that are significantly correlated with the given social support measure using a false discovery rate of 0.05 (adjustedpcrit=0.046).

Analyses of parameter space. We performed PCA on the entire parameter space; (a) shows the magnitude of expression of the first principal component. The point in parameter space for which the magnitude of expression of the first principal component is greatest isγ=0.15, ω=0.001, indicated by a black outline. (b) Shows the variability in node flexibility across parameter space computed across individual partitions at a given point in parameter space. The point in parameter space that maximizes variability in flexibility isγ=0.55, ω=0.0316, indicated by a black outline. In (c), for every point in parameter space, the similarity between every partition was computed using the adjusted Rand Index then averaged for each point across layers. For this calculation, the global maximum similarity is achieved by five points in parameter space:γ=0.9, ω=0.1; γ=0.85, ω=0.3162; γ=0.9, ω=0.3162; γ=0.85, ω=1andγ=0.9, ω=1, indicated by a black outline. We also calculated the similarity of each individual partition to the consensus communities at that point in parameter space and took the average across partitions and layers (d). For this calculation, the global maximum was reached by a set of five points in parameter space:γ=0.9, ω=0.1; γ=0.85, ω=0.3162; γ=0.9, ω=0.3162; γ=0.85, ω=1andγ=0.9, ω=1, indicated by a black outline. Panels (e) and (f) show seed-based analyses, where we computed the similarity of all partitions for each point in parameter space to the seed. For (e), we used as a seedγ=0.55, ω=0.0316, which maximizes the variability in flexibility. For (f), we used as a seedγ=0.35, ω=1, which is the point in parameter space whose communities are most similar to the Yeo systems. (g) Shows the extent to which different regions of the cortex and subcortex express the mode of flexibility present in the first principal component.

Multilayer consensus community structure. Panel (a) shows the consensus community structure across social support measures at the point in parameter space whose coefficient for principal component 1 has the largest magnitude (γ=0.15, ω=0.001). Colors indicate communities. There are 14 communities total, but only four communities are present across all layers (shown in dark blue, light blue, teal, and red). Panel (b) shows the number of communities in each layer according to the consensus partition. Panel (c) gives the average community size for each layer according to the consensus partition. The MNI coordinates for the five panels of subcortical communities are, from left to right: x = -23, x = 10, z = -3, x = 13, x = 27.

Flexibility. Panels (a), (b), (c), and (d) represent data from the point in parameter space (γ=0.15, ω=0.001) with the largest loading on the first principal component. Panel (a) shows the flexibility of cortical and subcortical regions. Panel (b) shows the distribution of flexibility across canonical systems, as defined byYeo et al. (2011);Schaefer et al. (2018), ordered by median flexibility. We note that because trends largely hold for both the system and subsystem scale (e.g., Control and Subcortex are more flexible than other systems, even though Limbic and DefaultC are more flexible than many Control or Subcortex subsystems), we present the results for subsystems here and then proceed with system-level reporting for the remainder of our results. Panel (c) shows the proportion of all nodes that change community affiliation between each pair of the seven layers: Friendship, Stress, Rejection, Hostility, Emotional Support, Loneliness, and Instrumental Support. Panel (d) shows how similar all systems are to each other in terms of the proportion of nodes changing community between each pair of layers. Panel (e) plots the correlation between each system and the global pattern of community change between layers. The MNI coordinates for the five panels of subcortical flexibility are, from left to right: x = -23, x = 10, z = -3, x = 13, x = 27.

Full text

Introduction

1

Research in health psychology and social neuroscience indicate an urgent need to understand how social support is related to whole-brain connectivity (Xiao et al., 2021). Brain networks (and their features) are increasingly used as biomarkers for mental illness and disease (Lydon-Staley & Bassett, 2018;Rashid & Calhoun, 2020;Song et al., 2018). At the same time, the particular patterns of brain network organization that are associated with poorer mental health and cognition depend on the social environment (Ellwood-Lowe et al., 2021;Rakesh, Seguin, et al., 2021;Ramphal et al., 2020;Xiao et al., 2021). Thus, before translating proposed biomarkers to diagnostic or healthcare settings, it is important to characterize how brain network organization varies across social contexts. This endeavor involves linking many strands of research; to date, studies have connected perceived social support to wellbeing (Cohen, 2004;Williams et al., 2018), experience in the environment to brain network organization (Chan et al., 2018;Ellwood-Lowe et al., 2021), and the social world to the brain (Krendl & Betzel, 2022;Morawetz et al., 2021;Rudolph et al., 2021;Schmälzle et al., 2017).

More than an arena for social cognition to play out, the social environment is home to resources important for human cognition (Gross & Medina-DeVilliers, 2020) and well-being (Coan et al., 2017;De Jaegher et al., 2010;Morawetz et al., 2021;Zaki & Williams, 2013). The quality of our social relationships and the support we receive from them has a large impact on our health (Beckes & Sbarra, 2022), including better prognoses related to disease (Cohen, 1988,2004;Cohen & Wills, 1985;Williams et al., 2018), increased antibody response to vaccines (Uchino et al., 2020), lower levels of inflammation (Uchino et al., 2018), reduced physiological impact of adversity (Wymbs et al., 2020), conserved energy in response to threats (Gonzalez et al., 2021;Saxbe et al., 2020), decreased risk of mortality (Berkman & Syme, 1979), and more (Chen & Miller, 2007;Hawkley & Cacioppo, 2007;House et al., 1988;Miller et al., 2009;Repetti et al., 2002;Seeman & McEwen, 1996;Uchino et al., 1996).

From an ecological neuroscientific perspective, social support can be seen as a resource, reliable access to which can dramatically alter health and behavioral outcomes. Indeed, ecological approaches to neuroscience have long investigated the link between social relationships and brain connectivity (Callaghan & Tottenham, 2016;Coan et al., 2017;DeCross et al., 2022;Gee et al., 2013;Gonzalez et al., 2015;McLaughlin et al., 2014;Morawetz et al., 2021;Sheridan & McLaughlin, 2014). While neuroscientists who study how brain organization and dynamics vary across individuals given the experiences during development, often have extensive data about individuals’ environments, brain analyses from an ecological perspective often focus on the connectivity between a few seed regions of interest (Callaghan and Tottenham, 2016;Gee et al., 2013;Gonzalez et al., 2015, but seeDeCross et al., 2022;Pelletier-Baldelli et al., 2022; for network approaches). Network neuroscience, on the other hand, provides tools for understanding environmentally influenced differential connectivity in the context of the entire brain. Although we are far from a thorough understanding of brain network reconfigurations in response to dynamic social environmental resources, integrating network and ecological perspectives serves to advance both fields.

Network techniques take advantage of whole-brain data by considering all interactions between all regions of the brain simultaneously. This holistic approach clarifies the contributions of regions of interest in the context of the whole brain. Importantly, whole-brain network approaches do not preclude using extensive data about the environment; indeed, in recent years, many studies have made connections between changes in brain networks over the course of development (Baum et al., 2017;Betzel et al., 2014;Dosenbach et al., 2010;Fair et al., 2007;Uddin et al., 2011) and variation in brain networks across different environmental experiences (Chan et al., 2018;Chen et al., 2016;Ellwood-Lowe et al., 2021;Rakesh, Zalesky, et al., 2021). There are many existing methodological approaches for linking brain network data to social environment (or behavioral/clinical/cognitive) data. Several studies have made this link using multivariate statistical approaches, including canonical correlation analysis (Xia et al., 2018), partial least squares (Mwilambwe-Tshilobo et al., 2019), and linear models (Ellwood-Lowe et al., 2021). By drawing from methods in network science as opposed to multivariate statistics, our approach makes no assumptions of linearity, involves no dimensionality reduction, and can quantify the extent to which brain-behavior associations are stable versus variable across different contexts (Mucha et al., 2010). Thus, more than just disentangling overlapping or distinct associations between brain networks and social environment (Rakesh, Zalesky, et al., 2021), network science methods can facilitate identification of brain systems that differentiate between social contexts.

Moreover, a variety of network tools (e.g., multilayer networks) exist to study brain connectivity over time or across contexts. When network neuroscientists have considered social information, the focus is often on social cognition and affect (Barrett & Satpute, 2013), social networks (Noonan et al., 2018), face perception (Cassidy et al., 2021), emotion regulation (Morawetz et al., 2020,2021), loneliness (Mwilambwe-Tshilobo et al., 2019,2023;Spreng et al., 2020), or socioeconomic status (Chan et al., 2018;Ellwood-Lowe et al., 2021). Understudied in this work, then, is a network neuroscientific approach tosocial support(but seeZajner et al. (2022)andSchurz et al. (2021)for recent examples of such an investigation). How do whole-brain patterns of connectivity vary across different experiences of perceived social support? Addressing such a question would contribute to ongoing work linking peoples’subjectiveperceptions to neural data, like fMRI (Mwilambwe-Tshilobo et al., 2019;Spreng et al., 2020), which, for social information, can yield different and additional insights than objective measures of social experience alone (Chen et al., 2004;Coan et al., 2017).

The present study builds on existing work linking the social environment to variation in brain network organization. Specifically, we use a well-established multilayer modeling framework (Mucha et al., 2010) to examine correlations of brain functional connectivity weights with a suite of social support measures. In a multilayer model, each layer is an individual network (e.g., a single functional connectivity network) which can be coupled to other layers to measure change in a network across the layers. Whereas past research has used the same framework to track variation in network structure across time (Bassett et al., 2011;Betzel et al., 2017;Finc et al., 2020), individuals (Betzel et al., 2019), or connectivity modalities (Puxeddu et al., 2022), we use this technique to track variation in brain-behavior correlation patterns. Our application of the multilayer framework presents a series of conceptual advances: (1) it preserves correlation profiles of individual social support measures, conferring a greater degree of interpretability, and (2) it identifies groups of brain regions whose correlation profiles are stableversusvariable across social support measures, that is, inflexibleversusflexible, providing a substantive contribution to the science of neural sensitivity to dimensions of social context. Put simply, we link the intricacies of social environment experience to the intricacies of brain network connectivity, without dimensionality reduction on either data source to make the connection.

Materials and Methods

2

Data

2.1

Our data included self-report social support data and resting-state fMRI (rsfMRI) data from the Human Connectome Project (HCP) (Van Essen et al., 2013). For the social data, we used the seven perceived social support measures from the NIH Toolbox included in HCP data:Friendship, Emotional Support, Instrumental Support, Hostility, Rejection, Loneliness, and Stress(Cyranowski et al., 2013). This set of measures was designed to capture various aspects of social support, companionship, and social distress. Specifically, Friendship refers to the availability of companions with whom to interact or affiliate; Emotional Support refers to the availability of people to listen to one’s problems with empathy, caring, and understanding; Instrumental Support refers to the perceived availability of people who can provide functional aid in completing daily tasks; Hostility captures the amount that others in one’s social environment argue, yell, or criticize them; Rejection measures how often others don’t like or make fun of one; Loneliness refers to the perception that one is lonely or socially isolated from others; and Stress captures general life stress (Cyranowski et al., 2013). All measures had eight items exceptLoneliness, which had five items. All items from all measures were rated on a five-point scale (ranging from “Never” to “Always”) and T-scored to fall between a range of 0 to 100, with 50 representing the population average. Within the present sample, the average score and standard deviation for each measure was: 51.11 $[eqn]$ 9.38 forFriendship; 50.39 $[eqn]$ 8.58 forEmotional Support, 49.79 $[eqn]$ 8.27 forInstrumental Support, 47.74 $[eqn]$ 8.27 forHostility, 51.91 $[eqn]$ 9.66 forRejection, 47.86 $[eqn]$ 8.38 forLoneliness, and 47.76 $[eqn]$ 7.87 forStress (seeSupplementary Fig. 1a, bfor score distributions and correlations between scores). Social measures were typically assessed on the same day as or within one day of rs-fMRI data collection.

The HCP dataset (Van Essen et al., 2013) consisted of resting-state functional magnetic resonance imaging (fMRI) scans (henceforth REST1 and REST2) from 100 unrelated adult subjects. These “100 Unrelated Subjects” were identified by the WU-Minn consortium as healthy, unrelated subjects with quality data, allowing researchers to forego controlling for family structure in analyses (WU-Minn, 2017). After excluding subjects based on data completeness and quality control (see Quality Control), the final subset utilized included 92 subjects (56% female, mean age = 29.29 $[eqn]$ 3.66, age range = 22–36, but seeSupplementary Fig. 2bfor data about how results generalize with different sample sizes). The study was approved by the Washington University Institutional Review Board, and informed consent was obtained from all subjects. A comprehensive description of the imaging parameters and image preprocessing can be found inGlasser et al. (2013). Images were collected on a 3 T Siemens Connectome Skyra with a 32-channel head coil. Subjects underwent two T1-weighted structural scans, which were averaged for each subject (TR = 2400 ms, TE = 2.14 ms, flip angle = 8°, 0.7 mm isotropic voxel resolution). Subjects underwent four resting-state fMRI scans over a two-day span, henceforth we refer to the rsfMRI scans from day 1 as REST1 and the rsfMRI scans from day 2 as REST2. (SeeSupplementary Fig. 2cfor comparisons of our findings across REST1 and REST2.). The fMRI data were acquired with a gradient-echo planar imaging sequence (TR = 720 ms, TE = 33.1 ms, flip angle = 52°, 2 mm isotropic voxel resolution, multiband factor = 8). Each resting-state run duration was 14:33 min, with eyes open and instructions to fixate on a cross.

Quality control

2.2

All preprocessed time series were visually inspected for visual artifact. Subject motion measurements during the fMRI scanning sessions were obtained from the HCP minimal preprocessing pipeline output directories (files: Movement_RelativeRMS.txt and eddy_unwarped_images.eddy_movement_rms). Across fMRI sessions and within an fMRI session, the mean and mean absolute deviation of the motion measurements were calculated, resulting in four summary motion measures per subject. Subjects without complete data or exceeding 1.5 times the inter-quartile range (in the adverse direction) of the measurement distribution for more than one of these summary measurements were excluded. We discarded frames $[eqn]$ 0.15 motion threshold. After these quality assurance steps, data from 92 subjects remained.

Image processing

2.3

Functional images of the HCP dataset were minimally preprocessed according to the description provided inGlasser et al. (2013). Briefly, T1w images were aligned to MNI space before undergoing FreeSurfer’s (version 5.3) cortical reconstruction workflow. fMRI images were corrected for gradient distortion, susceptibility distortion, and motion, and then aligned to the corresponding T1w with one spline interpolation step. This volume was further corrected for intensity bias and normalized to a mean of 10000. This volume was then projected to the32k_fs_LRmesh, excluding outliers, and aligned to a common space using a multi-modal surface registration (Robinson et al., 2014). The resultantciftifile for each HCP subject used in this study followed the file naming pattern: *REST{1,2}{LR,RL}_Atlas_MSMAll.dtseries.nii.

Parcellation

2.4

As HCP fMRI was provided in32k_fs_LRspace, these data could be parcellated based on the available Schaefer 400 parcellation (Schaefer et al., 2018) in theciftifile format. Additionally, a novel gradient-based subcortical parcellation was used to delineate nodes within the amygdala, hippocampus, thalamus, and striatum consisting of 41 regions per hemisphere (Scale III parcellation fromTian et al., 2020). Volumetric ROIs of subcortical and cerebellar (Diedrichsen et al., 2009) regions were defined in the MNI152 1 mm space, conforming to the cifti specification. We also compared coarser 200 cortical node and 100 cortical parcellations to this maximally fine-grained parcellation (seeSupplementary Fig. 2c).

Functional connectivity

2.5

Each preprocessed BOLD image was linearly detrended, band-pass filtered (0.008-0.08 Hz), confound regressed, and standardized using Nilearn’s signal.clean function, which removes confounds orthogonally to the temporal filters. The confound regression strategy included six motion estimates, mean signal from a white matter, cerebrospinal fluid, and whole-brain mask, derivatives of these previous nine regressors, and squares of these 18 terms. Spike regressors were not applied to the HCP data. The 36 parameter strategy (with and without spike regression) has been show to be a relatively effective option to reduce motion-related artifacts (Parkes et al., 2018). Following these preprocessing operations, the mean signal was taken at each node in the surface space. We report on global signal regressed data in the main text but provide comparisons to data without global signal regression inSupplementary Figure 1a.

Correlating functional connectivity and social support

2.6

For each resting scan, we computed functional connectivity using the Pearson correlation between all 482 nodes for each individual (seeFig. 1a–c). This gave us 92 matrices of dimension 482 by 482, whose entries $[eqn]$ represent the subject-specific edge weight between nodesiandj(seeFig. 1a). Each of the 482 nodes was assigned a functional system (and subsystem) byYeo et al. (2011);Schaefer et al. (2018), which we later use to describe the results at a meso-scale. We vectorized these subject-specific functional connectivity matrices to obtain a $[eqn]$ vector of edge weights. Subjects’ edge weight vectors were concatenated to create a $[eqn]$ edge weight matrix (seeFig. 1c). This allowed us to calculate the Spearman correlation between scores on each of the seven social support measures (seeFig. 1b) and all edge weights, resulting in a $[eqn]$ matrix.

A schematic of the method for computing the correlations between functional connectivity and social support. (a) shows the 92 (one per subject) functional connectivity matrices, where elementijindicates the Pearson correlation between nodesiandj. (b) represents all 92 subjects scores on seven multi-item measures of perceived social support:Friendship, Stress, Rejection, Hostility, Emotional Support, andInstrumental Support. The matrices in (a) are vectorized and stacked on top of each other in (b) to yield a 92 (number of subjects) by 115921 (number of edges on each functional connectivity matrix) matrix. We compute the correlation between social support measures (b) and edge weights (c), resulting in seven (one for each measure) correlation matrices (d), where elementijin matrixkgives the correlation between social support measurekand the edge weight between nodesiandj. Because modularity maximization treats negative weights as belonging to just one community, we transform each correlation matrix to a similarity matrix (e). We do this by computing the similarity between every pair of rows in the correlation matrix. Each of the seven similarity matrices becomes itself a layer in the multilayer model (f). The multilayer model is a 2D tensor (f) whose on-diagonal elements are the seven correlation matrices and whose off-diagonal elements are the coupling between layers. We perform community detection on this multilayer model, which yields a community assignment for every node for each layer (g). We compute for each node its flexibility, or how much it changes its community affiliation between layers (h).

Each of these vectorized correlations could be returned to the upper triangle of a $[eqn]$ matrix, resulting in seven such matrices, where the elementijof matrixkgave the correlation between the edge weight of nodesiandjand social support measurek(seeFig. 1d). These seven correlation matrices had a high correlation across REST1 and REST2 (seeSupplementary Fig. 1c), with and without global signal regression (seeSupplementary Fig. 2a), across sample sizes (seeSupplementary Fig. 2b), and across parcellations of cortex into 100, 200, and 400 regions (seeSupplementary Fig. 2c). The results of those analyses are reported inSupplementary Figure 1. The results presented in the remaining sections focus on correlation values averaged over both REST1 and REST2 scans based on data processed using a pipeline that includes global signal regression and a cortical parcellation of 400 nodes with 82 subcortical parcels.

Examining brain-behavior correlation patterns

2.7

Rather than examine the correlation patterns independently for each social support measure, we study them jointly, representing each as a layer in a multi-layer network. This approach has been used previously inBassett et al. (2011,2013),Finc et al. (2020), andBetzel et al. (2019), though to the the best of our knowledge this is the first time that the layers have corresponded to patterns of brain-behavior correlations, as opposed to functional connectivity. Multilayer networks provide a framework to study different versions of a network at once. In a multilayer model, each layer is an individual network (e.g., a single functional connectivity network) which can be linked to other layers to capture change in a network over time or in different contexts. This allows any network computation to be performed on all layers simultaneously, providing insight into how a network changes over time or across different contexts (Mucha et al., 2010).

We sought to quantify which patterns were consistent and which fluctuated across the social support measures. To do this, we used community detection. If we were to identify the community structure for each social support measure individually, there would be no principled way to compare structure across social support measures. We cannot know whether Community 1 for one measure is referring to the same set of nodes as Community 1 for another measure if we perform community detection on the layers independently, because the community detection algorithm involves some stochasticity. Furthermore, because the seven social support measures are strongly correlated with one another—some positively and some negatively—we want to take all measures into account when performing any analysis. Not doing so would provide redundant information. To accommodate this, we used multilayer modularity maximization (Newman & Girvan, 2004). Modularity maximization (both the single- and multi-layer versions) aims to identify groups of nodes whose connectivity to one another exceeds that of a chance model. Here, however, we aim to detect groups of brain regions with similar brain-behavior correlation profiles. We use a uniform null model and the Louvain algorithm to maximize the modularity function $[eqn]$ :

$$[eqn]$$

Nodes are linked to each other across layers through the coupling parameter $[eqn]$ . The value of $[eqn]$ influences the homogeneity of communities across layers (indicated by $[eqn]$ and $[eqn]$ ), such that small $[eqn]$ values emphasize layer-specific communities and large $[eqn]$ values identify communities shared across layers. $[eqn]$ and $[eqn]$ are the actual and expected weights of the link connecting nodes $[eqn]$ and $[eqn]$ . $[eqn]$ indicates to which cluster node $[eqn]$ belongs. The Kronecker delta function, $[eqn]$ takes a value of 1 when $[eqn]$ and a value of 0 otherwise. The parameter $[eqn]$ represents a spatial resolution weight that scales the influence of the null model. The optimization of $[eqn]$ returns for each layer a partition of the network into assortative communities (seeFig. 1g).

Multilayer model

2.7.1

Before we applied multilayer modularity maximization to the multilayer model, a transformation of the data was necessary. When a network includes negative weights, as is the case with correlation-based networks, modularity maximization will treat all negatively weighted items as part of one community. Since this could misrepresent patterns of correlations found in our data, we transform the seven correlation matrices into similarity matrices. We do this by computing the similarity between every pair of rows in a given correlation matrix, such that the entry $[eqn]$ corresponded to the Pearson correlation between the correlation profiles of regions $[eqn]$ and $[eqn]$ (seeFig. 1e). Each of the seven similarity matrices becomes itself a layer in the multilayer model, such that within a given layer we are given the relationships between scores on one social support measure and all edge weights for all subjects (seeFig. 1f). By using a multilayer approach, we can run community detection on all social support layers with varying amounts of coupling between layers and track which nodes change their community affiliation between which layers (seeFig. 1g). This particular formalization represents an innovation in applications of multilayer networks to neuroscience.

To build the multilayer model, we flatten our seven similarity matrices into a 2D tensor, such that each matrix (which is itself the similarity matrix of correlations between a social support measure and functional edge weights) is positioned on the diagonal of the tensor matrix (seeFig. 1f). Off-diagonal entries represent the coupling parameter $[eqn]$ . If $[eqn]$ is 0, each layer is independent in any calculation performed. As $[eqn]$ approaches 1, all the layers are increasingly treated as a single object where the individual layers are indistinguishable. By varying $[eqn]$ , then, we see how the community structure of the individual layers differs from each other and from the community structure of the layers taken as a whole. We tested a range of nine logarithmically spaced $[eqn]$ s from 0.0001 to 1.

Community detection

2.7.2

Since we are interested in the community structure of this multilayer model, we ran modularity maximization (Newman & Girvan, 2004) with a uniform null model and the Louvain algorithm to maximize the modularity function $[eqn]$ . By tuning the resolution parameter $[eqn]$ , we see differences in communities according to how large their average correlations are. Larger $[eqn]$ s result in communities for which the average correlation between members is high, yielding numerous small communities. Smaller $[eqn]$ s give communities for which the average correlation between members is low, yielding fewer and bigger communities. We tested a range of 19 linearly spaced $[eqn]$ s from 0 to 0.9. Using 100 iterations of the modularity maximization algorithm, we also computed consensus communities (seeFig. 4a).

Flexibility

2.8

Using the community assignments resulting from the community detection algorithm, we can calculate how much a given node changes its community affiliation between layers, or the flexibility $[eqn]$ of node $[eqn]$ :

$$[eqn]$$

where $[eqn]$ is the number of times node $[eqn]$ changes its community affiliation, and $[eqn]$ is the total number of changes possible (seeFig. 1h). When the coupling parameter $[eqn]$ is large, in general most nodes will have lower flexibility because the layers are so highly coupled that the amount a node can change is constrained. Similarly, when $[eqn]$ is small, most nodes will have higher flexibility.

Results

3

Our data included self-report social support data and resting-state fMRI (rsfMRI) data from 100 unrelated adult subjects from the Human Connectome Project (HCP;Van Essen et al., 2013; further description of sample selection inMethodssection). For the social data, we used the seven perceived social support measures from the NIH Toolbox included in HCP data:Friendship, Emotional Support, Instrumental Support, Hostility, Rejection, Loneliness,andStress(Cyranowski et al., 2013). For the brain data, we used four resting-state functional magnetic resonance imaging (fMRI) scans. The social and brain data were collected within one day of each other.

Permutation test of correlations between edge weights and perceived social support

3.1

To test which Spearman correlations between edge weights and social support measures had a greater magnitude than would be expected by chance, we used permutation testing. Here, edges refer to the average magnitude of co-fluctuations between nodes (brain regions). Edge weights are calculated by taking the Pearson correlation of two nodes’ time series for every pair of nodes. Computing Spearman correlations between edge weights and social support measures thus quantifies the association between a given social support measure and the co-fluctuation of a pair of nodes. We shuffled social support scores across subjects (i.e., we shuffled the order of the rows of the matrix inFig. 1b) 1000 times and recalculated the correlation between edge weights and social support scores each time. We counted the number of instances for which the observed correlation between edge weight and social support score had a greater magnitude than the permuted correlation. To obtain a $[eqn]$ -value for each edge weight-social support measure correlation, we divided this count by the number of permutations (1000) and used false-discovery rate of 0.05 to correct for multiple comparisons. This gave an adjusted significance criterion of $[eqn]$ .Figure 2shows which elements and what proportion of each canonical system survived permutation testing.

System-level correlations between edge weights and social support measures. We shuffled social support scores across subjects and recalculated the correlations with edge weights 1000 times. For each measure, we summarize the results of this permutation testing by showing the log proportions of edges of each canonical subsystem (e.g., ContA, ContB) that are significantly correlated with the given social support measure using a false discovery rate of 0.05 (adjustedpcrit=0.046).

Analyses of the parameter space of the multilayer model

3.2

How do system- and layer-specific behaviors change as community resolution $[eqn]$ and inter-layer coupling $[eqn]$ change? To answer this, we analyzed the parameter space of our multilayer model using several methods: applying principal components analysis to a measure of community flexibility to identify modes of community variation across parameter values, by tracking variability in regional flexibility patterns, by studying the partition landscape degeneracy (a measure of how similar detected partitions are to one another at a given point in parameter space), and by comparing select partition “seeds” to the rest of the parameter space. These analyses motivate our selection of points in parameter space for further study.

Principal components analysis

3.2.1

First, with the aim of arriving at a holistic understanding of how varying the parameter values is related to changes in the relationship between brain networks and social support, we performed a principal components analysis (PCA) on the flexibility of all nodes in all layers across all values in parameter space (seeSupplementary Fig. 3for amount of variance explained) (Betzel et al., 2019).Figure 3ashows the magnitude of expression of the first principal component, which explains about 35% of the variance, by each point in parameter space (see alsoSupplementary Fig. 6for the magnitudes of expression of the second and third principal components in parameter space and in brain space). Our results suggest that strong positive coefficients were associated with small values of the community resolution parameter $[eqn]$ , and strong negative coefficients with larger $[eqn]$ s.Figure 3gprojects this into brain space for both the cortex and subcortex. The point in parameter space whose coefficient for principal component 1 had the largest magnitude was at $[eqn]$ . At these parameter values, there were 14 communities with an average size of $[eqn]$ nodes.

Analyses of parameter space. We performed PCA on the entire parameter space; (a) shows the magnitude of expression of the first principal component. The point in parameter space for which the magnitude of expression of the first principal component is greatest isγ=0.15, ω=0.001, indicated by a black outline. (b) Shows the variability in node flexibility across parameter space computed across individual partitions at a given point in parameter space. The point in parameter space that maximizes variability in flexibility isγ=0.55, ω=0.0316, indicated by a black outline. In (c), for every point in parameter space, the similarity between every partition was computed using the adjusted Rand Index then averaged for each point across layers. For this calculation, the global maximum similarity is achieved by five points in parameter space:γ=0.9, ω=0.1; γ=0.85, ω=0.3162; γ=0.9, ω=0.3162; γ=0.85, ω=1andγ=0.9, ω=1, indicated by a black outline. We also calculated the similarity of each individual partition to the consensus communities at that point in parameter space and took the average across partitions and layers (d). For this calculation, the global maximum was reached by a set of five points in parameter space:γ=0.9, ω=0.1; γ=0.85, ω=0.3162; γ=0.9, ω=0.3162; γ=0.85, ω=1andγ=0.9, ω=1, indicated by a black outline. Panels (e) and (f) show seed-based analyses, where we computed the similarity of all partitions for each point in parameter space to the seed. For (e), we used as a seedγ=0.55, ω=0.0316, which maximizes the variability in flexibility. For (f), we used as a seedγ=0.35, ω=1, which is the point in parameter space whose communities are most similar to the Yeo systems. (g) Shows the extent to which different regions of the cortex and subcortex express the mode of flexibility present in the first principal component.

Variability in flexibility

3.2.2

Second, we computed for each point in parameter space the variability of flexibility for all nodes across all layers (seeFig. 3c). Previous work (Bassett et al., 2011) has selected the parameter values which maximize variability in flexibility for further analyses. In our case, this point was at $[eqn]$ . With this pair of parameter values, we detected 79 communities with an average size of $[eqn]$ nodes.

Partition landscape degeneracy

3.2.3

Given that we ran 100 iterations of community detection, how similar was each resulting partition to other runs with the same parameter values? That is, for which points in parameter space did the community detection algorithm consistently find the same communities across runs? We approached this two ways. First, for a given pair of parameter values, we calculated the similarity (using the adjusted Rand Index; ARI) between every pair of partitions then took the average (seeFig. 3d). In general, a larger ARI for a pair of partitions indicates that they are relatively similar to each other; smaller values mean greater dissimilarity between partitions. Second, for every pair of parameter values, we calculated the similarity of each partition to the consensus communities (seeFig. 3e). This second approach tells us more about how representative the consensus communities are for a given pair of parameter values. For both of these approaches, similarity was maximized at the highest values of the community resolution parameter $[eqn]$ (i.e., many small but internally very similar communities) and the coupling parameter $[eqn]$ (i.e., maximal inter-layer coupling).

Seed-based similarity

3.2.4

Next, we wanted to quantify how similar a community partition at one point in parameter space was to all other points. In other words, how representative of all of parameter space is the one pair of $[eqn]$ and $[eqn]$ values that maximize variability in flexibility? We approached this two ways. First, we selected the point in parameter space whose consensus communities maximized variability in flexibility ( $[eqn]$ ). Then, we used this point as a seed and computed the adjusted Rand Index to the seed for the consensus communities for all other points in parameter space (seeFig. 3e). Additionally, to link our results more clearly to canonical brain systems in network neuroscience, we also used as a seed the point in parameter space that maximized similarity (using ARI) to the Yeo systems (Yeo et al., 2011). This point was found at $[eqn]$ , when layers are maximally coupled (seeFig. 3f).

Selecting points in parameter space for further analysis

3.2.5

In principle, any of these approaches could be used to select a point or points in parameter space for further analysis. In general, we saw more variation across values of the community resolution parameter $[eqn]$ than the coupling parameter $[eqn]$ , in line with other studies (Betzel et al., 2019). Coefficients for the first principal component suggested different patterns are to be found for smaller $[eqn]$ s versus larger $[eqn]$ s. To select specific smaller and larger $[eqn]$ values, we opted for the two points in parameter space whose coefficients on the first and second principal components have the greatest magnitude. The parameter values for which the magnitude of the coefficient of the first principal component is maximized is $[eqn]$ . For principal component two these values are $[eqn]$ . By using this approach, we have identified points are maximally representative of the dominant modes of variation in parameter space.

Multilayer community structure

3.3

What does the multilayer community structure look like at the points in parameter space we have selected? When the magnitude of the coefficient for the first principal component is maximized ( $[eqn]$ ), we find a total of 14 communities, only four of which were present in every layer (seeFig. 4). Across the seven layers of perceived social support, a few patterns were consistent. Somatomotor and adjacent nodes tended to be in the same community in dark blue. For most layers, visual nodes were also in this same community. Frontal, temporal, medial, and subcortical nodes had more variety across layers.

Multilayer consensus community structure. Panel (a) shows the consensus community structure across social support measures at the point in parameter space whose coefficient for principal component 1 has the largest magnitude (γ=0.15, ω=0.001). Colors indicate communities. There are 14 communities total, but only four communities are present across all layers (shown in dark blue, light blue, teal, and red). Panel (b) shows the number of communities in each layer according to the consensus partition. Panel (c) gives the average community size for each layer according to the consensus partition. The MNI coordinates for the five panels of subcortical communities are, from left to right: x = -23, x = 10, z = -3, x = 13, x = 27.

When the magnitude of the coefficient of the second principal coefficient is maximized ( $[eqn]$ ), there were 25 communities, all of which were present in every layer (seeSupplementary Fig. 5). Somatomotor and visual nodes again tended to be in the same community, but with less consistency. While there were many communities represented by the subcortex, a single community (yellow) was most prominent across all layers, unlike when $[eqn]$ . Again, frontal, temporal, and medial nodes were associated with a variety of communities across layers. To understand in more depth which nodes change their community affiliation between which layers, we looked at flexibility.

Flexibility

3.4

Flexibility across the brain and across canonical systems

3.4.1

Which regions of the brain switch their community affiliation more (i.e., are more flexible)? Do flexible nodes tend to belong to the same canonical system? To address these questions, we examined the flexibility (i.e., the amount a node switches communities between layers) of nodes across the brain, across canonical systems, and across social support layers. We used theYeo et al. (2011);Schaefer et al. (2018)parcellation to define our (sub)systems. We found that, at the point in parameter space where the magnitude of the coefficient for the first principal component is maximized ( $[eqn]$ ), flexibility was widely distributed across the brain (seeFig. 5a) but unevenly distributed across the nodes of canonical systems ( $[eqn]$ using a Kruskal-Wallis rank-sum test; seeFig. 5b). The system with the most flexible nodes was ContC (median flexibility = 0.849), and the system with the least flexible nodes was SomMotA (median flexibility = 0.63).

Flexibility. Panels (a), (b), (c), and (d) represent data from the point in parameter space (γ=0.15, ω=0.001) with the largest loading on the first principal component. Panel (a) shows the flexibility of cortical and subcortical regions. Panel (b) shows the distribution of flexibility across canonical systems, as defined byYeo et al. (2011);Schaefer et al. (2018), ordered by median flexibility. We note that because trends largely hold for both the system and subsystem scale (e.g., Control and Subcortex are more flexible than other systems, even though Limbic and DefaultC are more flexible than many Control or Subcortex subsystems), we present the results for subsystems here and then proceed with system-level reporting for the remainder of our results. Panel (c) shows the proportion of all nodes that change community affiliation between each pair of the seven layers: Friendship, Stress, Rejection, Hostility, Emotional Support, Loneliness, and Instrumental Support. Panel (d) shows how similar all systems are to each other in terms of the proportion of nodes changing community between each pair of layers. Panel (e) plots the correlation between each system and the global pattern of community change between layers. The MNI coordinates for the five panels of subcortical flexibility are, from left to right: x = -23, x = 10, z = -3, x = 13, x = 27.

To compare the flexibility of nodes in different systems, we used permutation testing to construct a null model. We permuted all nodes’ flexibility 10000 times; for each permutation we compared the flexibility of the nodes for all pairs of systems. To compute $[eqn]$ -values, we counted the number of instances in which the direction of difference in flexibility between a pair of systems was greater in the observed data than in the permuted data, then divided this by the number of permutations (10000). ContC emerged as more flexible than DefaultA, DefaultB, DorsAttnA, DorsAttnB, SomMotA, SomMotB, TempPar, VisCent, and VisPeri (all $[eqn]$ s $[eqn]$ 0.0015). SomMotA was less flexible than all systems except SomMotB, TempPar, and Amygdala (all $[eqn]$ s $[eqn]$ 0.0037).

The Limbic system had the second most flexible nodes after ContC (median flexibility = 0.845), and was more flexible than DefaultA, DefaultB, DorsAttnA, SomMotA, SomMotB, TempPar, VisCent, and VisPeri (all $[eqn]$ s $[eqn]$ 0.005). Subcortical nuclei varied in their flexibility, with median flexibility ranging from 0.845 in the Cerebellum to 0.722 in the Amygdala. Ultimately, none of the subcortical nuclei were statistically distinguishable from one another. The only coarse-grained system that had distinguishable fine-grained systems was Default Mode: DefaultC was more flexible than DefaultA ( $[eqn]$ = 0.004).

The relative lack of flexibility in somatomotor, temporoparietal, and visual systems suggests that these systems are largely indifferent to the different measures used here. Given that our neuroimaging data were collected during isolated rest, these results do not preclude the importance of sensory regions during live social interaction. The systems with the most flexible nodes, on the other hand, (e.g., ContC and Limbic, and to a lesser extent Cerebellum and DefaultC) may show patterns of activity that are more sensitive to the differences between the seven measures of perceived social support used here. In other words, given the connectivity of these more flexible systems, the measures of perceived social support are more differentiable. It is possible but unlikely that this is explained by variability (seeSupplementary Fig. 8). Caution is warranted in further interpretation of these results, given the low test-retest reliability of, for example, subcortical nodes (Noble et al., 2017).

Flexibility across social support measures

3.4.2

Are there particular social support layers that are substantially different from the rest, forcing nodes to switch their community affiliation and therefore driving flexibility higher? To address this question, we looked at the proportion of nodes that switched their community affiliation between each pair of layers (seeFig. 5cfor switches between layers for all nodes and seeSupplementary Fig. 4for switches between layers broken down by canonical system). To construct a null model, for each partition and for each node we permuted community affiliations across layers 1000 times. We then computed the proportion of nodes that changed communities between every pair of layers and took the average across partitions. To obtain a $[eqn]$ -value, we divided this average by the number of permutations. The proportion of nodes that switched their communities between theFriendshiplayer and every other layer was significantly larger in the observed data than in the permuted data (all $[eqn]$ s $[eqn]$ 0.0001). No other layer had significantly different proportions of nodes changing foreverycomparison, butInstrumental Supporthad the next most. The proportions of nodes changing communities were significantly larger forInstrumental Supportthan all layersexcept Hostility(all $[eqn]$ s $[eqn]$ 0.03). To further examine whetherFriendshipdrove patterns of flexibility, we recalculated flexibility by systematically excluding each layer in turn then calculating the correlation between recalculated and original flexibility values (seeSupplementary Fig. 7). WhileFlexibility, when excluded, had the lowest correlation, all correlations were $[eqn]$ (and all $[eqn]$ s $[eqn]$ ). Taken together, this suggests that communities, and by extension the roles nodes play, in theFriendshiplayer were especially different from all the other layers, even if all layers contributed substantially to flexibility. Communities in theInstrumental Supportlayer also tended to be different from most other layers.

How similar are systems to this global pattern, in which theFriendshipandInstrumental Supportlayers contain relatively distinct communities? To determine this, we computed the proportions of nodes changing communities between each pair of layers both (1) for the whole brain and (2) for each system. We then correlated these lists of proportions. All systems except Somatomotor and Temporoparietal (which are both generally inflexible; seeFig. 5b) had significant correlations with the global pattern (all $[eqn]$ s $[eqn]$ 0.56, all $[eqn]$ s $[eqn]$ 0.0009; seeFig. 5e). DorsAttn, Control, and Vis were the systems most similar to the global pattern (all $[eqn]$ s $[eqn]$ 0.79); they were also most similar to each other (all $[eqn]$ s $[eqn]$ 0.81; seeFig. 5d), even though they varied greatly in the flexibility of their nodes (seeFig. 5b). Temporoparietal, on the other hand, was dissimilar from all other systems ( $[eqn]$ with Default = 0.47, all other $[eqn]$ s $[eqn]$ 0.22). Similarity to the global pattern of Friendship driving flexibility was not observed only in more flexible systems, like Control; less flexible systems, like Vis, also showed this pattern.

Discussion

4

We set out to link different experiences in the social environment to variation in brain network organization. We use the established methods of community detection (Newman & Girvan, 2004) and multilayer modeling (Mucha et al., 2010) in multiple innovative ways to make this connection. We highlight two findings that contribute to research on the neural sensitivity to social support. First, the subcortex and control systems are especially sensitive to the different constructs of perceived social support. Network nodes in these systems are highly flexible; their community affiliations differ at high rates across social support measures. Nodes in the somatomotor, temporoparietal, and visual systems, on the other hand, were relatively stable in their community affiliations. Second, theFriendshiplayer, and to a lesser extentInstrumental Support, drove flexibility; higher proportions of nodes across all systems switched their community affiliations betweenFriendship(andInstrumental Support) and all other layers.

When examining the stability of our results across processing pipelines, we find that correlational patterns between brain network edge weights and scores on seven measures of social support are consistent across resting fMRI scans and with and without global signal regression. While many edges do not survive permutation testing, notably many in the subcortex do. We examine the parameter space of the multilayer model in detail according to variability in flexibility; according to similarity of partitions to each other, to consensus communities, and to canonical systems; and according to modes of variation based on PCA. In general, there is more variation across the resolution parameter ( $[eqn]$ ), which tunes the size and number of communities, than the coupling parameter ( $[eqn]$ ), which tunes the strength of the relationship between layers, in line with previous network neuroscience studies that use multilayer models (Betzel et al., 2019).

We zoom into two points in parameter space: the parameter value pairs for which the magnitudes of the coefficients of the first and second principal components are maximized. The point for PC1 yield 14 total communities, while the point for PC2 yields 25 total communities. In both cases, Somatomotor, Default, and Visual nodes mostly have low flexibility while Control and Subcortex (and to a lesser extent Attention and Limbic) nodes have a wider distribution to their flexibility values. The increased flexibility of Control and Subcortex nodes suggests these systems may be better able to differentiate between the different social contexts (represented by the seven measures used here). Given that previous studies have found these systems to vary between groups from different socioeconomic statuses (Ellwood-Lowe et al., 2021;Rakesh, Seguin, et al., 2021), they may be key for supporting adaptation to social environments and the behaviors relevant for these environments. Indeed, myriad studies have noted the sensitivity of Control (Coan et al., 2017) and Subcortex (Morawetz et al., 2021;Pozzi et al., 2021;Rudolph et al., 2021) systems to social support. Given that subcortical nodes tend to have low test-retest reliability (Noble et al., 2017), caution is warranted in interpreting these results.

When we broke down flexibility across the brain and canonical systems (Yeo et al., 2011), we found that ContC had the most flexible nodes while SomMotA had the least flexible. Greater flexibility, in our case, suggests that the role(s) of nodes in a system differentiate between the measures of social support. Lower flexibility has two possible interpretations: (1) the nodes of a system do not distinguish between measures of social support because their activity is not driven by or correlated with social support or (2) the nodes of a system do not distinguish between measures of social support because their activity is similarly correlated with each measure. Disambiguating these two explanations might involve cross-referencing flexibility with system-level activation in social tasks. Given that previous studies have linked loneliness and social support with the Default network (Morawetz et al., 2021;Mwilambwe-Tshilobo et al., 2019;Pozzi et al., 2021,2023;Rudolph et al., 2021;Spreng et al., 2020;Xiao et al., 2021;Zajner et al., 2022), it might be the case that interpretation 2 applies to inflexible Default nodes. Interpretation 1 might then apply to inflexible Somatomotor and Visual nodes since we would expect these brain regions to have a weaker relationship with social support (but seeRakesh, Seguin, et al., 2021;Rakesh, Zalesky, et al., 2021for associations between somatomotor systems and socioeconomic status). Adjudication between these explanations should also consider the limited visual and sensory input subjects receive during rest. Lack of activation of sensory regions during rest may not preclude the importance of these regions during real life social interactions.

Just as we broke down flexibility across brain systems, we also broke down flexibility across the social support measures. We find that, for principal component 1, node roles (operationalized here as community affiliation) are substantially different in theFriendshiplayer beyond what would be predicted by measure dissimilarities. The same is true to a lesser extent for theInstrumental Supportlayer. Put another way, the associations between these measures and brain network edge weights are especially different from other layers, beyond the differences between the measures themselves. Given that social support influences bodily systems beyond the brain (Holt-Lunstad et al., 2008;Seeman & McEwen, 1996;Uchino et al., 1996), it is important to consider whether theFriendshiplayer would have comparably distinct associations between physiological and social support measures. Because integrated brain-body data is a relatively new development, it is difficult to comment with confidence on the generality of our results to other physiological measures, but this avenue of questioning will be critical for future work.

Notably, our layer-based analysis does not yield any results showing different patterns between layers indicating “positive” social experiences (e.g.,Friendship, Instrumental Support, Emotional Support) and “negative” social experiences (e.g.,Stress, Hostility, Rejection, Loneliness), even though subjects’ scores alone tend to cluster together in this way. Several neuroecological studies have proposed different dimensions of environmental stressors (McLaughlin & Sheridan, 2016;Sheridan & McLaughlin, 2014), including deprivation versus threat and social versus economic versus neighborhood. The indifference of community changes to “positive” versus “negative” experiences suggests that valence may not be an especially meaningful dimension in the context of these social support measures for the brain. A more salient distinction, according to previous literature, may be between emotional support and instrumental support (Morelli et al., 2015). While both measures capture aspects of social support, emotional support speaks more to empathy and instrumental support more to resource provision. This difference translates to varying risks for all-cause mortality (Holt-Lunstad et al., 2017;Uchino et al., 2018), varying trajectories with age (Mwilambwe-Tshilobo et al., 2023), and varying associations with brain organization and development (Morawetz et al., 2021;Mwilambwe-Tshilobo et al., 2019;Spreng et al., 2020;Xiao et al., 2021).

While the importance of control and subcortical nodes for social processes aligns with work in social neuroscience and neuroecology (Cacioppo et al., 2007;Coan et al., 2017;Morawetz et al., 2021;Pozzi et al., 2021;Rudolph et al., 2021;Steinberg, 2017;Xin & Lei, 2015), our approach provides rich detail that traditional region of interest analyses cannot. In particular, multilayer modeling allows us to take all social support variables into consideration, get detailed information about all nodes in the system, and make direct comparisons about nodes and systems across social support variables (seeRakesh, Zalesky, et al., 2021for a similar approach with linear models to socioeconomic status). This detailed information can be connected to canonical brain systems and social support, making the multilayer approach a bridge between subjective social measures and brain functional connectivity. To clarify what exactly the role of control and subcortical nodes might be, it is important to consider how the observed patterns of flexibility change across individuals according to their experiences in the social environment. That is, how do patterns of flexibility for control and subcortical nodes differ between individuals with consistent high quality social support and individuals who experience more rejection, hostility, loneliness, and stress in their social environment? More than that, how do different patterns of brain network organization relate to well-being within a particular environment? Several researchers have begun addressing these questions using measures of socioeconomic status (Ellwood-Lowe et al., 2021;Ramphal et al., 2020).

Our method offers innovation on commonplace network neuroscientific tools (Bassett et al., 2011;Vaiana & Muldoon, 2020) to provide details about associations between social context and whole-brain network organization. In this way, it offers more information than traditional neuroecological approaches, which often focus on particular brain regions in linear models. By drawing from methods in network science as opposed to multivariate statistics, our approach makes no assumptions of linearity, involves no dimensionality reduction, and can quantify the extent to which brain-behavior associations are stable versus variable across different contexts (Mucha et al., 2010). Moreover, this methodological innovation can be applied toanyset of measures, whether they are cognitive, behavioral, or demographic. We want to highlight that our layer-based analysis (which identified that more nodes change community affiliation betweenFriendshipand all other layers) is only feasible when there is all-to-all coupling between layers. In other words, applications of multilayer modeling to study time-varying functional connectivity cannot identify layer differentiation, because layers are often only coupled to adjacent layers. This layer-based analysis is thus a methodological innovation that helps shed light on how the role of brain systems might vary according to (in our case, social) context.

It is worth noting that our particular multilayer implementation, and by extension our layer-based analysis, can be applied to any set of interrelated measures. Our application of multilayer modeling to perceived social support is motivated by work in neuroecology (Coan et al., 2017;DeCross et al., 2022;Gonzalez et al., 2015;Morawetz et al., 2020;Pelletier-Baldelli et al., 2022), but the implementation’s agnosticism toward the measures used speaks to its generality: it is a powerful tool for a wide variety of research questions. This approach can be used to investigate how brain connectivity varies across experiences in the physical or economic environment or across different behavioral or cognitive measures.

Limitations and future directions

4.1

Within a developmental neuroecological framework, variation across individuals is meaningful and can point to important differences in environmental experience (Ellwood-Lowe et al., 2021;Ramphal et al., 2020). Critically, this study does not investigate individual or even group differences. By using group-averaged functional connectivity, we cannot speak directly to adaptations, individual differences, or precise neuroendophenotypes. Before such an approach can be taken, it was imperative to first show that network science tools could contribute to our understanding of neuroendophenotypes beyond multivariate statistics. That our results align with neuroecological findings (e.g.Coan et al., 2017;Gonzalez et al., 2015,2019;Ellwood-Lowe et al., 2021;Rakesh, Zalesky, et al., 2021;Ramphal et al., 2020) suggests network tools stand to add something to this kind of investigation. Therefore, future studies should explore how individual variation in brain networks relates to variation in experiences in the environment and how brain-environment and brain-behavior associations are related to well-being (seeEllwood-Lowe et al., 2021for such approaches using measures of socioeconomic status;Ramphal et al., 2020).

Because the HCP dataset includes only adults, we were limited in the developmental conclusions we can draw. Bigger effect sizes would be expected in adolescent populations in particular, since adolescent represents a critical period of rapid and high-stakes adaptation to the (especially social) environment (DeCross et al., 2022;Gonzalez et al., 2015;Pelletier-Baldelli et al., 2022;Rudolph et al., 2021). Applying methods similar to ours to a dataset like that of Adolescent Brain Cognitive Development study may be particularly illuminating.

Conclusions

5

The present study offers theoretical and methodological innovations, both of which establish a link between the fields of neuroecology and network neuroscience. Our theoretical contribution is motivated by neuroecology. Neuroecological perspectives hold that individual differences in brain features are meaningful and are often related to differences in experience in the environment over the course of development (Gottlieb, 1991). While there are myriad environmental measures that could be linked to brains, we focus on the social environment. Although network neuroscientists have explored links between the social environment and brain networks before (Barrett & Satpute, 2013;Cassidy et al., 2021;Chan et al., 2018;Ellwood-Lowe et al., 2021;Krendl & Betzel, 2022;Morawetz et al., 2020,2021;Mwilambwe-Tshilobo et al., 2019;Noonan et al., 2018;Spreng et al., 2020), the relationship between social support specifically and brain networks has been heretofore untouched. Our study addresses this gap explicitly and establishes a theoretical connection between neuroecology and network neuroscience by identifying how brain networks change in relation to social contexts.

Our methodological contribution is grounded in network science. Studying the relationship between social support and brain networks involves integrating rich environmental data with whole-brain data. In a way, this is structurally identical to the many brain-behavior association studies. We diverge from this resemblance by leveraging tools from network science to address neuroecological questions. We employ multilayer modeling to identify stability and variation in network community structure across seven measures of perceived social support. Because we construct the multilayer model to have all-to-all coupling, we are able to make pairwise comparisons between all layers to see how brain-environment relationships change across social support contexts. This methodology and our results align with and build on existing work in both network neuroscience (Bassett et al., 2011;Vaiana & Muldoon, 2020) and social neuroecology (Coan et al., 2017;Gonzalez et al., 2019), facilitating a rich connection between the two fields.

There are two key takeaways from this work: (1) dimensions of social environment are associated with variation in brain networks and (2) the methods used in this paper can be used to understand the interplay between brain and environment, behavior, cognition, or disease. This work serves as a foundation for future research at the intersection of network neuroscience and social neuroecology, or other fields.

Supplementary Material

Supplementary Material