Brain signaling becomes less integrated and more segregated with age

Aging is the number one risk factor for almost all neurodegenerative diseases (Kennedy et al., 2014). For every 5 years after the age of 65, the probability of acquiring Alzheimer’s disease doubles (Bermejo-Pareja et al., 2008). An influential conceptual approach to begin making sense of brain dynamics frames it in terms of a balance between integrated and segregated network states (Deco, Tononi, Boly, & Kringelbach, 2015; Friston, 2009; Sporns, 2010, 2013; Tononi, Sporns, & Edelman, 1994; Wig, 2017). On one hand, the brain faces functional pressure to have as many regions directly connected for quick communication. On the other hand, the brain is constrained to minimize metabolic energy consumption because it consumes 10 times more of the body’s energy than expected by mass (Raichle, 2006). Tuning the balance between extensive global signaling, referred to as integration, and limited local signaling, referred to as segregation, optimally compromises between functional and energetic constraints (Bullmore & Sporns, 2012; Cohen & D’Esposito, 2016; Manza et al., 2020; Wang et al., 2021). Although these constraints remain throughout life, age related glucose hypometabolism disrupts their balance.

Previous research has found mixed aging results, depending on the metrics used to measure integration and segregation (Chan, Park, Savalia, Petersen, & Wig, 2014; Chen et al., 2021; Onoda & Yamaguchi, 2013; Zhang et al., 2021). Although most in the literature use the system segregation metric (Chan et al., 2014), no consensus exists surrounding integration. In general, the problem facing the integration-segregation framework is that there is no one way to define the two states. Many graph theoretical metrics could potentially be used (Rubinov & Sporns, 2010), and it is unclear why one should take precedence over the other, particularly when their aging outcomes are mutually inconsistent. Thus, there is a need to more fundamentally define integration and segregation to transform it from a proxy to a physical quantity.

Here, we provide a physical foundation for the framework by applying the mean-field Ising model to treat integration and segregation as physical two-phase systems like magnets and liquids. After demonstrating that the Ising model can capture global brain dynamics as measured by functional MRI (fMRI) once the effective number of nodes is properly set, we proceed to calculate probabilities of being in the integrated or segregated states and find that younger and older brains are bounded by optimal and random signaling, respectively. We then explore diffusion and structural MRI data to ask if the age-related changes in signaling are due to changes in topological network connectivity.

We model human brain signaling patterns obtained from resting-state fMRI datasets. As in previous work (Weistuch et al., 2021), we capture those patterns with the Ising model, a widely used theoretical method for expressing macroscale behaviors in terms of interactions among many underlying microscale agents (Dill & Bromberg, 2010). We first transform the continuous fMRI data into a representation as discrete Ising spins via binarization of the data (Figure 1). That is, we reduce the state of the region as either −1 or 1 based on whether fMRI signaling is decreasing or increasing, respectively. Second, we calculate the synchrony by summing over all spins in a given time interval and dividing by the total number of spins (Figure 1). Synchronies are collected over the entirety of the scan to obtain a distribution. Based on the Ising model theory, the synchrony threshold delineating between integrated and segregated states is set such that P_int = P_seg = 1/2 at the Ising model’s critical point (see the Methods section). P_seg is the probability that the brain is in the segregated state and is defined as the relative number of time points for which the absolute value of synchrony is less than the synchrony threshold (Figure 1). P_int is defined as the relative number of time points for which the absolute value of synchrony is greater than the synchrony threshold and trivially relates to P_seg because P_int + P_seg = 1.

