# Experimental Comparison of Hardware-Amenable Spike Detection Algorithms for
iBMIs

Shoeb Shaikh1, Rosa So2, Camilo Libedinsky3 and Arindam Basu1 **This work was
supported through grant RG87/16 by MOE, Singapore1Authors are with Nanyang
Technological University, Singapore<EMAIL_ADDRESS>is with the
Institute for Infocomm Research, Singapore<EMAIL_ADDRESS>is
with SiNAPSE, National University of Singapore, Singapore<EMAIL_ADDRESS>
###### Abstract

This paper presents an experiment based comparison of absolute threshold (AT)
and non-linear energy operator (NEO) spike detection algorithms in Intra-
cortical Brain Machine Interfaces (iBMIs). Results show an average increase in
decoding performance of $\approx$ 5% in monkey A across 28 sessions recorded
over 6 days and $\approx$ 2% in monkey B across 35 sessions recorded over 8
days when using NEO over AT. To the best of our knowledge, this is the first
ever reported comparison of spike detection algorithms in an iBMI experimental
framework involving two monkeys. Based on the improvements observed in an
experimental setting backed by previously reported improvements in simulation
studies, we advocate switching from state of the art spike detection technique
- AT to NEO.

##  I INTRODUCTION

Intra-cortical Brain Machine Interfaces (iBMIs) have enabled patients
suffering from debilitating spinal cord afflictions to regain a sense of
control through demonstrations such as [1, 2] where the subjects are
successfully able to feed themselves and [3] where record communication speeds
have been achieved among other things [4, 5]. Spikes serve as an input to
these systems, where instantaneous firing rates are computed on every recorded
electrode and subsequently mapped to a behavioral co-variate to drive an
effector or stimulate a paralysed limb.

The process of spike detection entails identifying a spike amidst the recorded
background noise. Absolute Threshold (AT) method is the state of the art
technique currently employed for this purpose in both non-human primate (NHP)
[6, 4, 7] and human subject [3, 2] based studies. However, this technique is
sensitive to placement of an optimal threshold. If one sets it too high,
spikes are missed and if one sets it too low, false detections occur [8].
Secondly, hardware development in iBMIs is progressing towards including spike
detection in the implant itself to counter the problem of transmitting
broadband raw data rate [9, 10, 11, 12]. For such a design decision, it is not
feasible to implement AT on chip as it requires large memory and concomitantly
large area for computation of a median-based threshold [13].

Figure 1: Experimental setup for the iBMI experiment adapted from [4]

An alternate technique by the name of non-linear energy operator (NEO) [14]
has been shown to perform better in spike detection than AT in studies based
on semi-synthetic data such as [15, 16]. Apart from being better adept at
spike detection due to reasons which will be briefly touched upon in $Section$
$II$, NEO enjoys two specific advantages over AT. Firstly, it has been shown
to be generally less sensitive to the placement of threshold as compared to AT
[16]. Secondly, it is computationally friendly without any demanding memory
requirements and extreme low power implementations such as [17, 18] consuming
$40$ nW/channel, $50$ nW/channel respectively. Other spike detection
techniques such as stationary wavelet transform product [19], exponential
component - power component (EC-PC) [20] have been reported in literature, but
these were not considered in this analysis owing to their relatively greater
computational complexity [16, 21] than both AT and NEO.

Spike detection methods reported so far in literature are based on studies
performed on synthetic and semi-synthetic data with realistic assumptions [8].
The logical next step we believe is to apply these findings in a behavioral
task-based experimental framework to check for their validity. An analysis
along similar lines has been reported in [13]. However, this study lacks on
two fronts. Firstly, it is restricted to just one NHP instead of the standard
protocol of two NHPs. Secondly, the proposed method’s hardware implementation
consumes up to three orders of magnitude more power than the ones reported for
NEO [21, 17, 18] reflecting it’s computational complexity.

Thus, the object of this study is to apply the simulation based finding of
NEO’s preponderance over state of the art - AT in a simple experimental
setting for two NHPs with iBMI performance as a yardstick. The line of
reasoning is similar to [7], where spike sorting and thresholding were
compared in context of iBMI performance.

##  II Spike Detection Algorithms

In this section, we briefly present AT and NEO spike detection methods.

###  II-A Absolute threshold

This method involves applying a threshold to the digitally filtered raw data
$x[n]$ at time $n$ to identify a spike. Threshold is computed over $N$ samples
of filtered raw data as,

$\sigma_{N}=\frac{median|x[n]|}{0.6745}$ (1) $\displaystyle
Thr_{AT}=k\times\sigma$ (2)

where $\sigma_{N}$ is the standard deviation estimate over $N$ samples of
filtered raw data, $k$ is an integer typically between $2$ and $5$. Negative
threshold crossings of the filtered raw signal are considered as spikes.

###  II-B Non-linear Energy Operator

This method involves application of a pre-emphasis step before spike detection
[15]. The idea behind this approach is to accentuate high energy, high
frequency content amidst background noise thereby increasing spike to noise
separation. NEO threshold calculation involves the following steps,

$\displaystyle\psi[n]=x^{2}[n]-x[n-1]\times x[n+1]$ (3) $\displaystyle
Thr_{NEO}=l\times\frac{\sum_{n=1}^{N}\psi[n]}{N}$ (4)

where $\psi[n]$ is the NEO processed, pre-emphasised signal at time $n$.
$Thr_{NEO}$ is then computed as a multiple - $l$ times the average of
$\psi[n]$ over $N$ samples of filtered raw data. The value of $l$ is optimized
as per the experimental conditions [16, 21] and is not as sensitive as AT due
to better spike to noise separation. Other methods of automatic threshold
setting that are less sensitive to amplitude of the NEO output can also be
used [17]. Also, NEO can be applied directly to the analog signals after the
amplifier and before the ADC to save the ADC power dissipation [17].

(a)
(b)

Figure 2: Decoder results for (a) monkey A and (b) monkey B. Statistical
significance test is conducted only between best decoding performance of AT
and NEO, significant results are denoted by (*), p$<$0.05 (paired t-test).
Higher decoding accuracies are obtained using NEO in both cases. Threshold
multiples $k$ and $l$ for AT and NEO respectively are marked in the figure for
each bar.

##  III METHODS

Figure 3: NEO vs AT visualization for a segment of neural recording showing
improved detection using NEO.

###  III-A Neural Data Acquisition

Multiple microelectrode arrays were implanted in the hand/arm area of primary
motor cortex. There were a total of $96$ and $64$ implanted electrodes sensing
extracellular neural signals in monkeys A and B respectively. These signals
were then passed through front-end circuitry sampling processed raw data at
$13$ $k$Hz [22]. The digitized raw data so obtained was subjected to digital
bandpass filtering between $300$ and $3000$ Hz implemented on a PC to yield
$x[n]$ at time $n$.

###  III-B Experiment

The experiment protocol is described in detail in [4]. Briefly, a monkey was
trained to manipulate the position of a robotic platform through a hand-
controlled three-directional joystick. Four tasks - moving forward by $2$ m,
turning $90\degree$ right, turning $90\degree$ left, staying still for $5$
seconds constituted the experiment. Trial was considered successful if the
monkey reached the target within $15$ seconds for tasks involving movement and
stayed still for $5$ seconds during the stop task. A successful trial was
rewarded with a fruit reward by the trainer.

###  III-C Decoding Algorithm

We have used linear discriminant analysis (LDA) as the decoder of choice for
mainly two reasons. Firstly, it is popularly used in BMIs operated in a
discrete control fashion [4, 23, 24]. Secondly, LDA is a relatively simpler
algorithm compared to more complex algorithms such as support vector machine
and deep neural nets which have in some cases shown to outperform simple
linear algorithms [24]. This ensures that the gains in decoding accuracy are
primarily attributed to the spike detection method and not the superior
generalization capabilities of the classification algorithm [13].

##  IV RESULTS

The LDA based decoder was trained de novo on the first session’s data on each
day. Subsequent $3$-$4$ sessions were used as test sessions for offline
analysis. A total of $28$ sessions across $6$ days and $35$ sessions across
$8$ days were used in the analysis for monkeys A and B respectively.
Thresholds for both AT and NEO methods were set using the initial 30 seconds
of raw data from the first session. For a fairer comparison multipliers $k$
and $l$ described in equations (2), (4) were swept for both AT and NEO
respectively.

