# Harnessing Synthetic Active Particles for Physical Reservoir Computing

Xiangzun Wang1,2 and Frank Cichos1∗
1Peter Debye Institute for Soft Matter Physics, Leipzig University, 04103
Leipzig, Germany
2Center for Scalable Data Analytics and Artificial Intelligence (ScaDS.AI)
Dresden/Leipzig, 04105 Leipzig, Germany
∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>leipzig.de

###### Abstract

The processing of information is an indispensable property of living systems
realized by networks of active processes with enormous complexity. They have
inspired many variants of modern machine learning one of them being reservoir
computing, in which stimulating a network of nodes with fading memory enables
computations and complex predictions. Reservoirs are implemented on computer
hardware, but also on unconventional physical substrates such as mechanical
oscillators, spins, or bacteria often summarized as physical reservoir
computing. Here we demonstrate physical reservoir computing with a synthetic
active microparticle system that self-organizes from an active and passive
component into inherently noisy nonlinear dynamical units. The self-
organization and dynamical response of the unit is the result of a delayed
propulsion of the microswimmer to a passive target. A reservoir of such units
with a self-coupling via the delayed response can perform predictive tasks
despite the strong noise resulting from Brownian motion of the microswimmers.
To achieve efficient noise suppression, we introduce a special architecture
that uses historical reservoir states for output. Our results pave the way for
the study of information processing in synthetic self-organized active
particle systems.

## Introduction

Storing and processing of information is vital for living systems Tkacik.2014
. The detection of low amounts of chemicals by a bacterium to navigate
environments Wadhams2004 ; Celani.2010 , the feedback mechanisms controlling
and maintaining the function of organisms Cosentino.2019 , or the highly
sophisticated computations in large biological neural networks in the brain
Knudsen.1987  are intricate examples of this importance and created by
evolutionary development. All these processes with living matter as the
substrate of computation rely on its inherent activity, e.g., the microscopic
energy conversion to power the signalling cascades in the presence strong
thermal noise. They have inspired many computational models of machine
learning that are not executed on living matter, but on well-designed
electronic hardware using completely different information representation than
living matter Marković.2020 . Recurrent neural networks are a variant of such
mathematical algorithms with a fading memory that allow learning from
information sequences as in language or time series Abiodun2018 . Reservoir
computers employ sparsely and statically connected recurrent nodes Jaeger2001
; Maass2002 ; Verstraeten2007  or even a single node by using time-
multiplexing Paquot.2010 ; Appeltant2011  to create a high dimensional space.
Information can be injected into this space to spread over the many degrees of
freedom. Unlike the training of other neural networks, where the interactions
of all components are optimized, training reservoir computers is often only
restricted to finding how the desired information can be retrieved from the
node states using adjustable readout only Lukosevicius2009 ; Lukosevicius2012
.

As one of the main properties of recurrent nodes is the memory of past states,
reservoir computers also allow for a physical realization on unconventional
computational substrates Tanaka2019 ; Nakajima2020 ; Rafayelyan.2020  using
optoelectronic oscillators Larger2012 ; Paquot2012 ; VanDerSande2017 ,
mechanical oscillators Coulombe2017 ; Dion2018 , carbon nanotubes Wu.2021  or
passive soft bodies Tanaka2019  as excitable physical systems. It is thus
intruiging to close the loop and draw inspiration from active living systems
to explore microscopic reservoir computing in synthetic active microsystems,
where noise is omnipresent as well but a precise control over the shape and
the physics of the active system is possible. Motile synthetic active
particles have generated enormous interest as a model for self-propelled
systems far from equilbrium and emergent collective effects Bechinger.2016fj9
, that mimic, for example, the dynamics of swarms Lavergne.2019 ;
Baeuerle.2020 ; Wang2023 . Information processing Woodhouse.2017 ; Colen.2021
and learning Cichos2020  in experiment Muinos-Landin2021  and simulation
Liebchen.2019 ; Schneider.2019 ; Colabrese.201760u ; Colabrese.2018 ;
Gustavsson.2017  have also entered the field of synthetic active matter.
However, studies that extend the use of synthetic microparticles as
computational substrates are still rare or purely computational Lymburn2021 .

We demonstrate that motile self-propelled active microparticles can be used
for physical reservoir computing. An active particle self-organizes into a
nonlinear dynamical unit based on a retarded propulsion towards an immobile
target forming a noisy physical recurrent node. The node is perturbed by a
time-multiplexed input signal to form a network of virtual nodes with sparsely
connected topology. Multiple of these active units realize a high dimensional
space of our reservoir computer. Harnessing the physics and inherent dynamics
of the active particles, the reservoir computer is capable of predicting
chaotic signals despite the strong influence of Brownian motion of the active
particles. In particular, we find, that using historical reservoir states for
the output derivation effectively suppresses the intrinsic noise of the
reservoir opening new routes for reservoir computing in noisy systems.

## Results

### Active particle recurrent node

A reservoir computer (RC), as a paradigm derived from recurrent neural
networks, consists of recurrent nodes that nonlinearly process external signal
inputs as well as their previous outputs Jaeger2001 ; Verstraeten2007 . We
realize a simple recurrent node with the help of a single synthetic active
particle as a microscopic model for motile active matter Bechinger.2016fj9 .
The active particle is a polymer microbead of
$2.19\,$\mathrm{\SIUnitSymbolMicro m}$$ diameter with the surface decorated by
gold nanoparticles (about $8\,$\mathrm{nm}$$ diameter). It is immersed in a
thin layer of water bounded by two cover glasses and can move freely in two
dimensions. It is propelled with a speed $v_{0}$ Franzl2021  towards an
immobile target by partial heating of the gold nanoparticles using a focused
laser in a microscopy setup (see Fig. 1A, B and Methods section). The
continuous propulsion in a desired direction is realized by adjusting the
focused laser spot position in real time with the help of a feedback loop
using a spatial light modulator (SLM). A time delay $\delta t$ is added to the
intrinsic feedback latency to control the active particle Khadka2018 ;
Wang2023  and to introduce a finite reaction time, as is inherent in many
processes in living systems. Correspondingly, the attraction
$\hat{\mathbf{F}}(t)$ experienced by the particle at time $t$ is determined by
its position at previous time $t-\delta t$,

$\hat{\mathbf{F}}(t)=-\frac{\mathbf{r}(t-\delta t)}{|\mathbf{r}(t-\delta
t)|},$ (1)

where $\mathbf{r}$ denotes the location of the active particle with respect to
the immobile particle center. The consequence of this retardation is a
particle angular displacement $\theta(t)=\phi(t)-\phi(t-\delta t)$ during the
delay time, which results in a transient rotational motion of the active
particle around the immobile target (Fig. 1C), where $\phi$ denotes the
angular position of the particle. The dynamics of the angular position can
then be described by a nonlinear delay differential equation

