Criticality as It Could Be: organizational invariance as self-organized criticality in
embodied agents
Miguel Aguilera1,2,3 and Manuel G. Bedia1,2

arXiv:1704.05255v3 [nlin.AO] 24 May 2017

1

Dept. of Computer Science, Univ. of Zaragoza, Zaragoza (Spain)
2
Aragón Institute of Engineering Research, Zaragoza (Spain)
2
Institute for Cross-Disciplinary Physics and Complex Systems, Palma (Spain)
sci@maguilera.net
Abstract
This paper outlines a methodological approach for designing adaptive agents driving themselves near points of criticality. Using a synthetic approach we construct a conceptual
model that, instead of specifying mechanistic requirements to
generate criticality, exploits the maintenance of an organizational structure capable of reproducing critical behavior. Our
approach exploits the well-known principle of universality,
which classifies critical phenomena inside a few universality classes of systems independently of their specific mechanisms or topologies. In particular, we implement an artificial embodied agent controlled by a neural network maintaining a correlation structure randomly sampled from a lattice
Ising model at a critical point. We evaluate the agent in two
classical reinforcement learning scenarios: the Mountain Car
benchmark and the Acrobot double pendulum, finding that in
both cases the neural controller reaches a point of criticality,
which coincides with a transition point between two regimes
of the agent’s behaviour, maximizing the mutual information
between neurons and sensorimotor patterns. Finally, we discuss the possible applications of this synthetic approach to
the comprehension of deeper principles connected to the pervasive presence of criticality in biological and cognitive systems.

Introduction
Biological systems at a wide range of scales – from protein
families to brains – show signatures of criticality. These systems do not operate deep into one or other regime of activity
instead of they operate at or near critical points, poised at
transitions in their parameter space (Mora and Bialek, 2011).
In a nutshell, criticality refers to a distinctive set of properties found at the boundary separating regimes with different
dynamics: the transition between an ordered and a disordered phase. Some of these properties include power-law divergences of some quantities described by critical exponents
and fractal behaviour (Salinas, 2001). Signatures of criticality have been detected in neural cultures (Schneidman et al.,
2006), immune receptor proteins (Mora et al., 2010), the
network of genes controlling morphogenesis in fly embryos
(Krotov et al., 2014) or lipid membranes (Honerkamp-Smith
et al., 2009). As well, indicators of critical behaviour have

been observed in the human brain (Chialvo, 2014) and cognitive behavioural patterns (Van Orden et al., 2012). These
results hint at general theoretical principles underlying biological self-organization, compelling us to ask what type of
mechanisms are driving biological systems at a dauntingly
diverse range of levels of organization to operate near critical points of activity.
This question is largely unresolved. Firstly, because the
connection between experimental descriptive indicators of
criticality and models is often fragile (Wagenmakers et al.,
2012). On the other hand, models of criticality are generally
used at the level of analogy, based on specific mechanistic requirements reproducing criticality in a set of particular
cases, providing a myriad of different models but often failing to capture explanations of a few more universal classes
and properties. During the last decade a new generation of
experiments has attempted to go beyond analogies generating models directly inferred from biological data (Mora and
Bialek, 2011). However, the difficulty of inferring general
principles from the specific modelled mechanisms is still
patent in these models.
In order to progress in the comprehension of critical phenomena, there are some alternatives and suggestions that we
could explore with more emphasis. For example, we could
try to detect universal mechanisms able to generate critical
activity in a wide set of systems and contexts. This approach
– a kind of ‘Artificial Life route’ to self-organized criticality
– could promote the development of conceptual models explaining how organisms are driven towards criticality in an
abstract level, working as ‘proofs of concept’ (Barandiaran
and Chemero, 2009) and supporting existing and future experimental findings. Specifically, here we focus on the relation of criticality with organizational invariants of the model.
By this, we mean that instead of specifying specific mechanistic requirements to generate criticality, we explore the
possibility of generating criticality through the conservation
of certain invariants in the relations between the components
of the system (e.g. correlations between components), i.e.
imposing certain patterns in the system’s organization. This
organizational point of view, in contrast with focusing on