The averaged decoder results across all the test sessions are plotted in
figures 2(a) and (b) respectively. For monkey A, NEO with $l=15$ ($NEO_{15}$)
and AT with $k=3$ ($AT_{3}$) perform the best respectively, and $NEO_{15}$
outperforms $AT_{3}$ by $5.29$ %. Along similar lines for monkey B, NEO with
$l=13$ ($NEO_{13}$) and AT with $k=3$ ($AT_{3}$) perform the best
respectively, and $NEO_{15}$ outperforms $AT_{3}$ by $2.16$ %. Another point
to note is that in both monkey cases, standard deviations of decoding
accuracies across sessions for NEO are less than AT. For monkey A, the
standard deviation of decoding accuracies across sessions for $NEO_{13}$ is
$8.09$ %, whereas for $AT_{3}$ it is $10.72$ %. Similarly for monkey B,
standard deviation of decoding accuracies across sessions for $NEO_{13}$ is
$5.04$ %, whereas for $AT_{3}$ it is $6.88$ %. Thus, NEO performs better than
AT with greater average and lower standard deviation of decoding accuracies
across sessions for both monkeys.

##  V DISCUSSION

Initial iBMI studies resorted to spike sorting to map cortical neural firing
rates to behavioral variables [25, 26]. However, studies such as [7, 27] have
shown that perhaps it is unnecessary, with little or no loss in decoding
accuracies obtained by forsaking sorting for simple threshold crossings.
Furthermore, there have also been studies to provide a quantitative
understanding for using threshold crossing in lieu of sorted spikes [28].
Thus, there is pressing evidence to simply use threshold crossings to monitor
the population level response as a whole from intra-cortical microelectrode
recordings. We argue that, NEO by virtue of being better adept at spike
detection yields a more informative population response.

In figure 3, we have visualized filtered raw data and NEO processed signal
along with different thresholds for NEO, AT for a single channel as a way to
demonstrate the difference in spike detection between these two approaches.
Ground truth spikes are marked in the form of downward pointing arrows. One
can see that for AT with $k=4$, $k=5$, the threshold levels are relatively
higher and some spikes are missed being detected. One might be tempted to go
for a further lower threshold with $k=3$ to rein in the missed detections.
However, with $k=3$, some false detections can be seen to occur. Simple
eyeballing in the lower half of the plot for NEO with $l=13$ reveals that
spikes are comfortably detected with relatively lesser concern for false
detection. One can easily observe the increased spike to noise ratio advantage
offered by NEO over AT in figure. 3.

The color-coded firing rate plots in figures 4(a) and (b) also show evidence
of increased firing rates appearing at the electrodes for NEO with $l=13$ as
compared to AT with $k=4$. This goes to show that NEO helps in detecting
spikes from neurons which were perhaps missed by AT and helps in gleaning
increased population-level information.

(a)

(b)

Figure 4: Mean firing rates on electrodes (plotted in an $8\times 8$ array
fashion) with (a) AT ($k=4$) and (b) NEO ($l=13$) for Monkey B when it is
performing the forward trials in a given session

##  VI CONCLUSION

To conclude, we have shown NEO to be better than AT in context of iBMI
performance for two monkeys with $\approx 2\%$ and $\approx 5\%$ improvement
in decoding accuracy respectively. The underlying hypothesis is that NEO helps
in better spike detection which in turn leads to enhanced performance. Thus,
based on previous simulation studies [15, 16], current results and hardware
amenability [17] we advocate switching from AT to NEO as the de facto spike
detection method. Future work includes closed-loop trials and applying NEO in
conjunction with more sophisticated decoders such as neural networks and
Kalman filters.

## ACKNOWLEDGMENT

The authors would like to thank Clement Lim for helping in animal training and
data collection.

## References

  * [1] Collinger, Wodlinger, _et al._ , “High-performance neuroprosthetic control by an individual with tetraplegia,” _The Lancet_ , no. 9866, pp. 557–564, feb.
  * [2] Ajiboye, Willett, _et al._ , “Restoration of reaching and grasping movements through brain-controlled muscle stimulation in a person with tetraplegia: a proof-of-concept demonstration,” _The Lancet_ , no. 10081, pp. 1821–1830.
  * [3] Pandarinath, Nuyujukian, _et al._ , “High performance communication by people with paralysis using an intracortical brain-computer interface,” _eLife_ , p. e18554, feb.
  * [4] Libedinsky, So, _et al._ , “Independent mobility achieved through a wireless brain-machine interface,” _PLoS ONE_ , vol. 11, no. 11, pp. 1–13, 2016.
  * [5] Rajangam, Tseng, _et al._ , “Wireless Cortical Brain-Machine Interface for Whole-Body Navigation in Primates,” _Scientific Reports_ , p. 22170\.
  * [6] Gilja, Nuyujukian, _et al._ , “A high-performance neural prosthesis enabled by control algorithm design,” _Nature Neuroscience_ , vol. 15, no. 12, pp. 1752–1757, 2012.
  * [7] Christie, Tat, _et al._ , “Comparison of spike sorting and thresholding of voltage waveforms for intracortical brain-machine interface performance,” _Journal of Neural Engineering_ , no. 1, p. 16009.
  * [8] Rey, Pedreira, _et al._ , “Past, present and future of spike sorting techniques,” _Brain Research Bulletin_ , pp. 106–117.
  * [9] Chattopadhyay, Chang, _et al._ , _Emerging technology and architecture for big-data analytics_ , 2017.
  * [10] Chen, Yao, _et al._ , “A 128-Channel Extreme Learning Machine-Based Neural Decoder for Brain Machine Interfaces,” _IEEE Transactions on Biomedical Circuits and Systems_ , vol. 10, no. 3, pp. 679–692, 2016.
  * [11] Boi, Moraitis, _et al._ , “A Bidirectional Brain-Machine Interface Featuring a Neuromorphic Hardware Decoder,” _Frontiers in Neuroscience_ , no. December, pp. 1–15.
  * [12] Harrison, Watkins, _et al._ , “A low-power integrated circuit for a wireless 100-electrode neural recording system,” _IEEE Journal of Solid-State Circuits_ , vol. 42, no. 1, pp. 123–133, 2007.
  * [13] Wing-kin Tam, So, _et al._ , “EC-PC spike detection for high performance brain-computer interface,” _2015 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)_ , no. 5, pp. 5142–5145.
  * [14] Kaiser, “On a simple algorithm to calculate the ’energy’ of a signal,” in _International Conference on Acoustics, Speech, and Signal Processing_ , April 1990, pp. 381–384 vol.1.
  * [15] Mukhopadhyay and Ray, “A new interpretation of nonlinear energy operator and its efficacy in spike detection,” _IEEE Transactions on Biomedical Engineering_ , vol. 45, no. 2, pp. 180–187, 1998.
  * [16] Gibson, Judy, _et al._ , “Technology-aware algorithm design for neural spike detection, feature extraction, and dimensionality reduction.” _IEEE transactions on neural systems and rehabilitation engineering : a publication of the IEEE Engineering in Medicine and Biology Society_ , vol. 18, no. 5, pp. 469–478, 2010.
  * [17] Yao, Chen, _et al._ , “A 0.7 V, 40 nW Compact, Current-Mode Neural Spike Detector in 65 nm CMOS,” vol. 10, no. 2, pp. 309–318, 2016.
  * [18] Yang and Mason, “Hardware Efficient Automatic Thresholding for NEO-Based Neural Spike Detection,” _IEEE Transactions on Biomedical Engineering_ , vol. 64, no. 4, pp. 826–833, 2017.
  * [19] Kim and Kim, “A wavelet-based method for action potential detection from extracellular neural signal recording with low signal-to-noise ratio,” _IEEE Transactions on Biomedical Engineering_ , vol. 50, no. 8, pp. 999–1011, Aug 2003.
  * [20] Yang, Liu, _et al._ , “A new ec–pc threshold estimation method for in vivo neural spike detection,” _Journal of Neural Engineering_ , vol. 9, no. 4, p. 046017.
  * [21] Zhou, Wu, _et al._ , “On the robustness of ec–pc spike detection method for online neural recording,” _Journal of Neuroscience Methods_ , vol. 235, pp. 316 – 330, 2014.
  * [22] Zou, Liu, _et al._ , “A 100-channel 1-mw implantable neural recording ic,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 60, no. 10, pp. 2584–2596, Oct 2013.
  * [23] Lotte, Congedo, _et al._ , “A review of classification algorithms for eeg-based brain–computer interfaces,” _Journal of Neural Engineering_ , vol. 4, no. 2, p. R1.
  * [24] Lotte, Bougrain, _et al._ , “A review of classification algorithms for eeg-based brain–computer interfaces: a 10 year update,” _Journal of Neural Engineering_ , vol. 15, no. 3, p. 031005.
  * [25] Serruya, Hatsopoulos, _et al._ , “Instant neural control of a movement signal.” _Nature_ , no. 6877, pp. 141–2.
  * [26] Taylor, Tillery, _et al._ , “Direct cortical control of 3d neuroprosthetic devices,” vol. 296, no. 5574, pp. 1829–1832, 2002.
  * [27] Fraser, Chase, _et al._ , “Control of a brain-computer interface without spike sorting,” _Journal of Neural Engineering_ , vol. 6, no. 5, 2009.
  * [28] Trautmann, Stavisky, _et al._ , “Accurate estimation of neural population dynamics without spike sorting,” _bioRxiv_ , p. 229252.