Before proceeding to calculate P_seg, we first check whether the model can capture the experimental synchrony distributions. A mean-field Ising model only considering pairwise interactions has one quantity of interest, λ. The strength of connection λ between any two regions corresponds to the degree to which signals between any two brain regions are correlated. However, we find that a naive fit of λ based on the principle of maximum entropy (Dill & Bromberg, 2010; Schneidman, Berry, Segev, & Bialek, 2006; Weistuch et al., 2021) fails to capture the synchrony distribution from the fMRI data (Figure 2, orange). To improve upon a standard Ising model approach, here, we introduce a hyperparameter N_eff. Brain atlas parcellations provide N brain regions; however, those N regions must be identically distributed across time for the Ising model to apply. We find that when setting N to a lower value N_eff, fixed for all individuals within a dataset, the Ising model accurately captures synchrony distributions (Figure 2). The optimal value of N_eff = 40 is determined by scanning across N_eff multiples of 5 to find which best captures the next order moment not fit by our maximum entropy setup across all individuals (see the Methods section; Figure 6). For our particular preprocessing (see the Methods section), we find that N_eff = 40 performs best for individuals in the Cambridge Centre for Ageing and Neuroscience (Cam-CAN) (Taylor et al., 2017) and the Human Connectome Project (HCP) Aging (Harms et al., 2018). For the UK Biobank (UKB) (Alfaro-Almagro et al., 2018), N_eff = 30 performs best (Figure 6).

Based on identified N_eff hyperparameter values, brains act as if they have a few tens of functional units. If different preprocessing decisions are considered, such as atlas resolution, N_eff values are still within an order of magnitude. At the voxel level (N = 125,879), we obtain an N_eff value of 65 for Cam-CAN and 125 for HCP using the same procedure as for the Seitzman atlas (N = 300) considered in the previous paragraph (Supporting Information Figure S2). Future work will pinpoint how N_eff depends on preprocessing to enable a future study creating a physics-based parcellation of the brain.

We also tried an alternative fitting strategy by fitting N_eff per individual rather than having the same value for all individuals in a respective dataset. We show that individually fitted N_eff values trivially relate to λ as expected by theory (Supporting Information Figure S1). Moreover, individually fitted N_eff are not found to be related to global differences in anatomical brain connectivity (Supporting Information Figure S3).

With an appropriately determined N_eff, we can accurately set the same synchrony threshold s* for all individuals within a dataset to calculate P_seg. The value of s* is set such that at the Ising model’s critical point in connection strength λ, P_seg equals to 1/2 for the ideal synchrony distribution based on the Ising model theory (see the Methods section). This enables P_seg comparisons across datasets that may have different N_eff values. For Cam-CAN and HCP, the value of s* is s* = 0.33 because N_eff = 40 for both datasets. For UKB, s* = 0.36 (Supporting Information Table S1).

Across the three publicly available datasets, we find that the balance shifts toward more segregation at older ages (Figure 3). Note that if we plotted P_int rather than P_seg, Figure 3 would be horizontally flipped, where P_int goes from high to low values as a function of increasing age because P_seg + P_int = 1. There is a large variation among the subjects (Supporting Information Figure S4). However, the correlation between age and P_seg is significant with the largest coefficient being 0.40 for Cam-CAN, while the lowest being 0.08 for UKB. Discrepancies in study designs may explain correlation magnitude differences: Cam-CAN and HCP are designed to study healthy aging (Bookheimer et al., 2019; Shafto et al., 2014), while the goal of UKB is to identify early biomarkers for brain diseases (Sudlow et al., 2015).

To better highlight how P_seg changes across Cam-CAN’s large age range, we present violin plots for younger, middle age, and older individuals’ P_seg (Supporting Information Figure S5). We also investigate how P_seg varies across time for a given individual. In Supporting Information Figure S6, we show that per individual P_seg, standard deviations decrease across age for Cam-CAN and HCP individuals. Finally, we perform a multiple linear regression with sex and handedness as additional covariates and show that age still strongly explains increasing segregation (Supporting Information Tables S2– S4 and Figures S7–S9).

Informed by the Ising model, increases in segregation result from network reorganization to more local signaling because of weakened connection strength between regions. Interestingly, younger individuals exhibit segregation behavior closer to the Ising model’s critical point of connection strength (Supporting Information Figure S10). At the critical point, we define P_seg = 1/2 (see the Methods section) and find experimental P_seg values closer to 1/2 for younger individuals (Figure 3). Older individuals, on the other hand, approach P_seg = 1 on average. This limit corresponds to functionally uncoupled brain regions that are randomly activating. Our results support the critical brain hypothesis that healthy brains operate near a critical point (Beggs, 2022; Beggs & Plenz, 2003; Ponce-Alvarez, Kringelbach, & Deco, 2023; Tagliazucchi, Balenzuela, Fraiman, & Chialvo, 2012) and implicate aging as pushing brain dynamics further away from criticality.

