## 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 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 sex-matched 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.

## RESULTS

### Constructing a Function-by-Structure Embedding Using a Constrained Maximum Likelihood Estimation

*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*) ∈ {−1, +1}

^{N}. Note that

*t*

_{max}is determined as a result of the 8-min fMRI 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*(

**) = exp (−**

*s**βH*(

*s*))/

*Z*, where

*H*(

**) = −∑**

*s*_{<ij>}

*J*

_{i,j}

*s*

_{i}

*s*

_{j}, with

*i*,

*j*∈ [1, 2, …,

*k*] is the Hamiltonian function describing the energy of the system, and

*Z*= ∑

_{s}exp(−

*βH*(

**)) is the partition function. Here, the spin configuration**

*s***is defined as the column vector**

*s***s**= [

*s*

_{1},

*s*

_{2}, …,

*s*

_{N}]

^{tmax}, 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:

*W*

_{i,j}is the structural connectivity between pairs of ROIs. This ensures that in the pseudolikelihood estimation of

**, 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, 1977):**

*J**Pr*(

**) by the product of the conditional probabilities $p\u02dc$ =**

*s**Pr*(

*s*

_{i}(

*t*)|

**,**

*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:

_{pseudo}(

**,**

*J**β*) expands to the following:

*s*

_{i}(

*t*) is binary (one positive, and one negative). The likelihood function may be simplified by setting

*C*

_{i}(

*t*) = β $\u2211m=1k$

*J*

_{i,m}

*s*

_{m}(

*t*), resulting in the following formulation:

*J*

_{i,j}by computing the partial derivative of the log-pseudolikelihood as the following:

*A*= $\lambda \beta $. The updating scheme follows: $Ji,jn+1$ = $Ji,jn$ +

*γ*$\u2202\u2113\u2202Ji,j$|

_{n}. 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
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.

### 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
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 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 difference with *P* = 0.008. Here, we perform a
calculation on group difference by computing *delta* = 1 − $E/INCE/IAPOE$ 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 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 self-organize 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\beta $).
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 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, excitation-inhibition 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-of-the-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.

## MATERIALS AND METHODS

### 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.

**Table 1.**

. | ε4 carriers (N = 38)
. | non-ε4 carriers (N = 38)
. |
---|---|---|

Age (years) | 50.8 (0.99) | 50.9 (0.99) |

Sex (M:F) | 16:22 | 16:22 |

Education (years) | 15.4 (2.5) | 15.2 (2.4) |

DRS-2 (total) | 139.9 (2.3) | 139.9 (2.3) |

MMSE (total) | 28.5 (1.1) | 28.8 (1.3) |

GDS (total) | 1.8 (2.3) | 2.4 (2.7) |

. | ε4 carriers (N = 38)
. | non-ε4 carriers (N = 38)
. |
---|---|---|

Age (years) | 50.8 (0.99) | 50.9 (0.99) |

Sex (M:F) | 16:22 | 16:22 |

Education (years) | 15.4 (2.5) | 15.2 (2.4) |

DRS-2 (total) | 139.9 (2.3) | 139.9 (2.3) |

MMSE (total) | 28.5 (1.1) | 28.8 (1.3) |

GDS (total) | 1.8 (2.3) | 2.4 (2.7) |

All participants were screened for any contraindications to MRI. Imaging sessions
lasted 75 min. To determine the structural and functional connectivity maps,
multimodal imaging (including T1-weighted MRI, resting-state fMRI, and diffusion
weighted MRI) was performed. For structural MRI imaging, a
“spoiled-grass” (SPGR) sequence (axial acquisition: TR = 35 ms, TE
= 5 ms, flip angle = 45°, matrix = 256 × 256, FOV = 24 cm, NEX =
1) was obtained, followed by a T2*-weighted functional scan with an
echo-planar pulse imaging (EPI) sequence (28 axial slices, 20 × 20
cm^{2} FOV, 64 × 64 matrix, 3.125 mm × 3.125 mm
× 4 mm voxels, TE = 40 ms, TR = 2,000 ms). The 8-min rs-fMRI scan was
acquired while participants were under task-free conditions (i.e., resting
state). Additionally, a 3-min, 30-s diffusion tensor imaging sequence was
acquired with a spin echo single shot, echo-planar imaging sequence with
sensitivity (SENSE = 2.5) encoding (2.2 mm isotropic voxels, 212 × 212 mm
FOV, 96 × 96 acquired matrix), TR/TE = 6,338/69 ms, 60 slices for
whole-brain coverage. Diffusion gradients were applied along 32 noncollinear
directions at a b-factor of 700 s/mm^{2}, including one minimally
weighted image with b = 0 s/mm^{2}.

