Effective networks: a model to predict network structure and critical transitions from datasets

TL;DR

This paper introduces an effective network model that predicts the structure and critical transitions of complex systems like neuronal networks using data-driven reconstruction, even outside observed parameter ranges.

Contribution

It presents a novel approach to predict network behavior and critical transitions from limited data by combining local dynamics with statistical interaction models.

Findings

Successfully reconstructed neuronal network structure from data

Predicted critical transitions outside observed parameter ranges

Validated approach on cat cerebral cortex networks

Abstract

Real-world complex systems such as ecological communities and neuron networks are essential parts of our everyday lives. These systems are composed of units which interact through intricate networks. The ability to predict sudden changes in network behaviour, known as critical transitions, from data is important to avert disastrous consequences of major disruptions. Predicting such changes is a major challenge as it requires forecasting the behaviour for parameter ranges for which no data on the system is available. In this paper, we address this issue for networks with weak individual interactions and chaotic local dynamics. We do this by building a model network, termed an effective network, consisting of the underlying local dynamics at each node and a statistical description of their interactions. We illustrate this approach by reconstructing the dynamics and structure of realistic neuronal interaction networks of the cat cerebral cortex. We reconstruct the community structure by analysing the stochastic fluctuations generated by the network and predict critical transitions for coupling parameters outside the observed range.

Equations

$${\mathbf{x}}_i(t+1)$$

$${\mathbf{G}}_i({\mathbf{x}}_i)={\mathbf{F}}_i({\mathbf{x}}_i)+{\beta }_i\mathbf{V}({\mathbf{x}}_i(t))$$

$${\mathbf{x}}_i(t+1)={\mathbf{G}}_i({\mathbf{x}}_i(t))+{\mathbf{\xi }}_i(t)$$

$$y_i(t+1)=g_i(y(t))+{\xi }_i(t)$$

$$g_i=f_i(y_i(t))+{\beta }_{i}v(y_i(t))$$

$$r(t)e^{i\psi (t)}=\frac{1}{N_c}\sum _{j=1}^{N_c}{e}^{i{\theta }_j(t)}$$

$$\alpha \sum _{j=1}^N{A}_{ij}\mathbf{H}({\mathbf{x}}_j(t),{\mathbf{x}}_i(t))=\alpha k_i\mathbf{V}({\mathbf{x}}_i(t))+\alpha k_i{\mathbf{\xi }}_i(t)$$

$$\mathbf{V}(\mathbf{x})=\int \mathbf{H}(\mathbf{x},\mathbf{y})d\mu (\mathbf{y})$$

$${\mathbf{x}}_i(t+1)={\mathbf{G}}_i({\mathbf{x}}_i(t))+\alpha k_i{\mathbf{\xi }}_i(t)$$

$$y_i(t+1)=g_i(y_i(t))$$

$${\hat{\mu }}_{i\ell }=\frac{|\Gamma (i)\cap \Gamma (\ell )|}{|\Gamma (i)\cup \Gamma (\ell )|}$$

$${\hat{\mu }}_{i\ell }=\frac{(A+A^2{)}_{i\ell }}{k_i+k_{\ell }-(A+A^2{)}_{i\ell }}.$$

$$\text{mbox}Corr({\xi }_i^c,{\xi }_j^c)\propto {\hat{\mu }}_{ij}.$$

$$\text{mbox}Corr(y_i,y_j)\approx \frac{\text{mbox}Cov(\zeta ,\zeta )}{{\sigma }_{Y_1}^2{\sigma }_{Y_2}^2}$$

$$F_1(u,w)=\frac{\beta }{1+u^2}+w\text{mbox}and{F}_2(u,w)=w-\nu u-\sigma .$$

$$\Gamma (u_j)=1/(1+\exp \{\lambda (u_j-{\Theta }_s)\})$$

$$h(u_i,u_j)=u_j-u_i$$

$$\theta (t)=2\pi \left(\frac{t-t_n}{t_{n+1}-t_n}+t_n\right),t_n<t<t_{n+1}$$

$$\sum _j{A}_{ij}\mathbf{H}({\mathbf{x}}_i,{\mathbf{x}}_j)=\sum _{j\in C}{A}_{ij}\mathbf{H}({\mathbf{x}}_i,{\mathbf{x}}_j)+\sum _{j\notin C}{A}_{ij}\mathbf{H}({\mathbf{x}}_i,{\mathbf{x}}_j)$$

$$\sum _{j\notin C}{A}_{ij}\mathbf{H}({\mathbf{x}}_i,{\mathbf{x}}_j)=(1-\kappa )\Delta \int \mathbf{h}({\mathbf{x}}_i,\mathbf{y})d\mu (\mathbf{y})+{\mathbf{\xi }}_i^o(t)$$

$${\mathbf{x}}_i(t+1)={\mathbf{q}}_i(x_i(t))+\sum _{j\in C}{A}_{ij}\mathbf{H}({\mathbf{x}}_i,{\mathbf{x}}_j)+{\mathbf{\xi }}_i^o(t),$$

$$w(t+1)=w_0+\mu \sum _{n=0}^t(u(n)-1)$$

$$u_i(t+1)=F_{1,i}(u_i(t))+u_i(t)+\alpha \Delta [⟨u⟩-u_i(t)]+{\xi }_i(t)$$

$$\zeta (t)=u_i(t)-u_j(t)$$

$$\zeta (t+1)=D{F}_1(x_i(t))\zeta (t)+\mu \sum _{n=0}^t\zeta (n)-\alpha \Delta \zeta (t)$$

$$\eta (t)=\sum _{n=0}^{t}z(n)$$

$$\sum _{n=0}^{t}D{F}_1(x_i(n))\zeta (n)\approx \lambda \sum _{n=0}^t\zeta (n)$$

$$\eta (t+1)=(\lambda -\alpha \Delta )\eta (t)+\mu \sum _{n=0}^t\eta (n)$$

$$\displaystyle\Theta=\left(\begin{array}[]{cccc}\phi_{1}(x_{1}(t_{1}))&\phi_{1}(x_{2}({t_{1}}))&\cdots&\phi_{k}(x_{N}({t_{1}}),x_{N}({t_{1}}))\\ \phi_{1}(x_{1}(t_{2}))&\phi_{1}(x_{2}({t_{2}}))&\cdots&\phi_{k}(x_{N}({t_{2}}),x_{N}({t_{2}}))\\ \vdots&\vdots&\vdots&\vdots\\ \phi_{1}(x_{1}(t_{n-1}))&\phi_{1}(x_{2}({t_{n-1}}))&\cdots&\phi_{k}(x_{N}({t_{n-1}}),x_{N}({t_{n-1}}))\\ \end{array}\right)$$

$$\Xi ={\left({\zeta }_{1,1},{\zeta }_{1,2},\cdots ,{\zeta }_{1,N},{\zeta }_{2,1},{\zeta }_{2,2},\cdots {\zeta }_{kN}\right)}^{\ast }$$

Full text

• 1

Instituto de Ciências Matemáticas e Computação, Universidade de São Paulo, São Carlos, Brazil

• 2

Department of Mathematics, Imperial College London, London SW7 2AZ, UK

• 3

Department of Bioinformatics and Genetics, Kadir Has University, 34083 Istanbul, Turkey

• 4

Department of Mathematics and Statistics, University of Victoria, Victoria BC, CA

We are surrounded by a range of complex networks consisting of a large number of intricately coupled nodes. Neuron networks form an important class examples where the interaction structure is heterogeneous [1]. Because changes in the interaction can have massive ramifications on the system as a whole, it is desirable to predict such disturbances and thus enact precautionary measures to avert potential disasters. For instance, neurological disorders such as Parkinson’s disease, schizophrenia, and epilepsy, are thought to be associated with anomalous interaction structure among neurons [2]. Often, as in the case of neuron networks, it is impossible to directly determine the interaction structure. Therefore, a major scientific challenge is to develop techniques which use measurements of the time evolution of the nodes state to indirectly recover the network structure and predict the network’s behaviour when essential characteristics of the interactions change.

The literature on data-based network reconstruction is vast. Reconstruction methods can be classified into model-free methods and model-based methods. The former identify the presence and strength of a connection between two nodes by measuring the dependence between their time-series in terms of: correlations [3, 4], mutual information [5], maximum entropy distributions [6, 7], Granger causality, and causation entropy [8, 9]. Such methods alone do not provide information on the dynamics of the network, which is necessary to predict critical transitions. Model-based methods instead provide estimates (or assume a priori knowledge) of the dynamics and interactions in the system, and use this knowledge to reconstruct the network topology. When the interactions are sufficiently strong, the network structure can be recovered [11, 12, 10]. For a more extensive account of reconstruction (model-free and -based) methods see the reviews [10, 13, 14] and references therein.

In many real-world applications, the behaviour of isolated nodes is erratic (chaotic) and the interaction is weak [15, 16, 1]. Moreover, the network structure typically has community structures and hierarchical organisations such as the rich-club networks found in the brain [17]. As the interaction strength per connection is weak and the statistical behaviour of the nodes robust, the influence of a single link on the overall dynamics is negligible and only the cumulative contribution of many links matter. Furthermore, because of the chaotic dynamics and weak coupling, the influence of a single node on the rest of the network corresponds essentially to a random signal that is superposed to the randomness generated by the chaos in the local dynamics. The existing reconstruction techniques fail to reconstruct a model from the data, as they require the interaction to be of the same magnitude as the isolated dynamics.

In this paper, we introduce the notion of an effective network. Its aim is to make a model of a complex system from observations of the nodes evolution when the network has heterogeneous structure, the strength of interaction is small and local dynamics are highly erratic. This approach starts by reconstructing the local dynamics from observations of nodes with relatively few connections, and then recovers the interaction function from observations of the highly connected nodes whose dynamics are the most affected by the interactions as a result of the multitude of connections they receive from the rest of the networks[18, 19]. A key achievement is that we are able to identify community structures even in the presence of weak coupling. An outcome of the effective network is that it recovers enough information to forecast and anticipate the network behaviour, even in situations where the parameters of the system change into ranges that have not been previously encountered.

