- 1ETIS, ENSEA, CNRS, UMR8051, CY Cergy-Paris University, Cergy, France
- 2Laboratoire de Physique Théorique et Modélisation, UMR 8089, CNRS, Cergy-Pontoise, CY Cergy Paris Université, Cergy, France
- 3Simulation and Data Lab Neuroscience, Institute for Advanced Simulation, Jülich Supercomputing Centre (JSC), JARA, Forschungszentrum Jülich GmbH, Jülich, Germany
- 4Department of Biomedical Engineering, University of Illinois at Chicago, Chicago, IL, United States
- 5Department of Psychology, University of Wisconsin-Milwaukee, Milwaukee, WI, United States
- 6Department of Electrical and Computer Engineering, University of Pittsburgh, Pittsburgh, PA, United States
- 7Department of Psychiatry, University of Illinois at Chicago, Chicago, IL, United States
The human brain, composed of billions of neurons and synaptic connections, is an intricate network coordinating a sophisticated balance of excitatory and inhibitory activities between brain regions. The dynamical balance between excitation and inhibition is vital for adjusting neural input/output relationships in cortical networks and regulating the dynamic range of their responses to stimuli. To infer this balance using connectomics, we recently introduced a computational framework based on the Ising model, which was first developed to explain phase transitions in ferromagnets, and proposed a novel hybrid resting-state structural connectome (rsSC). Here, we show that a generative model based on the Kuramoto phase oscillator can be used to simulate static and dynamic functional connectomes (FC) with rsSC as the coupling weight coefficients, such that the simulated FC aligns well with the observed FC when compared with that simulated traditional structural connectome.
1 Introduction
The human brain is a complex neural network that self-organizes into different emergent states, which are crucial for its functions. Such states include spatiotemporal patterns of neural synchronization associated with cognitive processes (Bansal et al., 2019). Brain regions can be modeled as dynamically interacting nodes in a functional network on a 3D space (functional brain networks), which was coupled in a complex manner driven by the structure of these networks. Over the past years, interdisciplinary approaches, using concepts from non-linear dynamics, physics, biology, and medicine to name a few, allowed us to understand in more depth how the human brain functions and how certain brain disorders and their underlying mechanisms can be further studied using mathematical models. It is feasible to ameliorate even more the predictive performance of such models, since a vast amount of neuroimaging data e.g., electroencephalography (EEG), magnetoencephalography (MEG), and functional magnetic resonance blood-oxygen-level-dependent (BOLD) functional magnetic resonance imaging (fMRI) became available in the last two to three decades. Such data may provide information not only for healthy or pathological brain activity but can also be used to fingerprint functional connectomes by identifying individuals using brain connectivity patterns (Finn et al., 2015).
Together with extensive experimental study, mathematical/computational modeling of the whole brain dynamics has been an active research topic for years (Deco et al., 2008; Sanz-Leon et al., 2015; Jirsa et al., 2017; Murray et al., 2018; Young, 2020). In such a setting, one can model populations of neurons as nodes in a graph structure. Then, one can obtain information about relative connection weights (coupling strength) and communication lag (delay) between different nodes by diffusion-weighted magnetic resonance imaging (dwMRI) techniques (Ghosh et al., 2008; Hagmann et al., 2010; Deco et al., 2011). This is termed as the structural connectivity (SC) of the network, and it is, in general, subject-dependent with a certain degree of variability (gender/age/healthy vs. diseased). Furthermore, statistical analysis of BOLD time series inferred from fMRI can provide the functional relationships between different brain regions. It is usually calculated as the Pearson correlation coefficient of the activity between regions and results in the empirical functional connectivity (FC) matrix per brain recording and subject (Sporns et al., 2005; Horn et al., 2014).
By working in silico, one can seek for model parameters that are able to produce simulated time series and global dynamics that fairly resemble the empirical ones. One way in achieving that is to tune selected parameters which optimize the similarity between empirical FC and simulated FC (Cabral et al., 2011; Deco and Jirsa, 2012). Hence, these parameters can serve as dynamical biomarkers and predictors of different brain states and behavioral modes [see Popovych et al. (2019) for a recent review]. Along this direction, the virtual epileptic patient has been recently proposed, where medical-treatment approaches using personalized mathematical models for epileptic patients have been illustrated (Jirsa et al., 2017). Furthermore, the choice of the brain atlases, (i.e., the mapping of the different regions of interest (ROIs) based on functional or anatomical criteria using different parcellations) can affect the quality of model performance and its level of agreement with the empirical data [see Popovych et al. (2021) and references therein for more details].
In recent years, substantial research efforts have been directed toward understanding the brain (large-scale activity) using resting state fMRI (rs-fMRI), employing sophisticated mathematical and statistical tools to investigate the FC from rs-fMRI data (Biswal et al., 1997). So far, the mainstream approach is to consider SC to be static and the FC to be dynamic. However, this is not necessarily the case as the white matter tracts can be in use or engaged when the brain is performing certain tasks but inactive or disengaged during other tasks and hence not static. An altered and more sophisticated “functional connectivity-informed structural connectivity” has been introduced by Ajilore et al. (2013), employing information from fMRI to infer the underlying pattern of white matter engagement specific to the state of the brain. The resulting connectome, the so-called resting-state informed structural connectome (rsSC), encodes the structural network that underlies and facilitates the observed rs-fMRI correlation connectome which is able to detect altered rsSC community structure in diseased subjects relative to controls. In the original set up, there is no “directionality” inferred, i.e., whether the white matter tract of interest is of “excitatory” vs. “inhibitory” nature.
However, understanding the dynamical balance between excitation and inhibition, a concept termed E-I balance is vital for adjusting neural input/output relationships in cortical networks and regulating the dynamic range of their responses to stimuli (Kinouchi and Copelli, 2006) such that information capacity and transfer are maximized (Shew et al., 2011). This is the central thesis of the criticality hypothesis (Beggs and Plenz, 2003; Muñoz, 2018), i.e., that brain activity self-organize into a critical state (Wilting and Priesemann, 2019), a unique configuration likened to a phase transition in physical systems where a dynamical system transitions from order (balanced excitation-inhibition) to disorder (disrupted excitation-inhibition balance) (Sornette, 2004; Cocchi et al., 2017; Hahn et al., 2017; Tagliazucchi, 2017). Indeed, evidence supporting that the brain is operating near criticality has been reported in studies examining neuronal signaling (Beggs and Plenz, 2003; Shew et al., 2009; Hahn et al., 2017) and BOLD fMRI signals (Tagliazucchi et al., 2012; Haimovici et al., 2013; Lombardi et al., 2017; Rabuffo et al., 2021).
To incorporate co-activation (excitatory) or silencing (inhibitory) effects into our hybrid rsSC framework that would allow us to infer the E-I balance of the brain, in the study by Fortel et al. (2019), we then introduced an improved framework based on the Ising model representation of the brain as a dynamical system, wherein self-organized patterns are formed through the spontaneous fluctuations of random spins. This Ising spin-glass model has been previously used to successfully characterize complex microscale dynamics (Tkacik et al., 2015; Kadirvelu et al., 2017) and macroscale interactions (Schneidman et al., 2006; Marinazzo et al., 2014; Ezaki et al., 2017; Nghiem et al., 2018; Niu et al., 2019; Nuzzi et al., 2020) of the human brain and to accurately represent spatiotemporal co-activations in neuronal spike trains (Schneidman et al., 2006; Shlens et al., 2006; Roudi et al., 2009) and patterns of BOLD activity (Watanabe et al., 2013; Ashourvan et al., 2017; Cocco et al., 2017; Ezaki et al., 2020).
In this study, we use The Virtual Brain (TVB, Sanz-Leon et al., 2015), a whole-brain simulation platform part of the EBRAINS infrastructure (https://ebrains.eu/), to investigate the potential benefits in employing rsSC instead of the traditional SC for simulating whole-brain dynamical activity. For example, one major limitation of using traditional SC with certain dynamical models, such as the Kuramoto phase oscillators (Kuramoto, 2003) and the generic limit-cycle oscillators (Kuznetsov, 1998), to model the mean neural activity of each node using traditional SC connectomes, is that the resulting signals from different ROIs do not adequately produce negative correlations obtained in the empirical ones, even in the presence of delays in the system (Popovych et al., 2021). Here, we show that using rsSC such as dynamical systems succeed to produce simulated signals with both positive and negative correlations which sufficiently follow the trends of the empirical ones.
2 Methods and materials
2.1 Empirical data and signed resting state structural connectome
The structural and functional connectivity (resting state) for 38 cognitively normal APOE ε4 allele carriers with the detailed information on the imaging and processing steps can be found in the study by Korthauer et al. (2018). The algorithm to obtain rsSC can be found in the study by Fortel et al. (2019, 2020); Tang et al. (2021); Fortel et al. (2022), while in the Supplementary material, we provide a concise description with details on its implementation.
2.2 Models and simulated data
To produce simulated fMRI time series in the given connectomes, we employ the Kuramoto phase oscillator model (Kuramoto, 2003; Lee and Frangou, 2017; Popovych et al., 2021):
where θi are the phases, N is the number of oscillators, fj are the natural frequencies (Hz), cij and τij (ms) represent the individual coupling weight and propagation delay in the coupling, respectively, from oscillator j to oscillator i, while K is the global coupling parameter. The time t in the model and delay in coupling term are measured in ms.
For each individual subject, we produced a “personalized” model Equation (1) to simulate dynamics of the network and calculate time series. To this end, two cases of connectivity matrices were compared: (i) in the first one, the cij values are defined by simply counting the number of streamlines connecting regions i and j normalized to 1 and with zero diagonal (i.e., define cij as a normalized version of the empirical tractography-derived SC or eSC), leading to only excitatory interactions between ROIs; and (ii) in the second one, the cij values are assigned by the corresponding entries of the hybrid rsSC connectomes, leading to both excitatory and inhibitory interactions between ROIs. The delays τij were calculated as τij = Lij/V, where Lij (mm) is the average tract (path) length of the streamlines connecting regions i and j, and V (m/s) is an average velocity of signal propagation. In this particular dataset, the exact path lengths are not available; hence we used the euclidean distance between nodes in the Desikan atlas (Desikan et al., 2006) as proxies. The euclidean distance has been used in the literature in the construction of structural networks (Ercsey-Ravasz et al., 2013) and found to closely follow the trends obtained by anatomical tract-tracing studies. Furthermore in the study by Deco et al. (2021), the authors showed that such networks also strongly correlate with MRI tractography-based networks. The matrix L = Lij can thus be used to calculate the delays τij in the coupling, which can be expressed as τij = τ·Lij/〈Lij〉, where τ = 〈Lij〉/V is the global (or mean) delay. In Equation (1), the self-connections were excluded by setting the diagonal elements in the matrices eSC/rsSC and L to zero (i.e., cii = Lii = 0, respectively).
In Figure 1A, we show the empirical SC matrix (weights of the node-to- node connections) for a Non-Carrier subject from the dataset with 80 nodes (ROIs). Figure 1B shows the corresponding SC and subject rsSC matrix calculated as described earlier. Notably, the hybrid rsSC contains negative entry values as opposed to SC one that is restricted to having only positive values. Figure 1C) shows the tract length L matrix (in mm) that we used for the simulation of all subjects in the absence of the actual measured ones from a neuroimaging prepossessing pipeline. The phases in our model (Equation 1) were initialized randomly. We set the intrinsic frequencies to be uniformly distributed with mean = 60 Hz and SD = 1 Hz, corresponding to oscillations within the gamma frequency range [see Cabral et al. (2011); Mess et al. (2014); Váča et al. (2015); Lee and Frangou (2017) for more details and motivation], as gamma local field potential (LFP) power is coupled to the BOLD fMRI signal and is considered representative of the overall neuronal activity (Niessing et al., 2005; Nir et al., 2007; Miller et al., 2009; Schölvinck et al., 2010)].
Figure 1. Empirical connectivity matrices example (non-carrier subject). (A) Weights of SC matrix. (B) Hybrid rsSC matrix. (C) Tract length (mm) matrix L based on the euclidean distance of the nodes on the Desikan atlas (same for all simulations and subjects).
For our simulations, we used a TVB tailored version for the Kuramoto model and we made adjustments for efficient parallelization on CPUs using MPI on the supercomputer JUSUF located at the Jülich Supercomputing Center. Our model generates time series which correspond initially to electrical activity (fast oscillations) for each node, i.e., we register the observable xi = sin(θi) for each brain region. Then, in order to estimate the simulated BOLD signal, we use the built-in tool of TVB to calculate the hemodynamic response function kernel (i.e., “fMRI activity”) associated with a given neural activity time series, also known as the Balloon-Windkessel model (Friston et al., 2003). Our simulations ran for 500 s. The first 20 s were discarded to remove transient effects, resulting in T = 480 s (8 min), i.e., a time interval identical to the time-length of the empirical fMRI signals. We set the time-step at 0.1 ms, and we integrated the system with an Euler scheme. In this particular study, we did not consider the presence of noise.
3 Results
We numerically simulate BOLD time series varying two model parameters, namely, the global coupling strength K and the delay τ in Equation (1), with respective ranges K ∈ [1, 75] and τ ∈ [1, 33], resulting in a 32 × 32 grid. For each pair of parameters, we begin by producing the matrix of the simulated FC (sFC). The latter is measured by the Pearson Correlation Coefficient (CC) between the simulated BOLD signals xi, i = 1, 2, …, N from different ROIs [also referred to as Static Functional Connectivity in the literature, see Cabral et al. (2017)] as follows
Then, we compare each sFC with the eFC ones using the Pearson Correlation Coefficient; however, this time we calculate it for the two respective matrices (upper triangular parts) as follows:
The optimal match between sFC and eFC in the parameter space is acquired for (K, τ)-values where CCFC becomes maximal [see also Popovych et al. (2021) and references therein for more details and motivation].
In Figure 2, we present the first main result, namely, the superiority of hybrid rsSC over standard SC matrices in generating simulated BOLD time series with models such as Equation (1), which better approximate the empirical BOLD signals (shown here for one example healthy subject). The upper row refers to simulations performed using the standard SC matrix of the respective subject to define the coupling weights in the Kuramoto model. Figure 2A shows the parameter sweep exploration (PSE) for eFC vs. sFC and for the parameters (K, τ) and measured as CCFC = corr(sFC, eFC). The five white circles on the red regions indicate the highest correlations (sizes of larger circles correspond to larger CCFC values). In Figure 2B, we present the eFC calculated from the empirical BOLD signal while in Figure 2C the sFC matrix with the larger ccFC. We can observe that sFC did not capture adequately the negative correlations that are present in eFC (compareing the minimum values in the barplots of Figures 2B, C).
Figure 2. Parameter sweep exploration, correlation & Statistical analysis for eFC vs. sFC example. First row (using the respective subject's standard SC matrix to define the weights in model Equation (1): (A) The colormap depicts the CCFC = corr(sFC, eFC) for the parameters (K, τ). The 5 circles on the red regions indicate the highest correlations found (larger circles' sizes correspond to larger CCFC values). (B) The eFC calculated from the empirical BOLD signal. (C) The sFC matrix with the larger CCFC. Second row (using the respective subject's hybrid rsSC matrix to define the weights in model Equation (1). (D) The respective colormap for CCFC = corr(sFC, eFC). (E) The eFC calculated from the empirical BOLD signal [same as (B)]. (F) The sFC matrix with the larger CCFC. Note the ranges for the two colorbars in (A, D) are kept intact (i.e., no scaling) for visualization purposes. (G, H) Scatterplots between empirical (y−axis) and optimal simulated BOLD correlations (x−axis) aggregated across all entries in the corresponding FC matrices, i.e., (B, C). (G) eFC vs. sFC using rsSC, (H) eFC vs. sFC using standard SC. Both panels refer to same subject (blue lines indicate the linear regression model). (I, J) boxplots for the non-carriers and carriers dataset (38 subjects per type) correlation coefficients between eFC and sFC using SC and rsSC matrices for the simulated time series respectively. For each subject we considered the 5 maximum values (see circles in (A, D)). The difference in the respective mean values of the two datasets is statistically significant measured by the t-test with very small p−value (p ≤ 0.0001) for both non-carriers and carriers sets.
In the second row of Figures 2D–F, we perform similar simulations; however, we now use the hybrid rsSC matrix of the respective subject to define the coupling weights in the Kuramoto model. Notably, the significant improvement in the maximum value of the CCFC≈0.86 compared with the one found when using the standard SC matrix (CCFC≈0.33). Notably, also the better agreement is observed between the two FC matrices (empirical (E) and simulated (F)) and how better the sFC captures both positive and negative correlations (indicated by the range of the respective colorbars). We should stress that we did not opt to use the same range for the two colorbars in Figures 2A, D, as in this way it would be difficult to visually identify the PSE region in (A), depicting the optimal parameter values.
The respective scatterplots and CC values between empirical and optimal simulated FC matrices are presented in the third row of Figure 2 using rsSC (G) and standard SC (H) matrices. Here, we plot the empirical (y-axis) against the optimal simulated BOLD correlations (x-axis) aggregated across all entries in the corresponding FC matrices; thus, a perfect match between the two would place all the points along the line x = y. The higher CCFC(sFC, eFC) value (using hybrid rsSC matrices) is well reflected by a rather clear linear trend in the distribution of the points (Figure 2A). On the other hand, only a relatively weak linear trend is obtained using standard SC matrices (Figure 2B). Both panels refer to the same subject presented in Figure 2, with the lines indicating the corresponding linear fit in each case.
In the forth row of Figure 2, we present a statistical analysis for all subjects per category, i.e., 38 non-carriers (I) and 38 carriers (J). For each subject, we considered the five maximum values of correlation coefficients between eFC and sFC using SC and rsSC matrices for the simulated time series, respectively, (circles in (A) and (D)) and produced boxplots. We then used the t-test to measure the difference in the mean Pearson correlation value of each group (carriers or non-carriers) between the empirical FC vs. simulated FC (5 optimal cases) when using standard SC and rsSC ones.. The difference in the respective mean values of the two datasets is found to be statistically significant with very small p-value (p ≤ 0.0001) for both non-carriers. The respective BOLD time series (empirical, optimal simulated using SC and rsSC respectively) from this example can be found in Supplementary Material (Supplementary Figure S1). The separation of the two groups is indeed biologically important given the already demonstrated differences in E-I dynamics (Fortel et al., 2019, 2020, 2022, 2023). As this current study leveraged one of the several datasets that we previously used to demonstrate sex-by-ε4 hyperexcitation, by showing that the model fits between the two groups are equally optimal we further establish that differences in E-I dynamics are not an artifact secondary to differences in model fit.
Next, we conducted a statistical analysis of negative and positive correlations in empirical (eFC) and simulated (sFC) functional connectomes (Figure 3). We used boxplots to visualize correlation coefficients for eFC and sFC matrices informed by structural connectivity (SC) or resting-state SC (rsSC) for the non-carrier group. We considered sFC matrices produced by the five parameters (K, τ) that maximized eFC and sFC similarity per subject (marked with circles in Figures 2A, D). The leftmost three boxplots show minimum (negative) correlations in actual eFC (light green), which was simulated using rsSC (light blue) and using standard SC (red), respectively, while the rightmost three boxplots show maximum (positive) correlation values.
Figure 3. Statistical Analysis for negative and positive correlations in empirical and simulated FCs. Boxplots of the correlation coefficients for eFC and sFC obtained by using either SC or rsSC as the coupling coefficient matrices (non-carriers group) in the simulations.. For each subject we considered the parameters (K, τ) which correspond to the 5 maximum values that optimize the similarity between eFC and sFC matrices (indicated with circles in Figures 2 A, D)). The 3 leftmost boxplots indicate the negative correlations in the actual eFC (light green), simulated using rsSC (light blue) and using standard structural connectome SC (red). Similarly, the 3 rightmost boxplots compare the maximum (positive) correlation values. Similar results were also obtained for the carriers, see Supplementary Figure S2.
The simulated functional connectomes (sFC), generated using both rsSC (used here for the first time in simulating BOLD time series) and standard SC, sufficiently recovered the positive correlations observed in the empirical BOLD signals. However, the sFC derived from SC (left red boxplots) do not correctly recover the negative correlations present in the empirical data (left light green boxplots). In contrast, the BOLD signals generated with rsSC (left light blue boxplots) exhibit negative correlations much closer to those observed in the empirical data. This consistent trend holds true for both non-carrier and carrier datasets (see Supplementary Figure S2). A more detailed analysis on the role of positive and negative coupling coefficients in rsSC vs. SC connectommes in simulating high-fidelity fMRI correlations with a dynamical system can be found in the discussion and figures of Supplementary material.
Next, to further explore the advantage in using the hybrid rsSC matrices beyond Static Functional Connectivity matrices, we sought out to perform a similar PSE analysis for Dynamic Functional Connectivity, which allows us to capture switching trends in the resting state activity. To this end, we calculate the Phase Coherence Connectivity [see e.g., Cabral et al. (2017); Hancock et al. (2022) and references therein], which does not suffer from time-window length effects such as other similar techniques based on calculating successive FC(t) matrices using a sliding-window [see discussion in Cabral et al. (2017); Hancock et al. (2022)]. Hence, we use BOLD Phase Coherence Connectivity to measure time-resolved dynamic FC matrices (dFC), with size N × N × T, where N refers to the number of ROIs and T = 236 refers to the total number of recording frames. Then, we begin by estimating the phases from the BOLD time series (empirical and simulated) for all ROIs i (θ(i, t)) applying a Hilbert transform and we bandpass filter the parcellated fMRI time-series within 0.01 − 0.1Hz [see e.g., Popovych et al. (2021) and references therein] using a discrete Fourier transform computed with a fast Fourier transform. Then, the phase coherence between brain areas i and j at time t, dFC(i, j, t) is defined as follows:
When two ROIs have temporarily aligned BOLD signals their respective dFC(i, j, t)≈1 while their BOLD signals are orthogonal dFC(i, j, t)≈0. Notably, the matrix dFC serves as the foundation of Leading Eigenvector Dynamic Analysis (LEiDA) which has been used to detect subtle FC patterns that distinguish between healthy and diseased BOLD signals (Cabral et al., 2017; Hancock et al., 2022).
In Figure 4, we present the Phase Coherence Connectivity PSE for edFC vs. sdFC (in a similar way as in Figure 2). However, we now compare each simulated mean dFC (sdFC) calculated by Equation (5) with the empirical (edFC) ones using the Pearson Correlation Coefficient from the upper triangular section of the two respective matrices as follows:
The optimal match between sdFC and edFC in the parameter space is acquired for (K, τ)-values where CCdFC becomes maximal (in Figures 4A, D, we indicate five maximum values with white circles). The upper row shows the results when we use the standard SC matrix of the respective subject to define the weights in model Equation (1). Figure 2A shows CCdFC = corr(sdFC, edFC) for the parameters (K, τ), while Figure 4B shows that the edFC calculated from the empirical BOLD signal. Figure 4C shows the sdFC matrix obtained by the larger CCdFC. In the lower row, we show the same analysis using the the hybrid rsSC matrix of the respective subject. Once again, in the context of dynamical functional connectivity, we find that the use of hybrid rsSC yields substantial improvement in the best fit between empirical and simulated BOLD activities (CCdFC(rsSC)≈0.72 while CCdFC(SC)≈0.29). In Figures 4G, H, we present the corresponding correlation analysis and scatterplots and conclusions as those found earlier (Figure 4). Here, we have presented the output for the same example subject (as the one in previous figures). However, this conclusion holds for all subjects similar to what we show in Figure 2, in Figures 4I, J, we show the respective statistical analysis and boxplots.
Figure 4. Parameter sweep exploration for phase coherence connectivity (edFC vs. sdFC). First row (using the respective subject's standard SC matrix to define the weights in model Equation (1): (A) The colormap depicts the CCdFC = corr(sdFC, edFC) for the parameters (K, τ). (B) The edFC calculated from the empirical BOLD signal. (C) The sFC matrix with the larger CCdFC. Second row (using the respective subject's hybrid rsSC matrix to define the weights in the model. (D) The respective colormap for CCdFC = corr(sdFC, edFC). (E) The edFC calculated from the empirical BOLD signal (same as (B)). (E) The sdFC matrix with the larger CCdFC. Note the different ranges in the respective colorbars of (B–F) capturing the different temporal alignment of the phase of the respective BOLD signals. See text for more details. (G, H) Scatterplots between empirical (y−axis) and optimal simulated dFC correlations (x−axis) aggregated across all entries in the corresponding dFC matrices, i.e., (B, C). (G) edFC vs. sdFC using rsSC (H) edFC vs. sdFC using standard SC. Both panels refer to same subject (blue lines indicate the linear regression model). (I, J) Boxplots for the non-carriers (I) and carriers (J) dataset (38 subjects per type) correlation coefficients between edFC and sdFC using SC and rsSC matrices for the simulated time series respectively. For each subject we considered the 5 maximum values (see circles in (A, D). The difference in the respective mean values of the two datasets is statistically significant measured by the t-test with very small p−value (p ≤ 0.0001) for both non-carriers and carriers sets.
4 Discussion
In this study, we showed that a coupled Kuramoto oscillator system built on a novel brain connectome can yield simulated BOLD brain activities that strongly resemble actual BOLD signals observed during resting-state fMRI. We used the TVB computational platform with the Kuramoto model (Kuramoto, 2003) and generated simulated BOLD time series across a range of different model parameters (K, τ) (producing PSE colormaps like in Figure 2). This allowed us to optimize model parameters and tune generated synthetic BOLD signals that produce simulated functional connectivity (FC) most similar to actual observed FC. Overall, we found that there are important advantages in using hybrid rsSC as it can produce BOLD sequences and synthetic FC that follow the general trends of the empirical BOLD time series and empirical FC (Figures 2, 4).
Despite the fact that, in general, both sFC (simulated with rsSC/SC matrices) perform rather well in capturing the positive correlations observed in the empirical BOLD signals, only the rsSC ones can effectively produce negative correlations, closely matching those occurring in the empirical BOLD signals (Figure 3, see also Zhan et al., 2017).
Our study has a few limitations. First, we restricted ourselves to a specific frequency band during simulations, and thus, future studies should further explore different ranges of frequencies in the Kuramoto model, e.g., either in different Hz ranges or extracted directly from the empirical BOLD signals per node and per subject [see Lee and Frangou (2017); Popovych et al. (2021) and references therein]. Furthermore, one may validate these findings for different dynamical models or to consider additional relevant dynamical features such as noise or the use of neuroimaging data where the path lengths are also available. In the study by (Popovych et al., 2021), the authors compared BOLD simulated signals (with SC) obtained using different dynamical models, such as the Kuramoto phase oscillators and the Hopf limit-cycle oscillators. They reported that both models perform rather similarly and that the role of such a model is not crucial as well as differences in the quality of the simulated optimal BOLD signals when using different atlases (structural vs functional) and parcellations. In our study, we achieved a significantly better agreement between optimal sFC and eFC compared the ones reported in the literature (Popovych et al., 2021). Let us also stress that in this study, we do not seek to detect model parameter settings that could distinguish between carriers and non-carriers based on the presence or absence of the APOE ε4 gene or age and gender factors, which is a research direction we plan to take in the near future.
In summary, here, we showed that our recently proposed hybrid connectome rsSC can produce simulated synthetic BOLD signals that yield functional connectivity matrices, which are strikingly similar to those actually obtained during the resting state. Thus, we conclude by highlighting that existing publicly available open-source pipelines, such as the TVB platform, could be easily equipped to include an add-on module that incorporates rsSC for the neuroscientific community interested in the modeling of simulated fMRI BOLD time series.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.
Ethics statement
This study represents a secondary data analysis using imaging data obtained as part of a larger study. The latter one was carried out in accordance with guidelines set by the institutional review boards at the University of Wisconsin-Milwaukee and the Medical College of Wisconsin. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.
Author contributions
TM: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing—original draft, Writing—review & editing. SD-P: Conceptualization, Methodology, Resources, Software, Writing—review & editing. IF: Data curation, Software, Writing—review & editing. ID: Data curation, Formal analysis, Software, Writing—review & editing. LZ: Data curation, Software, Writing—review & editing. AL: Conceptualization, Formal analysis, Funding acquisition, Methodology, Software, Supervision, Writing—review & editing.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was partially funded by the Helmholtz Association through the Helmholtz Portfolio Theme “Supercomputing and Modeling for the Human Brain”. This project was also received funding from the European Union's Horizon 2020 Research and Innovation Program under grant agreement No. 945539 (Human Brain Project SGA3). TM, SD-P, and AL were also supported by the Labex MME DII (ANR-11-LBX-410 0023-01) French national funding program. Additionally, LZ and AL were partially supported by NIH RF1MH125928 and R01AG071243. Moreover, LZ was also partially supported by NSF IIS 2045848. Open access publication funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation-491111487).
Conflict of interest
AL is a consultant for Otsuka USA, and serves on the medical advisory board of Buoy Health. AL was a co-founder of Keywise AI.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.
Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2023.1295395/full#supplementary-material
References
Ajilore, O., Zhan, L., GadElkarim, J., Zhang, A., Feusner, J., Yang, S., et al. (2013). Constructing the resting state structural connectome. Front. Neuroinform. 7, 30. doi: 10.3389/fninf.2013.00030
Ashourvan, A., Gu, S., Mattar, M. G., Vettel, J. M., and Bassett, D. S. (2017). The energy landscape underpinning module dynamics in the human brain connectome. NeuroImage 157, 364–380. doi: 10.1016/j.neuroimage.2017.05.067
Bansal, K., Garcia, J. O., Tompson, S. H., Verstynen, T., Vettel, J. M., and Muldoon, S. F. (2019). Cognitive chimera states in human brain networks. Sci. Adv. 5, eaau8535. doi: 10.1126/sciadv.aau8535
Beggs, J. M., and Plenz, D. (2003). Neuronal avalanches in neocortical circuits. J. Neurosci. 23, 11167–11177. doi: 10.1523/JNEUROSCI.23-35-11167.2003
Biswal, B. B., Kylen, J. V., and Hyde, J. S. (1997). Simultaneous assessment of flow and bold signals in resting-state functional connectivity maps. NMR Biomed. 10, 165–170. doi: 10.1002/(SICI)1099-1492(199706/08)10:4/5<165::AID-NBM454>3.0.CO;2-7
Cabral, J., Hugues, E., Sporns, O., and Deco, G. (2011). Role of local network oscillations in resting-state functional connectivity. NeuroImage 57, 130–139. doi: 10.1016/j.neuroimage.2011.04.010
Cabral, J., Vidaurre, D., Marques, P., Magalh aes, R., Silva Moreira, P., Miguel Soares, J., et al. (2017). Cognitive performance in healthy older adults relates to spontaneous switching between states of functional connectivity during rest. Scient. Rep. 7, 5135. doi: 10.1038/s41598-017-05425-7
Cocchi, L., Gollo, L. L., Zalesky, A., and Breakspear, M. (2017). Criticality in the brain: a synthesis of neurobiology, models and cognition. Progr. Neurobiol. 158, 132–152. doi: 10.1016/j.pneurobio.2017.07.002
Cocco, S., Monasson, R., Posani, L., and Tavoni, G. (2017). Functional networks from inverse modeling of neural population activity. Curr. Opin. Syst. Biol. 3, 103–110. doi: 10.1016/j.coisb.2017.04.017
Deco, G., and Jirsa, V. K. (2012). Ongoing cortical activity at rest: criticality, multistability, and ghost attractors. J. Neurosci. 32, 3366–3375. doi: 10.1523/JNEUROSCI.2523-11.2012
Deco, G., Jirsa, V. K., and McIntosh, A. R. (2011). Emerging concepts for the dynamical organization of resting-state activity in the brain. Nat. Rev. Neurosci. 12, 43–56. doi: 10.1038/nrn2961
Deco, G., Jirsa, V. K., Robinson, P. A., Breakspear, M., and Friston, K. (2008). The dynamic brain: From spiking neurons to neural masses and cortical fields. PLoS Computat. Biol. 4, 1–35. doi: 10.1371/journal.pcbi.1000092
Deco, G., Sanz Perl, Y., Vuust, P., Tagliazucchi, E., Kennedy, H., and Kringelbach, M. L. (2021). Rare long-range cortical connections enhance human information processing. Curr. Biol. 31, 4436–4448.e5. doi: 10.1016/j.cub.2021.07.064
Desikan, R. S., Sgonne, F., Fischl, B., Quinn, B. T., Dickerson, B. C., Blacker, D., et al. (2006). An automated labeling system for subdividing the human cerebral cortex on mri scans into gyral based regions of interest. NeuroImage 31, 968–980. doi: 10.1016/j.neuroimage.2006.01.021
Ercsey-Ravasz, M., Markov, N., Lamy, C., Van Essen, D., Knoblauch, K., Toroczkai, Z., et al. (2013). A predictive network model of cerebral cortical connectivity based on a distance rule. Neuron 80, 184–197. doi: 10.1016/j.neuron.2013.07.036
Ezaki, T., Fonseca dos Reis, E., Watanabe, T., Sakaki, M., and Masuda, N. (2020). Closer to critical resting-state neural dynamics in individuals with higher fluid intelligence. Commun. Biol. 3, 52. doi: 10.1038/s42003-020-0774-y
Ezaki, T., Watanabe, T., Ohzeki, M., and Masuda, N. (2017). Energy landscape analysis of neuroimaging data. Philosoph. Trans. R. Soc. A. 375, 20160287. doi: 10.1098/rsta.2016.0287
Finn, E. S., Shen, X., Scheinost, D., Rosenberg, M. D., Huang, J., Chun, M. M., et al. (2015). Functional connectome fingerprinting: identifying individuals using patterns of brain connectivity. Nat. Neurosci. 18, 1664–1671. doi: 10.1038/nn.4135
Fortel, I., Butler, M., Korthauer, L. E., Zhan, L., Ajilore, O., Driscoll, I., et al. (2019). “Brain dynamics through the lens of statistical mechanics by unifying structure and function,” in Medical Image Computing and Computer Assisted Intervention-MICCAI 2019, eds. D. Shen, T. Liu, T. M. Peters, L. H. Staib, E. Carolineand, S. Zhou, et al. (Cham. Springer International Publishing), 503–511. doi: 10.1007/978-3-030-32254-0_56
Fortel, I., Butler, M., Korthauer, L. E., Zhan, L., Ajilore, O., Sidiropoulos, A., et al. (2022). Inferring excitation-inhibition dynamics using a maximum entropy model unifying brain structure and function. Netw. Neurosci. 6, 420–444. doi: 10.1162/netn_a_00220
Fortel, I., Korthauer, L. E., Morrissey, Z., Zhan, L., Ajilore, O., Wolfson, O., et al. (2020). Connectome Signatures of hyperexcitation in cognitively intact middle-aged female APOE-ε4 carriers. Cerebral Cortex 30, 6350–6362. doi: 10.1093/cercor/bhaa190
Fortel, I., Zhan, L., Ajilore, O., Wu, Y., Mackin, S., and Leow, A. (2023). Disrupted Excitation-Inhibition balance in cognitively normal individuals at risk of alzheimer's disease. J. Alzheimers Dis. 95, 1449–1467. doi: 10.3233/JAD-230035
Friston, K., Harrison, L., and Penny, W. (2003). Dynamic causal modelling. NeuroImage 19, 1273–1302. doi: 10.1016/S1053-8119(03)00202-7
Ghosh, A., Rho, Y., McIntosh, A. R., Kötter, R., and Jirsa, V. K. (2008). Noise during rest enables the exploration of the brain's dynamic repertoire. PLoS Comput. Biol. 4, 1–12. doi: 10.1371/journal.pcbi.1000196
Hagmann, P., Cammoun, L., Gigandet, X., Gerhard, S., Ellen Grant, P., Wedeen, V., et al. (2010). Mr connectomics: Principles and challenges. J. Neurosci. Method 194, 34–45. doi: 10.1016/j.jneumeth.2010.01.014
Hahn, G., Ponce-Alvarez, A., Monier, C., Benvenuti, G., Kumar, A., Chavane, F., et al. (2017). Spontaneous cortical activity is transiently poised close to criticality. PLoS Computat. Biol. 13, 1–29. doi: 10.1371/journal.pcbi.1005543
Haimovici, A., Tagliazucchi, E., Balenzuela, P., and Chialvo, D. R. (2013). Brain organization into resting state networks emerges at criticality on a model of the human connectome. Phys. Rev. Lett. 110, 178101. doi: 10.1103/PhysRevLett.110.178101
Hancock, F., Cabral, J., Luppi, A. I., Rosas, F. E., Mediano, P. A., Dipasquale, O., et al. (2022). Metastability, fractal scaling, and synergistic information processing: what phase relationships reveal about intrinsic brain activity. NeuroImage 259, 119433. doi: 10.1016/j.neuroimage.2022.119433
Horn, A., Ostwald, D., Reisert, M., and Blankenburg, F. (2014). The structural-functional connectome and the default mode network of the human brain. NeuroImage, 102, 142–151. doi: 10.1016/j.neuroimage.2013.09.069
Jirsa, V., Proix, T., Perdikis, D., Woodman, M., Wang, H., Gonzalez-Martinez, J., et al. (2017). The virtual epileptic patient: Individualized whole-brain models of epilepsy spread. NeuroImage 145, 377–388. doi: 10.1016/j.neuroimage.2016.04.049
Kadirvelu, B., Hayashi, Y., and Nasuto, S. J. (2017). Inferring structural connectivity using ising couplings in models of neuronal networks. Scient. Rep. 7, 8156. doi: 10.1038/s41598-017-05462-2
Kinouchi, O., and Copelli, M. (2006). Optimal dynamical range of excitable networks at criticality. Nat. Phys. 2, 348–351. doi: 10.1038/nphys289
Korthauer, L., Zhan, L., Ajilore, O., Leow, A., and Driscoll, I. (2018). Disrupted topology of the resting state structural connectome in middle-aged apoe ε4 carriers. NeuroImage 178, 295–305. doi: 10.1016/j.neuroimage.2018.05.052
Kuramoto, Y. (2003). Chemical Oscillations, Waves, and Turbulence. Dover Books on Chemistry Series. Mineola, NY: Dover Publications.
Kuznetsov, Y. A. (1998). Elements of Applied Bifurcation Theory. Number 112 in Applied Mathematical Sciences. Berlin: Springer.
Lee, W. H., and Frangou, S. (2017). Linking functional connectivity and dynamic properties of resting-state networks. Scient. Rep. 7, 16610. doi: 10.1038/s41598-017-16789-1
Lombardi, F., Herrmann, H. J., and de Arcangelis, L. (2017). Balance of excitation and inhibition determines 1/f power spectrum in neuronal networks. Chaos. 27, 047402. doi: 10.1063/1.4979043
Marinazzo, D., Pellicoro, M., Wu, G., Angelini, L., Corts, J. M., and Stramaglia, S. (2014). Information transfer and criticality in the ising model on the human connectome. PLoS ONE 9, 1–7. doi: 10.1371/journal.pone.0093616
Mess, A., Rudrauf, D., Benali, H., and Marrelec, G. (2014). Relating structure and function in the human brain: Relative contributions of anatomy, stationary dynamics, and non-stationarities. PLoS Comput. Biol. 10, 1–9. doi: 10.1371/journal.pcbi.1003530
Miller, K. J., Weaver, K. E., and Ojemann, J. G. (2009). Direct electrophysiological measurement of human default network areas. Proc. Natl. Acad. Sci. U S A. 106, 12174–12177. doi: 10.1073/pnas.0902071106
Muñoz, M. A. (2018). Colloquium: Criticality and dynamical scaling in living systems. Rev. Mod. Phys. 90, 031001. doi: 10.1103/RevModPhys.90.031001
Murray, J. D., Demirta, M., and Anticevic, A. (2018). Biophysical modeling of large-scale brain dynamics and applications for computational psychiatry. Biol. Psychiat. 3, 777–787. doi: 10.1016/j.bpsc.2018.07.004
Nghiem, T.-A., Telenczuk, B., Marre, O., Destexhe, A., and Ferrari, U. (2018). Maximum-entropy models reveal the excitatory and inhibitory correlation structures in cortical neuronal activity. Phys. Rev. E 98, 012402. doi: 10.1103/PhysRevE.98.012402
Niessing, J., Ebisch, B., Schmidt, K. E., Niessing, M., Singer, W., and Galuske, R. A. W. (2005). Hemodynamic signals correlate tightly with synchronized gamma oscillations. Science 309, 948–951. doi: 10.1126/science.1110948
Nir, Y., Fisch, L., Mukamel, R., Gelbard-Sagiv, H., Arieli, A., Fried, I., et al. (2007). Coupling between neuronal firing rate, gamma LfP, and bold fMRI is related to interneuronal correlations. Curr. Biol. 17, 1275–1285. doi: 10.1016/j.cub.2007.06.066
Niu, W., Huang, X., Xu, K., Jiang, T., and Yu, S. (2019). Pairwise interactions among brain regions organize large-scale functional connectivity during execution of various tasks. Neuroscience 412, 190–206. doi: 10.1016/j.neuroscience.2019.05.011
Nuzzi, D., Pellicoro, M., Angelini, L., Marinazzo, D., and Stramaglia, S. (2020). Synergistic information in a dynamical model implemented on the human structural connectome reveals spatially distinct associations with age. Netw. Neurosci. 4, 910–924. doi: 10.1162/netn_a_00146
Popovych, O. V., Jung, K., Manos, T., Diaz-Pier, S., Hoffstaedter, F., Schreiber, J., et al. (2021). Inter-subject and inter-parcellation variability of resting-state whole-brain dynamical modeling. NeuroImage 236, 118201. doi: 10.1016/j.neuroimage.2021.118201
Popovych, O. V., Manos, T., Hoffstaedter, F., and Eickhoff, S. B. (2019). What can computational models contribute to neuroimaging data analytics? Front. Syst. Neurosci. 12, 68. doi: 10.3389/fnsys.2018.00068
Rabuffo, G., Fousek, J., Bernard, C., and Jirsa, V. (2021). Neuronal cascades shape whole-brain functional dynamics at rest. eNeuro 8, 1523. doi: 10.1523/ENEURO.0283-21.2021
Roudi, Y., Tyrcha, J., and Hertz, J. (2009). Ising model for neural data: model quality and approximate methods for extracting functional connectivity. Phys. Rev. E 79, 051915. doi: 10.1103/PhysRevE.79.051915
Sanz-Leon, P., Knock, S. A., Spiegler, A., and Jirsa, V. K. (2015). Mathematical framework for large-scale brain network modeling in the virtual brain. NeuroImage 111, 385–430. doi: 10.1016/j.neuroimage.2015.01.002
Schneidman, E., Berry, M. J., Segev, R., and Bialek, W. (2006). Weak pairwise correlations imply strongly correlated network states in a neural population. Nature 440, 1007–1012. doi: 10.1038/nature04701
Schölvinck, M. L., Maier, A., Ye, F. Q., Duyn, J. H., and Leopold, D. A. (2010). Neural basis of global resting-state fMRI activity. Proc. Natl. Acad. Sci. U S A. 107, 10238–10243. doi: 10.1073/pnas.0913110107
Shew, W. L., Yang, H., Petermann, T., Roy, R., and Plenz, D. (2009). Neuronal avalanches imply maximum dynamic range in cortical networks at criticality. J. Neurosci. 29, 15595–15600. doi: 10.1523/JNEUROSCI.3864-09.2009
Shew, W. L., Yang, H., Yu, S., Roy, R., and Plenz, D. (2011). Information capacity and transmission are maximized in balanced cortical networks with neuronal avalanches. J. Neurosci. 31, 55–63. doi: 10.1523/JNEUROSCI.4637-10.2011
Shlens, J., Field, G. D., Gauthier, J. L., Grivich, M. I., Petrusca, D., Sher, A., et al. (2006). The structure of multi-neuron firing patterns in primate retina. J. Neurosci. 26, 8254–8266. doi: 10.1523/JNEUROSCI.1282-06.2006
Sornette, D. (2004). Critical Phenomena in Natural Sciences: Chaos, Fractals, Selforganization, and Disorder: Concepts and Tools. Berlin: Springer.
Sporns, O., Tononi, G., and Kötter, R. (2005). The human connectome: a structural description of the human brain. PLoS Comput. Biol. 1, e42. doi: 10.1371/journal.pcbi.0010042
Tagliazucchi, E. (2017). The signatures of conscious access and its phenomenology are consistent with large-scale brain communication at criticality. Consc. Cogn. 55, 136–147. doi: 10.1016/j.concog.2017.08.008
Tagliazucchi, E., Balenzuela, P., Fraiman, D., and Chialvo, D. (2012). Criticality in large-scale brain fmri dynamics unveiled by a novel point process analysis. Front. Physiol. 3, 15. doi: 10.3389/fphys.2012.00015
Tang, H., Ma, G., He, L., Huang, H., and Zhan, L. (2021). Commpool: an interpretable graph pooling framework for hierarchical graph representation learning. Neural Netw. 143, 669–677. doi: 10.1016/j.neunet.2021.07.028
Tkacik, G., Mora, T., Marre, O., Amodei, D., Palmer, S. E., Berry, M. J., et al. (2015). Thermodynamics and signatures of criticality in a network of neurons. Proc. Natl. Acad. Sci. 112, 11508–11513. doi: 10.1073/pnas.1514188112
Váča, F., Shanahan, M., Hellyer, P. J., Scott, G., Cabral, J., and Leech, R. (2015). Effects of lesions on synchrony and metastability in cortical networks. NeuroImage 118, 456–467. doi: 10.1016/j.neuroimage.2015.05.042
Watanabe, T., Hirose, S., Wada, H., Imai, Y., Machida, T., Shirouzu, I., et al. (2013). A pairwise maximum entropy model accurately describes resting-state human brain networks. Nat. Commun. 4, 1370. doi: 10.1038/ncomms2388
Wilting, J., and Priesemann, V. (2019). 25 years of criticality in neuroscience–established results, open controversies, novel concepts. Curr. Opin. Neurobiol. 58, 105–111. doi: 10.1016/j.conb.2019.08.002
Young, L.-S. (2020). Towards a mathematical model of the brain. J. Statist. Phys. 180, 612–629. doi: 10.1007/s10955-019-02483-1
Keywords: whole brain dynamics, resting-state brain dynamics, neuroimaging data, functional connectivity, resting-state informed structural connectome, Alzheimer's disease
Citation: Manos T, Diaz-Pier S, Fortel I, Driscoll I, Zhan L and Leow A (2023) Enhanced simulations of whole-brain dynamics using hybrid resting-state structural connectomes. Front. Comput. Neurosci. 17:1295395. doi: 10.3389/fncom.2023.1295395
Received: 16 September 2023; Accepted: 05 December 2023;
Published: 19 December 2023.
Edited by:
Alessandro E. P. Villa, Neuro-heuristic Research Group (NHRG), SwitzerlandReviewed by:
Richard Betzel, University of Pennsylvania, United StatesVictor Buendia, University of Granada, Spain
Copyright © 2023 Manos, Diaz-Pier, Fortel, Driscoll, Zhan and Leow. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Thanos Manos, dGhhbm9zLm1hbm9zJiN4MDAwNDA7Y3l1LmZy; Sandra Diaz-Pier, cy5kaWF6JiN4MDAwNDA7ZnotanVlbGljaC5kZQ==