Sex differences in multilayer functional network topology over the course of aging in 37543 UK Biobank participants

Although the human life span has increased from 50 to 80 years of age in the past two centuries, this has not been matched by an improvement in health span (Crimmins, 2015). In fact, age is one of the major risk factors for debilitating conditions such as cardiovascular and neurodegenerative diseases, having a large societal and economic impact (Dhingra & Vasan, 2012; Hou et al., 2019). However, not all individuals age in the same way. In particular, sex seems to be responsible for a substantial interindividual variability during aging, with women displaying a higher probability of developing certain age-related disorders such as Alzheimer’s disease (Mazure & Swendsen, 2016) and multiple sclerosis (Golden & Voskuhl, 2017), whereas men are more likely to develop Parkinson’s disease (K. M. Smith & Dahodwala, 2014). These differences in vulnerability to distinct diseases suggest that men and women have a distinct underlying brain network organization that might predispose them to develop specific pathological processes.

The functional network organization of the brain can be assessed using the correlations of spontaneous fluctuations in brain activity across brain regions by measuring the blood oxygen level–dependent signals on resting-state functional magnetic resonance imaging (rs-fMRI) (B. Biswal, Zerrin Yetkin, Haughton, & Hyde, 1995). Using this technique, several studies have identified highly reproducible resting-state networks in the brain such as the sensorimotor, dorsal attention, and default-mode networks, which play an important role in motor and cognitive functions (van den Heuvel & Pol, 2010). The communication between these networks is particularly important for brain function and has been shown to change during the course of aging (Betzel et al., 2014; Sala-Llonch, Bartrés-Faz, & Junqué, 2015), with older individuals showing a loss of anticorrelations (negative connections) (Keller et al., 2015; Spreng, Stevens, Viviano, & Schacter, 2016) and increases in positive correlations (positive connections) (Damoiseaux, 2017; Ferreira et al., 2016) between resting-state networks. These changes reflect the tendency of older individuals to over-recruit functional networks, needing to activate more brain networks than younger individuals, thus decreasing functional specialization and spending more neural resources (Goldstone et al., 2016).

However, the impact of sex on the communication between functional brain networks during aging is still not well understood (Ritchie et al., 2018; Scheinost et al., 2015), in part due to the small number of participants included in previous studies and their limited statistical power (Ritchie et al., 2018; Ruigrok et al., 2014). Assessing sex differences in functional connectivity is important for several reasons. For example, because functional networks are closely associated with cognitive and sensorimotor functions (van den Heuvel & Pol, 2010), understanding how their communication deteriorates with aging might provide important clues on why men and women are vulnerable to different diseases (Golden & Voskuhl, 2017; Mazure & Swendsen, 2016; K. M. Smith & Dahodwala, 2014) and why they show differences in other important health aspects such as brain structure, cardiovascular risk factors, and cognitive function (McCarrey, An, Kitner-Triolo, Ferrucci, & Resnick, 2016; Ramirez & Sullivan, 2018; Ritchie et al., 2018; Sachdev, Parslow, Wen, Anstey, & Easteal, 2009; Weiss, Kemmler, Deisenhammer, Fleischhacker, & Delazer, 2003).

From a methodological point of view, studies analyzing functional connectivity have mainly focused on positive connections (Chan, Park, Savalia, Petersen, & Wig, 2014; Tomasi & Volkow, 2012a). While this approach is more straightforward to assess the organization or topology of brain networks (Fornito, Zalesky, & Breakspear, 2013), negative connections are commonly found between brain networks or areas and seem to play an important role in brain communication (Fox et al., 2005; Hampson, Driesen, Roth, Gore, & Constable, 2010) and cognition (Barber, Caffo, Pekar, & Mostofsky, 2013). Since the negative connections or anticorrelations carry behaviorally relevant information (Barber et al., 2013), an integrative approach that incorporates information from negative correlations as well as positive correlations may reveal unique insights on sex differences throughout aging.

In this study, we developed this approach by combining the positive and negative functional connections between networks as separate layers in a complex multilayer network. We demonstrated that the multilayer measures provide novel insights about the impact of age and sex on the different functional connectivity trajectories in men and women between 47 and 79 years old in a large cross-sectional cohort of 37,543 individuals. Moreover, together with other measures of network organization, they were associated with structural brain imaging markers, cardiovascular risk factors, and cognitive functions, which typically differ between men and women (McCarrey et al., 2016; Ramirez & Sullivan, 2018; Ritchie et al., 2018; Sachdev et al., 2009; Weiss et al., 2003), as well as genes involved in physiological processes associated with aging. These findings open new avenues for the study of functional brain connectivity in aging by using multilayer network measures.

We included 19,975 women and 17,568 men with resting-state functional MRI from the UK Brain Biobank cohort (Miller et al., 2016) (ec12Methods, section ec13Participants and Supporting Information Figure S1). A subset of these individuals had available structural brain imaging data such as T1-weighted imaging and underwent a comprehensive assessment of cognitive functions and cardiovascular risk factors (ec12Methods, sections ec17Structural Imaging Preprocessing, ec19Cardiovascular Risk Factors, and ec20Cognitive Tests). Using permutation testing to compare the demographic characteristics between women and men, we found that men showed higher scores than women in the executive cognitive domain during middle ages and in the visuospatial cognitive domain across all ages (Supporting Information Figure S2C and D). In line with previous research showing that men are more susceptible to cardiovascular problems during middle adulthood (Anand et al., 2008; Gillis & Sullivan, 2016), we observed a significantly greater prevalence of high blood pressure, heart attack and white matter hyperintensities in men between the ages of 51 and 76 compared to women (Supporting Information Figure S2E, F, K). Finally, men had larger subcortical volumes across all ages (Supporting Information Figure S2I–K), consistent with previous studies (Ritchie et al., 2018; Wang, Xu, Luo, Hu, & Zuo, 2019). There were no significant differences between sexes in professional qualifications and years of education.

