Inferring excitation-inhibition dynamics using a maximum entropy model unifying brain structure and function

Abstract Neural activity coordinated across different scales from neuronal circuits to large-scale brain networks gives rise to complex cognitive functions. Bridging the gap between micro- and macroscale processes, we present a novel framework based on the maximum entropy model to infer a hybrid resting-state structural connectome, representing functional interactions constrained by structural connectivity. We demonstrate that the structurally informed network outperforms the unconstrained model in simulating brain dynamics, wherein by constraining the inference model with the network structure we may improve the estimation of pairwise BOLD signal interactions. Further, we simulate brain network dynamics using Monte Carlo simulations with the new hybrid connectome to probe connectome-level differences in excitation-inhibition balance between apolipoprotein E (APOE)-ε4 carriers and noncarriers. Our results reveal sex differences among APOE-ε4 carriers in functional dynamics at criticality; specifically, female carriers appear to exhibit a lower tolerance to network disruptions resulting from increased excitatory interactions. In sum, the new multimodal network explored here enables analysis of brain dynamics through the integration of structure and function, providing insight into the complex interactions underlying neural activity such as the balance of excitation and inhibition.


INTRODUCTION
The brain is a complex dynamical system whose functional properties are largely determined by the characteristics of its neurons and patterns of synaptic connectivity, resulting in a balance of excitatory (E) and inhibitory (I) interactions. For example, if the number of neurons that are coactivated from one signal is too high (increased excitation), the result is wide-scale activations and errant signal propagation across the brain's subnetworks. On the other hand, if the number of coactivated neurons is too low (increased inhibition), the propagation of the signal may diminish too quickly, limiting information transfer. The dynamical balance between excitation and inhibition is important for adjusting neural input/output relationships in cortical networks and regulating the dynamic range of their responses to stimuli (Kinouchi & Copelli, 2006) as well as the optimal dynamic range where information capacity and transfer are maximized (Shew et al., 2011). This is the central thesis of the criticality hypothesis, a phenomenon that suggests that neural networks and many aspects of brain activity self-organize into a unique configuration, sometimes called a critical state (Wilting & Priesemann, 2019). This state represents the transition of complex dynamical systems like the brain from order (balanced excitation-inhibition) to disorder (disrupted excitation-inhibition balance) and has found applications in many scientific domains, including neuroscience and clinical neurology (Cocchi et al., 2017;Hahn et al., 2017;Sornette, 2006;Tagliazucchi, 2017). Studies have demonstrated that the cortex operates near criticality during neuronal signaling (Beggs & Plenz, 2003;Hahn et al., 2017;Shew et al., 2009), as well as in studies utilizing blood oxygen level-dependent (BOLD) signals extracted from fMRI imaging (Haimovici et al., 2013;Lombardi et al., 2017;Rabuffo et al., 2021;Tagliazucchi et al., 2012). In fact, there is growing evidence from animal models and whole-cell recordings supporting the hypothesis that synaptic dysfunction leading to neuronal hyperexcitation may represent some of the earliest changes in the progression of neurodegenerative disease like Alzheimer's disease (AD; Busche & Konnerth, 2016;Palop et al., 2007;Petrache et al., 2019;Ren et al., 2018). However, the major challenge with early detection and intervention is that both normal aging and AD are associated with alterations to neural structure and function (McDonald et al., 2009;Schuff et al., 1999). This includes regional hypometabolism (Chételat et al., 2013;Curiati et al., 2011), white matter (WM) changes (Barrick et al., 2010;Michielse et al., 2010), Aβ deposition (Rodrigue et al., 2012;Rowe et al., 2010), and disrupted resting-state functional connectivity (Damoiseaux et al., 2008;Sheline et al., 2010;Wang et al., 2006). To improve our understanding of neurodegenerative diseases (accounting for major factors such as age, sex, or genetic phenotypes) and improve early detection, we investigate a model that can integrate microscale principles at a connectome level to bridge the gap between cell-to network-level degeneration. However, we acknowledge that some abstraction is required in this strategy; in models of large-scale effects, physiological information may be more abstract, and details of cellular processes potentially lost. While this may seem counterintuitive from a biological perspective, it is necessary for describing higher-level phenomena informed by MRI neuroimaging.
To this end, in this paper we introduce a method based on statistical physics to jointly model both brain structure and function via a pairwise maximum entropy model (pMEM). Our framework is inspired by the Ising model representation of brain dynamics whereby self-organized patterns of connectivity are formed through the spontaneous fluctuations of random spins (Reichl & Luscombe, 1999). This model has been used to characterize complex microscale dynamics of the human brain (Deco et al., 2008;Kadirvelu et al., 2017;Ostojic & Brunel, 2011;Tkačik et al., 2015), as well as macroscale interactions (Ezaki et al., 2017;Marinazzo et al., 2014;Nghiem et al., 2018;Niu et al., 2019;Nuzzi et al., 2020;Schneidman et al., 2006). Unconstrained maximum entropy models (MEM) have been shown to accurately represent spatiotemporal coactivations in neuronal spike trains (Roudi et al., 2009;Schneidman et al., 2006;Shlens et al., 2006) as well as patterns of BOLD activity (Ashourvan et al., 2017;Cocco et al., 2017;Ezaki et al., 2020;Watanabe et al., 2013). In fact, Zanoci et al. (2019) recently showed that the Ising model captures collective neuronal behavior during wakefulness, light sleep, and deep sleep when both excitatory (E) and inhibitory (I) neurons are modeled. Further, at the macroscale, Ashourvan et al. (2021) recently developed a maximum entropy-based framework that derives functional connectivity measures from Criticality: From physics, representing the state of a dynamical system between order and disorder.
Functional connectivity (FC): Undirected measure describing the statistical dependance between brain regions based on blood oxygen level-dependent (BOLD) signals from fMRI. Pairwise maximum entropy model (pMEM): A maximum entropy model that takes into account the average activity of a node as well as pairwise interactions.
Ising model: Network model in which node activity is represented with binary states and energy is computed based on pairwise interactions.

