# On Model Stability as a Function of Random Seed

Pranava Madhyastha
Department of Computing
Imperial College London
<EMAIL_ADDRESS>
Rishabh Jain
Bloomberg
London
<EMAIL_ADDRESS>
This work was conducted when the author was a student at Imperial College
London.

###### Abstract

In this paper, we focus on quantifying model stability as a function of random
seed by investigating the effects of the induced randomness on model
performance and the robustness of the model in general. We specifically
perform a controlled study on the effect of random seeds on the behaviour of
attention, gradient-based and surrogate model based (LIME) interpretations.
Our analysis suggests that random seeds can adversely affect the consistency
of models resulting in counterfactual interpretations. We propose a technique
called  _Aggressive Stochastic Weight Averaging (ASWA)_ and an extension
called  _Norm-filtered Aggressive Stochastic Weight Averaging (NASWA)_ which
improves the stability of models over random seeds. With our ASWA and NASWA
based optimization, we are able to improve the robustness of the original
model, on average reducing the standard deviation of the model’s performance
by $72\%$.

##  1 Introduction

There has been a tremendous growth in deep neural network based models that
achieve state-of-the-art performance. In fact, most recent end-to-end deep
learning models have surpassed the performance of careful human feature-
engineering based models in a variety of NLP tasks. However, deep neural
network based models are often brittle to various sources of randomness in the
training of the models. This could be attributed to several sources including,
but not limited to, random parameter initialization, random sampling of
examples during training and random dropping of neurons. It has been observed
that these models have, more often, a set of _random seeds_ that yield better
results than others. This has also lead to research suggesting random seeds as
an additional hyperparameter for tuning Bengio
(2012)111http://www.argmin.net/2018/02/26/nominal/. One possible explanation
for this behavior could be the existence of multiple local minima in the loss
surface. This is especially problematic as the loss surfaces are generally
non-convex and may have multiple saddle points making it difficult to achieve
model stability.

if high crimes were any more generic it would have a universal product code
instead of a title ($\Pr{(Y_{{negative}})}=0.99$)
---
if high crimes were any more generic it would have a universal product code
instead of a title ($\Pr{(Y_{{negative}})=0.98}$)

Figure 1: Importance based on attention probabilities for two runs of the same
model with same parameters and same hyperparameters, but with two different
random seeds (color magnitudes: pink$<$magenta$<$red)

Recently the NLP community has witnessed a resurgence in interpreting and
explaining deep neural network based models Jain et al. (2019); Jain and
Wallace (2019); Alvarez-Melis and Jaakkola (2017). Most of the interpretation
based methods involve one of the following ways of interpreting models: a)
sample oriented interpretations: where the interpretation is based on changes
in the prediction score with either upweighting or perturbing samples Jain et
al. (2019); Jain and Wallace (2019); Koh and Liang (2017); b) interpretations
based on feature attributions using attention or input perturbation or
gradient-based measures; Ghaeini et al. (2018); Feng et al. (2018); Bach et
al. (2015); c) interpretations using surrogate linear models Ribeiro et al.
(2016) – these methods can provide local interpretations based on input
samples or features. However, the presence of inherent randomness makes it
difficult to accurately interpret deep neural models among other forms of
pathologies Feng et al. (2018).

In this paper, we focus on the stability of deep neural models as a function
of random-seed based effects. We are especially interested in investigating
the hypothesis focusing on model stability: do neural network based models
under different random seeds allow for similar interpretations of their
decisions? We claim that for a given model which achieves a substantial
performance for a task, the factors responsible for any decisions over a
sample should be approximately consistent irrespective of the random seed. In
Figure 1, we show an illustration of this question where we visualize the
attention distributions of two CNN based binary classification models for
sentiment analysis, trained with the same settings and hyper-parameters, but
with different seeds. We observe that both models obtain the correct
prediction with significantly high confidence. However, we note that both the
models attend to completely different sets of words. This is problematic,
especially when interpreting these models under the influence of such
randomness. We observe that on average $40-60\%$ of the most important
interpretable units are different across different random seeds for the same
model. This phenomenon also leads us to the question on the exact nature of
interpretability – are the interpretations specific to an instantiation of the
model or are they general to a class of models?

We also provide a simple method that can, to a large extent, ameliorate this
inherent random behaviour. In Section 3.1, we propose an aggressive stochastic
weight averaging approach that helps in improving the stability of the models
at almost zero performance loss while still making the model robust to random-
seed based instability. We also propose an improvement to this model in
Section 3.2 which further improves the stability of the neural models. Our
proposals significantly improve the robustness of the model, on average by
$72\%$ relative to the original model and on Diabetes (MIMIC), a binary
classification dataset, by $89\%$ (relative improvement). All code for
reproducing and replicating our experiments is released in our
repository222https://github.com/rishj97/ModelStability.

##  2 Measuring Model Stability

In this section, we describe methods that we use to measure model stability,
specifically — prediction and interpretation stability.

###  2.1 Prediction Stability

We measure prediction stability using standard measures of the mean and the
standard deviations corresponding to the accuracy of the classification based
models on different datasets. We ensure that the models are run with exactly
the same configurations and hyper-parameters but with different random seeds.
This is a standard procedure that is used in the community to report the
performance of the model.

###  2.2 Interpretation Stability

For a given task, we train a set of models only differing with random-seeds.
For every given test sample, we obtain interpretations using different
instantiations of the models. We define a model to be stable if we obtain
similar interpretations regardless of different random-seed based
instantiations. We use the following metrics to quantify stability:

a) Relative Entropy quantification ($\mathcal{H}$): Given two distributions
over interpretations, for the same test case, from two different models, it
measures the relative entropy between the two probability distributions. Note
that, the higher the relative entropy the greater the dissimilarity between
the two distributions.

$\mathcal{H}=\sum_{i{\in}d}{{{\Pr}_{1}}\cdot\log{\frac{\Pr_{1}}{\Pr_{2}}}}$

where, $\Pr_{1}$ and $\Pr_{2}$ are two attention distributions of the same
sample from two different runs of the model and $d$ is the number of tokens in
the sample. Given $n$ differently seeded models, for each test instance, we
calculate the relative entropy obtained from the corresponding averaged
pairwise interpretation distributions.

b) Jaccard Distance ($\mathcal{J}$): It measures the dissimilarity between two
sets. Here higher values of $\mathcal{J}$ indicate larger variances. We
consider top-$n$ tokens which have the highest attention for comparison. Note
that, Jaccard distance is over sets of word indices and do not take into
account the attention probabilities explicitly. Jaccard distance is defined
as:

$\mathcal{J}=(1-\frac{A{\cap}B}{A{\cup}B})*100\%$

where, $A$ and $B$ are the sets of most relevant items. We specifically
decided to use ‘most’ relevant (top-$n$ items) as the tail of the distribution
mostly consists of values close to $0$.

#### Interpretation methods under study:

In this paper we study interpretation stability using the following three
interpretation methods:

  1. 1.

Attention based interpretation: We focus on attention probabilities as the
mode of interpretation and consider the model to be stable if different
instantiations of the model leads to similar attention distributions. Our
major focus in this paper is attention based interpretation. As we use Jain et
al. (2019) as a testbed for our investigation, we focus heavily on attention.
Also, as the attention layer has a linear relationship with the prediction, we
consider attention to be more indicative of the model stability.

  2. 2.

Gradient-based feature importance: Given a sample, we use the input gradients
of the model corresponding to each of the word representations and compute the
magnitude of the change as a local explanation. We refer the reader to
Baehrens et al. (2010) for a good introduction to gradient-based
interpretations. As all of our models are differentiable, we use this as an
alternative method for interpretation. We follow the standard procedure as
followed in Feng et al. (2018) and note that we do not follow Jain and Wallace
(2019) and do not disconnect the computational graph at the attention module.
We obtain probabilistic gradient scores by normalizing over the absolute
values of gradient values.

  3. 3.

LIME based interpretation: We use locally interpretable model-agnostic
interpretations Ribeiro et al. (2016) that learns a surrogate interpretable
model locally around the predictions of the deep neural based model. We obtain
LIME based interpretations for every instantiation of the models. We then use
Jaccard Distance to measure the divergence.

We note that, we observe similar patterns across the three interpretation
methods and the interpretations consistently differ with random seeds.

##  3 Reducing Model Instability with an Optimization Lens

We observe that different instantiations of the model can cause the model have
different starts on the optimization surface. Further, stochastic sampling
might result in different paths. Both of these factors can lead to different
local minimas potentially leading to different solutions. With this
observation as our background we propose two, closely related, methods to
ameliorate divergence: Agressive Stochastic Weight Averaging and Norm-filtered
Agressive Stochastic Weight Averaging. We describe these two in the following
subsections.

###  3.1 Aggressive Stochastic Weight Averaging (ASWA)

Stochastic weight averaging (SWA) Izmailov et al. (2018) works by averaging
the weights of multiple points in the trajectory of gradient descent based
optimizers. The algorithm typically uses modified learning rate schedules. SWA
is itself based on the idea of maintaining a running average of weights in
stochastic gradient descent based optimization techniques Ruppert (1988);
Polyak and Juditsky (1992). The principle idea in SWA is averaging the weights
that are maximally distant helps stabilize the gradient descent based
optimizer trajectory and improves generalization. Izmailov et al. (2018) use
the analysis of Mandt et al. (2017) to illustrate the stability arguments
where they show that, under certain convexity assumptions, SGD iterations can
be visualized as sampling from a Gaussian distribution centred at the _minima_
of the loss function. Samples from high-dimensional Gaussians are expected to
be concentrated _on the surface of the ellipse_ and not close to the _mean_.
Averaging iterations is shown to stabilize the trajectory and further improve
the width of the solutions to be closer to the _mean_.

In this paper, we focus on the stability of deep neural models as a function
of random-seeds. Our proposal is based on SWA, but we extend it to the
extremes and call it _Aggressive Stochastic Weight Averaging_. We assume that,
for small batch size, the loss surface is locally convex. We further relax the
conditions for the optimizer and assume that the optimizer is based on some
version of gradient descent — this means that our modification is valid even
for other pseudo-first-order optimization algorithms including Adam Kingma and
Ba (2014) and Adagrad Duchi et al. (2011).

We note that, Izmailov et al. (2018) suggest using SWA usually after
’pre-’training the model (at least until $75\%$ convergence) and followed by
sampling weights at different steps either using large constant or cyclical
learning rates. While, SWA is well defined for convex losses Polyak and
Juditsky (1992), Izmailov et al. (2018) connect SWA to non-convex losses by
suggesting that the loss surface is  _approximately_ convex after convergence.
In our setup, we investigate the utility of averaging weights over every
iteration (an iteration consists of one batch of the gradient descent).
Algorithm 1 shows the implementation pseudo-code for SWA. We note that, unlike
Izmailov et al. (2018), we average our weights at each batch update and assign
the ASWA parameters to the model at the end of each epoch. That is, we replace
the model’s weights for the next epoch with the averaged weights.

1:

2:$e=$ Epoch number

3:$m=$ Total epochs

4:$i=$ Iteration number

5:$n=$ Total iterations

6:$\alpha=$ Learning rate

7:$\mathcal{O}=$ Stochastic Gradient optimizer function

$e\leftarrow 0$;

while  _$e <m$_ do

$i\leftarrow 1$

while  _$i\leq n$_ do

$W_{swa}\leftarrow W_{swa}+\frac{(W-W_{swa})}{(e*n+i+1)}$;

$W\leftarrow W-\mathcal{O}(\alpha,W)$;

$i\leftarrow i+1$

$W\leftarrow W_{swa}$;

$e\leftarrow e+1$

Algorithm 1 Aggressive SWA algorithm

In Figure 2, we show an SGD optimizer (with momentum) and the same optimizer
_with SWA_ over a $3$-dimensional loss surface with a saddle point. We observe
that the original SGD reaches the desired minima, however, it almost reaches
the saddle point and does a course correction and reaches minima. On the other
hand, we observe that SGD with ASWA is very conservative, it repeatedly
restarts and reaches the minima without reaching the saddle point. We
empirically observe that this is a desired property for the stability of
models over runs of the same model that differ only over random
instantiations. The grey circles in Figure 2 highlight this conservative
behaviour of SGD with ASWA optimizer, especially when compared to the standard
SGD. Further, Polyak and Juditsky (1992) show that for convex losses,
averaging SGD proposals achieves the highest possible rate of convergence for
a variety of first-order SGD based algorithms.

(a) Trajectory for Stochastic Gradient Descent

(b) Trajectory for Stochastic Gradient Descent with ASWA

Figure 2: Trajectory for gradient descent algorithms with red and black arrows
on (b) indicating movements from consecutive epochs with restarts.
Conservative behaviour of ASWA algorithm helps avoid the saddle point without
ever reaching it.

###  3.2 Norm-filtered Aggressive Stochastic Weight Averaging (NASWA)

1:

2:$e=$ Epoch number

3:$m=$ Total epochs

4:$i=$ Iteration number

5:$n=$ Total iterations

6:$\alpha=$ Learning rate

7:$\mathcal{O}=$ Stochastic Gradient optimizer function

8:$N_{s}=$ List of previous iterations’ norm differences

$e\leftarrow 0$;

while  _$e <m$_ do

$i\leftarrow 1$

while  _$i\leq n$_ do

$N_{cur}\leftarrow\|W-W_{swa}\|_{1}$;

$N_{mean}\leftarrow\frac{\sum_{i=1}^{|N_{s}|}N_{s}[i]}{|N_{s}|}$;

if  _$N_{cur} >N_{mean}$_ then

$W_{swa}\leftarrow W_{swa}+\frac{(W-W_{swa})}{(e*n+i+1)}$;

$N_{s}\leftarrow[N_{cur}]$;

else

$N_{s}\leftarrow N_{s}+[N_{cur}]$;

$W\leftarrow W-\mathcal{O}(\alpha,W)$;

$i\leftarrow i+1$

$W\leftarrow W_{swa}$;

$e\leftarrow e+1$

Algorithm 2 Norm-filtered Aggressive SWA algorithm

We observe that the ASWA algorithm is especially beneficial when the norm
difference of the parameters of the model are high. We hypothesise that in
general, the norm difference indicates the divergence between optimizers’
steps and we observe that the larger the norm difference, the greater the
change in the trajectory. Therefore, we propose to maintain a list that stores
the norm differences of the previous iterations. If the norm difference of the
current iteration is greater than the average of the list, we update the ASWA
weights and reinitialize the list with the current norm difference. When the
norm difference, however, is less than the average of the list, we just append
the current norm difference to the list. After the completion of the epoch, we
assign the ASWA parameters to the model. This is shown in Algorithm 2. We call
this approach _Norm-filtered Aggressive Stochastic Weight Averaging_.

##  4 Experiments

We base our investigation on similar sets of models as Jain and Wallace
(2019). We also use the code provided by the authors for our empirical
investigations for consistency and empirical validation. We describe our
models and datasets used for the experiments below.

###  4.1 Models

We consider two sets of commonly used neural models for the tasks of binary
classification and multi-class natural language inference. We use CNN and bi-
directional LSTM based models with attention. We follow Jain and Wallace
(2019) and use similar attention mechanisms using a) additive attention
Bahdanau et al. (2014); and b) scaled dot product based attention Vaswani et
al. (2017). We jointly optimize all the parameters for the model, unlike Jain
and Wallace (2019) where the encoding layer, attention layer and the output
prediction layer are all optimized separately. We experiment with several
optimizers including Adam Kingma and Ba (2014), SGD and Adagrad Duchi et al.
(2011) but most results below are with Adam.

For our ASWA and NASWA based experiments, we use a constant learning rate for
our optimizer. Other model-specific settings are kept the same as Jain and
Wallace (2019) for consistency.

Dataset | Avg. Length | Train Size | Test size
---|---|---|---
IMDB | 179 | 12500 / 12500 | 2184 / 2172
Diabetes(MIMIC) | 1858 | 6381 / 1353 | 1295 / 319
SST | 19 | 3034 / 3321 | 652/653
Anemia(MIMIC) | 2188 | 1847 / 3251 | 460 / 802
AgNews | 36 | 30000 / 30000 | 1900 / 1900
ADR Tweets | 20 | 14446 / 1939 | 3636 / 487
SNLI | 14 | 182764 / 183187 / 183416 | 3219 / 3237 / 3368

Table 1: Dataset characteristics. Train size and test size show the
cardinality for each class. SNLI is a three-class dataset while the rest are
binary classification

###  4.2 Datasets

The datasets used in our experiments are listed in Table 1 with summary
statistics. We further pre-process and tokenize the datasets using the
standard procedure and follow Jain and Wallace (2019). We note that IMDB Maas
et al. (2011), Diabetes(MIMIC) Johnson et al. (2016), Anemia(MIMIC) Johnson et
al. (2016), AgNews Zhang et al. (2015), ADR Tweets Nikfarjam et al. (2015) and
SST Socher et al. (2013) are datasets for the binary classification setup.
SNLI Bowman et al. (2015) is a dataset for the multiclass classification
setup. All of the datasets are in English, however we expect the behavior to
persist regardless of the language.

###  4.3 Settings and Hyperparameters

We use a $300$-dimenstional embedding layer which is initialized with FastText
Joulin et al. (2016) based free-trained embeddings for both CNN and the bi-
directional LSTM based models. We use a $128$-dimensional hidden layer for the
bi-directional LSTM and a $32$-dimensional filter with kernels of size
$\\{1,3,5,7\\}$ for CNN. For others, we maintain the model settings to
resemble the models in Jain and Wallace (2019). We train all of our models for
20 Epochs with a constant batch size of 32. We use early stopping based on the
validation set using task-specific metrics (Binary Classification: using roc-
auc, Multiclass and question answering based dataset: using accuracy).

Dataset | CNN(%) | CNN+ASWA(%) | CNN+NASWA(%)
---|---|---|---
IMDB | 89.8 ($\pm$0.79) | 90.2 ($\pm$0.25) | 90.1 ($\pm$0.29)
Diabetes | 87.4 ($\pm$2.26) | 85.9 ($\pm$0.25) | 85.9 ($\pm$0.38)
SST | 82.0 ($\pm$1.01) | 82.5 ($\pm$0.39) | 82.5 ($\pm$0.39)
Anemia | 90.6 ($\pm$0.98) | 91.9 ($\pm$0.20) | 91.9 ($\pm$0.19)
AgNews | 95.5 ($\pm$0.23) | 96.0 ($\pm$0.11) | 96.0 ($\pm$0.07)
Tweet | 84.6 ($\pm$2.65) | 84.4 ($\pm$0.54) | 84.4 ($\pm$0.54)

Table 2: Performance statistics obtained from 10 differently seeded CNN based
models. Table compares accuracy and its standard deviation for the normally
trained CNN model against the ASWA and NASWA trained models, whose deviation
drops significantly, thus, indicating increased robustness.

##  5 Results

In this section, we summarize our findings for $10$ runs of the model with
$10$ different random seeds but with identical model settings.

###  5.1 Model Performance and Stability

We first report model performance and prediction stability. The results are
reported in Table 2.

Dataset | LSTM(%) | LSTM+ASWA(%) | LSTM+NASWA(%)
---|---|---|---
IMDB | 89.1 ($\pm$1.34) | 90.2 ($\pm$0.32) | 90.3 ($\pm$0.17)
Diabetes | 87.7 ($\pm$1.44) | 87.7 ($\pm$0.60) | 87.8 ($\pm$0.55)
SST | 81.9 ($\pm$1.11) | 82.0 ($\pm$0.60) | 82.1 ($\pm$0.57)
Anemia | 91.6 ($\pm$0.49) | 91.8 ($\pm$0.34) | 91.9 ($\pm$0.36)
AgNews | 95.5 ($\pm$0.32) | 96.1($\pm$0.17) | 96.1 ($\pm$0.10)
Tweet | 84.7 ($\pm$1.79) | 83.8 ($\pm$0.45) | 83.9 ($\pm$0.45)

Table 3: Performance statistics obtained from 10 differently seeded LSTM based
models.

We note that the original CNN based models, on an average, have a standard
deviation of $\pm 1.5\%$. Which seems standard, however, we note that ADR
Tweets dataset has a very high standard deviation of $\pm 2.65\%$. We observe
that ASWA and NASWA are almost always able to achieve higher performance with
a very low standard deviation. This suggests that both ASWA and NASWA are
extremely stable when compared to the standard model. They significantly
improve the robustness, on an average, by $72\%$ relative to the original
model and on Diabetes (MIMIC), a binary classification dataset, by $89\%$
(relative improvement). We observe similar results for the LSTM based models
in Table 3.

(a) CNN models

(b) LSTM models

Figure 3: Prediction’s standard deviation for CNN and LSTM based models for
all binary classification datasets under consideration. Predictions are
bucketed in intervals of size 0.1, starting from 0 (containing predictions
from 0 to 0.1), until 0.9

We further analyze the prediction score stability by computing the mean
standard deviation over the binned confidence intervals of the models in
Figure 3(a). We note that on an average, the standard deviations are on the
lower side. However, we observe that the mean standard deviation of the bins
close to $0.5$ is on the higher side as is expected given the high
uncertainty. On the other hand both, ASWA and NASWA based models are
relatively more stable than the standard CNN based model. We observe similar
behaviours for the LSTM based models in Figure 3(b). This suggests that our
proposed methods, ASWA and NASWA, are able to obtain relatively better
stability without any loss in performance. We also note that both ASWA and
NASWA had relatively similar performance over more than $10$ random seeds.

###  5.2 Attention Stability

(a) CNN models

(b) LSTM models

Figure 4: Average attention entropy against the bucketed predictions for CNN
and LSTM based models. Figure highlights the high entropy between attention
based distributions from differently seeded models (especially for the
Diabetes-MIMIC datatset), indicating towards model instability.

We now consider the stability of attention distributions as a function of
random seeds. We first plot the results of the experiments for _standard_ CNN
based binary classification models over uniformly binned prediction scores for
positive labels in Figure 4(a). We observe that, depending on the datasets,
the attention distributions can become extremely unstable (high entropy). We
specifically highlight the Diabetes(MIMIC) dataset’s entropy distribution. We
observe similar, but relatively worse results for the LSTM based models in
Figure 4(b). In general, we would expect the entropy distribution to be close
to zero however, this doesn’t seem to be the case. This means that using
attention distributions to interpret models may not be reliable and can lead
to misinterpretations.

(a) CNN models

(b) LSTM models

Figure 5: Jaccard distance highlighting instability in attention distributions
of CNN and LSTM based models.

(a) CNN+ASWA

(b) CNN+NASWA

(c) LSTM+ASWA

(d) LSTM+NASWA

Figure 6: Improved prediction stability from ASWA and NASWA for CNN and LSTM
based models

(a) Diabetes

(b) SST

Figure 7: Gradient based interpretations’ stability improvement from NASWA on
CNN based models. The Jaccard distance is calculated using the top 20%
attentive items.

(a) IMDB

(b) Diabetes

(c) SST

(d) Anemia

(e) AgNews

(f) ADR Tweets

Figure 8: Attention stability improvement from ASWA and NASWA on CNN based
models.

(a) Label 0 prediction vs entropy

(b) Label 1 prediction vs entropy

(c) Label 2 prediction vs entropy

Figure 9: Attention stability improvement from ASWA and NASWA on CNN based
model for the SNLI dataset.

We use the top $20\%$ of the most important items (indices) in the attention
distribution for each dataset over $10$ runs and plot the Jaccard distances
for CNN and LSTM based models in Figure 5(a) and Figure 5(b). We again notice
a similar trend of unstable attention distributions over both CNN and LSTM
based attention distribution.

In the following sections for space constraints, we focus on CNN based models
with additive attention. Our results on LSTM based models are provided in the
attached supplementary material. We note that the observations for LSTM models
are, in most cases, similar to the behaviour of the CNN based models. Scaled
dot-product based models are also provided in the supplementary material and
we notice a similar trend as the additive attention.