intrinsic properties of the components of the system, is supported by the existence of well-known universality classes
that provide a unified expression for families of systems operating under criticality (Kadanoff, 2009). Systems under
the same universality class, even if they are defined by very
different material parameters or physical properties, present
the same critical exponents characterizing diverging observables which are defined by the symmetries of the system.
In line with these ideas, the paper is structured as follows.
First, we propose a conceptual model based on statistical
mechanics to design an artificial agent that maintains organizational invariants in its structure, specifically a distribution of correlations between the components of the system. Operatively, the model is implemented as a ‘Boltzmann machine’ neural network reproducing a correlation
structure randomly sampled from a 2D lattice Ising model
at a critical point. Then, we test the model in two classical examples of learning and control: the Mountain Car and
the double pendulum. In both examples we find that agents
with no free parameters exploit the whole dynamic range of
available configurations being poised near critical transition
points of their parameter space between qualitatively different behavioural regimes. Finally, we discuss the possible
applications of this synthetic approach to contribute to the
comprehension of the deeper principles that governs biological and cognitive systems.

Organizational invariants of self-organized
criticality
The development of mechanistic models of criticality in biological systems – based on finding the particular properties and generative mechanism that give support this regime
of behaviour – typically assumes a particular way of modelling. Typically, we can find models of criticality exploiting clever local rules generating critical behaviour (e.g. Bak
et al., 1987; Bak and Sneppen, 1993), parameter tuning of
systems showing critical phase transitions (e.g. Beggs and
Plenz, 2003; Kitzbichler et al., 2009) and/or the use of specific structural properties of the underlying network of the
system (e.g. Rubinov et al., 2011; Haimovici et al., 2013).
However, as we have pointed out above, families of critical phenomena can be classified into universality classes determined only by a few properties of the system. One of
these properties is that all the critical exponents of all models within a given universality class are exactly the same.
For example, in all Ising models in 2D lattices (square, triangular, hexagonal and so forth) spin-spin correlations follow
the asymptotic form c(r) ∝ 1/rη , where η = 1/4 (Salinas,
2001). This surprising property provides a perspective about
criticality in terms of universal relations, suggesting that we
could model criticality using simple and non-specific mechanisms independently on the topology of the system.
Following this intuition, we propose a model not committed to a particular local behavioural rule or connectivity

topologies. Instead, our model is built on the maintenance of
an organizational structure for reproducing the global functional properties of a family of critical models. There exists some experimental evidence showing that, given a Ising
model near a critical point, one could build a family of models by learning correlations drawn at random from the original model, which will be poised near a similar critical point
(Tkacik et al., 2009). Based on this inspiration, we propose
to reproduce and support criticality by the invariant maintenance of an organizational structure of correlations following a 1/rη law. Thus, instead of restricting a model to a particular set of mechanisms, we find a general organizational
distribution which could be easily implemented by very different structures, driving a system to a regime of criticality.
We define our model as a neural network as an Ising
model or Boltzmann Machine (Ackley et al., 1985) following a maximum entropy distribution:
"
#
X
X
1
hi si +
Jij si sj
(1)
P (s) = exp β
Z
i
i<j
where the distribution follows a en exponential family
P (s) = Z1 e−βE(s) , Z is a normalization value, where the
energy E(s) of each state is defined in terms of the bias hi
and couplings Jij between pairs of units, and β = 1/(T kB ),
being kB Boltzmann’s constant and T the temperature of the
system. Units si can take discrete values of +1 or −1 and
the couplings and bias can take continuous values. Without
loss of generality we can set an operating working temperature such that β = 1. To simulate the network, we use
Glauber dynamics, by which the value of a unit is inverted
each time with probability:
#−1

"
β∆Ei (s)

Pi (s) = 1 + e

(2)

P
where ∆Ei (s) = 2(hi si + j Jij si sj ) is the energy difference between the inverted and original state. Throughout
the paper, we simulate the network updating its state sequentially by applying Glauber dynamics to all units in the network in a random order at each simulation step.
In order to induce criticality in its behaviour, our model
maintains a specific internal structure of correlations, conserving the statistical properties of a system that we know
it is at criticality. Since the size of our models will be far
from the thermodynamic limit, instead of directly using the
asymptotic form c(r) ∝ 1/rη , we approximate it by computing the correlation structure of a known model showing
this distribution in the thermodynamic limit. One of the few
specific cases where the Ising model presents an exact solution is a model with zero fields and 2D lattice connectivity,
in
√
which a critical point appears at J = βlog(1 + 2)/2 (Onsager, 1944). Specifically, we use a 20x20 2D lattice Ising
model at critical temperature with periodic boundary conditions. We simulate the model using Glauber Dynamics,

