(C) PLOS One
This story was originally published by PLOS One and is unaltered.
 . . . . . . . . . .



Strength-dependent perturbation of whole-brain model working in different regimes reveals the role of fluctuations in brain dynamics [1]

['Yonatan Sanz Perl', 'Center For Brain', 'Cognition', 'Computational Neuroscience Group', 'Department Of Information', 'Communication Technologies', 'Universitat Pompeu Fabra', 'Barcelona', 'Department Of Physics', 'University Of Buenos Aires']

Date: 2022-11 

Despite decades of research, there is still a lack of understanding of the role and generating mechanisms of the ubiquitous fluctuations and oscillations found in recordings of brain dynamics. Here, we used whole-brain computational models capable of presenting different dynamical regimes to reproduce empirical data’s turbulence level. We showed that the model’s fluctuations regime fitted to turbulence more faithfully reproduces the empirical functional connectivity compared to oscillatory and noise regimes. By applying global and local strength-dependent perturbations and subsequently measuring the responsiveness of the model, we revealed each regime’s computational capacity demonstrating that brain dynamics is shifted towards fluctuations to provide much-needed flexibility. Importantly, fluctuation regime stimulation in a brain region within a given resting state network modulates that network, aligned with previous empirical and computational studies. Furthermore, this framework generates specific, testable empirical predictions for human stimulation studies using strength-dependent rather than constant perturbation. Overall, the whole-brain models fitted to the level of empirical turbulence together with functional connectivity unveil that the fluctuation regime best captures empirical data, and the strength-dependent perturbative framework demonstrates how this regime provides maximal flexibility to the human brain.

How and why do complex, fluctuating, and oscillating dynamics characterise brain states? We combined a whole-brain model and strength-dependent perturbation frameworks to investigate the causal mechanistic explanation behind the human brain function. We demonstrated by fitting whole-brain models to the level of empirical turbulence together with functional connectivity that the fluctuation regime best captures empirical data. Furthermore, the strength-dependent perturbative approach allows us to assess the computational capabilities of different dynamical regimes. We showed that the fluctuations regime provides maximal flexibility to the human brain, a desirable property for brain dynamics to interact with the environment.

To address the question of distinguishing between noise, fluctuation and oscillations regimes in shaping global brain dynamics, we constructed and perturbed whole-brain models taking advantage of the recently proposed turbulence framework. First, we fitted the model to the level of empirical turbulence to precisely locate the optimal parameters in each regime, showing that the noise regime is suboptimal compared with fluctuations and oscillations that are equally good. To disentangle the generative roles of the fluctuation (subcritical) and oscillations (supercritical) models, we investigated the performance of each model regime in the corresponding optimal working point to fit empirical turbulence at fitting the functional connectivity (probably one of the major outcome measures evaluated thus far in the literature [ 31 , 32 ]). We found that the fluctuations regime overperforms the oscillatory regime showing that the conjunction of turbulence and FC fitting provides the information to unveil the best model regime that captures the empirical observables. We then created a global and local strength-dependent perturbational framework to characterise the computational capabilities in each regime in terms of the responsiveness to external stimulus. To do so, we quantified the evolution of three perturbative sensitivity measures, susceptibility, information capability (defined in the turbulence framework [ 10 ]) and the well-known PCI, as a function of the applied global sustained perturbation. The susceptibility is the ability of the system to be externally perturbed, based on the work of Hiroaki Daido [ 33 ], who defines the susceptibility of a large population of coupled oscillators as the variation in the system synchronisation under external perturbation. The information capability measures the ability of the system to encode external inputs, as such is closely related to automatic complexity evaluator (ACE), and synchrony coalition entropy (SCE) (used and defined in [ 34 ]). We found that the fluctuation regime provides maximal flexibility to brain dynamics, a desirable property for brain dynamics to interact with the environment. The stimulation in the fluctuation regime also provides a mechanistic explanation of experimentally reported brain dynamics showing the modulation of specific resting-state networks when they are perturbed. We show that using a combined whole-brain model fitting and strength-dependent perturbations—instead of the classical flat, constant perturbations–framework provides the means for disentangling alternative model-based hypotheses about the underlying empirical dynamics and characterising their computational capabilities. Our results show that the fluctuation regime more faithfully reproduces empirical properties such as the level of turbulence and functional connectivity, allowing maximal flexibility in information processing.

In a complementary direction, empirical perturbations have proved to be an excellent approach to provide insights into the complexity of brain dynamics, such as the one proposed by Massimini and colleagues. The authors used transcranial magnetic stimulation (TMS) with electroencephalography (EEG) to demonstrate perturbation-elicited changes in global brain activity in the perturbative complexity index (PCI) between different brain states (wakefulness, sleep, anaesthesia and coma) [ 25 – 27 ]. The results showed, for example, that non-REM sleep is accompanied by a breakdown in cortical effective connectivity, where the stimuli rapidly extinguish and do not propagate beyond the stimulation site [ 25 – 27 ]. Based on the same strategy, previous experimental research demonstrated that stimulation with TMS in specific brain regions can differentially modulates specific networks [ 28 , 29 ]. Computational approaches also demonstrated that constant perturbation can also provide clear insights into the complexity of brain dynamics (Deco et al., 2018; Goldman et al., 2020; Kunze et al., 2016).