Functional brain connectivity was assessed for each participant using the negative and positive correlations between 21 nodes that correspond to the resting-state functional MRI networks shown in Figure 1A and Supporting Information Figure S3 (ec12Methods, section ec16Functional Image Preprocessing). First, we computed classical single-layer connectivity measures, namely, the average connectivity for the whole correlation network, followed by the average negative connectivity, the average positive connectivity, and the number of negative correlations (ec12Methods, section ec21Connectivity Measures). Then, the negative and positive correlations of each functional network were separated (Figure 1B.1 and 1B.2) and analyzed as two independent layers. To evaluate the topological organization of each single layer, we used two measures: the clustering coefficient and the global efficiency. The clustering coefficient is a measure of segregation that increases with the number of local connections and represents the average clustered connectivity around all nodes in the network. The global efficiency is a measure of integration that increases as the paths connecting any two nodes in the network become shorter and estimates the average capability with which different nodes communicate with each other (Mijalkov et al., 2017) (ec12Methods, section ec23Single-Layer Network Measures). To assess the relationship between these two layers, we integrated them into a multiplex network (Figure 1C) and a multilayer network (Figure 1D). In the multiplex network approach (Figure 1C), each node in the positive layer was connected with the same node in the negative layer. We computed two multiplex measures: the multiplex clustering coefficient, a measure that increases with the local connections in neighboring nodes between the two layers, and the multiplex participation, which is a measure of integration that assesses how evenly a node is connected in the two layers (Battiston, Nicosia, & Latora, 2014) (ec12Methods, section ec26Multiplex Network Measures). A disadvantage of the multiplex approach is that the relation between the two layers is local being only allowed between the same nodes. To address this limitation, in the multilayer network approach we connected each node in one layer to every node in the other layer (Figure 1D). The strength of the relationship between the two layers can be changed by adjusting the weight of the interlayer connections, σ. For each multilayer network, we define σ as a fraction of the strongest functional connections in the corresponding network, and evaluate the measures’ ability to characterize sex differences across the wide range of σ values. We developed two new measures to assess the integration and segregation properties of these multilayer networks. Specifically, we calculated the multilayer global efficiency, which compares the global efficiency differences due to the intra- and interlayer connections. Similarly, we also calculated the multilayer clustering coefficient, which compares the clustering coefficients or triangles in all nodes between the two layers (ec12Methods, section ec29Multilayer Network Measures).

In a first step, to identify which simple connectivity measures showed the greatest differences between sexes over the course of aging, we compared men and women at all ages using a permutation test and used separate linear models that included whole-brain average connectivity, average negative connectivity, average positive connectivity, and the number of negative correlations as the outcome and age, sex, age², age × sex, and age² × sex interactions as predictors (the models were fit on the average data by sex and within each age; ec12Methods, section ec35Statistical Analysis). These models showed that women had significantly higher average functional connectivity than men (R² = 0.769; AIC = 244.651; MSE = 2.184; Figure 2A and Supporting Information Figure S4A), whereas men had a significantly higher number of negative connections compared to women (R² = 0.778; AIC = 501.106; MSE = 108.420; Figure 2D and Supporting Information Figure S4D) across a broad age range (50 to 71 years). These differences in connectivity between sexes diminished with increasing age, to the point where there were almost no differences in mean connectivity strength or number of negative connections between men and women after 75 years (age × sex interaction: average connectivity p = 0.01; number of negative connections p < 0.001). On the other hand, there were no significant differences between sexes in the average positive and negative connectivity strength (R² = 0.403; AIC = 281.229; MSE = 3.947 and R² = 0.408; AIC = 281.229; MSE = 2.013, respectively; Figure 2B and C and Supporting Information Figure S4B and C). These results suggest that women have higher connectivity strengths than men due to a lower number of negative connections, but these differences dissipate with increasing age.

To identify which single-layer topological measures showed the greatest differences between sexes over the course of aging, we repeated the above analysis by including them as dependent variables in separate linear models with age, sex, age², age × sex, and age² × sex interactions as predictors. These models showed that men had higher global efficiency than women in the positive and negative layers (R² = 0.924; AIC = −1.051; MSE = 0.055 and R² = 0.649; AIC = 88.603; MSE = 0.209, respectively; Figure 2G–H and Supporting Information Figure S4G–H), indicating that their functional connectomes were characterized by shorter paths in the networks with negative and positive connections. Interestingly, these sex differences remained constant across different ages, suggesting that they were independent of age (age × sex or age² × sex interaction not significant in global efficiency measures; Supporting Information: Linear Models). In contrast, no significant differences in the clustering coefficients in the positive and negative networks were observed between women and men (R² = 0.459; AIC = 226.729; MSE = 1.728 and R² = 0.533; AIC = −96.294; MSE = 0.013, respectively; Figure 2E and F and Supporting Information Figure S4E and F).

Our linear regression models showed that the multiplex participation coefficient was significantly higher in men than in women across all ages (R² = 0.891; AIC = −172.848; MSE = 0.004; Figure 2J, Supporting Information Figure S4J and Supporting Information: Linear Models), indicating that the men’s functional connectomes were characterized by nodes with a similar number of connections in the negative and positive connectivity layers. However, these differences decreased with age, similarly to the average connectivity and number of negative connections (age × sex interaction: p = 0.036). The multiplex clustering coefficient did not reveal any significant differences between men and women (R² = 0.096; AIC = −106.989; MSE = 0.282; Figure 2I and Supporting Information Figure S4I). We note that caution needs to be applied when considering multiplex clustering coefficients as their interpretation can be cumbersome for multiplex correlation networks (Masuda, Sakaki, Ezaki, & Watanabe, 2018; Zalesky, Fornito, & Bullmore, 2012).