$\dot{\phi}(t)=\frac{v_{0}}{R_{0}}\sin\big{(}\theta(t)+u(t)\big{)}+\frac{\sqrt{2D}}{R_{0}}w(t)$
(2)

assuming that the active and the immobile particle (radius of
$a_{\mathrm{act}}$ and $a_{\mathrm{imm}}$) are in physical contact, i.e.
$R_{0}=a_{\mathrm{act}}+a_{\mathrm{imm}}$, due to the delayed attraction. The
active particle in water is subject to Brownian motion as represented by the
noise term in Eq. 2 with $w(t)$ denoting Gaussian white noise. The diffusion
coefficient was determined from the experiment to be $D$ = 0.08 µm2s-1, giving
rise to a Pećlet number of $\mathrm{Pe}$ = $a_{\mathrm{act}}v_{0}/D$ = 38.7.
To make the physical recurrent node capable of receiving external inputs, we
introduce the $u(t)$ in Eq. 2 representing an angular deviation of the
particle propulsion direction from $\hat{\mathbf{F}}$ (Fig. 1C).

Figure 1:  Experimental realization. A Experimental setup (see Sec. Method for
detail). B Active particle recurrent unit consisting of a gold-nanoparticle
covered melamine resin particle (MF-Au) and an immobile target particle (PS).
For active particle propulsion a $532~{}$\mathrm{nm}$$ laser is focused to a
distance d $d$ from the particle center. The resulting heat and asymmetric
temperature induce a self-thermophoretic motion of the particle with a speed
of $v_{0}$ and a direction set the vector from the laser to the particle
center. C Top view of the active particle system. The active particle is
controlled to carry out a motion towards the immobile particle along
$\hat{\mathbf{F}}(t)$ with a time delay $\delta t$ in its reaction. The
direction of $\hat{\mathbf{F}}(t)$ (dashed arrow) is determined by the
previous active particle location $\mathbf{r}(t-\delta t)$. An additional
angle $u(t)$ between the particle propulsion direction (solid arrow) and
$\hat{\mathbf{F}}(t)$ represents an external input into the system. D
Darkfield microscopy image of the sample consisting of 10 active-immobile
particle pairs (larger circle is the immobile particle, smaller circle the
active particle) as physical nodes in the experiment. An additional calibrator
is used as an active particle swimming along a square route to measure the
propulsion speed $v_{0}$. The real-time video is provided in Supplementary
Movie 1.

For its function as a recurrent node, the dynamics of the angle $\theta(t)$ is
important. The dynamics described by Eq. 2 can be approximated by an
overdamped motion of a particle in a self-generated effective quartic
potential $U_{\rm eff}(\theta)$ Wang2023 , which resembles a generic Landau-
type description Goldenfeld.2018 ; Wang2023 . The shape of the potential is
controlled by a dimensionless parameter $\theta_{0}=v_{0}\delta t/R_{0}$. With
increasing $\theta_{0}$, $U_{\rm eff}(\theta)$ transitions from a near
parabolic shape with a minimum at $\theta=0$ to a symmetric double well shape
with minima at $\theta_{+}$ and $\theta_{-}$ in a pitchfork bifurcation (see
Fig. 2C and Sec. S2 of Supplementary Information).

A perturbation of $\theta$ (solid arrow in Fig. 2C), will result in a
relaxation with a dynamics determined by the control parameter $\theta_{0}$.
We have experimentally determined the response $\theta(t)$ to an impulsive
perturbation $u(t)=\delta(t)$ for different values of $\theta_{0}$ (Fig. 2A,
B). The individual experimental trajectories $\theta(t)$ (grey lines in Fig.
2A) strongly fluctuate due to the Brownian motion of the active particle.
However, the ensemble average of 500 trajectories of $\theta(t)$ (red lines in
Fig. 2A and solid lines in B) exhibits an asymptotic behavior, which nicely
reflects the response evaluated in deterministic simulations (Fig. 2B, dashed
lines). The characteristic relaxation time $\tau_{\theta}$ extracted from the
deterministic simulations reveals the expected strong increase around the
transition point. Tuning the control parameter $\theta_{0}$, i.e. activity
($v_{0}$) and/or delay ($\delta t$) therefore allows to manipulate the fading
memory of the active particle recurrent node, which is paramount for the
coupling and response of the recurrent nodes in our reservoir computer. Note,
that the nonlinear dynamics of the active particle system is the result of the
delayed propulsion towards the target.

### Reservoir computer with active particle nodes

The asymptotic relaxation of $\theta(t)$ demonstrated in Fig. 2B represents a
basic requirement for reservoir computing Boyd1985 ; Maass2002 ; Maass2004 ;
Dambre2012 . Based on the time-delay it also allows the setup of coupled
recurrent nodes. In the discrete-time setting of our experiment with the
sampling period $\Delta t$, Eq. 2 can be rewritten as

$\phi(T)=\phi(T-1)+\frac{v_{0}\Delta
t}{R_{0}}\sin\big{(}\theta(T-1)+u(T-1)\big{)}+\frac{\sqrt{2D\Delta
t}}{R_{0}}W(T-1),$ (3)

where $T$ = $t/\Delta t$ is an integer number representing the timestep. $W$
denotes Gaussian random numbers with zero mean and unit variance. The
evolution of $\theta$ then follows as

$\displaystyle\theta(T)=\theta(T-1)$ $\displaystyle+\frac{v_{0}\Delta
t}{R_{0}}\sin\big{(}\theta(T-1)+u(T-1)\big{)}$ (4)
$\displaystyle-\frac{v_{0}\Delta t}{R_{0}}\sin\big{(}\theta(T-\delta
T-1)+u(T-\delta T-1)\big{)}$ $\displaystyle+\frac{\sqrt{2D\Delta
t}}{R_{0}}\big{(}W(T-1)-W(T-\delta T-1)\big{)},$

with the discrete-time delay $\delta T$ = $\delta t/\Delta t$. Refering to the
concept of virtual nodes and time-multiplexing Appeltant2011 , we consider the
transient state $\theta(T)$ of a physical node at different timesteps as
virtual nodes constituting the reservoir. Each virtual node state $\theta(T)$
is, according to Eq. 4, coupled to its previous states $\theta(T-1)$ and
$\theta(T-\delta T-1)$. The virtual nodes thus reflect a topology with sparse
interconnections realized by the delay of the physical node (Fig. 3A). The
interconnections are inherently nonlinear due to the sine function in Eq. 4,
which originates from the physical interaction between the active and the
immobile particle and naturally serves as the activation function in our RC.