# Adaptation to criticality through organizational invariance in embodied
agents

Miguel Aguilera  Dept. of Computer Science, University of Zaragoza, Zaragoza,
Spain  IAS-Research Center for Life, Mind, and Society, University of the
Basque Country, Donostia-San Sebastián, Spain  Institute for Cross-
Disciplinary Physics and Complex Systems, Palma, Spain<EMAIL_ADDRESS>Manuel G. Bedia  Dept. of Computer Science, University of Zaragoza, Zaragoza,
Spain  Aragón Institute of Engineering Research, Zaragoza, Spain

###### Abstract

Many biological and cognitive systems do not operate deep within one or other
regime of activity. Instead, they are poised at critical points located at
phase transitions in their parameter space. The pervasiveness of criticality
suggests that there may be general principles inducing this behaviour, yet
there is no well-founded theory for understanding how criticality is generated
at a wide span of levels and contexts. In order to explore how criticality
might emerge from general adaptive mechanisms, we propose a simple learning
rule that maintains an internal organizational structure from a specific
family of systems at criticality. We implement the mechanism in artificial
embodied agents controlled by a neural network maintaining a correlation
structure randomly sampled from an Ising model at critical temperature. Agents
are evaluated in two classical reinforcement learning scenarios: the Mountain
Car and the Acrobot double pendulum. In both cases the neural controller
appears to reach a point of criticality, which coincides with a transition
point between two regimes of the agent’s behaviour. These results suggest that
adaptation to criticality could be used as a general adaptive mechanism in
some circumstances, providing an alternative explanation for the pervasive
presence of criticality in biological and cognitive systems.

## Introduction

Generally, the behaviour of biological and cognitive systems is not steadily
poised at one phase or another. Instead, living beings operate around points
of critical activity, at the boundary separating ordered and disordered
dynamics. Critical activity, or _criticality_ , refers to a distinctive set of
properties found at this transition. Some of these properties include the
presence of a wide range of scales of activity and maximal sensitivity to
external fluctuations [1, 2], facilitating systems at criticality to present
optimal responses when facing complex heterogeneous environments [3]. The
surprising fact is that, unlike unanimated matter where critical transitions
from order to disorder take place by fine-tuning of the parameters of the
system [1], living systems appear to be ubiquitously poised near critical
points [4]. For instance, signatures of criticality have been detected in
neural cultures [5], immune receptor proteins [6], the network of genes
controlling morphogenesis in fly embryos [7] and bacterial clustering [8].
Indicators of critical behaviour have also been observed in the brain [9] and
cognitive behavioural patterns [10].

Although these results suggest that general theoretical principles might
underlie biological self-organization, there is no well-founded theory yet for
understanding how living systems operate near critical points in a broad range
of contexts. This compels us to ask what type of mechanisms are driving
biological systems at a dauntingly diverse span of levels of organization to
operate near critical points of activity. During the last couple of decades,
the issue has been popularized as the ‘adaptation to the edge of chaos’ [11]
and different solutions have been tested through modelling approaches.
However, as we have stated, this question is still unresolved and a general
mathematical framework for understanding how living systems are driven to
criticality is yet lacking.

If we had to synthesize the approaches that have been proposed so far to
explain the presence of criticality spreading over such a wide range of
natural systems, we could broadly consider two general positions. On the one
hand, we find approaches assuming that biological systems are self-tuned,
either by learning or evolution, to regions of the parameter space displaying
optimal fitness, and that these optimal points are often placed near critical
points due to the functional advantages of critical behaviour[12, 3]. In other
words, criticality is understood as a by-product derived from the adaptive or
survival processes of living systems. On the other hand, we find different
views focused on the idea of self-organized criticality systems (SOC), in
which criticality emerges spontaneously from simple local interactions,
without fine tuning of the parameters of the model. Typically, SOC models
exploit clever local rules (e.g. in cellular automata) that produce critical
behaviour in a specific context [13, 14, 15].

Nevertheless both approaches present some weaknesses. Understanding
criticality in biological systems as an indirect consequence of adaptation to
external circumstances is not very explanatory and does not provide ways to
test alternatives. In general, explanations that simply assume that they are
the consequence of the ‘survival of the fittest’ are often unfalsifiable. With
regard to the second approach, it is well known that many SOC models are
highly idealized and, in many cases, they are not able to capture the basic
interactions of living systems, often failing to provide general explanations
of the ubiquitous emergence of criticality[16]. For that reason, although both
adaptive self-tuning and SOC models present interesting insights, they are
generally applicable in a relatively narrow set of contexts or under highly
idealized conditions.

In this paper, we explore an alternative approach to examine how a system can
display critical activity in a wide variety of situations. We propose a model
that, using only general local mechanisms, is aimed to adaptively maintain the
behaviour of the system around a critical point of its parameter space by
maintaining certain relational invariants. In other words, instead of thinking
about criticality as a by-product of adaptation to complex environments or a
spontaneous property of certain systems, we inquire into the possibility that
biological systems might are equipped with adaptive mechanisms aimed to
preserve an internal equilibrium near critical points while they interact with
their environment. Among other outcomes, the existence of cheap learning
mechanisms maintaining the parameters of a system around regions of
criticality could drastically reduce the cost of searching large parametric
spaces for finding fit solutions, or even generate interesting solutions in an
unsupervised manner.

The paper is structured as follows. First, we propose a novel method that
appears to be capable of driving a system near critical points by maintaining
certain relational invariants. In our case, these invariants are extracted
from the correlation structure of well-known model operating at a critical
point, where correlations scale with distance according to a power law
function. Second, we propose a simple learning rule maintaining this
correlation invariance, and hypothesize that it could be used for driving
systems in different contexts to operate near critical points. We test the
model using two classical examples of learning and control: the Mountain Car
and a double pendulum. We show evidence suggesting that the general rule
proposed here is able to drive adaptive agents with no free parameters towards
critical points of operation. At the same time, the agents themselves are
poised at points of behavioural transition, where they are able to exploit a
broad span of dynamic possibilities available in their environment, suggesting
a link between an internal search of critical points and the exploration of
external behavioural points of transition. Finally, we suggest further tests
of criticality and discuss the limitations and possible generalization of this
synthetic approach as a contribution towards understanding deeper principles
governing biological and cognitive systems.

## Model

Inspired by the ideas described above, we present a novel simple mechanism
designed to test whether a general adaptive system is able to drive a neural
controller near criticality by imposing certain patterns in the organizational
structure of the system. This focus on the system’s organization instead of
the mechanistic properties of its components is supported by the existence of
well-known universality classes that provide a unified expression for families
of systems operating under criticality [17].

In physics, the concept of universality allows to group a great variety of
different critical phenomena into a small number of universality classes in
such a way that all systems belonging to a given universality class are
essentially identical near the critical point. Thus, systems belonging to the
same universality class, even if defined by very different material parameters
or physical properties, have the same critical exponents characterizing
diverging observables. For example, in different spin and percolation models,
we find that the family of all bidimensional lattices (square, triangular,
hexagonal and so forth) spatial correlations follow the asymptotic form
$c(r)\propto 1/r^{\eta}$ near the phase transition point, where $\eta$ is the
same for all lattice structures of dimension 2 in a particular model.