The differences between men and women were robust for the complete range of interlayer weights σ for the multilayer clustering and for smaller values of σ in the case of multilayer global efficiency. However, the strongest differences were observed for the σ = 0.7 and σ = 0.2 in the case of multilayer clustering and multilayer global efficiency, respectively (Supporting Information Figures S5 and S6). Whereas the multilayer clustering coefficient differences between men and women remained stable with aging (R² = 0.799; AIC = −130.497; MSE = 0.008; age × sex interaction not significant; Figure 2K and Supporting Information: Linear Models), the differences in multilayer global efficiency decreased with aging (R² = 0.795; AIC = −82.061; MSE = 0.016; age × sex interaction: p < 0.001; Figure 2L and Supporting Information: Linear Models). Up to 80% of the variance in both multilayer measures was associated with sex differences over time, outperforming all previous models for simple connectivity, single-layer measures, and multiplex topological measures. For a more consistent comparison of the different measures, we have repeated our analysis when the functional connectivity metrics were fitted to a full model with age, sex, age², age × sex and age² × sex) as predictors (Supporting Information: Linear Models). These results further confirm that multilayer measures and multiplex participation are the best performing measures and suggest that complex network measures that account for the relationship between positive and negative functional connections are more sensitive to sex differences across aging.

Finally, we have conducted several analyses to assess the reproducibility of our findings. Specifically, we assessed the reproducibility of the multilayer measures as a function of sample size (Marek et al., 2022) (Supporting Information: Figures, Section VII), whether they were replicable at individual network densities (Supporting Information: Figures, Section VIII), as well as their test-retest reliability using the intraclass correlation coefficient (ICC) (McGraw & Wong, 1996) (Supporting Information: Test-Retest). Our findings showed that multilayer measures can uncover similar patterns of between-sex differences for different analysis parameters. Specifically, at most ages, they can detect the differences observed in the original sample in subsamples with sizes as low as 55%–60% of the original sample, as well as at individual network densities higher than 15%. Finally, their test-retest reliability was greatest in younger individuals and higher than that of test-retest reliability of single-layer and connection measures; their ICC values ranged from 0.4 to 0.65, suggesting a fair to good degree of clinical significance (Cicchetti & Sparrow, 1981).

Next, we assessed whether the observed sex differences in functional connectivity throughout aging were associated with the cognitive functions, structural brain measures and vascular risk factors that differed between men and women in our cohort (executive functions, visuospatial functions, blood pressure and heart attack prevalence, subcortical volumes, white matter hyperintensities; see ec2Results, section ec3Sample, and Supporting Information Figure S2). Due to the high collinearity between the functional connectivity measures (Supporting Information Figure S7), we examined these associations using partial least squares (PLS) regressions. The executive cognitive scores were best explained by the single-layer positive global efficiency, followed by the single-layer negative clustering, average positive connectivity and negative connectivity. The visuospatial cognitive scores were best explained by the multilayer functional connectivity measures, followed by the average connectivity, number of negative connections, and the positive single-layer global efficiency. Heart attack and high blood pressure prevalence were best predicted by multiplex participation and single-layer global efficiency of the positive connections, with blood pressure additionally being predicted by multilayer clustering coefficient, number of negative connections, average connectivity, and single-layer global efficiency of the positive network. Multilayer global efficiency was a significant predictor of subcortical brain volumes and white matter hyperintensities, in addition to the average connectivity, number of negative connections, multilayer clustering coefficient, multiplex participation, and single-layer clustering coefficient in the case of subcortical volumes and single-layer clustering coefficient in the case of white matter hyperintensities. The detailed VIP scores for each predictor and the PLS model cumulative explained variance for the predicted variables are shown in Figure 3 and Supporting Information: Linear Models. Altogether, these results indicate that, compared to the other measures, the multilayer and multiplex measures were the best predictors of cognitive functions, structural brain measures, and vascular risk factors that differed between men and women during aging.

To identify genetic variation associated with the 12 multilayer, single-layer and connectivity functional measures we examined 9,356,431 genetic variants with minor allele frequency > 1% in a total of 33,773 European participants (ec12Methods, section ec33Genetic Association Analyses). A total of 4 loci exceeded a genome-wide significance level among one or more of these traits (Figure 4 and Supporting Information Figure S8). A locus near Paired Box 8 (PAX8), a gene involved in sleep efficiency, diastolic blood pressure, development, and vulnerability to neurodegenerative diseases (Elliott et al., 2018; Foo et al., 2021b; Jones et al., 2019), showed an association with the majority of functional connectivity measures, including multilayer and multiplex clustering, as well as all single-layer measures. Colocalization analysis with eQTL summary statistics of 49 tissues in GTEx project and brain tissues in BrainSeq, ROSMAP, Braineac2, and CommonMind datasets (ec12Methods, section ec34Colocalization) (GTEx Consortium, 2020; Guelfi et al., 2020; Jaffe et al., 2018; Kerimov et al., 2021; Ng et al., 2017; Sieberts et al., 2020), suggests a consistent strong colocalization between this locus and gene expression patterns of Immunoglobulin Kappa Variable 1/OR2-108 (IGKV1OR2-108), COBW domain-containing protein 2 (CBWD2), and Forkhead Box D4 Like 1 (FOXD4L1) over the multilayer and multiplex clustering coefficients, average positive and negative connectivity, and all single-layer functional measures (Supporting Information: Genetic Results). Furthermore, a locus near Inositol Polyphosphate-5-Phosphatase A (INPP5A) gene showed significant associations with the single-layer positive global efficiency. Interestingly, the genome-wide association study (GWAS) lead variant at this locus colocalized with eQTL for INPP5A in brain dorsolateral prefrontal cortex (DLPFC) from CommonMind dataset. Finally, loci near Adenosine Deaminase RNA specific B2 (ADARB2) and D-Amino Acid Oxidase Activator (DAOA) genes, which have previously been linked with processes involved in normal memory functioning and ATP metabolism (Prata et al., 2012), were associated with multiplex clustering and multiplex participation coefficients respectively.

Complex network measures are becoming increasingly popular in the field of network science (van den Heuvel & Pol, 2010). In combination with large samples, these measures can improve our understanding of brain connectivity (Fornito et al., 2013), and, in particular, our ability to predict sex differences over the course of aging (Stumme, Jockwitz, Hoffstaedter, Amunts, & Caspers, 2020). Previous studies of brain connectivity and network topology have mainly focused on measures that exclude the negative connections or average the effects of positive and negative connections (Chan et al., 2014; Ferreira et al., 2016; Foo et al., 2021a; Tomasi & Volkow, 2012b; Zhang et al., 2016). However, changes in the balance between positive and negative connectivity could be a more sensitive marker of abnormalities that occur in men and women in middle and late adulthood. Here we show that measures that assess this balance are indeed better predictors of sex differences during aging. Furthermore, we also show that these measures are associated with genes implicated in aging-related physiological processes as well as cognition, brain structure, and vascular disease, which have been previously shown to differ between men and women (McCarrey et al., 2016; Ramirez & Sullivan, 2018; Ritchie et al., 2018; Sachdev et al., 2009; Weiss et al., 2003). Altogether these findings highlight the importance of integrating the information from positive and negative connections in a multiplex or multilayer network approach to provide a more holistic view of functional brain connectivity changes.

Increasing age is associated with lower physical fitness and worse cognitive abilities, being one of the greatest risk factors for the development of neurodegenerative diseases (Deary et al., 2009; Wyss-Coray, 2016). However, there are several factors that determine how an individual ages and his/her predisposition to develop certain diseases. Sex is one of these factors and it has recently received a lot of attention due to increasing recognition that precision medicine approaches should take into account biological sex to treat neurodegenerative diseases (Ferretti et al., 2018). Sex was found to affect the functional connectivity (Allen et al., 2011; B. B. Biswal et al., 2010; Ritchie et al., 2018); however, the location and nature of the functional connectivity differences between men and women vary across studies (Scheinost et al., 2015). Between-sex differences were most consistently observed in the default mode network (Allen et al., 2011), but also reported in other resting-state networks (Allen et al., 2011; B. B. Biswal et al., 2010; Filippi et al., 2013; Scheinost et al., 2015), as well as in the functional connectivity patterns between them (Goldstone et al., 2016; Satterthwaite et al., 2015; Stumme et al., 2020). These differences between men and women change with aging (Cosgrove, Mazure, & Staley, 2007), with studies showing different age-related trajectories in functional connectivity that can be either diverging or converging in different resting-state networks (Goldstone et al., 2016; Scheinost et al., 2015). However, other studies have reported that between-sex differences persist also in older individuals (Stumme et al., 2020). Such differing results suggest that the understanding of sex differences in the adult brain is still in its infancy and utilizing larger samples or more sensitive methods are needed to address these discrepancies.

Functional brain connectivity plays a crucial role in how brain networks communicate with each other, being closely associated with behavior and cognition (van den Heuvel & Pol, 2010). Studies have shown that aging is accompanied with an increase in the connections between these networks (Betzel et al., 2014; Chan et al., 2014; Damoiseaux, 2017; Ferreira et al., 2016). These increases in connectivity are thought to be due to a less efficient use of neural resources in older individuals, who tend to over-recruit brain networks to compensate for the detrimental effects of aging (Geerligs, Renken, Saliasi, Maurits, & Lorist, 2015; Goldstone et al., 2016; Park & Reuter-Lorenz, 2009). Here, we confirm these findings in a larger sample of middle-age and old adults by showing that aging is associated with increases in the average whole-brain connectivity and decreases in the number of negative connections. Moreover, our findings revealed that men displayed lower average connectivity and a greater number of negative connections than women at young ages, but these differences dissipated with increasing age, to the point that no differences between sexes were observed any longer at old ages. These findings agree with results from previous studies showing different functional trajectories between men and women (Satterthwaite et al., 2015; Scheinost et al., 2015). However, in contrast to other studies (Stumme et al., 2020), here we demonstrate that, although the functional connectomes of younger men and women are different, they become increasingly similar with older age, possibly due to a faster rate of functional changes observed in the brains of men (Foo et al., 2021a), which is in line with studies showing that men have lower resilience to age-related cognitive decline compared to women (McCarrey et al., 2016).

Regarding measures of network organization, we found that the positive global efficiency was higher in men compared to women, but these differences remained stable across different ages. Similar differences were found in the negative global efficiency; however, in this case, the differences were not stable and seemed to decrease with aging. The global efficiency is used to assess integration or the ability for an efficient processing to occur between distant brain networks (Fornito et al., 2013); an excessive integration is thought to impede the ability of the brain to process information in a meaningful way (Lord, Stevner, Deco, & Kringelbach, 2017). Furthermore, global efficiency is associated with the presence of long range connections between networks. The healthy brain connectome is characterized by a low number of such connections because long range connections are associated with higher metabolic costs and minimizing these costs is essential for evolution (Bullmore & Sporns, 2012). Therefore, these differences between sexes might provide clues on why men are more prone to develop specific diseases with age, for example, Parkinson’s disease and epileptic seizures, which are associated with abnormal organization in the functional connectivity networks as a result of increased integration (Lord et al., 2017; Mijalkov, Volpe, & Pereira, 2022).

To evaluate whether the relationship between the positive and negative connections can reveal additional insights into sex differences over the course of aging, we integrated these connections as two separate layers in a multiplex and a multilayer networks. This relationship estimates the brain’s ability to inhibit certain functional connections in order to modulate the coactivation between different resting-state networks. Because the positive and negative connections can arise from different neurovascular mechanisms (Goelman, Gordon, & Bonne, 2014), the presence of both connection types is necessary to reach a more balanced communication between resting-state networks and consequently a better cognitive performance (Saberi, Khosrowabadi, Khatibi, Misic, & Jafari, 2021); in fact, the lack of anticorrelation or negative connections has been associated with lower cognitive control and working memory performance during neurodevelopment (Chai, Ofen, Gabrieli, & Whitfield-Gabrieli, 2014). This balance is disrupted in older adults, which over-recruit the functional networks as a compensatory mechanism to maintain or improve function (Goldstone et al., 2016; Park & Reuter-Lorenz, 2009), for example, by engaging both the left and right hemispheres (HAROLD model) (Cabeza, Anderson, Locantore, & McIntosh, 2002) or showing more distributed activation when compared to younger adults in performing verbal working and long-term memory as well as category-learning tasks (Park & Reuter-Lorenz, 2009).