The working principle of our RC is now illustrated in Fig. 3B. For simplicity
of the discussion, we consider a RC with a single physical node, scalar input
$X_{n}$ and output $Y_{n}$ at the $n$-th computation step (for general case
see Sec. 4.1 of Supplementary Information). The input layer $\mathbf{u}_{n}$
is generated via a matrix of input weights
$\mathbb{W}_{\mathrm{in}}\in\mathbb{R}^{P_{\mathrm{in}}\times 2}$

$\mathbf{u}_{n}=\mathbb{W}_{\mathrm{in}}[b_{\mathrm{in}},\,X_{n}]^{\mathrm{T}},$
(5)

with the scalar $b_{\mathrm{in}}$ as the input bias. $\mathbf{u}_{n}$ is a
one-dimensional array containing $P_{\mathrm{in}}$ elements, which are
sequentially input into the physical node as the angle $u(T)$ of the active
particle in Eq. 3. This operation is equivalent to a time-multiplexing of the
input $X_{n}$ into $P_{\mathrm{in}}$ virtual nodes with
$\mathbb{W}_{\mathrm{in}}$ as the mask, as commonly applied in continuous-time
single node physical RC approaches Appeltant2011 ; Appeltant2014 . The output
$Y_{n}$ is derived as a linear combination of a scalar bias
$b_{\mathrm{out}}$, the input signal $X_{n}$, and the node states $\theta$ of
the past $P_{\mathrm{out}}$ timesteps using an output weight matrix
$\mathbb{W}_{\mathrm{out}}\in\mathbb{R}^{1\times(2+P_{\mathrm{out}})}$

$Y_{n}=\mathbb{W}_{\mathrm{out}}[b_{\mathrm{out}},\,X_{n},\,\theta(T=0),\,\theta(T=-1),\cdots,\,\theta(T=1-P_{\mathrm{out}})]^{\mathrm{T}},$
(6)

where $T=0$ denotes the time of the current computation step. The output
weight matrix $\mathbb{W}_{\mathrm{out}}$ is the only quantity to be trained
in the RC framework to tune the output towards the target signal. It is
calculated via a ridge regressions Lukosevicius2012  in each computation
cycle.

As compared to conventional time-multiplexed physical RCs where
$P_{\mathrm{in}}$ = $P_{\mathrm{out}}$, we explicitly allow these two
parameters to be independent and freely adjustable (see Fig. 3B). In
particular, we set $P_{\mathrm{out}}\gg P_{\mathrm{in}}$ for our RC. By this
means, each output is not only derived from the reservoir states of its
corresponding step, but also from previous $n_{\mathrm{hist}}$ =
$P_{\mathrm{out}}/P_{\mathrm{in}}-1$ steps. It will be demonstrated in the
next sections that this setting enables us to carry out the RC with a good
stability of the output and, most importantly, an effective reduction of the
impact of the intrinsic noise.

Figure 2:  Impulse response of teh active particle node. A Respons of
$\theta(t)$ to an impulsive input $u(t)$ = $\delta(t)$ measured in the
experiment with $v_{0}\delta t/R_{0}$ of 0.05 and 0.16. The grey curves denote
the measured $\theta(t)$ traces. They strongly fluctuate due to the Brownian
motion of the active particles. The red curves denote the means of 500
$\theta(t)$ traces in each case. B Comparison of the mean impulse responses of
$\theta(t)$ from the experiment (solid lines, same as in panel (A)) to the
ones obtained from the deterministic simulation (dashed lines). C Effective
potential $U_{\mathrm{eff}}(\theta)$ with different $v_{0}\delta t/R_{0}$.
$U_{\mathrm{eff}}(\theta)$ transitions from a single well to a double well
form at the transition point, with its local minima positions bifurcating from
zero to two opposite values ($\theta_{+}$, $\theta_{-}$). After perturbations
$\Delta\theta$ (solid arrows), $\theta$ relaxes to one of the minima (dashed
arrows). D Relaxation time $\tau_{\theta}$ of $\theta$ as function of
$v_{0}\delta t/R_{0}$ evaluated in deterministic simulations with $\delta
t=0.6\,$\mathrm{s}$$ varying $v_{0}$ from $0\,$\mathrm{\SIUnitSymbolMicro
m}\text{\,}{\mathrm{s}}^{-1}$$ to $2.74\,$\mathrm{\SIUnitSymbolMicro
m}\text{\,}{\mathrm{s}}^{-1}$$. $\theta$ is perturbed with $\Delta\theta$ from
one of the states $\theta_{+,-}$ at $t=0\,$\mathrm{s}$$, then relaxes to
$\theta(\tau_{\theta})-\theta_{+,-}$ = $\Delta\theta/10$ . $\tau_{\theta}$
diverges to infinity at the transition point of $U_{\rm eff}(\theta)$ (dash-
dotted line) at $v_{0}\delta t/R_{0}=0.46$.

This single physical node architecture can be further extended to multiple
physical nodes operated in parallel. In our experiment, we control
$N_{\mathrm{node}}=10$ independent physical nodes in one sample
simultaneously. The transient state $\theta$ of the 10 physical nodes together
constitutes the virtual nodes of the reservoir. Fig. 1D and Supplementary
Movie 1 show the real time image and video of the sample in the experiment.
Fig. 3C plots an exemplary trace of $\theta(t)$ of four physical nodes (blue
lines) driven by an external signal $X$ (red line) in the experiment.

The configuration of the RC is optimized in a simulation for the best
performance and then applied to the experiment. The input weights are selected
from a binary distribution $\mathbb{W}^{i}_{\mathrm{in}}$ $\in\\{-2,2\\}$,
which results in a better noise resistance of the RC than using
$\mathbb{W}^{i}_{\mathrm{in}}$ $\in\\{-1,1\\}$. The details of the RC
configuration are described in Sec. 4 of Supplementary Information.

### Chaotic series prediction

We test our RC with the free-running prediction of the chaotic Mackey-Glass
series (MGS). The MGS is generated by the delay differential equation

$\frac{\,\mathrm{d}S(n)}{\,\mathrm{d}n}=\alpha\frac{S(n-\tau)}{1+S(n-\tau)^{\beta}}-\gamma
S(n),$ (7)

which was introduced to model the complex dynamics of physiological feedback
systems Mackey1977 . It has been widely used as a benchmark task for series
forecasting Jaeger2004 ; Antonik2017b ; Fujii2017 ; Dong2020 . With parameters
$\alpha$ = 0.2, $\beta$ = 10, $\gamma$ = 0.1, and a delay parameter $\tau$ =
17, the MGS exhibits a chaotic behavior with the a Lyapunov exponent of around
0.006 Jaeger2004 . The performance of the prediction is evaluated by the
normalized root-mean-square error (NRMSE, see Sec. S4, S5 of Supplementary
Information for details). Fig. 4 shows the results of the MGS predictions by
our RC in experiments and simulations.