Further to these empirical observations, the relevance of whole-brain computational model within this analytical framework was demonstrated based on the fact that brain dynamics can be accurately modelled by a system of coupled non-linear Stuart-Landau working in a turbulent regime [ 8 , 10 ], which integrates anatomy and local dynamics [ 16 – 19 ]. As it happens, this whole-brain model is suited to resolve the question in hand, since it naturally describes the transitions between noise, fluctuation and oscillation. The whole-brain system can produce three radically different regimes, simply by varying the local bifurcation parameter: 1) Noise regime–when the parameter is much less than zero, 2) fluctuating subcritical regime–when the parameter is just below zero; and 3) oscillatory supercritical regime–when the parameter is larger than zero. Previous research has shown that the three different scenarios of noise [ 20 , 21 ], subcritical [ 16 , 21 – 23 ] and supercritical [ 24 ] are equally able to fit the empirical neuroimaging data in terms of functional connectivity. However, it is not clear when fitting the whole-brain model to the level of turbulence which regimes are able to fit the empirical data.

Deco & Kringelbach have proposed a novel framework based on recent results showing turbulence in the brain dynamics, which is based on quantifying the level of local synchronisation in whole-brain activity [ 8 – 10 ]. Briefly, turbulence dynamics in non-fluid systems can be defined by the coupled oscillator framework of Kuramoto [ 11 , 12 ]. In this sense, brain dynamics present turbulent behaviour regarding the amplitude variability of the local Kuramoto order parameter. Furthermore, turbulence has been shown to provide the optimal transmission of energy (which is closely related to information [ 13 ]), and at the core of this transmission are the scale-free mixing properties of turbulence. Kolmogorov’s seminal phenomenological statistical studies have shown that this transmission is highly efficient across scales within the turbulent regime [ 14 , 15 ]. In Deco & Kringelbach approach, it has also been shown that this efficient transmission, which is demonstrated in power laws relation between information and space, is also present in brain dynamics [ 8 ]. Note that for these power laws, the flow is not provided by billions of molecules in a fluid but by the flow of information coming from the interplay and mutual entrainment in billions of neurons underlying brain dynamics.

Already at the birth of neuroscience, a deep problem emerged: namely that local and global recordings from inside and outside the brain show very complex fluctuating and oscillating patterns of brain activity [ 1 – 5 ]. This gave rise to the fundamental question of the importance of synchronous or asynchronous local dynamics as the origin of the dynamical behaviour of brain states [ 6 , 7 ]. In global brain dynamics, a purely fluctuating scenario will give rises to patterns formed due to noise correlations, whereas a purely oscillatory regime would produce patterns arising mainly from cluster synchronisation. In both cases, the activity is shaped by the underlying brain anatomy but the generating principles are clearly different. Even more, the asynchronous, irregular background dynamics have been associated with conscious, responsive brain state, while synchronisation and regular dynamics have been linked with reduced states of conscious awareness [ 6 ].

Results

In this study, we were interested in revealing the underlying mechanisms of the different dynamical regimes available in the resting state [2,35]. We extended the well stablish perturbative approach to use strength-dependent, non-constant perturbation in a whole-brain model fitting the empirical data to provide a causal mechanistic explanation for disentangling fluctuating from oscillating regimes in the underlying empirical dynamics. Fig 1 shows the details of our framework, which has two key ingredients: 1) a model-based approach which is probed with 2) varying levels of strength-dependent perturbations. The whole-brain model is based on the recent of demonstrating turbulence (Fig 1A) in empirical neuroimaging data (Fig 1B). Turbulence is a property found in high-dimensional non-linear systems, where its mixing capability is crucial for giving rise to the efficient energy/information cascade, whereby large whirls turns into smaller whirls and eventually energy dissipation. Using this turbulence framework, we were able to determine the vorticity, i.e. the local level of synchronisation, allowing us to extend the standard global time-based measure of metastability to become a local-based measure of both space and time for capturing the essential features of the underlying brain dynamics.

PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 1. Overview of the framework. A) Turbulence provides a good description of the seemingly chaotic dynamics of fluids as first described by Leonardo da Vinci [9] (left panel drawing of turbulent whirls). The physical principles giving rise to turbulence are given by high-dimensional spacetime non-linear coupled systems. In turbulence, a fundamental property is its mixing capability which yields the energy cascade through turning large whirls into smaller whirls and eventually energy dissipation (middle panel). Furthermore, the turbulent energy cascade has been shown to be highly efficient across scales, as evidenced by a power law (right panel). B) Empirical brain dynamics was recently shown to exhibit turbulence [8]. The fMRI resting state analysis over 1000 healthy participants (left panel) shows the presence of highly variable, local synchronisation vortices across time and space (middle panel). Equally, the turbulent brain regime also gives rise to an efficient information cascade obeying a power law (right panel). C) Furthermore, Hopf whole-brain models [16] (left panel) were able to fit both turbulence and the empirical data at the same working point (right panel). D) We model brain activity as a system of non-linear Stuart Landau oscillators, coupled by known anatomical connectivity. E) The Stuart Landau equation (top panel) is suited for describing the transitions between noise and oscillation. By varying the local bifurcation parameter, a, the equation will produce three radically different regimes: Noise (a<<0), fluctuating subcritical regime (a<0 & a~0) and oscillatory supercritical regime (a>0) (bottom panel). F) We evaluated the fitting capacity of the three model regimes in terms of functional connectivity and turbulence (with the dashed line showing the empirical level of turbulence). G) However, it is well-known that physical systems can be more deeply probed by perturbing them. Therefore, we used strength-dependent perturbations to disentangle the generative roles of the fluctuation (subcritical) and oscillations (supercritical) models. We observed the evolution of two key perturbative measures, susceptibility and information capacity, as a function of the applied global sustained perturbation. H) Finally, in order to generate experimentally testable hypotheses, we used local strength-dependent, non-sustained perturbations and measured the elicited dynamics in terms of the empirical perturbative complexity index [25]. Specifically, we simulated 600 volumes with the perturbation active, and we then evaluated the evolution of the signals in the following 200 volumes without perturbation and computed the difference between the PCI after and before perturbation. https://doi.org/10.1371/journal.pcbi.1010662.g001

The Hopf whole-brain model can fit the complex spatiotemporal brain dynamics in terms of both functional connectivity and turbulence (Fig 1C). More generally, the Hopf whole-brain model integrates anatomical connections [36–38] with local dynamics to explain and fit the emergence of global dynamics in empirical data [16,32,39–41] (Fig 1D). For decades, brain signals have been recorded with a plethora of different techniques showing them to be combinations of at least three different regimes: noise, fluctuating, and oscillatory. The non-linear Stuart Landau oscillator is perfect for generating and testing these three regimes, given that the local bifurcation parameter in the equation governs the dynamics of each local brain region (Fig 1E). Indeed, by varying this parameter the Stuart-Landau equation will produce three radically different signals: 1) a noise signal resulting from Gaussian noise added to a fixed point when the parameter is much less than zero resulting in a pure noise signal; 2) a fluctuating stochastically structured signal when the parameter is just below zero, which allows the system to fluctuate between noise and oscillations; 3) an oscillatory signal when the parameter is larger than zero. Technically, these three regimes are termed noise, subcritical and supercritical, respectively.

In the following, we show that the subcritical fluctuating and supercritical oscillatory regimes are equally able to fit the empirical data in terms of functional connectivity and turbulence (Fig 1F). Crucially, however, our framework includes the second ingredient of strength-dependent perturbation, which, as shown below, has allowed us to distinguish between the two regimes. We probe the model in two ways using both global (Fig 1G) and local strength-dependent perturbations (Fig 1H) and measuring the sensitivity of the system through quantifying the elicited susceptibility, information capacity and perturbative complexity index.

Hopf whole-brain model of large-scale empirical neuroimaging data We first investigated the ability of the three regimes to fit empirical data. We fitted whole-brain models of Stuart Landau oscillators in the three regimes to the large-scale neuroimaging resting state fMRI data from 1003 healthy human participants in the Human Connectome Project (HCP) [42]. We extracted the timeseries in the Schaefer1000 parcellation [43], a fine-grained atlas that allowed quantified turbulence in empirical data [8]. Previous Hopf whole-brain models have successfully fitted functional neuroimaging data with different acquisition parameters from many different neuroimaging setups [32,44,45] using a fluctuating regime with a local bifurcation parameter close to the bifurcation point. Here we aim to fit both functional connectivity and turbulence of functional neuroimaging data. The functional connectivity defined as the Pearson correlation between all pairs of nodes signal and the turbulence defined as the amplitude variability of the local Kuramoto order parameters. In order to fit turbulence with the Stuart-Landau oscillator in the oscillatory supercritical regime, Kuramoto and colleagues [12] have shown that an extra parameter, the so-called shear parameter, is fundamental. This parameter is also called the “nonisochronicity parameter” and provides a different path to controls the synchronicity of coupled oscillators [13,46,47]. Specifically, this parameter can prevent the full synchronization of oscillators when are coupled [46], which is counterintuitive and necessary to fit brain dynamics. Following the analogy with fluids dynamics this parameter is related with the viscosity of the media [11–13]. Therefore, we extend the Hopf whole-brain model to use the appropriate formulation of the Stuart-Landau equation (see Methods) to be able to fit the data with the supercritical regime. We explored the parameter space of varying the global coupling (G) and the shear parameter (β). The coupling parameter (G) scales the local fibre densities of the anatomical structural connectivity (see Methods) to capture the effectivity of the coupling by assuming a single global conductivity parameter. The shear parameter (β) acts similar to viscosity in fluid dynamics [11] in that it is able to affect both the frequency and amplitude of the generated oscillations [12]. Importantly, for the structural connectivity in the whole-brain model, we used a combination of exponential distance rule, EDR [48] and long-range connections (LR), which improve the fit to the available dMRI tractography from humans [10](see Methods). In order to fit the whole-brain model, we used the following observables: 1) the empirical mean level of amplitude turbulence, as the standard deviation of the Kuramoto Local order parameter (D), in fine parcellation, and metastability, as the standard deviation of the Kuramoto Global order parameter (M), in coarse parcellation; and 2) the grand average functional connectivity (FC) from the neuroimaging empirical data (see Methods). For measuring the level of fitting for each: 1) for turbulence/metastability measure, we computed the error (eD = abs(D sim -D emp )/(eM = abs(M sim -M emp )), i.e. by the absolute difference between the simulated and empirical amplitude turbulence and 2) for the functional connectivity, we computed Euclidean distance (eFC) between the simulated and empirical FC.