A
hi ← hi + µ(mi − mm
i )

30

Jji ← Jji + µ(cij − cm
ij )

P (cij ) 20
10
0

0.50

0.55

0.60

0.65

cij

0.70

B
N=64
N=32
N=16
N=8
N=4

1.5

C(β)

1.0

0.5

0.0

6 × 10−1

2 × 100 3 × 100

100

β

C

D
Jij

0

0.40

Jij

0

0.8

0.35

0.6

0.30

20

0.25

i

20

0.4

i

0.20

40

0.15

0.2

40

0.0

0.10
60

0.05
0

25

50

j

0.00

60
0

25

50

−0.2

j

Figure 1: (A) Distribution of correlation values used for
learning, extracted from a 20x20 lattice Ising model at critical temperature. (B) Divergence of the heat capacity in 10
models learning random correlations sampled from Figure
1.A. Maximum and minimum values are shown by the grey
area. (C) Distribution of the coupling matrix of a 8x8 lattice
Ising model with periodic boundaries. (D) Distribution of
a 64 units Ising model learning correlations sampled from
Figure 1.A. The order of the nodes in the coupling matrix
has been rearranged using hierarchical clustering.
generating 106 samples, after initializing the model running
105 updates. From this simulation, we obtain the distribution of correlations in the system cij = hsi sj i observed in
Figure 1.A. Since the fields at all units are zero, the means
mi = hsi i of all units are also zero.
Once obtained a distribution of correlations, we will generate new models by assigning to each unit and pair of units
means and correlations randomly drawn at random from the
distribution of the 20x20 Ising model. At this point, the
problem is that it is not trivial finding which combination of
hi and Jij generates a specific combination of mj and cij .
This is known as the ‘inverse Ising problem’. This problem
was formulated by Ackley et al. (1985) in their discussion of
Boltzmann machines and can be solved by a simple gradient
descent rule:

(3)

where µ is a constant learning rate, mi and cij are the objective mean and correlations of the learning algorithm, and
m
mm
i and cij are the mean and correlations of the model for
the current values of hi and Jij . Computing each learn step
is computationally costly, since it requires to sum over all
possible states of s, although approximate methods as Monte
Carlo sampling are generally used to speed up learning. As a
demonstration, we apply the learning rule to different models assigning them objective means and correlations drawn
at random from the distribution found for the 20x20 lattice Ising model. For each model, we apply Equation 3 for
m
learning mi and cij , computing the actual mm
i and cij with
Glauber dynamics. 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 used a learning rate
µ = 0.01 and computed 1000N samples for each learning
step, being N the size of the system.
Using this method, we apply the learning rule to 10 different models for sizes N = 4, 8, 16, 32, 64. For each
model, we test if the models are at criticality by computing their heat
C(β) = β 2 hE 2 (s)i−hE(s)i2 , where
P capacityP
E(s) = − i hi si − i<j Jij si sj is the energy of the Ising
model. The divergence of the heat capacity is an indicator of
criticality. We simulate each model for 105 steps for different values of β, and we found that all the 10 models diverge
at the operating temperature of β = 1 (Figure 1.B), showing
similar values of heat capacity to the original lattice Ising
model with periodic boundaries. Nevertheless, although the
distribution of correlations is similar than the lattice Ising
model, the structure of the model is radically changed. Instead of the original ordered structure of a uniform lattice
(Figure 1.C), we have now a disordered distribution of couplings Jij (Figure 1.D), including both positive and negative couplings. Also, each random selection of correlations
yields a completely different arrangement of values of couplings Jij .
In the following section, we use this learning rule to drive
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 will not use an
Ising model as an environment but instead we will use two
classical examples from reinforced learning. 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 are universal

approximators (Montúfar, 2014) and any arbitrary environment could be approximated by an equivalent Ising model.
A
Sensor
neurons

Hidden
neurons
Motor
neurons

B

C