Complex networks of nonlinear systems. We consider networks with $N$ nodes with chaotic isolated dynamics and pairwise interactions satisfying some additional assumptions described below. The network is described by its adjacency matrix $\bm{A}$, whose entry $A_{ij}$ equals $1$ if node $i$ receives a connection from $j$ and equals [math] otherwise. We assume that the time evolution of the state $\bm{x}_{i}(t)$ of node $i$ at time $t$ is expressed as

[TABLE]

When performing reconstruction from data, the isolated local dynamics $\bm{F}_{i}\colon M\to M$, the coupling function $\bm{H}$, the coupling parameter $\alpha$ (that is small), the adjacency matrix $\bm{A}$, and even the dimension $k$ of the space $M$, are assumed to be unknown. This class of equations can be applied to modelling many important complex systems found in neuroscience (e.g. coupled equations in neuron networks) [20], engineering (e.g. smart grids) [21, 22], material sciences (e.g. superconductors) [23], and biology (e.g. cardiac pacemaker cells) [24], among other systems.

Our three assumptions are: $(a)$ the local dynamics are close to some unknown ergodic and chaotic map $\bm{F}$ (that is, $\|\bm{F}-\bm{F}_{i}\|\leq\delta$, which is often the case in applications [25, 26]). $(b)$ The network connectivity is heterogeneous, which means that the number of incoming connections at a node $i$ (given by its degree $k_{i}=\sum_{j}A_{ij}$) varies widely across the network. More precisely, $k_{i}$ is larger for a few nodes called hubs. $(c)$ $\alpha$ is such that, denoting by $\Delta=\max_{i}k_{i}$ the maximum number of connections $\alpha\Delta$ is of the same magnitude of $\bm{DF}$. Assumptions $(a)$ and $(c)$ imply that the the total effect of the network on a node $i$, which is of the order $\alpha k_{i}$, ranges from small (for nodes $i$ with few incoming connections) to of order 1 for the hub nodes. A prime example is the cat cerebral cortex which possesses inter-connected regions split into communities with a hierarchical organization. This network has heterogeneous connectivity, chaotic motion and weak coupling [27, 28, 29]. Other examples, includes the drosophila optic lobe network [30, 31]. For a given dataset, our effective network first tests whether the underlying system satisfies assumptions $(a)-(c)$, and, if so, reconstructs the model (see the Methods for more detail).

We assume the availability of a time series of observations $y_{i}(t)=\phi({\bm{x}}_{i}(t))$, where $\phi$ is a projection to a variable on which unit interactions depend. This situation occurs frequently in applications; for example, with measurements of membrane potentials in neurons or voltages in electronic circuits. Here, we demonstrate how to obtain a model for the system by constructing an effective network from the time series $y_{i}(t)$.

Effective networks recover both structure and dynamics

To obtain an effective network from observations of a complex system, we combine statistical analysis, machine-learning techniques, and dynamical systems theory for networks. An effective network provides (i) local evolution laws and averaged interactions for each unit that, in combination, closely approximate the unit dynamics, and (ii) a network with the same degree distribution and communities as the original system. We use the term “effective” because, even if it does not carry information on each link and interaction in the original system, the network gathers sufficient data to reproduce the behaviour of the original network and predict its critical transitions.

Using our assumptions for the network and local dynamics, we can show that the evolution at each node will have low-dimensional excursions over finite time scales. In particular, the approximated evolution rule at node $i$ has a low dimensional description over large time scales given by . More precisely, the evolution rule at node $i$ is given by

[TABLE]

where $\bm{F}_{i}\approx\bm{F}$ is the isolated dynamics, $\beta_{i}=\alpha k_{i}$ is the rescaled degree, and $\bm{V}$ is the averaged mean-field effect of the interaction function (effective coupling function) that takes into account the cumulative effect of interactions on node $i$. The true dynamics $\bm{x}_{i}(t)$ slightly deviates from this rule and is given by

[TABLE]

where $\bm{\xi}_{i}(t)$ is a fluctuation term that is small for an interval of time which is exponentially large in terms of the size of the network. This low-dimensional reduction has been rigorously established in important test cases over long time scales (see [19]). This approximation holds, roughly speaking, for two reasons. Firstly, a low degree node $i$ has a small number of connections compared to the maximum degree, and the interactions slightly perturb its statistical behaviour. Secondly, hub nodes receive a huge number of connections so that the sum of their interactions with the other nodes in the network, $\alpha\sum_{j}A_{ij}\bm{H}(\bm{x}_{i},\bm{x}_{j})$, can be approximated as the integral $\alpha k_{i}\int\bm{H}(\bm{x}_{i},\bm{y})\rho(\bm{y})\,d\bm{y}$ where $\rho(\bm{y})$ is the stationary distribution of typical orbits of the low degree nodes. This approximation holds up to a small fluctuation, $\bm{\xi}(t)$, which depends on the state of neighbours of the $i$th node. This fact will be key to detecting communities.

The approximation described above also applies to the measured state variable $y_{i}(t)$. Depending on the system, we need to preprocess the data (see Methods). The processed variable is still referred to as $y_{i}(t)$. The Takens reconstruction technique tells us that $y_{i}(t+1)$ is a nonlinear function of $k+1$ past points $y_{i}(t),\dots y_{i}(t-k)$, for a given number $k$ provided by the approach. Here, we focus on the case when $k=1$, which occurs in many real-world examples, and discuss cases with $k\geq 2$ in the Supplementary Materials. This means that

[TABLE]

where

[TABLE]

is a low-dimensional law that approximates the evolution of $y_{i}(t)$. We first employ Takens reconstruction to estimate $g_{i}$ and coarsely classify the nodes by their degrees. For instance, in rich-club motifs, the network has low-degree nodes organised in communities and high-degree in the rich-club. Since $\beta_{i}=\alpha k_{i}$ is small for low-degree nodes $i$, the evolution rules at such nodes will be similar to the isolated dynamics. Thus we identify which nodes are low-degree nodes and we can recover the local dynamics $f$. Next, we use Eq. (2) and a classification of nodes by their time-series, to obtain the coupling function and estimate the degree distribution $k_{i}$ and the coupling parameter $\alpha$. More precisely, an effective network is obtained in three main steps (see Methods for additional details):

• (i)

Reduced dynamics. Since the fluctuation term $\xi_{i}(t)$ is small, the rule $g_{i}$ can be estimated fitting the points $(y_{i}(t),y_{i}(t+1))$ using machine-learning techniques. Here, we estimate the function $g_{i}$ by decomposing it as a linear combination of basis functions, which is chosen according to the application. The parameters of the basis functions are obtained by performing a $10$-fold cross-validation with $90\%$ training and $10\%$ test [32, 33] (see Methods for details). The particular technique used in this step is not essential as the dynamics will be low-dimensional. Other techniques such as compressive sensing based techniques [34, 35, 36] could be also employed.

• (ii)

Isolated dynamics and effective coupling. Since the network is heterogeneous, it has many low degree nodes for which $g_{i}$ is close to the uncoupled dynamics $f_{i}\approx f$. We first run a model-free estimation to coarsely classify nodes into low-degree nodes and hubs and estimate $f$ (see Supplementary Material). We then use the expression of $g_{i}$ recovered at low degree nodes to obtain an approximation for $f$, while $g_{i}$ recovered at hub nodes allows to estimate $\beta_{i}v$. We estimate the parameter $\beta_{i}$ by using Bayesian inference (see Methods).

• (iii)

Network structure and communities. Since $\beta_{i}=\alpha k_{i}$, we can recover the network’s degree distribution from knowledge of $\beta_{i}$ at every node. To reconstruct the community structures, we use the fact that the term $\xi_{i}(t)$ at a node $i$ depends on the state of all nodes interacting with $i$. Thus the correlation between $\xi_{i}$ and $\xi_{j}$ is proportional to the matching index (number of common neighbors) of the nodes $i$ and $j$. Thus, we can create a network that has the same statistical properties of the actual one. Further details are provided in the next section.

Benchmark model for the isolated dynamics. We consider heterogeneous networks of neurons described by the Rulkov model (see Methods for details). The model has two variables, $u$ and $w$, evolving at different time scales. The fast variable $u$ describes the membrane potential and $w$ the slow currents. In the following discussion, we choose the parameters such that the isolated dynamics presents a tonic spiking behaviour. In the Supplementary Materials, we also show the results for bursting dynamics, i.e. when the parameters are such that the membrane potential oscillates between two phases: a bursting phase of a rapid spiking, and a quiescent phase. We also validated our methods for other types of isolate chaotic dynamics (also in the Supplementary materials).

Revealing community structure: the rich-club motif

We apply the effective network in a wide range of setups as discussed in the supplementary material. Here we focus on the network structure of the cat cerebral cortex for which the mesoscopic connectivity has been probed and detailed [29]. The regions and their connections were discovered by using datasets from tract-tracing experiments analyzed via nonmetric multidimensional scaling, an optimization approach that minimizes the distance between connected structures and maximizes between the unconnected ones [27, 28].

The network contains 53 meso-regions arranged in four communities that closely follow functional subdivisions; namely, visual (16 nodes), auditory (7 nodes), somatomotor (16 nodes) and frontolimbic (14 nodes), as shown in Fig 2 (a). In addition to these cortical regions, some cortical areas (hubs) form a hidden layer called a rich-club and are densely connected to each other and the communities. A set of nodes form a rich-club if their level of connectivity exceeds what would be expected by chance alone. The maximum number of connections in this network is $\Delta=37$. The network obtained is weighted [27, 28]. For simplicity and to improve the performance in detecting communities, we turn the network into an undirected simple graph [29]. We simulate each mesoregion as a neuron interacting via electrical synapses (see Methods for details) and obtain a multivariate data $\{y_{1}(t),y_{2}(t),\dots,y_{N}(t)\}$ for a time $T=5000$. For simplicity, we will denote $y_{i}=\{y_{i}(t)\}_{t=0}^{T}$, whenever there is no risk of confusion.

For comparison, we first test the recovery of the network using the functional network framework [39, 37, 38]. The intuition behind this approach is that nodes with similar time series perform analogous functions in the network and are assumed to have similar characteristics. The functional network can thus be constructed by the matrix of the similarities between nodes via statistical analysis [40, 41]. Here we employ a covariance analysis between the time series. The functional network alone does not detect communities from the data because the bursts are erratic (Fig. 2 (b)). Other similarity measures such as mutual information give no significative improvement.