### 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 *n*th-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*.

*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:

*β*is the inverse temperature, and

*Z*is the partition function:

*Z*= ∑

_{s}exp(−

*βH*(

**)).**

*s*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).

**can thus be derived by maximizing the pseudolikelihood function (Besag, 1975, 1977):**

*J***s**(t) by the product of the conditional probability $p\u02dc$ =

*Pr*(

*s*

_{i}(

*t*)|

**,**

*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 log-pseudolikelihood function as the following:

*s*

_{i}(

*t*) is binary. The likelihood function may be simplified further by setting

*C*

_{i}(

*t*) =

*β*$\u2211m=1k$

*J*

_{i,m}

*s*

_{m}(

*t*), resulting in the following:

*J*

_{i,j}by computing the partial derivative of the log-pseudolikelihood as the following:

*γ*$\u2202\u2113\u2202Ji,j$|

_{n}. 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*H*(*s*_{i}). - (3)
If

*H*(*s*_{i}) ≤ 0 or*rand*(0, 1) ≤ exp$Hsi\beta simulated$, 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**χ**). - (6)
Do this for all

*β*_{simulated}.

*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* = $\u2211i=1N$ <*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\beta $ (<*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 self-organized 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
|*J*_{i,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 = $4534$, 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.

## AUTHOR CONTRIBUTIONS

Igor Fortel: Conceptualization; Formal analysis; Investigation; Methodology; Software; Validation; Visualization; Writing – original draft; Writing – review & editing. Mitchell Butler: Conceptualization; Methodology; Writing – review & editing. Laura E. Korthauer: Data curation; Investigation; Methodology. Liang Zhan: Conceptualization; Formal analysis; Funding acquisition; Methodology. Olusola Ajilore: Conceptualization; Investigation. Anastasios Sidiropoulos: Conceptualization; Investigation. Yichao Wu: Conceptualization; Investigation. Ira Driscoll: Data curation; Investigation. Dan Schonfeld: Conceptualization; Methodology; Supervision; Validation. Alex Leow: Conceptualization; Formal analysis; Funding acquisition; Methodology; Supervision; Validation; Writing – review & editing.

## FUNDING INFORMATION

Alex Leow and Liang Zhan, National Institutes of Health (https://dx.doi.org/10.13039/100000002), Award ID: R01AG071243. Alex Leow, Yichao Wu, and Liang Zhan, National Institutes of Health (https://dx.doi.org/10.13039/100000002), Award ID: R01MH125928. Liang Zhan, National Institutes of Health (https://dx.doi.org/10.13039/100000002), Award ID: U01AG068057. Liang Zhan, National Science Foundation (https://dx.doi.org/10.13039/100000001), Award ID: IIS2045848. Liang Zhan, National Science Foundation (https://dx.doi.org/10.13039/100000001), Award ID: IIS1837956.

## TECHNICAL TERMS

- 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.

- 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).

- Hamiltonian:
Equation representing the total energy of the system.

- Boltzmann distribution:
Quantifies the probability that a system will be in a specific state based on energy (of that state) and temperature (system).

- Critical temperature:
Parameter describing the tolerance of a system to increasing perturbations (commonly used in physics to identify loss of magnetic properties).

## REFERENCES

*A modern course in statistical physics*, 2nd edition [Book review]

## Author notes

Competing Interests: The authors have declared that no competing interests exist.

Handling Editor: Olaf Sporns