Figure 2: (A) Structure of the embodied neural controller.
(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 double
pendulum
In order to evaluate the behaviour of the proposed learning
model, we test it in two embodied situations using the OpenAI Gym toolkit (Brockman et al., 2016). For doing so we
define a neural network consisting on an Ising model defined as in Equation 1, describing a network of 8 units, with
6 hidden units and 2 motor units, and a variable number of
sensor units (4 or 6, depending on the environment). Motor units will define the actions performed by the agents,
while sensor units are updated directly from the state of the
environment. Sensor units and motor units are only connected to hidden neurons, while hidden neurons are fully
connected to the rest of the system (Figure 2.A). All connections are assigned an objective cij value (selected at random from distribution P (cij ) at Figure 1.A), and all units
except sensor units are assigned an objective mean value of
mi = 0. During learning, the agent will apply the rule in
Equation 3 for adjusting its means and correlations to the
objective values. At each simulation step, the Ising model

is simulated by updating its units in uniform random order using Glauber dynamics. The first embodiment of the
network consists in the Mountain Car environment (Moore,
1990). This environment is a classical testbed in reinforcement learning depicting an under-powered car that must
drive up a steep hill (Figure 2.B). Since gravity is stronger
than the car’s engine, the vehicle must learn to leverage potential energy by driving to the opposite hill before the car
is able to make it to the goal at the top of the rightmost
hill. In this environment, the horizontal position x of the
car is limited to an interval of [−1.5π, 0.5π], and the vertical position of the car is defined as y = 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. The sensors of the neural network are defined as an array of four
units, which are fed the instantaneous velocity of the car,
discretized into an array of 4 bits. Each sensor unit is assigned a value of 1 if its corresponding bit is active and −1
otherwise.
The second embodiment consists in a double pendulum or
‘Acrobot’ (Sutton, 1996) which has to coordinate the movements of two connected links to lift its weight (Figure 2.C).
The position of the system is defined by the angles of both
pendulums θ1 and θ2 , whose behaviour is defined by the following system of differential equations:
θ̈1 = −(d2 θ̈2 + φ1 )/d1

d2
d22 −1
) (τ + φ1 − φ2 )
d1
d1
2
2
d1 = m1 lc1 + m2 (l12 + lc2
+

2
θ̈2 = (m2 lc2
+ I2 −

+2l1 lc2 cos(θ2 )) + I1 + I2
d2 =

2
m2 (lc2

+ l1 lc2 cos(theta2 )) + I2

(4)

φ2 = m2 lc2 g cos(θ1 + θ2 − π/2)

φ1 = −m2 l1 lc2 θ̇22 sin(θ2 )−
−2m2 l1 lc2 θ̇2 θ̇1 sin(θ2 )+

+(m1 lc1 + m2 l1 )g cos(θ1 − π/2) + φ2

where τ is the torque applied to the system which can be
either −1, 0, 1, m1 = m2 = m is the mass of the links, l1 =
l2 = 1 is the length of the links and lc1 = lc2 = 0.5 are the
lengths to the center of mass of the links, and I1 = I2 = 1
are the moments of inertia of the links and g = 9.8 is the
gravity.
The Acrobot embodiment was defined with 6 sensor units,
divided in two groups of 3. Each group encoded the binarized horizontal and vertical positions of the tip of the
second link, defined as x = sin(θ1 ) + sin(θ1 + θ2 ), y =
− cos(θ1 ) − cos(θ1 + θ2 ).
In order to make the tasks challenging, we set the maximum velocity allowed to the Mountain car to ±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 11.2%
for the Mountain Car and 7.7% for the Acrobot1 .
We train 10 agents for each embodiment applying the
learning rule from Equation 3, with η = 0.01. In order to avoid overfitting, we add an L2 regularization term
λ = 0.004. Note that agents during learning have no explicit
goal other than adjusting the correlations of the system to a
random sample extracted from the probability distribution
in Figure 1.A. Agents are initialized in the starting random
positions in the bottom of their environments (x ∈ [0.4, 0.6]
for the Mountain Car and θ1 , θ2 , θ̇1 , θ̇2 ∈ [−0.1, 01] for the
Acrobot). The state of the neural network is randomized and
the initial parameters hi and Jij are set to zero. Agents are
simulated for 1000 trials of 5000 steps, computing each trial
m
the values of mm
i and cij and applying Equation 3 at the end
of the trial . Note that agents are not reset at the end of the
trial.

Results
In this section, we analyze the behaviour of the neural controllers and the behavioural patterns of the agents respect to
the possibilities of their parameter space. The first striking
result is that all 10 agents present quite similar behaviour for
each embodiment, despite the fact that each one has learned
different values of cij and Jij . In order to compare the
agents with other behavioural possibilities, we explore the
parameter space by changing the parameter β of the agents.
Modifying the value of β is equivalent to a global rescaling of the parameters of the agent transforming hi ← β · hi
and Jij ← β · Jij , thus exploring the parameter space along
one specific direction. For 21 values of β logarithmically
distributed in the interval [10−1 , 101 ] we simulate the 10
agents for each embodiment during 106 simulation steps, after starting the agents from the random starting position and
an initial run of 104 simulation steps. We will use the results
of those simulations for all the situations in this section.

Signatures of criticality in the neural controller
First, in order to test whether the agents are being driven near
a critical point, we analyze signatures of critical behaviour
in the neural controller of the agent. Counting the occurrence of each possible state of all the neurons of the agents
1

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.

(including sensor, hidden and motor neurons) we can compute the probability distribution of the Ising model P (s).
We observe that all agents approximately follow a Zipf’s
law at β = 1 for the Mountain Car (Figure 3.A) and Acrobot embodiments (Figure 3.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.
Since the neural network is now connected to an environment that is driven by deterministic mechanics, the probability distribution of the Ising neural controller is no longer
described by equation 1. Thus, classical indicators of criticality as the divergence of the heat capacity are not directly
calculable from the energy of the model. However, we can
look for other indicators related to second order phase transitions, computed from the probability distribution P (s) calculated from simulation. One indicator can be the behaviour
of order parameters such as the entropy of the system, a
transition around the critical temperature. For example, if
we compute the entropy of the
P probability function of the
neural controller H(s) = − x P (x) log P (x) for different values of β we observe that the agent is near an orderdisorder transition (Figure 3.C,D) (as well, some Mountain
Car agents present a smaller transition at larger values of β,
showing that interesting behaviours can also arise in the ordered phase).
In lattice Ising models, critical points are also indicated
by other measures such as a peak in the mutual information
found in the system (Barnett et al., 2013). Mutual information is defined as:
X
P (x, y)
(5)
I(x, y) =
P (x, y) log
P (x)P (y)
x,y
In our case, we measure the mutual information between the
set of sensor units ss and the set containing the rest of the
neurons of the controller sn . The objective is to capture how
much information of the sensorimotor coordination of the
robot is captured by the system as a whole instead of being contained in the variables alone. In Figure 3.E,F we can
observe how mutual information between sensors and neurons has a peak very close to β = 1, suggesting that agents
are poised in a point of the parameter space maximizing the
exchange of information between the agent and its environment.
These results suggest that the agent’s neural controller is
operating near a point of criticality, resembling a second order phase transition. As we will see now, no only the agent’s
neural controller presents indicators of critical activity, but
also the behaviour of the agent as a whole.

Behavioural transitions in the parameter space
What does it imply for the agent to poise its neural controller
at a critical point? It should be remarked here that our agents
are given no explicit goal but they only tend to behavioural

A

B

A

100

P (s)10

P (s)10−2
10−3

0.02

−2

v

10−4
103

100

rank

102

0.02

v

0.00

−0.02

−1.0

−0.5

0.0

104

−0.04
−1.5

0.5

x

−1.0

−0.5

0.0

−1.5

0.5

x

−1.0

−0.5

0.0

0.5

2000

3000

4000

1

2

7500

10000

x

rank
0

0.0

C

x

−1

−1.0

8

0

x

x−0.5

D
8

0

1000

2000

3000

−1
0

4000

1000

2000

3000

4000

0

1000

t

t

6

H(s)

0.00

−0.02

−0.04
−1.5

102

0.04

0.02

v

0.00

−0.04

10−4

β=2

0.04

−0.02

101

C
β=1

0.04

10−1

100

B
β = 0.5

100

t

6

H(s)
4

D

2

2

2

0 −1
10

0 −1
10

4

100

101

β

100

101

β

E

F
0.6
0.4

I(ss; sn) 0.4

I(ss; sn)

E
β = 0.5

F
β=1

2

1

1

y 0

y 0

y 0

−1

−1

−1

−2
−2

−2
−2

−1

0

1

2

x

−1

0

1

2

x

−2
−2

0.3

2

2

2

0.2

y 0

y 0

y 0

0.2

β=2

2

1

−1

0

x

0.1
0.0 −1
10

100

β

101

0.0 −1
10

−2

100

101

0

2500

5000

t

7500

10000

−2

0

2500

5000

7500

t

10000

−2

0

2500

5000

t

β

Figure 3: Signatures of criticality for 10 different agents in
the Mountain Car (left) and Acrobot (right) embodiments.
(A-B) Ranked probability distribution function of the Ising
models. The real distribution is compared with a distribution
following Zipf’s law, (i.e. P (s) = 1/rank, dashed line).
We observe a good agreement between the model and Zipf’s
law, suggesting critical scaling. (C-D) Entropy of the neural system H(s). A transition is observed near the operating temperature. (E-F) Mutual information between sensors
and neurons I(ss ; sn ). In all figures the mean of all agents
is shown as a solid line, while the region between maximum and minimum values showed by the different agents
is shown as a grey area.

patterns maintaining a distribution of correlations randomly
sampled from the distribution in Figure 1.A. Thus, we explore the effects of transiting the critical point of the neural
controller observing the different behavioural modes of the
agent in the parameter space by changing the value of β. The
behaviour of the car can be described just by the position x
and speed v at different moments of time. As well, the position of the tip of the Acrobot’s links shows a good image of
the system’s behaviour.
In Figure 4.A-C we can observe the behaviour of the
Mountain Car for β = {0.5, 1, 2} respectively. We observe
that for values of β lower than the operating temperature, the
agents are not able to reach the top of the mountain. On the
other hand, when β is higher, the agents present more ‘rigid’
trajectories going form one top of the mountain to the other.

Figure 4: Transition in the behavioural regime of the agents.
We show the behaviour of two agents with different values
of β for the Mountain Car (A, B, C) and Acrobot (D, E, F)
embodiments. We observe that β = 1 is a transition point
between two modes of behaviour in both agents.
At β = 1 the agent is able to reach the top of the mountain
(note that the peaks of the mountain are located at x = −π/2
and x = π/6) while displaying larger behavioural diversity.
Similarly, in Figure 4.D-F, we observe that the Acrobot at
β = 1 displays a diverse range of behaviours, being able to
reach to the top of the plane while, when β is lowered or increased, it drifts to other behavioural modes in less diverse
regimes. Furthermore, if we compute the median value of
height ỹ for both agents at different values of β (Figures 5),
we observe that there is a transition in the parameter space
around β = 1, in which both agents maximize the diversity
of positions reached in their environment. This suggests that
there is a behavioural transition connected with the phase
transition of the neural Ising controller, in which the agent
maximizes the dynamic range of inputs of the neural controller.

Discussion
Recapitulating the main ideas presented so far, we have
tested how 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 is also near
a critical point of its parameter space. Moreover, imposing an embodied agent to maintain these correlations using

A

B
2
1.0

1

0.5

ỹ

ỹ

0.0

−1

−0.5
−1.0

0

10−1

100

β

101

−2

10−1

100

101

β

Figure 5: Behavioural transitions in the (A) Mountain Car
and (B) Acrobot environments, showing the median height
ỹ of the agent, together with the 5% and 95% percentiles
in the Mountain Car and the 5% and the 95% and 99.9%
percentiles in the Acrobot (dotted lines).

a simple learning rule – as a sort of organizational homeostasis – drives the agent to a critical point, which coincides
with behavioural transitions of its parameter space. This
suggests the possibility that criticality could be diagnosed
and even induced directly from the maintenance of a given
distribution of correlations rather than modelling a precise
mechanistic structure. Also, criticality could be caused 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,
as the combination of Hebbian and anti-Hebbian tendencies
in specific ratios.
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
(Bernard and Greene, 1957; Rosen, 1972; Ashby, 1962), and
it is also in line with notions of relational invariance as Piaget’s approach on functional invariants in cognitive development (Piaget, 1971) or Maturana and Varela’s ideas of
autopoietic machines, defined as homeostatic systems that
maintain constant its own organization as a network of relations between components (Maturana and Varela, 1980).
Assuming a similar systemic perspective, we have derived
learning rules for a system that drives itself near a critical
point by maintaining an invariant structure of correlations
roughly defined by a critical exponent 1/rη . This promotes
a different perspective on criticality. In our model, the distribution of correlations is not the consequence of criticality
in a specific topology, but the cause driving an indeterminate topology to a critical point. The question now could be
whether imposing connections derived from a 1/rη 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

(Greenwood and Yule, 1920), or directly as a natural consequence of certain geometries of space (as e.g. gravitation
laws, see Barrow, 2002).
In this way, our model only assumes that a system is going to adapt to preserve an internal network of relations.
It emphasizes the maintenance of organizational structures
capable of reproducing living systems behaviors being in
opposition with the ones relying on internal models of the
external source of sensory input. Thus, this contrasts with
other approaches which have focused in understanding criticality as a strategy to effectively represent a complex and
variable external world, for example discovering criticality
in predictive coding or deep learning architectures dealing
with complex inputs (Friston et al., 2012; Lin and Tegmark,
2016; Hidalgo et al., 2014). In those cases, an internalist
view is assumed, where the neural controller is representing
structures of the external world, whose complexity may be
the cause of criticality being present in the neural controller.
Instead, our approach is agnostic about what are 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 to explore the space of parameters and eventually to
achieve solutions of problems for which it was not designed.
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 stable and self-equilibrating regime of behaviours, that is not directly required from the environment
and without concrete goals (Di Paolo et al., 2010). A model
as the one presented here could be used for exploring lifelike autonomous behaviour without the need of explicit internal representations, goals, or rule-based behaviour. Conceptual models of critical activity based in the maintenance
of a system’s relational invariants could help developing a
synthetic route towards the exploration of adaptive and embodied criticality.

Acknowledgements
Research was supported in part by the Spanish National Programme for Fostering Excellence in Scientific and Technical Research project PSI2014-62092-EXP and by the project
TIN2016-80347-R funded by the Spanish Ministry of Economy and Competitiveness.

Materials and methods
The source code implementing the learning rule in the different examples is freely available at https://github.com/
MiguelAguilera/Criticality-as-It-Could-Be.

References
Ackley, D. H., Hinton, G. E., and Sejnowski, T. J. (1985). A learning algorithm for Boltzmann machines. Cognitive science,
9(1):147–169.

Kitzbichler, M. G., Smith, M. L., Christensen, S. R., and Bullmore, E. (2009). Broadband Criticality of Human Brain
Network Synchronization. PLOS Computational Biology,
5(3):e1000314.

Ashby, W. R. (1962). Principles of the self-organizing system.
Principles of Self-organization, pages 255–278.

Krotov, D., Dubuis, J. O., Gregor, T., and Bialek, W. (2014).
Morphogenesis at criticality. Proceedings of the National
Academy of Sciences, 111(10):3683–3688.

Bak, P. and Sneppen, K. (1993). Punctuated equilibrium and criticality in a simple model of evolution. Physical Review Letters, 71(24):4083–4086.

Lin, H. W. and Tegmark, M. (2016). Critical Behavior from Deep
Dynamics. arXiv:1606.06737. arXiv: 1606.06737.

Bak, P., Tang, C., and Wiesenfeld, K. (1987). Self-organized criticality. Physical Review Letters, 59(4):381–384.
Barandiaran, X. E. and Chemero, A. (2009). Animats in the Modeling Ecosystem. Adaptive Behavior, 17(4):287–292.
Barnett, L., Lizier, J. T., Harré, M., Seth, A. K., and Bossomaier, T.
(2013). Information flow in a kinetic Ising model peaks in the
disordered phase. Physical review letters, 111(17):177203.

Maturana, H. R. and Varela, F. J. (1980). Autopoiesis and Cognition: The Realization of the Living. D. Reidel Publishing
Company, 1st edition.
Montúfar, G. F. (2014). Universal Approximation Depth and Errors of Narrow Belief Networks with Discrete Units. Neural
Computation, 26(7):1386–1407.
Moore, A. W. (1990). Efficient memory-based learning for robot
control. Univ. of Cambridge.

Barrow, J. D. (2002). New Dimensions. In The Constants of Nature, pages 203–205. Random House, New York.

Mora, T. and Bialek, W. (2011). Are biological systems poised at
criticality? Journal of Statistical Physics, 144(2):268–302.

Beggs, J. M. and Plenz, D. (2003). Neuronal Avalanches in Neocortical Circuits. Journal of Neuroscience, 23(35):11167–
11177.

Mora, T., Walczak, A. M., Bialek, W., and Callan, C. G. (2010).
Maximum entropy models for antibody diversity. Proceedings of the National Academy of Sciences, 107(12):5405–
5410.

Bernard, C. and Greene, H. C. (1957). An introduction to the study
of experimental medicine. Courier Corporation.
Brockman, G., Cheung, V., Pettersson, L., Schneider, J., Schulman, J., Tang, J., and Zaremba, W. (2016). OpenAI gym.
arXiv:1606.01540.
Chialvo, D. R. (2014). Critical Brain Dynamics at Large Scale. In
Plenz, S. I. D. and Niebur, E., editors, Criticality in Neural
Systems, pages 43–66. Wiley-VCH Verlag GmbH & Co.
Di Paolo, E. A., Rohde, M., and De Jaegher, H. (2010). Horizons
for the enactive mind: Values, social interaction, and play.
Enaction, pages 33–87.
Friston, K., Breakspear, M., and Deco, G. (2012). Perception and
self-organized instability. Frontiers in Computational Neuroscience, 6.
Greenwood, M. and Yule, G. U. (1920). 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(2):255–279.
Haimovici, A., Tagliazucchi, E., Balenzuela, P., and Chialvo, D. R.
(2013). Brain Organization into Resting State Networks
Emerges at Criticality on a Model of the Human Connectome.
Physical Review Letters, 110(17):178101.
Hidalgo, J., Grilli, J., Suweis, S., Muñoz, M. A., Banavar, J. R.,
and Maritan, A. (2014). Information-based fitness and the
emergence of criticality in living systems. Proceedings of the
National Academy of Sciences, 111(28):10095–10100.
Honerkamp-Smith, A. R., Veatch, S. L., and Keller, S. L. (2009).
An Introduction to Critical Points for Biophysicists; Observations of Compositional Heterogeneity in Lipid Membranes.
Biochimica et biophysica acta, 1788(1):53–63.
Kadanoff, L. P. (2009). More is the same; mean field theory and
phase transitions. Journal of Statistical Physics, 137:777–
797.

Onsager, L. (1944). Crystal Statistics. Physical Review, 65(34):117–149.
Piaget, J. (1971). Biology and knowledge: An essay on the relations between organic regulations and cognitive processes.
Rosen, R. (1972). Some relational cell models: The metabolismrepair systems. In Foundations of Mathematical Biology,
pages 217–253. Academic Press. DOI: 10.1016/B978-0-12597202-4.50011-6.
Rubinov, M., Sporns, O., Thivierge, J.-P., and Breakspear, M.
(2011). Neurobiologically Realistic Determinants of SelfOrganized Criticality in Networks of Spiking Neurons. PLOS
Computational Biology, 7(6):e1002038.
Salinas, S. R. A. (2001). Scaling Theories and the Renormalization
Group. In Introduction to Statistical Physics, pages 277–304.
Springer New York.
Schneidman, E., Berry, M. J., Segev, R., and Bialek, W. (2006).
Weak pairwise correlations imply strongly correlated network states in a neural population. Nature, 440(7087):1007–
12.
Sutton, R. S. (1996). Generalization in Reinforcement Learning:
Successful Examples Using Sparse Coarse Coding. In Advances in NIPS, 8, pages 1038–1044. MIT Press.
Tkacik, G., Schneidman, E., Berry II, M. J., and Bialek, W.
(2009). Spin glass models for a network of real neurons.
arXiv:0912.5409.
Van Orden, G., Hollis, G., and Wallot, S. (2012). The blue-collar
brain. Fractal Physiology, 3:207.
Wagenmakers, E. J., van der Maas, H. L. J., and Farrell, S. (2012).
Abstract Concepts Require Concrete Models. Topics in Cognitive Science, 4(1):87–93.