Modelling results for fine-scale parcellation with 1000 regions Fig 2 shows the results of fitting the Hopf whole-brain model in the three different regimes (noise, fluctuating and oscillatory, see upper row) for the Schaefer1000 parcellation in terms of functional connectivity and turbulence. For each of these regimes, we defined a grid of the parameter space (G,β), where G is the coupling strength factor, i.e. the global scaling factor of regional connectivity and β, the shear parameter (see above and Methods). For each pair in the grid, the whole-brain dynamics were simulated 100 times, and we computed the level of fitting between amplitude turbulence (second row) and simulated and empirical FC (third row). PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 2. Model Schaefer1000 fitting of noise, fluctuations, and oscillatory for 1) Turbulence and 2) FC. A-C) We explored the bi-dimensional parameter space defined by β and G for noise, fluctuating and oscillatory regime (bifurcation parameter a = -1.3, a = -0.02 and a = 1.3, respectively, indicated in upper row). We computed the level of amplitude turbulence error as the absolute difference between the empirical and simulated turbulence. Yellow stars indicate the (β, G) combination that reaches the lowest turbulence error in each regime. D) The upper subpanel shows the model fitting scheme in fine Schaefer1000 parcellation (the render on a flatmap of the hemisphere stands for a scheme of brain regions considered in this parcellation). The bottom subpanel displays the barplot that indicates the statistical distribution of the level of amplitude turbulence obtained by simulating 20 trials with 100 subjects for each model regime with the parameters set at the corresponding working point. We also display the results of two model-based surrogates created by increasing the shear parameter of each model regime. The red dashed line indicates the empirical level of amplitude turbulence averaged across participants. The subcritical, supercritical and empirical level of turbulence are not statistically different (Wilcoxon test, P = 0.33), the rest of the comparison are statistically significant (Wilcoxon test, P<0.001). E-G) We explored the bi-dimensional parameter space defined by β and G for noise, fluctuating and oscillatory regime computed the FC fitting as Euclidean distance between the empirical and simulated FC. Yellow stars indicate the (β, G) combination that reaches the lower turbulence error in each regime (the optimal working point obtained in panels A-C). H) The barplot indicates the statistical distribution of the FC fitting obtained by simulating 20 trials with 100 subjects for each model regime at the corresponding working point defined as the minimum turbulence error. We also display the results for the model-based surrogates. All comparisons are statistically significant (Wilcoxon, P<0.001). I) Visualization of the change of the local Kuramoto order parameter, R, in space and time reflecting amplitude turbulence in a single simulation at the optimal working point of each regime (noise, fluctuating and oscillatory cases) and one participant (empirical). This can be appreciated from continuous snapshots for segments separated in time rendered on a flatmap of the hemisphere. https://doi.org/10.1371/journal.pcbi.1010662.g002 We found the optimal fitting for turbulence for each of the three regimes, indicated with a star in the second row of Fig 2A–2C. Fig 2A shows the best fit for the noise regime (a = –1.3) with optimal (G,β) = (1.8,0) as the absolute difference between the simulated and empirical turbulence D (here eD = 0.0473). Fig 2B shows the best fit for the fluctuating regime (a = –0.02) with optimal (G,β) = (1.2, 0.1), which produces a good fit with eD = 4x10-4. Fig 2C shows the best fit for the oscillatory regime (a = 1.3) with optimal (G,β) = (0.15, 2.2), which also produces an good fit with eD = 3x10-4. We also computed the grid fitting for the FC for all three regimes (Fig 2E–2G), with a star in the grid indicating the optimal fit of turbulence which is the criterium for selecting the optimal working point since this is a measure that favours the brain information transmission and its responsiveness [8,10]. Note that these points do not correspond to the optimal fitting with FC in the three regimes. We replicated this analysis by changing the connectivity between non-linear oscillators from the EDR-LR to the structural connectivity (SC) to assess how robust are the results (S1 Fig). We found that despite the levels of fitting and parameters space values are modified, both the fluctuating and oscillatory regimes are good for fitting the level of empirical turbulence in the same region of the (G,β) parameter space. In summary, both the fluctuating and oscillatory regimes are good for fitting the turbulence in the empirical data, while the noise regime is not. This is quantified in Fig 2D, which shows the statistical comparisons between optimal fitting (indicated with the stars in Fig 2A–2C) of the three regimes with turbulence (repeated 20 times 100 simulations) and a horizontal line of D = 0.1976 indicates the level of empirical turbulence. The results show that the best working point for the noise regime is only giving mean D = 0.1484, which is significantly worse than both fluctuating and oscillatory regimes (Wilcoxon p<0.001, compare noise and fluctuation; Wilcoxon p<0.001, compare noise with oscillatory in Fig 2D). On the other hand, the fluctuating and oscillatory regimes are good at fitting the data (mean D = 0.1972 and mean D = 0.1973 respectively) but not significantly different (Wilcoxon P = 0.33, comparing second with the third bar in Fig 2D). Further bolstering these findings, we also generated model-surrogates (see Methods) to compare with the corresponding optimal working point by setting the parameters of the model in the optimal working point but increasing β, which is known to suppress turbulence [8]. Hence, we produced two surrogate models: surr_fluct for the fluctuating model surrogates using a = -0.02 and (G,β) = (1.2, 6) and surr_osc for the oscillatory model surrogates using a = 1.3. and (G,β) = (0.15, 6). The results clearly show significant differences comparing with the level of turbulence fitting obtained by the optimal working point of the model in different regimes (Wilcoxon p<0.001, comparing fluctuation with surr_fluct and comparing oscillations with surr_osc). Similarly, we fitted the whole-brain model with the functional connectivity by means of the Euclidean distance with the empirical. In Fig 2E–2G, we show the fitting for (G,β) for the noise, fluctuating and oscillating regimes. We quantify the fit (using the optimal points from the turbulence fitting indicated with the stars in Fig 2A–2C) in Fig 2H, which shows the statistical comparisons of the three regimes with functional connectivity (see Methods). The results show that the best working point for the turbulence fitting for the noise regime is only giving a functional connectivity fitting mean ErrFC = 0.2594, which is significantly worse than both fluctuating and oscillatory regimes (Wilcoxon p<0.001, comparing noise with fluctuations [mean ErrFC = 0.1422]; Wilcoxon p<0.001, comparing noise with oscillations [mean ErrFC = 0.2063] in Fig 2H). On the other hand, the fluctuating and oscillatory regimes better fit the functional connectivity than the noise regime. However, in this case, the fluctuating regime is significantly better than the oscillatory regime (Wilcoxon p<0.001, comparing fluctuations with oscillations box in Fig 2H). We also evaluated the functional connectivity fitting for the same model-surrogates generated previously (see Methods) to compare with the corresponding optimal working point. The results clearly show significant differences with the obtained level of fitting with the optimal working point of the models (Wilcoxon p<0.001, comparing box fluctuations with surr_fluct and box oscillations with surr_osc). Finally, in Fig 2I, we demonstrate the amplitude turbulence (the local Kuramoto parameter, R, see Methods) at the optimal fitting point of the three whole-brain model regimes contrasted with the empirical data (right subpanel) by rendering continuous snapshots for segments separated in time rendered on a flatmap of a brain hemisphere.

