Search tips
Search criteria 


Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2017; 7: 44037.
Published online 2017 March 10. doi:  10.1038/srep44037
PMCID: PMC5345088

Multiplex visibility graphs to investigate recurrent neural network dynamics


A recurrent neural network (RNN) is a universal approximator of dynamical systems, whose performance often depends on sensitive hyperparameters. Tuning them properly may be difficult and, typically, based on a trial-and-error approach. In this work, we adopt a graph-based framework to interpret and characterize internal dynamics of a class of RNNs called echo state networks (ESNs). We design principled unsupervised methods to derive hyperparameters configurations yielding maximal ESN performance, expressed in terms of prediction error and memory capacity. In particular, we propose to model time series generated by each neuron activations with a horizontal visibility graph, whose topological properties have been shown to be related to the underlying system dynamics. Successively, horizontal visibility graphs associated with all neurons become layers of a larger structure called a multiplex. We show that topological properties of such a multiplex reflect important features of ESN dynamics that can be used to guide the tuning of its hyperparamers. Results obtained on several benchmarks and a real-world dataset of telephone call data records show the effectiveness of the proposed methods.

A current research trend aims at investigating complex time-variant systems through graph theory, by considering suitable features associated with vertices and edges1. Of particular interest are those systems that also perform a computation when driven by an external input signal. An example is that of artificial RNNs2,3,4, which are computational dynamical systems whose link with physics and neurosciences dates back to the ‘80 with some pioneering works from Jordan5 and Amit et al.6. Nowadays, RNNs are gaining renewed interest in neuroscience due to their biological plausibility7,8,9,10 and in computer science and engineering for their modeling ability11,12. RNNs are capable to generate complex dynamics and perform inference based on current inputs and internal state, the latter maintaining a vanishing memory of past inputs13,14.

Let us consider trajectories describing the evolution of a dynamical system in state space, e.g., the space containing all possible system states. As an example, in Fig. 1 we show the trajectories of a dynamical system operating in ordered (left) and chaotic (right) regimes, whose state is defined by the values of variables θ1, θ2, and θ3 at time t. Discriminating between order and chaos is of fundamental importance for investigating the properties of a dynamical system. It emerges that, such properties are manifested in the memory of the dynamical system and in the divergence rate of its state trajectories15. Fading memory is a desirable property of a dynamical system, and is characterized by ordered (contractive) dynamics. This is also referred to as the echo state property in the reservoir computing community16 and ensures that the current state/output of the system depends only on a finite number of past states/inputs17. At the same time, a high divergence rate between state trajectories (a property of chaotic dynamics) is also a desirable feature of RNNs. Results show that RNNs operating in a chaotic regime are able to produce meaningful patterns of activity18 and a balance must be struck in order to meet both properties. As a consequence, a (computational) dynamical system has to operate on the transition between order and chaos, in a region of the controlling parameter space called “edge of criticality”. On the edge of criticality, RNN internal dynamics becomes richer, meaning that neuron activations become heterogeneous15. Such a diversity both improves memory capacity of RNNs as well as their capability to reproduce complex target dynamics19,20. The notion of edge of criticality permeates several complex systems21,22,23, including random Boolean networks24 and populations of spiking neurons25. It is furthermore well-known that several complex systems spontaneously adapt to operate towards the edge of criticality, according to a mechanism known as self-organizing criticality in the statistical physics literature26.

Figure 1
Trajectories in the state space of a system described by the evolution of variables θ1, θ2, and θ3 over time.

In this work, we study ESNs27,28, a class of RNNs characterized by a large, untrained recurrent layer producing outputs by means of a linear (trained) readout layer. As such, ESNs are trained by optimizing readout weights only through a fast linear regression, rather than learning the weights of the recurrent layer by means of backpropagation through time, as in more general RNNs11. Please, note that more sophisticated readout layers, requiring longer training procedures, could be used as well29. ESNs trade the precision of gradient descent with the “brute force” redundancy of the random recurrent layer. This, inevitably, makes ESNs more sensitive to selection of the hyperparameters that control their internal dynamics. Therefore, the difficulty involved in finding the optimal hyperparameters is more crucial for ESNs than in other RNN models. The non-linear recurrent nature of ESNs hampers interpretability of internal state dynamics30. As such, ESNs (as any other RNN) are usually treated as black box models with hyperparameters tuned in a supervised manner through cross-validation. However, cross-validation requires evaluation of the computational performance on a validation set for each hyperparameter configuration. This might be a serious limitation in real-life applications, whenever the data are scarce and/or supervised information unavailable. Moreover, the model has to be re-trained for each hyperparameter configuration; this is an issue whenever the learning procedure is computationally demanding29.

Several unsupervised methods for ESNs tuning have been proposed in the literature, e.g., see refs 31,32. Although unsupervised tuning seldom achieves the same level of accuracy of the supervised counterpart, it offers some insights to the interpretation of the system behavior33. Typically, unsupervised learning is achieved by observing the dynamics of the states. These methods are usually based on a statistical framework, which requires reservoir outputs to be independent and identically distributed (i.i.d.). However, internal states are not identically distributed when ESNs are driven by non-stationary inputs34. Independence, instead, is always violated both by dependent inputs and by temporal dependencies introduced by recurrent connections. Therefore, even if the input signal is i.i.d., neurons activations may be identically distributed but never independent.

With this work, we propose two unsupervised graph-based methods for tuning ESN hyperparameters, with the aim to maximize (i) prediction accuracy and (ii) memory capacity. In our framework, temporal dependencies in time series of neurons activations (states) are converted into connections of a graph representation. This step allows to relax the i.i.d. assumption of statistical methods. Then, graphs associated with all neurons activations become layers of a structure called a multiplex. A multiplex graph35,36,37,38,39 is a special type of multilayer network, whose vertices are replicated through each layer and connected across layers only with their replicas. The topology (i.e., the edges) of the graph in each layer can be different. In particular, here we use a multiplex visibility graph40, whose layers consist of a particular class of graphs called horizontal visibility graphs (HVGs)41. HVGs are planar graphs built from real-valued time series, whose elements have one-to-one correspondence with the vertices. Important properties, linking the structure of HVGs with features (e.g., onset of chaotic behavior) of the dynamic system underlying the analyzed time series, have been recently studied42,43,44,45,46. This last aspect provides an important justification for using HVGs to model neuron activations for the purpose of hyperparameter tuning.