Maximum entropy model (MEM):
A statistical model of data representing the highest entropy that satisfies the constraint of prior knowledge.
intracranial EEG recordings; their findings suggest that structural connections in the brain give rise to large-scale patterns of functional connectivity by promoting coactivation between connected structures. Thus, MEM may be an ideal tool to model functional connectivity and ultimately link microscale interactions (such as excitation and inhibition in neuronal circuits) to the functional connectome (FC) captured through fMRI BOLD activity.
Described as a function-by-structure embedding (FSE), our model infers the organization of functional connectivity from global activity patterns (i.e., simultaneously considering the activity of more than two brain regions) constrained to the structural connectome. We present a robust numerical approach for our model, optimizing a constrained maximum likelihood estimation. The use of a structural connectome to inform the modeling of BOLD activity is motivated by a strong link between fMRI-based functional connectivity and white matter-based structural connectivity (Bettinardi et al., 2017;Honey et al., 2009;Shen et al., 2015). These studies suggest that models of functional dynamics should also be governed by the underlying structure to include direct and indirect connections between brain regions. Thus, if our model accurately describes large-scale brain activity patterns during rest, it will provide a much richer representation of functional interactions governing global dynamics that may give rise to hyperexcitation. With our framework we construct hybrid resting-state structural connectomes (rs-SC) for a group of 76 middle-aged and cognitively intact individuals. These unique structural networks are informed by a spin-glass-like Ising model, whose dynamics resemble that of traditional FC. We demonstrate that our new structurally informed networks can consistently and accurately reconstruct observed BOLD correlations. Investigating macroscale brain dynamics through the lens of statistical physics allows us to infer computationally the nature of resting-state activity (corresponding to inhibition or excitation) and probe potential disruptions to E/I balance that may lead to hyperexcitation and subsequent increased vulnerability to neurodegeneration. To evaluate this phenomenon, we create subgroups of 38 age-and sexmatched individuals based on whether one is a carrier of the apolipoprotein E (APOE)-ε4 allele, a well-known genetic risk factor of AD. Recent studies have shown that APOE-ε4 may contribute directly to early neuronal dysfunction, either directly via modification of the excitation/inhibition balance or linked with amyloid deposition (Bi et al., 2020;Koelewijn et al., 2019;Nuriel et al., 2017;Stargardt et al., 2015). Using our new hybrid rs-SC, we investigate the relationship between E/I balance and criticality in these two groups. We hypothesized that because of a shift in E/I balance towards hyperexcitation, the female APOE-ε4 carrier group would exhibit a lower tolerance to perturbations in the network when simulating brain dynamics using Monte Carlo simulations of the Ising model as compared with the female noncarrier group. Herein we aim to demonstrate that an increase in excitatory interactions at the connectome level, identified using our new hybrid connectome, may provide new evidence of vulnerability among females to AD neuropathology due to disruptions in E/I balance.

Constructing a Function-by-Structure Embedding Using a Constrained Maximum Likelihood Estimation
In constructing the function-by-structure embedding (FSE), we begin with the unconstrained pairwise maximum entropy model (pMEM) as described in the Methods section. The pMEM is sometimes referred to as the inverse Ising model, where the pairwise interactions (represented as J i,j , with i and j representing regions of interest, ROIs, in the brain network) are inferred from the observed data (BOLD time series). As the model assumes binary data, we binarized the resting-state fMRI signals obtained from the 76 cognitively intact middle-aged subjects. The binarized activity pattern of N = 80 ROIs at time t (t = 1, 2, …, t max ; t max = 236) is denoted s(t ) = s 1 (t ), s 2 (t ), …, s N (t ) 2 {−1, +1} N . Note that t max is determined as a result of the 8-min fMRI Structural connectome (SC): Network representation of physical connections in the brain. Nodes represent brain regions. Edges represent a measure of connectivity between them. APOE-ε4: An allele representing one of the strongest genetic risk factors for developing Alzheimer's disease (AD). scan time with TR = 2 s (see the Methods section). Here s 1 (t ) = ± 1 indicates that an ROI is either active (+1) or inactive (−1). First, the time series goes through a z-score normalization procedure, resulting in zero mean and unitary variance. To assess the sensitivity of our results to thresholding, we tested thresholds of 0 and ±1 SD. The results of this assessment will be presented in the section called Determining Parameters for Generating the Optimal Resting-State Structural Connectome. For the unconstrained pMEM, we fit the following probability distribution to all 76 subjects by maximizing a pseudolikelihood (see the Methods section): Pr(s) = exp (−βH(s))/Z, where H s ð Þ ¼ − P <i j> J i;j s i s j ; with i; j 2 1; 2; …; k ½ is the Hamiltonian function describing the energy of the system, and Z ¼ P Þis the partition function. Here, the spin configuration s is defined as the column vector s = [s 1 , s 2 , …, s N ] t max , where s i and s j are the spin states of region i and j, and J i,j represents a pairwise interaction between those regions. Traditionally, the Hamiltonian includes a term for external influences that we assume to be zero for resting-state data. We use the unconstrained pMEM as a control for comparison purposes. In our approach, we hypothesized that the interaction J i,j between two regions should be directly linked back to the diffusion MRI-derived structural connectivity between them as informed by tractography, so we add a constraint to the Hamiltonian function as follows: where W i,j is the structural connectivity between pairs of ROIs. This ensures that in the pseudolikelihood estimation of J, we constrain it with the structural connectivity (under the assumption that structural connectivity informs spin models governing brain dynamics). Thus, the optimal interaction matrix J is derived by maximizing the pseudolikelihood function as follows (Besag, 1975(Besag, , 1977: Pseudolikelihood substitutes Pr (s) by the product of the conditional probabilities e p ¼ Pr s i t ð Þj J;β; s −i t ð Þ À Á , observing one element s i (t) with all the other elements (denoted s −i(t ) ) fixed. To ensure that the magnitude of the coupling interactions is scaled relative to structural connectivity, the constraint is formulated as |J i,j | ∼ = μW i,j , where μ is a normalization constant and W i,j is the structural connectivity between ROI pairs. Without loss of generality, we assume that μ = 1 with appropriate normalization. We therefore present a penalty-based optimization scheme to maximize the constrained log-pseudolikelihood function as follows: The pseudolikelihood component 1 t max ln L pseudo J; β ð Þexpands to the following: Hamiltonian: Equation representing the total energy of the system.
We formulate the probability distribution based on the Boltzmann distribution under pseudolikelihood conditions. Thus, the numerator describes the energy of the system, while the denominator is the sum of all possible energies. Hence, there are only two terms in the denominator since s i (t ) is binary (one positive, and one negative). The likelihood function may be simplified by setting C i t ð Þ ¼ β P k m¼1 J i;m s m t ð Þ, resulting in the following formulation: Here we may construct the gradient ascent procedure with respect to J i,j by computing the partial derivative of the log-pseudolikelihood as the following: where A ¼ λ β . The updating scheme follows: J nþ1 Here, n is the iteration number and γ is the learning rate. In this way, the penalty function ensures that the inferred pairwise interaction is scaled relative to the estimated structure of the brain. This process is followed for all 76 subjects as part of first stage in constructing an optimized rs-SC as shown in the schematic of Figure 1. The next stage involves optimization of the two parameters that control scale and convergence, β and A.