This surprising property provides a perspective on criticality in terms of
universal relations, suggesting that we could model criticality using simple
and non-specific mechanisms independently of the individual parameters of the
system. Our hypothesis is the following: if all systems belonging to the same
universality class present the same distribution of correlations at
criticality, adjusting an arbitrary system to reproduce the same distribution
of correlations might drive the system to a similar critical point. If, as we
have said, in the neighborhood of critical points, critical exponents assume
the same universal values for a particular class, it could be enough to use in
our analysis a very simple (but nontrivial) model. Looking for generality, we
use the least structured statistical model (i.e. a maximum entropy model) of a
network of interacting units, constrained only by pairwise correlations
between them. This is known as the Ising model in physics or the Boltzmann
Machine in computer science [18]. The interest of using it is that it is also
one of the simplest models of criticality that that can be solved
analytically[19].

An Ising model can be specified as a neural network of $N$ binary variables
only constrained by pairwise correlations. Units can have a value of $+1$ or
$-1$ and are affected by local bias $h_{i}$ and couplings $J_{ij}$ between
pairs of units. These parameters take continuous values and we assume
couplings to be symmetric with $J_{ii}=0$. In order to simulate the behaviour
of the model, units are updated sequentially in a random order using Glauber
dynamics, by which each unit is activated with a probability that follows a
sigmoid function:

$P(s_{i}(t+1))=\Big{[}1+e^{-\beta 2H_{i}(t)s_{i}(t+1)}\Big{]}^{-1}$ (1)

where $H_{i}(s)=h_{i}+\sum_{j}J_{ij}s_{j}(t)$ is the effective field received
by neuron $i$ summing the internal field and inputs from other neurons (and
$2H_{i}(s)s_{i}(t+1)$ is the energy difference required to flip the sign of
unit $i$). The state of the model will be updated by sequentially applying
Glauber dynamics (i.e. Equation 1) each simulation step to all units of the
network in a random order. When updated sequentially, an Ising model with
symmetric couplings will reach an equilibrium maximum entropy distribution:

$P(s)=\frac{1}{Z}exp\Bigg{[}{\beta(\sum_{i}h_{i}s_{i}+\sum_{i<j}J_{ij}s_{i}s_{j}})\Bigg{]}$
(2)

where the distribution follows an exponential family
$P(s)=\frac{1}{Z}e^{-\beta E(s)}$ and $Z$ is a normalization value. The energy
$E(s)$ of each state of size $N$ is defined in terms of the bias $h_{i}$ and
couplings $J_{ij}$ between pairs of units, with $\beta=1/({Tk_{B}})$, being
$k_{B}$ the Boltzmann’s constant and $T$ the temperature of the system.
Without loss of generality we can set an operating working temperature such
that $\beta=1$.

Following the intuition introduced above about universality classes, we are
looking for a model that preserves certain structure in the correlations of
the system. There is some experimental evidence showing that, given an Ising
model near a critical point, one could build a family of models by learning
correlations drawn at random from the original system, which will be poised
near a critical point [20]. Inspired by this idea, we propose to reproduce and
support criticality by maintaining a distribution of correlations of a
particular universality class. Interestingly, the Ising model is a well-
studied example of a universality class. In the case of bidimensional
lattices, pairwise correlations follow the asymptotic form $c(r)\propto
1/r^{\eta}$, where $\eta=1/4$ [2], and $r$ is the distance between units.

However, instead of restricting a model to a particular set of mechanisms or a
given topology, we decide to design a learning rule that preserves a
distribution of correlations of a critical point belonging to a specific
universality class. This could capture some of the properties of that
universality class without choosing a specific topology or parametrical
configuration. The next goal would be to test whether this adaptive mechanism
has the capability of driving an arbitrary system to a critical point. And, in
case it were so, the system poised near a critical point should display
interesting features of adaptive behaviour as maximal sensitivity or a wide
dynamic range of behaviours.

In order to test this idea we design a simple learning rule to adjust the
parameters of an arbitrary Ising model to the desired distribution of
correlations. In a nutshell, the learning rule will operate as follows: 1) The
distribution of correlations of a finite square lattice Ising model is
calculated, 2) the new model is defined by assigning each neuron reference
correlation values for its synapses randomly sampled from the previous
distribution, 3) during learning, each neuron sorts its synapses by their
correlation strength, and adjusts these correlations to the reference values
assigned using an inverse Ising learning rule. We proceed now to explain these
points in detail.

First, since the size of the models employed here is far from the
thermodynamic limit, instead of directly using the diverging asymptotic form
$c(r)\propto 1/r^{\eta}$, we approximate it by computing the correlation
structure of a finite model operating at a known critical temperature. One of
the few cases where the Ising model presents an exact solution is a model with
zero fields and a bidimensional square lattice connectivity, in which a
critical point appears at $J=log{(1+\sqrt{2})}/(2\beta)$ [19]. Exploiting
this, we build a 20x20 square lattice Ising model operating at critical
temperature with periodic boundary conditions to generate reference
correlations to be used by the learning rule. We simulate the model using
Glauber Dynamics, generating $10^{6}$ samples, after an initial run of
$10^{5}$ updates starting from a random state. From this simulation, we obtain
the distribution of correlations in the system $P(c_{ij})$, where
$c_{ij}=\langle s_{i}s_{j}\rangle$, observed in Figure 1. Since the fields
$h_{i}$ of all units are zero, the means $m_{i}=\langle s_{i}\rangle$ of all
units are also zero.

Figure 1: Distribution of correlation values used for learning, generated from
a 20x20 lattice Ising model at critical temperature.

Second, once the distribution of correlations has been obtained, we generate
new models by adjusting their correlations to match the distribution
$P(c_{ij})$. For doing so applying only local information, we assign each
neuron a set of correlation values $c^{*}_{ik},k=1...N$ drawn at random from
$P(c_{ij})$. At each step, we will compute the actual correlations between a
neuron $i$ and its neighbours $j$ as $c^{m}_{ij}$, and generate reference
values for learning by sorting $c^{*}_{ik}$ to match the order of
$c^{m}_{ij}$. We will denote the sorted reference values as $c^{*}_{ij}$. All
the values of $m^{*}_{i}$ will be set equal to zero. The reason for sorting
the values of the correlations is to give more flexibility to the rule, since
we are only interested in maintaining a $c(r)\propto 1/r^{\eta}$ relation,
independently of which connection holds each value.

For the third step, the problem is that it is not trivial finding which
combination of $h_{i}$ and $J_{ij}$ generates a specific combination of
$m_{j}$ and $c_{ij}$. This is known as the ‘inverse Ising problem’, which can
be solved by using a simple gradient descent rule [18]:

$\begin{split}h_{i}\leftarrow h_{i}+\mu(m^{*}_{i}-m_{i}^{m})\\\
J_{ji}\leftarrow J_{ji}+\mu(c^{*}_{ij}-c_{ij}^{m})\end{split}$ (3)

where $\mu$ is a constant learning rate, $m^{*}_{i}$ and $c^{*}_{ij}$ are the
reference mean and correlations of the learning algorithm, and $m_{i}^{m}$ and
$c_{ij}^{m}$ are the mean and correlations of the model for the current values
of $h_{i}$ and $J_{ij}$. Generally, performing each learning step is
computationally expensive, since it requires summing over all possible states
of $s$, although approximate methods such as Monte Carlo sampling are
generally used to speed up learning. Similarly, we compute the approximate
correlations by simulating our networks using the Glauber dynamics in Equation
1 for a number of steps.

Now, before jumping to the main results, we will present a test of the model
under idealized conditions and describe the embodied agents we use for
evaluating the performance of the presented learning rule.

### Testing adaptation to criticality in isolated networks

As a demonstration of our method, we apply the learning rule to 10 different
networks for sizes $N=4,8,16,32,64$ assigning them means and correlations
drawn at random from the distribution found for the 20x20 lattice Ising model.
For each network, we apply Equation 3 for learning $m^{*}_{i}$ and
$c^{*}_{ij}$, estimating the actual $m_{i}^{m}$ and $c_{ij}^{m}$ with Glauber
dynamics.

In order to simplify the process, we have made some operative decisions. For
instance, since precision of learning is not important (the objective is to
capture the overall distribution), we do not wait for convergence of the
algorithm and simply update the learning rule $1000$ times. We use a learning
rate $\mu=0.01$ and compute $1000N$ samples for each learning step, being $N$
the size of the system. For simplicity also, instead of the critical
temperature of the lattice Ising model, we set an arbitrary inverse
temperature of $\beta=1$. Note that the choice of operating temperature is
irrelevant, since it only implies a rescaling of the parameters that the
algorithm will compensate.