The novelty introduced in our paper refers to the characterization of ESN dynamics through structural characteristics of the associated multiplex network. To the best of our knowledge, this is the first time ESNs are studied within the framework of multiplex networks. However, it is worth citing a loosely related paper by Zhang et al.47, where the authors set the basis for a graph-theoretical analysis of RNNs. Here, we represent the instantaneous state of an ESN through a set of vertex properties of the multiplex (e.g., degree, clustering coefficient). Edges in HVGs might cover a relevant time interval (contained within the longest period of the signal), while they are still local in terms of topology. Accordingly, the ESN state can be characterized also by information that is non-local in time. To find hyperparameters yielding highest prediction accuracy, we search hyperparameter configurations producing neuron activations as diverse as possible. This occurs when neurons dynamics is maximally heterogeneous (critical dynamics), a characteristic that, as we will show in the paper, is well-captured by the average entropy of vertex properties of the multiplex. Successively, to quantify the amount of memory in an ESN, we check the existence of neuron activations that are “similar” to different delayed versions of the input. In fact, memory in ESNs depends on the ability to reproduce past input sequences from information kept within some neuron activations. We describe dynamics of delayed inputs and of neuron activations through unsupervised graph-based measures. Then, by evaluating the agreement between such measurements, we quantify the memory capacity. We provide experimental evidence that our methods achieve performance comparable with supervised techniques for identifying hyperparamer configurations with high prediction accuracy and large memory capacity.


Echo state networks

A standard ESN architecture consists of a large, recurrent layer of non-linear neurons sparsely interconnected by edges with randomly generated weights, called a reservoir, and a linear, feedforward readout layer that is usually trained with regularized least-squares optimization27. De facto, the recurrent layer acts as a non-linear kernel48 that maps inputs to a high-dimensional space. The time-invariant difference equations describing the ESN state-update and output are, respectively, defined as

An external file that holds a picture, illustration, etc.
Object name is srep44037-m1.jpg
An external file that holds a picture, illustration, etc.
Object name is srep44037-m2.jpg

The reservoir contains Nr neurons whose transfer/activation function f(·) is usually a hyperbolic tangent and g(·) is usually the identity function. At time instant t, the ESN is driven by the input signal An external file that holds a picture, illustration, etc.
Object name is srep44037-m3.jpg and produces the output An external file that holds a picture, illustration, etc.
Object name is srep44037-m4.jpg, being Ni and No input and output dimensions, respectively. The vector An external file that holds a picture, illustration, etc.
Object name is srep44037-m5.jpg describes the ESN (instantaneous) state. The weight matrices An external file that holds a picture, illustration, etc.
Object name is srep44037-m6.jpg (reservoir connections), An external file that holds a picture, illustration, etc.
Object name is srep44037-m7.jpg (input-to-reservoir connections), and An external file that holds a picture, illustration, etc.
Object name is srep44037-m8.jpg (output-to-reservoir connections) contain real values randomly sampled from a uniform or Gaussian distribution. Weights in An external file that holds a picture, illustration, etc.
Object name is srep44037-m9.jpg and An external file that holds a picture, illustration, etc.
Object name is srep44037-m10.jpg, instead, are trained (typically by means of regularized least-squares algorithm) according to the task at hand.

Once given, the designer can control An external file that holds a picture, illustration, etc.
Object name is srep44037-m11.jpg, An external file that holds a picture, illustration, etc.
Object name is srep44037-m12.jpg, and An external file that holds a picture, illustration, etc.
Object name is srep44037-m13.jpg only through scaling coefficients. For instance An external file that holds a picture, illustration, etc.
Object name is srep44037-m14.jpg is scaled by a multiplicative constant ωo, which has the effect to tune the impact of the ESN output on the next state. In this study, the output feedback is removed by setting ωo = 0. Hence, the resulting state-update difference equation (1) becomes

An external file that holds a picture, illustration, etc.
Object name is srep44037-m15.jpg

while the output equation (2) remains unchanged.

The properties of the ESN described in (3) heavily depend on two hyperparameters. First, input weights in An external file that holds a picture, illustration, etc.
Object name is srep44037-m16.jpg are scaled by a multiplicative constant ωi, which controls the amount of non-linearity introduced by neurons. Large values for ωi tend to saturate the non-linear activation functions. The second hyperparameter is the spectral radius ρ (eigenvalue with largest absolute value) of An external file that holds a picture, illustration, etc.
Object name is srep44037-m17.jpg, which is related to the echo-state property, discussed earlier. For a detailed discussion on the relationships between stability, performance and ρ, we suggest to the interested reader ref. 30 and references therein. Here, it suffices to say that a widely adopted rule-of-thumb49 suggests to set ρ to a value slightly smaller than 1 (e.g., ρ = 0.99). However, to reach higher performance in some practical tasks, it could be necessary to pick a small value for ρ or to breach the aforementioned “safety” bound and push its value beyond unity. Note that in this latter case asymptotic stability of the ESN might still hold, even if some assumptions are locally violated.

In this work, we focus only on the tuning of the hyperparamers ρ and ωi, without affecting generality of the proposed methods that can be adopted to tune also other hyperparamers. The optimal values of ρ and ωi for the task at hand are typically identified with a cross-validation procedure, with the potential associated shortcomings mentioned above. For that reason, different unsupervised approaches have been proposed for their tuning30,31,33. For example, the effect of ρ and ωi on the ESN computational capability can be investigated through the maximal local Lyapunov exponent (MLLE), which measures the divergence rate in state space of trajectories with similar initial conditions. In autonomous (not input-driven) systems, chaos occurs when the maximal Lyapuanov exponent becomes positive, while in input-driven systems, like ESN, one typically relies on local first-order approximations of this quantity (see ref. 30 for details). Accordingly, the onset of criticality in ESNs can be detected by checking when MLLE crosses 0. Another quantity, which was shown to be more accurate in detecting criticality in dynamic systems and well-correlated with ESN performance, is the minimal singular value (in average, over time) of the reservoir Jacobian, denoted as λ33. λ is unimodal and in correspondence of its maximum the dynamical system is far from singularity, has many degrees of freedom, has a good excitability, and it separates well the input signals in state space33. By assuming a null input in Eq. (3), the Jacobian matrix of the reservoir at time t is given by An external file that holds a picture, illustration, etc.
Object name is srep44037-m18.jpg, where the diag(·) operator returns a diagonal matrix. In this paper, λ will serve as a baseline for comparison to our proposed graph-based unsupervised methods for improved ESN hyperparameter tuning, as discussed in the next section.

Horizontal visibility graph and multiplex network