We now focus on the effect of ASWA and NASWA on binary and multi-class CNN
based neural models separately.

#### Binary Classification

In Figure 8, we plot the results of the models with ASWA and NASWA. We observe
that both these algorithms significantly improve the model stability and
decrease the entropy between attention distributions. For example, in Figure
8(b), both ASWA and NASWA decrease the average entropy by about $60\%$. We
further notice that NASWA is slightly better performing in most of the runs.
This empirically validates the hypothesis that averaging the weights from
divergent weights (when the norm difference is higher than the average norm
difference) helps in stabilizing the model’s parameters, resulting in a more
robust model.

#### Multi-class Classification

In Figure 9, we plot the entropy between the attentions distributions of the
models for the SNLI dataset (CNN based model), separately for _each label_
(_neutral_ , _contradiction_ , and _entailment_). We notice, similar
observations as the binary classification models, the ASWA and NASWA
algorithms are able to significantly improve the entropy of the attention
distributions and increases the robustness of the model with random seeds.

###  5.3 Gradient-based Interpretations

We now look at an alternative method of interpreting deep neural models and
look into the consistency of the gradient-based interpretations to further
analyze the model’s instability. For this setup, we focus on binary classifier
and plot the results on the SST and the Diabetes dataset in particular since
they cover the low and the high end of the entropy spectrum (respectively). We
notice similar trends of instability in the gradient-based interpretations
from model inputs as we did for the attention distributions. Figure 7 shows
that the entropy between the gradient-based interpretations from differently
seeded models closely follows the same trend as the attention distributions.
This result further strengthens our claim on the importance of model stability
and shows that over different runs of the same model with different seeds, we
may get different interpretations using gradient-based feature importance.
Moreover, Figure 7 shows the impact of ASWA towards making the gradient-based
interpretations more consistent, thus, significantly increasing the stability.

###  5.4 LIME based Interpretations

We further evaluated the surrogate model based interpretability using LIME
Ribeiro et al. (2016). LIME obtains a locally linear approximation of the
model’s behaviour for a given sample by perturbing it and learning a sparse
linear model around it. We focus on AgNews and SST based datasets and obtain
interpretability estimates using LIME. Once again, we notice a similar pattern
of instability as the other two interpretability methods. In Figure 10 we
present our results from the LIME based interpretations with Jaccard distance
as the measure. Note that we measure the Jaccard distance over the top $20\%$
most influential items. We observe once again that NASWA helps in reducing the
instability and results in more consistent interpretations.

(a) AgNews

(b) SST

Figure 10: LIME based interpretations’ stability improvement from NASWA on CNN
based models. The Jaccard distance is calculated using the top 20% attentive
items.

In all our experiments, we find that a significant proportion of
interpretations are dependent on the instantiation of the model. We also note
that we perform experiments over $100$ random seeds for greater statistical
power and see similar patterns333These results are provided in the appendix..

##  6 Discussion

Recent advances in adversarial machine learning Neelakantan et al. (2015);
Zahavy et al. (2016) have investigated robustness to random initialization
based perturbations, however, to our knowledge, no previous study investigates
the effect of random-seeds and its connection on model interpretation. Our
study analyzed the inherent lack of robustness in deep neural models for NLP.
Recent studies cast doubt on the consistency and correlations of several types
of interpretations Doshi-Velez and Kim (2017); Jain and Wallace (2019); Feng
et al. (2018). We hypothesise that some of these issues are due to the
inherent instability of the deep neural models to random-seed base
perturbations. Our analysis (in Section 4) leads to the hypothesis that models
with different instantiations may use completely different optimization paths.
The issue of variance in all black-box interpretation methods over different
seeds will continue to persist until the models are fully robust to random-
seed based perturbations. Our work however, doesn’t provide insights into
instabilities of different layers of the models. We hypothesise that it might
further uncover the reasons for the relatively lower correlation between
different black-box interpretation methods as these are effectively based off
on different layers and granularity.

There has been some work on using noisy gradients Neelakantan et al. (2015)
and learning from adversarial and counter-factual examples Feng et al. (2018)
to increase the robustness of deep learning models. Feng et al. (2018) show
that neural models may use redundant features for prediction and also show
that most of the black-box interpretation methods may not be able to capture
these second-order effects. Our proposals show that aggressively averaging
weights leads to better optimization and the resultant models are more robust
to random-seed based perturbation. However, our research is limited to
increasing consistency in neural models. Our approach further uses first order
based signals to boost stability. We posit that second-order based signals can
further enhance consistency and increase the robustness.

##  7 Conclusions

In this paper, we study the inherent instability of deep neural models in NLP
as a function of random seed. We analyze model performance and robustness of
the model in the form of attention based interpretations, gradient-based
feature importance and LIME based interpretations across multiple runs of the
models with different random seeds. Our analysis strongly highlights the
problems with stability of models and its effects on black-box interpretation
methods leading to different interpretations for different random seeds. We
also propose a solution that makes use of weight averaging based optimization
technique and further extend it with norm-filtering. We show that our proposed
methods largely stabilize the model to random-seed based perturbations and, on
average, significantly reduce the standard deviations of the model performance
by $72\%$. We further show that our methods significantly reduce the entropy
in the attention distribution, the gradient-based feature importance measures
and LIME based interpretations across runs.

## Acknowledgments

We thank Panos Parpas and Emtiyaz Khan for their feedback on an earlier draft
of this paper. We thank the anonymous reviewers for their thorough reviews and
constructive comments. Pranava Madhyastha kindly acknowledges the support of
Amazon AWS Cloud Credits for Research Award, hardware grant from NVIDIA, Anne
O’Neill and the Imperial Corporate Partnership Programme.

## References

  * Alvarez-Melis and Jaakkola (2017) David Alvarez-Melis and Tommi S Jaakkola. 2017.  A causal framework for explaining the predictions of black-box sequence-to-sequence models.  _arXiv preprint arXiv:1707.01943_.
  * Bach et al. (2015) Sebastian Bach, Alexander Binder, Grégoire Montavon, Frederick Klauschen, Klaus-Robert Müller, and Wojciech Samek. 2015.  On pixel-wise explanations for non-linear classifier decisions by layer-wise relevance propagation.  _PloS one_ , 10(7):e0130140.
  * Baehrens et al. (2010) David Baehrens, Timon Schroeter, Stefan Harmeling, Motoaki Kawanabe, Katja Hansen, and Klaus-Robert MÃžller. 2010.  How to explain individual classification decisions.  _Journal of Machine Learning Research_ , 11(Jun):1803–1831.
  * Bahdanau et al. (2014) Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio. 2014.  Neural machine translation by jointly learning to align and translate.  _arXiv preprint arXiv:1409.0473_.
  * Bengio (2012) Yoshua Bengio. 2012.  Practical recommendations for gradient-based training of deep architectures.  In _Neural networks: Tricks of the trade_ , pages 437–478. Springer.
  * Bowman et al. (2015) Samuel R Bowman, Gabor Angeli, Christopher Potts, and Christopher D Manning. 2015\.  A large annotated corpus for learning natural language inference.  _arXiv preprint arXiv:1508.05326_.
  * Doshi-Velez and Kim (2017) Finale Doshi-Velez and Been Kim. 2017.  Towards a rigorous science of interpretable machine learning.  _arXiv preprint arXiv:1702.08608_.
  * Duchi et al. (2011) John Duchi, Elad Hazan, and Yoram Singer. 2011.  Adaptive subgradient methods for online learning and stochastic optimization.  _Journal of Machine Learning Research_ , 12(Jul):2121–2159.
  * Feng et al. (2018) Shi Feng, Eric Wallace, Alvin Grissom II, Pedro Rodriguez, Mohit Iyyer, and Jordan Boyd-Graber. 2018.  Pathologies of neural models make interpretation difficult.  In _Empirical Methods in Natural Language Processing_.
  * Ghaeini et al. (2018) Reza Ghaeini, Xiaoli Z Fern, and Prasad Tadepalli. 2018.  Interpreting recurrent and attention-based neural models: a case study on natural language inference.  _arXiv preprint arXiv:1808.03894_.
  * Izmailov et al. (2018) Pavel Izmailov, Dmitrii Podoprikhin, Timur Garipov, Dmitry Vetrov, and Andrew Gordon Wilson. 2018.  Averaging weights leads to wider optima and better generalization.  _arXiv preprint arXiv:1803.05407_.
  * Jain et al. (2019) Sarthak Jain, Ramin Mohammadi, and Byron C Wallace. 2019.  An analysis of attention over clinical notes for predictive tasks.  _arXiv preprint arXiv:1904.03244_.
  * Jain and Wallace (2019) Sarthak Jain and Byron C. Wallace. 2019.  Attention is not explanation.  _CoRR_ , abs/1902.10186.
  * Johnson et al. (2016) Alistair EW Johnson, Tom J Pollard, Lu Shen, H Lehman Li-wei, Mengling Feng, Mohammad Ghassemi, Benjamin Moody, Peter Szolovits, Leo Anthony Celi, and Roger G Mark. 2016.  Mimic-iii, a freely accessible critical care database.  _Scientific data_ , 3:160035.
  * Joulin et al. (2016) Armand Joulin, Edouard Grave, Piotr Bojanowski, Matthijs Douze, Hérve Jégou, and Tomas Mikolov. 2016.  Fasttext. zip: Compressing text classification models.  _arXiv preprint arXiv:1612.03651_.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_.
  * Koh and Liang (2017) Pang Wei Koh and Percy Liang. 2017.  Understanding black-box predictions via influence functions.  In _Proceedings of the 34th International Conference on Machine Learning-Volume 70_ , pages 1885–1894. JMLR. org.
  * Maas et al. (2011) Andrew L Maas, Raymond E Daly, Peter T Pham, Dan Huang, Andrew Y Ng, and Christopher Potts. 2011.  Learning word vectors for sentiment analysis.  In _Proceedings of the 49th annual meeting of the association for computational linguistics: Human language technologies-volume 1_ , pages 142–150. Association for Computational Linguistics.
  * Mandt et al. (2017) Stephan Mandt, Matthew D Hoffman, and David M Blei. 2017.  Stochastic gradient descent as approximate bayesian inference.  _The Journal of Machine Learning Research_ , 18(1):4873–4907.
  * Neelakantan et al. (2015) Arvind Neelakantan, Luke Vilnis, Quoc V Le, Ilya Sutskever, Lukasz Kaiser, Karol Kurach, and James Martens. 2015.  Adding gradient noise improves learning for very deep networks.  _arXiv preprint arXiv:1511.06807_.
  * Nikfarjam et al. (2015) Azadeh Nikfarjam, Abeed Sarker, Karen O’connor, Rachel Ginn, and Graciela Gonzalez. 2015.  Pharmacovigilance from social media: mining adverse drug reaction mentions using sequence labeling with word embedding cluster features.  _Journal of the American Medical Informatics Association_ , 22(3):671–681.
  * Polyak and Juditsky (1992) Boris T Polyak and Anatoli B Juditsky. 1992.  Acceleration of stochastic approximation by averaging.  _SIAM Journal on Control and Optimization_ , 30(4):838–855.
  * Ribeiro et al. (2016) Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. 2016.  Model-agnostic interpretability of machine learning.  _arXiv preprint arXiv:1606.05386_.
  * Ruppert (1988) David Ruppert. 1988.  Stochastic approximation.  Technical report, Cornell University Operations Research and Industrial Engineering.
  * Socher et al. (2013) Richard Socher, Alex Perelygin, Jean Wu, Jason Chuang, Christopher D Manning, Andrew Ng, and Christopher Potts. 2013.  Recursive deep models for semantic compositionality over a sentiment treebank.  In _Proceedings of the 2013 conference on empirical methods in natural language processing_ , pages 1631–1642.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  In _Advances in neural information processing systems_ , pages 5998–6008.
  * Zahavy et al. (2016) Tom Zahavy, Bingyi Kang, Alex Sivak, Jiashi Feng, Huan Xu, and Shie Mannor. 2016\.  Ensemble robustness and generalization of stochastic deep learning algorithms.  _arXiv preprint arXiv:1602.02389_.
  * Zhang et al. (2015) Xiang Zhang, Junbo Zhao, and Yann LeCun. 2015.  Character-level convolutional networks for text classification.  In _Advances in neural information processing systems_ , pages 649–657.

Figure 11: Entropy improvement for tanh Attention based CNN model for the SST
dataset using 100 different seeds.

Figure 12: Entropy improvement for tanh Attention based CNN model for the
AgNews dataset using 100 different seeds.

Figure 13: Entropy improvement for tanh Attention based CNN model for the IMDB
dataset using 100 different seeds.

(a) IMDB

(b) Diabetes

(c) SST

(d) Anemia

(e) AgNews

(f) ADR Tweets

Figure 14: Attention stability improvement from ASWA and NASWA on LSTM based
models.

##  Appendix A Jaccard Distance Experiments

In Figure 5(a) and Figure 5(b) we plot the Jaccard distance plots with CNN
models and LSTM models and note that ASWA consistently improves the stability
of the interpretations.

Figure 15: Jaccard Distance improvement for Diabetes.

Figure 16: Jaccard Distance improvement for Diabetes.

Figure 17: Entropy improvement for dot Attention based CNN model for SST
dataset.

Figure 18: Entropy improvement for dot Attention based CNN model for Diabetes
dataset.

##  Appendix B Results with 100 seeds

We perform attention stability based experiments as mentioned in the paper,
but now with 100 seeds, i.e, 100 models initialized with different seeds
instead of 10. Figures 11, 13, and 12 show the entropy of attention based
interpretations for different datasets. Experimenting with 100 seeds helps
strengthen our claims about the instability of the model and the effectiveness
of our proposed algorithms like ASWA and NASWA.

##  Appendix C Hyperparameter Settings

For training purposes, we use the same model settings for the models as
mentioned in the paper Jain and Wallace (2019) (or the Github
implementation444https://github.com/successar/AttentionExplanation), our port
of the code is made available at: https://github.com/rishj97/ModelStability.
Additional hyper-parameters for replication studies are:

  * •

Number of epochs: 20

  * •

Optimizer: Adam

  * •

Learning rate: 0.001

The exact seeds used for running the experiments can be found in our code
repository.

##  Appendix D Binary Classification with LSTM based Models

For LSTM based models, we notice (in Figure 14) similar trends as to the CNN
models in terms of the instability of the attention based interpretations.


# Complex topological features of reservoirs shape learning performances in
bio-inspired recurrent neural networks

Valeria d’Andrea<EMAIL_ADDRESS>Fondazione Bruno Kessler, Via Sommarive 18,
38123 Povo (TN), Italy  Michele Puppin  Dipartimento di Fisica e Astronomia
”Galileo Galilei”, Università di Padova, Via Marzolo 8, 35131 Padova, Italy
Manlio De Domenico<EMAIL_ADDRESS>Dipartimento di Fisica e
Astronomia ”Galileo Galilei”, Università di Padova, Via Marzolo 8, 35131
Padova, Italy

###### Abstract

Recurrent networks are a special class of artificial neural systems that use
their internal states to perform computing tasks for machine learning. One of
its state-of-the-art developments, i.e. reservoir computing (RC), uses the
internal structure – usually a static network with random structure – to map
an input signal into a nonlinear dynamical system defined in a higher
dimensional space. Reservoirs are characterized by nonlinear interactions
among their units and their ability to store information through recurrent
loops, allowing to train artificial systems to learn task-specific dynamics.
However, it is fundamentally unknown how the random topology of the reservoir
affects the learning performance. Here, we fill this gap by considering a
battery of synthetic networks – characterized by different topological
features – and 45 empirical connectomes – sampled from brain regions of
organisms belonging to 8 different species – to build the reservoir and
testing the learning performance against a prediction task with a variety of
complex input signals. We find nontrivial correlations between RC performances
and both the number of nodes and rank of the covariance matrix of activation
states, with performance depending on the nature – stochastic or deterministic
– of input signals. Remarkably, the modularity and the link density of the
reservoir are found to affect RC performances: these results cannot be
predicted by models only accounting for simple topological features of the
reservoir. Overall, our findings highlight that the complex topological
features characterizing biophysical computing systems such as connectomes can
be used to design efficient bio-inspired artificial neural networks.

Brain networks $|$ Bio-inspired Reservoir Computing $|$ Complex Networks

## Introduction

The physical and biological basis of cognition are embedded in the complex
temporally and spatially multiscale structure of the brain Bullmore and Sporns
[2012], Bassett and Gazzaniga [2011], de Vico Fallani et al. [2014], Presigny
and Fallani [2022]. Even if it is not yet fully understood how cognitive
functions can emerge from brain activity, it has been demonstrated that
structural properties of human brain networks play a key role in brain
functioning. For instance, in healthy individuals the configuration of brain
connectivity is built upon different aspects, such as behavioral variability
Bassett et al. [2009], cognitive ability Li et al. [2009], van den Heuvel et
al. [2009], age Meunier et al. [2008], Sifis et al. [2009] and gender Gong et
al. [2009]. Connectivity patterns exhibit plasticity: they are inherently
dynamical and reshaped as a consequence of experimental tasks Bassett et al.
[2007], De Vico Fallani et al. [2008] and drug treatment Achard and Bullmore
[2007], Schwarz et al. [2009]. Similarly, individuals affected by psychiatric
and neurodegenerative disorders present abnormal structural and functional
connectivity between brain regions Catani et al. [2005]. Some diseases
characterized by abnormal cortical structural/functional organization are
schizophrenia Lynall et al. [2010], De Domenico et al. [2016], Alzheimer’s
disease Stam et al. [2008], Guillon et al. [2017], epilepsy Horstmann et al.
[2010] and many others.

The characterization of the human brain architecture requires also
investigation from a dynamical perspective Deco et al. [2011], Churchland et
al. [2010], which has been successfully investigated with methods from
statistical physics Burioni et al. [2014], di Santo et al. [2018], Munoz
[2018]. In the quest to build computational models resembling the most salient
feature of the human brain, Reservoir Computing (RC) Lukoševičius and Jaeger
[2009], Gauthier et al. [2021] provides a suitable framework for its
transparency and its versatility for implementation on physical devices Paquot
et al. [2012], Van der Sande et al. [2017], Tanaka et al. [2019a]. RC
identifies a class of algorithms that are biologically-inspired computing
systems composed by a set of processing units called neurons, organized within
a network. A neuron receives a signal, processes it and propagates the output
signal to the connected neurons. The strength of the connection, analogously
to the synaptic strength, encodes how much a neuron influences another neuron.
By adapting network structure, i.e. finding the optimal set of link weights,
the machine can learn to perform desired tasks and to develop meaningful
internal representations of the environment. RC algorithms represent the
state-of-the-art algorithms to tackle tasks such as images recognition, text
comprehension, natural language processing, game playing and many others. In
these tasks, RC algorithms have reached super-human performances Storrs and
Kriegeskorte [2019], Murugan [2018].

From a network neuroscience point of view, RC-based approaches can be used to
gain new insights about the role of brain connectivity patterns in the
execution of complex cognitive tasks. In fact, at variance with other Machine
Learning methods, RC approaches are characterized by biological plausibility,
with several architectural and dynamical properties of mammalian brains used
in RC design Lukoševičius and Jaeger [2009]. RC can provide insights in how
brains perform accurate computations with an inaccurate physical scaffolding
Buonomano and Merzenich [1995], Häusler and Maass [2007] and it helps
explaining brain time encoding Karmarkar and Buonomano [2007], visual
information processing in primary visual cortex Stanley et al. [1999], Nikolić
et al. [2006] and the representation of sequential information. In addition,
RC offers a functional interpretation of the cerebellar circuitry Kistler and
De Zeeuw [2002], Yamazaki and Tanaka [2007]. The structure of RC networks
plays a fundamental role, allowing to approximate very complex functions and
it has been shown that real-world inspired connectivity such as Watts-Strogatz
and Barabási-Albert models has resulted in competitive performances compared
to optimized state-of-the-art architectures Xie et al. [2019]. Furthermore,
during training feedforward networks have proved to spontaneously form non-
random topologies also found in brains, such as modular structures Filan et
al. [2020]. This and other evidences support the notion that non-random
topologies can lead to desired performance of ANNs Damicelli et al. [2021].

In our work we probe the inverse approach by using real-world connectomes as
reservoirs, in order to investigate how RC performance in a prediction task is
affected by topological features occurring in empirical network structures. We
also investigate the dynamical input representation by studying the
dimensionality of elicited connectome internal states.

To verify if the RC behavior is affected by the topological features
identified in real-world connectomes, we test our results on three null models
of increasing complexity, namely an Erdős Rényi (ER) model, Configuration
Model (CM) and a Stochastic Block Model (SBM).

Finally, we test our results on multiple signals to verify if the RC
performance depends on the stochastic or deterministic nature of the input
signal.

## Results

In this study, we collected and use 45 connectomes providing a network
representation of the neuronal system of C. elegans, C. intestinalis,
Drosophila, Human, Mouse, Platynereis dumerilly, Rattus norvegicus and Rhesus
macaque in terms of adjacency matrices $A_{conn}$. For each connectome, we use
the corresponding adjacency matrix to build a reservoir.

Figure 1: Allometric scaling in empirical connectomes. Scatter plot describing
the number of neuronal links, $E$, as a function of the number of nodes, $N$,
for the connectome data sets ($n=45$ networks) considered in this study. Each
point is a single network, its color identifies the species and its size is
proportional to the network modularity (mod).

For each empirical network we calculate some topological indicators, such as
average degree and modularity, which are widely used to characterize network
structures and their mesoscale organization. Figure 1 shows that connectome
data span more than 2 orders of magnitude in the number of nodes and more than
3 orders of magnitude in number of edges; a complete list of connectomes is
shown in Tab. 2 along with the number of nodes and the average degree of the
network.

The performances of the connectome-based reservoirs are tested in a one-step-
ahead prediction task, where given a point of a temporal sequence as input the
model is asked to predict the next point of the sequence. To probe the ability
of our bio-inspired machines to deal with signals of different nature, we
consider four different time series as inputs: two stochastic signals (white
noise and fractional Brownian noise) and two deterministic chaotic signals
(Ikeda and Mackey-Glass signals).

Organism | Num. of | Num. of | Brain mass | Brain volume
---|---|---|---|---
| neurons | synapses |  |
|  |  | $(\mathbf{\text{g}})$ | $(\mathbf{\text{cm}^{3}})$
C. intestinalis | $213$ | $3105$ | - | -
C. elegans | $302$ | $7500$ | - | 1
Platynereis dumerilli | $2728$ | $1.14\cdot 10^{4}$ | - | -
Drosophila | $10^{5}$ | $10^{7}$ | $3.5\cdot 10^{-5}$ | $0.08$
Mouse | $7\cdot 10^{7}$ | $10^{12}$ | $0.36$ | $440$
Rattus norvegicus | $2\cdot 10^{8}$ | $5\cdot 10^{11}$ | $1.8$ | $1200$
Rhesus | $6\cdot 10^{9}$ | - | $75$ | $94968$
Human | $10^{11}$ | $10^{14}$ | $1300$ | $1.2\cdot 10^{6}$
Table 1: Neurophysiology of the empirical connectomes considered in this work.
Information about the number of neurons, number of synapses, brain mass (in
grams) and brain volume (in cm3) for C. elegans, C. intestinalis, Drosophila,
Human, Mouse, Platynereis dumerilly, Rattus norvegicus and Rhesus macaque.
Figure 2: Prediction errors depend on system’s size and the nature of input
signals. Correlations between normalized MSE and number of nodes $N$. a)
white, b) fractional Brownian, c) Ikeda and d) Mackey Glass. Each point
represents the average across all networks of the same organism with the same
number of nodes. The color darkness reflects the total number of neurons in
the organism’s brain (Tab. 1). Correlations are significant ($p-value\leq
0.05$) for all signals. Positive values are measured for white and fractional
Brownian signals, while negative values are measured for Ikeda and Mackey
Glass signals.

For each combination of connectome and input signal, the corresponding
reservoir is trained and the activation states are stored in matrix $X$. Then
we compute the rank of the covariance matrix of $X$, as it measures the number
of independent (uncorrelated) vector states and can be used to evaluate the
dimensionality of the reservoir as a dynamical system Carroll and Pecora
[2019]. Consequently, we evaluate how RC performances – here measures as
normalized MSE – vary as a function of different network measures, namely
number of nodes $N$, average degree $\langle k\rangle$, rank $r$ of the
covariance matrix of activation states and network modularity $Q$.
Correlations between performances and networks’ macroscopic properties are
quantified by the Spearman coefficient, whose p-value is estimated using a
permutation test, that is by resampling without replacements $N$ and the
normalized MSE, and then evaluating the Spearman correlation for the resampled
data. This procedure is independently repeated $n=10,000$ times.

Figure 2 shows RC performance as a function of the number of nodes $N$ of the
connectomes. On the one hand, results show that the RC is actually able to
learn the chaotic dynamics of Ikeda and Mackey-Glass deterministic signals and
to reproduce them over time with a good approximation, with errors that range
between $10^{-2}$ and $10^{-7}$ (arbitrary units). On the other hand, the
performance of RC is lower against white and fractional Brownian noisy
signals, where prediction errors range between $10^{-1}$ and $10^{2}$ as
expected due to the highly stochastic nature of those sequences.

Correlations are found to be all significant; positive values are observed for
white and fractional Brownian signals while negative values are observed for
Ikeda and Mackey-Glass signals. Thus, for stochastic signals, the error
increases as the number of nodes increases, while for deterministic chaotic
signals the error decreases as the number of nodes increases. No correlations
are found between RC performance and the total number of neurons in the
organism brain from which the connectome is obtained (represented in different
colors in Fig. 2, see also Table 1). Similar correlation values and
statistical significances are observed for the average degree (Supplementary
Fig. S1), the rank of the covariance matrix of activation states
(Supplementary Fig. S2) and the network modularity (Supplementary Fig. S3).

Then we test the hypothesis that the observed results are consistent with the
topology of a set of synthetic data generated with three different generative
models preserving one or more features of the original connectomes. With this
approach, we probe which topological features characterizing empirical
connectomes are mostly responsible for the measured performances of the
machine. In this work, we consider three different classes of null models: an
Erdős Rényi (ER) random network, that only preserves the number of nodes and
the number of edges of the real connectome; a Configuration Model (CM), that
reproduces the degree distribution but destroys degree-degree correlations and
other structures observed in the real networks; and a Stochastic Block Model
(SBM), where the block structure is the topological feature that is preserved
by the synthetic network.

For each connectome, 20 independent realizations of each generative model are
considered, and for each realization the RC algorithm is run for one random
initialization of the weights and the average performance over all the
realizations is computed. To compare empirical and synthetic measures we use a
Z-score: negative scores indicate that the measure computed with the empirical
network is significantly greater than the one of the corresponding synthetic
network; conversely, positive scores indicate the opposite behavior. The two
measures we compare in empirical and synthetic networks include the rank of
the covariance matrix of activation states and the normalized MSE.

The number of independent activation states of the network, here measured as
the rank of their covariance matrix, characterizes the dimension of the
network actually involved in the requested task, as it discounts all the nodes
that follow identical trajectories. Higher ranks have been found to correlate
with lower error prediction Carroll and Pecora [2019], Lukoševičius and Jaeger
[2009]. Figure 3 shows the heatmaps representing the differences in the rank
of the covariance matrix between empirical networks and the corresponding
synthetic networks. Positive Z-scores prevail regardless of the nature of the
input signal. At fixed connectome, higher values are found as the complexity
of the null model increases, that is, the highest rank is found in SBM model
while the lowest is found in ER model. This result suggests that, among the
network measures that characterize an empirical connectome, the block
structure is the one that mainly could affect prediction performance. Along
the same synthetic model, higher values of the Z-score are prevalent for human
connectomes and, in general, the number of independent nodes increases with
the number of neurons of the considered organism (see Table 1). Our results
attest that synthetic networks, and in particular SBM, have a higher number of
independent vector states involved in the input prediction task and this would
suggest lower prediction errors for all the considered signals.

To test this hypothesis, and to verify if the trends of prediction
performances we found in empirical connectomes actually depend on topological
features, we compare errors in RC performances between empirical and synthetic
connectomes. Results are shown in Fig. 4 and indicate that, when the input
signal is deterministic, as for Ikeda and Mackey Glass signals, as expected
the normalized MSE is greater for real connectome reservoirs with respect to
the case in which the reservoir is a synthetic model. For white and fractional
Brownian input signals, we find that the performances of the empirical
connectomes are significantly better than the ones of the synthetic null
models. These results could be explained by overfitting occurring, for
stochastic signals, when the number of independent network nodes increases.

Figure 3: Number of independent activation states increases with null model
complexity and organism brain volume. Heatmap of Z-scores between the rank of
the covariance matrix of activation states obtained in empirical networks and
in Erdős Rényi (ER), Configuration Model (CM) and Stochastic Block Model (SBM)
networks. Positive Z-scores (blue) indicate that the number of independent
nodes of the empirical network is significantly lower than the one of the
corresponding synthetic network; negative Z-scores (brown) indicate that the
number of independent vector states is significantly greater in real
connectomes. Light gray indicates non significant results (Z-score in
$[-2,2]$). Positive Z-scores prevail for all signals and, across different
synthetic networks, higher absolute values of the Z-score are found as the
complexity of the null model increases. High absolute values of the Z-score
are prevalent for Human connectomes. Figure 4: Topological features and input
dynamics affect prediction performance. Heatmap of Z-scores between normalized
MSE obtained in empirical networks and normalized MSE in Erdős Rényi (ER),
Configuration Model (CM) and Stochastic Block Model (SBM) networks. Positive
Z-score (blue) indicate that the performances of the empirical network are
significantly better than the ones of the corresponding synthetic network;
negative Z-scores (brown) indicate that the performances of the synthetic
network are significantly better than the ones of the empirical network. Light
gray indicates non significant results (Z-score in $[-2,2]$). Positive
Z-scores prevail for white and fractional Brownian signals. Negative values,
instead, prevail for Ikeda and Mackey-Glass signals. High absolute values of
the Z-score are prevalent for Human connectomes. Also, higher absolute values
of the Z-score are found as the complexity of the null model increases.

##  I Discussion

In the present work, we exploit the framework of RC to investigate the
relationship between structure and function in real brain connectomes. In
practice, we use that state-of-the-art machine learning technique to gain
insights about empirical neural systems in terms of how the underlying network
structure can influence the performances of a prediction task.

Previous studies used biologically-inspired structures to perform
computational tasks and showed that the network topology could relevantly
affect neural computation and achieved valuable performance in a requested
task Bray [1995], Damicelli et al. [2021], Suárez et al. [2021]. Our work
capitalizes on those findings and clarifies which topological features,
typical of real networks, are relevant for computation. To this purpose, we
used a wide data set of 45 connectomes extracted from 8 different species,
validating our results by means of different null synthetic models with
increasing topological complexity. We found that multiple features correlate
with RC performance and that, in particular, a modular structure has a
significant impact on prediction error. This result is in agreement with
studies that suggested how functional Meunier et al. [2009] and structural
Chen et al. [2008], Bassett et al. [2010] hierarchical modularity Meunier et
al. [2009], Werner [2010] is characteristic of human brain connectivity. Block
structure produces a compartmentalization that reduces the interdependence of
modules and enhances robustness Kirschner and Gerhart [1998], Félix and Wagner
[2008] and allows the formation of complex hierarchical architectures that
promote a high functional specificity of subsystems Bassett and Gazzaniga
[2011]. Furthermore, brain modular structure allows different regions to act
either as hubs of high connectivity or as peripheral nodes of local processing
Bassett and Gazzaniga [2011], Rodriguez et al. [2019]. In addition to
robustness enhancement and specificity, modularity facilitates behavioral
adaptation Kirschner and Gerhart [1998] because each subsystem can modify its
function without altering other subsystems functioning. The fact that real
connectomes achieve a lower performance than the corresponding modular
synthetic model can be explained by considering the wiring cost that real
brain networks have to minimize: synapses functioning contributes to most of
the energy consumption in the human brain Lennie [2003] thus hierarchical
modularity of brain structure needs to be compatible with an evolutionary
pressure to minimize energy consumption while building and maintaining brain
connectivity Bassett et al. [2010], Niven and Laughlin [2008]. Short-range
connections between neurons are predominant in real brain networks Horvát et
al. [2016], while wiring cost constraints are not affecting synthetic
generative models.

Another measure that we find to affect RC performances is the dimensionality
of the reservoir as a dynamical system, that is the number of independent
vector states. High reservoir dimensionality is fundamental to ensure the
input representation is rich enough for the readout training procedure to work
properly. It helps the separation of initially inseparable inputs when
performing classification tasks Tanaka et al. [2019b]. Moreover, it allows
learning spatio-temporal dependencies of inputs when performing prediction
tasks Tanaka et al. [2019b]. Usually, reservoirs contain a large number of
neurons with order ranging from tens to thousands Lukoševičius and Jaeger
[2009]. An interesting result of our work is that differences in ranks between
empirical and synthetic networks depend on the number of neurons of the
considered organisms, suggesting that wiring patterns change as a function of
the species’ complexity.

Moreover, we run our analysis for different input signals and we find that,
unlike what we could expect from the number of independent vector states, when
the signals that have to be predicted are noisy the performances for
reservoirs with modular structures are lower than the one evaluated with
empirical connectomes. Our finding could be computationally explained by an
overfitting effect and is in agreement with the results of previous works,
according to which biological networks have optimized topologies capable of
suppressing noisy fluctuations that consist in sparser and more hierarchically
organized structures Ronellenfitsch et al. [2018].

There are also some limitations of our study that have to be considered in the
interpretation of the results. First, we study the reservoir dynamics at a
fixed amount of memory and nonlinearity, while it has been shown that
performances could vary depending on the underlying dynamical regime Suárez et
al. [2021]. Second, we collected our connectome data set from several studies.
Even though this has the advantage to offer an exhaustive analysis across
several organisms and species, the different experimental methodologies that
have been used to reconstruct connectomes could affect our comparisons.
Finally, in this work, we focus on a specific task, namely the prediction of a
signal provided in input to the reservoir, that represents a narrow window on
the complexity of the cognitive functions performed by the brain. In order to
better understand how cognition and perception are implemented in the brain,
other processes should be investigated, including memory and classification as
well as dynamic regimes with critical transitions that are typical of the
brain.

## Material and Methods

### Connectome and synthetic networks data

For each empirical connectome analogous synthetic networks are generated as
null models of increasing complexity to verify if the behavior of the
empirical networks actually depends on topological correlations. Erdős Rényi
(ER) model, Configuration Model (CM) and Stochastic Block Model (SBM) are used
to generate random graphs with the same number of nodes of the corresponding
empirical network. With ER, CM and SBM it is possible to produce synthetic
networks that randomize the network topology while preserving respectively:

  * (ER) the average degree,

  * (CM) the degree distribution,

  * (SBM) the degree distribution and the block structure

of the empirical network.

For each connectome (45) and signal (8) the algorithm is run with 20 different
random initializations of the weights and the average performances over all
the realizations are computed. Results referring to networks of the same
organism with the same number of neurons are averaged as well if convenient.