Figure 3:  Architecture of the reservoir computer. A Topology of the reservoir
with a single physical node with a discrete-time delay $\delta T=1$ in this
example. The angle $\theta$ of the active particle at different timesteps $T$
is considered as a virtual node. Each virtual node state $\theta(T)$ is
nonlinearly coupled to previous states $\theta(T-1)$ and $\theta(T-\delta
T-1)$ through the physics of the system (Eq. 4). B Sketch of the information
processing in the single node RC with $P_{\mathrm{in}}=2$ and
$P_{\mathrm{out}}=4$ in this example. The external signal $X$ of each
computation step is multiplexed by a weight matrix $\mathbb{W}_{\mathrm{in}}$
to $P_{\mathrm{in}}$ elements, which are sequentially input into the node as
the perturbation $u(T)$ of the active particle (Fig. 1C). Each output is
linearly derived from the $\theta$ states of previous $P_{\mathrm{out}}$
timesteps using a weight matrix $\mathbb{W}_{\mathrm{out}}$ that is trained
using ridge regression. Biases $b_{\mathrm{in/out}}$ for the input and output
layers are not plotted for the sake of simplicity. C Examples of $\theta(t)$
traces of four physical nodes (blue lines) driven by an external signal $X(t)$
(red line, see Eq. 7) measured in the experiment.

#### Simulated prediction

The deterministic simulations have been carried out with a very small
reservoir with only $N_{\mathrm{node}}P_{\mathrm{in}}$ = 20 virtual nodes but
using a large $P_{\mathrm{out}}$ = 400, i.e., $n_{\mathrm{hist}}$ = 199
historical reservoir states for each output. The simulation result shows a
very good prediction of the target MGS up to around 900 steps (corresponding
to 5.4 Lyapunov time) with a NRMSE of 6.7$\times$10-2 (Fig. 4A). Similarly,
Sec. S7.2 of the Supplementary Information also describes the prediction of a
three-dimensional chaotic Lorenz series by our RC in a deterministic
simulation. These results signify the capability of our architecture for
chaotic systems predictions.

#### Experimental prediction

As compared to the simulations, the experimental RC performance is
significantly degraded as a result of the Brownian motion of the active
particles and the sensitivity of chaotic systems. The noise acts on the
particle angular position $\phi(T)$ (Eq. 3) and propagates to the virtual node
states $\theta(T)$ (Eq. 4). The signal-to-noise ratio (SNR) of the RC can
estimated by comparing the angular displacement of the active particle between
neighboring timesteps $\Delta\phi(T)=\phi(T)-\phi(T-1)$ in the simulation with
and without noise (see Sec. S6 of Supplementary Information) and reveals an
extremely low value (SNR = 1.9 (2.78 dB)). Fig. 4B depicts the experimental
results of 50 MGS predictions with identical configurations
($\mathbb{W}_{\mathrm{out}}$ trained for each repetition). The outputs $Y(n)$
(gray curves) from different experimental predictions exhibit a notable
influence of the noise. The mean of the outputs $\overline{Y}(n)$ (blue curve)
can reproduce the fundamental period of the target MGS for round 200 steps and
details of the series for about 50 steps, which is, considering the extremely
low SNR, still remarkable. These results obtained for the experimental RC are
further underscored by the very good prediction of a non-chaotic periodic
trigonometric series (see Sec. S7.1 of Supplementary Information).

Figure 4:  Results of free-running predictions of the Mackey-Glass series. A
Prediction (blue line) of the target Mackey-Glass series (MGS, red line) by
the RC in the deterministic simulation using
$N_{\mathrm{node}}P_{\mathrm{in}}$ = 20 virtual nodes and $n_{\mathrm{hist}}$
= 199 historical states (i.e. $P_{\mathrm{out}}$ = 400). The simulations use a
active particle speed of $v_{0}=2.78\,$\mathrm{\SIUnitSymbolMicro
m}\text{\,}{\mathrm{s}}^{-1}$$ and a delay of $\delta t=0.05\,$\mathrm{s}$$.
Other parameters of the RC are given in Sec. 4.3 of Supplementary Information.
B Experimental results of MGS predictions with the same RC configuration as in
panel (A). The gray curves denote the RC outputs from 50 repeated predictions
with strong fluctuations due to Brownian motion. The blue curve represents the
mean of the output traces. C–G Comparison of the experimental results with
$P_{\mathrm{out}}$ from 200 to 400 evaluated by 50 repeated predictions for
each $P_{\mathrm{out}}$. C Means (colored lines) and corresponding standard
deviations (colored areas) of the RC output traces. D The deviation between
the mean of predictions and the target MGS, and E the standard deviations
$\sigma$ of the predictions versus the step $n$. The curves are smoothed via
100-step moving average. F NRMSE of 200 steps predictions from experiments and
stochastic simulations versus $P_{\mathrm{out}}$. G Means of the output
weights $\mathbb{W}_{\mathrm{out}}$ trained in experiments. The elements of
$\mathbb{W}_{\mathrm{out}}$ are in turn the weight for the output bias
$b_{\mathrm{out}}$, the input signal $X$, and the historical $\theta$ states
of the physical nodes (see Sec. S4.1 of Supplementary Information). Each
$N_{\mathrm{node}}P_{\mathrm{in}}$ = 20 weights for $\theta$ correspond to one
computation step, which is denoted by $n$ on the top axis with the negative
values representing the past. The peak marked by the arrow indicates the high
contribution of the historical reservoir states of around $n$ = -17, which
reveals the property of the target MGS with the delay parameter $\tau$ = 17
(Eq. 7).

The performance under these very noisy conditions becomes possible by the
special architecture using a number of historical states $n_{\mathrm{hist}}$
for the output derivation. Fig. 4C–G demonstrate the experimental results of
the RC with $P_{\mathrm{out}}$ varying from 200 to 400, corresponding to
$n_{\mathrm{hist}}=99\ldots 199$ historical states. The accuracy of the mean
of the predictions is improved by increasing $P_{\mathrm{out}}$ (Fig. 4C,
colored lines). Changing $P_{\mathrm{out}}$ from 200 to 400 results in a
decrease of the difference between $\overline{Y}(n)$ and the target MGS by 22%
(evaluated by the root-mean-square (RMS) of 200 steps $|\overline{Y}$ \-
$Y^{\mathrm{target}}|$ in Fig. 4D). The standard deviation $\sigma_{Y}$ of the
prediction decreases by 32% (200 steps RMS of $\sigma_{Y}$ in Fig. 4E) and the
resulting NRMSE of 200 step predictions diminishes by 28%. This trend of lower
prediction error for the experimental system is also picked up by stochastic
simulations with an NRMSE improvement of 11% (Fig. 4F).

