# Evidence of Scaling Regimes in the Hopfield Dynamics of Whole Brain Model

Giorgio Gosti  Center for Life Nano- and Neuro-Science, Istituto Italiano di
Tecnologia, Viale Regina Elena 201, I-00161, Rome, Italy  Istituto di Scienze
del Patrimonio Culturale, Sede di Roma, Consiglio Nazionale delle Ricerche,
00010, Montelibretti (RM), Italy  Sauro Succi  Center for Life Nano- and
Neuro-Science, Istituto Italiano di Tecnologia, Viale Regina Elena 201,
I-00161, Rome, Italy  Istituto per le Applicazioni del Calcolo del Consiglio
Nazionale delle Ricerche, via dei Taurini 19, I-00185, Rome, Italy  Giancarlo
Ruocco  Center for Life Nano- and Neuro-Science, Istituto Italiano di
Tecnologia, Viale Regina Elena 201, I-00161, Rome, Italy  Dipartimento di
Fisica, Universitá di Roma “La Sapienza”, P.le Aldo Moro 5, I-00185, Roma,
Italy

###### Abstract

It is shown that a Hopfield recurrent neural network, informed by
experimentally derived brain topology, recovers the scaling picture recently
introduced by Deco Deco2020 , according to which the process of information
transfer within the human brain shows spatially correlated patterns
qualitatively similar to those displayed by turbulent flows. Although both
models employ a coupling strength which decays exponentially with the
euclidean distance between the nodes, their mathematical nature is widely
different, Hopf oscillators versus Hopfield neural network. Hence, their
convergence suggests a remarkable robustness of the aforementioned scaling
picture. Furthermore, the present analysis shows that the Hopfield model brain
remains functional by removing links above about five decay lengths,
corresponding to about one sixth of the size of the global brain. This
suggests that, in terms of connectivity decay length, the Hopfield brain
functions in a sort of intermediate “turbulent liquid”-like state, whose
essential connections are the intermediate ones between the connectivity decay
length and the global brain size. This ”turbulent-like liquid” appears to be
more spiky than actual turbulent fluids, with a scaling exponent around $2/5$
instead of $2/3$.

##  I Introduction

In a recent paper Deco2020 , it was argued that the information transfer
within the human brain may proceed in close analogy with the mechanisms which
govern mass and energy transport in turbulent fluids. The analogy refers to
the emergence of correlated spacetime patterns which would facilitate and
possibly optimize information transfer across the brain, the “same” way as
fluid turbulence enhances energy transport across fluids. Based on an extended
dynamical system of oscillators (Hopf model), supplemented by large-scale
neuro-imaging empirical data on the connections between different brain areas,
the authors of Deco2020  numerically showed that the statistical correlations
of human brain signals may exhibit a scaling regime similar to the one
observed in turbulent fluids.

At a closer scrutiny, though, the analogy does not appear to be quantitative,
meaning that the spatial scaling exponents of the second order structure
function $S_{2}(d)=<(O(d)-O(0))^{2}>\sim d^{\alpha}$, $O$ being a suitable
order parameter, are not the same, about $\alpha\sim 1/2$ for the Hopf model
of the human brain against $\alpha$=2/3 for fluid turbulence. Yet, a scaling
regime in space is found in either cases, which is intriguing per se,
regardless of the specific value of $\alpha$.

In this paper, we probe the robustness of the above analogy by performing a
similar analysis based on the same neuro-data, but using a very different
model, namely a Hopfield neural network, and measuring scaling exponents as a
function of the connectivity decay length $\delta$.

The recurrent Hopfield neural network model is a well-established model for
small-scale neuronal networks Amari1972 ; Little1974 ; Hopfield1982 . It is a
neural network model based on binary McCulloch-Pitts neurons Amit1985 ;
Brunel2016 ; Hillar2018 ; Hillar2021  in which every processing unit ($i$) is
connected to all the other ones ($j$) through a set of weights ($J_{ij}$).
Most of the literature considers fully connected and symmetric Hopfield
networks since this simplifies the mathematical analysis of the model,
although recently the investigation of diluted networks Stauffer2003 ;
Brunel2016 ; Kim2017 ; Folli2018 ; Leonetti2020  and asymmetric networks
Gopalsamy1994 ; Xu1996 ; Chen2001 ; Franca2000 ; Zheng2010 ; Szedlak2014 ;
Folli2018 ; Leonetti2020  was introduced.

Usually, a recurrent Hopfield neural network is trained to store patterns or
memories as steady state attractors with the Hebbian prescription. In this
case, the coupling matrix $J_{ij}$ assumes the form of the sum of diadic
forms, it is therefore symmetric and fully connected. The maximum storage
capacity of such a network is found to be linear with the number of nodes $N$.
The capacity can be increased by the process of learning (Hebbian learning):
in this process the simultaneous activation of two nodes will strengthen their
connection. This is inspired by biological neural networks where the
simultaneous activation of neurons leads to increments in synaptic strength
Hebb1949 , and it allows the network to store a larger number of patterns
Amit1985 , but this quantity still remains linear with $N$. Lately,
researchers have also investigated how unlearning methods can improve the
stability of the stored patterns Fachechi2019 ; Benedetti2022 ; Benedetti2023
.

Interestingly, it has been recently demonstrated both numerically Folli2018 ;
Gosti2019  and theoretically Hwang2019 ; Hwang2020  that random, non-Hebbian,
coupling matrices have a storage capacity much larger than the Hebbian ones.
Indeed the storage capacity $C$ of random matrices scales exponentially with
$N$: $C\sim e^{\Sigma N}$ being $\Sigma$ the complexity. The latter quantity
is found to be maximised for asymmetric and diluted matrices, with values that
match the connectivity found in the hippocampus areas.

The approach outlined in reference Deco2020 , as elaborated in the subsequent
chapter, diverges from both the Hebbian and random connectivity models. It
relies on a deterministic coupling matrix, which is derived from experimental
observations of the connections between human brain regions. This approach
assumes that the coupling strength is proportionally linked to a specific
function of the distances between nodes (see Eq. 3):
$J_{ij}=e^{-d_{ij}/\delta}$, $d_{ij}$ being the Euclidean distance between
nodes $i$ and $j$ in the network. At variance with the latter approach, where
the dynamics is based on the Hopf oscillator model, we here use the same
coupling matrix of Deco2020 , but employ the Hopfield dynamics that, to the
best of our knowledge, was never used to model whole brain activity dynamics.

Our main conclusions are summarized as follows:

i) We confirm the existence of a scaling regime of the neural activity in
space for each value of the connectivity decay length $\delta$ in the explored
range, 0$\leq$$\delta$$\leq$10 mm.

ii) At $\delta\sim 5.5$ mm we found $\alpha\sim 2/5$, very close to the one
reported by Deco _et al._ Deco2020  for the same value of $\delta$ 111In Ref.
Deco2020  the used parameter is the inverse of the decay scale,
$\lambda=1/\delta$..