In the previous subsection, we discussed the disruption of the integration and segregation balance from the perspective of phase transitions in physics. Here, we explore the physiological mechanism underlying increasing segregation in the aging brain. We consecutively simulate the Ising model on a hypothetically degrading brain structure and show that random removal of edges yields qualitatively similar results to those of fMRI (Figure 4). Note that Figure 4 is horizontally flipped from that of P_seg (Figure 3) because the average degree (relative number of edges) is on the x-axis. It is presumed that edges are lost as age increases. In Figure 4, edges are lost linearly in time; however, more complicated monotonic functions can be employed to yield a quantitative match with experimental data in Figure 3. We can also capture variability among individuals by assuming that connection strengths within an individual are drawn from a distribution, rather than all being equal (Supporting Information Figure S11). In the Supporting Information, we also demonstrate that similar qualitative trends are obtained when starting with other individuals’ structures, regardless of their age (Supporting Information Figure S12).

We now begin to investigate possible mechanisms of connection degradation. First, we find that our simulation is agnostic to the detailed mechanism of connection degeneration because connection strength is essentially modulated by the probability that a given edge exists (Supporting Information Figure S13). In other words, the simulation cannot inform whether connections are degraded based on some targeted property. Thus, we turn to structural MRI and diffusion MRI data from UKB to investigate possible properties being degraded with age. In Supporting Information Figure S15, we confirm that white matter volume decreases as a function of adult age, as previously reported (Bethlehem et al., 2022; Lawrence et al., 2021; Lebel et al., 2012). However, this decrease does not correspond to a loss of anatomical connections because we find that neither average degree, average tract length, nor average tract density monotonically decrease with age when analyzing diffusion MRI scans using the Q-Ball method (Supporting Information Figure S16). This seems to contradict previous findings, which report decreases (Betzel et al., 2014; Lim, Han, Uhlhaas, & Kaiser, 2015). However, previous results employed the more simple diffusion tension imaging (DTI) method, which is known to be less accurate at performing tractography (Garyfallidis et al., 2014; Jones, Knösche, & Turner, 2013; Rokem et al., 2015). When rerunning our analysis for DTI, we can reproduce previously reported tract properties’ anticorrelations with age (Supporting Information Figure S16). We also investigate a graph property that captures polysynaptic connectivity called communicability (Andreotti et al., 2014; Estrada & Hatano, 2008; Seguin, Sporns, & Zalesky, 2023) and find that it also does not decrease age when using Q-Ball-derived tract density (Supporting Information Figure S17).

We propose that observed white matter volume reduction (Supporting Information Figure S15) and brain dynamics change correspond to less myelin covering axons as functions of age. Despite rejecting anatomical connections as a possible mechanism in the previous paragraph, it remains inconclusive whether myelin underlies trends because we are not aware of such data being publicly available. Although axons are still physically present, myelin coverage disruption causes regions to no longer be functionally connected because signals do not arrive on time. Previously reported results from myelin water imaging confirm reduction in myelin at advanced ages (Arshad, Stanley, & Raz, 2016; Buyanova & Arsalidou, 2021). We also investigated whether degraded functional connections are likely to be longer than average with age, as previously reported for certain brain regions (Tomasi & Volkow, 2012). Although we indeed find that the average correlation of the 25% longest connections is slightly more strongly anticorrelated with age compared with the average correlation of the 25% shortest connections for Cam-CAN (Supporting Information Figure S18, left), we find the opposite trend for HCP (Supporting Information Figure S18, right). Thus, myelin reduction does not seem to have a stronger impact on longer connections and conclude that the loss of functional connections happens randomly with respect to length at the brain-wide scale.

We apply the mean-field Ising model to physically quantify integration and segregation at the emergent scale of the whole brain. From resting-state fMRI scans across three publicly available datasets, we find that brain dynamics steadily becomes more segregated with age. Physically, aging leads to brain dynamics moving further away from its optimal balance at the critical point. Physiologically, analyses of white matter properties point to random functional connection losses due to myelin degeneration as the possible culprit for more segregated dynamics. This expands upon our previous work finding metabolic dysfunction to underly brain aging (Weistuch et al., 2021), hinting that myelin may be especially vulnerable to energy imbalance.