Modelling results for less fine parcellation with 68 regions Following the precise results of fitting the whole-brain model to the empirical data using a fine parcellation, we turned our attention to showing the fitting a less fine parcellation. We found that this was also not able to distinguish between fluctuating subcritical and oscillatory supercritical regimes. Specifically, we found that the level of fitting the empirical metastability defined as the standard deviation of the global Kuramoto order parameters is the same for both regimes. We used during this second analysis a smaller brain parcellation, the Desikan-Killiany with 68 cortical regions of interest (ROIs), to be able to establish a node-level perturbative in silico protocol. We repeated the fitting procedure by exploring the parameter space (G,β) for the model in fluctuation supercritical and oscillatory subcritical regime but we also included the model considering the same absolute values of a than in fluctuation regime but with the opposite signs, that we called supercritical fluctuations. This parcellation is not suitable for computing amplitude turbulence, as is defined in Kawamura et al. [12] and Deco et al. [8], due to the lack of spatial resolution. We thus fitted the metastability, which is the most similar measure computable in coarser parcellation [16]. We found the pair (G,β) that minimizes the absolute difference between the empirical and simulated levels of metastability. Fig 3 shows the results of fitting the Hopf whole-brain model in the two different regimes (fluctuations and oscillations) but also for the supercritical fluctuations, a~0 and a>0 (see upper row panel A, B and C). We computed for the Desikan-Killiany parcellation (upper row panel G) the fitting in terms of functional connectivity and metastability. For each of these regimes, we defined a grid of the parameter space (G,β), and for each pair in the grid, we simulated 100 times the whole-brain dynamics, and we computed the level of fitting between the metastability (second row) and simulated and empirical FC (third row). PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 3. Model Desikan-Killiany fitting of fluctuations and oscillatory for 1) Metastability and 2) FC fitting in. A-C) We explored the bi-dimensional parameter space defined by β and G for fluctuating, supercritical fluctuations and oscillatory regime (bifurcation parameter a = -0.02, a = 0.02 and a = 1.3, respectively, indicated in the upper row) and computed the level of metastability error as the absolute difference between the empirical and simulated metastability. Yellow stars indicate the (β,G) combination that reaches the lowest metastability error in each regime. D-F) We explored the bi-dimensional parameter space defined by β and G for fluctuating, supercritical fluctuations and oscillatory regime computed the FC fitting as Euclidean distance between the empirical and simulated FC. Yellow stars indicate the (β,G) combination that reaches the lowest metastability error in each regime (the optimal working point obtained in panels A-C). G) The upper subpanel shows the model fitting scheme procedure in coarser Desikan-Killiany parcellation (the render on flatmap of the hemisphere stands for a scheme of brain regions considered in this parcellation). The bottom subpanel displays the barplot that indicates the statistical distribution of the metastability error obtained by simulating 20 trials with 100 subjects for each model regime with the parameters set at the corresponding working point. We also display the results of two model-based surrogates created by increasing the shear parameter of each model regime. The comparison between fluctuations, supercritical fluctuations and oscillations model’s regimes at fitting the metastability shows that the two regimes are equally good (Wilcoxon, fluctuations vs oscillations P = 0.21; fluctuations vs supercritical fluctuations P = 0.26; and supercritical fluctuations vs oscillations P = 0.46), while the rest of the comparisons are statistically significant (Wilcoxon, P<0.001). H) The barplot indicates the statistical distribution of the FC fitting obtained by simulating 20 trials with 100 subjects for each model regime at the corresponding working point defined as the minimum metastability error. We also display the results of two model-based surrogates created by increasing the shear parameter of each model. All comparisons are statistically significant (Wilcoxon, P<0.001). https://doi.org/10.1371/journal.pcbi.1010662.g003 We found the optimal fitting for the level of metastability for each of the three scenarios, indicated with a star in the second row of Fig 3A–3C. Fig 3A shows the best fit for the fluctuating regime (a = -0.02) with optimal (G,β) = (2.2,0) with minimal absolute difference between the simulated and empirical metastability M (here eM = 4x10-3). Fig 3B shows the best fit for the supercritical regime close to the bifurcation, supercritical fluctuations (a = 0.02) with optimal (G,β) = (0.2, 1.4), which fit the empirical data with eM = 7x10-3. Fig 3C shows the best fit for the oscillatory regime (a = 1.3) with optimal (G,β) = (0.4, 2.2), which also produces an good fit with eM = 1x10-3. We also computed the grid fitting for the FC, defined as the Euclidean distance between the simulated and empirical FC, for the tree scenarios (Fig 3D–3F), with a star in the grid indicating the optimal fit of metastability (which is the criterium for selecting the optimal working point since this is the most similar measure to turbulence). Note that these points do not correspond to the optimal fitting with FC. We found that the supercritical fluctuation regime behaves similar to the oscillatory regime showing an optimal working point with small coupling strength (G) and β greater than one but the oscillatory regime provides better fit (Fig 3B). We replicated this analysis by expanding the region in the (G,β) parameter space to the same grid for both regimes. In S2 Fig, we display the level of fitting of turbulence and FC for a grid of G = [0–3.4] and β = [0–2.4] in 0.2 steps in both dimensions equally for the three regimes (we included the noise regime). We found that fluctuations and oscillatory regimes present the optimal working point in the regions we have previously analysed (red squares in S2 Fig). The statistical comparison between the optimal working point of subcritical and supercritical model regimes (defined by the optimal fitting of the level of metastability, indicated with the stars in Fig 3A and 3C) is quantified in Fig 3G and 3H. We simulated 20 times the 100 repetitions of the whole-brain model at each regime working point and compared the level of metastability fitting and functional connectivity fitting (see Methods). The results show that the fluctuating, supercritical fluctuating and oscillatory regimes are good at fitting the metastability and not significantly different (Wilcoxon, fluctuations vs oscillations P = 0.21; fluctuations vs supercritical fluctuations P = 0.26; and supercritical fluctuations vs oscillations P = 0.46, first bars of Fig 3G). On the other hand, the fluctuating regimes better fit the functional connectivity than the supercritical fluctuations and oscillatory regime (Wilcoxon p<0.001, for all comparisons in Fig 3H). Note that the supercritical fluctuations regime is similar to the supercritical regime in terms of parameter space working point and metastability fitting, but it produces a worst fitting in terms of functional connectivity, thus in the following we consider only the oscillatory regime for the study. We also generated model-surrogates to compare with the corresponding optimal working point by setting the parameters of the model in the optimal working point but increasing β. Hence, we produced two surrogate models: surr_fluct for the fluctuating model surrogates using a = -0.02 and (G,β) = (2.2, 3) and surr_osc for the oscillatory model surrogates using a = 1.3. and (G,β) = (0.4, 3). The results clearly show significant differences comparing with the level of metastability fitting obtained by the optimal working point of the model in different regimes (Wilcoxon p<0.001, comparing fluctuation with surr_fluct box and oscillations with surr_osc boxes in the second row of Fig 3C). The same results were obtained for the fitting of the functional connectivity (Wilcoxon p<0.001, comparing fluctuation with surr_fluct box and oscillations with surr_osc boxes of Fig 3F). To complete the analysis, we computed within an equally and extended grid for the three regimes (noise, fluctuations and oscillations) other two metrics comparing the simulated and empirical FC, the traditional Pearson correlation and the structural similarity index (SSIM)[49]. This metric balances sensitivity to absolute (e.g. Euclidean distance) and relative (e.g. correlation distance) differences between the FC matrices because it is based on three different observables: the luminance, the contrast and the structure. We found that despite the correlation is not an optimal metric to constrain the models it shows that the fluctuation regime is performing better than the others (S3 Fig). In the same direction, the SSIM confirms that the fluctuation regime is the best model to fit FC and also is a good metric to define an optimal working point in the parameter space defined by [β,G]. For reference we display the yellow star representing the minimum found in metastability error exploration in Fig 3A and 3C. The results also show that it could not distinguish between fluctuating subcritical and oscillatory supercritical regimes in coarser parcellation in terms of fitting the empirical data. We then focused our analysis on the perturbation response as an approach to disentangle between both models.