Determining Parameters for Generating the Optimal Resting-State Structural Connectome
Within the framework for the FSE, parameters that need to be tuned are the constraint scale (parameter A) and the convergence parameter β. We note two points here: First, in the gradient ascent procedure, the influence of the β parameter is primarily in the hyperbolic tangent, which converges to 1 as β → ∞ and 0 as β → 0. Second, in a similar fashion, the influence of A is in the scale of the constraint. Thus, if as A → 0, then the model converges to an unconstrained pMEM, and if A → ∞, then the constraint will completely dominate the ascent and the system will converge to a pure structural connectome.
With that in mind, we develop two metrics to evaluate the accuracy and performance of our rs-SC, constructed with our FSE method. First is a similarity metric S m (β, A) that, as described in the Methods section, is simply the correlation between |rs-SC| and the structural connectome. This metric is used to gauge the quality of the constraint component in the framework. Second, we generate a correlation function f c (β, A) by simulating the Ising model with Markov chain Monte Carlo (MCMC) simulations (see Methods), computing a Pearson correlation between observed and simulated functional connectivity for all β simulated . As shown in Stage 2 of Figure 1, this results in a bell-like curve where each point represents a correlation value between the observed FC and reconstructed FC for each β simulated . Here we identify the max f c (β, A), which represents the maximum achieved correlation between observed and reconstructed FC for the parameters β, A. As described in the Methods section, we perform these computations for a range of empirically determined values. The results are presented in Figure 2. As expected, as the parameter values of β, A increase, the similarity metric S m evaluating the constraint increases to be almost perfectly correlated with the structural Boltzmann distribution: Quantifies the probability that a system will be in a specific state based on energy (of that state) and temperature (system).
connectome. However, the quality of FC reconstruction, evaluated using max f c (β, A), decreases as the values of β, A increase. This suggests that, while structural connectivity is important for reconstructing functional dynamics, a purely structural network is limited in the accuracy of FC reconstruction. Hence, using a grid-search technique shown in Stage 3 of Figure 1, we compute f (β, A) = max f c + S m , which intends to take equal weight between the underlying structure of the rs-SC and the accuracy with which it can reconstruct FC. As shown in Figure 2, this grid-search optimization results in the optimal values of β = 0.8 and A = 1.4 for a single subject (shown as an example). This process is computed individually for each subject, resulting in a unique parameter set for each individual. The range of parameter values is presented in Supplementary Figure 1, noting that the approximate mean value for A = 1.75 and β = 0.75. Using these parameters, we then reconstruct an optimized rs-SC for each of the 76 subjects to be used in the ensuing analyses.
Using this optimization process, we further tested different binarization strategies for the z-scored time series data to identify the optimal thresholding parameter. Here we tested 0, and ±1 SD. We used max f c (β, A) as a measure of performance to determine whether one threshold results in better FC reconstruction quality. The results are presented in Figure 3, noting that using zero as the threshold results in more consistent and accurate quality. We note the existence of a handful of outliers when binarizing ±1 SD. This is in part due to processing steps of the observed BOLD time series, where upwards of 40% of the TRs were excluded due to imaging artifacts for some subjects. As with all inference-based methods, accuracy of estimations increases or decreases with the amount of observed data. Future studies using this methodology will focus on cohorts with more consistent time series data across subjects. Figure 1. Schematic for the function-by-structure embedding (FSE) and ensuing parameter optimization strategy. The framework for constructing the hybrid resting-state structural connectome using the FSE is based on the principle of maximum entropy. Using a constrained maximum likelihood estimation where structural and functional connectivity are combined, we estimate both an edge strength in the network as well as a sign (±) representing excitatory or inhibitory interactions (Fortel et al., 2019). In Stage 2, two metrics are used to evaluate parameter quality, namely a similarity metric S m and the maximum of a functional correlation function max f c . As these values are dependent on parameter choices within the FSE framework, in Stage 3 a grid search is performed to find the optimal values for the tuned parameters that maximize f (β, A) = max f c + S m . These two metrics were computed for each subject individually, to identify the optimal parameters for constructing a hybrid resting-state structural connectome (rs-SC) for each subject.