The Ising model and integration-segregation frameworks are considered as the simplest approaches to capture dynamics in their respective fields. Thus, it is fitting to map segregated and integrated states in neuroscience to disordered and ordered Ising model phases in physics, respectively. One general challenge in applying graph theory to MRI-level data is identifying what constitutes a node (DeFelipe, 2010; Lacy & Robinson, 2020; Seung, 2012; Sporns, 2010; Wig, Schlaggar, & Petersen, 2011; Yeo & Eickhoff, 2016). We identify the best number of effective brain regions N_eff such that the Ising model accurately captures individuals’ synchrony distributions across the corresponding dataset, improving upon our original application of the Ising model, which lacked the N_eff hyperparameter (Weistuch et al., 2021). Future work will utilize N_eff calculations to guide the creation of a parcellation in which brain regions are constrained to be physically independent based on their collective functional activity.

The field is inundated with integration and segregation metrics that have different aging trends. We go beyond heuristic definitions, such as one that we previously proposed based on matrix decomposition (Weistuch et al., 2021), by self-consistently defining the two states within the Ising model framework. This makes our metric mechanistically based on the connection strength between regions and further stands out because P_seg and P_int are naturally at the emergent scale of the brain. Notably, we do not calculate a local property and then average over nodes to yield a brain-wide value (Wang et al., 2021’s metric also has this advantage). In addition, P_seg and P_int are directly related because P_seg + P_int = 1. Most integration and segregation metrics (Chan et al., 2014; Rubinov & Sporns, 2010; Tononi et al., 1994; Wang et al., 2021) are not defined to be anticorrelated. This could be advantageous because greater complexity can be captured (Sporns, 2010).

Taken together, it is not surprising that P_seg and P_int results are not consistent with some previous aging reports. For example, a property called system segregation, defined as the difference between inter- and intracorrelations among modules, was found to decrease with age (Chan et al., 2014). Although most report that segregation decreases with age, regardless of the specific metric (Chan et al., 2014; Damoiseaux, 2017; King et al., 2018; Zhang et al., 2021; see Chen et al., 2021 for an exception), integration trends are less clear. Global efficiency, taken from the graph theory, was found to increase with age (Chan et al., 2014; Yao et al., 2019); however, others found different integration metrics decreasing with age (Chong et al., 2019; Oschmann, Gawryluk, & The Alzheimer’s Disease Neuroimaging Initiative, 2020; Zhang et al., 2021), consistent with the results reported here.

The utility of the integration-segregation framework lies in its simplicity. However, its simplicity has led to various heuristic definitions that have qualitatively different aging trends. By physically defining integration and segregation based on connection strength between regions, we provide an interpretable foundation for more detailed studies going beyond the two-state approximation to investigate brain dynamics.

We access three publicly available resting-state fMRI datasets: Cam-CAN (Taylor et al., 2017), UKB (Alfaro-Almagro et al., 2018), and HCP (Harms et al., 2018). Acquisition details such as the field strength and repetition time can be found in Supporting Information Table S5. Demographic details can be found in Supporting Information Table S6.

UKB and HCP fMRI data are accessed in preprocessed form (for details, see Alfaro-Almagro et al., 2018, and Glasser et al., 2013, 2018, respectively). We preprocessed Cam-CAN data as done in our previous work (Weistuch et al., 2021). For all three datasets, the cleaned, voxel space time series are band-pass filtered to only include neuronal frequencies (0.01–0.1 Hz) and smoothed at a full width at half maximum of 5 mm. Finally, we parcellate into 300 regions of interest according to the Seitzman atlas (Seitzman et al., 2020). For our voxel-wide analysis presented in the Supporting Information, we do not perform parcellation and just consider gray mater voxels by masking.

Applying the Ising model requires the data to only take two possible values: −1 or 1. After performing the preprocessing outlined in the previous paragraph, we binarize the continuous signal for a given region based on the sign of the slope of subsequent time points (Weistuch et al., 2021). We previously showed that such binarization still yields similar functional connectivities as that of the continuous data (Weistuch et al., 2021).

Finally, we only consider brain scans that have the same number of measurements as the predominant number of individuals in the respective dataset (Supporting Information Table S5). If the fitted connection strength parameter λ is less than 0, reflecting a nonphysical value, we do not include that individual’s brain scan in our analysis. In the HCP dataset, we excluded individuals aged 90 years or older because their exact age, considered protected health information, is not available.