iii) The scaling exponent shows a decreasing trend with the decay length
$\delta$. Importantly, the functional relation $\alpha(\delta)$ shows a steep
increase in the region around $\delta\sim 5.55$ mm, the value indicated in
Deco2020  as the physiologically relevant one. This means that the specific
value of the decay length has a significant impact on the scaling exponent,
for instance a slight increase to $\delta=6$ mm would yield a scaling exponent
slightly above $2/3$, the value associated to homogenoeus incompressible
turbulence. Conversely, a slight decrease to $\delta=5$ mm, would yield
$\alpha\sim 0.1$, a much less organized regime. This pinpoints the importance
of fine-tuning the value of $\delta$.

iv) By systematically pruning the couplings below a given running threshold,
$J_{ij}<J_{th}$ (which correspond to setting to zero the connection above a
certain distance), it is found that the scaling exponent remains largely
unaffected up to we get to pruning 95$\%$ \- 98$\%$ of the couplings,
suggesting that nodes at a distance $d>\sim 5\delta$, do not significantly
partake to the physical mechanisms at the roots of the scaling regime.
Accordingly, the model brain loses functionality by removing connections below
$6\delta$, which is close to the geometric mean between the local scale
$\delta$ and the size of the global brain, about $25\delta$, circa $20$ cm.

##  II The turbulent-like whole-brain model

.

Deco _et al._ model the global human brain as a set of $N$ coupled Hopf
oscillators of the form

$\displaystyle\frac{dx_{i}}{dt}=a_{i}x_{i}$ $\displaystyle+$
$\displaystyle(\beta y_{i}-x_{i})(x_{i}^{2}+y_{i}^{2})-\omega_{i}y_{i}$ (1)
$\displaystyle+$ $\displaystyle
G\sum_{j=1}^{N}J_{ij}(x_{j}-x_{i})-\nu_{i}\xi_{i}$
$\displaystyle\frac{dy_{i}}{dt}=a_{i}y_{i}$ $\displaystyle+$
$\displaystyle(\beta x_{i}+y_{i})(x_{i}^{2}+y_{i}^{2})-\omega_{i}x_{i}$ (2)
$\displaystyle+$ $\displaystyle
G\sum_{j=1}^{N}J_{ij}(y_{j}-y_{i})-\nu_{i}\xi_{i}$

In the above, $(x_{i},y_{i})$ are the real and imaginary part of the BOLD
signal recorded at the $i$-th node of the network, $a_{i}=-0.02$ are decay
terms at the brink of the Hopf bifurcation, where scaling is best observed,
$\omega_{n}$ is the natural frequency of the oscillators, $\beta$ the so
called shear factor and $\xi$ a gaussian noise with unit variance. The global
coupling across the network is represented through the connectivity matrix
$J_{ij}$ which is assumed to obey the Exponential Distance Rule

$J_{ij}=e^{-d_{ij}/\delta}$ (3)

where $d_{ij}$ is the euclidean distance between nodes $i$ and $j$ and we have
set $\delta$ ($\equiv 1/\lambda$ in Deco’s notation). Deco2020  argues that
the exponential distance rule matches very well the empirical HCP dMRI
tractography of the human brain. In particular, it reproduces the decay of the
fiber densities between the parcel centroid pairs as a function of the
Euclidian distance. The authors Deco2020  further show that the optimal fit is
obtained for $\delta=5.55$ $\mathrm{m}\mathrm{m}$.

###  II.1 Amplitude of turbulence and scaling regime

In order to assess the scaling properties of the “turbulent brain” (TB), Deco
_et al._ introduce the following complex Kuramoto Order Parameter:

$\Psi_{i}(t)\doteq
R_{i}(t)e^{i\theta_{i}(t)}\doteq\sum_{j=1}^{N}W_{ij}e^{i\phi_{i}(t)}$ (4)

where $W_{ij}=J_{ij}/\sum_{k=1}^{N}J_{ik}$ and $\phi_{i}(t)$ is the phase of
the $i$-th oscillator obtain by Hilbert-transform of the (narrow band) BOLD
signal. The authors then measure the amplitude of turbulence as the variance
of the fluctuations of amplitude of the order parameter, namely:

$D\doteq<R^{2}>-<R>^{2}$ (5)

where bracket indicate average in space and time over the entire network and
over hundreds of noise realizations.

The scaling regime is assessed by inspecting the associated structure function

$S_{2}(d)=<(R(x+d)-R(x))^{2}>\sim d^{\alpha}$ (6)

where “$x$” and “$x+d$” indicate that the two network nodes are at distance
$d$ and the average is taken at statistical steady state, over all network
nodes and noise realizations. This observable is commonly used to assess the
statistical properties of fluid turbulence, and for the case of Kolmogorov
homogeneous incompressible turbulence one obtains

$\alpha_{T}=2/3,$

indicating that the change in the (longitudinal) flow field at a distance $d$,
$\delta u(d)=\langle|u(x+d)-u(x)|\rangle$, scales like $d^{1/3},$ hence the
turbulent flow field is not differentiable (sometimes $\alpha$ is also called
the roughness exponent) Frisch1995 ; Benzi1993 . For the whole-brain model,
Deco _et al._ Deco2020  report

$\alpha_{D}\sim 1/2,$

indicating that the correlation pattern of the brain signals is less regular
than those observed in homogeneous incompressible turbulence.

As we shall show, these values should be taken with extreme caution, as they
are very sensitive to a number of assumptions, primarily the structure and
functional dependence of the connectivity matrix on the internode distance.

##  III  Discrete-time Hopfield Recurrent Neural Networks

Figure 1:  A) 1000 parcel centroids in three dimensional space from Schaefer’s
cerebral cortical parcellation atlas Schaefer2018 . B) Frequency histogram of
the Euclidean distances $d_{ij}$ between each pair of nodes from Schaefer’s
cerebral cortical parcellation atlas Schaefer2018 . The linear dependence on
$d_{ij}$ indicates that the nodes are basically distributed on a two
dimensional surface. C) Matrix composed of the Euclidean distances $d_{ij}$
between each pair of nodes from Schaefer’s cerebral cortical parcellation
atlas Schaefer2018 . D) Matrix formed by the values such that $J_{ij}>0.4$,
where $J_{ij}=e^{-d_{ij}/\delta}$ with $\delta$=0.55.

Following Amari1972 ; Little1974 ; Hopfield1982 ; Amit1985 ; Folli2017 ;
Gosti2019 ; Leonetti2020 ; Gosti2024 , we consider a network of $N$-binary
neurons that can be either firing or silent.

The state of each neuron is a binary variable $s_{i}$, with $i=1,2,\ldots,N$,
such that $s_{i}=-1$ or $s_{i}=1$, so that the network state is represented by
a binary string $\mathbf{s}=(s_{0},s_{1},\ldots,s_{N})$. Neurons interact
through the connectivity matrix $\mathbf{J}$, with matrix elements $J_{ij}$
following the same exponential dependence $J_{ij}=e^{-d_{ij}/\delta}$ used in
Deco2020 . We assume a discrete time $t$, and, at each time step $t$, the
evolution of the neuron state $s_{i}(t)$ is given by the following non-linear
dynamic equation,