Evaluating the Quality and Performance of the Resting-State Structural Connectome (rs-SC)
With the now optimized rs-SC it is important to compare this combined structure-function network with our control (or null model) pMEM-based interaction network, and the traditional Pearson-correlated FC. To do this we first perform MCMC simulations of the Ising model with 64,000 runs (N × N × 10; N = 80 ROIs) over a range of β simulated (see the Methods section). In this we may compute the max f c correlation results for each participant using the pMEM-based network and FSE-based networks (rs-SC), comparing the performance of both in reconstructing observed functional connectivity. In evaluating the results, we also separate our 76 subjects into two groups of noncarriers (NC) and carriers (APOE) to determine whether there are any within-group performance differences. The results of these simulations are presented in Figure 3, with the violin plot showing the median and range of max f c correlation values.
The median max f c for (NC, APOE) using the pMEM-based interaction network is (0.67, 0.60) with a range of {(0.84, 0.32), (0.85, 0.28)}. The 95% CI for the NC group is 0.66 ± 0.036 and 0.6 ± 0.044 for the APOE group. Using the FSE-based rs-SC, the median correlation for (NC, APOE) is (0.89, 0.90) with a range of {(0.94, 0.85), (0.94, 0.84)}. The 95% CI for the NC group is 0.89 ± 0.015 and 0.87 ± 0.021 for the APOE group. Based on comparison with the pMEM-based network, the FSE-based network results in a more consistent and accurate Here, β True is the value of β used in the FSE algorithm. In the top left is the average S m , the correlation r (|J |, W ) used to evaluate the performance of the constraint. The maximum value achieved in this example is r = 0.979 when β True , A are maximized. In the top right is the average of max f c , computed using MCMC simulations for the Ising model as described in the Methods section. The maximum value achieved is r = 0.91 for β True = 0.6, A = 0.2. Given the inverse effect of these two metrics, we compute f (β, A) = max f c + S m , identifying a parameter set that maximizes both metrics. Thus, in the grid search the maximum value is achieved at f (β, A) = 1.78 where β = 0.8 and A = 1.4. Last, in the bottom right is the average β simulated (max f c ) during MCMC simulations. We note that as the parameters β, A → 0, β simulated (max f c ) converge to the unconstrained pMEM, and as β, A → ∞, β simulated (max f c ) → 1. reconstruction of FC for both the NC and the APOE groups. Further, we evaluate the quality of our constraint under the optimized parameters β, A for all subjects. Violin plots for both NC and APOE groups, evaluating the similarity metric S m , are presented in Figure 3. The median S m for both groups is 0.94 with a range of (0.96, 0.90) in the NC group and (0.96, 0.88) in the APOE group. Further, the 95% CI for the NC group is 0.93 ± 0.0017 and 0.93 ± 0.0003 in the APOE group. These results suggest strong and consistent performance of the structurally Figure 3. Violin plots evaluating optimal hybrid resting-state connectomes (rs-SC), generated using the FSE framework. The top left plot evaluates the reconstruction quality, that is, the ability of our new network to reconstruct traditional FC correlation patterns. We use the network estimated using an unconstrained pMEM as a control for comparison. Using these networks, we identify a max f c value, representing the maximum correlation between the observed FC and reconstructed FC using MCMC simulations of the Ising model as described in the Methods section. Presented here are results for the noncarrier (NC) and APOE groups based on the two estimation techniques. The median max f c for (NC, APOE) using the pMEM is (0.67, 0.60) with a range of {(0.84, 0.32), (0.85, 0.28)}. The 95% CI for the NC group is 0.66 ± 0.036 and 0.6 ± 0.044 for the APOE group. Using the FSE, the median correlation for (NC, APOE) is (0.89, 0.90) with a range of {(0.94, 0.85), (0.94, 0.84)}. The 95% CI for the NC group is 0.89 ± 0.015 and 0.87 ± 0.021 for the APOE group. Last, given that the FSE framework relies on the structural connectivity as a constraint on the network estimation, we evaluate the quality of the constraint using the similarity metric S m described in the Methods section. In the top right plot, the median S m for both groups is 0.94 with a range of (0.96, 0.90) in the NC group and (0.96, 0.88) in the APOE group. Further, the 95% CI for the NC group is 0.93 ± 0.0017 and 0.93 ± 0.0003 in the APOE group. These results suggest strong and consistent performance of the structurally informed rs-SC in reconstructing functional dynamics for both groups, as well as constraint on the estimated network. Last, the bottom plot presents an evaluation of binarization thresholds using the max f c value, representing the maximum correlation between the observed FC and reconstructed FC using MCMC simulations of the Ising model. The median value when binarizing about zero is 0.89, while the median value is 0.62 and 0.60 for ±1 SD, respectively. These results reveal outliers when binarizing the time series data using a value other than zero. This is most likely due to data quality questions related to fMRI processing resulting in time series with up to 40% of the TRs excluded due to imaging artifacts for some subjects.
informed rs-SC in reconstructing functional dynamics for both groups, as well as constraint on the estimated network.
Last, in previous studies no group differences could be identified between NC and APOE groups using traditional FC measures (Fortel et al., 2019;Korthauer et al., 2018). Here, we test whether an excitation-inhibition (E/I) ratio may be used to differentiate between the rs-SC network as well as the pMEM-derived interaction network. As defined in the Methods section, the E/I ratio is simply the sum of the positive edges, divided by the sum of negative edges (computed at either the whole-brain or ROI level). For both networks, the E/I ratio is computed for all ROIs and averaged for the NC and APOE groups. In Figure 4, the scatterplots display a weak association between the average E/I ratio for all ROIs between the NC and APOE groups using the pMEM-derived network with R 2 = 0.44, and a paired t test across all ROIs results in P >> 0.1, suggesting no statistically significant group differences in the pMEM-based networks. Conversely, performing a similar computation of the E/I ratio on the rs-SC networks results in a strong association between NC and APOE groups, with R 2 = 0.76, as well as a notable shift observed for all ROIs (i.e., globally). A paired t test between the groups results in P = 0.037, suggesting a statistically significant group difference in the rs-SC networks between NC and APOE groups, as evaluated with the E/I ratio that cannot be identified with the unconstrained model. In sum, the results presented in this section indicate that the novel rs-SC network constructed with the FSE framework can not only describe structural and functional dynamics, but also probe brain dynamics that may not be captured using a similar unconstrained methodology. Last, we compare the E/I ratio in a group comparison of males and females (NC versus APOE) and present the results in the scatterplot of Figure 5. Both males and females exhibit a positive association between groups with R 2 = 0.65 and R 2 = 0.75, respectively; however, only the female group has a statistically significant group . Group comparison of the excitation-inhibition ratio for each brain region based on the unconstrained pairwise maximum entropy model and the function-by-structure embedding. As described in the Methods section, the E/I ratio is simply the sum of positive edges divided by the sum of negative edges for each ROI. Here, we present a plot comparing the E/I ratio between the NC and APOE groups using the pMEMbased network and our rs-SC network, computed and averaged at the ROI level. This results in a weak association with R 2 = 0.44 for for the pMEM-based network, and R 2 = 0.76 for the rs-SC network, with paired t tests across all ROIs resulting in P >> 0.1 and P = 0.037, respectively. This suggests no statistically significant differences in E/I balance when using the unconstrained model; however, there is a statistically significant difference in the two groups when using our structurally informed model. We note that numerically, an increase in group-averaged E/I ratio would move a point (representing one ROI) above the x = y reference line, suggesting a shift in E/I balance towards hyperexcitation. A tabular version of these results is included in the Supporting Information. difference with P = 0.008. Here, we perform a calculation on group difference by computing to evaluate the average change between NC and APOE groups. The males have an average increase of 2.4% (averaged across all ROIs), while the females have an average increase of 5.9% in E/I ratio (approximately 2.4x higher than the male group).
The raw values for each brain region are shared in Supplementary Table 1 for reference.