Global strength-dependent perturbation distinguishes between fluctuating and oscillatory regimes We implemented a global strength-dependent sustained perturbation that allows us distinguishing between fluctuating and oscillatory regimes for both the fine and coarse parcellations (Fig 4). We generated an in silico stimulus by adding an external periodic force applied equally to all nodes at the optimal working point in both model regimes (see Methods). We varied the strength of the external forcing F 0 from 0 to 0.001 in steps of 0.0001, and for each amplitude, we simulated 100 times the perturbed and unperturbed model signals. We obtained the local and global Kuramoto order parameters (lKoP and gKoP) for the perturbed and unperturbed cases for the fine and coarse parcellation, respectively. We then computed the local and global Susceptibility and absolute Information Capacity as the mean and standard deviation of the subtraction between the perturbed and unperturbed lKoP and gKoP across trials (see Methods). We repeated this computation 20 times and Fig 4C–4F shows the mean and standard deviation across repetitions. The subcritical fluctuating regime shows a rapid increase of the level of local Susceptibility in the fine parcellation (Fig 4C) and the level of global Susceptibility in the coarse parcellation (Fig 4D). The global absolute Information Capability also rapidly increase while the forcing strength increases (Fig 4E and 4F dark colors). We have explored how the asymptotic values in both measures in the fluctuating regime depend on the bifurcation parameters (a). To do so we slightly changed the value of a to keep in the subcritical regime but close to the bifurcation point. We found that the asymptote changes with the value of a: while a is closer to the bifurcation point (a = 0) the system becomes less sensitive in terms of Susceptibility. The Information Capability behaves as opposite, while a is close to the bifurcation the system presents less Information Capability (S4 Fig). This could be interpreted that the fluctuation regime presents an optimal working point in terms of proximity of the bifurcation parameter to the bifurcation point (a = 0) that balances the level of both measures. PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 4. Global and sustained strength-dependent perturbation. A) We applied global strength-dependent, sustained perturbation in Schaefer1000 parcellation, and B) the same perturbation in Desikan-Killiany parcellation. C-D) The evolution of local and global Susceptibility (fine parcellation, panel C and coarse parcellation, panel D, respectively) as a function of perturbation strength. In dark purple is shown the response of the subcritical fluctuating regime, while in light purple, the behaviour of the supercritical oscillating regime. The subcritical regime is clearly more susceptible than the supercritical regime that is almost unaltered by the perturbation. E-F) The evolution of global absolute Information Capacity (fine parcellation, panel E and coarse parcellation, panel F, respectively) as a function of perturbation strength. In dark orange is shown the response of the subcritical fluctuating regime, while in light orange, the behaviour of the supercritical oscillating regime. The subcritical regime clearly changes the Information Capacity with the perturbation strength comparing with the supercritical regime that is almost unaltered by the perturbation. https://doi.org/10.1371/journal.pcbi.1010662.g004 In the supercritical oscillatory regime, the global Susceptibility and absolute Information Capability are constant in both parcellation along F 0 (Fig 4C–4F light colours). It is remarkable that the level of these measurements in this regime keep almost zero for all the strength forcing range, showing that the model in that regime do not respond under this global perturbation.

Local strength-dependent perturbation also distinguishes between fluctuating and oscillatory regimes We then used local strength-dependent sustained and non-sustained perturbations that also allows us distinguishing between fluctuating and oscillatory regimes (Fig 5). This is demonstrated using the less fine parcellation. This reduction of the number of regions allowed us to define a node-by-node perturbative approach. Firstly, we explored the model’s regime response by applying a sustained perturbation, and then we quantified the response to non-sustained external perturbation by the PCI. PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 5. Local and Sustained/non-sustained strength-dependent perturbation. A) The evolution of Susceptibility (first row) and the absolute Information Capability (second row) as a function of the perturbation strength and the perturbed pairs of homotopic nodes. The middle left panel displays the results for the subcritical regime (fluctuations, a = -0.02, G = 2.2 and β = 0), and the middle right panel shows the response of the supercritical regime (oscillations, a = 1.3, G = 0.4 and β = 2.2). The right panels present the perturbative node hierarchy rendered onto the brain cortex for both measures (first and second row) for the case of a perturbation strength of 0.01 indicated with a box in middle left panel. B) Non-sustained PCI: The PCI is obtained by perturbing by pairs of homotopic nodes and different forcing amplitude. In the left column, the PCI results are obtained by perturbing the subcritical model in its corresponding working point with an external periodic force applied by pairs of homotopic nodes as a function of the amplitude of this forcing. In the right column, the same measurement is displayed but, in this case, for the supercritical model in its corresponding working point. The right panel shows the node-perturbative hierarchy in terms of PCI of each region for the maximum value of the forcing amplitude (indicated with black box in the middle-left panel) rendered onto a brain cortex. https://doi.org/10.1371/journal.pcbi.1010662.g005