In the multiplex network, where only connections between the same nodes in the two layers are allowed, we measured the balance between the positive and negative connections by the multiplex participation coefficient, which quantifies the heterogeneity of the connectivity patterns of a given node in both layers. We found that men had higher multiplex participation coefficients compared to women, indicating that they had a higher balance between positive and negative connections in the two layers that decreased with aging at a faster rate compared to women. Consistent with earlier findings, our results suggest that the less efficient connectomes reported in older persons may be related to an increase in the number of positive connections at the expense of negative connections (Ferreira et al., 2016), thereby disrupting the balance between coactivation and inhibition connectivity between functional networks.

The multilayer network approach, where connections are allowed between all nodes in the two layers, extends the multiplex approach by providing an estimation of the degree to which the between-layer balance affects the topological organization of the network. We found that men had higher multilayer clustering and multilayer global efficiency compared to women, which both decreased with aging. These findings indicate that the relation between positive and negative layers is stronger in men than in women. This might come at the expense of greater neural resources and metabolic costs, which could predispose the male connectome to the effects of increased oxidative stress and poor antioxidant defense mechanisms, which have been suggested to accompany higher brain connectivity (Griffa, Baumann, Thiran, & Hagmann, 2013) and could potentially lead to steeper rates of cognitive decline (McCarrey et al., 2016).

When the different connectivity and topological network measures were compared to each other, we found that the multiplex and multilayer network measures were the variables that were best predicted by sex differences over aging. In particular, age and sex were able to predict multiplex participation coefficient, multilayer clustering and multilayer global efficiency by explaining up to 89.1%, 79.9%, and 79.5% of the variance. These findings indicate that the integration of positive and negative connections as separate layers in a complex network approach is sensitive to important age- and sex-related variability not captured by conventional measures. This approach could thus be used to understand why men and women age differently. For instance, we found that cardiovascular risk factors such as hypertension and heart attack prevalence, which were higher in men between 51 and 79 ages, were best predicted by these measures. In addition, differences between sexes that remained stable over aging such as lower visuospatial cognitive scores and lower subcortical volumes in women compared to men were also best explained by multilayer measures.

Regarding the results of the GWAS analysis, a locus near the PAX8 gene showed an association with the majority of functional connectivity measures. Interestingly, this locus has been associated with sleep efficiency, diastolic blood pressure, insomnia, and sleep length (Jones et al., 2019). These data are consistent with the aging process resulting in changes in sleep habits and high blood pressure (Ancoli-Israel, 2009; Ramirez & Sullivan, 2018), both of which have been linked to functional connectivity changes (Foo et al., 2021b; Neitzel et al., 2021; Tagliazucchi et al., 2012). The PAX8 protein has also been associated to the regulation of multiple genes involved in thyroid hormone synthesis, which is necessary for brain development and function, for example, through processes such as neuronal differentiation, synaptogenesis, and dendritic proliferation (Foo et al., 2021b). All functional measures associated with the PAX8 gene demonstrated substantial colocalization with the FoxD4L1 gene, which is similarly implicated in processes that promote the onset of neural differentiation (Klein et al., 2013). They also demonstrated colocalization with CBWD2, which has been linked to sleep duration (Foo et al., 2021b), and IGKV1OR2-108, which has been found to be elevated in the livers of type 2 diabetes patients (Li, Pan, & Yang, 2019) and can lead to abnormal functional connectivity (Chen et al., 2014).

The rs4309079 locus, associated with single-layer positive global efficiency, has been linked to functional connectivity measures in previous research (Elliott et al., 2018). It is located adjacent to the INPP5A gene, which is involved in calcium signaling. As a fundamental cellular mechanism, inositol calcium signaling is expected to play a role in a range of neurological pathways that underlie functional connectivity (Neitzel et al., 2021). The INPP5A gene has been further associated with brain age explained by changes in functional connectivity and decreased system segregation (Neitzel et al., 2021; S. M. Smith et al., 2020), which is consistent with our current findings.

Multiplex participation showed a strong association in a locus near the DAOA gene, which is known to play an important role in the control of glutamatergic transmission (Prata et al., 2012). Glutamate is the most common excitatory neurotransmitter in mammals, and glutamate’s activation of NMDA receptors is critical for normal memory function. Glutamatergic antagonists (e.g., ketamine) have been demonstrated to lower performance on tests of declarative memory, verbal fluency, and problem solving, all of which have been linked to aging (Prata et al., 2012). We nevertheless could not find any striking colocalization evidence in the examined gene expression QTL datasets. Finally, for multiplex clustering, we found a locus near the ADARB2 (index variant: rs2152237) gene, which is involved in the ATP/ITP metabolic pathway. As a result, our findings add to our understanding of the genetic influences on functional connectivity and provide a link between the functional and genetic architecture of the brain, which might be relevant in explaining the changes in a variety of biological processes throughout healthy aging.