Network label | Num. of nodes | $\mathbf{\langle k\rangle}$ | Ref.
---|---|---|---
C. elegans herm. | 279 | 16.4 | git [a]
C. elegans male | 269 | 21.6 | Cook et al. [2019]
C. intestinalis 1 | 222 | 20.3 | Ryan et al. [2016]
C. intestinalis 2 | 138 | 6 | Ryan et al. [2016]
Drosophila medulla | 1770 | 10.1 | Takemura et al. [2013]
Drosophila con. mat. | 365 | 32.6 | Neu
Drosophila adj. | 49 | 43 | Neu
Human DS01216 1 | 1121 | 32.4 | Neu
Human Talairach 1 | 959 | 44.7 | Neu
Human DS01216 2 | 1121 | 30.9 | Neu
Human Talairach 2 | 959 | 43.2 | Neu
Human DS01216 3 | 1121 | 32.2 | Neu
Human Talairach 3 | 959 | 45 | Neu
Human DS01216 4 | 1121 | 29.6 | Neu
Human Talairach 4 | 959 | 40.2 | Neu
Human DS01216 5 | 1121 | 30.1 | Neu
Human Talairach 5 | 959 | 42.5 | Neu
Human Rockland 1 | 188 | 13.2 | Hagmann et al. [2008]
Human Rockland 2 | 188 | 14.8 | Hagmann et al. [2008]
Human Rockland 3 | 188 | 14.4 | Hagmann et al. [2008]
Human Rockland 4 | 188 | 12.8 | Hagmann et al. [2008]
Human Rockland 5 | 188 | 16.7 | Hagmann et al. [2008]
Human Rockland 6 | 188 | 16.6 | Hagmann et al. [2008]
Human Rockland 7 | 188 | 16.0 | Hagmann et al. [2008]
Human Rockland 8 | 188 | 13.5 | Hagmann et al. [2008]
Human Rockland 9 | 188 | 15.5 | Hagmann et al. [2008]
Human Rockland 10 | 188 | 18.1 | Hagmann et al. [2008]
Human Hagmann 1 | 998 | 29.2 | Hagmann et al. [2008]
Human Hagmann 2 | 998 | 29.5 | Hagmann et al. [2008]
Human Hagmann 3 | 998 | 29.4 | Hagmann et al. [2008]
Human Hagmann 4 | 998 | 27.3 | Hagmann et al. [2008]
Human Hagmann 5 | 998 | 28.6 | Hagmann et al. [2008]
Mouse kasthuri | 987 | 3.1 | Neu
Mouse brain | 213 | 151.1 | Neu
Mouse retina | 1076 | 168.8 | git [b]
Mouse visual cortex 1 | 29 | 3 | Bock et al. [2011]
Mouse visual cortex 2 | 193 | 2.2 | Bock et al. [2011]
Platynereis dumerilli | 2728 | 7.9 | Verasztó et al. [2020]
Rattus norvegicus 1 | 503 | 91.6 | Neu
Rattus norvegicus 2 | 503 | 98.2 | Neu
Rattus norvegicus 3 | 503 | 105.4 | Neu
Rhesus brain 1 | 242 | 25.2 | Harriger et al. [2012]
Rhesus brain 2 | 91 | 12.8 | Neu
Rhesus cerebral cortex | 91 | 30.8 | Neu
Rhesus interareal cort. net. | 93 | 48.6 | Neu
Table 2: List of connectomes of C. elegans, C. intestinalis, Drosophila,
Human, Mouse, Platynereis dumerilly, Rattus norvegicus and Rhesus macaque. For
each connectome the label of the network, the number of nodes of the network
and the average degree $\langle k\rangle$ of the network are shown.

### Connectome-based RC

RC is a Recurrent Neural Network passively excited by the input signal and
that maintains, in its internal state, a non-linear transformation of the
input history. The output signal is then generated as a linear combination of
the neuron’s signals using a set of readout units that are trained to produce
the desired target output. The key point is that the high-dimensional
reservoir contains a representation of the input that is sufficiently rich to
obtain the target output with a cheap linear learning procedure. Let’s
consider a training data set of ordered input vectors
$\vec{u}(n)\in\mathbb{R}^{N_{u}}$, $n=1,\dots,T$ each with its corresponding
target vector $\vec{y}_{target}(n)\in\mathbb{R}^{N_{y}}$. The goal is to learn
the function

$\vec{y}(n)=y(\dots,\vec{u}(n-1),\vec{u}(n))$ (1)

that maps the input vector into the output vector minimizing an error function
$E(\vec{y},\vec{y}_{target})$.

Firstly, the input vector is nonlinearly mapped into the reservoir, which is
composed of $N$ neurons. The activation of each neuron determines the
reservoir state vector $\vec{x}(n)\in\mathbb{R}^{N}$. The input vector
expansion can be written as:

$\vec{x}(n)=x(\dots,\vec{u}(n-1),\vec{u}(n))=x(\vec{x}(n-1),\vec{u}(n))$ (2)

with the expansion function that has memory, namely it depends on previous
input vectors. In this work, the following form for the expansion function is
chosen

$\vec{x}(n)=(1-\alpha)\vec{x}(n-1)+\alpha\tanh(W_{in}\vec{u^{\prime}}(n)+W\vec{x}(n-1)+\vec{b})$
(3)

where $\alpha$ is the leakage rate,
$\vec{u^{\prime}}(n)=(1,u_{0}(n),\dots,u_{N_{u}}(n))^{T}$ is the $n$-th biased
input vector and $\vec{b}\in\mathbb{R}^{N}$ is a bias. Hyperbolic tangent
function is used as neuron activation function. Secondly, the state vector is
mapped into the output with a linear transformation

$\vec{y}(n)=W_{out}\vec{x}(n)$ (4)

where $W_{out}\in\mathbb{R}^{N_{y}\times N}$ is the output weights matrix.
Learning procedure consists in finding the optimal values of the output
weights matrix. The vector states are sequentially computed for all points in
the training data set and are stored in the matrix
$X=(\vec{x}(1),\dots,\vec{x}(T))$. Given the target vectors
$\vec{Y}_{target}=(\vec{y}_{target}(1),\dots,\vec{y}_{target}(T))$, the
optimal output weights matrix can be computed as:

$W_{out}=\vec{Y}_{target}X^{T}(XX^{T}+\beta^{2}\mathbb{I})^{-1}$ (5)

where $\beta$ is a regularization parameter. In this work, RC networks are
connectomes provided as edge lists, namely lists of two indexes specifying the
connected nodes. The weight of the connection is not considered. Connectome
data are pre-processed to obtain undirected unweighted connected graphs
without multiple edges and self-loops. From each edge list the adjacency
matrix $A_{conn}$ of the network is computed. Given a connectome with $N$
neurons described by the adjacency matrix $A_{conn}\in\mathbb{R}^{N\times N}$,
the weight matrix $W$ of the reservoir is computed as element-wise
multiplication

$W=A_{conn}W_{weight},$ (6)

where $W_{weight}$ is drawn from a uniform distribution in $[-1,1]$.
Eventually, $W$ is rescaled as follows:

$W=W\frac{\rho_{exp}}{\rho(W)}$ (7)

with $\rho_{exp}=1.25$ to set the value of the spectral radius to obtain the
desired amount of memory and nonlinearity required.

The input matrix $W_{in}$ is drawn from a uniform distribution in
$\big{[}-\frac{1}{2},\frac{1}{2}\big{]}$. Moreover, the leakage rate is set to
$\alpha=0.3$ and the regularization parameter is set to be $\beta=10^{-8}$.
The bias vector $\vec{b}$ elements are sampled from a uniform distribution in
$[-\frac{1}{2},\frac{1}{2}]$ and multiplied by a bias intensity factor of 0.1.

### RC training and testing

Each signal is a series $\vec{U}=(u(1),\dots,u(T))$ of one dimensional vectors
$u(n)\in\mathbb{R}$. A transient period is of $L_{trans}=100$ points of the
signal is considered

$\vec{U}_{trans}=(u(1),\dots,u(L_{trans}));$ (8)

the reservoir dynamics runs, as in Eq. 3, to initialize the activation pattern
and no state vector is stored in this phase.

Then, the system runs with $L_{train}=2500$ points of the signal

$\vec{U}_{train}=(u(L_{trans}+1),\dots,u(L_{trans}+L_{train}))$ (9)

and activation states are stored in matrix $X$.

The columns of $X$ may be correlated with each other so, in order to evaluate
the number of uncorrelated vector states, we compute the rank of the
covariance matrix of $X$ as

$\Gamma=\text{rank}(X^{T}X)$ (10)

The output weight matrix $W_{out}$ can now be trained for the one-step-ahead
prediction task as in Eq. 4, using as target vector the input vector shifted
forward of one unit

$\vec{Y}_{target}=(u(L_{trans}+2),\dots,u(L_{trans}+L_{train}+1)).$ (11)

The trained RC is assessed by analyzing the performances on a test set of
$L_{test}=5000$ points of the signal

$\vec{U}_{test}=(u(L_{trans}+L_{train}+1),\dots,u(L_{trans}+L_{train}+L_{test})).$
(12)

The predictions are stored in the vector $\vec{Y}_{pred}$ and the performances
are estimated based on the Mean Square Error (MSE) computed over a number
$\tau$ of points as in the following equation:

$MSE(\tau)=\frac{\sum_{i=1}^{\tau}(\vec{Y}_{pred,i}-\vec{Y}_{target,i})^{2}}{\tau};$
(13)

and in particular the normalized MSE:

$MSE^{*}=\frac{MSE}{\text{std}(\vec{Y}_{target})}.$ (14)

For each connectome and input signal, the algorithm is run with 20 different
random initializations of the weights and the average performances over all
the realizations are computed.

Differences in performances between empirical connectomes and synthetic null
models are evaluated using the Z-score as in the following equation:

$Z=\frac{\overline{MSE}_{syn}-\overline{MSE}_{emp}}{\sqrt{\sigma_{MSE_{syn}}^{2}+\sigma_{MSE_{emp}}^{2}}}$
(15)

where $\overline{MSE}$ is the the mean value of the normalized MSE over
different realizations and $\sigma_{MSE}$ is the standard deviation.

### Input signals

Different time series (signals) are used as inputs. Each signal is a series
$\vec{U}=(u(1),\dots,u(T))$ of one dimensional vectors $u(n)\in\mathbb{R}$.
Four signals are considered: two stochastic signals (white noise, fractional
Brownian noise) and 2 deterministic chaotic signals (Ikeda and Mackey-Glass).

The fractional Brownian motion is a generalization of the Brownian motion that
consists in a continuous-time Gaussian process where the increments are not
independent Mandelbrot and Van Ness [1968]. The Brownian motion describes the
random motion of particles suspended in a medium at thermal equilibrium,
namely with no preferential direction of flow. The motion is described by the
following equation:

$\begin{split}&B_{H}(t)=B_{H}(0)\\\
+&\frac{1}{\Gamma(H+1/2)}\bigg{\\{}\int_{-\infty}^{0}\left[(t-s)^{H-1/2}-(-s)^{H-1/2}\right]dB(s)\\\
+&\int_{0}^{t}(t-s)^{H-1/2}dB(s)\bigg{\\}}\end{split}$ (16)

where the integration is with respect to the white noise measure $dB(s)$. $H$
is a parameter satisfying $0<H<1$; note that with $H=1/2$ the Brownian motion
is recovered. The increment process can be written as

$X(t)=B_{H}(t+1)-B_{H}(t).$ (17)

Deterministic chaotic signals are generated by dynamical systems which appear
to show a completely chaotic random behavior. Though, their behavior is
deterministic because it can be fully determined by time evolution equations
and initial conditions alone Lynch [2018]. Chaotic systems are usually
characterized by a long-term aperiodic behavior and a strong sensitivity to
initial conditions Lynch [2018]. The following systems are considered to
generate the signals.

  * •

Ikeda map is derived as a model of the light path across a non-linear optical
resonator Ikeda [1979]. The equations are:

$\begin{split}x_{n+1}&=A+B\left(x_{n}\cos\lvert
E_{n}\rvert^{2}-y_{n}\sin\lvert E_{n}\rvert^{2}\right)\\\
y_{n+1}&=B\left(x_{n}\sin\lvert E_{n}\rvert^{2}+y_{n}\cos\lvert
E_{n}\rvert^{2}\right)\\\ E_{n}&=x_{n}+iy_{n}\end{split}$ (18)