We also compared our method with a reconstruction by the sparse recovery technique [34, 10]. There, the dynamics is written as linear combination of chosen basis functions and the unknowns are the coefficients. Sparse recovery works well when many coefficients are exactly zero. The presence of a link is determined when any coefficient of the corresponding interaction is nonzero. In our setting, each link provides a small contribution to the individual dynamics and only their cumulative effect can influence the dynamics. Therefore, the coefficients to be recovered are close to zero, and cannot be distinguished from zero terms. In this setup, sparse recovery can also uncover a model for the isolated dynamics, however, it misidentifies the network structure. We implemented the sparse recovery method to our benchmark method and could assign only $20\%$ of the nodes to the right cluster. These results are discussed in the Methods.

Remarkably, however, the effective network is able to recover the community structures (Fig. 2 (c)). To every pair of $y_{i}$ and $y_{j}$ we assign a Pearson distance $s_{ij}\geq 0$ which measures the difference between the observed dynamics at nodes $i$ and $j$ so that $s_{ij}\approx 0$ if the attractors of $i$ and $j$ are similar and $s_{ij}\approx 1$ if they are reasonably distinguishable (see Methods for details). The higher the number of nodes with behaviour different from $i$, the larger the intensity $S_{i}=\sum_{j}s_{ij}$. The analysis of intensities $S_{i}$ allows us to distinguish nodes in the rich-club. Next, we find the local low-dimensional rules $g_{i}$ for each time series $y_{i}$ by applying Step (i). These rules are similar and display a chaotic evolution, thus, confirming that the similarity analysis by the functional networks would provide little insight into the reconstruction. Having the local rules $g_{i}$, we can decompose the time series in terms of a low-dimensional deterministic part and a fluctuation term $\xi_{i}$, which depends on the neighbours of the $i$th node. One of our main key points is that these fluctuations allow us to reveal the community structure.

Indeed, if nodes $i$ and $j$ interact with the same nodes, they are subject to the same inputs and the covariance Cov$(\xi_{j},\xi_{i})$ is high. If not, Cov$(\xi_{i},\xi_{j})$ is nearly zero due to the decay of correlations in the deterministic part. Therefore, Cov$(\xi_{j},\xi_{i})$ is high when nodes $i$ and $j$ have high matching index (high fraction of common connections). It is crucial that the correlation analysis is restricted to the dynamical fluctuations $\xi_{i}$. Indeed, since the variance of the deterministic part of $y_{i}$ is larger than that of the small fluctuations $\xi_{i}$, performing a direct correlation analysis between $y_{i}$ and $y_{j}$ hides all the contributions to the correlations coming from the covariance between $\xi_{i}$ and $\xi_{j}$. Consequently, the correlation of the deterministic part is close to zero due to the chaotic nature of the dynamics (see Supplementary Materials). So, remarkably, the noise terms $\xi_{i}$ are needed to access the community structure.

Using steps (i)-(iii) we obtained a model for the isolated dynamics, coupling function, and distribution of degrees. To apply current methods of community detection, we obtain an adjacency matrix with entries equal to zero or one from the matrix of correlations Fig. 12 (c). This is done by thresholding the entries and considering a link only when the correlations were greater than $0.5$. We tested threshold values ranging from $0.3$ to $0.6$ and obtained the same results. Indeed, the distribution of the entries of the matrix of correlations is unimodal and has a peak near 0.5. We computed rich-club coefficients for each node [42] (see details in Methods), and nodes with the highest coefficients formed the rich-club displayed at the centre of the network. By applying the community detection method [43] to the thresholded effective matrix, we recovered the communities according to their function as shown in Figure 2 (d).

There is no reason to expect that the links in the effective network correspond to links in the actual network. However, since every node makes most of its interactions within a cluster, two nodes with highly correlated fluctuations $\xi(t)$ are likely to belong to the same community, and this can be enforced in the effective network by adding a connection between them. Further results for a variety of rich-club motifs are shown in the Supplementary Materials.

Performance of the communities reconstruction.

To quantify the effectiveness of community reconstruction, we compute $PE=\frac{m}{N}$, where $N$ is the total number of nodes and $m$ is the number of nodes assigned to the wrong community. We compute $PE$ for different values for the coupling strength $\Delta\alpha$ between $0.05$ and $0.4$. For each value of $\alpha$, we considered 50 different simulations of the dynamics obtained by choosing different initial conditions. The figure shows the plot of the mean of $PE$ and a shaded region corresponding to the standard deviation. For $\Delta\alpha$ values larger than $0.4$, the reconstructing procedure cannot identify the communities correctly. This is due to the synchronization rich club which appears around this value.

Predicting critical transitions in rich-clubs

The ability to reconstruct the network and dynamics from data can be exploited to predict critical transitions that may occur when the coupling strength. In the cat brain, for example, a transition to collective dynamics in the rich-club has drastic repercussions for the functionality of the network [29, 44].

Our goal now is to obtain data when the network is far from a collective dynamics and predict the onset of collective motion in the rich-club. The effective network can predict the onset of such collective dynamics based on a single multivariate time series for fixed coupling strength in a regime far from the synchronized state. We analyze time-series obtained simulating the dynamics at a value of the coupling strength $\Delta\alpha=0.3$, and apply our reconstruction procedure to obtain the network structure and a model of the isolated dynamics.

Estimating the transition to synchronization in the rich-club. Transitions to synchronization between the bursts is possible while the fast spikes remain out of synchronization [45]. To estimate the transition to burst synchronization, we obtain the slow variable as a filter over the membrane potential (fast variable). Indeed, since we measure the membrane potential $y_{i}(t)=u_{i}(t)$, the slow variable is given as $z_{i}(t)=\mu\sum_{k=1}^{t}(y_{i}(k)-\sigma)$ and for a good choice of parameters $\mu$ and $\sigma$ this can be identified with the slow variable of the model $w$. In the methods section, we obtain an equation for the slow variable and analyze the effect of the network of the dynamics of the slow variable.

We can use the data on the network and the dynamics recovered from the time-series recorded at $\Delta\alpha=0.3$ to predict that at the value $\Delta\alpha\approx 0.42$ the rich-club will develop a burst synchronization (details in the Methods).

To capture such a transition in the bursts of the rich club, we introduce a phase $\theta_{j}(t)$ for the slow variable, the definition of $\theta_{j}(t)$ can be found in the Methods section. Once we have the phase we compute the order parameter

[TABLE]

Loosely speaking, a small value of the order parameter, $r\approx 0$, means that no collective state is present in the system, whereas $r(t)\approx 1$ means that the bursts are synchronized. Figure 4 shows that behaviour of the order parameter as a function of the coupling. The rich-club undergoes a transition to burst synchronization at $\Delta\alpha\approx 0.4$ that corresponds to an increase of roughly $40\%$ of the coupling strength and is close to the predicted value $\Delta\alpha\approx 0.42$. In the Supplementary Materials, we present other examples where the local dynamics is fully chaotic.

Obtaining a statistical description of the network

The effective network can also provide detailed information on the structure of each community and a statistical description of the network. To illustrate this, we reconstruct the statistical properties of scale-free networks.

Scale-free networks of coupled bursting neurons. We consider coupled bursting neurons with excitatory synapses [45] in scale-free networks (see Methods for details on simulation of scale-free networks and formulation of the model). Our techniques work equally well for spiking dynamics.

We generate a scale-free network with $N=10^{4}$ nodes such that the probability of having a node of degree $k$ is proportional to $k^{-\gamma}$, where $\gamma=2.53$. For this reconstruction we only need 2000 data points for each node. Again, to every pair of time series $y_{i}$ and $y_{j}$ we assign a Pearson distance $s_{ij}\geq 0$ and the node intensity $S_{i}=\sum_{j}s_{ij}$. The empirical distribution of the intensities $S_{i}$ approximates the degree distribution of the network, see Fig 5(d). In the example here, the estimated structural exponent from the distribution of $S_{i}$ is $\gamma_{\rm est}=3.1$, which yields a relative error of nearly 25% with respect to the true value of $\gamma$ (see also the plots in Figure 6 a)). The functional network therefore overestimates $\gamma$, which has drastic consequences for the predicted character of the network. For example, the number of connections of a hub for a scale-free network is concentrated at $k_{\rm max}\sim N^{1/(\gamma-1)}$, so the relative inaccuracy for the estimate $k_{\rm est}$ of the maximal degree is $k_{\rm max}/k_{\rm est}=N^{1/\gamma-1/\gamma_{\rm est}}$, which is about 500%. Such inaccuracy has important repercussions for the ability to predict the emergence of collective behaviour [46, 19].

The statistical measure used for the construction of functional network typically depend in a nonlinear way on the degrees, thus causing a distortion in the statistics. We will discuss the case of Pearson distance. Suppose that the signals $\{(y_{i}(t),y_{i}(t+1))\}$ are purely deterministic, $y_{i}(t+1)=g_{i}(y_{i}(t))$, and thus perfectly fall on the graphs of the functions $g_{i}$, that are determined by $\beta_{i}$ and so by the degree $k_{i}$ of the node. The Pearson distance $s_{ij}$ between the signal at $i$ and $j$ is going to be a number between 0 and 1, depending on how close these graphs are. Unfortunately this distance depends nonlinearly on the degrees $k_{i}$ and $k_{j}$. Devising another distance $s^{\prime}_{ij}$ without the knowledge of the interaction, in general, would still carry the nonlinear dependence on the degrees. Moreover, once fluctuations from the network are included as $y_{i}(t+1)=g_{i}(y_{i}(t))+\xi_{i}(t)$ the differences between time-series $\{(y_{i}(t),y_{i}(t+1))\}$ for different $i$s can be due to fluctuations rather than differences in the degrees. So the decomposition of the rules in terms of interactions and fluctuations is essential to recover precise details of the real degree distribution.

The effective network provides a better statistical description of the network structure. To compare with the functional network approach described above, we constructed an effective network of the same system tested for the functional network approach. The estimate for $\gamma$ from the effective network is $\gamma_{\rm est}=2.55$, which has an error of only $1\%$ ( Fig. 5 (d)). We repeat the analysis on a different network with different parameters $\gamma$ in the degree distribution. The estimated $\gamma_{\rm est}$ values are shown in Fig. 6 (c) as a function of the true parameter $\gamma$. The relative error on the estimated exponent is within $2\%$.