This study has some limitations. First, a longitudinal study design would have been more appropriate to assess age-related differences between sexes in functional brain connectivity. In particular, recent studies have shown that results regarding aging-related changes obtained from cross-sectional and longitudinal designs can be different from each other (Vidal-Piñeiro et al., 2021; Xing, 2021); therefore, our findings should be interpreted with this limitation in mind. However, multiple efforts to collect large longitudinal samples are currently underway (Liu et al., 2021; Miller et al., 2016), which opens the opportunity for future studies to overcome this limitation. Second, there was a considerable overlap between the values for individual men and women in functional connectivity measures. This overlap has also been observed in previous studies (Joel & Fausto-Sterling, 2016; Ritchie et al., 2018; Stumme et al., 2020; Zhang et al., 2016; Zhang, Dougherty, Baum, White, & Michael, 2018), indicating that the prediction of an individual’s sex from a small number of functional connectivity measures is challenging (Zhang et al., 2018). Nonetheless, our findings revealed consistent group differences between men and women across a wide age range (47–79 ages), suggesting that, although there is considerable variability, some changes seem to be quite robust (Chekroud, Ward, Rosenberg, & Holmes, 2016). Finally, since in this study the nodes of the networks corresponded to 21 resting-state networks derived from independent component analyses (Miller et al., 2016), the negative connections are related to anticorrelations between resting-state networks, which have been consistently observed and demonstrated to have a neurophysiological basis (Chai, Castañón, Öngür, & Whitfield-Gabrieli, 2012; Fox et al., 2005; Fox, Zhang, Snyder, & Raichle, 2009). Thus, these connections should be interpreted differently than the negative correlations between brain regions analyzed in studies using brain areas as nodes, which are still not clearly interpretable (Chai et al., 2012; Fornito et al., 2013).

To summarize, in this study we developed novel multilayer connectivity measures in order to assess the connectivity patterns and topological architecture between resting-state networks in a large cohort of middle age and old adults. We showed that these multilayer measures are superior at capturing sex-related effects during aging when compared to simpler connectivity measures that do not account for the relationship between positive and negative connections. The multilayer measures were also significant predictors of sex differences in cognitive, structural, and cardiovascular measures, and they were associated with genes that have previously been implicated in aging-related processes. Thus, our findings highlight the importance of studying the balance between positive and negative functional connections to understand the effects of sex over aging, which should be included in future studies.

