Changes in brain balance of excitation and inhibition by LSD increase overall brain synchrony

In this study, Kuramoto Order Parameter (OP) was employed to assess the alignment of BOLD signal phases across the entire brain under LSD and placebo conditions, with the aim of determining the impact of LSD on overall brain synchrony. Results from paired-sample t-tests indicated that the OP during the LSD condition (0.2708 ± 0.0540) was significantly higher than that observed during the placebo condition (0.2367 ± 0.0305) as shown in Fig 2, with a p-value of 0.0101 (threshold set at 0.05), and the effect size (Cohen's d) was 0.7671. This finding suggests that LSD significantly enhances global brain synchrony, aligning with previous research [4].

Additionally, the standard deviation of the OP, STD(OP), was analyzed to evaluate the dynamic stability of the system, providing insights into the complexity of brain networks under different conditions. The results revealed that the STD(OP) in the LSD condition (0.1210 ± 0.0216) was significantly greater than that in the placebo condition (0.1079 ± 0.0106), with a p-value of 0.0140, and the effect size (Cohen's d) was 0.7246. This indicates that under LSD, the brain network may exhibit transitions between multiple metastable states, reflecting more complex cognitive processes and richer neural activity patterns.

To further elucidate the precise mechanisms underlying the increased synchrony of brain activity induced by LSD, LEiDA was employed. This method characterizes the phase synchronization of BOLD signals across different brain regions at each repetition time (TR). By applying k-means clustering to these synchronization patterns, we derived a repertoire of recurrent phase-locking states (PL states) across all subjects under both LSD and placebo conditions. The results for k = 6, where k represents the pre-defined number of clusters, are illustrated in Fig 3. These PL states provide a compact representation of how the brain dynamically organizes into distinct large-scale functional configurations during rest. We then quantified key dynamic properties of these states—including their probability of occurrence and dwell time—and examined the transitions between them. This approach enables a granular investigation into how LSD alters the brain's dynamic landscape, thereby revealing the specific mechanisms behind the observed increase in synchrony. Fig 3A displays bar plots of each cluster centroid vector (i.e., PL state), where each element of the vector represents the projection of the BOLD phase of each brain region onto the leading eigenvector, indicating the contribution of that brain region to the respective PL state. Notably, PL state 3 exhibits predominantly negative values, corresponding to a global mode in which all BOLD signals change in a unified direction (i.e., all signals project towards the same direction on the leading eigenvector), consistent with findings from previous studies utilizing the LEiDA approach [20,21,27]. In addition, several non-global modes were also identified, and the transitions between these modes influenced the degree of phase alignment of the BOLD signals. The rendering of each PL state based on its values projected onto the cortical surface is depicted in Fig 3B. The matrix representation of the PL states, shown in Fig 3C, is obtained by calculating the outer product of V_c × V_c^T.

To further explore the mechanism of LSD's effects, we examined the spatial correspondence between the six extracted PL states and the seven canonical resting-state networks (RSNs) defined by Yeo et al. [28]. Specifically, we computed Pearson correlation coefficients between the spatial pattern of each PL state (represented by the cluster centroid vector V_c, where negative elements were set to zero) and the spatial distribution of each RSN (represented by a 210-element vector encoding the contribution of each Brainnetome area to the network). This analysis quantifies the degree of spatial overlap between the functional architecture captured by each PL state and the established RSNs, as illustrated in Fig 4.

The spatial distribution of PL state 1 showed significant positive correlations with the spatial distributions of the Default Mode Network (DMN) and the Frontoparietal Network (FPN), both of which are transmodal association networks. Conversely, PL state 1 exhibited significant negative spatial correlations with the Somatomotor (SOM) and Ventral Attention (VAT) networks, which represent unimodal sensorimotor regions. Similarly, the spatial pattern of PL state 2 demonstrated significant negative correlations with the DMN and FPN networks—networks positioned at one extreme of the principal gradient—while showing positive spatial correlations with the SOM and VAT networks located at the opposite extreme. PL state 3 displayed low correlation coefficients (all < 0.12) with all RSNs, and none reached statistical significance, consistent with its globally synchronized yet functionally non-modular pattern as indicated by the nearly uniform negative values in its centroid vector. PL state 4 was significantly positively correlated with the FPN and negatively correlated with the SOM and Visual (VIS) networks in spatial distribution. Finally, both PL state 5 and PL state 6 showed significant positive spatial correlations exclusively with the DMN. All correlation coefficients and corresponding FDR corrected p-values are provided in. S1 Table

The paired-sample t-test, after applying FDR correction for multiple comparisons, revealed that the probability of occurrence in PL state 3 under the LSD condition was significantly higher than that under the placebo condition (p = 0.0168, Cohen's d = 0.8661), as shown in the first row of Fig 5A. No significant group differences were found for the other PL states (see S2 Table for detailed p-values and effect sizes).

In terms of dwell time, after FDR correction, only PL state 3 showed a significant between-group difference (p = 0.0264, Cohen's d = 0.8386). Additionally, although not surviving FDR correction, PL state 6 exhibited a trend towards shorter dwell time under the LSD condition compared to the placebo condition (uncorrected p = 0.0044, Cohen's d = -0.6657), as shown in the second row of Fig 5A. No significant group differences were found for the other PL states (see S3 Table for FDR corrected p-values and effect sizes).

Considering the increased global brain activity synchronization mentioned earlier, it can be speculated that the increase in the probability of occurrence of PL state 3, which represents a global connectivity pattern, under the LSD condition may contribute to the overall higher synchronization of brain activity. After applying FDR correction, the correlation between the probability of occurrence of the different PL states and OP revealed a highly significant positive correlation between PL state 3 and OP (r = 0.8632, p < 0.0001), which supports this hypothesis. The probability of occurrence for both PL state 2 and PL state 6 was significantly negatively correlated with OP, with correlation coefficients of r = -0.4997 (p = 0.0147) and r = -0.4796 (p = 0.0147), respectively. Correlation results for the remaining PL states can be found in. Furthermore, significant correlations were found between the probability of occurrence of PL state 3, and PL state 6 and STD(OP). These changes in the probability of occurrence of the three PL states also contribute to the increased dynamic flexibility of brain activity. S4 Table

To better understand the dynamic trajectories of the PL states, Fig 5B presents the probability matrices for transitions between a given PL state (rows) and other PL states (columns) under LSD (upper panel) and placebo conditions (lower panel). The results show that, compared to the placebo condition, the probabilities of transitioning from PL state 3 to other PL states were all reduced under LSD. Notably, a nominal significance (i.e., uncorrected p < 0.05) was observed for the difference in the probability of transitioning from PL state 3 to itself (uncorrected p = 0.0338, Cohen's d = 0.7273), and for the transition to PL state 5 (uncorrected p = 0.0066, Cohen's d = -0.6535). In contrast, the probabilities of transitioning from other PL states to PL state 3 were all increased, with the transition from PL state 4 to PL state 3 also showing a nominally significant difference (uncorrected p = 0.0376, Cohen's d = 0.5201). However, it is important to note that after applying False Discovery Rate (FDR) correction for multiple comparisons across all state transitions, none of these differences remained statistically significant (all pFDR > 0.05). Despite this, the consistent direction of effects across related transitions and the moderate effect sizes suggest that these patterns may hold biological relevance and warrant further investigation in larger studies. To enhance clarity, Fig 5C displays this result with arrows in different colors. Results of the other paired-sample t-tests can be found in S5 Table.

To determine whether the results were influenced by the value of k, we identified the PL states with the most significant group differences in probability of occurrence and dwell time between LSD and placebo conditions across various partitioning models. For partitioning models with k values of 2, 3, and 4, no significant group differences were found for any PL state in terms of probability of occurrence or dwell time. However, for partitioning models with k values ranging from 5 to 10, significant group differences were observed for certain PL states (see S6 Table for detailed p-values). The PL states with the most significant differences are shown in Fig 6. The bar charts in S1A Fig clearly show that for models with k values from 5 to 10, the cluster center vectors corresponding to significant PL states at each k value are predominantly negative. These PL states also exhibit very similar cortical distribution patterns, all representing global connectivity, as shown in S1B Fig. Additionally, Pearson correlations were calculated between these states, revealing that, with the exception of the PL state with the most significant differences at k = 5, which showed the lowest correlation (r = 0.78) with the most significant PL states at other k values, all other k values showed high correlations (r ≥ 0.95) between the most significantly differing PL states, as shown in S2 Fig. This likely occurs because at a lower k-value (k = 5), the clustering algorithm merges slightly distinct "global states" that are separated at finer partitions (k > 5), leading to a centroid that differs subtly from the "purer" global state identified for k ≥ 6. Beyond this technical observation, the stable emergence of this specific "global state" across multiple partition scales (k-values) suggests it may represent a fundamental or robust dynamical attractor towards which the brain is driven under LSD.

Under the optimal parameters, the Pearson correlation coefficient (r) between the empirical and simulated FC for the LSD condition in the validation set was 0.34827, with a p-value less than 0.0001. For the placebo condition, the corresponding r value was 0.30049, with a p-value also less than 0.0001. The absolute difference (d) between the empirical and simulated FC was 0.0193 for the LSD condition and 0.0442 for the placebo condition. Additionally, the KS distance between the empirical and simulated FCD was 0.0482 for the LSD condition and 0.0719 for the placebo condition. These results suggest that the pFIC model generated realistic FC and FCD metrics.

Intuitively, regions with a negative E/I ratio difference in the cortical areas were primarily concentrated in SOM, as shown in Fig 6A. After FDR correction for multiple comparisons, the Pearson correlation results between the E/I ratio difference and the 7 RSNs revealed a significant negative correlation with SOM (r = -0.6286, p < 0.0001), and significant positive correlations with DMN (r = 0.2410, p = 0.0009), FPN (r = 0.2012, p = 0.0079), and LIM (r = 0.1770, p = 0.0178). The remaining correlation results are presented in S7 Table. As shown in Fig 6B, the distribution of the E/I ratio difference data for the LSD condition was more compact compared to the placebo condition, with closely aligned peak data densities for both conditions. When examining the E/I ratio differences for each of the 7 RSNs separately, the data distribution in the LSD condition was consistently more compact across all RSNs. Specifically, as shown in Fig 6C, the peak data densities for VIS, SOM, DAT, and VAT in the LSD condition were located to the left (i.e., smaller values) of those in the placebo condition, whereas the peak data densities for LIM, FPN, and DMN were located to the right (i.e., larger values) of the placebo condition distributions. The Wilcoxon signed-rank test revealed that the E/I ratio differences for all brain networks were significantly greater than zero for VIS, SOM, DAT, and VAT, or smaller than zero for LIM, FPN, and DMN. The p-value for LIM was 0.001, and all others were less than 0.001.

A key mechanism of LSD's effects is its ability to stimulate the 5-hydroxytryptamine-2A (5-HT_2A) receptors, which in turn activates glutamatergic transmission in the prefrontal cortex, ultimately influencing downstream subcortical regions and altering sensory and cognitive processing [29]. To further investigate this relationship, we calculated the Pearson correlations between the E/I ratio difference and the densities of 5-HT_2A and glutamate receptors. Results revealed significant positive correlations between the E/I ratio difference and both 5-HT_2A receptor density (r = 0.4593, p < 0.0001) and glutamate receptor density (r = 0.1967, p = 0.0042), as shown in Fig 6D.

The experimental protocol received approval from the West London Research Ethics Committee of the UK National Health Service. The experiments adhered to the revised Declaration of Helsinki (2000), the International Council for Harmonisation Good Clinical Practice guidelines, and the National Health Service Research Governance Framework. Twenty participants were recruited via word of mouth and provided written informed consent after a thorough study briefing and screening for physical and mental health.

The present analysis utilized neuroimaging data from a publicly available dataset on OpenNeuro, with detailed information provided in the original publication [4]. The screening included an electrocardiogram (ECG), routine blood tests, and a urine test to check for recent drug use and pregnancy. Participants also underwent a psychiatric interview to disclose their complete drug use history. Key exclusion criteria included being under 21 years of age, personal history of diagnosed psychiatric illness, immediate family history of psychotic disorders, lack of prior experience with classic psychedelic drugs (such as LSD, mescaline, psilocybin, or DMT), any psychedelic drug use within six weeks prior to the first scanning day, pregnancy, problematic alcohol consumption (more than 40 units per week), or any medically significant conditions that would render the volunteer unsuitable for the study.

Following the initial exclusion criteria, 20 healthy participants underwent resting-state fMRI Blood Oxygen Level Dependent (BOLD) data acquisition under both LSD and placebo conditions, with a minimum interval of 14 days between the two conditions. The order of administration of LSD or placebo was counterbalanced across subjects to control for potential order effects. Participants were blinded to the order of administration, while the researchers were informed of it. On each scanning day, each participant received an intravenous injection of 75 µg of LSD dissolved in 10 mL of saline for the LSD condition, or 10 mL of saline for the placebo condition, 70 minutes prior to the MRI scan. The acquisition protocol is elaborated upon in greater detail in the original publication [4]. The fMRI BOLD data were recorded using a specific gradient echo planar imaging sequence and underwent detailed preprocessing steps as described in the original publication [31].

BRANT version 3.35, a MATLAB toolbox specifically designed for the batch preprocessing of fMRI data, was utilized in this study [63]. The preprocessing procedures were comprehensively detailed by [64]. Following a series of preprocessing steps, including band-pass filtering (0.01–0.08 Hz) and motion correction, BOLD time series were extracted for each brain region. Specifically, for each participant at each condition, the fMRI data were averaged within the 210 regions of interest (ROIs) defined by the Brainnetome Atlas, resulting in a 210 (ROIs) × 216 (time points) matrix.

These matrices were then used to compute functional connectivity (FC) matrices (210 × 210) by calculating the pairwise correlations between the time courses of all ROIs. The group-level FC matrix was computed by averaging the participant-level FC matrices separately for the training (n = 7) and validation (n = 8) sets, with 15 participants randomly assigned to each set. The calculation of the Functional Connectivity Dynamics (FCD) matrix was carried out as follows. A sliding window of 60 seconds (equivalent to 30 TRs, given a TR of 2 seconds) was employed, consistent with recommendations from prior research [26]. This window was shifted from the first frame to the 187th frame of the BOLD time series, resulting in a total of 187 sliding windows. Each FC matrix corresponding to a sliding window was then vectorized by focusing solely on the upper triangular entries, and the vectorized FCs were subsequently correlated with one another, generating a 187 × 187 FCD matrix.

The diffusion tensor imaging (DTI) data from 1065 participants in the Human Connectome Project (HCP) S1200 release were utilized in this study. Detailed information regarding data collection is available in previous literature [65]. The participants were randomly divided into a training set (N = 532) and a validation set (N = 533).

To construct structural connectivity (SC) matrices for each participant, whole-brain tractography was performed using the fiber orientation distribution (iFOD2) algorithm [66] provided by MRtrix3 [67]. Each SC matrix entry represented the number of streamlines between two regions of interest (ROIs). A thresholding procedure was applied to eliminate false positives, setting entries to zero if fewer than 50% of participants exhibited a non-zero value. The number of streamlines was averaged across participants with non-zero values and log-transformed [67], while the main diagonal entries were set to zero. Group-level SC matrices were generated by averaging the individual-level SCs for both the training and validation sets, with the maximum value normalized to 0.02 following the previous study [26].

To evaluate the synchrony and metastability of brain network dynamics, we utilized the Kuramoto order parameter (OP), a metric commonly employed to quantify the collective dynamics of coupled oscillators, especially in systems consisting of multiple interacting units [68 –71]. This parameter provides data-driven framework for measuring the extent of phase coherence or the degree of synchrony among the oscillators, here, the brain regions, wherein the phase represents a large-scale hemodynamic rhythm rather than the instantaneous phase of neuronal spiking. Specifically, the Kuramoto order parameter at time point t is defined as:

where N is the number of oscillators, corresponding to the 210 brain regions in this study, and θ(n, t) represents the phase of the region n at time point t.

In the study, the BOLD signal from brain region n was treated as a continuous oscillatory signal x(n, t). The Hilbert transform takes this real-valued signal x(n, t) and creates a 90-degree phase-shifted version x^(n, t). Then the original signal x(n, t) and its Hilbert transform x^(n, t) were combined to form the analytic signal z(n, t)=x(n, t)+ix^(n, t), which can be expressed as A(n, t)eiθ(n, t), where A(n, t) is the instantaneous amplitude and θ(n, t) is the instantaneous phase [20,72]. The phase θ(t) of each brain region was then used to compute the Kuramoto order parameter at each time point.

The mean OP across time provides a measure of the system's overall synchrony. A value close to 0 indicates that the system remains largely desynchronized, with the brain regions operating independently and showing minimal phase coherence. A value between 0 and 1 suggests partial synchrony, where certain periods exhibit some degree of coordination among brain regions. A value approaching 1 indicates full synchrony, with most brain regions oscillating in phase, reflecting coordinated network activity. The standard deviation of the OP, STD(OP), was used to quantify the metastability of the system [73,74]. The STD(OP) reflects the temporal variability of synchrony: when OP remains relatively stable over time, it suggests that the system is in a stable equilibrium, whether synchronized or desynchronized. On the other hand, a higher standard deviation indicates that the system is dynamically shifting between different states of synchrony, characteristic of a metastable system, and implies that the system is exploring a broader range of functional metastable states. Finally, a paired-samples t-test was conducted to evaluate the differences in the OP and STD(OP) between the LSD and placebo conditions.

To extract the repertoire of large-scale functional networks, the phase-locking matrix was computed for each subject across 216 time points:

where θ(n,t) and θ(p,t) denote the instantaneous phase of brain region n and p, as shown in Fig 1A. This resulted in a three-dimensional matrix of size N × N × T for each subject, where N = 210 as shown in Fig 1B. Then, the leading eigenvector V₁(t) of the phase-locking matrix for each subject at each time point t was computed, which is an N × 1 vector serving as a low-dimensional representation of the BOLD phase-locking patterns, encapsulating the primary orientation of BOLD phases across all brain areas. Each element of V₁(t) reflects the projection of the BOLD phase of a given brain area into this leading eigenvector, indicating the degree to which each area aligns with the overall network mode. When all elements of V₁(t) share the same sign (either all positive or all negative), it indicates a state of full synchrony, where all brain regions oscillate in phase with respect to V₁(t). Conversely, differing signs among the elements suggest the presence of distinct clusters of brain areas, each exhibiting different phase relationships with respect to the leading eigenvector V₁(t).

To identify recurrent phase-locking states, k-means clustering was applied to the collection of leading eigenvectors V₁(t) across all subjects, time points, and conditions. A range of cluster numbers k from 2 to 10 was explored, allowing for the capture of varying levels of functional network granularity. In the current analysis, we identified the cluster solution k = 6 as the most favorable option by considering four evaluation metrics—Dunn score, distortion, silhouette score, and Davies–Bouldin score—while also aiming to optimize the statistical significance of the differences observed between patient groups. The resulting cluster centroids, denoted as V_c, represent distinct phase-locking states (PL states), capturing recurrent patterns of phase synchrony across the brain, which can be interpreted as metastable brain network configurations. Each PL state can be visualized as a network in cortical space, as shown in Fig 1C, where the magnitude of V₁(t) (the n-th element of the cluster center vector) is used to scale the color of each brain area, highlighting the relative contributions of different regions to the network.

Given that both V and -V occupy the same one-dimensional subspace, we adopted a convention where the majority of elements are negative. This choice is supported by previous study [27], which indicates that the smallest subset of regions, whose BOLD phases diverge from the global mode, corresponds to significant functional brain networks, specifically the canonical resting-state networks.

In this study, Pearson correlation was employed to assess the spatial overlap between the PL states and the seven resting state networks (RSNs) [27]. Initially, the seven RSNs defined by Yeo et al. (2011) [28] were transformed from the 1 mm³ MNI space into the Brainnetome Atlas space. This was achieved by counting the number of 1 mm³ MNI voxels within each Brainnetome area that belong to each of the seven networks, resulting in a vector with 210 elements that reflects the contribution of each Brainnetome area to the corresponding RSN. Subsequently, the Pearson correlation coefficients were computed between each pair of the 6 PL states and the seven 210-element vectors. The resulting p-values were corrected for multiple comparisons using the false discovery rate (FDR) procedure (q < 0.05). During this process, all negative elements in the V_c vector were set to zero, retaining only the positive values [27].

The clustering algorithm assigns each TR to the closest centroid VC, thereby designating a specific PL state for each time point of each subject, facilitating the calculation of the probability of occurrence for each PL state, dwell time, as well as the state transition probability matrices. The probability of occurrence for each PL state was computed by dividing the number of epochs assigned to that PL state by the total number of epochs within each scanning session for each partition model (ranging from k = 2 to k = 10), as shown in Fig 1D for k = 6. Additionally, the dwell time was obtained by calculating the average duration spent continuously in a specific PL state throughout the entire scan.

To assess the differences in probabilities of occurrence between LSD and placebo conditions, a permutation-based paired t-test was employed. This non-parametric test generates the null distribution by randomly permuting the group labels (LSD and placebo) for each subject. For each of 1000 permutations, a paired t-test was applied to compare the populations, resulting in corresponding p-values.

To further investigate the dynamics between different PL states, switching probability matrices were constructed to represent the mean probabilities of being in a given PL state (rows) and transitioning to other PL states (columns) under both LSD and placebo conditions. The switching probabilities were determined for each subject under both conditions, and inter-group differences were statistically evaluated using a permutation-based paired t-test, with 1000 permutations.

The large-scale spontaneous brain activity is believed to arise from the interactions of intrinsic dynamics within local circuits over long-range anatomical connections [54,75,76]. To understand the mechanisms by which LSD affects brain networks from a neurobiological perspective, we simulated the BOLD signals and estimated E/I ratio distribution under both LSD and placebo conditions. Previous studies have demonstrated that parameterizing local circuit properties using anatomical and functional gradients can yield more realistic static and dynamic resting-state FC [26,54]. This heterogeneity of local circuits can be informed by in vivo structural and functional neuroimaging measurements. The Parametric Feedback Inhibition Control (pFIC) Model incorporates estimates of cortical myelin derived from T1-weighted/T2-weighted (T1w/T2w) MRI and the principal resting-state FC gradient, both of which have been shown to reflect anatomical and functional hierarchies [77,78].

The derivation of the FIC model has been extensively detailed in a prior study [79]. In this section, we aim to offer some insights into the FIC model. The neuronal dynamics of the j-th cortical region are governed by the nonlinear differential equations presented below:

where I_j, r_j, and S_j represent the synaptic currents, firing rates, and synaptic gating variables of the excitatory (E) or inhibitory (I) neuronal populations in region j, respectively [22]. The first component of I_j, denoted as WI0, represents the overall effective external input, where W_E = 1, W_I = 0.7, and I₀ = 0.382nA. The second component accounts for the intra-regional E-to-E currents, with ω_EE and ω_EI representing the weights of recurrent excitation and excitatory-to-inhibitory connections, respectively. The third component captures inter-regional inputs, where G is a global coupling factor which are unknown parameters to be estimated by fitting to empirical fMRI data and uniformly scales all E-to-E connections across regions, as determined by the structural connectivity (SC) matrix, represented in C_jk (the connection between region j and k). The fourth component describes the local negative feedback from the inhibitory population, governed by the strength of the inhibitory-to-excitatory connection ω_IE (calculated analytically to maintain a typical noisy spontaneous activity with low firing rate of 3Hz observed in electrophysiology experiments to emulate resting-state conditions) [79 –81] and the strength of the inhibitory-to-inhibitory recurrent connection ω_II (set to 1).

ϕ denotes the neuronal response function [82] that converts the input current I_j into the firing rate r_j. The parameter a serves as the gain factor that determines the slope of ϕ, with a_E = 310n/C and a_I = 615n/C. b represents the threshold current, above which the firing rate r_j increases linearly with the input current, where b_E = 125Hz and b_I = 177Hz. d is a constant that determines the curvature of ϕ around the threshold current, with d_E = 0.16s and d_I = 0.087s.

The synaptic gating variable S_j^(E) of the excitatory pool is regulated by N-methyl-D-aspartate (NMDA) receptors with a decay time constant of τ_E = 0.1s and γ = 0.641. In contrast, the synaptic gating variable S_j^(I) of the inhibitory pool primarily depends on γ-aminobutyric acid (GABA), with a decay time constant of τ_I = 0.01s. Additionally, v_j represents an uncorrelated standard Gaussian noise with an amplitude of σ = 0.01nA. All parameter settings were based on previous studies [83]. In this study, the regional E/I ratio is defined as the ratio of the temporal average of the excitatory synaptic gating variable S_j^(E) to the inhibitory synaptic gating variable S_j^(I) of this region. Using the simulated excitatory synaptic gating variables S_j^(E) as inputs for the Balloon-Windkessel hemodynamic model [79,84,85] to simulate fMRI BOLD signals allows for the subsequent generation of simulated static FC and FCD.

As mentioned earlier, spatial heterogeneity in the excitatory-to-excitatory recurrent connection strength ω_EE, excitatory-to-inhibitory connection strength ω_EI, and regional noise amplitude σ must be taken into account in order to simulate more realistic brain activity. The use of Brainnetome Atlas with 210 cortical regions results in 210 × 3 + 1 (where the additional 1 represents the global coupling constant G) unknown parameters, posing a significant computational challenge. In this study, the pFIC model integrates estimates of cortical myelin derived from T1w/T2w MRI and the principal resting-state functional connectivity (FC) gradient by parameterizing ω_EE, ω_EI, and σ as a linear combination of the myelin map and the FC gradient, reducing the number of free parameters to 3 × 3 + 1 while still simulating more realistic brain activity, as follows:

Following previous studies [26,54], the objective of the fitting is to minimize the function (1−r)+d+KS, where r is the Pearson correlation coefficient between the empirical and simulated FC matrices, d represents the absolute difference between the average empirical FC matrices and simulated FC matrices, the inclusion of which accounts for scale differences between the empirical and simulated FC matrices. KS denotes the Kolmogorov-Smirnov (KS) distance. Specifically, since there is no temporal correspondence between the empirical and simulated FCD, it would be inappropriate to directly compute the Euclidean distance between them [26]. Instead, the difference between the empirical and simulated FCD matrices is quantified by the maximum distance between their cumulative distribution functions (CDFs), calculated using the Kolmogorov-Smirnov (KS) distance.

To minimize the objective function, Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [54] was employed. During the training phase, CMA-ES was initialized randomly five times, with each initialization undergoing 100 iterations, resulting in a total of 500 candidate parameter sets. For each set of parameters, the simulation was repeated 20 times on the validation set, and the mean FC and FCD were computed from these 20 simulations, which were then used to optimize the objective function, ultimately yielding the best parameter set. The E/I ratio differences were calculated as: