Universal and non-universal neural dynamics on small world connectomes: a finite size scaling analysis

TL;DR

This study investigates neural dynamics on small-world networks resembling the human connectome, revealing both universal and non-universal critical behaviors and phase transitions influenced by network topology.

Contribution

It demonstrates how topological parameters affect the criticality and universality class of neural dynamics on small-world networks, bridging the gap in understanding brain criticality.

Findings

Varying topological parameters induces mean field or non-universal critical behavior.

Identifies regions with discontinuous, non-critical phase transitions.

Highlights the role of network topology in brain-like critical dynamics.

Abstract

Evidence of critical dynamics has been recently found in both experiments and models of large scale brain dynamics. The understanding of the nature and features of such critical regime is hampered by the relatively small size of the available connectome, which prevent among other things to determine its associated universality class. To circumvent that, here we study a neural model defined on a class of small-world network that share some topological features with the human connectome. We found that varying the topological parameters can give rise to a scale-invariant behavior belonging either to mean field percolation universality class or having non universal critical exponents. In addition, we found certain regions of the topological parameters space where the system presents a discontinuous (i.e., non critical) dynamical phase transition into a percolated state. Overall these results shed light on the interplay of dynamical and topological roots of the complex brain dynamics.

Equations

$$⟨s⟩=\frac{{\sum }_s^{\prime }{s}^2{N}_s}{{\sum }_s^{\prime }s{N}_s}$$

Full text

I Introduction

The study of brain functional activity has revealed the existence of correlated fluctuations and scale invariance similar to those observed in critical phenomena. Such evidence prompted the conjecture that the large scale organisation of the brain emerges at criticalitybeggs ; chialvo ; mora . Despite the relevance of many of these findings, the small size of the available connectomes Haimovici2013 don’t allow to determine whether the expected finite-size scaling behavior at criticality holds, as well as the associated universality class. In this work we study the dynamical properties of a model defined on a class of small world networks that share some topological features with the human connectome, for a wide range of values of the parameters that define the associated topology and for different system sizes. A finite-size scaling analysis allowed us not only for a robust characterization of criticality, but also to estimate critical exponents and therefore to identify the universality class of the critical dynamic claimed to be relevant for the emergence of the large scale organization of brain activity.

The paper is organized as follows. In section II we describe the model as well as the simulation and finite size scaling methods used. The results are presented in section III and discuss the relevance of the main findings in section IV.

II Model and Methods

II.1 The model

This work uses an adaptation of the neural model presented in Ref.Haimovici2013 running over a small-world network with a weighted adjacency matrix $w_{ij}$. To mimic the weights distribution of the human connectomeHaimovici2013 ; Hangmann2008 , the non null $w_{ij}$ were chosen randomly from an exponential distribution $p(w)=\lambda\,e^{-\lambda\,w}$ with $\lambda=12.5$.

The node dynamics of the neural model respond to the Greenberg-Hastings cellular automatonGreenberg-Hastings , in which each node $i$ of the network has associated a three state variable $x_{i}=0,1,2$, corresponding to the following dynamical states: quiescent ($x_{i}=0$), excited ($x_{i}=1$) and refractory ($x_{i}=2$). The transition rules are the following: if a node at the discrete time $t$ is in the quiescent state $x_{i}(t)=0$ it can make a transition to the excited state $x_{i}(t+1)=1$ with a small probability $r_{1}$ or if $\sum_{j}w_{ji}\,\delta(x_{j}(t),1)>T$, where $T$ is a threshold and $\delta(x,y)$ is a Kroenecker delta function; otherwise, $x_{i}(t+1)=0$. If it is excited $x_{i}(t)=1$ then it becomes refractory $x_{i}(t+1)=2$ always. If it is refractory $x_{i}(t)=2$ then it becomes quiescent $x_{i}(t+1)=0$ with probability $r_{2}$ and remains refractory $x_{i}(t+1)=2$ with probability $1-r_{2}$. Following Ref.Haimovici2013 we set $r_{1}=10^{-3}$ and $r_{2}=0.3$111We made several checks for for different values of $r_{2}\sim 10^{-1}$ and no qualitative differences were observed. We kept these values fixed in all the calculations.

Each node interacts with the others according with the topology of a small-world network, constructed as the usual Watts-Strogatz modelwatts-strogatz . That is, we start from a ring of $N$ nodes in which each node is connected symmetrically to its $2m$ nearest neighbors. Then, for each node each vertex connected to a clockwise neighbor is rewired to a random node with a probability $\pi$ and preserved with probability $1-\pi$, so the average degree $\langle k\rangle=2m$ is preservedBarrat .

II.2 Statistical properties

The present analysis focuses on the dynamical clusters of coherent activity, namely, groups of nodes simultaneously activated ($x_{i}=1$) which are linked through non zero weights $w_{ij}$, for different values of the parameters $(\pi,\langle k\rangle)$ and different network sizes $N$. Each simulation started from a random distribution of activated sites and we let the system to run $100$ time steps before starting to collect data. We found that time interval to be enough for the system to reach a stationary state for any system size and for any value of the network parameters. Then we recorded data every 5 time steps, to avoid artifactual correlation effects. For each data set we computed several measures to describe a percolation like critical phenomenon as a function of the control parameter $T$Haimovici2013 . Specifically, we calculated the average sizes of the largest (i.e., giant) cluster $\overline{S}_{1}$ which can be considered the order parameter and also we compute the average size of the second largest cluster $\overline{S}_{2}$, together with the average cluster size,

[TABLE]

where the primed sum runs over all cluster sizes except the giant one and $N_{s}$ is the number of clusters of size $s$ Barrat ; Margolina . In some cases we also computed the cumulative complementary distribution function (CCDF) of clusters sizes. Depending of the case, averages were computed both over data obtained along the same simulation run and over different networks.

If the system presents a percolation like critical point, both $\langle s\rangle$ and $\overline{S}_{2}$ are expected to exhibit a (size dependent) maximum for certain pseudo critical value of the control parameter (the threshold $T$ in the present case), that scale with the system size asStauffer $\langle s\rangle\sim N^{\gamma/\nu d}$, $\overline{S}_{2}\sim N^{d_{f}/d}$. $\gamma$ and $\nu$ are the standard susceptibility and correlation length critical exponents, $d$ is the effective dimension of the system and $d_{f}$ the fractal dimension of the percolating cluster. Also at the critical point it is expected thatStauffer $P(s)\sim s^{-\tau}exp({-s/S^{*}})$, where $S^{*}\propto\overline{S}_{2}$.

III Results

As already mentioned, the distribution of cluster sizes of activity can be informative of the different dynamical regimes of the model. We computed such measures as a function of the threshold $T$ for models of increasing sizes $N$, and different topologies by varying average degree $\langle k\rangle$ and rewiring $\pi$. The model dynamics can be explored by inspecting the clusters’ finite-size scaling behaviour as a function of these three parameters. Figures 1, 2 and 3 correspond to the three typical behaviour found at well defined regions in parameter space which we will describe in detaill now.

At relatively low average degree $\langle k\rangle$, the example of Fig.1 is typical. We see that $\overline{S}_{1}/N$ goes to zero for any value of $T$ as $N\to\infty$. Also both $\overline{S}_{2}/N$ and $\langle s\rangle$ display a monotonous decrease (no peak) for any value of the system size $N$. Thus, for those small world topologies only a few small and uncorrelated clusters are present at any time for any value of the threshold $T$, resembling what TononiTononi describes as segregated regime, which per se is incompatible with normal brain functioning.

Increasing the value of $\langle k\rangle$ result in a change of the model dynamics. For such case, as shown in Fig.2, the clusters analysis reveals the characteristic finite-size scaling behavior exhibited by percolation phenomenaStauffer . There, we observe the existence of a critical threshold $T_{c}$ below which the order parameter converges to a finite value when $N\to\infty$, while it goes to zero above the critical threshold (Fig.2a). This means the existence of a dynamical phase transition, where the system develops long range order in the form of macroscopic correlated cluster of simultaneously active nodes, which resembles the correlations associated with the brain RSN described in [Haimovici2013 ] . Also both $\overline{S}_{2}$ and $\langle s\rangle$ exhibit a sharp peak around $T=T_{c}$ (See Fig.2b and c) that diverges as a power law when $N\to\infty$, thus exhibiting well defined critical exponents. Hence, for the associated topologies, the dynamical phase transition is indeed a second order or critical one.

A very different transition happens for large enough values of both parameters $(\pi,\langle k\rangle)$. Indeed, in that region both the order parameter and the fluctuations measures $\overline{S}_{2}$ and $\langle s\rangle$ exhibit an apparently discontinuous behavior at the transition point, as illustrated in Fig.3.

By analyzing the finite size scaling of any of the three quantities $\overline{S}_{1}/N$, $\overline{S}_{2}$ and $\langle s\rangle$ in the large $N$ limit we estimated the frontiers in the $(\pi,\langle k\rangle)$ space between the regions where the system exhibit or o not a dynamical phase transition, as well as the nature of the transition (see Fig.4). The frontier between the segregated and the Mean Field Critical (MFC) regimes was obtained by analizing whether there is a $T_{c}>0$ such that $\overline{S}_{1}/N$ converges to a finite value for $T<T_{c}$ when $N\to\infty$. The frontier between MFC and the Non Universal Critical (NUC) regime was identified by analyzing whether the difference between critical exponents $\gamma/\nu d$ and $d_{f}/d$ (obtained from the finite size scaling extrapolation) and the corresponding percolation mean field values exceeds or not the associated error bars. The frontier between NUC and the oscillatory regime was estimated by analyzing the presence or not of a discontinuity in the order parameter at the transition in the large $N$ limit.

At relatively low average degree $\langle k\rangle$ there is no transition (i.e., the region labeled NT in Fig.4). When present, the nature of the dynamical phase transition can be different depending on the topology of the network, as depicted in Fig.4.

For relatively small values of $\langle k\rangle$ and/or small enough values of $\pi$, the behavior is critical and the corresponding critical exponents are consistent with the universality class of mean field (or Bethe) classical percolation, namely $\gamma/\nu d\approx 1/3$ and $d_{f}/d=2/3$ (corresponding to $d=6$, $d_{f}=4$, $\nu=1/2$ and $\gamma=1$ ). An example is the case already shown in Fig.2.

To get further insight about the nature of the dynamical phase transition in the last mentioned region, we analyzed the associated behavior of the cluster size distribution (see Fig.5). For the critical dynamics, consistently, the cumulative cluster size distribution exhibits the expected behavior $P_{c}(s)\equiv\sum^{\prime}_{s^{\prime}\geq s}P(s^{\prime})\sim s^{-(\tau-1)}exp({-s/S^{*}})$, with an exponent $\tau\approx 5/2$ and $S^{*}\propto\overline{S}_{2}$(thus satisfying the scaling law $\tau=d/d_{f}+1$) as shown in Fig.4, corresponding to the region in the $(\pi,\langle k\rangle)$ space denoted as MCF (Mean Field Critical).

For even larger values of $\pi$ and/or $\langle k\rangle$ a wide region of the parameter space results in critical exponents changing continuously with the topological parameters, until they seems to saturate (see Fig.6). On the other hand, the power law increase of values of $\overline{S}_{2}$ (see Fig.3B) at both sides of the discontinuity with an exponent close to one third implies a roughly two dimensional percolating cluster ($d_{f}\approx 2$), assuming an effective dimension $d=6$.

For comparison, we also calculated the clusters size distribution for a set of values in the NUC region, as shown in Fig.5B. The power law distribution is indeed consistent with critical behavior. On the other hand, the data on Fig.5C shows that the cluster size distribution for large values of $(\pi,\langle k\rangle)$, for large enough system sizes ($N>10^{4}$), exhibits an exponential decay, which exclude critical behavior. Hence, we can conclude that in such region the system exhibits a first order or discontinuous dynamical percolation transition. In this region the analysis of the autocorrelation function of the mean activity (see Inset Fig. 3a) reveals that the dynamic exhibits a finite time scale corresponding to an oscillation whose period is a function of the topological parameters $\pi$ and $<k>$ as well as $T$.

The results of the scaling analysis are summarized in Fig.6 which show the computed critical exponents both for the maximum of $\overline{S}_{2}$ in panel A and for the maximum of $\langle s\rangle$ in panel B. The data shows the two behavior already discussed: the universal (mean field, denoted with filled circles) as well as for the cases in which there is a continuously changing exponent depending on the $\pi$ value.

IV Discussion

The main results of the present work illustrate how the large scale correlated activity depend on the topological details of the underlying structural network. The results show that the correlated patterns exhibiting scale invariance seen in the brain RSN need certain minimum of connectivity conditions. In other words, using Tononi’s terminologyTononi , for relatively low average degree $\langle k\rangle$ and fraction of network shortcuts $\pi$ the neuronal activity is too segregated. In the other extreme, for very large values of $\langle k\rangle$ and $\pi$ the activity is too integrated, exhibiting a first order or discontinuous phase transition. This is most unexpected, because in conventional static percolation on small world networks the transition is second order, belonging to the mean field universality class for any value of the rewiring probability $\pi$ Moore . Hence, the present effect is due entirely to dynamical correlations. In between these two extremes we found a regime that seems compatible with the fMRI brain data. There, two different dynamical regimes are found which are characterized by different kinds of dynamical phase transitions as the threshold controlling the neuronal activation is varied. For certain values the phase transition is universal and for others, relatively larger, while the scale invariance persist the universality is lost,i.e. the exponents change continuously. Such regions could correspond, with some differences, to the equivalent to the Griffiths phasesrare1 already described in another contextrare2 ; moretti . Please notice some recent results in the context of the brain connectome odor1 ; odor2

From the neuroscience point of view, the relevant finding is that the critical regime in this class of networks spans a wide region in the parameters space, corresponding to intermediate values of networks topologies, those that are too connected or too disconnected are not able to exhibit critical dynamics, regardless of the values of the other two parameters (system size and node’s threshold). In other words, from the point of view of neural architecture a minimum of connectivity needs to be predetermined (perhaps via evolution) in order for the dynamics to achieve the dynamical features of criticality (via modulation of excitability or threshold). On the other hand, the discontinuous percolation that appears only for extremely connected networks, might constitute pathological conditions observed in real nervous systems. Hence, critical fluctuations emerge as a robust characteristic towards variations in the anatomical network topology which is consistent with the expectations of criticality described for the spontaneous fluctuations of brain dynamics.

Acknowledgements.

This work was partially supported by CONICET (Argentina) through grant PIP 11220150100285, Ministerio de Ciencia y Tecnología de la Provincia de Córdoba (Argentina) through grant PID2018 and SeCyT (Universidad Nacional de Córdoba, Argentina) and by the NIH (USA) Grant 1U19NS107464-01. MZ is a recipient of a Doctoral Fellowship from CONICET (Argentina). This work used Mendieta Cluster from CCAD-UNC, which is part of SNCAD-MinCyT, Argentina.