Criticality and Hyperexcitation in Female APOE-ε4 Carriers
In this study, our subjects are separated into two age-and sex-matched groups (NC and APOE). One aspect of the link between APOE-ε4 and AD that has often been overlooked is that Figure 5. Gender-based comparison of critical behavior and E/I balance. As described in the Methods section, the E/I ratio is simply the sum of positive edges divided by the sum of negative edges for each ROI. In the top panels, we present plots comparing the E/I ratios between the NC and APOE groups for males and females, computed and averaged at the ROI level. This results in a strong association with R 2 = 0.65 for males, and R 2 = 0.75 for females with paired t tests across all ROIs resulting in P = 0.19 for males, and P = 0.008 for females. This suggests no statistically significant differences in E/I balance for males; however, there is a statistically significant difference for females. We note that numerically, an increase in group-averaged E/I ratio would move a point (representing one ROI) above the x = y reference line, suggesting a shift in E/I balance towards hyperexcitation with increased risk of chaotic activity. Thus, for each ROI, we can quantify the shift in E/I balance by computing delta to evaluate the average change between NC and APOE groups; this yields a shift of 5.94% in the female group between carriers and noncarriers, while in the male group it is 2.46% (approximately 2.4× difference between sexes). A tabular version of these results is included in the Supporting Information.
females with at least one ε4 allele are four times more likely to develop AD than males (Bretsky et al., 1999;Jack et al., 2015;Payami et al., 1994). Thus, we use our framework to evaluate not just group differences in criticality, but sex differences as well (22F/16M in each group). As previously mentioned, the brain criticality hypothesis suggests that neural networks selforganize into a unique configuration between order and disorder. In the context of statistical physics and the Ising model, this unique configuration occurs at some critical point (β critical ). Here, we again utilized MCMC simulations to generate a series of state configurations (±1) resulting in an N × t matrix, where N = 80 ROIs, and t = 100,000 runs (see Methods). In the previous section we used these states to compute a correlation between brain regions; however, in this case we will evaluate the critical dynamics elucidated from the rs-SC networks. Specifically, we are interested in the phase transitions based on the positive edges of the networks. The Ising model can be modified to model spin-glass behavior (full signed network); however, this can lead to "frustration" in the simulations. Frustration describes a scenario in which it is impossible to simultaneously minimize all the terms in the Hamiltonian. As a result, this generally leads to complex energy landscapes with many local minima. At low β, the system can get stuck in the local minima without ever reaching a true equilibrium. In future work we can investigate thermodynamic properties using the full signed network (spin-glass), but here we proceed in evaluating the ferromagnetic phase transitions.
For each β, we compute the order parameter (magnetization) and the variance (susceptibility) with respect to β simulated . Here, we again compute average values over NC and APOE groups to investigate potential group differences in critical behavior. As described in the methods, β is the inverse temperature (T ) parameter used in the Boltzmann distribution, and thus when simulating dynamics to identify phase transitions, we interpret temperature as a tolerance of the system when increased randomness is introduced. Performing Monte Carlo simulations of the Ising model using our hybrid network for a range of temperatures is used to identify a critical point, such that the system transitions from a hypoactive regime to a chaotic regime. Hence, the critical temperature is a measure of how much tolerance the system has to increased perturbations. We present the phase diagrams for susceptibility in Figure 5 for males and females, highlighting T critical for both groups (evaluated by the peak of susceptibility). It should be noted that β was simulated from 0.2 to 3.0 at increments of 0.05 (then plotted against T ¼ 1 β Þ. We identified a more pronounced deviation between NC and APOE females with T critical = 0.65 for the female APOE group as compared with T critical = 0.87 in the NC group. Conversely, T critical = 0.80 in the male APOE group as compared with T critical = 0.83 in the NC group. This suggests that the critical dynamics within the male group between NC and APOE are more similar in nature than the dynamics observed within the female group between NC and APOE. A lower critical temperature in the female carrier group suggests a lower tolerance to network dysfunction as a result of an increase in excitatory interactions, increasing vulnerability to chaotic activity. In sum, these results suggest that there is a link between brain criticality and excitation-inhibition balance that can be identified via our new connectome, demonstrating a disruption to this balance in APOE carriers (with a larger effect in females).
Further, presented here are plots demonstrating a global evaluation of critical brain dynamics. In the bottom panels, ferromagnetic susceptibility is shown for males and females, with the dashed lines representing the critical point T critical = T simulated (max χ). These charts demonstrate a more pronounced deviation between NC and APOE females with T critical = 0.65 for the female APOE group as compared with T critical = 0.87 in the NC group. Conversely, T critical = 0.80 in the male APOE group as compared with T critical = 0.83 in the NC group. This suggests that as the E/I balance shifts at global scale, the critical point also decreases because of an Critical temperature: Parameter describing the tolerance of a system to increasing perturbations (commonly used in physics to identify loss of magnetic properties). increase in excitatory interactions. As described in the Methods section, a lower critical temperature indicates a lower tolerance to network dysfunction, increasing vulnerability to chaotic activity.

DISCUSSION
Using a constrained maximum entropy model for our function-by-structure embedding (FSE), we have developed here a novel resting-state structural connectome (rs-SC), unifying connectome-level structure and function into a new spatiotemporal network. We constructed rs-SC networks for 76 cognitively intact participants with a grid-search parameter optimization scheme. Hence, we demonstrate two important results: First, the underlying structure of the rs-SC is as expected, strongly correlated with the empirical structural connectome (r > 0.9) due to it being used as a constraint in the FSE framework. Second, and more importantly, we demonstrated that it is possible to model the resting-state functional connectome based on a model of spin products, accounting for indirect or higher order structural connectivity. We acknowledge that when Ising dynamics are used to model neural firing patterns, these activations may amount to the collective behavior of a few neurons, and at the macro level of fMRI imaging used in this study each voxel may be providing information as a result of thousands of interacting neurons. However, simulation and empirical studies have demonstrated that increases in excitatory neuronal activity amplify oscillations associated with the transient BOLD response, while increasing inhibitory activity evokes an overall decrease in the BOLD signal (Aksenov et al., 2019;Krishnan et al., 2018;Sotero & Trujillo-Barreto, 2007;Sten et al., 2017). By grounding our macroscale methodology with models of microscale dynamics, we bridge the gap between the two, hereby inferring the nature (excitatory or inhibitory) of structural connectivity at rest. Further, the rs-SC can be used to simulate functional dynamics using Monte Carlo simulations, reconstructing traditional functional correlation patterns (r avg = 0.9). Beyond model quality and performance, we have also demonstrated that our rs-SC can distinguish between female noncarriers and APOE-ε4 carriers (age-and sex-matched) using our excitation-inhibition (E/I) ratio. Our results demonstrate that modeling with the rs-SC reveals a global shift of E/I balance for the APOE-ε4 carrier group. Given that APOE-ε4 carriers are at an elevated risk for AD, the observed shift in E/I balance in this sample may be a result of disease pathology. In many studies of AD, one critical feature that is often overlooked is that females with at least one ε4 allele are four times more likely to develop AD than are males (Jack et al., 2015). A comparison of group-averaged E/I ratio at the ROI level for each sex using the rs-SC (with new optimization strategy) yielded a global shift in E/I balance towards hyperexcitation, in line with our previous work (Fortel et al., 2020) and prior studies on sex differences related to the APOE genotype (Aboud et al., 2013;Bi et al., 2020;Jiménez-Balado & Eich, 2021;Leung et al., 2012). In future work, we may investigate in depth the relationship of our hybrid connectome with traditional measures of structural and functional connectivity in a larger cohort (with increased age range), to investigate known sex differences and further evaluate our method.
Further, in this study, we observe significant differences in critical behavior between a group of cognitively intact individuals with a genetic predisposition for late onset AD as compared with age-and sex-matched noncarriers. Traditional structural and functional connectivity based on BOLD correlations were unable to separate the two groups (Fortel et al., 2020;Korthauer et al., 2018). These results suggest that using a multimodal framework to unify structure and function can reveal underlying patters in brain dynamics that would otherwise not be captured using traditional methods. Further, we endeavored to identify a link between E/I balance and criticality. As a result of increased positive interactions (increased deviation from an E/I balance) in the hybrid connectome, simulations of brain dynamics using Monte Carlo simulations revealed a shift in criticality for female carriers compared with noncarriers of APOE-ε4 that may suggest an increased vulnerability to AD neuropathology in female APOE-ε4 carriers. We describe the critical temperature as a measure of tolerance in our modeled system that we simulate in dynamical regimes spanning from highly ordered (i.e., hypoactive) to highly disordered. This is in line with studies of preclinical neural models that have shown that networks operating at criticality exhibit an E/I balance as compared with networks that have been over excited or overinhibited by a controlled chemical stimulus (Heiney et al., 2019;Shew et al., 2011). In fact, many of the in vivo studies that have investigated the criticality hypothesis and excitation-inhibition balance in neurodegenerative disorders have relied on electroencephalography (EEG) or magnetoencephalography (MEG) recordings (Bruining et al., 2020;Montez et al., 2009;Rajkumar et al., 2021;Stam et al., 2005), which have inherent challenges with spatial resolution. By defining our activity states using both structural and functional connectivity together, we are capable of analyzing patterns of activity across both temporal and spatial scales, thereby improving the network inference and mitigating many challenges observed in unimodal and traditional analyses.
The results presented herein regarding E/I balance, criticality, and the APOE-ε4 genotype also coincide with the current understanding of the microscale mechanisms underlying AD pathology. A recent review article by Najm and colleagues explored the relationship among APOE-ε4, loss of GABAergic interneurons, and dysfunctional brain networks in the context of AD (Najm et al., 2019). In short, neurons responding to different factors (e.g., normal aging, injury, or stress) break down APOE-ε4 proteins and produce fragments that trigger phosphorylation of tau; this in turn disrupts mitochondrial function, leading to cell death. Destruction of inhibitory neurons in this way can alter network activity and produce hyperexcitability in neural circuits long before clinically identifiable symptoms arise. This may help explain the known associations of APOE-ε4 with memory deficits and severe epilepsy. Indeed, several in vitro and preclinical in vivo studies (cited by Najm et al., 2019) have demonstrated that intracellular APOE-ε4 is toxic to GABAergic interneurons, particularly in the hippocampus.
Moreover, other authors have recently suggested that neuronal hyperexcitability may be considered to be both a causal factor and a risk factor in the disease progression, even in the preclinical phase (Hijazi et al., 2020;Paterno et al., 2020;Tok et al., 2021). While significant structural and functional degeneration is well established in AD (DeTure & Dickson, 2019), our framework incorporates both structural and functional connectivity in order to provide a new multimodal perspective of connectome-level interactions in a preclinical group of individuals predisposed to AD. We acknowledge that our methodology is limited to insights that may be gained from macroscale BOLD activity as opposed to direct measurements of neuronal processes. That said, we reached a similar conclusion to independent studies of underlying neural mechanisms in AD: Individuals with the APOE-ε4 allele (females in particular) have a higher risk of neurodegeneration due to an increase of excitatory activity in neural circuits (Jiménez-Balado & Eich, 2021;Koutsodendris et al., 2022;Y. Li et al., 2016).
We note several limitations of this study. First, this study investigated only a small cross section of healthy middle-aged individuals at increased risk of developing AD. Further, the parcellation used in the processing used an atlas with 80 brain regions, which may be considered too coarse. Additional research with a longitudinal cohort and higher resolution parcellation would help improve the generalizability of results, providing important validation regarding within-subject variability, as well as broadening our understanding of longitudinal alterations in brain dynamics. Second, when interpreted as a strictly nodal property, excitationinhibition balance may be best measured at a regional level using FDG PET or phosphorous imaging. However, as conceptualized in this study, the concept of E/I balance may directly relate to this notion of "criticality" in brain dynamics. Further, in this group of participants, measurements of well-known biomarkers of Aβ and tau were not included in the protocol, and thus we could not add this layer of validation. Future studies comparing additional imaging modalities and biomarkers for validation and correlation purposes may be used to strengthen the results and methodology presented in this study (in addition to more state-ofthe-art diffusion tensor imaging and fMRI imaging protocols). Further, in this study as we are working with resting-state data processed with global signal regression (accounting for background and nonneural physiological noise), we model the BOLD activity assuming no external influences; future work can incorporate external influences in the framework to account for different interference scenarios.
It remains unclear whether the difference in criticality observed between the NC and APOE groups is because the NC group (on average) contains more inhibitory interactions or whether the APOE group has more excitatory interactions. Since we do not identify directionality in this study, this question is left for future work. Additionally, we have not performed an assessment herein on the potential relationships between traditional structural and functional connectivity measures, and metrics obtained with our rs-SC. This may be explored in detail with future investigations. Further, at the coarser spatial scale of human fMRI, there is evidence that the strength of functional connectivity between regions is greatest for region pairs separated by short physical distance and that connectivity strength decays rapidly as the Euclidean distance between brain regions increases (Alexander-Bloch et al., 2013). Likewise, the extent of white matter tract connectivity as measured with diffusion imaging also decays with distance. However, the inverse relationship between fMRI-based connectivity and distance is significant even after controlling for the strong association between anatomical connectivity and functional connectivity (Honey et al., 2009). In the future, the role of distance related to excitatory and inhibitory interactions should be explored in greater depth. Further utilizing thermodynamic principles, it should be investigated whether the rs-SC decays algebraically with a distance d (i.e., J(d ) ∝ d −α ) as well as what, if any, effect this distance decay would have on critical brain dynamics. Given the complex inner workings of the brain, it is entirely plausible that dynamics between brain regions at or near criticality rely on a balance between long-and short-range interactions. Again, this suggests that functional brain dynamics are governed by the underlying structure of the networks. Thus, after decades of research studying the brain's individual components, from neurons to neuronal ensembles and large-scale brain regions, conclusive evidence demonstrates the need for maps and models that incorporate interactions among these components in order to better understand the brain's ensemble dynamics, circuit function, and emergent behavior.

Participants and MRI Data Acquisition
The cohort used in this work has been described in a previous study (Korthauer et al., 2018). Participants (N = 76; all Caucasian) were selected based on APOE genotype from a larger sample of 150 adults aged 40-60 (age = 49.9 ± 6.0 in years; 60 men). The University of Wisconsin-Madison Biotechnology Center conducted the sequencing of the single nucleotide polymorphisms (SNPs; rs7412, rs429358) making up the common ε2, ε3, and ε4 APOE genotypes. Thirty-eight individuals out of the larger sample were APOE-ε4 carriers (either ε3/ε4 or ε4/ε4). Hence, a subset of noncarriers (ε3/ε3 or e2/ε3) were age-and sex-matched, creating equal groups (N = 38, 22 female) of carriers (APOE) and noncarriers (NC). The following exclusion criteria were used: (a) self-reported cognitive or memory complaints; (b) Mini-Mental State Exam (MMSE; Folstein et al., 1975) score ≤ 24; (c) Mattis Dementia Rating Scale Second Edition (DRS-2; Johnson-Greene, 2004) score ≤ 135; (d) Geriatric Depression Scale (GDS; Yesavage et al., 1982) > 10; (e) history of central nervous system disease (e.g., dementia, stroke, Parkinson's disease, epilepsy, other neurological disease); (f ) history of severe cardiac disease (e.g., myocardial infarction, coronary bypass surgery, angioplasty); (g) history of metastatic cancer; (h) history of serious psychiatric disorder or substance use disorder; and (i) any contraindication to MRI. MRI imaging was conducted on a GE Signa 3T scanner ( Waukesha, WI) with quad split quadrature transmit/receive head coil. All participants provided written informed consent, and were compensated financially for their participation; the imaging collection was carried out in accordance with the guidelines set by the institutional review boards of the University of Wisconsin-Milwaukee and Medical College of Wisconsin (Korthauer et al., 2018). Demographic characteristics and screening measures for each group are presented in Table 1.

Processing of fMRI and Diffusion Tensor Imaging
Preprocessing of rs-fMRI images was performed using Analysis of Functional NeuroImages (AFNI; Cox, 1996) and FMRIB Software Library (FSL; Smith et al., 2004) based on the rs-fMRI preprocessing pipeline from the Human Connectome Project (HCP; Smith et al., 2013). Detailed processing information steps can be found in prior work (Korthauer et al., 2018). Diffusion tensor imaging data processing was carried out with the FSL. The B0 image was skull-stripped using the brain extraction tool (Smith, 2002), with the resulting mask applied to the other images. Eddy current-induced distortions and subject movements were corrected using FSL's eddy tool (Andersson & Sotiropoulos, 2016). A probability distribution for fiber direction was generated at each voxel using BEDPOSTX (Behrens et al., 2007;Behrens et al., 2003), which was then used in probabilistic tractography. For individual subjects, FreeSurfer cortical parcellation and subcortical segmentation was used, defining the 80 ROIs (Dale et al., 1999;Fischl et al., 2002;Fischl et al., 2004). Affine registration with 6 degrees of freedom using FLIRT registered the ROIs to MNI and diffusion space (Jenkinson et al., 2002). For each ROI, the mean time course from the BOLD signal was extracted using global signal regression (GSR) from the preprocessed rs-fMRI data prior to constructing the functional connectivity matrix. The resulting zero-mean time courses for each ROI were then correlated using Pearson correlations to generate a traditional functional connectivity matrix. Probabilistic tractography was performed between pairs of ROIs using Probtrackx for estimating the structural connectivity. The resulting matrix was then further normalized by dividing each matrix row by the way-total for its corresponding seed ROI (Behrens et al., 2007;Behrens et al., 2003).

The Unconstrained Pairwise Maximum Entropy Model (pMEM)
This maximum entropy approach provides a way of quantifying the goodness of fit in models that include varying degrees of correlations (Schneidman et al., 2006). At a microscale level for example, a first-order model seeks to fit only the average firing rate of all neurons recorded in the ensemble. A second-order model would seek to fit the average firing rate and all pairwise correlations, with an nth-order model fitting all correlations up to and including those between all n-tuples of neurons in the ensemble. At macroscale, this amounts to fitting the average BOLD activation rate of a brain region and all pairwise correlations. Here, the observed bold activation rate is determined through a binarization of the BOLD time course. Thus, we construct unbiased predictions for the probabilities of functional brain states by fitting a pairwise maximum entropy model (pMEM). Here, in estimating the probability distribution, it is necessary to use the distribution that maximizes the uncertainty (e.g., entropy). To fit the pMEM, we must tune the first-and second-order interaction parameters between ROIs such that the predicted activation and coactivation rates match the observed data (the BOLD time series). An accurately fitted pMEM suggests that patterns of functional activity can be estimated from each ROI's independent activation rate combined with the joint activation rates. Thus, the pMEM represents a model of fMRI BOLD dynamics as a probabilistic process defined by underlying pairwise relationships between ROIs. In constructing this model, we leverage the Ising model, a special case of a Markov random field in which each ROI can exhibit two possible states s = ± 1. In this work, we first convert our BOLD time series to z-scores, ensuring that our BOLD date is represented as zero-mean with unitary variance, without altering the correlations between brain regions. As maximum entropy models of neural activity are developed based on Ising dynamics, studies investigating pairwise interactions using BOLD time course data are binarized to define activation states (either +1 for active, or −1 for inactive) in both simulated and empirical fMRI-based studies (Ashourvan et al., 2021;Cofré et al., 2019;Ezaki et al., 2020;Ezaki et al., 2017;Gu et al., 2018;Nghiem et al., 2018;Niu et al., 2019;Watanabe et al., 2013). We will show how the binarization strategy may be validated using Monte Carlo simulations, whereby we use the inferred interaction networks to reconstruct functional correlations. Our results will also show that for our network construction methodology, binarizing the z-scored time series at zero provides better inference of functional interactions than ±1 SD.
We first begin by modeling the neural system using an energy-based formulation, namely the Hamiltonian, as follows: Here, the spin configuration s is defined as the column vector s = [s 1 , s 2 , …s k ] T , k is the number of regions, s i and s j are the spin states of region i and j, and J i,j represents a pairwise interaction between ROIs. Conceptually, if two regions are co-active or co-inactive, the pairwise interaction is likely positive (excitatory), and if one region is active while the other is inactive, the pairwise interaction is likely negative (inhibitory). Here we assume that there is no external influence (i.e., resting state). Further, unless otherwise stated, the summations in this manuscript are for i < j to avoid double counting and exclude self-connections. The probability of observing a specific configuration is given as the following Boltzmann distribution: where β is the inverse temperature, and Z is the partition function: The summation in the partition function is over all possible configurations of states. Similar to other studies fitting pairwise models to neuronal firing data, a gradient ascent updating scheme is used (Watanabe et al., 2013;Yeh et al., 2010). Estimating a parameter set that minimizes the Kullback-Leibler (K-L) divergence between modeled and observed probability distributions is equivalent to maximizing a log likelihood of the observed data (the empirical BOLD time series). We note that a brute-force application of the maximum likelihood estimation requires heavy computational costs with calculations over all 2 N possible spin configurations for the partition function (Nguyen et al., 2017). To overcome the intractability of the partition function Z, we utilize a pseudolikelihood estimation method (Ezaki et al., 2017). Pseudolikelihood estimation has been shown to converge to a maximum likelihood estimator for large sample sizes (Besag, 1975).
The optimal interaction matrix J can thus be derived by maximizing the pseudolikelihood function (Besag, 1975(Besag, , 1977: Pseudolikelihood substitutes the probability of observing the state vector s(t) by the product of the conditional probability e p ¼ Pr s i t ð Þj J; β; s −i t ð Þ À Á of observing a single element s i (t ) while all the other elements, denoted s −i(t) , are fixed. Thus, we maximize the following logpseudolikelihood function as the following: This probability distribution is derived based on the Boltzmann distribution under pseudolikelihood conditions. The numerator describes the energy of the system, while the denominator is the sum of all possible energies. Hence, there are only two terms in the denominator, one positive and one negative since s i (t ) is binary. The likelihood function may be simplified further by setting C i t ð Þ ¼ β P k m¼1 J i;m s m t ð Þ, resulting in the following: The gradient ascent procedure can now be constructed with respect to J i,j by computing the partial derivative of the log-pseudolikelihood as the following: The updating scheme follows: J nþ1 Here, n is the iteration number and γ is the learning rate.

Monte Carlo Simulations for the Ising Model
All scripts were developed and executed in MATLAB R2018a on a Windows 10 machine with Intel i7 CPU@ 2.8 GHz and 16 GB of RAM. We used a Markov chain Monte Carlo (MCMC) method based on the metropolis algorithm to calculate the observables of the Ising model using the networks inferred from the pMEM and FSE. Here we present the simulations performed step by step: (1) Define the parameters J (network inferred with pMEM or FSE), the number of runs t, and a range of β simulated .
(2) For each run randomly fix an s i from the configuration and compute the Hamiltonian , flip the state. Note: The command rand(0, 1) generates a random value between 0 and 1. Complete this for all elements in the configuration.
(4) The final configuration of states is then used as the input for the next run.
(5) Concatenate all runs into an N × t array and compute the averages of the observables (i.e., Pearson correlation < s i s j >, magnetization |M|, susceptibility χ). Because of the computational cost, when performing MCMC simulations for the grid-search parameter optimization we used t = 2,000 runs and β simulated from 0.2 to 3.0 with increments of 0.2. For the control case based on pMEM, we used t = N × N × 10 runs with β simulated from 1 to 20 with increments of 0.5. Last, when evaluating the thermodynamic properties magnetization |M|, susceptibility χ using the rs-SC network, we use t = 100,000 with β simulated from 0.2 to 3.0 with increments of 0.05. The number of runs as well as range and increments of β simulated were selected based on the task performed to maximize algorithmic performance and to minimize processing time. The upper and lower bound of these values was first empirically determined to be containing the optimal range by simulations.

Phase Transition and Biological Motivation
The simplicity of the Ising model enables the prediction of cooperative behavior among a system of biological elements wherein each element has two states, and the energy of the system depends only on the state of each element and its neighbors. Moreover, the model parameters and representative physical properties are readily amenable to biological interpretation in the context of various complex systems. For example, a four-dimensional cellular automaton-like Ising model has been previously developed to investigate transitions between normal, proliferative, hypoxic, and necrotic states in the tumorigenesis processes (Durrett, 2013;Torquato, 2011). Ising-like models have also been implemented to estimate information transfer between spins occurring on the human connectome (Marinazzo et al., 2014) or to assess differentially expressed genes in cancer patients (X. Li et al., 2011), and even to model the joint expression profiles of genes to reconstruct E. coli gene interaction pathways (Santhanam et al., 2009). Hence, when we discuss a "phase transition," it is a result of the interactions among many elements, not from the specific nature of the individual units (be they ferromagnetic materials or biological elements like neurons, protein chains, or genes).
To evaluate these transitions, we look to the average of activations over the whole network (termed magnetization), which determines the ordering of the system. Magnetic susceptibility is simply the variance of the magnetization. If all the binary spin states are aligned in the same direction, a magnetization of ±1 corresponds to a configuration of complete order. The magnetization per site, defined as M = P N i¼1 < s i >; where < : > represents the ensemble average and quantifies the mean tendency that s i = 1 as opposed to s i = −1, is taken across the brain regions. The magnetic susceptibility is defined as χ ¼ 1 β <M 2 > − <M> 2 ð Þ (Landau, 2009).
Here, we consider brain networks positioned near a critical point between complete inactivity (i.e., neuronal death) and random activity (as in epilepsy, for example). In a less extreme sense, simulations of Ising dynamics can reveal a transition from a hypoactive state towards a more chaotic state. As described in Equation 9, the behavior of the modeled system depends on temperature. However, for a network of neurons or brain regions, there is no real concept of "temperature." Hence, when performing Monte Carlo simulations of the Ising model, we may describe temperature (T ) as a "tolerance" of the system in the sense that the effect of the T parameter injects additional randomness into the simulated dynamics of the system. Thus, for very low T (T < T critical ), spontaneous MCMC spin flips are less probable, with the spins in each configuration mostly aligned to contribute the minimum energy of the system. For very high T (T > T critical ), the magnetic ordering is completely lost as a result of a high number of spontaneous spin flips; thus, the magnetization tends to zero, which can be used to characterize the disordered (or chaotic) phase. In the intermediate range of T where selforganized criticality and second-order phase transitions occur, there is a point of maximal fluctuations in the magnetization at T = T critical that corresponds to a peak in the magnetic susceptibility (Chialvo, 2010). Thus, a system with lower critical temperature is suggestive of a lower tolerance to perturbations in the network as determined via Monte Carlo simulations of brain dynamics than is a higher critical temperature that would suggest a higher tolerance.

Parameter Optimization Using a Similarity Metric and Correlation Function
In this work, we use a grid-search optimization scheme to find the optimal parameters {β, A}. The parameters are evaluated from 0.2 to 3.0 with 0.2 increments for all 76 participants. With the FSE, J, we generate a correlation function max f c ( β, A) by simulating the Ising model with Monte Carlo simulations, computing a Pearson correlation between observed and simulated functional connectivity for all β simulated (from 0.2 to 3.0 with 0.2 increments). Further, we compute a similarity metric S m ( β, A) via the correlation r (|J i,j |, W i,j ) ∀ i, j to ensure that jJ i,j j ∝ W i,j , the structural connectome. To identify the optimal parameters, we find A, β such that f ( β, A) = max f c + S m is maximized.

Excitation-Inhibition (E/I) Ratio
It is important to note that in using the terminology connectome-level excitation-inhibition balance and hyperexcitation, we are not necessarily inferring directionality of these interactions or measuring processes at a neuronal level. Rather, we used such terminology to bridge the gap between microscale interactions (such as excitation and inhibition of neuronal circuits) and the connectome-level changes that may occur because of such processes. Note that similar terminologies have previously been adopted in several seminal studies that investigated neuronal firing patterns using the Ising model (Schneidman et al., 2006;Tkačik et al., 2013). To be clear, from a connectomics perspective, if several brain regions are identified to have an increase in positive edges in the rs-SC, collectively, that would suggest a wider spread pattern of coupling (i.e., more likely to exhibit a pattern of global coupling) that may subserve hyperexcitation. It is in this context that we conceptualize the excitation-inhibition (E/I) ratio, a global (whole-brain) or local (ROI-level) estimation of E/I balance, computed as the sum of positive edges divided by the sum of negative edges. For example, if an ROI in the network has 45 positive edges and 34 negative edges, then the E/I ratio = 45 34 , or 1.32 (a value of 1 indicates perfect E/I balance).

SUPPORTING INFORMATION
Supporting information for this article is available at https://doi.org/10.1162/netn_a_00220.