This finding is striking as increasing $P_{\mathrm{out}}$ does not add more
information to the RC, nor increases the reservoir dimensionality, which is
determined by the number of independent variables Dambre2012  of the
reservoir. During the computation of the RC, the node states of each step are
nonlinearly transformed and mapped into the states of the following steps
Dambre2012 , while they attenuate in magnitude due to the fading memory.
Conventional RCs derive the output from the reservoir state of its
corresponding step, which implicitly contains the information of the
historical states. Whereas in our RC, the historical reservoir states can
directly contribute to the output.

The actual contributions of the historical states are determined via the
training of the output weights $\mathbb{W}_{\mathrm{out}}$. The states
correlating more to the current output obtain higher weight magnitudes. Fig.
4G plots the $\mathbb{W}_{\mathrm{out}}$ trained in experiments with
$P_{\mathrm{out}}$ of 200 and 400. The first 2000 elements of
$\mathbb{W}_{\mathrm{out}}$ corresponding to the near past reservoir states
($n$ from 0 to -99) have similar structures in both cases. As compared to the
near past, the far past states ($n$ from -100 to -199 for $P_{\mathrm{out}}$ =
400) correspond to smaller weights in amplitude, indicating the weaker
correlations to the current output. The highest magnitude of
$\mathbb{W}_{\mathrm{out}}$ appears at around $n$ = $-$17, which coincides
with the delay parameter $\tau$ = 17 of the target MGS (Eq. 7). The trained
$\mathbb{W}_{\mathrm{out}}$ thereby partly reveal the property of the target
signal.

#### Impact of noise

These results suggests, that there might be an optimal relation between the
used number of historical states and the error of the RC to stabilize the
prediction and to reduce the impact of the inherent noise due to Brownian
motion. To investigate this interelation we refer to deterministic and
stochastic simulations. Fig. 5A, B display the RC performance as measured by
the NRMSE of the predictions as functions of $P_{\mathrm{in}}$ and
$P_{\mathrm{out}}$. Without noise (Fig. 5A), the results indicate a larger
error with low $P_{\mathrm{in}}$ or $P_{\mathrm{out}}$ and poor performances
for $P_{\mathrm{in}}<$ 2 and $P_{\mathrm{out}}<$ 50\. The performance is
improving for increasing $P_{\mathrm{in/out}}$ and with an approximate linear
correlation to obtain optimal results Fig. 5A. Best results for the
deterministic system are obtained with
$n_{\mathrm{hist}}=P_{\mathrm{out}}/P_{\mathrm{in}}-1=60$.

For the noisy active particles nodes (Fig. 5B), the RC outputs are unstable
for $P_{\mathrm{out}}<100$. The free-running prediction has a high probability
to yield fast diverging outputs to yield a large NRMSE (bottom subplot). With
$P_{\mathrm{out}}\geq 100$ a stable output is achieved. A trend towards lower
NRMSE with higher $n_{\mathrm{hist}}$ can be observed in agreement with our
experimental results (Fig. 4F). The results again suggest a linear relation
between $P_{out}$ and $P_{in}$ to obtain a minimum error. Thus, both
experiments and simulations confirm that using historical states
$n_{\mathrm{hist}}$ for the prediction of our RC improves the stability and
the quality of the prediction even under extremely noisy conditions.

Figure 5:  Simulation results of the RC performance as function of the RC
configuration. NRMSE of 200 steps MGS prediction versus $P_{\mathrm{in}}$ and
$P_{\mathrm{out}}$ (A–B), $v_{0}/R_{0}$ and $\delta t$ (C–E).
$\mathbb{W}_{\mathrm{in}}$ and $b_{\mathrm{in}}$ are optimized (see Sec. S4.3
of supplementary Information) for each grid point. A Results of RC without
noise. The white dashed lines represent the contour lines of
$n_{\mathrm{hist}}=P_{\mathrm{out}}/P_{\mathrm{in}}-1$. B Results of RC with
noises and the same parameters as in panel (A). Each grid point is evaluated
by 50 repeated predictions. The outputs of RC with $P_{\mathrm{out}}<100$ are
unstable, and result in large NRMSE as plotted separately in the bottom
subplot. C, D Results of RC without noise. The elements of
$\mathbb{W}_{\mathrm{in}}$ are selected from $\\{-1,1\\}$ (C) and $\\{-2,2\\}$
(D) respectively. The white dash-dotted curves denote the contour lines of
$v_{0}\delta t/R_{0}$. The yellow dashed curves denote the transition points
of $U_{\mathrm{eff}}(\theta)$ from single well to double well form. An
instrumental feedback latency $\delta t_{F}$ = 0.125 s is considered, thereby
the transition curves do not coincide with the $v_{0}\delta t/R_{0}$ contour
lines. For more details see Sec. S1–S2 of Supplementary Information. The
subplot on the right illustrates the corresponding $U_{\mathrm{eff}}(\theta)$,
also $\theta$ perturbations and relaxations. E Comparison of NRMSE results
from panel (C) and (D) as function of $v_{0}\delta t/R_{0}$ with
$v_{0}/R_{0}=1.01$. The dashed line denotes the transition point of
$U_{\mathrm{eff}}$.

#### Impact of the dynamical properties of the nodes

Our active particle recurrent node also provides the opportunity to tune the
dynamical response of the physical node across the transition of the pitchfork
bifurcation with either particle speed $v_{0}$ or delay $\delta t$. The tuning
varies the memory of each node as indicated in Fig. 2D by the relaxation time
$\tau_{\theta}$ and thereby also the coupling of the virtual nodes as stated
by Eq. 4.

The impact of this change on the performance of the RC is depicted in Fig.
5C–E with the NRMSE for the MGS predictions as function of $v_{0}/R_{0}$ and
$\delta t$ resulting from deterministic simulations. As discussed above,
$\theta_{0}$ determines the asymptotic behavior after an impulsive
perturbation of $\theta$. The optimum around $\theta_{0}\approx 0.2$ and the
correlation with $\theta_{0}={\rm const.}$ indicate that the dynamical
properties of the physical node, in particular the fading memory as
characterized by the relaxation time $\tau_{\theta}$ (Fig. 2B, D) are indeed a
key factor of the RC performance for this task. The optimal performance of the
deterministic system is found below the transition point (yellow line) where
the approximate effective potential is largely determined by the parabolic
part (Fig. 5C, right) Wang2023 . While the relaxation time is largest at the
transition point and an input $u$ perisists for the longest time, the NRMSE is
found to be largest considerbly beyond the transition point as indicated by
the yellow dashed line in (Fig. 5 C). In this region, the effective potential
is a double well potential allowing for an inconsistent response, i.e. an
input perturbation may result in a relaxation into either of the minima. If
the barrier of the double well potential increases further the NRMSE appears
to decrease again (Fig. 5E) presumably due to the fact that active particle
node is largely residing in one of the two potential wells only with a short
relaxation time.