$s_{i}(t+1)=sign\left[\sum_{j=1}^{N}J_{ij}s_{j}(t)\right],$ (7)

where the activation function $sign(x)$ is defined such as $sign(x)=1$ for
$x\geq 0$, and $sign(x)=-1$ otherwise. Consequently, whenever the summation of
the inputs on neuron $i$ is above zero, the neuron is active (“fires”),
otherwise, it remains silent.

###  III.1 The structural connectivity

We determined the Hopfield Brain connectivity using Schaefer’s cerebral
cortical parcellation atlas. This parcellation was obtained from the analysis
of a large dataset with a sample size of 1489 Schaefer2018 . This dataset is
publicly available, we used Schaefer’s 1000 parcellation with parcel centroids
mapped onto the MNI152 volumetric space222Github repository dataset url
Schaefer2018 :
https://github.com/ThomasYeoLab/CBIG/tree/master/stable_projects/brain_parcellation/Schaefer2018_LocalGlobal/Parcellations/MNI/Centroid_coordinates.
Fig. 1A shows the spacial distribution of the centroids in the MNI152
volumetric space. From the centroids of these areas, we computed the pairwise
distances. Fig. 1B shows the histogram of the pair distances $d_{ij}$ between
each couple of centroids. Figures 1C and D show respectively the matrix
composed of the Euclidean distances $d_{ij}$, the matrix obtained with the
threshold $J_{ij}>0.4$ were $J_{ij}=e^{-d_{ij}/\delta}$ with $\delta$=0.55.

##  IV The Structure Factor $S_{2}(d)$

Figure 2:  Structure factor $S_{2}(d)$ for $N_{r}=40$ realizations
respectively at $\delta=5.55$ mm (A), and $\delta=5.88$ mm (B). The values
$S_{2}(d)$ were binned over equal sized intervals. The dashed lines indicate
different $\alpha$ values: $\alpha=2/3$ which corresponds to turbulence,
$\alpha=1/2$, the value obtained in Deco2020 , and $\alpha=2/5$. We estimate
the scaling exponent $\alpha$ with a linear regression of
$\log(S(d))\sim\alpha\log(d)+\beta$ in the interval range $[2.7-33.1]$ mm, as
used in Deco2020 , and we obtain $\alpha\approx 2/5$ for $\delta=5.55$ mm, and
$\alpha\approx 2/3$ for $\delta=5.88$ mm.

We used Deco2020  adaptation of Kolmogorov’s concept of structure functions.
In turbulence, the Kolmogorov’s structure factor applies to the longitudinal
velocity but in our case it represents the activity of the parcel which as
measured through the BOLD signal in a fMRI experiment.

To compute the structure factor we use the auxiliary function $B(d)$ defined
as:

$B(d)=\frac{1}{T}\sum_{t=1}^{T}\frac{1}{{\cal{N}}(d)}\sum_{d_{ij}\in
d}s_{i}(t)s_{j}(t)$ (8)

where the function ${\cal{N}}(d)$ represents the number of $i,j$ pairs at
distance $d$. The notation $d_{ij}\in d$ indicates that the distance between
the nodes must be “close enough” to $d$. To this aim, it should be kept in
mind that in Schaefer’s cerebral cortical parcellation atlas, $d$ was recorded
with $2mm$ precision, thus there are multiple $i,j$ pairs with the same $d$.
Therefore we can bin the existing nodes’ distance and map them to the natural
numbers, so as to use of the Kronecker delta $\Delta(ij|d)$, which takes value
$1$ if $d_{ij}$ is in the bin identified by $d$, and is $0$ otherwise.

As a result, we have:

$B(d)=\frac{1}{T}\sum_{t=1}^{T}\frac{1}{{\cal{N}}(d)}\sum_{ij}s_{i}(t)s_{j}(t)\Delta(ij|d)$
(9)

and

${\cal{N}}(d)=\sum_{ij}\Delta(ij|d)$ (10)

The function $B(d)$ relates to the structure factor as follows:

$\displaystyle S_{2}(d)$ $\displaystyle\doteq$
$\displaystyle\frac{1}{T}\sum_{t=1}^{T}\frac{1}{{\cal{N}}(d)}\sum_{ij}(s_{i}(t)-s_{j}(t))^{2}\Delta(ij|d)$
(11) $\displaystyle=$ $\displaystyle 2\left[B(0)-B(d)\right]$

The Hopfield dyamics is deterministic, and since the phase space of the nodes’
states is finite, the asymptotic trajectories must necessarily be periodic.
Therefore, as also discussed in Folli2017 , Hopfield Recurrent Neural
Networks, given an initial state, consistently ends up into a specific steady
state (“attractor”). The steady states are deterministic attractors composed
of either limit cycles formed by repeating sequence of states, or a single
stationary state. Under the parametrizations considered here, we numerically
found that the steady states were always composed of a single fixed point, as
expected in the case of symmetric couplings $J_{ij}$. In the absence of a
further external input (“external stimuli”), the network (the brain) is stuck
in the attractor and the whole activity is constant in time, hence no time
summation was taken in the previous equations. This is of course a very crude
representation of the brain dynamics, which receives continuously “internal”
and “external” inputs, but we can think at the whole set of attractors as a
reasonable representation of the dynamical correlations within the whole brain
model.

Given different values of the decay length $\delta$, we have run the Hopfied
model for a system composed by $N=1000$ neurons (with pairwise distances
consistent with the Schaefer’s cerebral cortical parcellation atlas
Schaefer2018 ). The structure functions have then been calculated by following
the system until the fixed points were reached, and by averaging over
$N_{r}=40$ randomly chosen initial states. Figure 2 shows the structure factor
$S_{2}(d)$ for respectively $\delta=5.55$ mm (Fig. 2A) and $\delta=5.88$ mm
(Fig. 2B). The black points with the error bars represent the avarege $S(d)$
computed over 40 different realizations. The dashed green line corresponds to
$\alpha=2/3$, the orange line corresponds to $\alpha=1/2$, and the blue line
corresponds to $\alpha=2/5$. The $\alpha=1/2$ value (green line) is the value
estimated in Deco2020 . The $\alpha=2/3$ is the exponent associated with a
turbulece. Due to the significant level of statistical fluctuations, we have
estimated the average exponent in two different ways, first as the correct
ensemble average $<S_{2}(d)>\sim d^{\alpha}$ and also as an arithmetic average
of the single realization exponents,
$\overline{\alpha}=N_{r}^{-1}\sum_{r=1}^{N_{r}}\alpha_{r}$. Even though the
latter is not correct, the fact that both yield very similar values attests to
the robustness of the average exponent, notwithstanding sizeable fluctuations.
We estimate the scaling exponent $\alpha$ in the interval $[2.7-33.1]$ mm as
in Deco2020  and we obtained $\alpha\approx 2/5$ for $\delta=5.55$ mm, and
$\alpha\approx 2/3$ for $\delta=5.88$ mm.