The HVG41 associated with a finite univariate time series An external file that holds a picture, illustration, etc.
Object name is srep44037-m19.jpg, is constructed by assigning a vertex vt to each datum x[t]. The adjacency matrix A characterizes the graph: two vertices vi and vj, i  j are connected by an edge (A[i, j] = 1) iff the corresponding data fulfill the criterion x[ti], x[tj] > x[tp], ∀ ti < tp < tj. In a multivariate scenario, the data stream is composed of Nr different time series An external file that holds a picture, illustration, etc.
Object name is srep44037-m20.jpg, of equal length tmax. A multivariate time series can be mapped into a multiplex visibility graph An external file that holds a picture, illustration, etc.
Object name is srep44037-m21.jpg with Nr layers40. Specifically, the lth multiplex layer is defined by the HVG Gl constructed from xl. In the multiplex, a vertex is replicated on all layers and such replicas are linked by inter-layer connections, while intra-layer connections might change in each layer. From now on, we denote vl[t] to be the vertex of Gl in layer l associated with time interval t.

In this paper we introduce a weighted HVG (wHVG), with edge values defined as An external file that holds a picture, illustration, etc.
Object name is srep44037-m22.jpg An external file that holds a picture, illustration, etc.
Object name is srep44037-m23.jpg. Since self-loops are forbidden in HVGs, edge weights are always well-defined (i.e., finite). The use of weights permits to capture additional information, as it accounts for distance in time (j  i) and amplitude differences (x[i]  x[j]) of two data points connected by the visibility rule. This weighting scheme is motivated by our need to characterize and exploit the instantaneous state by means of a suitable measure of heterogeneity (discussed in the next section). To distinguish the original adjacency matrix from the one of the wHVG, we refer the former as binary adjacency matrix and the latter as weighted adjacency matrix. Algorithm 2 delivers the pseudo-code for constructing a HVG (and a wHVG) from time series x. The worst case complexity of this algorithm is An external file that holds a picture, illustration, etc.
Object name is srep44037-m24.jpg, which occurs when values in x are monotonically decreasing. Instead, the best case complexity is An external file that holds a picture, illustration, etc.
Object name is srep44037-m25.jpg and arises in correspondence of monotonically increasing values in x.

An external file that holds a picture, illustration, etc.
Object name is srep44037-m91.jpg

Vertex properties

Let us consider an ESN with a reservoir of size Nr, driven by an input of length tmax. According to (3), the ESN generates Nr time series An external file that holds a picture, illustration, etc.
Object name is srep44037-m30.jpg, of neuron activations (state). The multivariate time series An external file that holds a picture, illustration, etc.
Object name is srep44037-m31.jpg is represented here by a multiplex An external file that holds a picture, illustration, etc.
Object name is srep44037-m32.jpg. Indexes of the vertices on each layer l of An external file that holds a picture, illustration, etc.
Object name is srep44037-m33.jpg are associated one-to-one with the time index of the original time series. Hence, the ESN state An external file that holds a picture, illustration, etc.
Object name is srep44037-m34.jpg at time t is represented by the vertices An external file that holds a picture, illustration, etc.
Object name is srep44037-m35.jpg and by vertex properties An external file that holds a picture, illustration, etc.
Object name is srep44037-m36.jpg. In Table 1 we introduce four indexes: vertex degree, clustering coefficient50, betweenness and closeness centrality51.

Table 1
Vertex characteristics used to generate different instances of the vertex property vector ϕt.

Heterogeneity of neurons dynamics

The capability of an ESN to reproduce the dynamics of a target system, hence to predict its trajectory in state space, is maximized on the edge of criticality, where the internal dynamical patterns of an ESN become sufficiently rich. In the literature, such a “richness” is usually expressed in terms of diversity of connection weights52, entropy or rank of the matrix of neuron activations15,5,32.

In the same spirit, in order to find hyperparameters yielding maximal prediction accuracy, here we look for those hyperparameter configurations giving rise to neuron activations that are as heterogeneous as possible. Figure 2 provides a visual example of the concept herein discussed. As an illustration, we consider an ESN driven by a sinusoid. We depict the neuron activations and the corresponding HVGs. We select a time step (marked by a black square in the picture) and we show the correspondence between the element in each time series of the neuron activations and the vertices in the associated HVGs. When the ESN operates with a contractive dynamic (2a), neuron activations weakly depend on previous ESN states and they are all very similar to the input signal. This results in a lack of diversity among activations. Accordingly, the corresponding HVGs in the multiplex contain vertices with similar properties, e.g., similar degree and clustering coefficient. When the ESN approaches the edge of criticality, neuron activities highly depend on previous internal states, which encode information of past inputs and internal structure of neuron connections. As we can see from (2b), on the edge of criticality, the activations have different frequencies and their phases are shifted. Such a heterogeneity in the ESN instantaneous state is captured by vertex properties of the multiplex. Finally, if pushed beyond the edge of criticality, the ESN transits into a chaotic regime as shown in (2c). In this case, neuron activations become noise-like and disordered oscillations generate HVGs with vertex properties very different from previous configurations. However, the diversity of the patterns in different time series disappears, hence heterogeneity is again lost. This lack of variety is also highlighted by recurrence of similar motifs in the corresponding HVGs.

Figure 2
Heterogeneity of the graph-based representation of the ESN instantaneous state.

To determine the heterogeneity of neuron activations, we consider the entropy of the related vertex property distribution. In particular, heterogeneity is computed as follows:

  1. At each time t, evaluate the vertex properties An external file that holds a picture, illustration, etc.
Object name is srep44037-m37.jpg using one of the properties listed in Table 1.
  2. Estimate distribution p*[t]) using a histogram with b bins.
  3. Compute instantaneous entropy: Ht = H(p*[t])), where H(·) is the Shannon entropy.
  4. Define heterogeneity as average entropy over time instants: An external file that holds a picture, illustration, etc.
Object name is srep44037-m38.jpg.

The procedure described above, is visually represented in Fig. 3. By referring to the figure, we observe that, at given time t, the ESN state is represented by those vertices in the multiplex tagged with the same label across the layers. For example, the ones in red describe the instantaneous state of the ESN at time step t4 and a vector of vertex properties Φ*[t4] is computed for this set of vertices. Then, we estimate the distribution p*[t4]) and compute the entropy An external file that holds a picture, illustration, etc.
Object name is srep44037-m39.jpg. In order to compute An external file that holds a picture, illustration, etc.
Object name is srep44037-m40.jpg, the procedure is repeated for each time step. Finally, An external file that holds a picture, illustration, etc.
Object name is srep44037-m41.jpg is used to characterize the current hyperparameter configuration. The average entropy depends on the specific vertex property chosen for the analysis (we taken into account four properties, e.g., vertex degree, clustering coefficient, betweenness and closeness centrality – see Table 1). These vertex properties lead to four different entropy values An external file that holds a picture, illustration, etc.
Object name is srep44037-m42.jpg, An external file that holds a picture, illustration, etc.
Object name is srep44037-m43.jpg, An external file that holds a picture, illustration, etc.
Object name is srep44037-m44.jpg, and An external file that holds a picture, illustration, etc.
Object name is srep44037-m45.jpg.

