Thanos Manos
email: thanos.manos@cyu.fr

Sandra Diaz-Pier

Igor Fortel

Ira Driscoll

Liang Zhan

Alex Leow

Enhanced simulations of whole-brain dynamics using hybrid resting-state structural connectomes

Keywords: Whole brain dynamics, Resting-state brain dynamics, neuroimaging data, functional connectivity, resting-state informed structural connectome, Alzheimer's disease

The human brain, composed of billions of neurons and synaptic connections, is an intricate network coordinating a sophisticated balance of excitatory and inhibitory activity 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, 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 well aligns with the observed FC when compared to that simulated with traditional structural connectome.

INTRODUCTION

The human brain is a complex neural network that self-organises into different emergent states, crucial for its functions. Such states include spatiotemporal patterns of neural synchronisation associated with cognitive processes [START_REF] Bansal | Cognitive chimera states in human brain networks[END_REF]. Brain regions can be modelled as dynamically interacting nodes in a functional network on a 3D space (functional brain networks), coupled in a complex manner driven by the structure of these networks. Over the past years, interdisciplinary approaches using concepts from nonlinear 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 2 to 3 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 [START_REF] Finn | Functional connectome fingerprinting: identifying individuals using patterns of brain connectivity[END_REF].

Together with extensive experimental work, mathematical/computational modelling of the whole brain dynamics has been an active research topic for years (e.g., [START_REF] Deco | The dynamic brain: From spiking neurons to neural masses and cortical fields[END_REF][START_REF] Sanz-Leon | Mathematical framework for large-scale brain network modeling in the virtual brain[END_REF][START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF][START_REF] Murray | Biophysical modeling of large-scale brain dynamics and applications for computational psychiatry[END_REF]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 (see e.g., [START_REF] Ghosh | Noise during rest enables the exploration of the brain's dynamic repertoire[END_REF][START_REF] Hagmann | Mr connectomics: Principles and challenges[END_REF][START_REF] Deco | Emerging concepts for the dynamical organization of resting-state activity in the brain[END_REF]). 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 etc.). 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 (see e.g., [START_REF] Sporns | The human connectome: A structural description of the human brain[END_REF][START_REF] Horn | The structural-functional connectome and the default mode network of the human brain[END_REF].

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 with the simulated FC (see e.g., [START_REF] Cabral | Role of local network oscillations in resting-state functional connectivity[END_REF][START_REF] Deco | Ongoing cortical activity at rest: Criticality, multistability, and ghost attractors[END_REF]). Hence, these parameters can serve as dynamical biomarkers and predictors of different brain states and behavioural modes (see [START_REF] Popovych | What can computational models contribute to neuroimaging data analytics[END_REF] for a recent review). Along this direction, the virtual epileptic patient has been recently proposed, where medicaltreatment approaches using personalised mathematical models for epileptic patients have been illustrated (see e.g., [START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF]). 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 [START_REF] Popovych | Intersubject and inter-parcellation variability of resting-state whole-brain dynamical modeling[END_REF] and references therein for more details).

In recent years, substantial research efforts have been directed toward understanding the brain (largescale activity) using resting state fMRI (rs-fMRI) employing sophisticated mathematical and statistical tools to investigate the FC from rs-fMRI data [START_REF] Biswal | Simultaneous assessment of flow and bold signals in resting-state functional connectivity maps[END_REF]. So far, the mainstream approach is to consider SC to be static and the FC one dynamic. However, this is not necessarily the case as 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 connectivityinformed structural connectivity" has been introduced in [START_REF] Ajilore | Constructing the resting state structural connectome[END_REF] employing information from fMRI to infer the underlying pattern of white matter engagement specific to the brain's state. 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 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" versus "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 [START_REF] Kinouchi | Optimal dynamical range of excitable networks at criticality[END_REF]) such that information capacity and transfer are maximized [START_REF] Shew | Information capacity and transmission are maximized in balanced cortical networks with neuronal avalanches[END_REF]. This is the central thesis of the criticality hypothesis [START_REF] Beggs | Neuronal avalanches in neocortical circuits[END_REF][START_REF] Muñoz | Colloquium: Criticality and dynamical scaling in living systems[END_REF], i.e., that brain activity self-organize into a critical state [START_REF] Wilting | 25 years of criticality in neuroscience -established results, open controversies, novel concepts[END_REF], a unique configuration likened to a phase transition in physical systems where a dynamical system transitions from order (balanced excitation-inhibition) to disorder (disrupted excitationinhibition balance) [START_REF] Cocchi | Criticality in the brain: A synthesis of neurobiology, models and cognition[END_REF][START_REF] Hahn | Spontaneous cortical activity is transiently poised close to criticality[END_REF][START_REF] Sornette | Critical phenomena in natural sciences : chaos, fractals, selforganization, and disorder : concepts and tools[END_REF][START_REF] Tagliazucchi | The signatures of conscious access and its phenomenology are consistent with large-scale brain communication at criticality[END_REF]. Indeed, evidence supporting that the brain is operating near criticality has been reported in studies examining neuronal signaling [START_REF] Hahn | Spontaneous cortical activity is transiently poised close to criticality[END_REF][START_REF] Beggs | Neuronal avalanches in neocortical circuits[END_REF][START_REF] Shew | Neuronal avalanches imply maximum dynamic range in cortical networks at criticality[END_REF] as well as BOLD fMRI signals [START_REF] Haimovici | Brain organization into resting state networks emerges at criticality on a model of the human connectome[END_REF][START_REF] Lombardi | Balance of excitation and inhibition determines 1/f power spectrum in neuronal networks[END_REF][START_REF] Rabuffo | Neuronal cascades shape whole-brain functional dynamics at rest[END_REF][START_REF] Tagliazucchi | Criticality in large-scale brain fmri dynamics unveiled by a novel point process analysis[END_REF].

To incorporate co-activation (excitatory) or silencing (inhibitory) effects into our hybrid rsSC framework that would allow us to infer the brain's E-I balance, in [START_REF] Fortel | Brain dynamics through the lens of statistical mechanics by unifying structure and function[END_REF] we then introduced an improved framework based on the Ising model representation of the brain as a dynamical system, wherein selforganized 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 [START_REF] Kadirvelu | Inferring structural connectivity using ising couplings in models of neuronal networks[END_REF][START_REF] Tkačik | Thermodynamics and signatures of criticality in a network of neurons[END_REF] and macroscale interactions [START_REF] Marinazzo | Information transfer and criticality in the ising model on the human connectome[END_REF][START_REF] Ezaki | Energy landscape analysis of neuroimaging data[END_REF][START_REF] Nghiem | Maximum-entropy models reveal the excitatory and inhibitory correlation structures in cortical neuronal activity[END_REF][START_REF] Niu | Pairwise interactions among brain regions organize large-scale functional connectivity during execution of various tasks[END_REF][START_REF] Nuzzi | Synergistic information in a dynamical model implemented on the human structural connectome reveals spatially distinct associations with age[END_REF][START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF] of the human brain, and to accurately represent spatiotemporal co-activations in neuronal spike trains [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF][START_REF] Roudi | Ising model for neural data: Model quality and approximate methods for extracting functional connectivity[END_REF][START_REF] Shlens | The structure of multi-neuron firing patterns in primate retina[END_REF] and patterns of BOLD activity [START_REF] Ashourvan | The energy landscape underpinning module dynamics in the human brain connectome[END_REF][START_REF] Cocco | Functional networks from inverse modeling of neural population activity[END_REF][START_REF] Ezaki | Closer to critical resting-state neural dynamics in individuals with higher fluid intelligence[END_REF][START_REF] Watanabe | A pairwise maximum entropy model accurately describes resting-state human brain networks[END_REF].

In this paper, 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 [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF] and the generic limit-cycle oscillators [START_REF] Kuznetsov | Elements of applied bifurcation theory[END_REF], to model each node's mean neural activity 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 (see e.g. [START_REF] Popovych | Intersubject and inter-parcellation variability of resting-state whole-brain dynamical modeling[END_REF]). We here show that by using rsSC such dynamical systems succeed to produce simulated signals with both positive and negative correlations which sufficiently follow the trends of the empirical ones.

METHODS AND MATERIALS

Empirical data & 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 [START_REF] Korthauer | Disrupted topology of the resting state structural connectome in middle-aged apoe ε4 carriers[END_REF]. The algorithm to obtain rsSC can be found in [START_REF] Fortel | Brain dynamics through the lens of statistical mechanics by unifying structure and function[END_REF][START_REF] Fortel | Connectome Signatures of Hyperexcitation in Cognitively Intact Middle-Aged Female APOE-ε4 Carriers[END_REF][START_REF] Tang | Commpool: An interpretable graph pooling framework for hierarchical graph representation learning[END_REF][START_REF] Fortel | Inferring excitation-inhibition dynamics using a maximum entropy model unifying brain structure and function[END_REF] while in the Supplementary Material we provide a concise description with details on its implementation.

Models and simulated data

In order to produce simulated fMRI time series in the given connectomes, we employ the Kuramoto phase oscillator model [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF][START_REF] Lee | Linking functional connectivity and dynamic properties of restingstate networks[END_REF][START_REF] Popovych | Intersubject and inter-parcellation variability of resting-state whole-brain dynamical modeling[END_REF]:

θi (t) = 2πf i + K N N j=1 c ij sin[θ j (t -τ ij ) -θ i (t)], i = 1, 2, ..., N (1)
where θ i are the phases, N is the number of oscillators, f j are the natural frequencies (Hz), c ij 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 Eq. ( 1) to simulate the network's dynamics and to calculate time series. To this end, two cases of connectivity matrices were compared: (i) in the first one the c ij 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 c ij 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 c ij 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 = L ij /V , where L ij (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 instead the euclidean distance between nodes in the Desikan atlas [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on mri scans into gyral based regions of interest[END_REF] as proxies. The euclidean distance has been used in the literature in the construction of structural networks (see e.g. [START_REF] Ercsey-Ravasz | A predictive network model of cerebral cortical connectivity based on a distance rule[END_REF]) and found to closely follow the trends obtained by anatomical tract-tracing studies. Furthermore in [START_REF] Deco | Rare long-range cortical connections enhance human information processing[END_REF], the authors showed that such networks also strongly correlate with MRI tractography-based networks. The matrix L = L ij can thus be used to calculate the delays τ ij in the coupling, which can be expressed as

τ ij = τ • L ij /⟨L ij ⟩,
where τ = ⟨L ij ⟩/V is the global (or mean) delay. In Eq. (1) the self-connections were excluded by setting the diagonal elements in the matrices eSC/rsSC and L to zero (i.e. c ii = L ii = 0 respectively).

In Figure 1(A), 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 1(B) depicts the corresponding to this SC and subject rsSC matrix calculated as described earlier. Note that the hybrid rsSC contains negative entry values as opposed to SC one that is restricted to having only positive values. Figure 1(C) shows the tract length L matrix (in mm) that we used for all subjects' simulation in the absence of the actual measured ones from a neuroimaging prepossessing pipeline. The phases in our model (Eq. 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 e.g. [START_REF] Cabral | Role of local network oscillations in resting-state functional connectivity[END_REF][START_REF] Messé | Relating structure and function in the human brain: Relative contributions of anatomy, stationary dynamics, and non-stationarities[END_REF][START_REF] Váša | Effects of lesions on synchrony and metastability in cortical networks[END_REF][START_REF] Lee | Linking functional connectivity and dynamic properties of restingstate networks[END_REF] 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 (see also [START_REF] Niessing | Hemodynamic signals correlate tightly with synchronized gamma oscillations[END_REF][START_REF] Nir | Coupling between neuronal firing rate, gamma lfp, and bold fmri is related to interneuronal correlations[END_REF][START_REF] Miller | Direct electrophysiological measurement of human default network areas[END_REF][START_REF] Schölvinck | Neural basis of global resting-state fMRI activity[END_REF]).

For our simulations, we used a TVB tailored version for the Kuramoto model and we made adjustments for efficient parallelisation on CPUs using MPI on the supercomputer JUSUF located at the Jülich Supercomputing Centre. Our model generates time series which correspond initially to electrical activity (fast oscillations) for each node, i.e. we register the observable x i = sin(θ i ) for each brain region. Then, in order to estimate the simulated BOLD signal, we use the TVB's built-in tool 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 [START_REF] Friston | Dynamic causal modelling[END_REF]. Our simulations ran for 500 seconds in total. The first 20 seconds were discarded to remove transient effects, resulting in T = 480 seconds (8 minutes), i.e. a time interval identical to the time-length of the empirical fMRI signals. We set the timestep 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.

RESULTS

We numerically simulate BOLD time series varying two model parameters, namely the global coupling strength K and the delay τ in Eq. ( 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 x i , i = 1, 2, ..., N from different ROIs (also referred to as Static Functional Connectivity in the literature, see e.g. [START_REF] Cabral | Cognitive performance in healthy older adults relates to spontaneous switching between states of functional connectivity during rest[END_REF]), namely:

CC BOLD = corr(x i , x j ).
(2) Then, we compare each sFC with the eFC ones using, again, the Pearson Correlation Coefficient, however this time we calculate it for the two respective matrices (upper triangular parts), i.e.: CC FC = corr(sFC,eFC).

(3)

The optimal match between sFC and eFC in the parameter space is acquired for (K, τ )-values where CC FC becomes maximal (see also [START_REF] Popovych | Intersubject and inter-parcellation variability of resting-state whole-brain dynamical modeling[END_REF] 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 like Eq. ( 1) which better approximate the empirical BOLD signals (shown here for one example healthy subject). The upper row refers to simulations performed using the respective subject's standard SC matrix to define the coupling weights in Kuramoto model. Figure 2(A) shows the parameter sweep exploration (PSE) for eFC vs. sFC and for the parameters (K, τ ) and measured as CC FC = corr(sFC,eFC). The 5 white circles on the red regions indicate the highest correlations found (larger circles' sizes correspond to larger CC FC values). In Figure 2(B) we present the eFC calculated from the empirical BOLD signal while in Figure 2(C) the sFC matrix with the larger cc FC . We can observe that sFC did not capture adequately the negative correlations that are present in eFC (compare the minimum values in the barplots of panel (B) and (C)).

In the second row of Figure 2 (panels (D),(E),(F)), we perform similar simulations, however we now use the respective subject's hybrid rsSC matrix to define the coupling weights in the Kuramoto model). Note the significant improvement in the maximum value of the CC FC ≈ 0.86 compared to the one found when using the standard SC matrix (CC FC ≈ 0.33). Note also the better agreement 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 Figure 2(A) and (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 CC FC (sFC,eFC) value (using hybrid rsSC matrices) is well reflected by a rather clear linear trend in the distribution of the points (panel (A)), On the other hand, only a relatively weak linear trend is obtained using standard SC matrices (panel (B)). Both panels refer to 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 noncarriers (I) and 38 carriers (J). For each subject we considered the 5 maximum values of correlation coefficients between eFC and sFC using SC and rsSC matrices for the simulated time series respectively (circles in (A),(D)) and produced boxplots. We then used the t-test to measure the difference in each group's (carriers or non-carriers) mean Pearson correlation value 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 the Supplementary Material (Figure S1). The separation of the 2 groups is indeed biologically important given the already demonstrated differences in E-I dynamics (see e.g. [START_REF] Fortel | Brain dynamics through the lens of statistical mechanics by unifying structure and function[END_REF][START_REF] Fortel | Connectome Signatures of Hyperexcitation in Cognitively Intact Middle-Aged Female APOE-ε4 Carriers[END_REF][START_REF] Fortel | Inferring excitation-inhibition dynamics using a maximum entropy model unifying brain structure and function[END_REF][START_REF] Fortel | Disrupted Excitation-Inhibition balance in cognitively normal individuals at risk of alzheimer's disease[END_REF]. 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 noncarriers group. We considered sFC matrices produced by the 5 parameters (K, τ ) that maximized eFC and sFC similarity per subject, (marked with circles in Figure 2(A),(D)). The leftmost three boxplots depict minimum (negative) correlations in actual eFC (light green), simulated using rsSC (light blue) and using standard SC (red) respectively, while the rightmost three boxplots depict maximum (positive) correlation values.

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-carriers and carriers datasets (see Figure S2 in the Supplementary Material). 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 the 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. (see e.g. [START_REF] Cabral | Cognitive performance in healthy older adults relates to spontaneous switching between states of functional connectivity during rest[END_REF][START_REF] Hancock | Metastability, fractal scaling, and synergistic information processing: What phase relationships reveal about intrinsic brain activity[END_REF] and references therein) which does not suffer from time-window length effects like other similar techniques based on calculating successive FC(t) matrices using a sliding-window (see discussion in [START_REF] Cabral | Cognitive performance in healthy older adults relates to spontaneous switching between states of functional connectivity during rest[END_REF][START_REF] Hancock | Metastability, fractal scaling, and synergistic information processing: What phase relationships reveal about intrinsic brain activity[END_REF]). 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 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. [START_REF] Popovych | Intersubject and inter-parcellation variability of resting-state whole-brain dynamical modeling[END_REF] 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: dFC(i, j, t) = cos(θ(i, t) -θ(j, t)).

(4)

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. Note the matrix dFC serves as the foundation of Leading Eigenvector Dynamic Analysis (LEiDA) which has been used to detect subtle FC patterns that distinguish healthy versus diseased BOLD signals (see e.g. [START_REF] Cabral | Cognitive performance in healthy older adults relates to spontaneous switching between states of functional connectivity during rest[END_REF][START_REF] Hancock | Metastability, fractal scaling, and synergistic information processing: What phase relationships reveal about intrinsic brain activity[END_REF]).

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 Eq. ( 5) with the empirical (edFC) ones using the Pearson Correlation Coefficient, from the upper triangular section of the two respective matrices, i.e.:

CC dFC = corr(sdFC,edFC).

(5)

The optimal match between sdFC and edFC in the parameter space is acquired for (K, τ )-values where CC dFC becomes maximal (in panels (A) and (D) we indicate 5 maximum values with white circles).

The upper row shows the results when we use the respective subject's standard SC matrix to define the weights in model Eq. (1). Figure 4(A) depicts the CC dFC = corr(sdFC,edFC) for the parameters (K, τ ), while Figure 4(B) the edFC calculated from the empirical BOLD signal. Figure 4(C) shows the sdFC matrix obtained by the larger CC dFC . In the lower row, we show the same analysis using now the respective subject's hybrid rsSC matrix. Once again, we find that the use of hybrid rsSC yields substantial improvement in the best fit between empirical and simulated BOLD activity (CC dFC (rsSC) ≈ 0.72 while CC dFC (SC) ≈ 0.29), this time in the context of dynamical functional connectivity. In Figures 4(G-H), we present the corresponding correlation analysis and scatterplots and conclusions as those found earlier (Figure 2). Here, we have presented the output for the same example subject (as the one in previous

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 [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF] 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 well the general trends of the empirical BOLD time series and empirical FC (Figures 2 and4).

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 [START_REF] Zhan | The significance of negative correlations in brain connectivity[END_REF]).

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 e.g. [START_REF] Lee | Linking functional connectivity and dynamic properties of restingstate networks[END_REF][START_REF] Popovych | Intersubject and inter-parcellation variability of resting-state whole-brain dynamical modeling[END_REF] and references therein). Furthermore one may validate these findings for different dynamical models, or to furthermore consider additional relevant dynamical features such as noise or the use of neuroimaging data where the path lengths is also available. In [START_REF] Popovych | Intersubject and inter-parcellation variability of resting-state whole-brain dynamical modeling[END_REF], the authors compared BOLD simulated signals (with SC) obtained using different dynamical models, namely 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 (e.g., [START_REF] Popovych | Intersubject and inter-parcellation variability of resting-state whole-brain dynamical modeling[END_REF]). Let us also stress that in this work we do not seek to detect model parameters settings that could distinguish between carriers and non-carriers based on the presence or not 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 strikingly similar to those actually obtained during the resting-state. Thus, we conclude by highlighting that existing publicly available opensource 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.

SUPPLEMENTARY MATERIAL Empirical data

Structural and functional connectivity for 38 cognitively normal APOE ε4 allele carriers aged 40-60 (µ = 50.8, σ = 0.99 ) are compared with 38 age (µ = 50.9, σ = 0.99) and sex-matched (16 male/22 female) non-carriers (control -non-carriers). Resting state functional MRI (rs-fMRI)-A T2*-weighted functional scan was obtained 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-minute rs-fMRI scan was acquired under a task-free condition (i.e., resting state): subjects were instructed to relax with eyes closed and to "not think about anything in particular". Imaging included T1-weighted MRI, resting state fMRI and diffusion weighted MRI. Freesurfer cortical parcellation and sub-cortical segmentation was performed to derive 80 regions-of-interest (ROIs) registered on the Desikan atlas [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on mri scans into gyral based regions of interest[END_REF]. The mean time-course was extracted from the pre-processed rs-fMRI data. Probabilistic tractography was used to create the structural connectome matrices, and normalized by the way-total of the corresponding seed ROIs. The detailed information on the imaging and processing steps can be found in [START_REF] Korthauer | Disrupted topology of the resting state structural connectome in middle-aged apoe ε4 carriers[END_REF].

Signed resting state structural connectome

In constructing a signed resting state structural connectome, we use a novel approach introduced in [START_REF] Ajilore | Constructing the resting state structural connectome[END_REF] and has already been used in several studies (see e.g. [START_REF] Fortel | Brain dynamics through the lens of statistical mechanics by unifying structure and function[END_REF][START_REF] Fortel | Connectome Signatures of Hyperexcitation in Cognitively Intact Middle-Aged Female APOE-ε4 Carriers[END_REF][START_REF] Fortel | Inferring excitation-inhibition dynamics using a maximum entropy model unifying brain structure and function[END_REF][START_REF] Fortel | Disrupted Excitation-Inhibition balance in cognitively normal individuals at risk of alzheimer's disease[END_REF] which takes into account both structural connectivity and functional time series to form a signed coupling interaction network or "signed resting state structural connectome" (signed rsSC) to describe neural excitation and inhibition. To this end, an energy representation of neural activity based on the Ising model from statistical mechanics which ultimately bypasses traditional BOLD correlations. The spin model is a function of a coupling interaction (with positive or negative values) and spin-states of paired brain regions. Observed functional time series represent brain states over time. A maximum pseudolikelihood with a constraint is used to estimate the coupling interaction. The constraint is introduced as a penalty function such that the learned interactions are scaled relative to structural connectivity; the sign of the interactions may infer inhibition or excitation over an underlying structure. The efficiency of this approach was validated in comparing a group of healthy APOE-ε4 carriers (associated with genetic risk factor for Alzheimer's disease with a control (healthy) group of non APOE-ε4 subjects.

Here, we briefly describe the computational aspect of this approach. First, we adapted the Ising model, a well-known spin-glass model from statistical physics in which the states, also referred to as "spin configurations", of interacting units -in our case brain regions connected by white matter edges -are constrained to be either 1 ("active") or -1 ("inactive"). As described in [START_REF] Fortel | Inferring excitation-inhibition dynamics using a maximum entropy model unifying brain structure and function[END_REF][START_REF] Tang | Commpool: An interpretable graph pooling framework for hierarchical graph representation learning[END_REF] we construct a function-by-structure embedding (FSE) using a constrained pseudolikelihood estimation technique wherein pairwise interaction coefficients (represented as (J i,j ), with i and j representing ROIs in the brain network) are inferred from the observed data (BOLD time series). As the model assumes binary data, we binarize the resting-state fMRI signals. The binarized activity pattern of all ROIs at time t

(t = 1, 2, . . . , t max ) is (s(t) = s 1 (t), s 2 (t), . . . s N (t) ∈ {-1, +1} N ).
The time series data was thresholded at zero (after global signal regression), yielding a binarized sequence of -1 or +1. Note that our procedure follows the same procedure as previously proposed by other groups in this general research area that leverages the Ising model (see [START_REF] Fortel | Inferring excitation-inhibition dynamics using a maximum entropy model unifying brain structure and function[END_REF][START_REF] Fortel | Disrupted Excitation-Inhibition balance in cognitively normal individuals at risk of alzheimer's disease[END_REF] for more details).

Note that t max is determined as a result of the fMRI scan time. 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. 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:

H(s) = - i<j J i,j s i s j , (6)
such that |J i,j | ∝ W i,j , where (W i,j ) is the structural connectivity between pairs of ROIs, and the external force or bias terms are dropped in the case of resting-state. This ensures that in the pseudolikelihood estimation of (J), we constrain it with the structural connectivity (under the assumption that structural connectivity informs spin models governing brain dynamics). Thus, the optimal interaction matrix (J) is derived by maximizing the pseudo-likelihood function as:

L pseudo (J,β) = tmax t=1 k i=1
Pr s i (t) J,β, s -i(t) .

Pseudolikelihood substitutes Pr(s) by the product of the conditional probabilities p = 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:

ℓ(J,β) = 1 t max lnL pseudo (J,β) - λ 2 i<j (J i,j -sgn(J i,j )W i,j ) 2 . ( 8
)
And the pseudolikelihood component expands as follows:

1 t max ln L pseudo (J,β) = = 1 t max tmax t=1 N i=1 ln( exp(β N k=1 J i,k s i (t)s k (t)) exp(β N k=1 J i,k s k (t)) + exp(-β N k=1 J i,k s k (t))
).

(9)

Our formulation here is based on the Boltzmann distribution under pseudolikelihood conditions. Thus, the numerator describes the energy of the system, while the denominator is the sum of all possible energies. Hence, there are only two terms in the denominator since s i (t) is binary (one positive, and one negative).

The likelihood function may be simplified by setting C i (t) = β k m=1 J i,m s m (t), resulting in:

ℓ(J,β) = 1 t max tmax t=1 N i=1 C i (t)s i (t) -ln(exp(C i (t)) + exp(-C i (t)))- - λ 2 i<j (J i,j -sgn(J i,j W i,j )) 2 . ( 10
)
Here we may construct the gradient ascent procedure with respect to J i,j by computing the partial derivative of the log-pseudolikelihood as:

∂ℓ ∂J i,j = 1 t max tmax t=1 β{s i (t)s j (t) -s j (t) tanh(C i (t))} -λ(J i,j -sgn(J i,j )W i,j ). ( 11
)
The updating scheme follows:

J n+1 i,j = J n i,j + γ ∂ℓ ∂J i,j n . ( 12
)
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 procedure is followed for all subjects in constructing an optimized J matrix per subject, which we term the resting-state structural connectome or rsSC.

In Figure 5 we show the respective BOLD time series, In more detail, in Figure 5(A) we show the time evolution of the empirical BOLD signal.

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

The inferred matrix rsSC with c ij entries from the inverse Ising model (via a pseudo-likelihood maximization procedure that is further constrained by the structural connectivity, see [START_REF] Fortel | Inferring excitation-inhibition dynamics using a maximum entropy model unifying brain structure and function[END_REF]) encodes information that reveals excitation vs inhibitory relationship between brain regions. Thus, by deploying the c ij as coupling coefficients in a dynamical system (here coupled Kuramoto oscillators) we generate high-fidelity positive and negative fMRI correlations.

Our contribution is not in any way artificial, as the whole point of the paper is to outline and validate a Kuramoto-oscillator generative model with a novel brain connectome as the coupling coefficients for accurate simulations of resting-state fMRI dynamics. No prior attempt to our best knowledge has been able to match our reported accuracy in any artificial or non-artificial way. Next, we also studied what happens when some random weights of the SC are made negative, to check whether the presence of negative weights alone allows for more expressive dynamical regimes.

SC matrices

To this end, we modified the SC matrices as follows: We used the same original SC from the non-carrier subject in Figures 1-4 of the main manuscript. Then 10% of the positive connections were selected randomly and their signs were reversed to negative (all connections in the original SC matrix are nonnegative). This produced the first matrix. Similarly, we repeated the same procedure to randomly select and sign-flip an additional 20% of the non-zero positive connections to produce the second matrix. In total, this produced the 10%, 30%, 50%, 70% and 90% sign-reversed SC matrices Figures 9 and 10 shows the result for 10% and 30% respectively. Indeed, with a low fraction (10%) of sign-reversal the agreement between empirical FC and simulated FC (after fitting our proposed Kuramoto model) improved, suggesting that the presence of negative-valued coupling in the Kuramoto model is indeed necessary to sufficiently recover negative correlations in the functional connectivity. However, as expected with higher-fraction (50%, 70% and 90%, not shown here) sign-reversals of SC the agreement between empirical and simulated FC decreased again, supporting the notion of an optimal Excitation-Inhibition balance as we argue in this manuscript.

rsSC matrices

We have also performed a similar experiment by modifying the rsSC matrix of the same representative subject as follows: We first calculated the total amount of negative connections in rsSC. Then, we randomly selected 10% of these negative connections and turned them into positive by taking their absolute values. This produced the first sign-flipped rsSC matrix. We repeated the same process as before, yielding a total of 5 sign-flipped rsSC matrices at various fractions (10%, 30%, 50%, 70% and 90%) with randomly reversed signs from negative to positive. The general trend of the rsSC findings here is complementary to those in the previous SC sign-reversal experiments. Here, as the rsSC sign-reversal Taken together, given that structural connectivity is always non-negative, our simulations demonstrated that with structural connectivity as the coupling coefficients of a Kuramoto-oscillator model the recovery of negative correlations in FC is suboptimal even with the consideration of global delay. Thus, the main contribution of our manuscript is a principled biophysically-informed generative framework that leads to a signed coupling matrix (i.e. the rsSC) capable of generating substantially more accurate functional connectivity. Indeed, our model exhibits more expressive dynamical regimes and hence a better fit with the data.

Figure 1 .

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

Figure 2 .

 2 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 Eq. (1): (A) The colormap depicts the CC FC = 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 CC FC values). (B) The eFC calculated from the empirical BOLD signal. (C) The sFC matrix with the larger CC FC . Second row (using the respective subject's hybrid rsSC matrix to define the weights in model Eq. (1). (D) The respective colormap for CC FC = corr(sFC,eFC). (E) The eFC calculated from the empirical BOLD signal (same as (B)). (F) The sFC matrix with the larger CC FC . Note the ranges for the two colorbars in (A) and (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., panels (B) and (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 panels (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.

Figure 3 .

 3 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 Figures2(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 FigureS2in the Supplementary Material.

Figure 4 .

 4 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 Eq. (1): (A) The colormap depicts the CC dFC = corr(sdFC,edFC) for the parameters (K, τ ). (B) The edFC calculated from the empirical BOLD signal. (C) The sFC matrix with the larger CC dFC . Second row (using the respective subject's hybrid rsSC matrix to define the weights in the model. (D) The respective colormap for CC dFC = corr(sdFC,edFC). (E) The edFC calculated from the empirical BOLD signal (same as (B)). (E) The sdFC matrix with the larger CC dFC . Note the different ranges in the respective colorbars of panels (B), (C), (E) and (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., panels (B) and (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.10

  Figure 5(B) depicts the simulated BOLD signal with the parameters (K, τ ) ≈ (27, 20) found when optimizing CC FC using the respective rsSC matrix while in Figure 5(C) we plot the respective simulated BOLD signal with the parameters (K, τ ) ≈ (39, 18) using the SC matrix. All BOLD signals in all three panels are scaled to range in [-1, 1].

Figure 7

 7 Figure7depicts a visual comparison between the structural connectome (panel A), the c ij matrix of the resting-state structural connectome (panel B), and functional connectome (panel C) presented in the manuscript. One can clearly observe that, at least visually, the rsSC matrix (panel B) actually looks more similar to the structural connectivity (panel A) than to the function connectivity (panel C).

Figure 5 .

 5 Figure 5. Bold signals. (A) Empirical BOLD signals. (B) Simulated BOLD signals with the parameters (K, τ ) found when optimizing CC FC using the respective rsSC matrix. (C) Simulated BOLD signals with the parameters (K, τ ) found as in (B) but using the respective SC matrix. Note that all BOLD signals in all three panels are scaled to the range [-1, 1]. The small figures on top of the main panels show the respective zoomed areas for the first 20 seconds.

Figure 6 .

 6 Figure 6. 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 Figures2 (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).

Figure 7 .

 7 Figure 7. A visual comparison. Structural connectome (A), resting-state structural connectome (B) and functional connectome (C).

Figure 8 .

 8 Figure 8. Pearson correlation coefficient boxplots. Empirical (i) SC and rsSC, (ii) SC and FC and (iii) rsSC and FC matrices for all (A) non-carriers and (B) carriers datasets.

Figure 9 .

 9 Figure 9. A SC matrix with 10% random negative entries. Upper row: we selected 10% of the positive connections randomly and reversed their signs to negative. All connections in the original SC matrix are non-negative. Bottom row: parameter sweep exploration performing a similar analysis to the one presented in Figures 2 & 4 of the manuscript.

Figure 10

 10 Figure 10. A SC matrix with 30% random negative entries. Upper row: we selected 30% of the positive connections randomly and reversed their signs to negative. All connections in the original SC matrix are non-negative. Bottom row: parameter sweep exploration performing a similar analysis to the one presented in Figures 2 & 4 of the manuscript.

Figure 11

 11 Figure 11. A rsSC matrix with 10% random negative entries Upper row: After having calculated the total amount of negative connections in rsSC, we randomly selected 10% of these negative connections and turned them into positive by taking their absolute values. Bottom row: parameter sweep exploration performing a similar analysis to the one presented in Figures 2 & 4 of the manuscript.

Figure 12

 12 Figure 12. A rsSC matrix with 30% random negative entries Upper row: After having calculated the total amount of negative connections in rsSC, we randomly selected 30% of these negative connections and turned them into positive by taking their absolute values. Bottom row: parameter sweep exploration performing a similar analysis to the one presented in Figures 2 & 4 of the manuscript.

ACKNOWLEDGMENTS

We acknowledge the use of Fenix Infrastructure resources, which are partially funded from the European Union's Horizon 2020 research and innovation programme through the ICEI project under the grant agreement No. 800858. In particular, we acknowledge the access to the JUSUF supercomputer at the Jülich Supercomputer Centre. TM would like to thank the Institut Henri Poincaré for its support through the 'Research in Paris' programme and hospitality in the period during which part of this work took place.

FUNDING

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

DATA AVAILABILITY STATEMENT

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author. The TVB code used for the simulation of the brain network dynamics can be found in this repository: https://gitlab.jsc.fz-juelich.de/metaopt/rssc_ simulations_tvb.git. The scripts for generating rsSC are available here: https://github. com/iforte2/hybrid-connectome. Regarding the statistical tests used in the manuscript, standard python scripts were used (written by the authors).

CONFLICT OF INTEREST STATEMENT

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 rest of the authors declare that they have no conflicts of interest.

AUTHOR CONTRIBUTIONS

Thanos Manos: Conceptualization, Methodology, Software, Resources, Formal analysis, Writingoriginal draft, review and editing, Visualization, Supervision, Project administration, Funding acquisition. Sandra Diaz-Pier: Conceptualization, Methodology, Software, Resources, Writing review and editing. Igor Fortel: Software, Data curation, Writing review and editing. Ira Driscoll: Software, Formal analysis, Data curation. Liang Zhan: Software, Data curation, Writing review and editing. Alex Leow: Conceptualization, Methodology, Software, Formal analysis, Writing review and editing, Supervision, Funding acquisition.

ETHICS STATEMENT

This study represents a secondary data analysis using imaging data obtained as part of a larger study. See [START_REF] Korthauer | Disrupted topology of the resting state structural connectome in middle-aged apoe ε4 carriers[END_REF] for additional information on the study and the relevant institutional approval.