##  I Discussion

We have demonstrated above the realization of a physical reservoir computer in
an experiment using self-propelled active microparticles. A retarded
propulsion towards an immobilized target particle creates a self-organized
non-linear dynamical system that is, despite the strong Brownian noise,
capable of predicting chaotic time series such as Mackey-Glass and Lorenz
series when used as a physical node in a reservoir computer. The key element
of this physical recurrent node is a time-delay realizing a retarded
interaction Mijalkov.2016 ; Wang2023 ; Chen.2023  that creates the fading
memory as a basic requirement for reservoir computing Jaeger2001 . It also
provides the coupling of the active particle dynamics to its past allowing to
implement virtual nodes living on a single physical active particle system via
time-multiplexing Appeltant2011 ; Appeltant2014 . The information processing
that is provided by such a single node may therefore extend the rare
simulation work on reservoir computing with active particle swarms with
interparticle coupling Lymburn2021 . Additionally, the nonlinearity that is
required for computations is an intrinsic physical property of our active
particle system and requires no extra treatment of the output signal
Cucchi2022 . Future work could introduce direct physical coupling between the
isolated physical nodes in our configuration, e.g., through hydrodynamic or
other interactions, to obtain more complex dynamical networks of interacting
synthetic active particles.

The dynamics of the physical recurrent node that is the basic unit of our
reservoir computer is controlled by a parameter containing the product of
activity (active particle speed) and time-delay and can be understood with a
simple Landau-like self-induced quartic effective potential, which exhibits a
bifurcation to a double well potential Wang2023 . The system thereby allows to
address the relation of reservoir performance and node relaxation time, were
we found a clear indication for an optimal performance for small delays below
the bifurcation point in deterministic simulations. While the fading memory
becomes extremely long at the bifurcation and one might expect the worst
performance of the system, it is observed that a double well potential with a
small barrier leads to inconsistencies in the relaxation dynamics that are
more severe even in the deterministic system. Interestingly, such double well
potentials have recently been discussed as non-linear stochastic-resonance-
based activation functions in an attempt to provide better stability of Echo-
State-Networks against noise Guo2018 ; Liao2021 .

Noise is an inherent property of our information processing units, as Brownian
motion causes strong fluctuations of the node state. Such noises are
inevitable at the smallest scales also in the context of biological
information processing Tsimring2014 , for instance in neurons Faisal2008 ,
both with positive and negative effect Stein2005 ; Guo2018 . In physical RC
approaches, noise is commonly a major limiting factor for the performances
Dambre2012 ; Soriano2013 ; Antonik2017b  although subtle noises are reported
to be beneficial as well Jaeger2001 ; Jaeger2004 ; Sussillo2009 ; Paquot2012 ;
Estebanez2019 ). Yet, general strategies for the noise suppression are unclear
except increasing the reservoir size Alata2020 , which is normally costly for
physical RCs. While the performance of our reservoir computer for the chaotic
system prediction is highly degraded by the noise due to the sensitivity of
chaotic systems, we have introduced the output generation including historical
node states providing remarkable stability and noise reduction evenunder low
signal to noise ratios (see Sec. S7.1 of Supplementary Information). This
architecture is not increasing the dimensionality of the reservoir, but gives
different weight to contributions of reservoir states from the past for the
current output and could be potentially useful in future reservoir computing
studies.

In summary, simple retarded interactions in synthetic active microparticle
systems can give rise to non-linear self-driven dynamics that form a basis for
information processing with active matter. Our reservoir computer highlights
this connection between information processing, machine learning and active
matter on the microscale and paves the way for new studies on noise in
reservoir computing. While we so far referred to isolated active recurrent
units, we envision that the high level of control of synthetic active matter
design will yield new emergent physical collective states that may leverage
the field of active synthetic dynamical systems for information processing.

## Methods

###  I.1 Sample preparation

The sample used in experiments contains two kinds of micro particles:
polystyrene (PS) particles (microParticles GmbH) of
$3~{}$\mathrm{\SIUnitSymbolMicro m}$$ diameter and the melamine formaldehyde
(MF) (microParticles GmbH) particles of $2.19~{}$\mathrm{\SIUnitSymbolMicro
m}$$ diameter, suspended in a water solution (Fig. 1B). Gold (Au) nano-
particles of around $8~{}$\mathrm{nm}$$ diameter uniformly distributed on the
surface of the MF particle, cover about 10% of the total surface area of the
latter. Two glass coverslips ($20\,\times\,20~{}$\mathrm{mm}$^{2}$ and
$24\,\times\,24~{}$\mathrm{mm}$^{2}$) confine a
$3~{}$\mathrm{\SIUnitSymbolMicro m}$$ thick sample layer in between. Due to
the surface tension of water, the PS particles are compressed and immobilized
on the coverslips serving as spacers to define the sample thickness. The PS
and MF-Au particles are separately added into two 2% Pluronic F-127 solutions.
After 30 minutes, the Pluronic concentration of both solutions is decreased to
0.02% by diluting twice. Each dilution is followed by an centrifugation then a
removal of part of the solution to keep the particle concentration.
$0.3~{}$\mathrm{\SIUnitSymbolMicro l}$$ sample of PS particles is pipetted on
one of the coverslips, then $0.3~{}$\mathrm{\SIUnitSymbolMicro l}$$ sample of
MF-Au particles is pipetted in the droplet of PS particles. The sample is then
covered carefully with a second coverslip. The edges of the sample are sealed
by polydimethylsiloxane (PDMS) to prevent leakage and evaporation, also to
relieve the liquid flow inside the sample. The experiment starts about one
hour after sample preparation to wait until residual liquid flows in the
sample ceased.

###  I.2 Experimental setup

The experimental setup is illustrated in Fig. 1A. The micro MF-Au particles
are heated by a focused, continuous-wave laser with a wavelength of
$532~{}$\mathrm{nm}$$. The light from the laser module (CNI, MGL-H-532-1W) is
expanded by two tube lenses ($35,150~{}$\mathrm{mm}$$ focal lengths) in the
beam size, then guided by mirrors to a high-speed reflective Spacial Light
Modulator (SLM, Meadowlark Optics, HSP512-532), which modulates the phase of
the reflected laser. The reflected laser is then guided through two lenses
($500,300~{}$\mathrm{mm}$$ focal lengths) to an inverted microscope (Olympus,
IX73). A small opaque dot on a glass window located at the focal point of the
$500~{}$\mathrm{mm}$$ lens, serves as a mask to block the unmodulated laser
reflected by the SLM top surface. The laser in the microscope is reflected by
a dichroic beam splitter (Omega Optical, 560DRLP), then focused by an
objective lens (100x, Olympus, UPlanFL N x100/1.30, Oil, Iris, NA. 0.6–1.3) on
the sample plane with a beam width at half maximum about
$0.6~{}$\mathrm{\SIUnitSymbolMicro m}$$.