Figure 3
High-level illustration of the proposed methodology.

To summarize, given an input signal, we select hyperparameter configurations that maximize such entropy values. This criterion is inspired by the aforementioned observation linking performance of a computational dynamical system (i.e., prediction accuracy and memory) with heterogeneity of its dynamics (critical dynamics). Accordingly, in this case this information is exploited in order to derive, in an unsupervised way, the configuration yielding highest prediction accuracy.

Other multiplex complexity measures

Recently40, two measures have been proposed in order to characterize the dynamics of a system observed through a multivariate time series and represented as a multiplex composed of HVGs. Here, we consider these measures in order to evaluate whether they are useful for identifying the hyperparameter configurations yielding maximum accuracy performances or not.

The Average Edge Overlap (AEO) computes the expected number of layers of the multiplex on which an edge is present. For binary HVGs, it is defined as

An external file that holds a picture, illustration, etc.
Object name is srep44037-m46.jpg

where Al[i, j] equals 1 if vertices i and j are connected in layer l. The second measure is the Interlayer Mutual Information (IMI), which quantifies the correlations between the degree distributions of two different layers li and lj. It is defined, for binary HVGs, as

An external file that holds a picture, illustration, etc.
Object name is srep44037-m47.jpg

where An external file that holds a picture, illustration, etc.
Object name is srep44037-m48.jpg is the joint probability to find a vertex with degree An external file that holds a picture, illustration, etc.
Object name is srep44037-m49.jpg in layer li and degree An external file that holds a picture, illustration, etc.
Object name is srep44037-m50.jpg in layer lj, respectively. In the experiments, we use the average IMI between all pairs of multiplex layers.

Memory measures

The memory of an ESN is quantified by the capability to retain information about past inputs in its transient dynamics53. A supervised measure called memory capacity (MC) is usually adopted to quantify ESN memory27. When computing MC, the reservoir topology and weights are kept fixed and several readout layers are trained in order to reproduce delayed versions of the input at different time lags An external file that holds a picture, illustration, etc.
Object name is srep44037-m51.jpg. The MC is computed as the squared correlation coefficient between desired (delayed input) and computed outputs,

An external file that holds a picture, illustration, etc.
Object name is srep44037-m52.jpg

In order to properly evaluate the capability of ESNs to introduce memory through recurrent connections, x[t] is defined as a stationary uncorrelated noisy signal.

In the following, we propose an unsupervised graph-based method to identify hyperparameter configurations for which an ESN achieves large memory capacity. Given the input time series An external file that holds a picture, illustration, etc.
Object name is srep44037-m53.jpg, we determine if there exists a subset of neuron activations, which is correlated with a past input sequence xa,b = {x[t  τa], …, x[t  τb]}, with τa > τb and τa, τb [set membership] [1, tmax  1]. Being Gl the HVG representing the lth layer of the multiplex, An external file that holds a picture, illustration, etc.
Object name is srep44037-m54.jpg is the sequence of its vertex degrees ordered according to the time index (not to be confused with An external file that holds a picture, illustration, etc.
Object name is srep44037-m55.jpg, the degrees of vertices relative to time t across the different layers). With Gx, instead, we refer to the HVG constructed over the input x[t], while An external file that holds a picture, illustration, etc.
Object name is srep44037-m56.jpg is the vector of its vertex degrees.

First, we define a measure of maximum agreement between An external file that holds a picture, illustration, etc.
Object name is srep44037-m57.jpg and each sequence An external file that holds a picture, illustration, etc.
Object name is srep44037-m58.jpg as An external file that holds a picture, illustration, etc.
Object name is srep44037-m59.jpg. κ*(·,·) is a similarity measure between sequences: in this paper we consider the Pearson correlation κPC(·,·), the Spearman correlation κSC(·,·), and the mutual information κMI(·,·).

A second measure of agreement is defined on the adjacency matrix Ql = Ax  Al, where Ql[i, j] = 1 if Ax[i, j] = 1 and Al[i, j] = 1; Ql[i, j] = 0 otherwise. The agreement is then computed as the largest number of edges among all possible intersection graphs, An external file that holds a picture, illustration, etc.
Object name is srep44037-m60.jpg, where |·| counts the number of ones in the adjacency matrix.

Finally, one might consider the similarity between original input x and neuron activations hl, An external file that holds a picture, illustration, etc.
Object name is srep44037-m61.jpg. In this case, κ*(·,·) is directly evaluated on the time series rather than on the sequences of vertex degrees. This last measure of agreement is taken into account in order to quantitatively show (in the experiments) the benefits of using HVGs for representing neuron activations.

By referring to the illustrative example in Fig. 4, we generate the HVGs Gx, An external file that holds a picture, illustration, etc.
Object name is srep44037-m62.jpg, relative to delayed input x15,10 and activations of the Nr reservoir neurons, respectively. Successively, we evaluate their similarities by means of an agreement measure (δDG in this example). The similarity measure κ* is chosen among the three previously proposed measures: κPC, κSC, or κMI.

Figure 4
We extract a subsequence from the history of the input signal, for example x15,10, and compute the related HVG Gx.


In the following, we perform two experiments in order to evaluate the two proposed unsupervised methods for, respectively, finding hyperparameter configurations giving rise to ESNs with high prediction accuracy and large memory capacity. In the first experiment, we show that, on different real and synthetic tasks, (supervised) prediction accuracy is maximized for the same hyperparameter configurations that yields the largest heterogeneity for the vertex properties of the multiplex. In the second experiment, we show the reliability of the graph-based memory measures in identifying hyperparameters where the (supervised) MC is maximized.

Test for prediction accuracy

In this experiment, we consider several prediction tasks and, for each of them, we set the forecast step τf > 0 to be the smallest time-lag that guarantees the measurements in a time window of size τf to be uncorrelated (e.g., the first zero in the autocorrelation function of the input signal). Prediction error is evaluated by Normalized Root Mean Squared Error, An external file that holds a picture, illustration, etc.
Object name is srep44037-m63.jpg where An external file that holds a picture, illustration, etc.
Object name is srep44037-m64.jpg is the prediction provided by the ESN and y[t] is the desired/teacher output. The Prediction accuracy is defined as γ = max{1  NRMSE, 0}.

In the following, we describe the datasets used in this experimental campaign.

Sinusoidal input

We feed an ESN with a sinusoid y(t) = sin(ψt) and we predict future input values with a forecast step τf = 2π/ψ.

Mackey-Glass time series