Real-world scale-free networks: the optic lobe of Drosophila melanogaster. Effective networks can be used to provide estimates for real-world systems. As one example, we applied our method to data simulated from the neuronal network in the D. melanogaster optic lobe, which constitutes $>$50% of the total brain volume and contains 1781 nodes [30]. The degree distribution has a power-law tail [31]. We used spiking neurons with chemical coupling to simulate the multivariate time series, from which we constructed an effective model capable of estimating the degree distribution with great accuracy (Fig. 6 (b) and Methods).

Additional applications of our procedure for different types of chaotic systems (doubling map, logistic maps, spiking neurons with electrical synapses, Rössler systems, and Hénon maps) can be found in the Supplementary Materials, together with results demonstrating that effective networks are robust against noise.

Conclusions

We have introduced an effective network obtained from time-series of a complex network (by observing the dynamics at each node). Our method complements the existing ones in two ways. Firstly, it encompasses the case of isolated chaotic dynamics. Secondly it deals with weak coupling among the nodes. Both cases are commonly found in applications. Key to the success of the reconstruction is the heterogeneity of the network which allows us to separate the contribution from the local dynamics and the interactions and thus perform a multi-level reduction.

We have compared our procedure with methodologies most relevant for the systems we considered. We have thus excluded results tailored to specific experimental setups or dynamics (binary dynamics coupled on networks [47], and see [10] for a review). We also did not consider methods that rely on measurements obtained intervening on the system with controlled inputs (see e.g. [13]) and restrict our attention to the case where the time-series is recorded under constant conditions. When the coupling is strong, sparse recovery can be applied [34]. On the other hand, when the coupling is weak sparse recovery cannot distinguish small parameters from those that are identically zero thus misidentifying connections between nodes as it can be observed in Figure 7. Also model-free methods are ill-suited to analyse weakly coupled chaotic dynamical systems studied here. Indeed, the influence of a single pairwise interaction on the time-series is so weak that its effect on the correlation of the measurements can hardly be detected by direct measurements.

Our reconstruction works when no collective dynamics is induced by the weak coupling. When the rich-club synchronizes. A synchronous rich-club induces correlation between the fluctuations in the communities. This happens because two nodes in different communities but coupled to the rich-club would receive a similar forcing. Since we use the fluctuation to identify the community structure, the method would identify these nodes as belonging to the same community. While our method works well for various isolated dynamical systems, in two scenarios it didn’t provide a good reconstruction. First, when the local dynamics is the Lorenz with classical parameters. In this case, the passage near a fixed points smears the fluctuations. Second, in the bursting dynamics, when the quiescent state is too long. This failure seems to be related as well to the long passage near a fixed point. Our method works with large networks, up to $\leavevmode\nobreak\ 10^{5}-10^{6}$ nodes after which the computation time can be long. In fact, the number of operations to fit the the dynamics at all nodes, as in point (i), is proportional to Number of nodes $\times$ (Length of the time series) β, where the power $\beta$ depends on the optimization algorithm adopted. The number of operations to determine the correlations between the fluctuations, point (iii) above, is of the order of (Number of nodes) 2 $\times$ Length of the time series. The effective network exploits the network heterogeneous structure. This allows to single out different dynamical behaviours across the network and use this information to reconstruct the local dynamics and the effective interaction function. To obtain the community structures we use that certain noise terms associated with the time series at two nodes in the same community are correlated. By collecting data when the network is far from critical transitions, our method can reconstruct the effective network and enables predictions regarding critical transitions occurring in a network.

Methods

Random scale-free networks

A scale-free network has degree distribution $P(k)=Ck^{-\gamma}$, where $\gamma>0$ is the characteristic exponent and $C$ is a normalising constant. We use a random network model $\mathcal{G}(\bm{w})$ which is an extension of the Erdös-Rényi model for random graphs with a general degree distribution [48]. Given ${\bm{w}}=(w_{1},w_{2},\cdots,w_{n})$, each potential edge between $i$ and $j$ is chosen with probability $p_{ij}=w_{i}w_{j}\rho,$ and where $\rho=(\sum_{i=1}^{n}w_{i})^{-1}.$ To ensure that $p_{ij}\leq 1$, we consider $\max_{i}w_{i}\rho\leq 1.$ Then $w_{i}$ is the expected degree of the $i$th node. Taking $w_{i}=ci^{-1/(\gamma-1)}$, where $\gamma\geq 2$ and $c=\frac{\gamma-2}{\gamma-1}dn^{1/\gamma-1}$ with $d$ denoting the mean degree, the distribution of expected degrees follows a power law, and so $P(w)\propto w^{-\gamma}$ is the expected exponent of power-law distribution [48, 49]. If the random graph generated is not connected, we use only the largest connected component (the giant component).

Finding power-law distribution parameters from empirical data. For estimation of the power-law distribution parameters, we use the maximum likelihood estimator first introduced by Muniruzzaman [50], which is equivalent to the well-known Hill estimator [51]. After that, we test the reliability between the data and the power law by using the goodness-of-fit method. If the resulting $p$-value is larger than 0.1, the power-law estimation is an appropriate hypothesis for the data. A complete procedure for the analysis of power-law data can be found in Ref. [52].

Functional networks

For networks composed of chaotic oscillators, building the functional network from the standard Pearson correlation between time series gives no meaningful results because of the decay of correlation intrinsic to the chaotic dynamics. Hence, functional networks are built using a Pearson distance $s_{ij}\geq 0$, which describes the proximity of the dynamics (i.e. the ‘attractor’) at two nodes $i$ and $j$. To do this, we consider the time series ${z}_{i}(t):=({y}_{i}(t),{y}_{i}(t+1))$, $t=0,\dots,T-1$ reordered in $z_{i}^{\rm lex}(t)$ according to the lexicon order; that is, according to the magnitude of the first component of $z_{i}(t)$. Then, let $r_{ij}$ be the Pearson correlation, $r_{ij}=$ Cor$(z_{i}^{\rm lex},z_{j}^{\rm lex})$, so that $r_{ij}=1$ indicates that the attractors at nodes $i$ and $j$ fully agree. Finally, define the Pearson distance $s_{ij}=1-|r_{ij}|$ so that $s_{ij}=0$ indicates full agreement of the dynamics and $s_{ij}>0$ measures the difference between the attractors.

The intensity $S_{i}=\sum_{j}s_{ij}$ approximates how many nodes have a dynamical rule different from $i$ and helps to distinguish between poorly connected nodes and hubs. Since most of the network is composed of poorly connected nodes with similar evolution rules, they exhibit a smaller $S_{i}$ than high-degree nodes, which are scarcer and have different dynamics from the low-degree nodes. In the case of a scale-free network where the degree distribution is monotonic, the empirical distribution of $S_{i}$ is also expected to have the same trend.

Effective network

Dimensional reduction over finite time scales. To construct an effective network, we combine statistical estimation with mathematical results for high-dimensional dynamical systems that are based on theorems describing the evolution of dynamical systems coupled on a heterogeneous networks.

The assumptions we make are: (a) the dynamics determined by $\bm{F}$ is ergodic, (b) the network is heterogeneous in the sense that most nodes are low degree and a few hubs make many connections, and (c) the coupling strength $\alpha$ is small so that $\alpha\max_{i}k_{i}=O(1)$. These assumptions play a role in the reduction process. Indeed, (c) means that the contribution of the low-degree nodes to the dynamics of the hub nodes is an average, (b) means that the ergodic behaviour of low-degree nodes is similar to that of $\bm{F}$, and finally, (a) means that we obtain a hierarchy of reductions with different dynamics (depending on connectivity) from which the structure of the network can be inferred.

In heterogeneous networks, and up to a small error (relative to the network size), the dynamics at each node evolves according to a low-dimensional dynamical system for a certain time scale where the effect of the interaction of one node on the entire network is averaged to a give net interaction. This is rigorously proven in [19] under some additional technical assumptions.

One can write the interaction term in equation (1) at node $i$ as

[TABLE]

where

[TABLE]

is the averaged interaction in terms of the invariant measure $\mu$ for the isolated dynamics and $\beta_{i}=\alpha k_{i}$ where $k_{i}$ is the degree of node $i$. The term $\bm{\xi}_{i}(t)$ is determined by the states of the $i$th node neighbours. The low-dimensional system, called the reduced system, plus a fluctuation term thus reads as

[TABLE]

where ${\bm{G}_{i}}(\bm{x})=\bm{F}_{i}(\bm{x})+\alpha k_{i}\bm{V}(\bm{x}),$ for $i=1,\dots,N$. Using the chaotic and ergodic properties of the dynamics and concentration inequalities, we can show that $|\alpha k_{i}\bm{\xi}_{i}(t)|$ is small over exponentially large time scales, in terms of the system size, and for most initial conditions. Not all structures of the graph will follow this reduced equation. If a fully connected graph is weakly coupled to a scale-free network, then it can develop dynamics different from those of the reduced model. In this case, we cannot use the physical measure $\mu$ in the expectation; however, the reduction is still valid whenever $\mu$ is substituted with a measure that satisfies self-consistency relationships (see below).

k-Fold Cross-Validation

Cross-validation is a resampling method used to evaluate dynamical models on a limited data sample. The parameter $k$ refers to the number of groups that a given data is to be divided into. We use $k=10$, since it is known that this value produces the lowest bias in several tests, so that the method can be called 10-fold cross-validation. The procedure is the following: First, randomly shuffle the data and divide it into k groups. Hold few groups to test the system, use the remaining ones as training data. In our case, we used 1 group to test the system and 9 groups to train our model. We fit a model on the training set and evaluate it on the test one. According to the evaluation scores, we summarize the dynamical properties of the model. Further details on how this standard procedure works can be found in [33].

Reconstruction of the effective network representation from data

Step 1 : Reduced dynamics. Given the multivariate time series $\{y_{1}(t),\dots,y_{N}(t)\}$, we consider each time series separately and then perform a Takens reconstruction of the attractor for each time series $y_{i}(t)$. The reduction guarantees with high probability that the dynamics is low dimensional. If the time series is high dimensional, we discard it, and otherwise, once we are in the appropriate dimension, we continue to the next step. Many cases can be captured by a one-dimensional map and Takens reconstruction yields the return map