The UK Biobank cohort is a large population-based study with more than 500,000 participants from the United Kingdom (https://www.ukbiobank.ac.uk/). Following an initial visit for collection of medical and other clinical information, 37,704 individuals underwent MRI. To achieve robust group comparisons, we limited our analyses to age groups with at least 50 participants, resulting in a sample size of 37,543 people in the age range 47–79 years (17,568 men: mean age = 64.75; SD = 7.58 and 19,975 women: mean age = 63.49; SD = 7.33).

The functional MRI scans were performed on a standard Siemens Skyra 3T scanner using an echo-planar imaging (EPI) sequence with the following parameters: duration ∼ 6 min; 490 time points; repetition time = 735 ms; echo time = 39 ms; field of view = 88 × 88 × 64; voxel size = 2.4 mm³; flip angle = 52°. Participants were instructed to relax and think of nothing in particular while focusing their eyes on a crosshair during the scan. Regarding structural MRI, the T1-weighted images were obtained using a 3D magnetization-prepared rapid gradient-echo imaging sequence with the following parameters: 208 slices, echo time = 880 ms; repetition time = 2,000 ms; field of view = 208 × 256 × 256; voxel size = 1 mm³ (Alfaro-Almagro et al., 2018).

All images obtained by the different imaging modalities were analyzed using an image-preprocessing pipeline run by the UK Biobank imaging core, who also performed quality assessment of the images (Alfaro-Almagro et al., 2018). In the following, we summarize the procedure for the different modalities; full details of the preprocessing and quality control are available elsewhere (Alfaro-Almagro et al., 2018; Miller et al., 2016).

After data acquisition, a number of preprocessing steps were carried out using FSL, including motion correction using MCFLIRT, grand-mean intensity normalization by a single multiplicative factor, high-pass temporal filtering with a Gaussian-weighted least-squares straight line fitting (σ = 50.0 s), EPI unwarping by using a field map obtained before data collection, gradient distortion correction (GDC) unwarping, and removal of all artifacts by an ICA-based X-noiseifier. Finally, all datasets underwent temporal demeaning and variance normalization.

The preprocessed data of 4,100 participants was used for a Group-ICA analysis. Using FSL’s MELODIC tool and FSLNets toolbox (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki), a spatial-ICA with a dimensionality of 25 components was applied and the resulting ICA maps were mapped onto each subject’s resting-state fMRI time series data to generate one representative timeseries per ICA component. During this procedure, four networks were identified as artifacts and were discarded from further analysis. This resulted in a 21 × 21 connectivity matrix for each participant, where the functional connectivity between each pair of ICA spatial maps is characterized by full normalized temporal correlation. Analyzed networks included the default mode, fronto-parietal, sensorimotor, visual, attention, subcortical, and temporal networks (Figure 1A); the group-ICA spatial maps can be found at https://biobank.ctsu.ox.ac.uk/crystal/refer.cgi?id=9028, and the average functional connectomes for several representative age groups are shown in Supporting Information Figure S2. The complete description of the preprocessing procedures can be found elsewhere (Alfaro-Almagro et al., 2018; Miller et al., 2016) (online documentation: https://biobank.ctsu.ox.ac.uk/crystal/docs/brain_mri.pdf).

T1-weighted scans were preprocessed using the standard procedures of the FreeSurfer pipeline (version 6.0; https://surfer.nmr.mgh.harvard.edu/). We calculated the cortical thickness for each individual by averaging the regional cortical thicknesses from 68 regions from the Desikan–Killiany atlas (Desikan et al., 2006). In addition, we calculated the subcortical volumes for each individual by averaging the volumes of all subcortical gray matter structures (cerebellar cortex, thalamus, caudate, putamen, pallidum, hippocampus, amygdala, accumbens), which were corrected for total intracranial volume using a regression approach (O’Brien et al., 2011). This data was available for a subsample of 17,317 men and 19,842 women.

The white matter hyperintensities were segmented from a combined T1 and T2-FLAIR images, using a fully automated supervised method based on the k-nearest neighbors algorithm (Griffanti et al., 2016).

At their initial visit, individuals were questioned about their medical history, including whether or not they had a high blood pressure (subsample of 17,517 men and 19,911 women), a heart attack (subsample of 15,204 men and 16,918 women), a heart attack–related angina episode (subsample of 14,224 men and 15,520 women), or a stroke (subsample of 14,224 men and 15,520 women). We calculated the percentage of men and women diagnosed with the previous conditions at all age groups, and used these percentages as dependent variables in further analyses. Since the diagnosis of cardiovascular diseases was performed at a different time point than the brain scanning, we included the time difference between diagnosis and scanning as a covariate in the analyses with cardiovascular variables. Furthermore, as there were several age groups that did not have any individuals with heart attack and hypertension, the models for these variables were corrected for age prior to the PLS regression.

The cognitive assessments were administered on a touch screen and took place at the same visit as the brain scans. We averaged the standardized z-scores of 10 cognitive tests (Fawns-Ritchie & Deary, 2020) into four different cognitive domains: attention/psychomotor speed (reaction time, trail making - numeric path, symbol digit substitution tests; subsample of 11,531 men and 13,125 women), memory (numeric memory, paired associative learning, prospective memory, pairs matching tests; subsample of 11,657 men and 13,284 women), executive (fluid intelligence/reasoning, trail making - alphanumeric path tests; subsample of 11,535 men and 13,159 women), and visuospatial/visuoconstructional (matrix pattern completion, tower rearranging tests; subsample of 11,452 men and 13,013 women). More details about the cognitive tests are available at https://biobank.ndph.ox.ac.uk/showcase/label.cgi?id=100026.

The connectivity measures were calculated on the weighted connectivity networks consisting of positive and negative connections (shown in Figure 1A). For each weighted connectivity network, we calculated a corresponding binary network in which the individual connections retained their weight if they exceeded a certain threshold and were set to zero if they did not. In this process, the absolute value of each connection was compared to the threshold; however, their sign was preserved in the resulting binary matrix (i.e., negative connections in the weighted network remained negative in the binary network). As there are multiple thresholding approaches and there is currently no consensus as to which network density should be used (Fornito et al., 2013), we performed the thresholding at a density range of 6% to 33%, in steps of 1%. For densities below 6%, the networks became largely disconnected with fewer edges than nodes, whereas 33% was the maximum density that could be reached by all men’s and women’s networks of positive and negative connections (Supporting Information Figure S12).

After calculating all measures at each density within the complete density range, we evaluated the corresponding area under the curve (AUC) value, which was used to assess the between-sex differences (Supporting Information Figure S13). The AUC value was obtained by numerically integrating the measure values over the density range; this procedure resulted in a single numerical value for each network measure across the range of densities. As this analysis takes into account the entire density range, it is considered to be less susceptible to the thresholding process (Fornito et al., 2013).

We split the weighted connectivity networks into networks consisting of positive connections and negative connections (Figure 1B). As network measures are not defined for negative weights, the weights in the network of negative connections were substituted by their absolute values. Therefore, the topological analysis was equivalent for both positive and negative connections networks; in the following we omit the subscript (pos/neg) for clarity. Both networks were independently binarized at the density range 6%–33%, and the AUC value for all measures was obtained across this range, as outlined above. The topology of these networks was evaluated by calculating the global efficiency (SL-Ge) and clustering coefficient (SL-CC).

For each individual, we built a multiplex network with two layers, one with positive connections and the other with negative connections (Figure 1C). The multiplex networks were calculated at densities ranging from 6% to 33% by combining the two single-layer networks at the respective densities. At each density the multiplex networks can be represented by a supra-adjacency matrix, W, which consists of the intralayer adjacency matrices on the main diagonal and the interlayer connections in the off-diagonal entries. The interlayer connections are formed only among the node’s replicas, that is, $Wi∈αj∈β$ = 1 ⇔ i = j, where α and β denote layer 1 and layer 2, respectively. We evaluated the topology of these networks by calculating the participation and clustering coefficients and evaluating the AUC over the complete density range.

PLS performs a linear decomposition of the predictor and predicted variable matrix into latent variables (LVs) optimized so that the covariance between the resulting predictor and predicted matrix components (called factors and loadings, respectively) is maximal (Abdi & Williams, 2013). We fit a PLS regression model for each predicted variable (cognitive tests, structural measures, and prevalence of cardiovascular diseases) independently. The optimal component number for the decomposition was determined in each case by the Bayesian information criterion (BIC) based on the estimated degree of freedom of the PLS model (Krämer & Sugiyama, 2011). In particular, predictor and predicted variable matrices were decomposed into 2–12 LVs, with the optimal number of LVs based on BIC scores being defined as follows: attention = 5, executive = 6, memory = 2, thickness = 12, visuospatial = 6, volumes = 5, WMH = 5, heart attack = 3, hypertension = 4. The contribution of each variable to the prediction was quantified via variable importance in the projection (VIP) scores, calculated as the sum of PLS weights over LVs, weighted by the variance explained by each LV. Variables were defined as significant contributors to the prediction based on a VIP score of > 1 (Chong & Jun, 2005).

We examined the association between 9 million high-quality imputed common variants (MAF > 1%, imputation INFO score > 0.8 and Hardy–Weinberg equilibrium P > 10⁻¹⁰) from UK Biobank and each derived phenotype using linear mixed-effects models as implemented in REGENIE (Jamialahmadi et al., 2021; Mbatchou et al., 2021). All traits were rank-based inverse normal transformed before the analysis, and adjusted for age at MRI, sex, age², age * sex, age² * sex, body mass index (BMI), the first 10 principal components of ancestry, and genotyping array. For step 1, we used a subset of high-quality directly genotyped variants as described before (Jamialahmadi et al., 2021). We next performed a stringent physical LD clumping (PLINK parameters: –clump-p1 5e-8 –clump-r2 0.05 –clump-kb 1000, after excluding individuals with third-degree or closer relatives) (Bycroft et al., 2018; Chang et al., 2015), followed by an approximate step-wise model selection in conditional and joint multiple-SNP analysis (COJO-GCTA) (Yang et al., 2012), with a window of 1 Mb and using 50,000 randomly selected unrelated Europeans (in-sample LD structure) as described before (Jamialahmadi et al., 2021).

Colocalization was performed between independent genetic loci from COJO-GCTA and summary statistics of gene expression quantitative trait loci (eQTL) of 49 tissues in GTEx v8, BrainSeq, ROSMAP, Braineac2 and CommonMind from eQTL cataolgue release 4 (GTEx Consortium, 2020; Guelfi et al., 2020; Jaffe et al., 2018; Kerimov et al., 2021; Ng et al., 2017; Sieberts et al., 2020). All genes with at least one significant association (FDR adjusted p value < 0.1) within a window of 1 Mb around each index variant were tested using coloc R package (Giambartolomei et al., 2014) with default priors, and H4 posterior probability (PP) > 0.8 was considered as a strong evidence that both traits share the same causal variant.

The statistical significance of the differences between men and women was assessed by performing nonparametric permutation tests with 10,000 permutations, which were considered significant for a two-tailed test of the null hypothesis at p < 0.05. These results were adjusted for multiple comparisons by applying false discovery rate (FDR) corrections at q < 0.05 using the Benjamini–Hochberg procedure (Benjamini & Hochberg, 1995). All results remained statistically significant after adjusting for average cortical thickness and subcortical volumes, suggesting that the observed between-sex differences were independent of brain atrophy (data not shown). To detect sex differences over aging, we used linear regression models where all functional connectivity measures for men and women at each age were summarized by their average value due to the high overlap between brain measures across different ages. Then, they were included as the dependent variables in the linear regression models, whereas age, sex, age², age * sex, and age² * sex were included as independent variables. The best model for each measure (Supporting Information: Linear Models) consisted of a combination of predictors that resulted in the minimum value of the Akaike information criterion (AIC) for that model. The significance of the overall model and the independent coefficients was evaluated by an F test, which was considered significant at p < 0.05. To compare the performance of the different models, we used the AIC and the mean squared error of all models.

The authors did not participate in data collection. The primary data source used in the study was the UK Biobank, which requires an application for access (https://www.ukbiobank.ac.uk/). This study was conducted under application number 37142. All data from eQTL Catalogue project are freely available at https://www.ebi.ac.uk/eqtl/.

“UK Biobank has approval from the North West Multi-centre Research Ethics Committee (MREC) as a Research Tissue Bank (RTB) approval.” (from https://www.ukbiobank.ac.uk/learn-more-about-uk-biobank/about-us/ethics). Informed consent was obtained from all UK Biobank participants. The current study is covered by this approval as we did not use any additional data or re-contacted the participants.

We thank the following: Swedish Research Council; Swedish Alzheimer Foundation; Swedish Brain Foundation; Strategic Research Area Neuroscience (StratNeuro); Center for Medical Innovation (CIMED); Foundation for Geriatric Diseases at Karolinska Institutet; Gamla Tjänarinnor; Stohnes Foundation; Lars Hierta Memorial Foundation.

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00286.

Mite Mijalkov: Conceptualization; Formal analysis; Funding acquisition; Methodology; Visualization; Writing – original draft; Writing – review & editing. Dániel Veréb: Formal analysis; Methodology; Writing – original draft; Writing – review & editing. Oveis Jamialahmadi: Formal analysis; Writing – original draft; Writing – review & editing. Anna Canal-Garcia: Software; Visualization; Writing – review & editing. Emiliano Gómez-Ruiz: Software; Writing – review & editing. Didac Vidal-Piñeiro: Writing – review & editing. Stefano Romeo: Formal analysis; Writing – original draft; Writing – review & editing. Giovanni Volpe: Conceptualization; Methodology; Supervision; Writing – original draft; Writing – review & editing. Joana Pereira: Conceptualization; Funding acquisition; Methodology; Supervision; Writing – original draft; Writing – review & editing.

Joana Pereira, Vetenskapsrådet (https://dx.doi.org/10.13039/501100004359). Joana Pereira, Alzheimerfonden (https://dx.doi.org/10.13039/501100008599). Joana Pereira, Hjärnfonden (https://dx.doi.org/10.13039/501100003792). Joana Pereira, Strategic Research Area Neuroscience (StratNeuro). Joana Pereira, Center for Medical Innovation (CIMED). Mite Mijalkov, Foundation for Geriatric Diseases at Karolinska Institutet. Mite Mijalkov, Stiftelsen för Gamla Tjänarinnor (https://dx.doi.org/10.13039/100010815). Mite Mijalkov, Gun och Bertil Stohnes Stiftelse (https://dx.doi.org/10.13039/100009673). Mite Mijalkov, Stiftelsen Lars Hiertas Minne (https://dx.doi.org/10.13039/501100004722).

Negative connections:

Anticorrelations (negative associations) between the nodes in the whole-brain connectome.

Positive connections:

Correlations (positive associations) between the nodes in the whole- brain connectome.

Connectivity measures:

Measures summarizing the average functional connectivity strength of the whole-brain connectome.

Multiplex network:

A two-layer multilayer network where interlayer connections are allowed only between corresponding nodes in the positive and negative layers.

Multilayer network:

A two-layer network where interlayer connections are allowed between all pairs of nodes in the two layers.

σ:

The fraction of the strongest functional connection in the corresponding multilayer network that is used to calculate the weight of the interlayer connections.

Multilayer global efficiency:

Measure that assesses the level of functional integration in a multilayer network.

Multilayer clustering coefficient:

Measure that assesses the level of functional segregation in a multilayer network.

Single-layer measures:

Classical network measures of clustering coefficient and global efficiency calculated for the networks of positive and negative connections separately.

Partial least squares (PLS) regression:

A multivariate regression analysis that projects sets of dependent and explanatory variables to a lower dimensional subspace consisting of latent variables that are maximally correlated.