The Mackey-Glass (MG) system is commonly used as a benchmark in chaotic time series prediction. The input signal is generated from the MG time-delay differential equation

An external file that holds a picture, illustration, etc.
Object name is srep44037-m65.jpg

We adopt the standard parameters τMG = 17, α = 0.2, β = 0.1, initial condition x(0) = 1.2, and integration step equal to 0.1. The forecast step here is τf = 6.

Multiple superimposed oscillator

Prediction of superimposed sinusoidal waves with incommensurable frequencies is a hard forecasting exercise, due to the extension of the wavelength54. The ESN is fed with the multiple superimposed oscillator (MSO)

An external file that holds a picture, illustration, etc.
Object name is srep44037-m66.jpg

For this task, the ESN is trained to predict future input values, with forecast horizon τf = 16.


The chosen Non-Linear Auto-Regressive Moving Average (NARMA) tas55 consists in modeling the output of the r-order system:

An external file that holds a picture, illustration, etc.
Object name is srep44037-m67.jpg

x[t] is a uniform random noise signal in [0, 1] and is the input of the ESN, which is trained to reproduce y[t + 1]. The NARMA task is known to require a memory of at least r past time-steps, since the output is determined by inputs and outputs from the last r time-steps. For this task we set r = 20 and τf = 15.

Polynomial task

The ESN is fed with uniform noise in [−1, 1] and is trained to reproduce the following output

An external file that holds a picture, illustration, etc.
Object name is srep44037-m68.jpg

where ci,j is uniformly distributed in [0, 1]56. The difficulty of prediction can be controlled by varying the polynomial degree p and the time delay d. For this task we set p = 7, d = 10, and τf = d.

Telephone call load time series

As a last test, we consider a real-world dataset relative to the load of phone calls registered over a mobile network. The data comes from the Orange telephone dataset, published in the Data for Development (D4D) challenge57. D4D is a collection of call data records, containing anonymized events of Orange’s mobile phone users in the Ivory Coast, in a period spanning from December 1, 2011 to April 28, 2012. The dataset consists of 6 time series consisting of: number and volume of incoming calls, number and volume of outgoing calls, day and time (1 hour resolution) when the telephone activity was registered. More detailed information is available at the related website58. All 6 time series are fed into the ESN as inputs; the goal is to predict 6 hours ahead the volume of incoming calls – the profile of this latter time series is depicted in Fig. 5.

Figure 5
D4D dataset – load profile of incoming calls volume for the first 300 time intervals.

In each test, we evaluate the correlation between the average entropy of vertex properties in the multiplex and the prediction accuracy γ, as we vary the hyperparameters ρ and ωi. Multiplexes are generated using both the binary and weighted version of the HVG adjacency matrix. To appreciate the effectiveness of our methodology, we also estimate the correlations of γ with λ, the minimum singular value of the Jacobian of the reservoir (see Methods). Additionally, we consider the correlation of γ with the two layer-based measures IMI and AEO (see Methods). Correlations are evaluated as follows. For each configuration (ρ = k, ωi = j), we have the prediction accuracy γk,j, the entropy An external file that holds a picture, illustration, etc.
Object name is srep44037-m69.jpg, and so on. The values assumed by these quantities by varying ρ and ωi generate a two-dimensional manifold. The point-wise linear (Pearson) correlation between the manifolds relative to γ and the other considered measures is the result we are interested in. The configurations of the hyperparameters that are examined are generated by varying ρ in [0.5, 1.3] (20 different values) and ωi in [0.2, 0.9] (10 different values). A total of 200 configurations are evaluated. Due to the stochastic nature of the ESN initialization, for each configuration (ρ = k, ωi = j) we compute γk,j and all the other measures 15 different times. Successively, we compute the correlations among their average values. We use a reservoir with Nr = 100 neurons and sparsity of the internal connectivity equal to 25%. The readout is trained by standard ridge regression with regularization parameter set to 0.05. The distributions of the vertex properties are estimated by histograms with b = 50 bins. We comment that the reservoir size Nr can be increased/decreased without affecting the applicability of the proposed methodology (only the number of bins used to estimate the vertex properties distribution might be modified – see Methods). However, the number of neurons has an impact on the overall performance of the network and on the time and space complexities of the proposed method.

Figure 6 depicts the values assumed by γ and four graph-based measures in the case of the MSO prediction task. As can be seen, high correlation emerges between the (average) entropy of the vertex properties and the prediction accuracy. Since our approach is fully unsupervised, the proposed graph-based measures can approximate well the accuracy γ, regardless of the task learned by the readout (prediction, function approximation, reconstruction of past inputs, etc). In Table 2, we report the average correlation values and their statistical significance (expressed by p-values) on all tasks. As we can see, the highest (and statistically significant) correlation is achieved by using one of the four average entropy measures of vertex properties. In particular, the measure based on vertex cluster coefficient distribution, An external file that holds a picture, illustration, etc.
Object name is srep44037-m70.jpg, achieves the best results in 5 of the 6 tasks. For what concerns the D4D time series, we observe that λ achieves high correlation with γ, but still lower than the one achieved by An external file that holds a picture, illustration, etc.
Object name is srep44037-m71.jpg. This demonstrates the effectiveness of the proposed methodology, also in the case of a real-world application. In SIN and MSO tasks, the graph-based quantifiers estimated on the weighted HVG achieve a higher degree of accuracy, with respect to the binary counterpart. In these cases, additional qualitative information relative to temporal and amplitude differences in the connected data allows to better represent the dynamics of the system. Finally, it is worth noting that IMI takes high, yet negative correlation values on both MG and POLY tasks. In such cases, results are close to the ones achieved with our approach.

Figure 6
In each panel, the two-dimensional gray manifold represents the values of the prediction accuracy γ in the MSO task, for different configurations of ρ and ωi. The colored manifolds are the average entropy of the distributions of ...
Table 2
Correlations and related p-values (in brackets) of the proposed graph-based measures with accuracy γ in different prediction tasks, as ρ and ωi change.

Test for memory capacity

The performed experiment consists in generating 100 different random reservoirs, each one characterized by an increasing value of spectral radius ρ in the [0.1, 2] interval. As ρ varies, we evaluate the MC by training four readouts in order to reproduce different time-lagged versions of input signal x10,5, …, x25,20. Then, on the output of each reservoir, we evaluate the similarities δTS, δDG, and δAND, which are high if there exists at least one series of activations that is similar to the considered past input sequence. This is evaluated in such as way that the measure κ* taken into account. Even if some neurons retain dynamics of previous input sequences, the reservoir introduces shifts in the phase and the amplitude of the input signal. To filter out these differences, in this test we consider only HVGs defined by binary adjacency matrices, which do not account for differences in the amplitude of the connected values. To evaluate the effectiveness of the proposed unsupervised memory measures, we compute the correlation between the supervised MC and δTS, δDG, δAND, as ρ varies within the chosen interval. Note that we only monitor the effect of ρ on the dynamics, since it is the hyperparameter that mostly affects the memory capacity32. We kept the input scaling fixed, ωi = 0.7, while the remaining hyperparamers are configured as in the previous experiment. As before, we repeated each experiment 15 times with different and independent random initializations. In Table 3, we show the mean correlation values, along with the standard deviations, between the MC and the proposed unsupervised measures of memory capacity.