For each network, we test if it is near a critical point by computing its heat
capacity. A divergence in the heat capacity is a sufficient indicator of a
continuous phase transition indicating the presence of a critical point with
maximal sensitivity to parametrical changes [1]. In our model, the heat
capacity is represented by $C(\beta)=-\beta\frac{\partial
H}{\partial\beta}=\beta^{2}(\langle E^{2}(s)\rangle-\langle E(s)\rangle^{2})$,
where $E(s)=-\sum_{i}h_{i}s_{i}-\sum_{i<j}J_{ij}s_{i}s_{j}$ is the energy of
the Ising model and $H=-\sum_{s}P(s)\log{P(s)}$ is the entropy of the system.
So we test if a system is at criticality by using the entropy $H$ of the
system as an order parameter and looking for continuous phase transitions
associated with critical points. We detect a continuous phase transition if
the entropy presents a sharp but continuous transition in which the derivative
of the entropy (the heat capacity) diverges as the system size increases.

Figure 2: (A) Divergence of the heat capacity in 10 models after learning
random correlations sampled from Figure 1. Maximum and minimum values are
shown by the grey area. (B) Distribution of the values of the coupling matrix
$J_{ij}$ of a 64 units Ising model after learning correlations sampled from
Figure 1. (the order of the nodes in the coupling matrix has been arranged
using hierarchical clustering). Values of $J_{ij}$ correspond to the values in
the colour bar.

We simulate each network for $10^{5}$ steps for different values of $\beta$,
and we find that all the heat capacity of the 10 networks diverges at the
operating temperature $\beta=1$ (Figure 2.A), showing values similar to those
of the original lattice Ising model with periodic boundaries. Let us point out
that although the distribution of correlations is similar to the lattice Ising
model, the structure of the network is radically changed. Instead of the
original ordered structure of a uniform lattice, we now have a disordered
distribution of couplings $J_{ij}$ (Figure 2.B), including both positive and
negative values. Also, each execution of the learning algorithm yields a
completely different arrangement of values of couplings $J_{ij}$.

In the following section, we test the capacity of this learning rule for
driving the neural controller of an embodied agent towards a critical point.
In order to do so, we need to take into account the environment during
learning. If we consider two interconnected Ising models (one being the neural
controller and other being the environment) Equation 3 holds perfectly if we
only apply it to the values of $i$ and $j$ corresponding to units of the
neural controller. In our case, we do not use an Ising model as an environment
but instead we use two classical examples in reinforcement learning with the
goal of testing a more realist scenario. Therefore, our learning rule will be
valid as long as the statistics of the environment can be approximated by an
Ising model with an arbitrary number of units. Luckily, Ising models in the
form of Boltzmann machines are universal approximators [21] and the stationary
distribution of any arbitrary environment can be approximated by an equivalent
Ising model.

Figure 3: (A) Structure of the embodied neural controller for a model with
$N=12$ units and $N_{h}=6$ hidden units. (B) Mountain Car environment: an
under-powered car that must drive up a steep hill by balancing itself to gain
momentum. (C) Acrobot environment: an agent has to balance a double pendulum
to reach the high part of the environment.

### Embodied model: Mountain Car and Acrobot

In order to evaluate the behaviour of the proposed learning rule, we test it
in two embodied situations using the OpenAI Gym toolkit [22]. We define a
neural network consisting of an Ising model defined as in Equation 2,
describing a network of $N=6+H_{h}$ units, with $N_{h}$ hidden units, $2$
motor units and $4$ sensor units. Motor units define the actions performed by
the agents. In sensor units, the magnetic field of the unit is not a fixed
parameter but it is be updated with the value of an external input
$h_{i}=I_{i}$. Sensor units and motor units are only connected to hidden
neurons, while hidden neurons are connected to all other neurons (Figure 3.A).
We choose this configuration because it is widely used in neural networks and
allowed recurrence in the connectivity of hidden units, although the
architecture choice is not restrictive. All units are assigned an reference
$c^{*}_{ik}$ values (selected at random from distribution $P(c_{ij})$ shown in
Figure 1), and all units except sensor units are assigned an objective mean
value $m^{*}_{i}=0$. During learning, the agent applies the rule in Equation 3
for adjusting its means and correlations to the assigned values. Each
simulation step, the units of the Ising model are updated in a sequential
random order using Glauber dynamics (Equation 1).

The first embodiment of the network consists of the Mountain Car environment
[23]. This environment is a classical testbed in reinforcement learning
depicting an under-powered car that must drive up a steep hill (Figure 3.B).
Since gravity is stronger than the car’s engine, the vehicle must learn to
gain potential energy by driving to one hill before the car is able to make it
to the goal at the top of the opposite hill (see Methods). The neural network
receives the speed of the system as an input and controls the force of the
car’s engine as an output. The second embodiment consists of a double pendulum
or ‘Acrobot’ [24], which has to coordinate the movements of two connected
links to lift its weight (Figure 3.C, see Methods). The neural network
receives the speed of the first pendulum and controls the torque applied on
the joint between the two pendulums.

For sizes $N_{h}=1,2,4,8,16,32,64$, we train 10 agents in both environments,
applying the learning rule from Equation 3, with $\eta=0.01$. Note that agents
during learning have no other explicit goal other than adjusting the
correlations of the system to a random sample extracted from the probability
distribution in Figure 1. In Supplementary Videos S1 and S2 we can observe an
example of the behaviour of agents with $N=64$ after training. In Figure
S2.C-D we can see the distribution of correlations of an agent with $N_{h}=64$
(the result is similar for all agents and sizes) and Figure S2.E-F shows the
error between the distribution in Figures S2.C-D and 1 for this agents. In
Figures S2.C-D we can observe how agents after learning display a correlation
structure consistent with the correlation distribution $c(r)\propto
1/r^{\eta}$ of the lattice Ising model at criticality.

In the following section, we analyse whether the agents designed in this way
are able to operate in their environment in a regime of criticality.

## Results

In this section, we analyse the behaviour of the neural controllers and the
behavioural patterns of the agents with respect to the possibilities of their
parameter space. The goal is to test if the learning rule proposed here is
effective for driving the agent near critical points. As we observe below, the
10 agents display quite similar behaviour for each environment, despite the
fact that each one has learned different values of $c^{*}_{ij}$ and
$J^{*}_{ij}$. We are interested in analysing if there is anything special
about the configuration reached by the agents after learning.

In order to compare the agents with other behavioural possibilities, we
explore the parameter space by changing the parameter $\beta$ of the agents.
Modifying the value of $\beta$ is equivalent to a global rescaling of the
parameters of the agent transforming $h_{i}\leftarrow\beta\cdot h_{i}$ and
$J_{ij}\leftarrow\beta\cdot J_{ij}$, thus exploring the parameter space along
one specific direction. That is, changing $\beta$ is just a way of testing one
dimension in the parameter space of possible models. First, we assess the
presence of criticality in the neural controller of the agent. Specifically,
we look for the presence of a continuous phase transition in which an order
parameter of the system (the entropy) presents a sharp transition displaying a
divergence of its derivative (the heat capacity). Second, we analyze not only
the behaviour of the neural controller but the agent as a whole, looking for
behavioural transitions of the agents and divergences of the susceptibility of
the agent’s behaviour to parametrical changes.

### Signatures of criticality in the neural controller

In order to test whether the agents are being driven towards a critical point,
we analyse signatures of critical behaviour in the neural controller of the
agent. As we mentioned before, a sufficient indicator for criticality is the
presence of a divergence of the heat capacity of the system (as in Figure
2.A). A divergence in the heat capacity indicates the presence of a continuous
phase transition presenting a critical point in which the system is maximally
sensitive to parametrical changes. Unfortunately, when the neural controller
is embodied in the Mountain Car and Acrobot environments, we can no longer
access a formal description of the probability distribution of the agent-
environment system, thus we cannot directly compute from the energy of the
system values as the entropy or heat capacity of the system. Nevertheless, we
can still directly compute the entropy $H(x)=-\sum_{x}P(x)\log{P(x)}$ of any
variable of the system by estimating its probability distribution $P(x)$
through simulations. In order to compute the entropy and the heat capacity of
different variables, we simulate the agent’s behaviour for $101$ values of
$\beta$, logarithmically distributed in the interval $[10^{-1},10^{1}]$. We
run the 10 agents for each embodiment during $10^{6}$ simulation steps,
reseting the agent’s position and state every $5\cdot 10^{4}$ simulation
steps.