###  IV.1 Scaling exponents as a function of the decay length

Figure 3:  A) The blue line shows the average slope of $\alpha=\log S(d)/logd$
for different $\delta$s, given the connectivity matrix
$J_{ij}=\exp{(-d_{i,j}/\delta)}$. The horizontal dotted lines correspond to
$\alpha_{RW}=1$ and $\alpha_{T}=2/3$, corresponding to random walk and
homogeneous incompressible turbulence, respectively, while the dot-dashed
horizontal line corresponds to Deco _et al._ “brain” value $\alpha_{D}=1/2$.
The orange line shows the function $\alpha(\delta)$ with a connectivity matrix
obtained by randomly shuffling all the weights $J_{ij}$, which results into a
basically uncorrelated signal, indicating the importance of the metric
structure of the weights. The plot stops at $\delta=10$, a region where the
signal is highly correlated ($\alpha\sim 1.4$), but still below the smooth
regime marked by $\alpha=2$. A visual extrapolation indicates that the scaling
regime would remain non-smooth ($\alpha<2$) even in the limit of an ”infinite-
brain”, $\delta\to\infty$), the global brain ($d\sim 25\delta$) being pretty
close to the infinite-brain limit. Note however that in this region $B(d)\sim
B(0)$, so that the numerical calculation of the exponent becomes less and less
accurate (a fully correlated signal with $B(d)=B(0)$ yields formally
$\alpha\to\infty$.) Conversely, a brain below $4$ mm, would exhibit ”chaotic”
behaviour, with scaling exponents close to zero, i.e. no scaling regime. The
bottom panel shows the neuron activity mapped on the Schaefer’s cerebral
cortical parcellation atlas coordinates for (B) $\delta=4$ (“chaotic”), (C)
$\delta=5.99$, and (D) $\delta=6.66$ (random walk), respectively.

We repeated the previously depicted procedure for a sequence of values of the
decay length in the range $[3-10]$ mm and measured the associated scaling
exponents. The values of $\alpha(d)$ for different realizations are shown in
Fig.3.

Several comments are in order. First, we observe a clear increasing trend of
$\alpha$ with $\delta$, meaning that small decay lengths, _i.e._ localised
connectivity, promotes irregular patterns. For instance, below $\delta\sim 4$
mm, the scaling exponent is basically zero, corresponding to complete
randomness. This was checked independently, by running simulations in which
the couplings $J_{ij}$ were randomized though keeping the same distribution of
the elements of the matrix $J$ (flat orange line at the bottom). We were able
to simultaneously randomize the couplings $J_{ij}$, and keep the same value
distribution by randomly shuffling the $J_{ij}$ values. At the opposite end,
with $\delta=10$ mm, we measure $\alpha\sim 1.4$, which is still well below
the smoothness threshold $\alpha_{S}=2$, but above the random walk value
$\alpha_{R}=1$, and much above the Kolmogorov turbulence value
$\alpha_{T}=2/3$, which is in turn larger than $\alpha_{D}\sim 1/2$, the value
obtained by Deco _et al._ at $\delta\sim 5.55$ mm. Importantly, in the
intermediate regime, around the physiological value $\delta=5.55$ mm, the
present Hopfield model delivers $\alpha_{H}\sim 2/5$, different but still
close to Deco _et al._ value.

Figure 4 shows how the pair correlation for different values of $\delta$. For
$\delta<5$ mm, the pair correlation shows a basically uncorrelated activity,
except at short range. For $\delta=5.26$ mm, long-range pair correlations
start to emerge and finally, for larger $\delta$, the long-range pair
correlations become even stronger. In the limit of very large $\delta$,
stationary states emerge in which the system freezes into a single fully-
ordered state, with all nodes states either $+1$ or $-1$.

Our analysis confirms that there is indeed a scaling regime in the
organization of the brain patterns which reminds fluid turbulence, yet with a
higher degree of spikiness. Indeed, by sticking to the physiological value
$\delta=5.55$ mm, the Hopfield brain appears more irregular (lower scaling
exponent) than a random walk and also of homogeneous incompressible
turbulence. The most important feature, though, is the sigmoidal dependence of
the relation $\alpha=\alpha(\delta)$ with the rampup region between about $4$
and $7$ mm, centered around $\delta\sim 5.9$. This is indeed a physiologically
relevant scale of local regions in the brain. Furthermore, this rampup region
unveils the major sensitivity of the scaling exponent to small changes of the
decay length, indicating that small changes of the decay length lead to fairly
different scaling regimes. In particular, for $\delta$ below $4$ mm, the
activity appears to be nearly chaotic, whereas for $\delta$ above $8$ mm, it
becomes kind of “rigid”, i.e. globally correlated.

A word of caution must be spent on the statistical accuracy of the scaling
exponents, which for the present case of $40$ realizations, is around $25$
percent.

Figure 4:  Pair-correlation functions $B(d)$ for different $\delta$s. At small
$\delta$ the signal is largely uncorrelated, while at increasing $\delta$,
correlations start to emerge, although affected by significant amount of
statistical fluctuations.

##  V Thresholding the connections

Figure 5:  The average slope of $\log S(d)$ for different dilution levels
$\alpha(\rho)$, normalized to its threshold-free value $\alpha(\rho)$ at that
given $\delta$. The dilution is measured as the fraction of disconected pairs
$i,j$, $J_{ij}$ elements with value 0. The second $y$-axis on top shows the
cutoff distance in mm for wich at a certain dilution $\rho$ the pairs at a
larger or euql distance are disconected. The differently colored curves
indicate simulations with different $\delta$ values. The graph shows that
within the measurement error, the system is not affected by the dilution of
the edges, up to the removal of more than the $95\%$ of the connections, in
descending order of distance. This shows that, in view of the exponential
connectivity $J_{ij}=e^{-d_{ij}/\delta}$, the onset of collective patterns is
sustained mostly by the close connections, within about $4\delta$.

Having assessed that delocalization (large $\delta$s) promotes a correlated
response, it is natural to inquire the role of short versus long-range
couplings in promoting the correlated patterns sustaining the scaling regime.

To this purpose, we performed a series of simulations by progressively
removing pairs of nodes whose coupling strengths lie below a given threshold,
namely $J_{ij}<J_{th}$, where the coupling strengths vary in the range
$[0,1]$. Since $J_{ij}$ decays exponentially with the pair distance, it is
clear that the running threshold removes first the distant nodes and then the
less distant ones, upon increasing the threshold, thereby implementing a
“correlated dilution” process. More precisely, a threshold $J_{th}$ excludes
pairs beyond a distance

$d_{th}=\delta ln{(1/J_{th})},$ (12)