Table 3
Mean correlations and standard deviations of MC with the unsupervised memory quantifiers δ TS, δ DG, and δ AND.

From the table, we observe that the best agreement with the MC is achieved by the measures derived from the HVGs. In particular, δDG configured with the Spearman rank κSC is always highly correlated with MC and, in three of the four delayed input sequences taken into account, is the best performing one. In each setup, δDG works better if configured with κSC rather than κPC. Instead, results obtained with κMI are significantly worse in all cases. δAND achieves the best results only for the first time lag taken into account, while the agreement with MC is lower in the remaining cases. Interestingly, several measures show a high degree of correlation with the MC as the size of the delay increases. δTS, the unsupervised measure computed directly on the input time series and neuron activations, shows positive correlations with the MC, but the agreement is always lower with respect to the graph-based measures. For δTS, the setting with κPC works better than κSC. Finally, also in this case by using κMI we obtain the worst performance. In Fig. 7, we show an example of the values of MC, δTS (configured with κPC), δDG (configured with κSC), and δAND, as ρ is varied within the [0.2, 2] interval.

Figure 7
Values of supervised MC and of two selected unsupervised memory measures, when the input sequence x20,15 is taken into account.


Experimental results show satisfactory correlations for average entropy of vertex properties with respect to prediction accuracy. Moreover, the two unsupervised graph-based memory measures that we proposed (δDG and δAND) correlate well with the supervised measure of memory capacity.

We first discuss the results of the prediction accuracy test, where we analyzed topological properties of vertices in the multiplex, representing the ESN instantaneous state. On all tests taken into account, we observed a remarkable correlation between γ and the average entropy of the clustering coefficient distribution An external file that holds a picture, illustration, etc.
Object name is srep44037-m72.jpg, hence suggesting that the clustering coefficient is able to describe well the heterogeneity of the activations. To explain this result, it is necessary to elaborate on the properties of the clustering coefficient CL(·). In the HVG literature, CL(·) its behavior has been analyzed for time series characterized by different Hurst exponent59. Additionally, an upper bound (CL(v) [set membership] [0, 2/DG(v)]) is provided for HVGs derived from random time series41. In the following, we present an in-depth interpretation of the results by accounting for geometrical properties of the clustering coefficient.

In a HVG, CL(v) measures the inter-visibility among neighbors of v. For convex functions, it is possible to connect any two points with a straight line. This feature is also (partially) captured by the HVG. If v is contained in a convex part of the related time series, there is a high degree of intervisibility among the neighbor vertices to which v is connected, hence CL(v) is high. Additionally, moving along the same convex part of the time series results only in minor changes of the clustering coefficient in the associated HVG vertices. Instead, if v is a local maximum of a concave part, then it is connected to points belonging to two different basins, which do not have reciprocal visibility. In this case, CL(v) is low and its value rapidly changes as one moves away from the maximum. This results in great losses of visual information. Therefore, large values of CL(·) indicate the presence of dominating convexities, while low values characterize concavities60. Accordingly, CL(·) can be used to measure the length of a convex (concave) part of the time series and how fast the convexity is changing, which is a measure of the fluctuations in the time series61. In a regime characterized by contractive dynamics, convexity changes at the same (slow) rate in different neuron activations and this results in a low entropy value of the clustering coefficient distribution among vertices in different layers. On the edge of criticality, instead, convex and concave parts in the time series of activations are characterized by heterogeneous lengths and they change at different rates. This corresponds to a high degree of clustering coefficient diversity of the same HVG vertex, replicated at different layers in the multiplex. Finally, in the chaotic regime, all time series fluctuate very rapidly and their convexity changes every few time steps. In this case, in each time series of activations the lengths of convex and concave parts are always very short and hence the desired heterogeneity is again lost.

For what concerns experiments on memory, the best overall results in terms of agreement with the supervised MC are achieved by the graph-based measure δDG. As previously discussed, such a measure evaluates the maximum similarity between the sequence of vertex degrees on the input HVG Gx and the HVG Gl of neuron activations. This measure is closely related with the degree distribution P(k), whose importance is known in the HVG literature41. For example, it has been shown that for time series generated from an i.i.d. process, P(k) follows P(k) = (1/3)(2/3)k−2 and the mean degree is left angle bracketkright angle bracket = 4. As the correlations in the time series increase, the i.i.d. assumption is lost and P(k) decays faster. Furthermore, vertex degrees are key parameters to describe dynamic processes on the graph, such as synchronization of coupled oscillators, percolation, epidemic spreading, and linear stability of equilibrium in networked coupled systems62. Their role has been studied also in the HVG framework63. HVGs have been studied in the context of time series related to processes with power-law correlations59. In our case, the time series of neuron activations have short-term correlations and increments in the correlation coefficients can have opposite signs at consecutive time lags. For these cases, we are not aware of any previous study in terms of HVGs.

In networks which are inherently degree disassortative, the range of degree values increases with network size, with a consequent decrease of the assortativity value64. In such networks, the Spearman rank correlation provides a more suitable choice with respect to calculating degree-degree Pearson correlations. It is important to notice that the rank is computed through a non-linear rescaling, which is data dependent. The information on the actual values of the data is discarded as only its inherent ordering (rank) is preserved. We argue that HVGs convey the same type of information captured by the Spearman correlation. Hence, the latter should be preferred to Pearson correlation to characterize the characteristics of the vertices and related topological properties in HVGs. This fact justifies the higher agreement with memory capacity achieved by means of δDG when configured with κSC, which accounts for Spearman correlations between sequences of vertex degrees in the HVGs related to the input signal and the neuron activations.