[TABLE]

where $g_{i}:\mathbb{R}\rightarrow\mathbb{R}$. In fact, even when the local dynamics is three dimensional, by introducing a Poincarè section, we can sometimes reduce the analysis to one dimension (see Supplementary Materials for details).

Step 2: Isolated dynamics and effective coupling function. We first obtain a model for the isolated dynamics. Here we describe the procedure when the reduced dynamics can be modelled by a one-dimensional map, since the method works similarly in the higher-dimensional setting as shown in the Supplementary Materials.

To estimate the deterministic rule $g_{i}$ for each time series, we use a $10$-fold cross-validation with $90\%$ of the time series for training and $10\%$ for testing (see k-Fold Cross-Validation above). To obtain a model for the isolated dynamics we first build a functional network (see Functional Networks section in Methods). For the low-degree nodes, $\alpha k_{i}v$ is negligible and the dynamics at the low-degree nodes are close to $f$. We identify these nodes by analysing the distribution of $S_{i}$. We use the top $N_{\rm top}$ nodes of the highest intensity nodes to obtain a proxy for the isolated dynamics. We then average the rules obtained at those nodes to get $\langle g\rangle\approx f$. The choice of $N_{\rm top}$ is not fixed and depends upon the number of nodes and the fluctuation $\sigma_{g}^{2}=\langle(g_{i}-\langle g\rangle)^{2}\rangle$, which is averaged from the $N_{top}$ highest intensity nodes. A good heuristic to choose $N_{\rm top}$ is $\sigma_{g}^{2}/N^{1/2}_{\rm top}\ll 1$, as we have used here.

The effective coupling function. Since $\alpha k_{i}v=g_{i}-\langle g\rangle$, analysing the family $\{g_{i}-\langle g\rangle\}_{i=1}^{N}$ can yield the shape of $v$ up to a multiplicative constant via a nonlinear regression by imposing that $g_{i}-\langle g\rangle$ and $g_{j}-\langle g\rangle$ are linearly dependent. The choice of the base function for the fitting is supervised and depends on the particular application (see Supplementary Material for additional examples).

Step 3. Network Structural Statistics. After selecting a $v$ that satisfactorily approximates $g_{i}-\langle g\rangle$ up to a multiplicative constant over all indices $i$, we estimate the parameter $\beta_{i}$ using a dynamic Bayesian inference. Because the fluctuations $\xi_{i}(t)$ are close to Gaussian we use a Gaussian likelihood function and a Gaussian prior for the distribution of the values of $\beta_{i}$, and hence obtain equations for the mean and variance. We split the data into epochs of $200$ points and update the mean and variance iteratively.

Matching index: Recall that the degree of node $i$ is $k_{i}=\sum_{j}^{N}A_{ij}$ and count the number of neighbours it has. We now define the matching index of a graph [29]. Consider the neighbourhood of node $i$ as $\Gamma(i)=\{j\in\{1,\dots,N\}\,|\,A_{ij}=1\}$. This is the set of nodes that shares an edge with the node $i$. The matching index of nodes $i$ and $\ell$ is the cardinality of the overlap of their neighbourhoods $\mu_{i\ell}=|\Gamma(i)\cap\Gamma(\ell)|$. We are interested in the normalised matching index:

[TABLE]

or equivalently in terms of the adjacency matrix

[TABLE]

Clearly $\widehat{\mu}_{i\ell}=1$ if and only if $i$ and $l$ are connected to exactly the same nodes, i.e., $\Gamma(i)=\Gamma(\ell)$; whereas $\widehat{\mu}_{i\ell}=0$ if nodes $i$ and $\ell$ have no common neighbours.

Community structures. Once we have obtained the rules $g_{i}$, we filter the deterministic part of the time series $y_{i}$ and access the fluctuations $\xi_{i}$. We decompose the fluctuations $\xi_{i}=\xi_{i}^{\rm c}+\xi_{i}^{\rm o}$ where $\xi_{i}^{\rm c}$ is the fluctuation of the local mean field from nodes in the cluster containing $i$, and $\xi_{i}^{\rm o}$ is the contribution from outside the cluster. Since a node makes most of its connections within its cluster, $\xi_{i}^{c}\gg\xi_{i}^{o}$ with high probability, and thus if $i$ and $j$ belong to the same cluster $\mbox{Corr}(\xi_{i},\xi_{j})=\mbox{Corr}(\xi^{c}_{i},\xi^{c}_{j})$. Next, we notice that the common noise is generated by the common connections between nodes $i$ and $j$. In fact, for fixed isolated dynamics and coupling function

[TABLE]

It is well known that in the cat cerebral cortex nodes in the same community has a high matching index while nodes are distinct communities has a low matching index. In fact, this tends to be typically in modular networks [29]. Therefore, for nodes in distinct clusters the component $\xi_{i}^{\rm c}\approx 0$, so $\mbox{Corr}(\xi_{i},\xi_{j})\approx 0.$ Hence, we can recover the network structure by performing a covariance analysis of the noise.

Filtering out the deterministic part. Filtering out the deterministic part plays a major role in recovering community structures. Suppose we have two signals of the form $y_{i}(t)={Y}_{i}(t)+\zeta(t),\quad i=1,2$, where ${Y}_{i}$ is independent of $i$ and $\zeta(t)$ is a common noise term. $Y_{i}$ represents the superposition of the deterministic chaos and the independent fluctuations of the network setting. For the correlation, we have

[TABLE]

Hence, the large values of the variance of the time series ($\sigma_{y_{i}}\approx\sigma_{Y_{i}}\gg\sigma_{\zeta}$) suppress the contribution of the common noise, and an analysis solely based on the the original time series $y_{i}$ will overlook the common noise contribution. In fact, in the real system, we cannot assume ${Y}_{i}$ to be independently distributed. However, if the system is sufficiently chaotic, it will exhibit fast decay of correlations.

Generating the connectivity structure

To generate the adjacency matrix from a given degree distribution, we use the configuration model [53]. We know the number of stubs (half-links) for each node and the model randomly connects these stubs (half-links) to each other and generates an adjacency matrix that respects the given distribution.

Dynamical Models

We tested the effective network on a wide range of dynamics and networks satisfying assumptions (a)-(c). We focused on scale free networks and on rich club networks to evaluate the performance of the reconstruction. As a benchmark model, we restricted our attention to a neuron dynamics. Further models can be found in Supplementary Material.

Neuron Model

We use Rulkov maps to model spiking and bursting neurons [45]. In this setting, the variable $\bm{x}$ at each node is two dimensional and we will denote $\bm{x}=(u,w)$. The local uncoupled dynamics is $\bm{F}(\bm{x})=(F_{1}(u,w),F_{2}(u,w))$ with

[TABLE]

The variable $u$ represents the membrane potential of a neuron, and, in this case, $u$ is the state variable measured by the observed time series $y_{i}(t)$; that is, $\phi(\bm{x})=u$. Different combinations of parameters $\sigma$ and $\beta$ give rise to different dynamical states of the neuron, such as resting, tonic spiking, and chaotic bursts. To test our procedure we considered two cases where $\sigma=\nu=0.001$ and $\beta=5.9$, which correspond to tonic spiking and $\beta=4.4$ for bursting.

Chemical synapsis. Synaptic coupling occurs only through the membrane potential $u$. That is, $\bm{H}(\bm{x}_{i},\bm{x}_{j})=(h(u_{i},u_{j}),0)$ with $h(u_{i},u_{j})=(u_{i}-V_{s})\Gamma(u_{j})$, where

[TABLE]

and $V_{s}$ is a parameter called reverse potential. Choosing $V_{s}>u_{i}(t)$, the synaptic connection is excitatory. We take $V_{s}=20$, $\Theta_{s}=-0.25$, and $\lambda=10$. Applying the dimensional reduction in (3), we obtain a map with averaged value of the interactions $\bar{F}_{1}(u)=F_{1}(u)+\alpha k_{i}\langle\Gamma\rangle(u-V_{s}).$ One can a posteriori verify that $\langle\Gamma\rangle=\int\Gamma\mu(dx)$, where $\mu$ is the physical measure for the map $F_{1}$ and can be estimated from the dynamics of the low-degree nodes.

Electrical synapsis. Again, the synaptic coupling occurs only through the membrane potential $u$. That is, $\bm{H}(\bm{x}_{i},\bm{x}_{j})=(h(u_{i},u_{j}),0)$ with

[TABLE]

so the coupling depends only on the difference of states.

The cat cerebral cortex

The dataset for connections in the cat cerebral cortex consists of 53 geographic cortical areas investigated in detail with data-mining methods and was taken from [27, 28]. The network is organised into four functional areas: visual, auditory, somatosensory-motor, and frontolimbic.

Functional network. Several methods can be used to create a functional network to recover the clusters. The approach is to compute pairwise correlations of all nodes and to consider the similarities as interactions between nodes. This method is not appropriate for weakly coupled chaotic systems since they exhibit exponential decay of correlations.

Effective network representation. The connectivity matrix can be estimated using the correlation matrix $\rho_{ij}=\mbox{Corr}(\xi_{i},\xi_{j})$. The adjacency matrix $\bm{A}$ can be computed by thresholding the correlation matrix as $A_{ij}=\Theta(\rho_{ij}>\tau)$, where $\Theta$ is the Heaviside step function and $\tau=0.5$ is an arbitrary threshold. Note that we tested different thresholds with values from 0.3 to 0.6 and obtained similar results.

Finding communities. Given the reconstructed connectivity matrix $\bm{A}$, we used the modularity-based Louvain method [43] to detect communities.

Finding rich-club members. Colizza et al. developed an algorithm to detect members of the rich-club [42]. The algorithm provides a rich-club coefficient $\phi(k)\in[0,1]$ for each degree value in the network. We considered nodes with $\phi(k)>0.8$ to be the members of the rich-club.

Predicting critical transitions

Once we reconstruct the relevant information, we can perform numerical analysis of the recovered equations to explore the dynamics of the systems outside the observed range of parameters. Here we explain how to gather the information for a theoretical prediction of the critical transition.

Burst synchronization. We introduce the slow variable as discussed in the main body of the manuscript. To introduce a phase variable $\theta$ for the dynamics $y$ we first smooth the time series [54]. Then, we find the time $t_{n}$ of local maxima as the $n$th maximum point of the slow variable. We do this after performing a Gaussian filter to denoise slightly the data. Finally, we introduce the phase variable $\theta$ as