In Figure 2, our maximum entropy fit (orange line) fails to qualitatively capture the synchrony distribution for an arbitrary individual. To rescue the fit, we replace N with N_eff (Supporting Information Equation S1). In the right plot of Figure 5, we demonstrate that a mean-field Ising model with N_eff = 40 accurately captures the fourth moment of synchrony 〈s⁴〉 across all individuals in Cam-CAN preprocessed under the Seitzman atlas. Note that N_eff is not a parameter like λ; rather, it is a hyperparameter because it takes the same value across all individuals within the dataset. N_eff is necessary because the Ising model systematically underestimates 〈s⁴〉 when Λ > 0 (left plot of Figure 5). Note that Λ corresponds to rescaling λ such that Λ = 0 is at the critical point (Supporting Information Equation S13).

To identify N_eff = 40 as the best value, we perform a parameter scan over multiples of 5 and identify the N_eff at which the root mean square error (RMSE) between 〈s⁴〉_exp and 〈s⁴〉_model is minimized (Figure 6). We choose the fourth moment because it is the next order moment that our maximum entropy fit does not constrain. It is not the third moment because the distribution is assumed to be an even function as indicated by our prior (Supporting Information Equation S1).

In the second line, the mean-field Ising model’s P(n) is inserted (Supporting Information Equation S2). Z corresponds to the partition function and ensures that P(n) is normalized. The constant s* is the synchrony threshold for which segregated and integrated microstates are delineated. We set s* such that P_seg = 1/2 when Λ = 0 according to theory. More specifically, we numerically calculate P_seg(Λ = 0) for a given N_eff and extrapolate to find s* (Supporting Information Figure S19). Proper calibration ensures that the theory is accurate and enables apples to apples P_seg comparisons across datasets with different N_eff. The list of s* values for the three publicly available datasets studied can be found in Supporting Information Table S1.

We simulate the Ising model on an initial structure informed by diffusion MRI under the Harvard-Oxford atlas (Makris et al., 2006) (64 regions) for an arbitrarily chosen UKB individual (subject ID: 6025360). If no edge exists between two regions, then the regions are uncoupled. If an edge does exist, then regions i and j are coupled and contribute λ * σ_i * σ_j to the system’s energy, where λ corresponds to the connection strength and σ corresponds to the spin state of the corresponding region (−1 or 1). Under the standard notation of the Ising model, λ = J/T, where J corresponds to the coupling constant and T is the temperature of the bath. The starting λ is set to 34.4, which is above λ’s critical point (starting P_seg ≈ 0.2). By definition, N_eff = N = 64 in the simulations. Based on the atlas resolution, simulating the Harvard-Oxford atlas provides an N_eff similar to those found for the experimental data (N_eff = 40 for Cam-CAN and HCP; N_eff = 30 for UKB).

The simulation for a given structure starts by randomly assigning the 64 nodes up or down spins. Then, for each time step, we attempt 10 spin flips 64 times, for a total of 2,500 time steps. Spin flips are accepted according to the Metropolis-Hastings algorithm (Metropolis, Rosenbluth, Rosenbluth, Teller, & Teller, 1953). The exact number of spin flip attempts or total time points does not matter, as long as equilibrium is reached. For example, we find that for λ values larger than those presented in the text, synchrony distributions become asymmetric and exhibit only one of the two peaks corresponding to the integrated state because of the high kinetic barrier of going from all down spins to all up spins.

Although the starting structure is informed by diffusion MRI (dMRI), resulting structures after computational edge removals are based on the posited removal strategy. Edges informed by dMRI are undirected and removal maintains undirectedness. Effectively two times as many edges are removed because both forward and backward edges are concurrently eliminated. In Supporting Information Figure S14, we demonstrate how synchrony distributions change as edges are computationally removed for a UKB individual (subject ID: 6025360), with a starting λ = 34.4.