The sample is illuminated by white light from a LED lamp (Thorlabs, SOLIS-3C)
through an oil-immersion dark-field condenser (Olympus, U-DCW, NA 1.2–1.4).
The image of the sample is projected by the objective lens and a tube lens
($180~{}$\mathrm{mm}$$ focal length) inside the microscopy stand as well as
two additional lenses ($100,150~{}$\mathrm{mm}$$ focal lengths) outside the
microscopy stand to a camera (Hamamatsu digital sCMOS, C11440-22CU). The
numerical aperture (NA) of the objective is set to a value below the minimal
NA of the dark-field condenser. Two filters (EKSMA Optics 246-2506-532,
Thorlabs FESH0800) in front of the camera block the back reflections of the
laser from the microscope. A desktop PC (Intel(R) Core™ i7-7700K CPU
@4$\times$4.20 GHz, NVIDIA GeForce GTX 1050Ti) with a LabVIEW program (v.
2019) analyses the images, records data, and manipulates the active particles
by controlling the laser through the phase pattern on the SLM. More details
are given in Sec. S1 of Supplementary Information.

## Acknowledgement

This work is supported by Center for Scalable Data Analytics and Artificial
Intelligence (Scads.AI) Dresden/Leipzig.

## Competing interests

The authors have no competing interests.

## Author Contributions

## Data availability

All data in support of this work is available in the manuscript or the
supplementary materials. Further data and materials are available from the
corresponding author upon request.

## References

  * (1) Tkačik, G. & Bialek, W.  Information Processing in Living Systems.  _Annu Rev Condens_ 7, 1–29 (2014).
  * (2) Wadhams, G. H. & Armitage, J. P.  Making sense of it all: bacterial chemotaxis.  _Nat Rev Mol Cell Biol_ 5, 1024–1037 (2004).
  * (3) Celani, A. & Vergassola, M.  Bacterial strategies for chemotaxis response.  _Proc Natl Acad Sci USA_ 107, 1391–1396 (2010).
  * (4) Cosentino, C. & Bates, D.  _Feedback Control in Systems Biology_ (CRC Press, Boca Raton, 2019).
  * (5) Knudsen, E. I., Lac, S. & Esterly, S. D.  Computational Maps in the Brain.  _Annu Rev Neurosci_ 10, 41–65 (1987).
  * (6) Marković, D., Mizrahi, A., Querlioz, D. & Grollier, J.  Physics for neuromorphic computing.  _Nat Rev Phys_ 2, 499–510 (2020).
  * (7) Abiodun, O. I. _et al._ State-of-the-art in artificial neural network applications: A survey.  _Heliyon_ 4, e00938 (2018).
  * (8) Jaeger, H.  The "echo state" approach to analysing and training recurrent neural networks-with an erratum note’.  _Bonn, Germany: German National Research Center for Information Technology