thereby retaining only a fraction $1-J_{th}$ of the full set of interacting
nodes. The result of this progressive dilution on the scaling exponent are
reported in Figure 5. Where the dilution $\rho$ is defined as the fraction of
disconnected nodes $i,j$, $J_{ij}=0$, and the resolt on the scaling exponent
is masured comparing the undiluted exponent $\alpha(0)$ with the exponent
$\alpha(\rho)$ at that dilution. Since the dilution is obtaned thresholding in
order the $J_{ij}$ values untill we get a certain fraction of disconnected
nodes, and since eq. 12 we get that to a certain dilution corresponds a
certain cutoff distance above which all couples are disconnected.

From this figure, a smooth trend at increasing $\rho$ is observed, with the
scaling exponent starting to display significant decrease above
$\rho=\left[0.95,0.98\right]$, when just 2 or 5 percent of the nodes are left,
corresponding to $d_{ij}>4.6\delta$. With $J_{th}=0.10$, corresponding to
$d_{ij}>-\delta ln(0.1)\sim 0.46\delta$, the scaling exponent is basically
halved, indicating a major loss of correlation. This is expected since
$J_{th}=0.10$ cuts out all but shortest range interactions.

The conclusion of this dilution analysis is that the brain response is carried
almost entirely by interactions up to about five decay lengths, namely about
$1/5$ of the size of the global brain. This sounds reasonable, as it strikes a
plausible compromise between short-range ($d<4\delta$) and long-range
($d>6\delta$) interactions. It is interesting to observe that a “solid”
(strongly correlated) brain, such as the one that results for $\delta>10$ mm,
would be too homogeneous to be able to store and process the amount of
information required to function properly.

On the opposite side, for too short connectivity length, the brain would
behave like a “gas” of disconnected neurons, hence incapable of collective
behavior, which is key to its proper functioning

It appears like the physiological length is achieving an optimal compromise
between these two opposite, order-disorder, trends. With a daring but
conducive metaphor, we could speculate that the brain works in a sort of
turbulent liquid-like state, although more spiky than actual turbulent fluids.

##  VI Conclusions

It is shown that an Hopfield recurrent neural network, informed by
experimentally derived brain topology, recovers the scaling picture recently
proposed by Deco Deco2020 , according to which the process of information
exchange within the human brain shows spatially correlated patterns
qualitatively similar to those displayed by turbulent flows.

Our analysis confirms the initial finding by Deco2020  in near-quantitative
form, predicting very similar scaling exponent for the same value of the
connectivity decay length. This is pretty remarkable, given the very different
nature of the two mathematical models and provides a strong hint at the
robustness of the scaling picture.

It is further observed that the scaling exponents are smaller than for the
case of turbulence, indicating that the collective activity of the brain is
more irregular (spiky) than for a homogeneous incompressible turbulent fluid.

Importantly, the scaling exponents show a steep dependence on the spatial
range of neural connections, hence the specific value of the decay length has
a major impact on the degree of irregularity (spikiness) of the collective
patterns.

This provides a strong indication that there is an optimal connectivity length
at which the brain shows non-smooth functional behavior. This appears to
happen at a scale around $5,5$ mm, hence well below the global size of the
brain Longer/shorter connectivities would lead to excessive
smoothness/roughness respectively. Interestingly, a dilution analysis shows
that the Hopfield model brain remains functional by removing links above about
six decay lengths, namely $25-30$ mm, corresponding to about one fifth of the
size of the global brain. This indicates that, in terms of connectivity decay
length, the brain seems to function in a sort of intermediate ”turbulent
liquid-like” state , whose essential connections lie in near geometrical mean
between the decay length and the global size of the brain. This turbulent
liquid-like state appears to be more spiky than ordinary homogeneous
incompressible turbulent fluids.

##  VII Acknowledgments

One of the authors (SS) wishes to acknowledge financial support from the
Bertarelli Foundation during his stay at the Neurobiology Department of the
Harvard Medical School. The authors wish to acknowledge illuminating
discussions with John A. Assad, Haim Sompolinsky, and Bernardo Sabattini. GG
contributes to code writing and simulations. GG, SS, and GR conceived the
Idea, designed the experiments, interpreted the data, and wrote the
manuscript.

##  VIII Bibliography