First, we compute the entropy of the probability function of hidden neurons in
the controller. Due to computational constraints (the number of states of the
probability distribution increases with $2^{N_{h}}$) we only compute the
entropy in agents up to $N_{h}=16$ hidden units. Displaying the entropy of the
neurons for different values of $\beta$ we observe that the agents are near an
order-disorder transition (Figure 4.A,B). Larger sizes make the transition
sharper and closer to the operating temperature $\beta=1$. From the entropy of
a variable, its heat capacity can be computed as
$C(\beta)=-\beta\frac{\partial H(\beta)}{\partial\beta}$. From the computed
$101$ values of entropy, $H(\beta)$ is estimated by fitting a curve using
cubic B-splines [25] as indicated in the Methods section. In Figure 4.C,D we
observe how the system displays a similar divergence of the heat capacity of
the neural network as the Ising model, suggesting that the robot’s neural
controller is at a critical point.

In order to confirm a divergence in the observed transition larger systems
should be evaluated. Since evaluating the entropy of the hidden neurons is
computationally infeasible for large sizes, we repeat the analysis for the
state of the 4 sensor units of the network for a larger amount of hidden
units. In Figure 5 we show the entropy and heat capacity of the sensor units
for sizes up to $N_{h}=64$ hidden neurons, where we can observe a similar
picture than in Figure 4, suggesting that the heat capacity diverges, and a
second order transition takes place in the neural controller of the agent.
These results suggest that the agent’s neural controller is operating near a
point of criticality, resembling a continuous phase transition, indicating
that the neural controller self-organizes to present maximal sensitivity to
changes in its parameter space.

Figure 4: (A, B) Entropy of the agent’s neural controller for different sizes
up to $N=22$ ($N_{h}=16$), for 10 different agents and 101 values of $\beta$.
(C, D) Heat capacity of the agent’s neurons, computed as
$C(\beta)=-\beta\frac{\partial H(\beta)}{\partial\beta}$, where $H(\beta)$ is
estimated using B-splines. Maximum and minimum values are shown by the grey
area. The figures suggest that in both the Mountain Car (left) and Acrobot
(right) embodiments the model presents a divergence of its heat capacity as
the number of neurons $N$ increases, suggesting that the neural controller of
the system is near a continuous phase transition. Figure 5: (A, B) Entropy of
the sensor units for different sizes up to $N=70$ ($N_{h}=64$), for 10
different agents and 101 values of $\beta$. (C, D) Heat capacity of the
agent’s neurons, computed as $C(\beta)=\beta\frac{\partial
H(\beta)}{\partial\beta}$, where $H(\beta)$ is estimated using B-splines.
Maximum and minimum values are shown by the grey area. The figures suggest
that in both the Mountain Car (left) and Acrobot (right) embodiments the
system presents a divergence of its heat capacity as the number of neurons $N$
increases, suggesting that the neural controller of the system is near a
second order phase transition.

In addition to the presence of a continuous phase transition, a classical
signature of criticality is the presence of power law distributions in the
statistical descriptions of the states of a system. Unlike an isolated Ising
model, our neural network is connected to an environment and the probability
distribution of the Ising neural controller is no longer described by Equation
2. Thus, we compute the probability distribution of the system by counting the
occurrence of each state of the units $s$ to compute the probability
distribution of the Ising model $P(s)$. We simulate the system at $\beta=1$
for $10^{8}$ simulation steps. As in training, the agent’s position and state
are reset every $5\cdot 10^{4}$ simulation steps. We observe that all agents
approximately follow Zipf’s law for the Mountain Car (Figure 6.A) and Acrobot
embodiments (Figure 6.B), with error bars in a very narrow range. The power-
law distribution of neural activation patterns suggests that the neural
controller of the agents is operating near a critical point. We have to note
that the sole occurrence of a power law is generally insufficient to assess
the presence of criticality and may arise naturally in some non-equilibrium
conditions. Nevertheless, together with the apparent divergence of the heat
capacity it supports the idea that the neural controller of the agents might
be poised near a critical state.

Figure 6: Ranked probability distribution function of the neural network for
10 different agents and different sizes in (A) the Mountain Car and (B) the
Acrobot embodiments. The real distribution is compared with a distribution
following Zipf’s law, (i.e. $P(s)=1/rank$, dash-dotted line). We observe a
good agreement between the model and Zipf’s law, suggesting critical scaling.

### Behavioural transitions in the parameter space

What does it imply for the agent to adapt to be poised near a critical point?
It should be remarked here that agents are given no explicit goal. They only
tend to adapt to behavioural patterns maintaining a distribution of
correlations randomly sampled from the distribution shown in Figure 1. To
explore this issue, we examine the different behavioural modes of the agent
while exploring its parameter space by changing the value of $\beta$. The
behaviour of the Mountain Car can be described just by its horizontal position
$x$ and speed $v$ at different moments of time. As well, the horizontal and
vertical positions of the tip of the Acrobot’s links is a good description of
its behaviour.

In Figure 7.A-C we can observe the behaviour of a Mountain Car agent with
$N_{h}=64$ for $\beta=\\{0.75,1,1.2\\}$, respectively. We observe that for
values of $\beta$ lower than the operating temperature, the agents are not
able to reach the top of the mountain. On the other hand, when $\beta$ is
higher, the agents present more ‘rigid’ trajectories going from one mountain
peak to the other. At $\beta=1$ the agent is able to reach the top of the
mountain (note that the peaks of the mountain are located at $x=-\pi/2$ and
$x=\pi/6$) while displaying larger behavioural diversity. Similarly, in Figure
7.D-F, we observe that the Acrobot agent with $N_{h}=64$ at $\beta=1$ displays
a diverse range of behaviours, being able to reach the top of the plane while,
when $\beta$ is lowered or increased, it drifts to other behavioural modes in
less diverse regimes. Although only one agent is represented for each
environment, the results of Figure 7 are similar in all agents and sizes.

Figure 7: Transition in the behavioural regime of the agents with $N_{h}=64$.
We show the behaviour of two individual agents with different values of
$\beta$ for the Mountain Car (A, B, C) and Acrobot (D, E, F) embodiments. We
observe that $\beta=1$ is a transition point between two modes of behaviour in
both agents.

To get a more general picture of different agents and sizes, we can analyse
the behavioural transitions in relevant variables of the agent environment
systems. Furthermore, we analyse whether the behaviour of the agent, and not
only its neural controller, is near a critical point. For inspecting this, we
calculate the mean height of agents $\langle y\rangle$ in our simulations
(Figure 8.A,B), and compute the susceptibility of this height value as
$\chi_{y}(\beta)=\beta\frac{\partial\langle y\rangle}{\partial\beta}$ (Figure
8.C,D). In this case, the susceptibility appears to increase monotonically
with size when size is doubled, even if these increases are not as uniform as
in Figures 4 and 5. Although further tests of criticality could validate if
criticality is also found in the whole agent-environment system, the figures
suggest that the continuous phase transition of the neural controller
corresponds with a sharp transition in the agent’s behaviour.

Figure 8: (A, B) Mean height $\langle y\rangle$ of the agents for different
sizes up to $N=70$ ($N_{h}=64$), for 10 different agents and 101 values of
$\beta$. (C, D) Susceptibility the agent’s behaviour, computed as
$\chi_{y}(\beta)=\beta\frac{\partial\langle y\rangle}{\partial\beta}$, where
$\langle y\rangle$ is estimated using B-splines. Maximum and minimum values
are shown by the grey area. The figures suggest that in both the Mountain Car
(left) and Acrobot (right) embodiments the behaviour of the agent presents a
sharp transition around $\beta=1$.

## Discussion

