Equalization of the network sensitivity enhances predictive power at criticality in a wholebrain model
Abstract
Understanding the relationship between structural and functional connectivity remains a crucial issue in modern neuroscience. Here we investigate such relationship through a stochastic neural wholebrain model, using as input the structural connectivity matrix defined by the human connectome. Importantly, we address the structurefunction relationship both at the group and individualbased levels. In this work, we show that equalization of the nodes excitatory input improves the correlation between simulated neural patterns of the model and various brain functional data. We find that the best model performance is achieved when the model control parameter is in its critical value and that equalization minimizes both the variability of the critical points and neuronal activity patterns among subjects. Our results open new perspectives in personalized brain modelling with potential applications to investigate the deviation from criticality due to structural lesions (e.g. stroke) or brain disorders.
Introduction
The human brain constitutes an impressively complex system characterized by many spatiotemporal scales. At the largescale, white matter pathways derived from diffusion tensor or diffusion spectrum imaging (DTI/DSI) defines the socalled human connectome [1], i.e., a structural network of hardwired interconnections among mesoscopic brain regions. On the other hand, largescale brain activity can be accessed, among other techniques, through functional magnetic resonance imaging [2] (fMRI) which is a fourdimensional and noninvasive imaging technique that measures changes in the blood oxygen level dependent (BOLD) over time. BOLD timeseries behave as spontaneous lowfrequency ( Hz) fluctuations that have been shown to be highly correlated, across different brain areas, at rest or during a given cognitive task [3, 4, 5, 6]. Structural connections are the main input of whole brain models [7, 8] aimed at understanding how these functional patterns of correlated region activities, also called functional connectivity (FC), emerge. Understanding the relationship between structural and functional connectivity remains a crucial issue in modern neuroscience, and many studies focus on developing methods to increase the similarity between simulated and empirical functional activity, using the connectome as input [9, 10, 11, 12, 13].
From a theoretical point of view, statistical physics has decisively contributed to highlight the potential advantage that a brain may have in a critical state and also provide a quantitative description of brain activities though minimalist mesoscopic models [14, 15]. Systems consisting of many microscopic components (e.g. neurons) may exhibit rather diverse types of macroscopic collective behavior with different levels of internal organization (e.g. brain activity). Moreover, slight changes in external stimuli (e.g. auditory, visual , etc.) or in the strength of interactions themselves may induce dramatic structural rearrangements, i.e. phase transitions. It is thus tempting to hypothesize that biological states might be manifestations of similar collective phases and that shifts between them could correspond to phase transitions.
The emerging and fascinating hypothesis is that living systems, or parts of them, like the brain, are spontaneously driven close to a critical phase transition [16, 17, 18], thus conferring upon them the emergent features of critical systems like the lack of spatial and temporal scales and the high responsiveness to external perturbations. These characteristics would translate into the ability of the brain, through a large spatial and temporal scale activity, to promptly react to external stimuli by generating a coordinated global behavior [19], to maximize information transmission [20, 21], sensitivity to sensory stimuli [22] and storage of information [23].
These ideas have been particularly investigated in the last fifteen years in neuroscience and the hypothesis that the brain is poised near a critical state (in statistical mechanics sensu) is gaining consensus in the neuroscience community [14, 24, 25, 18]. In brain systems, the concept of criticality is mainly supported by the following two experimental findings: (i) the discovery of scalefree neural avalanches [20], as described by powerlaw distributions for the size and duration of the spontaneous bursts of activity in the cortex; (ii) the presence of longrange temporal correlations in the amplitude fluctuations of neural oscillations [27, 26]. Further studies reported the universality of the powerlaw exponents originally found in [20] among different species, for instance, rat [28]; nonhuman primate [29, 30] and humans using diverse techniques, such as MEG [31, 32, 33]; EEG [34] and fMRI [35, 15].
Also from a theoretical point of view, many whole brain models maximally describe realneuronal activities when they are poised at a critical point [20, 36, 15, 37, 38, 35]. Recently, a wholebrain mesoscopic model (which we call here HTC model), proposed by Haimovici et. al. [15], which is a variant of the GreenbergHastings cellular automata [39], predicts a phase transition between the subcritical regime with low activity, and the supercritical regime of high activations. When poised at the critical point, the HTC model [15] is able to capture, at the group level, the emergence of experimental spatiotemporal patterns, such as the functional connectivity (FC), the brain organization into resting state networks (RSNs) and the scaling law of the correlation length, among others. Typically these studies have been used to investigate healthy brain activity at the group level (using a single averaged functional and structural matrix from a cohort of healthy subjects) while little attention has been given to unhealthy brains. In particular, personalized brain modelling (which uses single individual DTI and fMRI as model input) has been largely unexplored for both healthy and unhealthy brains.
Recent experimental findings suggest that brain diseases (e.g., injuries, disorders) could promote a departure from the critical regime, as reported in studies of anesthesia [40], slow wave sleep [41] and epilepsy [42], where fundamentally deviations from healthy conditions promote a loss in longrange correlations and powerlaw distributions for the (spatiotemporal) neural avalanches. From a theoretical point of view, a recent work by Haimovici et. al [43] has quantified the way synthetic lesions (damages) may impact the largescale dynamical signatures of the (HTC) critical dynamics at a group level. Synthetic lesions are able to push the system out of the critical state towards a subcritical state, which is characterized by decreased levels of neural fluctuations. This result is in agreement with the experimental findings mentioned above. Indeed, Haimovici’s et. al. work is a first contribution in the direction of quantifying the effects of anatomical lesions in a whole brain model based on a critical neuronal dynamics. However, such an approach only holds for averaged groups, because the critical point in the HTC model is subject dependent resulting in a very high intersubject variability, thus inhibiting the possibility to distinguish markers of neural activity and functional patterns between healthy and injured brains. These fundamental issues have been largely unexplored so far.
In this work, we investigate the implications of the connectome topology and its degree of heterogeneity [44, 45] on the variability and the predictive power of some spatiotemporal patterns of neural [15] activity other than the neural avalanches as often investigated [46]. We address these issues for healthy brains both at the group and individual levels. In particular, we present a variant of the HTC model, where we introduce a suitable normalization of the structural connectivity matrix promoting an equalization of the nodes excitatory input, i.e. it maintains the original topology while rescaling the weights of existing connections. This equalization improves the simulated neural patterns of the model with respect to various brain functional data, such as the functional connectivity (FC), resting state networks (RSNs), and the powerlaw distribution for the sizes of active clusters in the cortex. An important result of the proposed framework is that we are able to reduce the intersubject variability within a given class of brain (e.g. healthy). In particular, we show that network equalization collapses the model state variables, i.e. neural activity patterns, of healthy subjects into universal curves, opening up a potential application on personalized brain modelling.
Theoretical Framework
The HTC model [15] consists of a discrete threestate cellular automaton in a network of nodes (i.e. cortical brain regions) linked with symmetric and weighted connections obtained from DTI/DSI scans of the white matter fiber tracts [47] and described by a matrix . The diagonal elements of (i.e., selfconnections) are all set to zero. At any given time step, each node can be in one of the three possible states: active (), inactive (), and refractory (). The state variable of a given node , , is set to if the node is active and otherwise. The temporal activity of the th node is governed by the following transition probabilities between pair of states: (i) either with a fixed probability or with probability if the sum of the connections weights of the active neighbors , , is greater than a given threshold , i.e., , otherwise , (ii) with probability , and (iii) with a fixed probability . The state of each node is overwritten only after the whole network is updated. The two parameters and controls the time scale of selfactivation and recovery of the excited state, while sets the rate of induced activity due to active nearest neighbors [15].
A mean field version of the dynamics is easily obtained in terms of the probability of node to be active, , or quiescent, , or refractory, (not to be confused with and which are the model parameters):
(1)  
(2) 
where is the Heaviside unit step function. In Eq. (1) we assumed the neighbors of node being excited as independent events. As discussed earlier [48], this approximation yields good results even when the network has a nonnegligible amount of short loops, which is the case of DTI/DSI networks considered in this study. Analytical solutions for and are difficult to be obtained. However, under suitable considerations one can obtain an approximate solution for the critical point (see Methods section), that explains its high variability within subjects. In addition, as we show in the Methods section, Eqs. (1)(2) correctly predict a collapse of across subjects when the normalized version of the input matrix is used to simulate the dynamics.
The complex behavior of the functional activities of the human brain is thought to emerge by the underlying architecture of the anatomical brain connections, as given by the binary adjacency matrix of the human connectome. It has been suggested that normalization (or equalization) of the network input sensitivity could be one of the basic ingredient [46, 49] that might lead to some simplification. Previous results [46] have shown that simulated mesoscopic neuronal dynamics is dominated by the central nodes, i.e. hubs with high indegree strength . Following [46, 50] we here propose a variant of the HTC model, by normalizing locally each entry of the structural matrix according to the following normalization rule:
(3) 
In simple terms, the above normalization fixes the indegree of all nodes to . As explained in the original reference, this normalization (or equalization) ensure that each node has at the mesoscopic level, a similar contribution on regulating the simulated brain activity. In the numerical simulation we have discretized time in steps . We set the total simulation timesteps , so to recover the length of typical (fMRI) BOLD experimental timeseries (about 520 minutes). In order to characterize simulated brain activity we have analyzed some standard quantities (see Methods section): the mean network activity (), the standard deviation of ( and the sizes of the averaged clusters, the largest and the second largest .
The simulated dynamics displays a phase transition as varies while keeping and fixed. For small values of the activation threshold , the activity is overresponsive and the signal from an active node will spread all over its first neighbors. We refer to this phase as a supercritical phase, which is characterized by sustained spontaneous activity with fast and temporally uncorrelated fluctuations. On the other hand, high values of leads to a subcritical phase, which is characterized by regular, short propagating and not selfsustained brain activity. In this phase, only those nodes with the strongest connections will determine the signal flow in the network. In between of these two phases a phase transition occurs at where brain activities have oscillatory behaviors, and longrange temporal correlations in their envelope [15, 43]. As shown in [15] the size of the second largest cluster is a suitable quantity to characterize the phase transition and it happens at the corresponding value of where is maximal. In addition to the second cluster size, the peak in the standard deviation, , may also be used to infer the critical transition [43].
To address the effects of network sensitivity in wholebrains, we thus performed our analysis using as input both structural matrices, and in its normalized counterpart . We here show that our approach is able to capture, at the critical point, the emergence of functional connectivity at rest, resting state networks (RSNs), among others. In particular, we find that the HTC model performs better when the normalization is considered and, in this case, it can also be successfully applied to individual personalized brain analysis.
Results
Group brain modelling: using average connectome as model input
Hagmann et al. dataset. We first compare the output of the presented wholebrain model on a lowresolution structural network with cortical regions obtained as an average connectome of 5 individuals [51]. The advantage of working with this dataset is that we have both average structural and functional networks, and we have a reference template for the resting state networks (see Methods). On the other hand, the DTI/fMRI matrix for each single individuals is missing, and therefore, in this case, we cannot perform an individual brain modelling.
We fixed the model parameters to the following values: , , timesteps with time discretized in seconds. We arbitrarily chose this parameterization in order to keep the ratio similar to that used originally in [15]. We then computed , , and , as a function of the threshold in the interval . We used black (red) color to represent the input matrix ().
In Fig. 1 we show (solid lines) and (dots) as a function of the rescaled threshold , where corresponds to the maximum of . Interestingly, the major global effect of normalization (for a fixed ) is not to increase the mean network activity , which in turn remains almost unchanged for small values of , but to distribute the activity more evenly across the network. Accordingly, we observe an overall increase in the strength of spatiotemporal neural variability as revealed by the peaks of and . In fact, both peaks become more pronounced after the network normalization. Finite size systems show a smooth behavior in correspondence of a phase transition in the infinite size system. Thus the divergence of a susceptibility at a critical point of an infinite system becomes a smooth peak for the corresponding finite system [52]. These peaks of maximum variability happen for different values of , but this effect is due to finite size effects (small ) affecting the position of the critical point in the system. The critical thresholds for the two quantities are expected to converge to the same as and this expectation is confirmed when using large enough (see the next section).
Another signature of criticality in the system is the cluster sizes distribution, . As shown in [35, 15] the brain forms activity clusters whose sizes follows a truncated powerlaw distribution, i.e. , with the exponential cutoff due to finite size () and a powerlaw slope of consistent with the hallmark exponent of neuronal avalanches [20]. Notice that should diverge at when . In order to identify the critical point of the system for both and in the case of small , we compute the distribution of cluster sizes for some values of , including those corresponding to the peaks of and (see insets of Fig. 1). The cluster size distributions are computed from single simulations lasting time steps. Indeed, at the peak of the second largest cluster, we find a truncated power lawdistribution, with an exponent consistent with and a cutoff which depends on the mean activity level [35, 15]. Furthermore, the normalized dynamics with , shows a better powerlaw distribution, extending up to cluster sizes , than the nonnormalized dynamics, , where it extends up to . On the other hand, scaling (if any) is less visible for corresponding to the peak of (see inset of Fig 1 (a)). Therefore from now on, we will define the critical point at the maximum of the second largest cluster. We finally note that the average size of the largest cluster (Fig. 1 (b)  solid lines) is almost indistinguishable from the average activity (Fig. 1 (a)  solid lines) in the full range of thresholds considered. This result shows that most of the time the active regions form patterns organized in a single giant component. Indeed, the largest cluster is almost two orders of magnitude larger than the second largest cluster.
We now investigate the consequences of these dynamical features on the simulated functional connectivity matrices. We employ the Pearson correlation and the distance to quantify the quality of our simulated averaged matrices (see Methods). The first index simply gives a linear correlation between the matrix elements, while the second one measures the distance between the two probability distribution functions. As usually done, we transform the model and empirical functional matrices (setting all diagonal elements to zero) in vectors, and respectively, and the Pearson correlation between both vectors, , is computed. The chisquared distance is then calculated from the corresponding (normalized) probability distribution functions and (see Methods). In Fig. 2 (a) we plot as a function of for both and its normalized counterpart . The equalization of the excitatory input leads to drastic effects on the simulated functional matrices. Indeed, it enhances the correlation with the empirical data by a factor of at with respect to the nonnormalized dynamics (see Figure 2). We find that in both cases the best model performances and occur near the critical point . How already stressed, due to finite size effects, for this database we have deviations from the thermodynamics () critical point. It is interesting to compare the performance of our critical wholebrain model with a previous work by Deco et. al. [54], where a mean field approach has been employed to study the emergence of functional connectivities. Using the same structural input (averaged DSI matrix), but different functional data, we obtain (at using the balanced matrix ) against of Deco et. al. best matching (see Fig. 3 in Ref. [54]).
For the purpose of visual inspection, we show in Figs. 2 (d)(f) the empirical connectivity matrix and the simulated ones at the corresponding minimum of , that is, and for the non normalized and normalized networks, respectively. Notice that for these particular values of both distributions have approximately . The connectivity patterns observed in the empirical matrix are best captured at criticality when using the normalized matrix , further suggesting that the normalized model is indeed better than its nonnormalized counterpart.
To gain a deep understanding of the effects of the normalization of the network nodes sensitivity we now analyze the brain organization into resting state networks (RSNs). These are a set of areas in the resting brain, i.e., when the brain is not performing any specific cognitive, language, or motor tasks, displaying BOLD fluctuations that are correlated and synchronous within the same network [6]. It has been found that RSNs are closely related to brain activation patterns seen during a given task execution, for instance, sensory (visual, auditory), cognitive, and motor etc. These spatiotemporal patterns can be obtained through spatial independent component analysis (sICA), that is the common statistical tool employed to extract RSNs from the BOLD activity (see Methods section).
In Fig. 2 (b) we show the overall match between simulated RSNs using spatial ICA and a template of wellestablished RSNs (taken from [53]) computed for the nonnormalized () and the normalized () networks as a function of . We use the Cohen’s Kappa index as a measure of similarity (see Methods). We again observe that at the normalization performed on the structural connectivity matrix significantly increases the overall match with the empirical resting state networks. This result is consistent with the previous result shown in Fig. 2 (a). In addition, some simulated RSNs are more affected by the normalization, thus showing increased similarity, as is the case of the default mode, the frontoparietal, the somatosensory and the visual networks (see the fourth column in Fig. 3). In contrast, the auditory, the cinguloopercular and the “other” networks are not significantly affected by the normalization and both input structural matrices present similar values. In Fig. 3 we also show a three dimensional projection of the templates and simulated RSNs at . It is encouraging that a quite reasonable quality of the simulated resting state networks are obtainable using a relatively lowresolution network.
Rudie et al. Dataset. We extend now our analysis to other largescale (openaccess) structural and functional dataset of Rudie et. al [55]. Different from the previous case, the actual dataset contains a large number () of individual DTI/fMRI matrices, all of them parcellated in largescale regions. By employing this dataset we can further investigate the issue of finite size effects stressed in the previous section and, most importantly, we can quantify the intersubject variability presented by the neural patterns of brain activity. The only limitation is that we do not have a reference template to the RSNs, and therefore we have not considered them in our analysis. Here we present the results for the group level DTI/fMRI matrices (averaged over 43 healthy individuals). We let the analysis of personalized brain modelling to the next section.
We fixed the two input parameters as in the previous dataset, i.e., and . The total simulation time was set to timesteps, with time discretized in seconds. Our results are shown in Fig. 4 (we maintain all the previous conventions about lines and colors as in Figs. 12). Although the network size is times larger than the one in the previous dataset, we still find finite size effects. Nevertheless the critical point for and are now closer. The nonnormalized system presents a rather broad peak of at , while in the normalized case the peak is much sharper. Similar results hold also for the mean activity and its standard deviation with peaks located at (not shown).
As in the previous section, the distribution of cluster sizes for the normalized input matrix displays a truncated powerlaw behavior in proximity and at the critical point with an exponent given by (see inset of Fig. 4 (a)).
We finally investigate whether the equalization of the excitatory input increases the correlation between simulated and empirical functional connectivity matrices. Our results are shown in Fig. 4 (b). We find again that model performance is maximized near the critical point (due to finite size effect the maximum is not exactly at ) and the normalization of the network weights caused a substantial improvement of the simulated functional connectivity matrices. At the critical point, the Pearson correlation increases by a factor of with respect to the nonnormalized input matrix. Regarding the chisquared distance, eq.(10), both systems present almost the same performance at the critical point ().
Personalized brain modelling.
By exploiting all the information in the Rudie et. al dataset, we can quantify the variability of the critical points and neural activity patterns for different individuals and their dependence on the topological properties of the underlying individual connectomes. In order to address the above issues, we have simulated the stochastic dynamics for each individual in the dataset ( healthy subjects) and calculated both the mean, the standard deviations and average clusters size of activity patterns. Simulations have been performed using the same model parameters of the group case.
We first analyzed the nonnormalized networks. Fig. 5 (a)(b) shows the behavior of , , and vs for each of the participants. The heavy lines correspond to the average curve, e.g. , where is the total number of participants. Consistent with the previous results, each of the four quantities displays a smooth behavior as a function of . Furthermore, in the nonnormalized case, each individual has its own critical threshold, according to the mean field prediction (see Methods). At the same time, consistent with the theory, the ratio between the critical threshold and the average strength for each individual connectome, does not change among individuals, i.e., , with a constant for . In the inset of Fig. 5 (b) we show the ratio for each of the participants. Except for few individuals, almost all points are peaked around , i.e. the dependence of on is correctly captured by our meanfield approach.
Next we considered the normalized networks (see Figs. 5 (d)(e)). The first striking result is the almost perfect collapse of , and (the latter shows a small variability in the peak heights) of all analyzed subjects. Regarding the second largest cluster, although the individual curves do not perfectly collapse into each other, their peaks are sharply distributed around . Thus importantly, as predicted by our meanfield approach (see Methods), for the normalized connectivity structure, the critical point of the dynamics is the same for all individuals in the dataset.
In order to test if the peak in the curve for each individual is a marker of criticality, we have analyzed the distribution of cluster sizes at the corresponding critical points, using as input both matrices and . Figures 5 (c) and 5 (f) exemplify our results for some representative individuals (results do not change for other individuals). We obtain truncated power law distributions with a relatively small cutoff () with an exponent for the nonnormalized matrix, , whereas for the normalized networks, , the power law extends about three times more (cutoff ) with an exponent , as already found for the group level.
Finally, we have analyzed the performance, at the individual level, of the normalized HTC model with respect to the functional connectivity matrix (FC). Individual FC matrices are computed at the corresponding critical points, both for and , for some representative individuals (). In all analyzed networks (results do not change for other individuals), the normalized HTC model shows a better performance. Indeed, we find an increased correlation between each simulated and empirical FCs by almost a factor of 1.5 as compared to the original HTC model. The mean of the two distributions is and (value) for the nonnormalized and normalized networks respectively. Taken together these results suggest that the normalized model is indeed better than its nonnormalized counterpart.
Discussion
One of the most compelling challenges in theoretical neuroscience is personalized brain modelling [56, 57, 58]. Indeed, when the complexity of the brain dynamics is investigated at the group level, it may hinder the interpretation of the diagnostic data leading to incomplete information for the individual subject. A key role in this challenge is played by whole brain models, which are grounded on the study of the human brain as a dynamical, complex and selforganized networked structure. Indeed, ideally, the brain activity derived from whole brain models should predict functional recovery in patients who have suffered brain damage (e.g. due to stroke). However, this attempt is strongly limited by the fact that the strength of the (cor)relations between model activity and data at the level of the individual subject [59] is very modest and the predictions can be very inaccurate. For this reason, most often whole brain models have their parameters tuned so at to replicate certain functional indicators at the population level, e.g. functional connectivity computed using the correlation between observed time series [60] or spatiotemporal patterns of local synchronization [61]. In particular, most of the models have developed indexes that are able to distinguish between different groups of subjects (e.g. healthy vs. stroked) using as input an average connectome obtained from many individuals and comparing the model output with average group properties [43, 62]. However, a prerequisite for theoretical models to be significant in terms of translational neuroscience, and thus possibly informative for therapeutic intervention, is to provide individuallevel markers (or “brain signatures”) that reliably predict cognitive and behavioral performance not only at the group level but also be adapted and tailored to the specific patient. Different markers have been recently proposed, for instance, information capacity [63] (an information theoretical measure that cannot be obtained empirically), integration [63] (a graph theoretical measure obtained from functional connectivity) and entropy [64] (a measure of repertoire diversity). All such quantities showed decreased values in stroke patients [63, 64], while information capacity and integration were additionally correlated with measures of behavioral impairment [63]. The relation between these measures and criticality will be tackled in a future work.
In this paper, we have shown that by introducing a network normalization on the HTC model, we are able not only to generate realistic brain dynamics at the group level but also to provide individualbased markers that reliably predict neuronal state activity, i.e. the criticality of the brain. In each of the analyzed connectomes of the Rudies et al. dataset, we, in fact, have found that the critical point of the generalized HTC model is located at the same . Moreover, both and have a maximum at this parameter value, and the system displays longrange correlations. A recent paper by Haimovici et. al. [43] investigated the effect of artificial lesions on the signatures of criticality of the network dynamics. Lesions were simulated by removing nodes/links of an averaged group level structural empirical matrix, targeting nodes/links according to a given network centrality and also in a random way. They found that for a stroked simulated brain there was a shift of the critical point with respect the healthy reference point, and other theoretical studies reported a loss in longrange correlations in the neural activity during anesthesia [40], slow wave sleep [41] and epilepsy [42]. However, such an approach would not be applicable in models where the critical point is individual dependent. Our results for the normalized HTC model show that introducing an equalization of the excitatory input in the simulated brain dynamics minimize the variability of the neural activity patterns and the critical point of the HTC model for different (healthy) subjects, allowing the opportunity of statistical comparison among model outputs for single individuals.
In summary, network equalization is useful in increasing the spatiotemporal variability of the brain dynamics at the individual level, which in turn increases the correlation between models outputs and empirical data. When applied at individual connectomes, the model collapses the state variables of healthy subjects into universal curves. A natural follow up of this work will be to develop an individuallevel marker based on criticality (calculated as for example) that reliably predict cognitive and behavioral performance (like in [63, 64]) as well as its evolution following therapeutic intervention. Therefore our results open new perspectives in modelling personalized brain patients with particular applications to investigate the deviation from criticality due to structural lesions (like stroke) or brain disorders.
Methods
Empirical datasets of structural connectivity and functional networks
We analyzed two different datasets of healthy subjects, consisting on both functional (fMRI) and structural (DTI/DSI) data. Specific details on the data acquisition and preprocessing can be found in the original studies.
The Hagmann et. al dataset consists of a group level (DSI) structural matrix averaged over five healthy subjects [51], and parcellated in cortical regions. The entries of the connectivity matrix represents the number of connecting fibers between a given pair of regions divided by the average area and by the average fiber length between the two regions. Functional data corresponds to BOLD timeseries measured from a cohort of other healthy subjects taken from Corbetta et. al [44]. Each subject performed two scanning runs of 10 minutes at rest. We used the Pearson correlation, eq. (9), to compute the FC matrix of each subject/scan. Then we averaged all FC matrices to obtain the group level FC.
The Rudie et. al dataset [55] consists of structural (DTI) and resting state functional matrices (obtained with BOLD fMRI) parcellated in cortical regions from a cohort of healthy typically developing individuals (13.1 2.4 years). In this case, the entries of the connectivity matrix represents the total number of fibers connecting a given pair of regions. To obtain the group level structural (functional) matrix we computed the average over the entire group of participants.
During our numerical experiments, we have analyzed the two datasets in different ways. In the first two parts of our analysis, we have investigated the effects of network equalization on the neural patterns at the group level, and thus we have simulated dynamics with an averaged structural network. Then we have compared the model output with the corresponding empirical group level FC. Finally, in the third part, we have investigated the intersubject variability simulating the neural dynamics on each individual structural matrix of the Rudie et. al dataset, then showing the feasibility of personalized modelling using our approach.
Characterization of simulated brain activity
In order to characterize the simulated brain activity through the generalized HTC model as a function of the control parameter , we have considered the following standard quantities (for simplicity we will refer to them as state variables):

the mean network activity,
(4) 
the standard deviation of ,
(5) where , is the total number of nodes and is the simulated total time;

the sizes of the averaged clusters, the largest and the second largest . Clusters were defined as ensembles of nodes that are structurally connected to each other and simultaneously active.
During simulations we kept fixed the model parameters of , and . Then we updated the network states, starting from random configurations of and states, for a total of timesteps. For each value of the threshold we computed the state variables, , , and . Throughout this study, unless stated otherwise, the final numerical results presented were averages over initial random configurations.
MeanField prediction of the critical point
Analytical solutions for and of Eqs. (1)(2) are very difficult to be obtained. However, by studying the steadstate solution () of the meanfield approximation we are able to explain the intersubject variability of the critical points. Indeed in the stationary state by setting and in Eqs. (1)(2) and approximating , where is the average network strength, one obtains, after straightforward manipulations,
(6) 
as the critical point. One finds that when and when . Notice that as one would expect, i.e. the activity is large when the threshold is low. This expression is only an approximation of the exact critical threshold. However, the meanfield solution accounts, with a reasonable agreement, for the variability of the critical points with the network strength observed numerically. Importantly, we see that for the nonnormalized dynamics, depends on the specific individuals, as is different among different brains. On the other hand, when using normalized structural connectivity then for all individuals and therefore is universal.
Model Validation
From the Model Output to BOLD Signal. Experimentally, brain activity at rest can be accessed through fMRI. In fMRI what is measured is the variation of the bloodoxygenlevel dependent (BOLD) signal. Moreover, following [15] we simulate BOLD timeseries of each node convolving the node variable with a canonical doublegamma hemodynamic response function (HRF),
(7) 
with,
(8) 
where is the BOLD signal of the th node. The free parameters in (8) were fixed according to values found in [65], i.e., , , , , and . Finally, the convolved timeseries, , were filtered with a zero lag finite impulse response band pass filter in the frequency range of . Although complicated, these steps are part of a standard procedure to transform model output in BOLD functional signals.
From the generated BOLD signal we can finally extract the functional connectivity (FC) networks. In fact, the FC matrix is defined through Pearson correlation:
(9) 
where is the standard deviation and is the temporal average of the BOLD time series.
To access the quality of our results we need to compare the generated FC matrix with the functional networks obtained from the fMRI data. In particular, we employed two distinct statistical measures to quantify the similarity between simulated and empirical functional matrices: (i) the Pearson correlation and, (ii) the chisquared distance (). As usually done, we transform the model and empirical functional matrices (setting all diagonal elements to zero) in vectors, and respectively, and the Pearson correlation between both vectors, , is computed. The distance is then calculated from the (normalized) probability distribution functions and ,
(10) 
where is the number of bins used to calculate both histograms.
Resting State Networks. The rest brain activity displays coherent spatiotemporal activation patterns which have been consistently found in healthy subjects [5, 6]. These spatiotemporal maps reflect regions that are functionally connected, i.e., with a similar BOLD activity, although they may be anatomically disconnected. The brain organization into resting state networks have been vastly extracted using spatial and temporal independent component analysis (sICA/tICA) [66, 67]. Here we applied the spatial ICA to extract RSNs from the BOLD activity. sICA decomposes a set of BOLD timeseries into a number of independent components (specified a priori) which are spatial maps associated with the time courses of the signal sources. In matrix notation it reads,
(11) 
where is the () raw matrix containing in its columns the simulated timeseries (of length ). Spatial maps are encoded in the rows of (of order ) and the corresponding time courses of the signal sources in (of order ).
We employed the fastICA algorithm in R (opensource platform) to estimate the independent components (ICs). After that ICs maps were ztransformed, i.e., for . For each value of we repeated such procedure times with distinct initial random configurations. At the end, we end up with a pool of ICs maps following a Gaussian distribution with zero mean and unity standard deviation. We finally threshold and binarize ICs. In particular, we set to 1 all elements such that () and zero otherwise.
We access the quality of our simulated spatial maps (ICs) by computing the Cohen’s kappa similarity index [68], , with a template of wellestablished human RSN taken from [53]. The template contains the name of the anatomical brain regions mainly involved in a given RSN network. We used it to match each node of our network belongs to a given empirical RSN. Such procedure resulted in a total of binary RSN template vectors, namely, auditory (6 matched nodes), cinguloopercular (12 nodes), default mode (8 nodes), frontoparietal (10 nodes), somatosensory (6 nodes), visual (10 nodes) and “other” (10 nodes) resting state networks. We omitted the ventral and dorsal attention templates because they comprised a small number of matched nodes (). Following [15, 69], we performed a best match approach to assign each simulated IC to be belonging to a given RSNs. Indeed, we computed between each with all template vectors, assigning the RSN with highest and averaging the corresponding values across the pool of ICs to obtain the overall average match (see Fig. 2 (b)). We also computed for each RSNs by simply averaging those ICs assigned to be closest to a given template vector (see Fig. 3).
We finally fixed the free parameters, namely, the threshold and the number of independent components , in a data drivenway. Accordingly, we applied the above framework to the Corbetta et. al. dataset [44] (48 empirical BOLD timeseries) and then fixed the parameters in such a way to maximize the overall match . By varying a twodimensional parameter space we found a maximum at and corresponding to the 92th percentile of the entire pool of ICs (see Supporting Figures, Figs. S1, S2 and S3).
Acknowledgements
R.P.R thank the University of Padova, Physics and Astronomy Department, and INFN; the financial support from the Brazilian agency CNPq (Grant No.201241/20153) and Francesco D’Angelo for useful discussions.
References
 Olaf Sporns, Giulio Tononi, Rolf Kötter. The human connectome: A structural description of the human brain. PLoS Comput Biol 1(4): e42 (2005).
 Martijn P. van den Heuvel, Hilleke E. Hulshoff Pol. Exploring the brain network: A review on restingstate fMRI functional connectivity. European Neuropsychopharmacology 20, 519–534 (2010).
 Ed Bullmore and Olaf Sporns. Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. 10, 186 (2009).
 H. J. Park and K. Friston. Structural and functional brain networks: from connections to cognition. Science 342, 1238411 (2013).
 Michael D. Fox, Abraham Z. Snyder, Justin L. Vincent, Maurizio Corbetta, David C. Van Essen, and Marcus E. Raichle. The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proc. Natl. Acad. Sci. USA. 102, 9673 (2005).
 J. S. Damoiseaux, S. A. R. B. Rombouts, F. Barkhof, P. Scheltens, C. J. Stam, S. M. Smith, and C. F. Beckmann. Consistent restingstate networks across healthy subjects. Proc. Natl. Acad. Sci. USA. 103, 13848 (2006).
 Joana Cabral, Morten L. Kringelbach, Gustavo Deco. Exploring the network dynamics underlying brain activity during rest. Progress in Neurobiology 114, 102 (2014).
 Joana Cabral, Morten L. Kringelbach, Gustavo Deco. Functional connectivity dynamically evolves on multiple timescales over a static structural connectome: Models and mechanisms. NeuroImage 160, 84 (2017).
 C. J. Honey, O. Sporns, L. Cammoun, X. Gigandet, J. P. Thiran, R. Meuli, and P. Hagmann. Predicting human restingstate functional connectivity from structural connectivity. Proc. Natl. Acad. Sci. USA. 106, 2035 (2009).
 Farras Abdelnour, Henning U. Voss, Ashish Raj. Network diffusion accurately models the relationship between structural and functional brain connectivity networks. NeuroImage 90, 335 (2014).
 Gustavo Deco, Anthony R. McIntosh, Kelly Shen, R. Matthew Hutchison, Ravi S. Menon, Stefan Everling, Patric Hagmann and Viktor K. Jirsa. Identification of Optimal Structural Connectivity Using Functional Connectivity and Neural Modeling. J. Neurosci. 34, 7910 (2014).
 Arnaud Messé, David Rudrauf, Alain Giron, Guillaume Marrelec. Predicting functional connectivity from structural connectivity via computational models using MRI: an extensive comparison study. NeuroImage 111, 65 (2015).
 Maria Luisa Saggio, Petra Ritter, Viktor K. Jirsa. Analytical Operations Relate Structural and Functional Connectivity in the Brain. PLoS ONE 11 (8): e0157292 (2016).
 D. R. Chialvo. Emergent complex neural dynamics. Nat. Phys. 6, 744 (2010).
 Ariel Haimovici, Enzo Tagliazucchi, Pablo Balenzuela, and Dante R. Chialvo. Brain Organization into Resting State Networks Emerges at Criticality on a Model of the Human Connectome. Phys. Rev. Lett. 110, 178101 (2013).
 Strictly speaking phase transitions exist only for systems with an infinite number of degrees of freedom, which at best are good approximation of large, but finite, systems like a brain.
 Jorge Hidalgo, Jacopo Grilli, Samir Suweis, Miguel A. Muñoz, Jayanth R. Banavar and Amos Maritan. Informationbased fitness and the emergence of criticality in living systems. Proc Natl Acad Sci USA 111 (28) 1009510100 (2014).
 Miguel A. Muñoz. Criticality and dynamical scaling in living systems. arXiv:1712.04499 (2017).
 Elad Schneidman, Michael J. Berry, Ronen Segev and William Bialek. Weak pairwise correlations imply strongly correlated network states in a neural population. Nature 440, 1007–1012 (2006).
 J. M. Beggs and D. Plenz. Neuronal avalanches in neocortical circuits. J. Neurosci. 23, 11167 (2003).
 Bertha VázquezRodríguez, Andrea AvenaKoenigsberger, Olaf Sporns, Alessandra Griffa, Patric Hagmann, and Hernán Larralde, Stochastic resonance at criticality in a network model of the human cortex. Sci. Reports 7,13020 (2017).
 O. Kinouchi and M. Copelli. Optimal dynamical range of excitable networks at criticality. Nature Phys. 2, 348 (2006).
 C. Haldeman and J. Beggs. Critical Branching Captures Activity in Living Neural Networks and Maximizes the Number of Metastable States. Phys. Rev. Lett. 94, 058101 (2005).
 Paolo Massobrio, Lucilla de Arcangelis, Valentina Pasquale, Henrik J. Jensen and Dietmar Plenz. Criticality as a signature of healthy neural systems. Frontiers in Systems Neuroscience 15, 22 (2015).
 J. Hesse and T. Gross. Selforganized criticality as a fundamental property of neural systems. Front. Syst. Neurosci. 8, 166 (2014).
 R. Hardstone, SimonShlomo Poil, G. Schiavone, R. Jansen, V. V. Nikulin, H. D. Mansvelder and K. LinkenkaerHansen. Detrended Fluctuation Analysis: A ScaleFree View on Neuronal Oscillations. Front. Physiol. 3, 450 (2012).
 K. LinkenkaerHansen, V. V. Nikulin, J. M. Palva and R. J. Ilmoniemi. Longrange temporal correlations and scaling behavior in human brain oscillations. J. Neurosci. 21, 1370 (2001).
 Gireesh, E. D., and Plenz, D. (2008). Neuronal avalanches organize as nested theta and beta/gamma oscillations during development of cortical layer 2/3. Proc. Natl. Acad. Sci. U.S.A. 105, 7576–7581. doi: 10.1073/pnas.0800537105
 Petermann, T., Thiagarajan, T. C., Lebedev, M. A., Nicolelis, M. A., Chialvo, D. R., and Plenz, D. (2009). Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. Proc. Natl. Acad. Sci. U.S.A. 106, 15921–15926 (2009).
 Yu, S., Yang, H., Nakahara, H., Santos, G. S., Nikolic, D., and Plenz, D. Higherorder interactions characterized in cortical activity. J. Neurosci. 31, 17514–17526 (2011).
 Poil, S.S., Hardstone, R., Mansvelder, H. D., and LinkenkaerHansen, K. Criticalstate dynamics of avalanches and oscillations jointly emerge from balanced excitation/inhibition in neuronal networks. J. Neurosci. 32, 9817–9823 (2012).
 Palva, J. M., Zhigalov, A., Hirvonen, J., Korhonen, O., LinkenkaerHansen, K., and Palva, S. Neuronal longrange temporal correlations and avalanche dynamics are correlated with behavioral scaling laws. Proc. Natl. Acad. Sci. U.S.A. 110, 3585–3590 (2013).
 Shriki, O., Alstott, J., Carver, F., Holroyd, T., Henson, R. N. A., Smith, M. L., et al. Neuronal avalanches in the resting MEG of the human brain. J. Neurosci. 33, 7079–7090. (2013).
 Meisel, C., Olbrich, E., Shriki, O., and Achermann, P. Fading signatures of critical brain dynamics during sustained wakefulness in humans. J. Neurosci. 33, 17363–17372 (2013).
 Enzo Tagliazucchi, Pablo Balenzuela, Daniel Fraiman and Dante R. Chialvo. Criticality in LargeScale Brain fMRI Dynamics Unveiled by a Novel Point Process Analysis. Frontiers in Physiology 3, 15 (2012).
 Daniel Fraiman, Pablo Balenzuela, Jennifer Foss, and Dante R. Chialvo. Isinglike dynamics in largescale functional brain networks. Phys. Rev. E 79, 061922 (2009).
 Gustavo Deco and Viktor K. Jirsa. Ongoing cortical activity at rest: criticality, multistability, and ghost attractors. J Neurosci. 32(10):33663375 (2012).
 Gustavo Deco, Anthony R. McIntosh, Kelly Shen, R. Matthew Hutchison, Ravi S. Menon, Stefan Everling, Patric Hagmann, and Viktor K. Jirsa. Identification of Optimal Structural Connectivity Using Functional Connectivity and Neural Modeling. J. Neurosci. 34(23):79107916 (2014).
 J. M. Greenberg and S. P. Hastings. Spatial Patterns for Discrete Models of Diffusion in Excitable Media. SIAM Journal on Applied Mathematics, 34 (3), 515–523 (1978).
 Scott G, Fagerholm ED, Mutoh H, Leech R, Sharp DJ, Shew WL, Knöpfel T. Voltage imaging of waking mouse cortex reveals emergence of critical neuronal dynamics. J. Neurosci. 34, 1661116620 (2014).
 Priesemann V, Valderrama M, Wibral M, Le Van Quyen M. Neuronal avalanches differ from wakefulness to deep sleep: evidence from intracranial depth recordings in humans. PLoS Comp. Biol. 9(3), e1002985 (2013).
 Meisel C, Storch A, HallmeyerElgner S, Bullmore E, Gross T. Failure of adaptive selforganized criticality during epileptic seizure attacks. PLos Comp. Biol. 8(1), e1002312 (2012).
 Ariel Haimovici, Pablo Balenzuela, and Enzo Tagliazucchi. Dynamical signatures of structural connectivity damage to a model of the brain posed at criticality. Brain Connectivity, 6(10): 759771 (2016).
 A. PonceAlvarez, Gustavo Deco, Patric Hagmann, G. Luca Romani, Dante Mantini, Maurizio Corbetta. RestingState Temporal Synchronization Networks Emerge from Connectivity Topology and Heterogeneity. PLoS Comput Biol 11(2): e1004100 (2015).
 Kaustubh Manchanda, Amitabha Bose and Ramakrishna Ramaswamy. Collective dynamics in heterogeneous networks of neuronal cellular automata. Physica A 487, 111 (2017).
 Géza Ódor. Critical dynamics on a large human Open Connectome network. Phys. Rev. E 94, 062411 (2016).
 Yaniv Assaf and Ofer Pasternak. Diffusion Tensor Imaging (DTI)based White Matter Mapping in Brain Research: A Review. J Mol Neurosci 34, 51 (2008).
 Daniel B. Larremore, Woodrow L. Shew, and Juan G. Restrepo. Predicting Criticality and Dynamic Range in Complex Networks: Effects of Topology. Phys. Rev. Lett. 106, 058101 (2011).
 Rony Azouz and Charles M. Gray. Dynamic spike threshold reveals a mechanism for synaptic coincidence detection in cortical neurons in vivo. Proc Natl Acad Sci USA 97, (14) 81108115 (2000).
 M.T. Hütt, M. K. Jain, C. C. Hilgetag, and A. Lesne, Stochastic resonance in discrete excitable dynamics on graphs, Chaos, Solitions Fractals 45, 611 (2012).
 Hagmann P, Cammoun L, Gigandet X, Meuli R, Honey CJ, et al. Mapping the structural core of human cerebral cortex. PLoS Biol 6: e159 (2008).
 Notice that a peak in a finite system does not necessarily is indicative of a transition in the corresponding infinite system. However in our case we make the assumption that the observed peaks in both and , in the normalized dynamics, correspond indeed to a critical phase transition in the infinite system.
 Ann E. Sizemore, Chad Giusti, Ari Kahn, Jean M. Vettel, Richard F. Betzel, Danielle S. Bassett. Closures and cavities in the human connectome. J. Comput. Neurosci. 44, 115145 (2018).
 Gustavo Deco, Adrián PonceAlvarez, Dante Mantini, Gian Luca Romani, Patric Hagmann, and Maurizio Corbetta. Restingstate functional connectivity emerges from structurally and dynamically shaped slow linear fluctuations. The Journal of Neuroscience 33, 11239 (2013).
 J. D. Rudie, J. A. Brown, D. BeckPancer, L. M. Hernandez, E. L. Dennis, P. M. Thompson, S. Y. Bookheimer and M. Dapretto. Altered functional and structural brain network organization in autism. Neuroimage Clin.2, 7994 (2012).
 Maria I. Falcon, Viktor Jirsa, and Ana Solodkin. A new neuroinformatics approach to personalized medicine in neurology: The Virtual Brain. Curr. Opin. Neurol. 29(4): 429–436 (2016).
 V.K. Jirsa, T. Proix, D. Perdikis, M.M. Woodman, H. Wang, J. GonzalezMartinez, C. Bernard, C. Bénar, M. Guye, P. Chauvel, F. Bartolomei. The Virtual Epileptic Patient: Individualized wholebrain models of epilepsy spread. NeuroImage 145, 377–388 (2017).
 Kanika Bansal, Johan Nakuci, Sarah Feldt Muldoon. Personalized brain network models for assessing structurefunction relationships. arXiv:1802.00473 (2018).
 Gustavo Deco and Morten L. Kringelbach. Great Expectations: Using WholeBrain Computational Connectomics for Understanding Neuropsychiatric Disorders. Neuron 84, 892905 (2014).
 Friston et al. Functional topography: multidimensional scaling and functional connectivity in the brain, Cereb Cortex 6:15664 (1996).
 Gustavo Deco et al. The role of rhythmic neural synchronization in rest and task conditions. Frontiers in human neuroscience 5: 4 (2011).
 G. Deco, G. Tononi, M. Boly, and M. L. Kringelbach. Rethinking segregation and integration: contributions of wholebrain modelling. Nature Reviews Neuroscience, 16(7), 430 (2015).
 Mohit H. Adhikari, Carl D. Hacker, Josh S. Siegel, Alessandra Griffa, Patric Hagmann, Gustavo Deco, and Maurizio Corbetta. Decreased integration and information capacity in stroke measured by whole brain models of resting state activity. Brain 140, 1068–1085 (2017).
 Victor M. Saenger, Adrián PonceAlvarez, Mohit Adhikari, Patric Hagmann, Gustavo Deco, and Maurizio Corbetta. Linking Entropy at Rest with the Underlying Structural Connectivity in the Healthy and Lesioned Brain. Cerebral Cortex, 2017, 111.
 Glover G. Deconvolution of Impulse Response in EventRelated BOLD fMRI. NeuroImage 9, 416–429 (1999).
 Calhoun VD, Adali T, Pearlson GD, Pekar JJ. A method for making group inferences from functional MRI data using independent component analysis. Hum Brain Mapp 14, 140–151 (2001).
 Mantini D, Corbetta M, Romani GL, Orban GA, Vanduffel W. Evolutionarily novel functional networks in the human brain? J Neurosci 33, 3259–3275 (2013).
 Mary L. McHugh. Interrater reliability: the kappa statistic. Biochem Med (Zagreb) 22, 276–282 (2012).
 Katharina Glomb, Adrian PonceAlvarez, Matthieu Gilson, Petra Ritter, Gustavo Deco. Resting state networks in empirical and simulated dynamic functional connectivity. NeuroImage 159, 388–402 (2017).