## References

  * (1) G. Deco and M. L. Kringelbach, “Turbulent-like dynamics in the human brain,” Cell Reports, vol. 33, p. 108471, 12 2020.
  * (2) S. I. Amari, “Learning patterns and pattern sequences by self-organizing nets of threshold elements,” IEEE Transactions on Computers, vol. C-21, pp. 1197–1206, 1972.
  * (3) W. A. Little, “The existence of persistent states in the brain,” Mathematical Biosciences, vol. 19, pp. 101–120, 2 1974.
  * (4) J. J. Hopfield, “Neural networks and physical systems with emergent collective computational abilities.,” Proceedings of the National Academy of Sciences of the United States of America, vol. 79, pp. 2554–2558, apr 1982.
  * (5) D. J. Amit, H. Gutfreund, and H. Sompolinsky, “Storing Infinite Numbers of Patterns in a Spin-Glass Model of Neural Networks,” Physical Review Letters, vol. 55, pp. 1530–1533, sep 1985.
  * (6) N. Brunel, “Is cortical connectivity optimized for storing information?,” Nature Neuroscience, vol. 19, pp. 749–755, may 2016.
  * (7) C. J. Hillar and N. M. Tran, “Robust Exponential Memory in Hopfield Networks,” Journal of Mathematical Neuroscience, vol. 8, no. 1, 2018.
  * (8) C. Hillar, T. Chan, R. Taubman, and D. Rolnick, “Hidden Hypergraphs, Error-Correcting Codes, and Critical Learning in Hopfield Networks,” Entropy 2021, Vol. 23, Page 1494, vol. 23, p. 1494, nov 2021.
  * (9) D. Stauffer, A. Aharony, L. D. F. Costa, and J. Adler, “Efficient hopfield pattern recognition on a scale-free neural network,” The European Physical Journal B - Condensed Matter and Complex Systems 2003 32:3, vol. 32, pp. 395–399, 2003.
  * (10) D.-H. Kim, J. Park, and B. Kahng, “Enhanced storage capacity with errors in scale-free hopfield neural networks: An analytical study,” PLOS ONE, vol. 12, p. e0184683, 10 2017.
  * (11) V. Folli, G. Gosti, M. Leonetti, and G. Ruocco, “Effect of dilution in asymmetric recurrent neural networks,” Neural Networks, vol. 104, pp. 50–59, 8 2018.
  * (12) M. Leonetti, V. Folli, E. Milanetti, G. Ruocco, and G. Gosti, “Network dilution and asymmetry in an efficient brain,” Philosophical Magazine, vol. 100, pp. 2544–2555, oct 2020.
  * (13) K. Gopalsamy and X. zhong He, “Stability in asymmetric hopfield nets with transmission delays,” Physica D: Nonlinear Phenomena, vol. 76, pp. 344–358, 9 1994.
  * (14) Z.-B. Xu, G.-Q. Hu, and C.-P. Kwong, “Asymmetric hopfield-type networks: Theory and applications,” Neural Networks, vol. 9, pp. 483–501, 4 1996\.
  * (15) T. Chen and S. I. Amari, “Stability of asymmetric hopfield networks,” IEEE Transactions on Neural Networks, vol. 12, pp. 159–163, 1 2001.
  * (16) F. M. Franca and Z. Yang, “Building artificial cpgs with asymmetric hopfield networks,” Proceedings of the International Joint Conference on Neural Networks, vol. 4, pp. 290–295, 2000.
  * (17) P. Zheng, J. Zhang, and W. Tang, “Analysis and design of asymmetric hopfield networks with discrete-time dynamics,” Biological Cybernetics, vol. 103, pp. 79–85, 7 2010.
  * (18) A. Szedlak, G. Paternostro, and C. Piermarocchi, “Control of asymmetric hopfield networks and application to cancer attractors,” PLOS ONE, vol. 9, p. e105842, 8 2014.
  * (19) G. Gosti, E. Milanetti, V. Folli, F. de Pasquale, M. Leonetti, M. Corbetta, G. Ruocco, and S. Della Penna, “A recurrent hopfield network for estimating meso-scale effective connectivity in meg,” Neural Networks, vol. 170, pp. 72–93, 2024.
  * (20) D. O. Hebb, The Organization of Behavior; A Neuropsychological Theory.  Wiley, 10 1949.
  * (21) A. Fachechi, E. Agliari, and A. Barra, “Dreaming neural networks: Forgetting spurious memories and reinforcing pure ones,” Neural Networks, vol. 112, 2019.
  * (22) M. Benedetti, E. Ventura, E. Marinari, G. Ruocco, and F. Zamponi, “Supervised perceptron learning vs unsupervised hebbian unlearning: Approaching optimal memory retrieval in hopfield-like networks,” Journal of Chemical Physics, vol. 156, 2022.
  * (23) M. Benedetti and E. Ventura, “Training neural networks with structured noise improves classification and generalization,” arXiv, vol. 2302.13417v4, 2 2023.
  * (24) G. Gosti, V. Folli, M. Leonetti, and G. Ruocco, “Beyond the Maximum Storage Capacity Limit in Hopfield Recurrent Neural Networks,” Entropy, vol. 21, p. 726, jul 2019.
  * (25) S. Hwang, V. Folli, E. Lanza, G. Parisi, G. Ruocco, and F. Zamponi, “On the number of limit cycles in asymmetric neural networks,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2019, p. 053402, 5 2019.
  * (26) S. Hwang, E. Lanza, G. Parisi, J. Rocchi, G. Ruocco, and F. Zamponi, “On the number of limit cycles in diluted neural networks,” Journal of Statistical Physics, vol. 181, pp. 2304–2321, 12 2020.
  * (27) U. Frisch, Turbulence.  Cambridge University Press, 11 1995.
  * (28) R. Benzi, S. Ciliberto, R. Tripiccione, C. Baudet, F. Massaioli, and S. Succi, “Extended self-similarity in turbulent flows,” Physical Review E, vol. 48, p. R29, 7 1993.
  * (29) A. Schaefer, R. Kong, E. M. Gordon, T. O. Laumann, X.-N. Zuo, A. J. Holmes, S. B. Eickhoff, and B. T. T. Yeo, “Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity mri,” Cerebral Cortex, vol. 28, pp. 3095–3114, 9 2018.
  * (30) V. Folli, M. Leonetti, and G. Ruocco, “On the Maximum Storage Capacity of the Hopfield Model,” Frontiers in Computational Neuroscience, vol. 10, p. 144, jan 2017.


# Physics of Smart Matter: integrating active matter and control to gain
insights into living systems

Herbert Levine1 and Daniel I. Goldman2 1 Departments of Physics and
Bioengineering, Northeastern University, Boston, MA, USA  2 School of Physics,
Georgia Institute of Technology, Atlanta, Georgia 30332, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We offer our opinion on the benefits of integration of insights from active
matter physics with principles of regulatory interactions and control to
develop a field we term “smart matter”. This field can provide insight into
important principles in living systems as well as aid engineering of
responsive, robust and functional collectives.

††preprint: APS/123-QED

Recent years have seen the dawning of a new vibrant subfield of physics, the
physics of living systems. An excellent discussion of the history and promise
of this research area was presented in a recent National Academy report [1].
The research agenda of this field encompasses systems ranging from
biomolecules to ecosystems, from building models to building robots [2], all
in service of a quantitative understanding of systems whose behavior
transcends what we have come to expect from experience with inanimate matter.
Here we offer our opinion as to what might be for a physical scientist some of
the necessary ingredients to consider a system to be living. We will argue
that a path toward addressing such questions and characterizing these systems
will be in developing ideas of “smart matter” (which we will define below)
that combine ideas of the exploding field of active matter with those of the
less appreciated (by physicists) regulatory interactions (e.g. feedback
control [3, 4]).

Our starting point is the idea of active matter [5], a field that many
physicists have come to believe underlies the secrets of life. Active matter
refers to physical organizations of interacting constituents that each have
their own access to energy sources. These constituents can be living, as is
the case in bacterial colonies [6], ant rafts [7] or bird flocks [8],
completely abiotic as in colloids propelled by catalyzed chemical reactions or
motor driven robots [9, 10], or “in-between” as in the beautiful dynamical
structures created in vitro by biopolymers activated by molecular motors [11]
or even biohybrid robots composed of soft materials and living cells [12]. The
study of active matter in the physics community took off with the seminal work
of Ben-Jacob, Vicsek and collaborators [13] who showed that these systems can
self-organize in ways that circumvent many of the restrictions exhibited by
“normal matter”, operating close to equilibrium. Clearly, any living system is
active, using stored energy and functioning far away from any thermal
equilibrium state.

So, is every active system alive? Clearly not. But, what does it take to go
from active matter to a living system? The first step involves the
predominance of regulatory interactions [14, 15] . A comparison of laboratory
preparations involving active biopolymers [16] and the actual situation that
prevails in the cytoskeleton of a living cell makes the point. In the latter,
there are many dozens of proteins that regulate all aspects of the polymer
chemistry and couple reactions to cellular conditions. Thereby, actin
polymerization is restricted to the front of a moving cell [17], microtubules
attach in a highly controlled manner to segregating chromosomes [18], and
intermediate filaments such as Vimentin arrange themselves geometrically to
cushion against nuclear deformation [19]. Active matter physics is slaved to
needed functionality, being necessary but not sufficient. Coupling active
matter to controlling regulators leads to a new and often qualitatively
different class of objects that we will refer to as “smart matter” [20, 21].