We also investigate other individuals’ structures in the UKB to test the robustness of our qualitative results. We arbitrarily chose the following six individuals to widely sample different ages: subject IDs 6025360 (51 years old), 4712851 (57 years old), 3081886 (61 years old), 1471888 (65 years old), 4380337 (72 years old), and 1003054 (74 years old) (Supporting Information Figure S12). To ensure that the starting λ are comparable despite differing in the probability that two regions are connected (p_edge), we set λ₀ = 86.0 for all simulations such that λ = λ₀ * p_edge. For example, for subject ID 6025360, p_edge = 0.40, thus, the starting λ = 34.4.

Diffusion MRI processing to obtain structural information such as tract length and streamline count, which we call tract density, is outlined in our previous work (Razban, Pachter, Dill, & Mujica-Parodi, 2023). Briefly, we take preprocessed dMRI scans from the UKB (Sudlow et al., 2015) and calculate connectivity matrices using the Diffusion Imaging in the Python software (Garyfallidis et al., 2014). We input the Talairach atlas (Lancaster et al., 2000) to distinguish between white and gray matter. We perform deterministic tractography and reconstruct the orientation distribution function using Constant Solid Angle (Q-Ball) with a spherical harmonic order of 6 (Aganj et al., 2010). For Supporting Information Figure S16, we also do reconstruction using diffusion tensor imaging (Garyfallidis et al., 2014). To generate the starting structure for Ising model simulations, we input the Harvard-Oxford atlas for tractography because it parcellates the brain into fewer regions, making it more computationally tractable to carry out simulations, and is closer to N_eff values found for experimental data.

We thank Ying-Jen Yang, Anthony Chesebro, Charles Kocher, Jonathan Pachter, and Lakshman Verma for insightful discussions. The research described in this paper is funded by the White House Brain Research Through Advancing Innovative Technologies Initiative (NSFNCS-FR 1926781 to L.R.M.-P. and K.A.D.) and the Stony Brook University Laufer Center for Physical and Quantitative Biology (K.A.D.). Data collection and sharing for this project was provided by the Cam-CAN. Cam-CAN funding was provided by the UK Biotechnology and Biological Sciences Research Council (grant number BB/H008217/1), together with support from the UK Medical Research Council and University of Cambridge, UK. This research has been conducted using the UKB Resource under Application Number 37462. Research reported in this publication was supported by the National Institute On Aging of the National Institutes of Health under Award Number U01AG052564 and by funds provided by the McDonnell Center for Systems Neuroscience at Washington University in St. Louis. The HCP-Aging 2.0 Release data used in this report came from DOI: 10.15154/1520707. The HCP data repository grows and changes over time. The HCP data used in this report came from NIMH Data Archive DOI: 10.15154/1526427.

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00389. The following references are cited in the supporting information: Friedli and Velenik (2017), Kochmański et al. (2013), and Landau (1937).

Rostam M. Razban: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Visualization; Writing – original draft; Writing – review & editing. Botond B. Antal: Conceptualization; Data curation; Investigation; Software; Writing – original draft; Writing – review & editing. Ken A. Dill: Funding acquisition; Investigation; Methodology; Writing – original draft; Writing – review & editing. Lilianne R. Mujica-Parodi: Funding acquisition; Supervision; Visualization; Writing – original draft; Writing – review & editing.

Lilianne Mujica-Parodi, Division of Information and Intelligent Systems (https://dx.doi.org/10.13039/100000145), Award ID: 1926781.

Scripts necessary to reproduce figures and conclusions reached in the text can be found at github.com/lcneuro/2state_brain. Please refer to the respective publicly available dataset to access previously published data (Cam-CAN, UKB, and HCP) (Alfaro-Almagro et al., 2018; Harms et al., 2018; Taylor et al., 2017).

State:

A particular combination of physical properties. Here, we assume that brain networks can only occupy either the integrated or segregated state.

Integration:

A network state composed of global signaling.

Segregation:

A network state limited to local signaling.

Ising model:

A classic model in physics that was first applied to ferromagnetism. It includes pairwise interactions between binary spin states.

Phase:

Interchangeable with the word “state” for the purposes of this text.

Critical point:

The point where two phases coexist. In this text, it is where the synchrony distribution dramatically changes from bimodal (primarily integrated) to unimodal (primarily segregated).

Principle of maximum entropy:

A fitting strategy that satisfies user-defined constraints in the most agnostic way.

White matter:

Bundles of axons connecting brain regions.