Modeling ESN dynamics through a multiplex network allowed us to connect two seemingly different research fields, thus fostering multidisciplinary research in the context of recurrent neural networks. By converting a temporal problem into a topological one, we handled temporal dependencies introduced by ESNs (as well as by other types of RNNs), hence overcoming technical limitations of statistical approaches that require independence of samples. We performed and discussed several experiments that provided empirical evidence that our methodology achieves performance higher than other unsupervised methods and comparable to cross-validation techniques. These results suggest to allocate efforts to further improve the effectiveness of unsupervised learning methods in the context of ESNs and RNNs. Finally, we would like to stress that, while this paper is primarily focused on network structures in machine learning, our results might suggest new ideas for theoretical understanding of recurrent structures in biological models of neuronal networks7,65,66. In particular, we believe that it is possible to identify emergent structural patterns in the developed graph-based representations of network dynamics. This would allow to further explore and analyze the route to chaos in input-driven neural models by exploiting the language of graph theory, which is an already established framework within the neuroscience field67.

Additional Information

How to cite this article: Bianchi, F. M. et al. Multiplex visibility graphs to investigate recurrent neural networks dynamics. Sci. Rep. 7, 44037; doi: 10.1038/srep44037 (2017).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors declare no competing financial interests.

Author Contributions L.L. outlined the research ideas. F.M.B. and L.L. conceived methods. F.M.B. performed experiments and generated figures. C.A. and R.J. contributed to the technical discussion. All authors took part to the paper writing and approved the final manuscript.


  • Barzel B. & Barabási A.-L. Universality in network dynamics. Nature Physics 9, 673–681 (2013).
  • Hammer B., Micheli A., Sperduti A. & Strickert M. Recursive self-organizing network models. Neural Networks 17, 1061–1085 (2004). [PubMed]
  • Maass W., Joshi P. & Sontag E. D. Computational aspects of feedback in neural circuits. PLoS Computational Biology 3, e165 (2007). [PMC free article] [PubMed]
  • Reinhart R. F. & Steil J. J. Regularization and stability in reservoir networks with output feedback. Neurocomputing 90, 96–105 (2012).
  • Jordan M. I. Serial order: A parallel distributed processing approach. In Donahoe J. W. & Dorsel V. P. (eds) Neural-Network Models of Cognition: Biobehavioral Foundations, vol. 121 of Advances in Psychology, 471–495 (North-Holland, 1997).
  • Amit D. J., Gutfreund H. & Sompolinsky H. Spin-glass models of neural networks. Physical Review A 32, 1007–1018 (1985). [PubMed]
  • Enel P., Procyk E., Quilodran R. & Dominey P. F. Reservoir computing properties of neural dynamics in prefrontal cortex. PLoS Computational Biology 12, e1004967 (2016). [PMC free article] [PubMed]
  • Barak O., Sussillo D., Romo R., Tsodyks M. & Abbott L. F. From fixed points to chaos: three models of delayed discrimination. Progress in Neurobiology 103, 214–222 (2013). [PMC free article] [PubMed]
  • Rajan K., Abbott L. F. & Sompolinsky H. Stimulus-dependent suppression of chaos in recurrent neural networks. Physical Review E 82, 011903 (2010). [PubMed]
  • Fusi S., Miller E. K. & Rigotti M. Why neurons mix: high dimensionality for higher cognition. Current Opinion in Neurobiology 37, 66–74 (2016). [PubMed]
  • Schmidhuber J. Deep learning in neural networks: An overview. Neural Networks 61, 85–117 (2015). [PubMed]
  • Barra A., Bernacchia A., Santucci E. & Contucci P. On the equivalence of Hopfield networks and Boltzmann machines. Neural Networks 34, 1–9 (2012). [PubMed]
  • Charles A., Yin D. & Rozell C. Distributed sequence memory of multidimensional inputs in recurrent networks. arXiv preprint arXiv:1605.08346 (2016).
  • Tiňo P. & Rodan A. Short term memory in input-driven linear dynamical systems. Neurocomputing 112, 58–63 (2013).
  • Legenstein R. & Maass W. What makes a dynamical system computationally powerful? In Haykin S., Principe J., Sejnowski & McWhirter (eds) New Directions in Statistical Signal Processing: From Systems to Brain, 127–154 (MIT Press, Cambridge, 2007).
  • Yildiz I. B., Jaeger H. & Kiebel S. J. Re-visiting the echo state property. Neural Networks 35, 1–9 (2012). [PubMed]
  • Dambre J., Verstraeten D., Schrauwen B. & Massar S. Information processing capacity of dynamical systems. Scientific Reports 2 (2012). [PMC free article] [PubMed]
  • Lajoie G., Lin K. K., Thivierge J.-P. & Shea-Brown E. Encoding in balanced networks: Revisiting spike patterns and chaos in stimulus-driven systems. PLoS Computational Biology 12, 1–30 (2016). [PMC free article] [PubMed]
  • Livi L., Bianchi F. M. & Alippi C. Determination of the edge of criticality in echo state networks through Fisher information maximization. IEEE Transactions on Neural Networks and Learning Systems 1–12, doi: (2017).10.1109/TNNLS.2016.2644268 [PubMed] [Cross Ref]
  • Mayer N. M. Input-anticipating critical reservoirs show power law forgetting of unexpected input events. Neural Computation 27, 1102–1119 (2015). [PubMed]
  • Moretti P. & Muñoz M. A. Griffiths phases and the stretching of criticality in brain networks. Nature Communications 4 (2013). [PubMed]
  • Mora T. & Bialek W. Are biological systems poised at criticality? Journal of Statistical Physics 144, 268–302 (2011).
  • Scheffer M. et al. . Anticipating critical transitions. Science 338, 344–348 (2012). [PubMed]
  • Wang X., Lizier J. & Prokopenko M. Fisher information at the edge of chaos in random Boolean networks. Artificial Life 17, 315–329 (2011). [PubMed]
  • Tkačik G. et al. . Thermodynamics and signatures of criticality in a network of neurons. Proceedings of the National Academy of Sciences 112, 11508–11513 (2015). [PubMed]
  • Marković D. & Gros C. Power laws and self-organized criticality in theory and nature. Physics Reports 536, 41–74 (2014).
  • Lukoševičius M. & Jaeger H. Reservoir computing approaches to recurrent neural network training. Computer Science Review 3, 127–149 (2009).
  • Grigoryeva L., Henriques J., Larger L. & Ortega J.-P. Optimal nonlinear information processing capacity in delay-based reservoir computers. Scientific Reports 5 (2015). [PMC free article] [PubMed]
  • Bianchi F. M., Scardapane S., Uncini A., Rizzi A. & Sadeghian A. Prediction of telephone calls load using echo state network with exogenous variables. Neural Networks 71, 204–213 (2015). [PubMed]
  • Bianchi F. M., Livi L. & Alippi C. Investigating echo state networks dynamics by means of recurrence analysis. IEEE Transactions on Neural Networks and Learning Systems 1–13, doi: (2016).10.1109/TNNLS.2016.2630802 [PubMed] [Cross Ref]
  • Boedecker J., Obst O., Lizier J. T., Mayer N. M. & Asada M. Information processing in echo state networks at the edge of chaos. Theory in Biosciences 131, 205–213 (2012). [PubMed]
  • Ozturk M. C., Xu D. & Prncipe J. C. Analysis and design of echo state networks. Neural Computation 19, 111–138 (2007). [PubMed]
  • Verstraeten D. & Schrauwen B. On the quantification of dynamics in reservoir computing. In Artificial Neural Networks-ICANN 2009 985–994 (Springer: Berlin Heidelberg, 2009).
  • Leisch F., Trapletti A. & Hornik K. Stationarity and stability of autoregressive neural network processes. In Kearns M. J., Solla S. A. & Cohn D. A. (eds) Advances in Neural Information Processing Systems 11, 267–273 (MIT Press, 1999).
  • Lee K.-M., Min B. & Goh K.-I. Towards real-world complexity: an introduction to multiplex networks. The European Physical Journal B 88, 1–20 (2015).
  • Boccaletti S. et al. . The structure and dynamics of multilayer networks. Physics Reports 544, 1–122 (2014).
  • Menichetti G., Remondini D., Panzarasa P., Mondragón R. J. & Bianconi G. Weighted multiplex networks. PloS ONE 9, e97857 (2014). [PMC free article] [PubMed]
  • Kivelä M. et al. . Multilayer networks. Journal of Complex Networks 2, 203–271 (2014).
  • De Domenico M. et al. . Mathematical formulation of multilayer networks. Physical Review X 3, 041022 (2013).
  • Lacasa L., Nicosia V. & Latora V. Network structure of multivariate time series. Scientific Reports 5 (2015). [PMC free article] [PubMed]
  • Luque B., Lacasa L., Ballesteros F. & Luque J. Horizontal visibility graphs: Exact results for random time series. Physical Review E 80, 046103 (2009). [PubMed]
  • Luque B., Lacasa L., Ballesteros F. J. & Robledo A. Feigenbaum graphs: A complex network perspective of chaos. PLoS One 6, e22411 (2011). [PMC free article] [PubMed]
  • Luque B., Cordero-Gracia M., Gómez M. & Robledo A. Quasiperiodic graphs at the onset of chaos. Physical Review E 88, 062918 (2013). [PubMed]
  • Lacasa L. & Toral R. Description of stochastic and chaotic series using visibility graphs. Physical Review E 82, 036120 (2010). [PubMed]
  • Ravetti M. G., Carpi L. C., Gonçalves B. A., Frery A. C. & Rosso O. A. Distinguishing noise from chaos: Objective versus subjective criteria using horizontal visibility graph. PloS ONE 9, e108004 (2014). [PMC free article] [PubMed]
  • Luque B., Lacasa L., Ballesteros F. J. & Robledo A. Analytical properties of horizontal visibility graphs in the Feigenbaum scenario. Chaos: An Interdisciplinary Journal of Nonlinear Science 22, 013109 (2012). [PubMed]
  • Zhang S. et al. . Architectural complexity measures of recurrent neural networks. arXiv preprint arXiv:1602.08210 (2016).
  • Hermans M. & Schrauwen B. Recurrent kernel machines: Computing with infinite echo state networks. Neural Computation 24, 104–133 (2012). [PubMed]
  • Lukoševičius M. A Practical Guide to Applying Echo State Networks, 659–686 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012).
  • Saramäki J., Kivelä M., Onnela J.-P., Kaski K. & Kertesz J. Generalizations of the clustering coefficient to weighted complex networks. Physical Review E 75, 027105 (2007). [PubMed]
  • Boccaletti S., Latora V., Moreno Y., Chavez M. & Hwang D. Complex networks: Structure and dynamics. Physics Reports 424, 175–308 (2006).
  • Bertschinger N. & Natschläger T. Real-time computation at the edge of chaos in recurrent neural networks. Neural Computation 16, 1413–1436 (2004). [PubMed]
  • Pascanu R. & Jaeger H. A neurodynamical model for working memory. Neural Networks 24, 199–207 (2011). [PubMed]
  • Jaeger H. & Haas H. Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication. Science 304, 78–80 (2004). [PubMed]
  • Jaeger H. Adaptive nonlinear system identification with echo state networks. In Advances in Neural Information Processing Systems, 593–600 (MIT Press, 2002).
  • Butcher J., Verstraeten D., Schrauwen B., Day C. & Haycock P. Reservoir computing and extreme learning machines for non-linear time-series data analysis. Neural Networks 38, 76–89 (2013). [PubMed]
  • Blondel V. D. et al. . Data for Development: the D4D Challenge on Mobile Phone Data. ArXiv preprint arXiv:1210.0137 (2012).
  • Orange d4d challenge. Accessed: 2016-09-22.
  • Xie W.-J. & Zhou W.-X. Horizontal visibility graphs transformed from fractional Brownian motions: Topological properties versus the Hurst index. Physica A: Statistical Mechanics and its Applications 390, 3592–3601 (2011).
  • Costa L. d. F., Rodrigues F. A., Travieso G. & Villas Boas P. R. Characterization of complex networks: A survey of measurements. Advances in Physics 56, 167–242 (2007).
  • Turner A., Doxa M., O’Sullivan D. & Penn A. From isovists to visibility graphs: A methodology for the analysis of architectural space. Environment and Planning B: Planning and Design 28, 103–121 (2001).
  • Restrepo J. G., Ott E. & Hunt B. R. Approximating the largest eigenvalue of network adjacency matrices. Physical Review E 76, 056119 (2007). [PubMed]
  • Fioriti V., Tofani A. & Di Pietro A. Discriminating chaotic time series with visibility graph eigenvalues. Complex Systems 21 (2012).
  • Newman M. E. J. Assortative mixing in networks. Physical Review Letters 89, 208701 (2002). [PubMed]
  • Marblestone A. H., Wayne G. & Kording K. P. Toward an integration of deep learning and neuroscience. Frontiers in Computational Neuroscience 10 (2016). [PMC free article] [PubMed]
  • Sussillo D., Churchland M. M., Kaufman M. T. & Shenoy K. V. A neural network that finds a naturalistic solution for the production of muscle activity. Nature Neuroscience 18, 1025–1033 (2015). [PMC free article] [PubMed]
  • Sporns O. Networks of the Brain (MIT press, Cambridge, MA, USA, 2011).
  • Lopez-Fernandez L., Robles G. & Gonzalez-Barahona J. M. Applying social network analysis to the information in CVS repositories. In Proceedings of the International Workshop on Mining Software Repositories 101–105 (Edinburgh, UK, 2004).

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group