Recapitulating the main ideas presented so far, we have tested how, by taking
a set of correlations chosen at random from a distribution generated by a
lattice Ising model at a critical point, we can construct a new model that
appears to be also near a critical point in its parameter space. Moreover, if
an embodied agent maintains these correlations using a simple learning rule
while interacting with its environment – as a sort of organizational
homeostasis – the agent neural controller seems to be driven to a critical
point, which coincides with behavioural transitions in the agent-environment
parameter space. Due to computational limitations for estimating the
probability distribution of the embodied Ising networks, we only calculate the
entropy and heat capacity of neural controllers up to $N_{h}=16$ hidden
neurons (and $N=22$ total neurons) and approximate larger sizes computing the
entropy and heat capacity of $N_{s}=4$ sensor neurons for neural controllers
up to $N=70$ total neurons. These results are still far from the thermodynamic
limit and further tests should confirm the results presented here for larger
sizes. Nevertheless, in all cases we observe a clear diverging tendency of the
heat capacity every time the size of the system is doubled. Tests for larger
sizes could be performed by designing environments that can be described by a
Gibbs distribution, avoiding the combinatorial explosion in computing the heat
capacity of the system by calculating it to be directly from the energy of the
whole agent-environment system.

These results suggest the possibility that criticality could be diagnosed and
induced directly from the maintenance of a given distribution of correlations
rather than modelling a precise mechanistic structure. Also, our results show
that criticality could be generated by quite simple mechanisms only relying on
local information, maintaining specific correlations around a given value.
Here we have implemented the mechanism as a simple Boltzmann Learning process,
but other rules could have the same effect, such as the combination of Hebbian
and anti-Hebbian tendencies in specific ratios or other simple mechanisms. In
our model, we only require the system to maintain a distribution of relations
between the components of the system. This connects with systemic approaches
to biology interested not in specific or intrinsic components of biological
systems but in the networks of relations and processes [26, 27, 28]. It is
also in line with notions of relational invariance in Piaget’s approach to
functional invariants in cognitive development [29] or Maturana and Varela’s
ideas of autopoietic machines, defined as homeostatic systems that maintain
their own organization constant as a network of relations between components
[30].

Assuming a similar systemic perspective, we have derived learning rule
intended to drive a system towards critical points by maintaining an invariant
structure of correlations roughly defined by a critical exponent $1/r^{\eta}$.
Our approach assumes a different point of view on self-organized criticality,
in which the distribution of correlations is not the consequence of
criticality in a specific topology but the cause driving an indeterminate
topology to what appears to be a critical point. The question now is whether
imposing connections derived from a $1/r^{\eta}$ function is a strong
assumption or implies particularly exigent circumstances. We do not think so,
since power law functions can be naturally generated by simple rules of
preferential attachment favouring ‘rich-get-richer’ cumulative inequalities
[31], or directly as a natural consequence of certain geometries of space
(e.g. gravitation laws[32]).

Our model only assumes that a system is going to adapt in order to preserve an
internal network of relations. It emphasizes the maintenance of organizational
structures capable of reproducing the behaviour of living systems without
relying on internal models of the external source of sensory input. This
contrasts with other approaches which have focused on understanding
criticality as a strategy to effectively represent a complex and variable
external world [3], for example studying criticality in predictive coding or
deep learning architectures dealing with complex inputs [33, 34]. In those
cases, an internalist view is assumed, where the neural controller represents
structures of an external world, whose complexity may be the cause of critical
activity in the neural controller. Instead, our approach is agnostic in terms
of the inputs or the external world of an organism, and deals only with how an
agent rearranges its internal structures facing different environments.

The agents presented here are not specifically designed for a particular
problem. In simple terms, our agents generate (preserving the same internal
neural organization) a wide variability and richness of behaviours (avoiding
both disorder and explosive and indiscriminate propagation) that permits them
to explore the space of parameters and eventually achieve solutions that they
were not designed to find. The empirical evidence of experiments shown here
supports this idea. A parallel could be established with the concept of play,
which can be understood as a ‘rule-breaker’ activity of the constraints of a
stable and self-equilibrating regime of behaviours which has no concrete goals
[35]. A model as the one presented here could be used for exploring life-like
autonomous behaviour without the need for explicit internal representations,
goals, or rule-based behaviour. Conceptual models of critical activity based
on the maintenance of a system’s relational invariants could help in the
development of a synthetic route towards the exploration of adaptive and
embodied criticality.

## Methods

#### Mountain Car.

This environment consists in a car with mass $m$ moving along a one-
dimensional environment. In this environment, the agent moves its position in
an horizontal axis $x$, limited to an interval of $[-1.5\pi/3,0.5\pi/3]$. Each
horizontal position represents a point in an environment with two mountains,
whose height is defined as $y=0.55+0.45\sin(3x)$. The velocity in the
horizontal axis is updated each time step as
$v(t+1)=v(t)+0.001a-0.0025\cos(3x)$, where $a$ is the action of the motor
which can be either ${-1,0,1}$, impulsing the car with a force $F=ma$. The
inputs $I_{i}$ fed to the sensor units are defined as an array of 4 units,
which are assigned the instantaneous velocity of the car discretized into an
array of 4 bits. Each input $I_{i}$ is assigned a value of $1$ if its
corresponding bit is active and $-1$ otherwise.

#### Acrobot.

The Acrobot is a two-link planar robot composed of two pendulums joined at
their tip, with a motor applying a torque in clockwise or counterclockwise
directions on the joint between the two links. The position of the system is
defined by the angles of both pendulums $\theta_{1}$ and $\theta_{2}$, whose
behaviour is defined by the following system of differential equations:

$\begin{split}\ddot{\theta}_{1}=-(d_{2}\ddot{\theta}_{2}+\phi_{1})/d_{1}\\\
\ddot{\theta}_{2}=(m_{2}l_{c2}^{2}+I_{2}-\frac{d_{2}^{2}}{d_{1}})^{-1}(\tau+\frac{d_{2}}{d_{1}}\phi_{1}-\phi_{2})\\\
d_{1}=m_{1}l_{c1}^{2}+m_{2}(l_{1}^{2}+l_{c2}^{2}+\\\
+2l_{1}l_{c2}\cos(\theta_{2}))+I_{1}+I_{2}\\\
d_{2}=m_{2}(l_{c2}^{2}+l_{1}l_{c2}\cos(\theta_{2}))+I_{2}\\\
\phi_{2}=m_{2}l_{c2}g\cos(\theta_{1}+\theta_{2}-\pi/2)\\\
\phi_{1}=-m_{2}l_{1}l_{c2}\dot{\theta}_{2}^{2}\sin(\theta_{2})-\\\
-2m_{2}l_{1}l_{c2}\dot{\theta}_{2}\dot{\theta}_{1}\sin(\theta_{2})+\\\
+(m_{1}l_{c1}+m_{2}l_{1})g\cos(\theta_{1}-\pi/2)+\phi_{2}\end{split}$ (4)

where $\tau$ is the torque applied to the system which can be either
${-1,0,1}$, $m_{1}=m_{2}=m$ is the mass of the links, $l_{1}=l_{2}=1$ is the
length of the links and $l_{c1}=l_{c2}=0.5$ are the lengths to the center of
mass of the links, $I_{1}=I_{2}=1$ are the moments of inertia of the links and
$g=9.8$ is the gravity. As well, variables $d_{1}$ and $d_{2}$ are the total
moments of inertia of each link, and $\phi_{1}$ and $\phi_{2}$ are linked to
the potential energy of the system

Similarly to the Mountain Car, the inputs fed to the sensor units in the
Acrobot embodiment are defined as an array defined with 4 sensor units,
encoding the angular speed of the first link with binary values (encoding
active and inactive bits as +1 and -1).

#### Task difficulty.

In order to make the tasks challenging, we set the maximum velocity allowed to
the Mountain car to $\pm 0.045$ (typically is set to $0.07$) and the mass of
the Acrobot links to $m=1.75$ (typically $m=1$). These parameters are designed
to make it difficult for agents controlled by neural networks with random
parameters solve the task (reaching the top of the environment), having
success rates of $6.1\%$ for the Mountain Car and $3.1\%$ for the Acrobot.
Success rates were evaluated by simulating 1000 neural controllers with random
parameters (sampled from a uniform distribution in the range $[-2,2]$). The
Mountain Car was simulated for 1000 simulation steps starting from a random
position in the valley between $[0.4,0.6]$, and was considered successful when
reached the maximum position at least once. The Acrobot was simulated for 5000
simulation steps from the bottom position (angles and angular speeds between
$[-0.1,01]$) and was consider successful if reached a vertical position higher
than $1.8$.

#### Training.