The operating principles behind regulatory interactions concern information
flow, by which we mean that details regarding the external world are used to
modulate active behavior. As the relevant environmental inputs vary, the
active system responds by realigning its dynamics accordingly; the more
complex and often, the more energy-intensive, a regulatory system is, the
better it can create useful correlations between the environment and active
matter behavior. The active matter can now behave intelligently. And, living
systems are clearly very smart in ways we are only beginning to understand.
Biology over the past half century has revealed an astounding complexity in
the control of all important processes. One can assume with confidence that
every step of a biological process will be regulated in multiple ways and over
multiple timescales. One can also assume with confidence that individual
components of the underlying active matter will evolve to become more flexible
in their ability to respond to regulatory input - a case in point is the
fascinating story of the evolution of the mammalian synapse [22], where the
sophistication of the molecular machinery has grown even as the neural systems
it serves have become larger. Perhaps these facts will eventually validate the
musings of Schrödinger in his famous “What is Life” text [23] on the need for
coming up with new concepts of physics to accommodate the workings of smart
matter, workings that just do not occur naturally in the abiotic world.

Of course, the boundary between active matter versus smart matter can
sometimes be rather fuzzy. Let us take for example the automobile traffic in a
big city freeway system. Until quite recently, this could serve as an obvious
example of an active matter system that had no overarching regulatory dynamics
shaping the local interactions of the active motorists. Even though individual
motorists clearly want to minimize their travel time and avoid collisions,
there was no goal driving the dynamics of the system as a whole. But, an
argument can be made that the advent of tools such as Waze and Google Maps has
indeed provided the regulatory feedback missing in the active matter paradigm.
Now, motorists do modulate their interactions and decisions based on freeway
conditions on a variety of length scales, and the network as a whole does
attempt to optimize transportation functionality. In fact, the city of Boston
has an ongoing partnership with Waze in which data is fed to the city’s
traffic management center in order to adjust traffic signals, explicitly meant
to optimize transportation efficiency. As surprising as that may seem to those
of us who live there, Boston traffic may be becoming intelligent.
Parenthetically, this example serves to emphasize that we do not restrict the
use of the word ”matter” to tangible physical ”stuff”; matter to us just means
a substrate upon which the actions of interactive components takes place. This
generalization might become particularly important in the future, when we are
forced to confront the question of living artificial intelligences.

Should we consider any smart matter system to be alive? The question of what
it means to be alive is fraught with millennia of philosophical debate that is
hard to place in a physical science context. Nonetheless, living systems might
require something beyond being simply smart matter possessing efficient
information flow governing active medium response. Let us turn to another
example in which smart management of environmental interactions is critical
for function. The engineering community has devoted considerable effort to
building multi-legged robots that can navigate effectively over difficult
terrain [24, 25, 26]. To do this, the robots are typically equipped with a
complex suite of sensors that provide input for controlling sophisticated
actuators and/or leveraging embodied (aka mechanical) ”intelligence” [27] in a
by now familiar smart matter pattern. However, despite the impressively life-
like agility and performance increasingly being displayed by such devices, it
is unlikely that many investigators would consider such a robot to be alive,
in the same sense that an insect navigating the same terrain would be. Why? We
can extend the question by imagining that we endow the robot with a battery
sensor such that when it detects that its power level is getting low, it stops
what it is doing and goes off in search of ”electric food” from the nearest
charging station. Is the robot now closer to being a living system? What if
there is in addition a sensor detecting potentially hazardous weather
conditions and the robot “knows” to seek shelter. What if the robot can
decide, based on the rate of progress it is making on an assigned task, that
it needs more copies of itself and can arrange to have that happen by ordering
from a factory with which it is in contact. Thus, the question posed to
physics of living systems researchers is whether the difference between living
and smart is just one of degree of systems integration and semantics, or
alternatively involves a true phase transition leading to new capabilities in
a discontinuous manner. There is no real hint at present as to what might
cause such a transition, or if the transition/bifurcation concept is even
relevant. The same question arises of course in the microscopic realm,
concerning the ancient origin of life and modern attempts to artificially
synthesize living cells and to characterize what a minimum cell must consist
of [28]. And, this is becoming ever more critical, as advances in astronomy,
such as the exoplanet revolution, have brought to the fore questions of how
best to search for indications of life elsewhere in the cosmos.

We note that even without discussing if such systems are “alive” we can
utilize smart matter systems as models to discover principles by which living
systems achieve robust function [29]. Indeed, physicists have a long history
of being interested in deep questions but not letting these get in the way of
making tangible progress; any history of quantum theory will clearly attest to
this useful duality. So, there is much work to done in figuring out how to
best couple active matter to smart controllers to enable the accomplishment of
various tasks. In this regard, we will need to work directly on all manner of
living systems and with all manner of biologists. In some parts of this
endeavor, there is a need for mutual re-education. While physicists are most
comfortable with the active matter paradigm, modern biological research often
stresses the regulatory aspects of living systems at the expense of working
back down to physical processes that interact with the environment. Many
papers focus almost exclusively on gene expression (aka transcriptomics) and
protein abundance as the ultimate in defining cell states and cell physiology.
This perspective is becoming even more entrenched as spatial transcriptomics
[30] and technologies such as tissue CyTOF [31] begin to define tissues and
organs solely in terms of omics profiles. It imagines real-world action as
something which can automatically flow from the information processing level -
there’s always “a gene for that”. An example of the problem arose quite a long
time ago in a paper [32] claiming to have found the gene responsible for
creating the chemical wave field responsible for guiding the motion of
hundreds of thousands of Dictyostelium amoebae towards aggregation centers, as
part of their survival strategy in the face of starvation. A gene and its
single protein product cannot make a millimeter scale chemical pattern;
instead a gene can help control a physiochemical system of interacting
components that are capable of doing the necessary spadework.

We therefore posit that the physics of living systems can benefit from a
synergistic merging of these two insufficient worlds views; creation of the
field of “smart matter” can provide researchers a way to frame such an
integration and develop new models of living systems. In our opinion, this
merger is essential. Ignoring the constraints placed on living systems by the
need to get the molecules, cells, tissues and organs to actually accomplish
the needed tasks will miss essential constraints on behavior. Assuming that
active matter systems are all we need to focus on as we move forward dismisses
out of hand many of the performance aspects of systems that allow us to
consider them living. Molecular and cellular biologists should realize genes
are not magic wands that can wish physical effects into existence. Active and
soft matter physicists need to take to heart a quote from Alan Turing
regarding patterning of the zebra, that “..the stripes are easy, but what
about the horse part?”.

###### Acknowledgements.

H.L. thanks the NSF Physics Frontier Center program for support. D.I.G. thanks
numerous colleagues for helpful discussion over the years and the NSF Physics
of Living Systems program, the Army Research Office MURI on algorithmic
matter, and the Dunn Family Professorship.