Susceptibility and information capability after local strength-dependent sustained perturbations We systematically perturbed the model in each regime optimal working point by adding an external periodic force (see Fig 5A and Methods). We performed this in silico stimulation approach by forcing the 34 pairs of homotopic nodes (in the parcellation with 68 nodes) with forcing strength ranging from 0 to 0.02 in 0.001 steps. For each combination of nodes and amplitude, we ran 50 trials with 100 simulations, each computing the global Kuramoto Order parameter (gKoP) for the perturbed and unperturbed case (see Methods). We defined the node-level global Susceptibility and Information Capability as the mean and standard deviation across simulations of the subtraction between the node-perturbed and unperturbed gKoP, and we then averaged across trials. Fig 5A (second and third rows) shows the results for oscillatory (supercritical) and fluctuating (subcritical) regimes for both measurements. As in the global perturbation experiment, we noticed that the supercritical regime shows almost non-response under the perturbation, while the subcritical case presents variations across nodes and forcing amplitude. In this node-level perturbative approach, we can determine a hierarchy of perturbative effects by assessing node-by-node perturbation effect while the forcing amplitude increases. Fig 5A left panels shows a render onto brain cortex for both measurements in the subcritical regimen for 0.01 of forcing strength.

PCI after local strength-dependent, non-sustained perturbations We slightly modified our perturbative approach to bring the in silico stimulation protocol closer to in vivo experiments, as proposed by Massimini and colleagues [25]. We simulated the external perturbation as an additive external force by pairs of homotopic nodes as in previous sections, but in this case, we focused on the response after the perturbation ends. Specifically, we simulated 600 volumes with the perturbation active, and we then evaluated the evolution of the signals in the following 200 volumes without perturbation. To investigate the behaviour of both model regimes, we adapted the PCI as is defined in Massimini et al. [25] to be applied on simulated BOLD signals (see Methods). This index gauges the amount of information contained in the integrated response to an external perturbation. Fig 5B displays the evolution of the PCI, computed as the normalised perturbed algorithmic complexity ( ) minus the background algorithmic complexity , for both model regimes, for each pair of nodes and forcing strength. We found that in the oscillatory regime, the behaviour of the system after the perturbation is almost the same as the behaviour of the system without perturbation ( ), for all nodes and amplitudes (Fig 5B middle right panel). On the other hand, assessing the perturbation of the subcritical regime unveils a node hierarchy of the response under external perturbations (Fig 5B middle left panel). These local responses under perturbations rendered onto the brain cortex for the maximal forcing amplitude is displayed in Fig 5B (right panel). It is remarkable that in the subcritical case, a set of nodes present the strongest response in terms of intensity (low values of PCI) and sensitivity (for lower forcing amplitudes). Most nodes present a moderate response for perturbations with forcing amplitude higher than 0.01, and other nodes remain unaltered.

Regional heterogeneity and node-hierarchy perturbative organization As shown in Fig 6, we also investigated how the node-hierarchy established in the previous section can be related to other sources of regional heterogeneity. We used different external sources of local heterogeneity, the T1w:T2w ratio and the principal component of transcriptional activity of an extensive set of specific brain genes (see Methods). We also compared with the anatomical and functional connectivity strength of each region, computed as and (well-known as Global brain connectivity, GBC), respectively, where C is the anatomical structural connectivity, and FC is the functional connectivity (see Methods). Finally, we compared the PCI node-hierarchy with the one found with global Susceptibility and Information Capability. We observed that the PCI hierarchical organisation is highly correlated with the other two perturbative measures obtained in the study. It is remarkable that this PCI node-hierarchy measure is computed after the perturbation ends, while the Susceptibility and Information Capability are computed during the perturbation. PPT PowerPoint slide

PNG larger image

TIFF original image Download: Fig 6. The correlation between the node-level PCI and different sources of regional-level heterogeneity. A-B) The correlation between node-level PCI and the node-level Susceptibility and Information Capacity are computed with significant negative correlation. C) The correlation between the node-level PCI and the first principal component of genes expression node information was computed. No correlation was found between variables. D) The same occurs in the correlation computed between the node-level PCI and the ratio between the T1/T2 MRI. E) The correlation between the node-level PCI and the node anatomical strength is computed obtaining a significant level of negative correlation. F) The correlation between the node-level PCI and the node functional connectivity strength (GBC) is computed obtaining a significant level of negative correlation. https://doi.org/10.1371/journal.pcbi.1010662.g006

[END]
---
[1] Url: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1010662

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/