During training, agents are initialized in the starting random positions in
the bottom of their environments ($x\in[0.4,0.6]$ and $v=0$ for the Mountain
Car and $\theta_{1},\theta_{2},\dot{\theta}_{1},\dot{\theta}_{2}\in[-0.1,01]$
for the Acrobot). The state of the neural network is randomized and the
initial parameters $h_{i}$ and $J_{ij}$ are set to zero. Agents are simulated
for $1000$ trials of $5000$ steps, computing each trial the values of
$m_{i}^{m}$ and $c_{ij}^{m}$ and applying Equation 3 at the end of the trial.
The agent’s position and state are reset every $5\cdot 10^{4}$ simulation
steps

#### Code availability.

The source code implementing the learning rule in the different examples is
freely available at https://github.com/MiguelAguilera/Adaptation-to-
criticality-through-organizational-invariance.

#### Curve interpolation.

The heat capacity of hidden and sensor units was computed by interpolating
entropy curves using cubic B-splines [25] using _scipy_ function _splrep_ with
as smoothing coeficient of 1. Different coeficients were tested with similar
results.

## References

  * [1] Salinas, S. R. A.  Phase Transitions and Critical Phenomena: Classical Theories.  In _Introduction to Statistical Physics_ , Graduate Texts in Contemporary Physics, 235–256 (Springer New York, 2001).  DOI: 10.1007/978-1-4757-3508-6_12.
  * [2] Salinas, S. R. A.  Scaling Theories and the Renormalization Group.  In _Introduction to Statistical Physics_ , 277–304 (Springer New York, 2001).
  * [3] Hidalgo, J. _et al._ Information-based fitness and the emergence of criticality in living systems.  _Proceedings of the National Academy of Sciences_ 111, 10095–10100 (2014).  DOI 10.1073/pnas.1319166111.
  * [4] Mora, T. & Bialek, W.  Are biological systems poised at criticality?  _Journal of Statistical Physics_ 144, 268–302 (2011).
  * [5] Schneidman, E., Berry, M. J., Segev, R. & Bialek, W.  Weak pairwise correlations imply strongly correlated network states in a neural population.  _Nature_ 440, 1007–12 (2006).  DOI 10.1038/nature04701.
  * [6] Mora, T., Walczak, A. M., Bialek, W. & Callan, C. G.  Maximum entropy models for antibody diversity.  _Proceedings of the National Academy of Sciences_ 107, 5405–5410 (2010).  DOI 10.1073/pnas.1001705107.
  * [7] Krotov, D., Dubuis, J. O., Gregor, T. & Bialek, W.  Morphogenesis at criticality.  _Proceedings of the National Academy of Sciences_ 111, 3683–3688 (2014).  DOI 10.1073/pnas.1324186111.
  * [8] Chen, X., Dong, X., Be’er, A., Swinney, H. L. & Zhang, H. P.  Scale-Invariant Correlations in Dynamic Bacterial Clusters.  _Physical Review Letters_ 108, 148101 (2012).  DOI 10.1103/PhysRevLett.108.148101.
  * [9] Chialvo, D. R.  Critical Brain Dynamics at Large Scale.  In Plenz, S. I. D. & Niebur, E. (eds.) _Criticality in Neural Systems_ , 43–66 (Wiley-VCH Verlag GmbH & Co., 2014).
  * [10] Van Orden, G., Hollis, G. & Wallot, S.  The blue-collar brain.  _Fractal Physiology_ 3, 207 (2012).  DOI 10.3389/fphys.2012.00207.
  * [11] Kauffman, S. A.  _The origins of order: Self-organization and selection in evolution_ (Oxford University Press, USA, 1993).
  * [12] Beggs, J. M. & Plenz, D.  Neuronal Avalanches in Neocortical Circuits.  _Journal of Neuroscience_ 23, 11167–11177 (2003).
  * [13] Bak, P., Tang, C. & Wiesenfeld, K.  Self-organized criticality.  _Physical Review Letters_ 59, 381–384 (1987).  DOI 10.1103/PhysRevLett.59.381.
  * [14] Bak, P. & Sneppen, K.  Punctuated equilibrium and criticality in a simple model of evolution.  _Physical Review Letters_ 71, 4083–4086 (1993).  DOI 10.1103/PhysRevLett.71.4083.
  * [15] Jensen, H. J.  _Self-organized criticality: emergent complex behavior in physical and biological systems_ (Cambridge Univ. Press, Cambridge [u.a., 1998).
  * [16] Watkins, N. W., Pruessner, G., Chapman, S. C., Crosby, N. B. & Jensen, H. J.  25 Years of Self-organized Criticality: Concepts and Controversies.  _Space Science Reviews_ 198, 3–44 (2016).  DOI 10.1007/s11214-015-0155-x.
  * [17] Kadanoff, L. P.  More is the same; mean field theory and phase transitions.  _Journal of Statistical Physics_ 137, 777–797 (2009).
  * [18] Ackley, D. H., Hinton, G. E. & Sejnowski, T. J.  A learning algorithm for Boltzmann machines.  _Cognitive science_ 9, 147–169 (1985).
  * [19] Onsager, L.  Crystal Statistics.  _Physical Review_ 65, 117–149 (1944).  DOI 10.1103/PhysRev.65.117.
  * [20] Tkacik, G., Schneidman, E., Berry II, M. J. & Bialek, W.  Spin glass models for a network of real neurons.  _arXiv:0912.5409_ (2009).
  * [21] Montúfar, G. F.  Universal Approximation Depth and Errors of Narrow Belief Networks with Discrete Units.  _Neural Computation_ 26, 1386–1407 (2014).
  * [22] Brockman, G. _et al._ OpenAI gym.  _arXiv:1606.01540_ (2016).
  * [23] Moore, A. W.  Efficient memory-based learning for robot control (Univ. of Cambridge, 1990).
  * [24] Sutton, R. S.  Generalization in Reinforcement Learning: Successful Examples Using Sparse Coarse Coding.  In _Advances in NIPS, 8_ , 1038–1044 (MIT Press, 1996).
  * [25] Dierckx, P.  _Curve and Surface Fitting with Splines_ (Oxford University Press, Inc., New York, NY, USA, 1993).
  * [26] Bernard, C. & Greene, H. C.  _An introduction to the study of experimental medicine_ (Courier Corporation, 1957).
  * [27] Rosen, R.  Some relational cell models: The metabolism-repair systems.  In _Foundations of Mathematical Biology_ , 217–253 (Academic Press, 1972).  DOI: 10.1016/B978-0-12-597202-4.50011-6.
  * [28] Ashby, W. R.  Principles of the self-organizing system.  _Principles of Self-organization_ 255–278 (1962).
  * [29] Piaget, J.  _Biology and knowledge: an essay on the relations between organic regulations and cognitive processes_ (University of Chicago Press, 1974).  Google-Books-ID: sXMSHUIBDNsC.
  * [30] Maturana, H. R. & Varela, F. J.  _Autopoiesis and Cognition: The Realization of the Living_ (D. Reidel Publishing Company, 1980), 1st edn.
  * [31] Greenwood, M. & Yule, G. U.  An inquiry into the nature of frequency distributions representative of multiple happenings with particular reference to the occurrence of multiple attacks of disease or of repeated accidents.  _Journal of the Royal statistical society_ 83, 255–279 (1920).
  * [32] Barrow, J. D.  New Dimensions.  In _The Constants of Nature_ , 203–205 (Random House, New York, 2002).
  * [33] Friston, K., Breakspear, M. & Deco, G.  Perception and self-organized instability.  _Frontiers in Computational Neuroscience_ 6 (2012).  DOI 10.3389/fncom.2012.00044.
  * [34] Lin, H. W. & Tegmark, M.  Critical Behavior from Deep Dynamics.  _arXiv:1606.06737_ (2016).  ArXiv: 1606.06737.
  * [35] Di Paolo, E. A., Rohde, M. & De Jaegher, H.  Horizons for the enactive mind: Values, social interaction, and play.  _Enaction_ 33–87 (2010).

## Acknowledgements

Research was supported in part by the Spanish National Programme for Fostering
Excellence in Scientific and Technical Research project PSI2014-62092-EXP, by
project TIN2016-80347-R funded by the Spanish Ministry of Economy and
Competitiveness and the UPV/EHU post-doctoral training program ESPDOC17/17.

## Author contributions statement

M.A. conceived and conducted the experiments. M.A. and M.B. analysed the
results and wrote the manuscript.

## Competing interests

The authors declare no competing interests.