GMD Technical Report_ 148 (2001).

  * (9) Maass, W., Natschläger, T. & Markram, H.  Real-Time Computing Without Stable States: A New Framework for Neural Computation Based on Perturbations.  _Neural Computation_ 14, 2531–2560 (2002).
  * (10) Verstraeten, D., Schrauwen, B., D’Haene, M. & Stroobandt, D.  An experimental unification of reservoir computing methods.  _Neural Netw_ 20, 391–403 (2007).
  * (11) Paquot, Y., Dambre, J., Schrauwen, B., Haelterman, M. & Massar, S.  Reservoir computing: a photonic neural network for information processing.  _Nonlinear Optics and Applications IV_ 77280B–77280B–12 (2010).
  * (12) Appeltant, L. _et al._ Information processing using a single dynamical node as complex system.  _Nat Commun_ 2, 468 (2011).
  * (13) Lukoševičius, M. & Jaeger, H.  Reservoir computing approaches to recurrent neural network training.  _Comput Sci Rev_ 3, 127–149 (2009).
  * (14) Lukoševičius, M.  A Practical Guide to Applying Echo State Networks.  In Montavon, G., Orr, G. B. & Müller, K.-R. (eds.) _Neural Networks: Tricks of the Trade_ , vol. 7700, 659–686 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012).
  * (15) Tanaka, G. _et al._ Recent advances in physical reservoir computing: A review.  _Neural Netw_ 115, 100–123 (2019).
  * (16) Nakajima, K.  Physical reservoir computing—an introductory perspective.  _Jpn J Appl Phys_ 59, 060501 (2020).
  * (17) Rafayelyan, M., Dong, J., Tan, Y., Krzakala, F. & Gigan, S.  Large-Scale Optical Reservoir Computing for Spatiotemporal Chaotic Systems Prediction.  _Phys Rev X_ 10, 041037 (2020).  eprint 2001.09131.
  * (18) Larger, L. _et al._ Photonic information processing beyond Turing: An optoelectronic implementation of reservoir computing.  _Opt Express_ 20, 3241–3249 (2012).
  * (19) Paquot, Y. _et al._ Optoelectronic Reservoir Computing.  _Sci Rep_ 2, 287 (2012).
  * (20) Van der Sande, G., Brunner, D. & Soriano, M. C.  Advances in photonic reservoir computing.  _Nanophotonics_ 6, 561–576 (2017).
  * (21) Coulombe, J. C., York, M. C. A. & Sylvestre, J.  Computing with networks of nonlinear mechanical oscillators.  _PLoS ONE_ 12, e0178663 (2017).
  * (22) Dion, G., Mejaouri, S. & Sylvestre, J.  Reservoir computing with a single delay-coupled non-linear mechanical oscillator.  _J Appl Phys_ 124, 152132 (2018).
  * (23) Wu, S., Zhou, W., Wen, K., Li, C. & Gong, Q.  Improved Reservoir Computing by Carbon Nanotube Network with Polyoxometalate Decoration.  _2021 IEEE 16th International Conference on Nano/Micro Engineered and Molecular Systems (NEMS)_ 00, 994–997 (2021).
  * (24) Bechinger, C. _et al._ Active particles in complex and crowded environments.  _Rev Mod Phys_ 88, 045006 (2016).  eprint 1602.00081.
  * (25) Lavergne, F. A., Wendehenne, H., Bäuerle, T. & Bechinger, C.  Group formation and cohesion of active particles with visual perception–dependent motility.  _Science_ 364, 70–74 (2019).
  * (26) Bäuerle, T., Löffler, R. C. & Bechinger, C.  Formation of stable and responsive collective states in suspensions of active colloids.  _Nat Comm_ 11, 2547 (2020).
  * (27) Wang, X., Chen, P.-C., Kroy, K., Holubec, V. & Cichos, F.  Spontaneous vortex formation by microswimmers with retarded attractions.  _Nat Commun_ 14, 56 (2023).
  * (28) Woodhouse, F. G. & Dunkel, J.  Active matter logic for autonomous microfluidics.  _Nat Comm_ 8, 15169 (2017).  eprint 1610.05515.
  * (29) Colen, J. _et al._ Machine learning active-nematic hydrodynamics.  _Proc Natl Acad Sci USA_ 118, e2016708118 (2021).  eprint 2006.13203.
  * (30) Cichos, F., Gustavsson, K., Mehlig, B. & Volpe, G.  Machine learning for active matter.  _Nat Mach Intel_ 2, 94–103 (2020).
  * (31) Muiños-Landin, S., Fischer, A., Holubec, V. & Cichos, F.  Reinforcement learning with artificial microswimmers.  _Sci Robot_ 6, eabd9285 (2021).
  * (32) Liebchen, B. & Löwen, H.  Optimal navigation strategies for active particles.  _EPL_ 127, 34003 (2019).
  * (33) Schneider, E. & Stark, H.  Optimal steering of a smart active particle.  _EPL_ 127, 64003 (2019).  eprint 1909.03243.
  * (34) Colabrese, S., Gustavsson, K., Celani, A. & Biferale, L.  Flow Navigation by Smart Microswimmers via Reinforcement Learning.  _Phys Rev Lett_ 118, 158004 (2017).  eprint 1701.08848.
  * (35) Colabrese, S., Gustavsson, K., Celani, A. & Biferale, L.  Smart inertial particles.  _Phys Rev Fluids_ 3, 084301 (2018).  eprint 1711.05853.
  * (36) Gustavsson, K., Biferale, L., Celani, A. & Colabrese, S.  Finding efficient swimming strategies in a three-dimensional chaotic flow by reinforcement learning.  _Eur Phys J E_ 40, 110 (2017).  eprint 1711.05826.
  * (37) Lymburn, T., Algar, S. D., Small, M. & Jüngling, T.  Reservoir computing with swarms.  _Chaos_ 31, 033121 (2021).
  * (38) Fränzl, M., Muiños-Landin, S., Holubec, V. & Cichos, F.  Fully Steerable Symmetric Thermoplasmonic Microswimmers.  _ACS Nano_ 15, 3434–3440 (2021).
  * (39) Khadka, U., Holubec, V., Yang, H. & Cichos, F.  Active particles bound by information flows.  _Nat Commun_ 9, 3864 (2018).
  * (40) Goldenfeld, N.  Lectures on Phase Transitions and the Renormalization Group.  In _Anomalous Dimensions_ , 189–199 (2018).
  * (41) Boyd, S. & Chua, L.  Fading memory and the problem of approximating nonlinear operators with Volterra series.  _IEEE Trans Circuits Syst_ 32, 1150–1161 (1985).
  * (42) Maass, W., Natschläger, T. & Markram, H.  Fading memory and kernel properties of generic cortical microcircuit models.  _J Physiol Paris_ 98, 315–330 (2004).
  * (43) Dambre, J., Verstraeten, D., Schrauwen, B. & Massar, S.  Information Processing Capacity of Dynamical Systems.  _Sci Rep_ 2, 514 (2012).
  * (44) Appeltant, L., Van der Sande, G., Danckaert, J. & Fischer, I.  Constructing optimized binary masks for reservoir computing with delay systems.  _Sci Rep_ 4, 3629 (2014).
  * (45) Mackey, M. C. & Glass, L.  Oscillation and chaos in physiological control systems.  _Science_ 197, 287–289 (1977).
  * (46) Jaeger, H. & Haas, H.  Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy in Wireless Communication.  _Science_ 304, 78–80 (2004).
  * (47) Antonik, P., Haelterman, M. & Massar, S.  Brain-Inspired Photonic Signal Processor for Generating Periodic Patterns and Emulating Chaotic Systems.  _Phys Rev Appl_ 7, 054014 (2017).
  * (48) Fujii, K. & Nakajima, K.  Harnessing Disordered-Ensemble Quantum Dynamics for Machine Learning.  _Phys. Rev. Appl._ 8, 024030 (2017).
  * (49) Dong, J., Rafayelyan, M., Krzakala, F. & Gigan, S.  Optical Reservoir Computing Using Multiple Light Scattering for Chaotic Systems Prediction.  _IEEE J Sel Top Quantum Electron_ 26, 1–12 (2020).
  * (50) Mijalkov, M., McDaniel, A., Wehr, J. & Volpe, G.  Engineering sensorial delay to control phototaxis and emergent collective behaviors.  _Phys Rev X_ 6, 1—16 (2016).  eprint 1511.04528.
  * (51) Chen, P.-C., Kroy, K., Cichos, F., Wang, X. & Holubec, V.  Active particles with delayed attractions form quaking crystallites (a).  _EPL_ 142, 67003 (2023).
  * (52) Cucchi, M., Abreu, S., Ciccone, G., Brunner, D. & Kleemann, H.  Hands-on reservoir computing: A tutorial for practical implementation.  _Neuromorph Comput Eng_ 2, 032002 (2022).
  * (53) Guo, D., Perc, M., Liu, T. & Yao, D.  Functional importance of noise in neuronal information processing.  _EPL_ 124, 50001 (2018).
  * (54) Liao, Z., Wang, Z., Yamahara, H. & Tabata, H.  Echo state network activation function based on bistable stochastic resonance.  _Chaos, Solitons & Fractals_ 153, 111503 (2021).
  * (55) Tsimring, L. S.  Noise in biology.  _Rep Prog Phys_ 77, 026601 (2014).
  * (56) Faisal, A. A., Selen, L. P. J. & Wolpert, D. M.  Noise in the nervous system.  _Nat Rev Neurosci_ 9, 292–303 (2008).
  * (57) Stein, R. B., Gossen, E. R. & Jones, K. E.  Neuronal variability: Noise or part of the signal?  _Nat Rev Neurosci_ 6, 389–397 (2005).
  * (58) Soriano, M. C. _et al._ Optoelectronic reservoir computing: Tackling noise-induced performance degradation.  _Opt Express_ 21, 12 (2013).
  * (59) Sussillo, D. & Abbott, L.  Generating Coherent Patterns of Activity from Chaotic Neural Networks.  _Neuron_ 63, 544–557 (2009).
  * (60) Estébanez, I., Fischer, I. & Soriano, M. C.  Constructive Role of Noise for High-Quality Replication of Chaotic Attractor Dynamics Using a Hardware-Based Reservoir Computer.  _Phys Rev Appl_ 12, 034058 (2019).
  * (61) Alata, R., Pauwels, J., Haelterman, M. & Massar, S.  Phase Noise Robustness of a Coherent Spatially Parallel Optical Reservoir.  _IEEE J Sel Top Quantum Electron_ 26, 1–10 (2020).