[TABLE]

as discussed in Refs. [55]

Reduction in the rich-club. Nodes in the rich-club have degrees of approximately $\Delta$ and make $\kappa\Delta$ connections inside the rich-club and $(1-\kappa)\Delta$ connections to the rest of the network. Following our reduction scheme, the interactions within and outside the rich-club can be described by the expected value of the interactions with respect to the invariant measure associated with each of them. Let $C$ denote the set of nodes in the rich-club, then the coupling term for the $i$th node in the rich-club is

[TABLE]

However,

[TABLE]

where $\mu$ is the invariant measure for the nodes outside the rich-club. Hence, for the rich-club we obtain

[TABLE]

where $\bm{q}_{i}(\bm{x}_{i}(t))=\bm{F}_{i}(\bm{x}_{i}(t))+(1-\kappa)\Delta\alpha\int\bm{H}(\bm{x}_{i},\bm{y})d\mu(\bm{y})$.

**Predicting the transition to collective behaviour ** Let us recall that when isolated $u_{i}(t+1)=F_{1,i}(u_{i}(t))+w_{i}(t)$ where $F_{1,i}\approx F_{1}$ and $w_{i}(t+1)=w_{i}(t)+\mu(w_{i}(t)-1)$. After some algebra we obtain

[TABLE]

Using the reduction Eq. (Predicting critical transitions), in the network we obtain

[TABLE]

where $i$ denotes the $i$th nodes in the rich-club, $\langle u\rangle$ is the mean membrane potential in the rich-club and $\xi_{i}$ are fluctuations. We fix two nodes $y_{i}=u_{i}$ and $y_{j}=u_{j}$ in the rich-club. Let us introduce the disturbance

[TABLE]

Applying the mean value theorem and using that $F_{1,i}\approx F_{1}$ we obtain

[TABLE]

and introducing a proxy for the dynamics of the slow variables

[TABLE]

and considering

[TABLE]

where we used that $\sum_{n=0}^{t}\zeta(n)$ is a slow variable. We obtain

[TABLE]

Recall that for the cat cerebral cortex $\Delta=37$. We given the time series $\{y_{i}\}$ for $\alpha\Delta=0.3$, we estimate $F_{1}$ using our method $i$ as the slow variables are constants over short time scales, and the obtain slow variables as a filter over the fast variables. Then, the parameters of the rich-club such as the mean number of connections among the rich-club in terms of the reduction and external connections in terms of strength of fluctuations (We will provide a throughly discussion about the estimation of these parameters for a fully chaotic system in the Methods). From the data and using Eq. (Predicting critical transitions) we estimate $\lambda=1.42$ and thus we obtain $\alpha\Delta=0.42.$ At this critical value the slow variables tend the stay together due to the contraction in the dynamics. This is related to the onset of synchronization in the bursts, which is captured via a phase variable through the order parameter.

Sparse Recovery. The method gives a way of recovering the equations from data by considering the evolution map as a linear combination of a well chosen basis of functions, called library. The main assumption is that many coefficients will be zero. That is, the vector of coefficients will be sparse [34]. Since we have $N$ nodes in our network, we consider the library $\mathcal{L}=[\phi_{1}(x_{1}),\phi_{1}(x_{2}),\cdots,\phi_{1}(x_{N}),\phi_{2}(x_{1},x_{N}),\cdots,\phi_{k}(x_{N},x_{N})]$ as the set of basis functions. And denote the $X=\left(x_{1}({t_{2}}),x_{2}({t_{2}}),\cdots,x_{N}(t_{1}),x_{1}(t_{2}),\cdots,x_{N}(t_{n})\right)^{*}$ where ∗ denotes the transpose and we introduce

[TABLE]

We then look for a solution of the system $X=\Theta\Xi$ where

[TABLE]

is the vector of coefficients. The sparse recovery technique then solves the linear equation for $\Xi$ iteratively enforcing the sparsity of $\Xi$ by introducing $\sigma$ such that if $|\zeta_{ij}|\leq\sigma$ we set such entry to zero [34].

In our case, we will consider the scenario where the synapsis between neurons are electric and we generate the multivariate data for the cat cerebral cortex. We will assume we have knowledge of the coupling function so we can easily read the network structure from the sparse recovery. We choose a library of polynomials as $\phi_{j}(x_{i})=x_{i}^{j}$ for $j=\{1,\cdots,5\}$ and the pairwise $\phi_{j}(x_{i},x_{k})$ to be homogeneous polynomials of degree two in the variables $x_{i}$ and $x_{k}$ for $i,k\in\{1,\cdots,N\}$.

We solve the minimization problem and group dependencies of the network through the coefficient vector as the network structure [34]. The strength of each connection is of order $\alpha\approx 0.015$. Hence we have chosen values of $\sigma$ close to this value. The reconstructed network does not identity the clusters correctly as can be seen by comparing the blue and red markers in Figure 7. Moreover, employing the technique of community detection shows that the probability of identifying a node in the correct cluster is around 45%. This happens as each single connection is not strong enough. In fact, since the dynamics of the neurons are chaotic each connection plays a role of noise which leads to the misidentification. Thus we have shown that the sparse recovery method is not appropriate in our setting.

Acknowledgements We would like to thank Tomislav Stankovski, Chiranjit Mitra, Mauro Copelli, Dmitry Turaev and Jeroen Lamb for enlightening discussions. This work was supported in part by the Center for Research in Mathematics Applied to Industry (FAPESP Cemeai grant 2013/07375-0), the European Research Council (ERC AdG grant number 339523 RGDD), and the Serrapilheira Institute.

Data Availability. The connection matrices of cat cortex can be found at https://sites.google.com/site/bctnet/datasets. Connectivity matrix of Drosophila Medulla can be found at https://neurodata.io/project/connectomes/.

1 Dimensional Reduction in Heterogeneous Networks

Using the notation in the main body of the paper, we present an informal statement of the theoretical result that suggests the reconstruction procedure. For a precise statement and details of the general setting see [1]. The theorem has three main assumptions concerning the uncoupled dynamics $\bm{F}$, the network encoded in its adjacency matrix $A_{ij}$, and the reduced dynamics $\bm{G}_{j}$. We discuss a particular case and use it to illustrate the reconstruction step-by-step.

The local dynamics must increase the distance between points by a constant factor. A paradigmatic one-dimensional example is $\bm{F}(x)=2x$ mod $1$ that presents some of the most interesting characteristics of chaotic maps.

• 2)

The networks are heterogeneous, having nodes with disparate degrees. This assumption is made precise in [1] by a set of conditions involving the size of the network $N$, the number of hub nodes $M$, the maximum degree of a hub node $\Delta$, and the maximum degree of a low degree node $\delta$. One instance in which these conditions are satisfied and that includes many situations of interest is the following. Suppose $M$ is constant or slowly increasing with the size of the network, $M\sim\log N$, $\Delta$ grows faster than the square root of the system size, $\Delta\sim N^{\frac{1}{2}+\epsilon}$, and $\delta$ increases slowly, but polynomially in the system size, $\delta\sim N^{\frac{\epsilon}{2}}$. Then for a sufficiently large $N$, the network satisfies the conditions.

• 3)

The reduced dynamics must be hyperbolic. This means that we assume the maps $\bm{G}_{j}$ to be either expanding (possibly by a non-constant factor), or to have a finite number of attracting periodic orbits. Every map can be perturbed by an arbitrarily small amount to obtain such an hyperbolic map [2].

Under these assumptions, we can prove the following approximation result

Theorem 1** ([1])**

For every hub node $j$, the dynamics at the hub is given by

[TABLE]

where $|\bm{\xi}_{j}(t)|<\xi$ for time $T$ with $1\leq T\leq\exp[C\xi^{2}\Delta]$, and a set of initial condition of measure $1-T/\exp[C\xi^{2}\Delta].$, where $C$ is constant in $\Delta$, and $\xi$.

Notice that one can pick the time scale $T$ exponentially large, but such that $T/\exp[C\xi^{2}\Delta]$ is very small so that, for large $\Delta$, the approximation result holds for very long time and for a large set of initial conditions. The reduced dynamics given by $\bm{G}_{j}$ depends on the hub’s state only. The term $\bm{\xi}(t)$ gives the fluctuations of the sum of all interactions from their expectation with respect to the physical measure of the unperturbed dynamics.

2 Reconstruction of degree distributions

Scale-free networks. We create ensembles of scale-free network of $N=6000$ nodes with distinct structural exponents $\gamma$ following the same technique as in the main body of the manuscript. As in the main body, we only consider the largest connected component of the network (giant component). For each network realisation we compute the largest degree $\Delta$ and, for simplicity, we fix the coupling strength $\alpha/\Delta=0.5$ throughout the rest of the section.

2.1 Doubling map

Let us now apply our approach when $\bm{F}_{i}$ is a perturbed version of the doubling map. Since the dynamics is one dimensional, we denote $\bm{x}=x$ and $\bm{F}_{i}(\bm{x})=f_{i}(x)$ with

[TABLE]

and where we take $\varepsilon_{i}$ to be i.i.d. random variables (i.e. independent over $i$) uniformly distributed on $[0,10^{-3}]$. Likewise we write $\bm{H}=h$ with

[TABLE]

We fixed $\alpha=10^{-2}.$ Since the unique absolutely continuous invariant measure for the doubling map is the Lebesgue measure $m$, we have

[TABLE]

yielding $v(x)=-\sin 2\pi x$ and the reduced dynamics takes the form

[TABLE]

We aim to recover the reduced dynamics and in particular $f$ and $v$ from the data of a multivariate time series with $T=2000$ time steps.

From data to Model. We assume not to have access to the network structure and only measure the time-series $\{x_{i}(t)\}$ at each node, as illustrated in Figure 8 Data. By performing the steps described in the main body of the manuscript, we can recover an effective network. We now illustrate in detail our procedure in the case above of the doubling map.

(i) Reduced dynamics. From the time series observed at each node, we construct the attractor for that particular node. Computing the embedding dimension we notice that the dynamic at each site (node) is well described by a one dimensional map, $g_{i}(x)=f_{i}(x)-\alpha k_{i}v(x)$, up to a fluctuation $\xi_{i}$. Hence, to reconstruct the attractor it suffices to obtain the return map. Return maps at different nodes are shown in Figure 8 (i).