where $A$, $B$ are constants.

  * •

Mackey-Glass equations are delay differential equations designed to mimic both
healthy and pathological behavior in some biological systems such as the
variation in the relative quantity of mature cells in the blood Mackey and
Glass [1977]. The equation is:

$\frac{dx}{dt}=\frac{\beta x(t-\tau)}{1+x(t-\tau)^{n}-\delta x(t)}$ (19)

with $\tau$ and $n$ constants and $x$ represents the density of cells over
time.

white random noise is added to the signals to provide the RC input vectors
with a predefined Signal to Noise Ratio (SNR) of 100. Moreover, each signal is
normalized in the interval $\big{[}-\frac{1}{2},\frac{1}{2}\big{]}$

## References

  * [1] Neurodata.  https://neurodata.io/project/connectomes/.
  * git [a] Anterior herm. c. elegans connectome.  https://github.com/ericmjonas/circuitdata/tree/master/celegans_herm, a.
  * git [b] Mouse retina connectome.  https://github.com/ericmjonas/circuitdata/tree/master/mouseretina, b.
  * Achard and Bullmore [2007] S. Achard and E. Bullmore.  Efficiency and cost of economical brain functional networks.  _PLoS computational biology_ , 3:e17, 03 2007.  doi: 10.1371/journal.pcbi.0030017.
  * Bassett and Gazzaniga [2011] D. Bassett and M. Gazzaniga.  Understanding complexity in the human brain.  _Trends in cognitive sciences_ , 15:200–9, 05 2011.  doi: 10.1016/j.tics.2011.03.006.
  * Bassett et al. [2007] D. Bassett, A. Meyer-Lindenberg, S. Achard, T. Duke, and E. Bullmore.  Adaptive reconfiguration of fractal small-world human brain functional networks.  _Proceedings of the National Academy of Sciences of the United States of America_ , 103:19518–23, 01 2007.  doi: 10.1073/pnas.0606005103.
  * Bassett et al. [2009] D. Bassett, E. Bullmore, A. Meyer-Lindenberg, J. Apud, D. Weinberger, and R. Coppola.  Cognitive fitness of cost-efficient brain functional networks.  _Proceedings of the National Academy of Sciences of the United States of America_ , 106:11747–52, 07 2009.  doi: 10.1073/pnas.0903641106.
  * Bassett et al. [2010] D. Bassett, D. Greenfield, A. Meyer-Lindenberg, D. Weinberger, S. Moore, and E. Bullmore.  Efficient physical embedding of topologically complex information processing networks in brains and computer circuits.  _PLoS computational biology_ , 6:e1000748, 04 2010.  doi: 10.1371/journal.pcbi.1000748.
  * Bock et al. [2011] D. D. Bock, W.-C. A. Lee, A. M. Kerlin, M. L. Andermann, G. Hood, A. W. Wetzel, S. Yurgenson, E. R. Soucy, H. S. Kim, and R. C. Reid.  Network anatomy and in vivo physiology of visual cortical neurons.  _Nature_ , 471(7337):177–182, Mar 2011.  ISSN 1476-4687.  doi: 10.1038/nature09802.  URL https://doi.org/10.1038/nature09802.
  * Bray [1995] D. Bray.  Protein molecules as computational elements in living cells.  _Nature_ , 376(6538):307–312, 1995.  ISSN 1476-4687.  doi: 10.1038/376307a0.  URL https://doi.org/10.1038/376307a0.
  * Bullmore and Sporns [2012] E. Bullmore and O. Sporns.  The economy of brain network organization.  _Nature reviews neuroscience_ , 13(5):336–349, 2012.
  * Buonomano and Merzenich [1995] D. V. Buonomano and M. Merzenich.  Temporal information transformed into a spatial code by a neural network with realistic properties.  _Science_ , 267:1028 – 1030, 1995.
  * Burioni et al. [2014] R. Burioni, M. Casartelli, M. di Volo, R. Livi, and A. Vezzani.  Average synaptic activity and neural networks topology: a global inverse problem.  _Sci. Rep._ , 4(1):4336, Mar. 2014.
  * Carroll and Pecora [2019] T. L. Carroll and L. M. Pecora.  Network structure effects in reservoir computers.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 29(8):083130, 2019.  doi: 10.1063/1.5097686.  URL https://doi.org/10.1063/1.5097686.
  * Catani et al. [2005] M. Catani, D. Jones, and D. ffytche.  Perisylvian language networks of the human brain.  _Annals of neurology_ , 57:8–16, 01 2005.  doi: 10.1002/ana.20319.
  * Chen et al. [2008] Z. Chen, Y. He, P. Rosa, J. Germann, and A. Evans.  Revealing modular architecture of human brain structural networks by using cortical thickness from mri.  _Cerebral cortex (New York, N.Y. : 1991)_ , 18:2374–81, 03 2008.  doi: 10.1093/cercor/bhn003.
  * Churchland et al. [2010] M. Churchland, J. Cunningham, M. Kaufman, S. Ryu, and K. Shenoy.  Cortical preparatory activity: Representation of movement or first cog in a dynamical machine?  _Neuron_ , 68:387–400, 11 2010.  doi: 10.1016/j.neuron.2010.09.015.
  * Cook et al. [2019] S. J. Cook, T. A. Jarrell, C. A. Brittin, Y. Wang, A. E. Bloniarz, M. A. Yakovlev, K. C. Q. Nguyen, L. T.-H. Tang, E. A. Bayer, J. S. Duerr, H. E. Bülow, O. Hobert, D. H. Hall, and S. W. Emmons.  Whole-animal connectomes of both caenorhabditis elegans sexes.  _Nature_ , 571(7763):63–71, Jul 2019.  ISSN 1476-4687.  doi: 10.1038/s41586-019-1352-7.  URL https://doi.org/10.1038/s41586-019-1352-7.
  * Damicelli et al. [2021] F. Damicelli, C. C. Hilgetag, and A. Goulas.  Brain Connectivity meets Reservoir Computing, 2021.  URL https://www.biorxiv.org/content/early/2021/01/23/2021.01.22.427750.
  * De Domenico et al. [2016] M. De Domenico, S. Sasai, and A. Arenas.  Mapping multiplex hubs in human functional brain networks.  _Frontiers in neuroscience_ , 10:326, 2016.
  * De Vico Fallani et al. [2008] F. De Vico Fallani, L. Astolfi, F. Cincotti, D. Mattia, M. Marciani, S. Gao, S. Salinari, R. Soranzo, and F. Babiloni.  Structure of the cortical networks during successful memory encoding in tv commercials.  _Clinical neurophysiology : official journal of the International Federation of Clinical Neurophysiology_ , 119:2231–7, 09 2008.  doi: 10.1016/j.clinph.2008.06.018.
  * de Vico Fallani et al. [2014] F. de Vico Fallani, J. Richiardi, M. Chavez, and S. Achard.  Graph analysis of functional brain networks: practical issues in translational neuroscience.  _Philosophical Transactions of the Royal Society B: Biological Sciences_ , 369(1653):20130521, 2014.
  * Deco et al. [2011] G. Deco, V. Jirsa, and A. McIntosh.  Emerging concepts for the dynamical organization of resting-state activity in the brain.  _Nature reviews. Neuroscience_ , 12:43–56, 01 2011.  doi: 10.1038/nrn2961.
  * di Santo et al. [2018] S. di Santo, P. Villegas, R. Burioni, and M. A. Muñoz.  Landau-Ginzburg theory of cortex dynamics: Scale-free avalanches emerge at the edge of synchronization.  _Proc. Natl. Acad. Sci. U. S. A._ , 115(7):E1356–E1365, Feb. 2018.
  * Filan et al. [2020] D. Filan, S. Hod, C. Wild, A. Critch, and S. Russell.  Neural networks are surprisingly modular.  _CoRR_ , abs/2003.04881, 2020.  URL https://arxiv.org/abs/2003.04881.
  * Félix and Wagner [2008] M.-A. Félix and A. Wagner.  Robustness and evolution: Concepts, insights and challenges from a developmental model system.  _Heredity_ , 100:132–40, 03 2008.  doi: 10.1038/sj.hdy.6800915.
  * Gauthier et al. [2021] D. J. Gauthier, E. Bollt, A. Griffith, and W. A. Barbosa.  Next generation reservoir computing.  _Nature communications_ , 12(1):1–8, 2021.
  * Gong et al. [2009] G. Gong, P. Rosa-Neto, F. Carbonell, Z. J. Chen, Y. He, and A. C. Evans.  Age- and gender-related differences in the cortical anatomical network.  _Journal of Neuroscience_ , 29(50):15684–15693, 2009.  ISSN 0270-6474.  doi: 10.1523/JNEUROSCI.2308-09.2009.  URL https://www.jneurosci.org/content/29/50/15684.
  * Guillon et al. [2017] J. Guillon, Y. Attal, O. Colliot, V. La Corte, B. Dubois, D. Schwartz, M. Chavez, and F. de Vico Fallani.  Loss of brain inter-frequency hubs in alzheimer’s disease.  _Scientific reports_ , 7(1):1–13, 2017.
  * Hagmann et al. [2008] P. Hagmann, L. Cammoun, X. Gigandet, R. Meuli, C. Honey, V. Wedeen, and O. Sporns.  Mapping the structural core of human cerebral cortex.  _PLoS biology_ , 6:e159, 08 2008.  doi: 10.1371/journal.pbio.0060159.
  * Harriger et al. [2012] L. Harriger, M. Heuvel, and O. Sporns.  Rich club organization of macaque cerebral cortex and its role in network communication.  _PloS one_ , 7:e46497, 09 2012.  doi: 10.1371/journal.pone.0046497.
  * Horstmann et al. [2010] M.-T. Horstmann, S. Bialonski, N. Noennig, H. Mai, J. Prusseit, J. Wellmer, H. Hinrichs, and K. Lehnertz.  State dependent properties of epileptic brain networks: Comparative graph-theoretical analyses of simultaneously recorded eeg and meg.  _Clinical neurophysiology : official journal of the International Federation of Clinical Neurophysiology_ , 121:172–85, 02 2010.  doi: 10.1016/j.clinph.2009.10.013.
  * Horvát et al. [2016] S. Horvát, R. Gămănuț, M. Ercsey-Ravasz, L. Magrou, B. Gămănuț, D. C. Van Essen, A. Burkhalter, K. Knoblauch, Z. Toroczkai, and H. Kennedy.  Spatial Embedding and Wiring Cost Constrain the Functional Layout of the Cortical Network of Rodents and Primates.  _PLOS Biol._ , 14(7):1–30, 2016.  doi: 10.1371/journal.pbio.1002512.  URL https://doi.org/10.1371/journal.pbio.1002512.
  * Häusler and Maass [2007] S. Häusler and W. Maass.  A statistical analysis of information-processing properties of lamina-specific cortical microcircuit models.  _Cerebral cortex (New York, N.Y. : 1991)_ , 17:149–62, 02 2007.  doi: 10.1093/cercor/bhj132.
  * Ikeda [1979] K. Ikeda.  Multiple-valued stationary state and its instability of the transmitted light by a ring cavity system.  _Optics Communications_ , 30(2):257–261, 1979\.  ISSN 0030-4018.  doi: https://doi.org/10.1016/0030-4018(79)90090-7.  URL https://www.sciencedirect.com/science/article/pii/0030401879900907.
  * Karmarkar and Buonomano [2007] U. Karmarkar and D. Buonomano.  Timing in the absence of clocks: Encoding time in neural network states.  _Neuron_ , 53:427–38, 03 2007.  doi: 10.1016/j.neuron.2007.01.006.
  * Kirschner and Gerhart [1998] M. Kirschner and J. Gerhart.  Evolvability.  _Proceedings of the National Academy of Sciences_ , 95(15):8420–8427, 1998.  ISSN 0027-8424.  doi: 10.1073/pnas.95.15.8420.  URL https://www.pnas.org/content/95/15/8420.
  * Kistler and De Zeeuw [2002] W. Kistler and C. De Zeeuw.  Dynamical working memory and timed responses: The role of reverberating loops in the olivo-cerebellar system.  _Neural computation_ , 14:2597–626, 12 2002.  doi: 10.1162/089976602760407991.
  * Lennie [2003] P. Lennie.  The cost of cortical computation.  _Current biology : CB_ , 13:493–7, 04 2003.  doi: 10.1016/S0960-9822(03)00135-0.
  * Li et al. [2009] J. Li, B. E. Carlson, and A. A. Lacis.  A study on the temporal and spatial variability of absorbing aerosols using total ozone mapping spectrometer and ozone monitoring instrument aerosol index data.  _J. Geophys. Res._ , 114:D09213, 2009.  doi: 10.1029/2008JD011278.
  * Lukoševičius and Jaeger [2009] M. Lukoševičius and H. Jaeger.  Reservoir computing approaches to recurrent neural network training.  _Computer Science Review_ , 3(3):127–149, 2009\.  ISSN 1574-0137.  doi: https://doi.org/10.1016/j.cosrev.2009.03.005.  URL https://www.sciencedirect.com/science/article/pii/S1574013709000173.
  * Lynall et al. [2010] M.-E. Lynall, D. Bassett, R. Kerwin, P. McKenna, M. Kitzbichler, U. Müller-Sedgwick, and E. Bullmore.  Functional connectivity and brain networks in schizophrenia.  _The Journal of neuroscience : the official journal of the Society for Neuroscience_ , 30:9477–87, 07 2010.  doi: 10.1523/JNEUROSCI.0333-10.2010.
  * Lynch [2018] S. Lynch.  _Dynamical Systems with Applications using Python (Information)_.  10 2018.  ISBN 978-3-319-78145-7.
  * Mackey and Glass [1977] M. C. Mackey and L. Glass.  Oscillation and chaos in physiological control systems.  _Science_ , 197(4300):287–289, 1977.  doi: 10.1126/science.267326.  URL https://www.science.org/doi/abs/10.1126/science.267326.
  * Mandelbrot and Van Ness [1968] B. B. Mandelbrot and J. W. Van Ness.  Fractional brownian motions, fractional noises and applications.  _SIAM Review_ , 10(4):422–437, 1968.  doi: 10.1137/1010093.  URL https://doi.org/10.1137/1010093.
  * Meunier et al. [2008] D. Meunier, S. Achard, A. Morcom, and E. Bullmore.  Age-related changes in modular organization of human brain functional networks.  _NeuroImage_ , 44:715–23, 12 2008.  doi: 10.1016/j.neuroimage.2008.09.062.
  * Meunier et al. [2009] D. Meunier, R. Lambiotte, A. Fornito, K. Ersche, and E. Bullmore.  Hierarchical modularity in human brain functional networks.  _Frontiers in neuroinformatics_ , 3:37, 10 2009.  doi: 10.3389/neuro.11.037.2009.
  * Munoz [2018] M. A. Munoz.  Colloquium: Criticality and dynamical scaling in living systems.  _Reviews of Modern Physics_ , 90(3):031001, 2018\.
  * Murugan [2018] P. Murugan.  Learning the sequential temporal information with recurrent neural networks.  _CoRR_ , abs/1807.02857, 2018.  URL http://arxiv.org/abs/1807.02857.
  * Nikolić et al. [2006] D. Nikolić, S. Häusler, W. Singer, and W. Maass.  Temporal dynamics of information content carried by neurons in the primary visual cortex.  pages 1041–1048, 01 2006.
  * Niven and Laughlin [2008] J. Niven and S. Laughlin.  Energy limitation as a selective force on the evolution of sensory systems.  _The Journal of experimental biology_ , 211:1792–804, 07 2008.  doi: 10.1242/jeb.017574.
  * Paquot et al. [2012] Y. Paquot, F. Duport, A. Smerieri, J. Dambre, B. Schrauwen, M. Haelterman, and S. Massar.  Optoelectronic reservoir computing.  _Scientific reports_ , 2(1):1–6, 2012.
  * Presigny and Fallani [2022] C. Presigny and F. D. V. Fallani.  Colloquium: Multiscale modeling of brain network organization.  _Reviews of Modern Physics_ , 94(3):031002, 2022\.
  * Rodriguez et al. [2019] N. Rodriguez, E. Izquierdo, and Y.-Y. Ahn.  Optimal modularity and memory capacity of neural reservoirs.  _Netw. Neurosci. (Cambridge, Mass.)_ , 3(2):551–566, 2019.  ISSN 2472-1751 (Electronic).  doi: 10.1162/netn˙a˙00082.
  * Ronellenfitsch et al. [2018] H. Ronellenfitsch, J. Dunkel, and M. Wilczek.  Optimal noise-canceling networks.  _Phys. Rev. Lett._ , 121:208301, Nov 2018.  doi: 10.1103/PhysRevLett.121.208301.  URL https://link.aps.org/doi/10.1103/PhysRevLett.121.208301.
  * Ryan et al. [2016] K. Ryan, Z. Lu, and I. Meinertzhagen.  The cns connectome of a tadpole larva of ciona intestinalis (l.) highlights sidedness in the brain of a chordate sibling.  _eLife_ , 5, 12 2016.  doi: 10.7554/eLife.16962.
  * Schwarz et al. [2009] E. Schwarz, M. Wdowiak, D. Almer-Jarz, and R. Breitenecker.  The effects of attitudes and perceived environment conditions on students’ entrepreneurial intent.  _Education + Training_ , 51:272–291, 05 2009.  doi: 10.1108/00400910910964566.
  * Sifis et al. [2009] M. Sifis, M. Vourkas, V. Tsirka, E. Karakonstantaki, K. Kanatsouli, and C. Stam.  The influence of ageing on complex brain networks: A graph theoretical analysis.  _Human brain mapping_ , 30:200–8, 01 2009.  doi: 10.1002/hbm.20492.
  * Stam et al. [2008] C. Stam, W. Haan, A. Daffertshofer, B. Jones, I. Manshanden, A.-M. van Cappellen van Walsum, T. Montez, J. Verbunt, J. De Munck, B. Dijk, H. Berendse, and P. Scheltens.  Graph theoretical analysis of magnetoencephalographic functional connectivity in alzheimer’s disease.  _Brain : a journal of neurology_ , 132:213–24, 11 2008\.  doi: 10.1093/brain/awn262.
  * Stanley et al. [1999] G. B. Stanley, F. F. Li, and Y. Dan.  Reconstruction of Natural Scenes from Ensemble Responses in the Lateral Geniculate Nucleus.  _J. Neurosci._ , 19(18):8036–8042, 1999.  ISSN 0270-6474.  doi: 10.1523/JNEUROSCI.19-18-08036.1999.  URL https://www.jneurosci.org/content/19/18/8036.
  * Storrs and Kriegeskorte [2019] K. R. Storrs and N. Kriegeskorte.  Deep learning for cognitive neuroscience, 2019.
  * Suárez et al. [2021] L. E. Suárez, B. A. Richards, G. Lajoie, and B. Misic.  Learning function from structure in neuromorphic networks.  _Nat. Mach. Intell._ , 3(9):771–786, 2021.  ISSN 2522-5839.  doi: 10.1038/s42256-021-00376-1.  URL https://doi.org/10.1038/s42256-021-00376-1.
  * Takemura et al. [2013] S.-y. Takemura, A. Bharioke, Z. Lu, A. Nern, S. Vitaladevuni, P. K. Rivlin, W. T. Katz, D. J. Olbris, S. M. Plaza, P. Winston, T. Zhao, J. A. Horne, R. D. Fetter, S. Takemura, K. Blazek, L.-A. Chang, O. Ogundeyi, M. A. Saunders, V. Shapiro, C. Sigmund, G. M. Rubin, L. K. Scheffer, I. A. Meinertzhagen, and D. B. Chklovskii.  A visual motion detection circuit suggested by drosophila connectomics.  _Nature_ , 500(7461):175–181, Aug 2013.  ISSN 1476-4687.  doi: 10.1038/nature12450.  URL https://doi.org/10.1038/nature12450.
  * Tanaka et al. [2019a] G. Tanaka, T. Yamane, J. B. Héroux, R. Nakane, N. Kanazawa, S. Takeda, H. Numata, D. Nakano, and A. Hirose.  Recent advances in physical reservoir computing: A review.  _Neural Networks_ , 115:100–123, 2019a.
  * Tanaka et al. [2019b] G. Tanaka, T. Yamane, J. B. Héroux, R. Nakane, N. Kanazawa, S. Takeda, H. Numata, D. Nakano, and A. Hirose.  Recent advances in physical reservoir computing: A review.  _Neural Networks_ , 115:100–123, 2019b.  ISSN 0893-6080.  doi: https://doi.org/10.1016/j.neunet.2019.03.005.  URL https://www.sciencedirect.com/science/article/pii/S0893608019300784.
  * van den Heuvel et al. [2009] M. P. van den Heuvel, C. J. Stam, R. S. Kahn, and H. E. Hulshoff Pol.  Efficiency of functional brain networks and intellectual performance.  _Journal of Neuroscience_ , 29(23):7619–7624, 2009.  ISSN 0270-6474.  doi: 10.1523/JNEUROSCI.1443-09.2009.  URL https://www.jneurosci.org/content/29/23/7619.
  * Van der Sande et al. [2017] G. Van der Sande, D. Brunner, and M. C. Soriano.  Advances in photonic reservoir computing.  _Nanophotonics_ , 6(3):561–576, 2017.
  * Verasztó et al. [2020] C. Verasztó, S. Jasek, M. Gühmann, R. Shahidi, N. Ueda, J. D. Beard, S. Mendes, K. Heinz, L. A. Bezares-Calderón, E. Williams, and G. Jékely.  Whole-animal connectome and cell-type complement of the three-segmented platynereis dumerilii larva, 2020.  URL https://www.biorxiv.org/content/early/2020/08/23/2020.08.21.260984.
  * Werner [2010] G. Werner.  Fractals in the nervous system: Conceptual implications for theoretical neuroscience.  _Frontiers in physiology_ , 1:15, 07 2010.  doi: 10.3389/fphys.2010.00015.
  * Xie et al. [2019] S. Xie, A. Kirillov, R. B. Girshick, and K. He.  Exploring randomly wired neural networks for image recognition.  _CoRR_ , abs/1904.01569, 2019.  URL http://arxiv.org/abs/1904.01569.
  * Yamazaki and Tanaka [2007] T. Yamazaki and S. Tanaka.  The cerebellum as a liquid state machine.  _Neural networks : the official journal of the International Neural Network Society_ , 20:290–7, 05 2007.  doi: 10.1016/j.neunet.2007.04.004.