## References

  * National Academies of Sciences  _et al._ [2022] E. National Academies of Sciences, Medicine, _et al._ ,  _Physics of Life_ (National Academies Press, 2022).
  * Aguilar  _et al._ [2016] J. Aguilar, T. Zhang, F. Qian, M. Kingsbury, B. McInroe, N. Mazouchova, C. Li, R. Maladen, C. Gong, M. Travers, _et al._ , A review on locomotion robophysics: the study of movement at the intersection of robotics, soft matter and dynamical systems, Reports on Progress in Physics 79, 110001 (2016).
  * Åström and Murray [2021] K. J. Åström and R. M. Murray,  _Feedback systems: an introduction for scientists and engineers_ (Princeton university press, 2021).
  * Cowan  _et al._ [2014] N. J. Cowan, M. M. Ankarali, J. P. Dyhr, M. S. Madhav, E. Roth, S. Sefati, S. Sponberg, S. A. Stamper, E. S. Fortune, and T. L. Daniel, Feedback control as a framework for understanding tradeoffs in biology, American Zoologist 54, 223 (2014).
  * Marchetti  _et al._ [2013] M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Hydrodynamics of soft active matter, Reviews of modern physics 85, 1143 (2013).
  * Ben-Jacob  _et al._ [2000] E. Ben-Jacob, I. Cohen, and H. Levine, Cooperative self-organization of microorganisms, Advances in Physics 49, 395 (2000).
  * Mlot  _et al._ [2011] N. J. Mlot, C. A. Tovey, and D. L. Hu, Fire ants self-assemble into waterproof rafts to survive floods, Proceedings of the National Academy of Sciences 108, 7669 (2011).
  * Cavagna and Giardina [2014] A. Cavagna and I. Giardina, Bird flocks as condensed matter, Annu. Rev. Condens. Matter Phys. 5, 183 (2014).
  * Li  _et al._ [2021] S. Li, B. Dutta, S. Cannon, J. J. Daymude, R. Avinery, E. Aydin, A. W. Richa, D. I. Goldman, and D. Randall, Programming active cohesive granular matter with mechanically induced phase changes, Science Advances 7, eabe8494 (2021).
  * Aguilar  _et al._ [2018] J. Aguilar, D. Monaenkova, V. Linevich, W. Savoie, B. Dutta, H.-S. Kuan, M. Betterton, M. Goodisman, and D. Goldman, Collective clog control: Optimizing traffic flow in confined biological and robophysical excavation, Science 361, 672 (2018).
  * Ndlec  _et al._ [1997] F. Ndlec, T. Surrey, A. C. Maggs, and S. Leibler, Self-organization of microtubules and motors, Nature 389, 305 (1997).
  * Filippi  _et al._ [2022] M. Filippi, O. Yasa, R. D. Kamm, R. Raman, and R. K. Katzschmann, Will microfluidics enable functionally integrated biohybrid robots?, Proceedings of the National Academy of Sciences 119, e2200741119 (2022).
  * Vicsek  _et al._ [1995] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel type of phase transition in a system of self-driven particles, Physical review letters 75, 1226 (1995).
  * Moat  _et al._ [2002] A. G. Moat, J. W. Foster, and M. P. Spector,  _Microbial physiology_ (John Wiley & Sons, 2002).
  * Schmidt-Nielsen [1997] K. Schmidt-Nielsen,  _Animal physiology: adaptation and environment_ (Cambridge university press, 1997).
  * Koenderink  _et al._ [2009] G. H. Koenderink, Z. Dogic, F. Nakamura, P. M. Bendix, F. C. MacKintosh, J. H. Hartwig, T. P. Stossel, and D. A. Weitz, An active biopolymer network controlled by molecular motors, Proceedings of the National Academy of Sciences 106, 15192 (2009).
  * Weiner [2002] O. D. Weiner, Regulation of cell polarity during eukaryotic chemotaxis: the chemotactic compass, Current opinion in cell biology 14, 196 (2002).
  * Hawkins  _et al._ [2010] T. Hawkins, M. Mirigian, M. S. Yasar, and J. L. Ross, Mechanics of microtubules, Journal of biomechanics 43, 23 (2010).
  * Pogoda  _et al._ [2022] K. Pogoda, F. Byfield, P. Deptuła, M. Cieśluk, Ł. Suprewicz, K. Skłodowski, J. L. Shivers, A. van Oosten, K. Cruz, E. Tarasovetc, _et al._ , Unique role of vimentin networks in compression stiffening of cells and protection of nuclei from compressive stress, Nano Letters (2022).
  * Savoie  _et al._ [2019] W. Savoie, T. A. Berrueta, Z. Jackson, A. Pervan, R. Warkentin, S. Li, T. D. Murphey, K. Wiesenfeld, and D. I. Goldman, A robot made of robots: Emergent transport and control of a smarticle ensemble, Science Robotics 4, eaax4316 (2019).
  * Ozkan-Aydin  _et al._ [2021] Y. Ozkan-Aydin, D. I. Goldman, and M. S. Bhamla, Collective dynamics in entangled worm and robot blobs, Proceedings of the National Academy of Sciences 118, e2010542118 (2021).
  * Ryan and Grant [2009] T. J. Ryan and S. G. Grant, The origin and evolution of synapses, Nature Reviews Neuroscience 10, 701 (2009).
  * Schrödinger  _et al._ [1992] E. Schrödinger _et al._ ,  _What is life?: With mind and matter and autobiographical sketches_ (Cambridge University Press, 1992).
  * Guizzo [2019] E. Guizzo, By leaps and bounds: An exclusive look at how boston dynamics is redefining robot agility, IEEE Spectrum 56, 34 (2019).
  * Saranli  _et al._ [2001] U. Saranli, M. Buehler, and D. E. Koditschek, Rhex: A simple and highly mobile hexapod robot, The International Journal of Robotics Research 20, 616 (2001).
  * Raibert [1986] M. H. Raibert,  _Legged robots that balance_ (MIT press, 1986).
  * Pfeifer and Bongard [2006] R. Pfeifer and J. Bongard,  _How the body shapes the way we think: a new view of intelligence_ (MIT press, 2006).
  * Glass  _et al._ [2017] J. I. Glass, C. Merryman, K. S. Wise, C. A. Hutchison, and H. O. Smith, Minimal cells—real and imagined, Cold Spring Harbor Perspectives in Biology 9, a023861 (2017).
  * Aina  _et al._ [2022] K. O. Aina, R. Avinery, H.-S. Kuan, M. D. Betterton, M. A. Goodisman, and D. I. Goldman, Toward task capable active matter: Learning to avoid clogging in confined collectives via collisions, Frontiers in Physics , 467 (2022).
  * Burgess [2019] D. J. Burgess, Spatial transcriptomics coming of age, Nature Reviews Genetics 20, 317 (2019).
  * Bodenmiller [2016] B. Bodenmiller, Multiplexed epitope-based tissue imaging for discovery and healthcare applications, Cell systems 2, 225 (2016).
  * Pálsson  _et al._ [1997] E. Pálsson, K. J. Lee, R. E. Goldstein, J. Franke, R. H. Kessin, and E. C. Cox, Selection for spiral waves in the social amoebae dictyostelium, Proceedings of the National Academy of Sciences 94, 13719 (1997).