(ii) Isolated dynamics and effective coupling. We start by introducing a similarity measure.

Similarity of the time series. Two time series are considered almost the same, if whenever $y_{i}(t)$ and $y_{j}(t^{\prime})$ are close then $y_{i}(t+1)$ and $y_{j}(t^{\prime}+1)$ are also close. This means that one considers the new times series ${z}_{i}(t):=({y}_{i}(t),{y}_{i}(t+1))$, $z_{j}(t):=({y}_{j}(t),{y}_{j}(t+1))$, $t=0,\dots,T-1$, reshuffle them according to the lexicographical ordering (i.e. according to the magnitude of the first coordinate), and then take $s_{ij}$ to be the Pearson correlation distance between these reshuffled sequences. So if the time series at $i$ and $j$ are just out of phase, then, thanks to the reordering, their distance $s_{ij}=0$. We then calculate the intensities $S_{i}=\sum_{j}s_{ij}$ and use this to classify the nodes. The similarity matrix is shown in Figure 8(ii).

Using the correlation distance described above, we determine for which nodes $S_{i}$ is minimal and in this way we identify the low degree nodes. For the low degree nodes, the sequence $\{(y_{i}(t),y_{i}(t+1))\}_{t=0,\dots,T-1}$ lies close to the graph of the return map for $f$, and thus we obtain a good estimate for the isolated dynamics $f$, see Figure 8(ii). Next, we apply the reduction to estimate the coupling function.

To obtain the effective coupling $v$ from data, take a hub $j$ and consider the sequence $\{(y_{j}(t),y_{j}(t+1))\}_{t=1}^{2000}$. The resulting time series approximates the graph of $g_{i}$ and subtracting $f$ gives an approximation of the function $v$ up to a multiplicative constant (shown in Figure 8(ii)).

(iii) Network Structural Statistics. There are two ways to obtain information about the statistics of the degrees $k_{i}$. The first one uses the noise variance. In fact, the size of the fluctuations $\xi_{i}$ depends on the number of connections the node $i$ makes, and with good approximation Var$(\xi_{i})\propto k_{i}$. The second one uses $f$ and $v$ recovered at the previous step. For every node $i$, choose $\beta_{i}$ to fit the time series $\{(y_{i}(t),y_{i}(t+1))\}_{t=0,\dots,T-1}$ with the map

[TABLE]

The value of $\beta_{i}$ is obtained by Bayesian inference on the return maps is shown in the first panel of Figure 8(i) as $f$ and $v$. The distribution of $\beta_{i}$ will be linear proportional to the degree distribution and so we can use it to obtain the the structural parameter. At this point we can construct an effective network with degree distribution close to that of the real network can be constructed with the configuration model as discussed in the main body of the manuscript. We can also check whether there are communities in the network by analyzing the covariance of the fluctuations $\xi_{i}$.

In Figure 10 a) we reconstruct from the distribution of $\beta_{i}$ the structural exponent $\gamma$ for 1000 scale-free network with exponents ranging from $2.4$ to $3.6$

2.1.1 Robustness of the reconstruction under noise

Adding some small independent noise to the dynamics does not influence much the reconstruction procedure for the doubling map. This is a consequence of stochastic stability of the local dynamics together with the persistence of the reduction [1]. On the other hand, when the fluctuations become too large the reconstruction will underestimate the network structure. To illustrate these effects, we consider the randomly perturbed doubling maps

[TABLE]

where the random variables $\eta_{i}(t)$ are independent over $i$ and $t$, and identically distributed uniformly in the interval $[-\eta_{0},\eta_{0}]$. Intuitively, as long as $\eta_{0}<\alpha\min k_{i}$ the reconstruction will go through as the noise fluctuation will not compete with the coupling term. Notice that we normalize $\|v\|=1$. This is illustrated in Figure 9. In inset a) the noise has a large support $\varepsilon=0.1$, and in particular larger than the coupling $\alpha\min k_{i}=10^{-2}$. As a consequence, the reconstruction understimate the number of low degree nodes. In inset b) the noise has a small support $\varepsilon=10^{-3}$.

2.2 Logistic Maps on Scale-Free Networks

Consider $M=[0,1]$, $\bm{x}=x$ and $\bm{F}_{i}(\bm{x})=f_{i}(x)$ where

[TABLE]

It is known that for such parameter, the map $f$ has an absolutely continuous invariant probability measure. In contrast to the doubling map, a small perturbation of the logistic map might result in a map that has a non-smooth physical measure, namely the measure is a distribution given by the sum of Dirac delta masses at a finite number of points that constitute the periodic orbits of the system. This compromises the validity of the rigorous result analysis that we used in the case of the doubling map, for which small perturbations produced small changes in the invariant measure [2]. Nonetheless, as we will show, the reconstruction analysis gives good results in this case as well. We believe that this happens because our reconstruction only needs a finite number of points and gives a representation of the dynamics over finite time-scales. As a coupling function, we consider

[TABLE]

and the reduction is given by $x_{i}(t+1)=f(x_{i}(t))+\alpha k_{i}v(x_{i}(t)),$ where $v$ is defined as above. Again, we consider the reduced model $g_{i}(x)=f(x)+\beta_{i}v(x)$ in terms of the free parameter $\beta_{i}=\alpha k_{i}$. Analysing the distribution of $\beta_{i}$ we access the distribution of the degrees. Results are presented in Figure10b).

2.3 Spiking Neurons with Electrical Synapses

We use the same models for the local dynamics of the spiking neurons which are described in the section Methods. We consider here a different coupling function. Denoting $\bm{x}=(u,w)$, the coupling is

[TABLE]

Results of the reconstruction are presented in Figure 10c).

2.4 Henon Maps on Scale-Free Networks

Using the same notation, $\bm{x}=(u,w)$, the coupled Hénon maps are given by

[TABLE]

We observe the dynamics of the first component of the Hénon map, that is, $y=\phi(\bm{x})=u$. In this case, the reconstruction starts by determining the dimension of the reduced system. Takens embedding reveals that the dimension is two for large time excursions. Hence, we aim at learning a function of the kind

[TABLE]

We use polynomial functions for the fitting and perform a 10-fold cross-validation. Just our theory implies that

[TABLE]

where $f$ models the isolated dynamics and $h$ the coupling. Again we obtain $f$ from the low-degree nodes via a similarity analysis. We learn the function $h$ by

[TABLE]

In our case, $v(y_{i}(t),y_{i}(t-1))\propto y_{i}(t-1)$ and we can obtain the degree distributions, Fig. 10d).

2.5 Coupled Roessler Oscillators

Assume that the local dynamics is modelled by a Rössler oscillator [3]. The dynamics is now in continuous time and our method can also be applied by using a suitable Poincaré section. This gives an induced map that describes the dynamics of the system at specific instants of time (when the system hits a selected subset of phase space). Denoting $\bm{x}=(x,y,z)^{*}$ the vector field is given by $\bm{F}(\bm{x})=(y-z,x+0.2y,0.2+z(x-9))$ and the coupling function, assumed to be diffusive, is given by $\bm{H}(\bm{x}_{i},\bm{x}_{j})=\bm{E}(\bm{x}_{j}-\bm{x}_{i})$, where $\bm{E}$ projects to the first component, i.e., $\bm{E}(x,y,z)=(x,0,0)$. So, our main equation reads as

[TABLE]

We perform a numerical integration of the equations on the Rich-Club network using a 4th order Runge-Kutta with integration step $10^{-4}$ and get the data $\{\bm{x}_{i}(t)\}_{t\geq 0}$. Using a statistical analysis of the time-series of the state variables, we are not able to reveal the connectivity structure.

The data is phase coherent, that is, taking a Hilbert Transform we can decompose the time series in terms of amplitude and phase we conclude that the spread in the phase variable is small and thus the return time to a given section is nearly constant. So, we consider the Poincaré section defined by the maxima $w_{i}$ of the time series $x_{i}(t)$. This gives us a time series $\{w_{i}(n)\}$ indexing all maxima. We then apply all the steps of the reconstruction procedure to this time series. Because of the coherent dynamics of the phase, the coupling form is preserved in the Poincaré section. The results of the network structure estimation are presented in Figure 11.

3 Reconstruction of Rich-Club Mofits

We report the performance of the method in the setting of a network of $100$ nodes having five clusters of $20$ nodes each. Four of these clusters are modelled as Erdös-Renyi random graph with connection probability $p=0.3$. The remaining cluster is the integrating clusters with connection probability $p=0.8$. The network resembles Fig. 12.

Filtering out the deterministic chaos. We need to filter the contributions of the deterministic parts to reconstruct the community structure. Indeed, for two nodes $i$ and $j$ in the same cluster, the signals have the form $x_{i}(t)=\tilde{X}_{i}(t)+\zeta(t),$ and $x_{j}(t)=\tilde{X}_{j}(t)+\zeta(t)$ where $\tilde{X}_{i}$ and $\tilde{X}_{j}$ is a superposition of the deterministic chaos depending on the variable at the node and independent fluctuations coming from the rest of the network, while $\zeta$ is common noise. $\tilde{X}_{i}$ and $\tilde{X}_{j}$ have fast decay of correlations, depends on different sets of variables, and for the sake of the following argument, can be assumed to be independent between each other and with $\zeta$. Under these assumptions, $\mbox{Corr}(\tilde{X}_{i}(t)+\zeta(t),\tilde{X}_{j}(t)+\zeta(t))=\mbox{Var}(\zeta)/\sigma_{x_{i}}\sigma_{x_{j}}$. Hence, the large values of the variance of the time series leads to strong suppression of the correlation coming from the small common noise $\zeta$.

3.1 Condition for recovering the community

In general, the coupling function is a sum of terms $h(x,y)=u(x)v(y).$ This leads to noise terms

[TABLE]

where $\mu$ is the physical measure of the local dynamics. Given $i$ and $j$ the sum can be split into connections commons the $i$ and $j$ and to the independent connections. So, for such $i$ and $j$ we can write

[TABLE]