Figure S1: Correlations between normalized MSE and average degree. a) White,
b) Fractional Brownian, c) Ikeda and d) Mackey Glass. The color darkness of
markers is proportional to the total number of neurons in the organisms brain.
Correlations are significant ($p-value\leq 0.05$) for all signals. Positive
values are observed for White and Fractional Brownian signals while negative
values are observed for Ikeda and Mackey Glass signals.  Figure S2:
Correlations between normalized MSE and rank of the covariance of matrix of
activation states.. Fa) White, b) Fractional Brownian, c) Ikeda and d) Mackey
Glass. Each dot represents the average across all networks of the same
organism with the same number of nodes.The color darkness of markers is
proportional to the total number of neurons in the organisms brain.
Correlations are found significant ($p-value\leq 0.05$) for all signals,
except for the fractional Brownian signal. Positive values are observed for
white and fractional Brownian signals while negative values are observed for
Ikeda and Mackey Glass signals.  Figure S3: Correlations between normalized
MSE and modularity. a) White, b) Fractional Brownian, c) Ikeda and d) Mackey
Glass. The color darkness of markers is proportional to the total number of
neurons in the organisms brain. Correlations are significant ($p-value\leq
0.05$) for all signals. Positive values are observed for white and fractional
Brownian signals while negative values are observed for Ikeda and Mackey Glass
signals.