where $w$ is the noise due to the common connections. Notice that $w$ has zero mean. Let’s estimate the covariance of the component $\xi_{i}$ and $\xi_{j}$. By abuse of notation we will omit the time index $t$ and write the covariance and by the previous computations $\mbox{Cov}(\xi_{i},\xi_{i})=\mathbb{E}[(u(x_{i})w)(u(x_{j})w)].$ After some manipulations, we obtain

[TABLE]

so, if $\int u(x)d\mu(x)=0,$ the correlation between the noise will vanish even though they have a common term. Thus, the generic condition in order for the above scheme to be able to recover communities is that $\langle v\rangle\not=0$. If this condition is not met, the network reconstruction via the $g_{i}$’s is also not possible. We recall that $\langle v\rangle=0$ is not a generic condition and this is why the effective network approach works in most cases.

3.2 Doubling Map on Rich Club Network

Using the doubling map with the diffusive coupling described in the above section on the reconstruction for scale-free networks, we test the community detection from time-series using the correlation between time series and correlations between the noise $\xi_{i}$. Then we apply our technique to recover the network structure. We show the results in Fig. 13

3.3 Logistic Maps on a Rich Club Network

The dynamic of the logistic map and the coupling are described in the above section on the reconstruction for scale-free networks. We test the community detection approach on time-series using the correlation between time series and correlations between the noise. To this end, we fix the coupling strength $\alpha=10^{-4}$ and simulate the network dynamics. Then we apply our technique to recover the network structure. We show the results in Fig. 14.

3.4 Spiking neurons coupled with Electric Synapsis on a Rich-Club

Again, we consider the spiking neurons and the electrical coupling described above on the same rich-club motif as before. The reconstruction and analysis of the noise correlation is able to detect the rich-club clusters. To this end we fix the coupling strength $\alpha=5\times 10^{-4}$ and simulate the network dynamics. Then we apply our technique to recover the network structure. We show the results in Fig. 15.

3.5 Bursting neurons coupled with Electric Synapsis on a Rich-Club

Our technique applies equally well then the local dynamics has multiple time-scales such as a bursting neuron. Our extensive numerical investigation reveals that when the resting time is not much larger then the total bursting time the reduced dynamics is capable of extracting the relevant information of the time series. Thus, we fixed the neuron parameter $\beta=4.4$ to obtain a bursting dynamics. The reconstruction and analysis of the noise correlation is able to detect the rich-club clusters. To this end we fix the coupling strength $\alpha=10^{-3}$ and simulate the network dynamics. Then we apply our technique to recover the network structure. We show the results in Fig. 16

4 Predicting Critical Transitions in Rich-Club Networks

We present another example of how to use the effective network methodology to predict critical transitions. In this case we choose the doubling map for the local dynamics coupled on a rich-club network with clusters sampled as Erdös-Renyi graphs. Consider a rich club network of 2200 elements with 5 clusters. Four of the clusters are made of $N_{\ell}=500$ nodes with small degrees, and one cluster (called integrating cluster or rich club) has $N_{I}=200$ nodes which are connected with most of the network. The edges within a cluster of low degree nodes are assigned as described above.

As a model for the isolated dynamics, we use the doubling map $f(x)=2x$ mod 1 with the diffusive coupling $\bm{H}(\bm{x}_{i},\bm{x}_{j})=h(x_{j},x_{i})=\varphi(x_{j})-\varphi(x_{i})$ where we picked $\varphi(x)=\sin(2\pi x)$.

We use the function

[TABLE]

as an empirical measure of the synchronization level at time $t$. From a single multivariate time-series, when the coupling is fixed at $\alpha_{0}\Delta=0.2$ (red marker Figure 17), the analysis of $\langle E\rangle$ gives no sign of critical transitions. A statistical analysis shows that the variance of the averaged synchronization error is not amplified, and the extreme value statistics fails to reveal a transition as it also does an analysis of dynamical correlations.

As described above, we can use the effective network methodology to recover the local dynamics $f$ and the effective coupling $v$. In particular we recover the isolated dynamics $f(x)=2x$ within $4\%$ accuracy for the poorly connected nodes. The unique physical equilibrium measure for this map is the uniform distribution $m$ on $[0,1)$. Furthermore, we recover the effective coupling

[TABLE]

At this point, we know that the approximate evolution of any node $i$ is given by $g_{i}=f_{i}+\alpha k_{i}v_{i}$, and we can study numerically or analitically this rule to find those values of $\alpha$ for which the integrating cluster exhibits synchrony. The analysis shows that excursions towards synchronization will start once the coupling is increased by $15\%$ (see Figure 17). Below we provide the details on how to recover this critical value of $\alpha$ for which a transition arises.

Reduction in the Integrating Cluster. Nodes in the integrating cluster have roughly degree $\Delta$ and make $\kappa\Delta$ connections inside the integrating cluster and $(1-\kappa)\Delta$ to the rest of the network.The interactions felt by a node in the rich-club can be split in those coming from nodes in other clusters and those coming from nodes within rich-club itself:

[TABLE]

But,

[TABLE]

where $\mu$ is the invariant measure for the nodes outside the integrating cluster 111In the example above $\mu$ is the Lebesgue measure $m$.. Hence, the equation of the integrating cluster can be written as

[TABLE]

where $q_{i}(x_{i}(t))=f_{i}(x_{i}(t))+(1-\kappa)\Delta\alpha\int h(x_{i},y)d\mu(y)$.

We can estimate $\mu$ empirically analyzing each cluster. Assuming that $h(x,y)=\varphi(y)-\varphi(x)$, we can recover $\varphi$ from the analysis of $q_{i}$ using the reconstruction techniques. In fact, for this particular choice of $h$, $\int h(x,y)d\mu(y)$ is equal to $-\varphi(x)$ plus a constant. Now if $\nu$ is the measure that describes the behaviour of a node in the integrating cluster, then the interaction within the rich-club can be written as

[TABLE]

Putting the two equations together one has that

[TABLE]

where $g_{i}=f_{i}-\alpha\Delta\varphi$ models the reduced dynamics and the $c(\mu,\nu)=(1-\kappa)\Delta\alpha\int\varphi(s)d\mu(s)+\kappa\Delta\alpha\int\varphi(s)d\nu(s)$ is the mean contribution from all interactions (inside and outside the integrating cluster). Finally $\zeta_{i}=\xi^{o}_{i}(t)+\xi^{c}(t)$ combines the effect of the fluctuations.

Model from Data. From a single multivariate time series at a given coupling parameter, shown in Figure 17 as a red dot, we reconstruct the model. First, we obtain the the rule $g_{i}$ which we uncover to be $g_{i}(x)=2x-\beta_{i}\sin 2\pi x$ mod $1$, hence $\beta_{i}=0.169$ and this number is nearly independent of the node in the integrating cluster. Hence, we obtain an estimate $(\alpha\Delta)_{\rm est}=0.169$. We also obtain that the dynamics of nodes in the communities is well approximated by $f(x)=2x$ mod $1$.

Next, we need to estimate $\kappa$ to construct a model for the connectivity of the integrating cluster. From the recovered local dynamics $f$, one knows that $\mu$ is the Lebesgue measure and, since $\int\varphi d\mu=0$, we obtain $c(\mu,\nu)=\kappa\int\varphi d\nu$. In the regime of parameters where the measurements have been made, $\nu$ can be obtained empirically and this allows to recover $\kappa$ since

[TABLE]

where we evaluate $\int\phi d\nu$ with respect to the empirically retrieved $\nu$, we substitute $\alpha\Delta$ with the estimated value above, and $\langle\cdot\rangle$ denotes the time average. We obtain that $\kappa=0.47$. The data analysis reveals that such $\kappa$ value is nearly independent of the node in the integrating cluster. Moreover, the analysis of the covariance of the noise in the integrating clusters shows a lack of communities and the functional analysis network indicates that the integrating cluster of is a random network 200 nodes with $p=0.83$. From this analysis, we obtain an estimate for $\Delta\approx p\times 200/\kappa=353=\Delta_{est}$. We are now able to estimate $\alpha$ by $\alpha_{est}=(\alpha\Delta)_{est}/\Delta_{est}=5\times 10^{-4}$.

Synchronization Prediction. With the reconstructed data $(f,v,A)$, we obtain that for a node in the integrating cluster, (11) reads as

[TABLE]

When $1/2\pi<\alpha\Delta<3/2\pi$, the map $2x_{i}(t)-\alpha\Delta\sin(2\pi x_{i}(t))$ mod $1$ has an attracting fixed point at 0. In this new regime, $\nu=\delta_{0}$ is a self-consistent measure for the integrating cluster, meaning that the measure $\nu$ gives rise to an approximated dynamical rule $g_{i}$ that has $\nu$ as equilibrium measure. This can be easily verified since $\int\sin(2\pi s)d\delta_{0}(s)=0$. We can then conclude that by selecting a value of the coupling strength such that $\alpha\Delta$ is in the range above, one expects all the states at the nodes in the integrating cluster to evolve towards the point [math] and fluctuate around this point by $\xi(t)$.

To obtain the range of $\alpha\Delta$ such that the fluctuations of $E$ are around twice $\max\|\xi\|$, we notice that since $\nu=\delta_{0}$ the stability properties are given by the linear stability around $x=0$. Let $u$ be a small displacement around [math] and let us denote at $J(\alpha)=Dg_{i}(0)$ the Jacobian of the map $g_{i}$ at [math]. Then we obtain that $u(t+1)=\sum_{i=0}^{t}J(\alpha)^{i}\xi_{i}$ and so $\|u\|\leq\max\|\xi\|/(1-J(\alpha))$, which yields that $1/(1-J(\alpha))<2$. Hence, we obtain $3/4\pi<\alpha\Delta<5/4\pi$. Because we measured $\alpha\Delta$ as $0.169$ in the data given we predict that a high quality coherent state in the rich-clubwill appear then $\alpha\Delta$ is increased by $40\%$. This is a agreement with the experiments.

5 Cat cerebral cortex

5.1 Bursting neurons with Chemical Synapesis

We simulated each mesoregion of the cat cerebral cortex network with bursting Rulkov oscillators coupled through chemical synapsis. For such dynamics the parameters are given as $\beta=4.4$ and $\Delta\alpha=0.05$ where there is no synchrony between oscillators Fig. 18. We use the simulated data to reconstruct the network structure as shown in Figure 18.