Harnessing synthetic active particles for physical reservoir computing

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 are the results 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 the 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 [1].The detection of low amounts of chemicals by a bacterium to navigate environments [2,3], the feedback mechanisms controlling and maintaining the function of organisms [4], or the highly sophisticated computations in large biological neural networks in the brain [5] 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 [? ].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 [7].Reservoir computers employ sparsely and statically connected recurrent nodes [8][9][10] or even a single node by using time-multiplexing [11,12] 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 [13,14].
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 [15][16][17] using optoelectronic oscillators [18][19][20], mechanical oscillators [21,22], carbon nanotubes [23] or passive soft bodies [15] 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 [24], that mimic, for example, the dynamics of swarms [25][26][27].Information processing [28,29] and learning [30] in experiment [31] and simulation [32][33][34][35][36] 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 [37].
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.

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 [8,10].We realize a simple recurrent node with the help of a single synthetic active particle as a microscopic model for motile active matter [24].The active particle is a polymer microbead of 2.19 µm diameter with the surface decorated by gold nanoparticles (about 8 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 [38] 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 δt is added to the intrinsic feedback latency to control the active particle [27,39] and to introduce a finite reaction time, as is inherent in many processes in living systems.Correspondingly, the attraction F(t) experienced by the particle at time t is determined by its position at previous time t − δt, where 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 θ(t) = φ(t) − φ(t − δt) during the delay time, which results in a transient rotational motion of the active particle around the immobile target (Fig. 1C), where φ denotes the angular position of the particle.The dynamics of the angular position can then be described by a nonlinear delay differential equation assuming that the active and the immobile particle (radius of a act and a imm ) are in physical contact, i.e.R 0 = a act + a 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 µm 2 s −1 , giving rise to a Pećlet number of Pe = a 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 F (Fig. 1C).For its function as a recurrent node, the dynamics of the angle θ(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 eff (θ) [27], which resembles a generic Landau-type description [27,40].The shape of the potential is controlled by a dimensionless parameter θ 0 = v 0 δt/R 0 .With increasing θ 0 , U eff (θ) transitions from a near parabolic shape with a minimum at θ = 0 to a symmetric double well shape with minima at θ + and θ − in a pitchfork bifurcation (see Fig. 2C and Sec.S2 of Supplementary Information).
A perturbation of θ (solid arrow in Fig. 2C), will result in a relaxation with a dynamics determined by the control parameter θ 0 .We have experimentally determined the response θ(t) to an impulsive perturbation u(t) = δ(t) for different values of θ 0 (Fig. 2A, B).The individual experimental trajectories θ(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 θ(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 τ θ extracted from the deterministic simulations reveals the expected strong increase around the transition point.
Tuning the control parameter θ 0 , i.e. activity (v 0 ) and/or delay (δ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 θ(t) demonstrated in Fig. 2B represents a basic requirement for reservoir computing [9,[41][42][43].Based on the time-delay it also allows the setup of coupled recurrent where D 0 ≈ 0.0642 μm 2 s −1 denotes the translational diffusion coefficient of the active particle and η(t) white noise.To solve this equation, we approximated _ ϕðtÞδt by θ(t) and expanded the sinðϕðtÞ À ϕðt À δtÞÞ in a Taylor series around δt = 0 up to the third order in δt.We dropped the term proportional to ϕ ⃛ (t) to secure the stability of the resulting equation 38 (for details, see Sec. 3 of Supplementary Information).The resulting noise term ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 8D 0 =ðω 0 δtRÞ 2 p turned out to be inaccurate compared to experimental and simulation data.We, therefore, introduce an effective diffusion coefficient D θ as a free parameter in the noise term in Eq. ( 4) to describe the rotation of the active particle around the target as the angular Brownian motion with Eq. ( 4) yields the stationary solutions 0 and θ ± with the bifurcation point ω 0 δt = 1, for the transition from a non-rotational to a rotational state.The data points in Fig. 3B display the experimentally obtained maxima of the histograms p(θ) of the propulsion angle (see Fig. 2C) as a function of ω 0 δt.The transition points in the experiments are located at lower values of the control parameter ω 0 δt, due to the mentioned instrumental delay Δt in the feedback loop of the experimental setup.This instrumental delay between the most recent exposure to the camera and the laser positioning affects the motion direction beyond the programmed delay δt 34,39 , causing an earlier onset of the transition to a stable rotation.The dashed line in Fig. 3B shows the theoretical prediction, which includes both the instrumental delay Δt and the programmed delay δt, as detailed in the Supplementary Information (Eq.( 11)).
The Langevin equation (4) can be interpreted as a dynamical equation for the position θ of an overdamped Brownian particle with diffusion coefficient D θ in a quartic potential (see derivation in Sec. 3 of Supplementary Information), which allows to classify the observed instability of the isotropic state as a normal supercritical pitchfork bifurcation 40 .The potential can also directly be extracted from the experimental data (Fig. 3B) by fitting the histogram p(θ) with a (normalized) Boltzmann distribution expðÀUðθÞ=D θ Þ=Z at the effective temperature D θ .The effective temperature thus links the measured potential of mean force ÀD θ log pðθÞ to Eq. ( 6).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 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 F(t) with a time delay δt in its reaction.The direction of F(t) (dashed arrow) is determined by the previous active particle location r(t − δt).An additional angle u(t) between the particle propulsion direction (solid arrow) and 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.
nodes.In the discrete-time setting of our experiment with the sampling period ∆t, Eq. 2 can be rewritten as where T = t/∆t is an integer number representing the timestep.W denotes Gaussian random numbers with zero mean and unit variance.The evolution of θ then follows as with the discrete-time delay δT = δt/∆t.Refering to the concept of virtual nodes and timemultiplexing [12], we consider the transient state θ(T ) of a physical node at different timesteps as virtual nodes constituting the reservoir.Each virtual node state θ(T ) is, according to Eq. 4, coupled to its previous states θ(T − 1) and θ(T − δ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 u n is generated via a matrix of input weights with the scalar b in as the input bias.u n is a one-dimensional array containing P 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 in virtual nodes with W in as the mask, as commonly applied in continuous-time single node physical RC approaches [12,44].
The output Y n is derived as a linear combination of a scalar bias b out , the input signal X n , and the node states θ of the past P out timesteps using an output weight matrix where T = 0 denotes the time of the current computation step.The output weight matrix W 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 [14] in each computation cycle.
As compared to conventional time-multiplexed physical RCs where P in = P out , we explicitly allow these two parameters to be independent and freely adjustable (see Fig. 3B).In particular, we set P out P 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 hist = P out /P 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.This where D 0 ≈ 0.0642 μm 2 s −1 denotes the translational diffusion coefficient of the active particle and η(t) white noise.To solve this equation, we approximated _ ϕðtÞδt by θ(t) and expanded the points in the experiments are located at lower values of the control parameter ω 0 δt, due to the mentioned instrumental delay Δt in the feedback loop of the experimental setup.This instrumental delay  2 and the linked small color plots show the corresponding potentials of mean force, determined from the propulsion angle histograms in Fig. 2C, together with a fit of the refined analytical model, including the instrumental delay Δt (see Sec. 2.2 and 3 of Supplementary Information).The only free parameter for fitting is the effective temperature of the system.C Relaxation time τ of a single active particle as determined experimentally from the autocorrelation of the propulsion angle fluctuations (Eq.( 8), data points).The solid lines correspond to the refined version of the theoretical prediction (Eq.( 7)), including the instrumental delay Δt (see Sec.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 W i in ∈ {−2, 2}, which results in a better noise resistance of the RC than using W i 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 which was introduced to model the complex dynamics of physiological feedback systems [45].It has been widely used as a benchmark task for series forecasting [46][47][48][49].With parameters α = 0.2, β = 10, γ = 0.1, and a delay parameter τ = 17, the MGS exhibits a chaotic behavior with the a Lyapunov exponent of around 0.006 [46].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.

Simulated prediction
The deterministic simulations have been carried out with a very small reservoir with only N node P in = 20 virtual nodes but using a large P out = 400, i.e., n 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×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 φ(T ) (Eq. 3) and propagates to the virtual node states θ(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 ∆φ(T ) = φ(T )−φ(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 where D 0 ≈ 0.0642 μm 2 s −1 denotes the translational diffusion coefficient of the active particle and η(t) white noise.To solve this equation, we approximated _ ϕðtÞδt by θ(t) and expanded the sinðϕðtÞ À ϕðt À δtÞÞ in a Taylor series around δt = 0 up to the third order in δt.We dropped the term proportional to ϕ ⃛ (t) to secure the stability of the resulting equation 38 (for details, see Sec. 3 of Supplementary Information).The resulting noise term ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 8D 0 =ðω 0 δtRÞ 2 p turned out to be inaccurate compared to experimental and simulation data.We, therefore, introduce an effective diffusion coefficient D θ as a free parameter in the noise term in Eq. ( 4) to describe the rotation of the active particle around the target as the angular Brownian motion with Eq. ( 4) yields the stationary solutions 0 and θ ± with the bifurcation point ω 0 δt = 1, for the transition from a non-rotational to a rotational state.The data points in Fig. 3B display the experimentally obtained maxima of the histograms p(θ) of the propulsion angle (see Fig. 2C) as a function of ω 0 δt.The transition points in the experiments are located at lower values of the control parameter ω 0 δt, due to the mentioned instrumental delay Δt in the feedback loop of the experimental setup.This instrumental delay between the most recent exposure to the camera and the laser positioning affects the motion direction beyond the programmed delay δt 34,39 , causing an earlier onset of the transition to a stable rotation.The dashed line in Fig. 3B shows the theoretical prediction, which includes both the instrumental delay Δt and the programmed delay δt, as detailed in the Supplementary Information (Eq.( 11)).
The Langevin equation ( 4) can be interpreted as a dynamical equation for the position θ of an overdamped Brownian particle with diffusion coefficient D θ in a quartic potential (see derivation in Sec. 3 of Supplementary Information), which allows to classify the observed instability of the isotropic state as a normal supercritical pitchfork bifurcation 40 .The potential can also directly be extracted from the experimental data (Fig. 3B) by fitting the histogram p(θ) with a (normalized) Boltzmann distribution expðÀUðθÞ=D θ Þ=Z at the effective temperature D θ .The effective temperature thus links the measured potential of mean force ÀD θ log pðθÞ to Eq. ( 6).A Graphical construction of condition (3) for a transition from an non-rotational state (red-shaded region) to a rotational state (green-shaded region).The red line (sin θ) and the black dashed line with slope 1/(ω0δt) intersect at several θ.The solution θ = θ+ in the green region and its chirally inverse image θ− in the third quadrant (not shown) correspond to co-and counter-clockwise rotation.B Experimentally measured propulsion angles (maxima of the histograms in Fig. 2C) as a function of ω0δt, exhibiting a bifurcation at ω0δt ≈ 0.76.The dashed line corresponds to the analytical prediction of the theoretical model ( 5), neglecting the inevitable instrumental delay Δt.The solid line shows the solution of the refined theoretical model, which includes the instrumental delay Δt = 64 ms of our setup in addition to the programmed delay δt.The colored dots indicate the control parameter values studied in Fig. 2 and the linked small color plots show the corresponding potentials of mean force, determined from the propulsion angle histograms in Fig. 2C, together with a fit of the refined analytical model, including the instrumental delay Δt (see Sec. 2.2 and 3 of Supplementary Information).The only free parameter for fitting is the effective temperature of the system.C Relaxation time τ of a single active particle as determined experimentally from the autocorrelation of the propulsion angle fluctuations (Eq.( 8), data points).The solid lines correspond to the refined version of the theoretical prediction (Eq.( 7)), including the instrumental delay Δt (see Sec. 2.2 of Supplementary Information for details).The colored dots have the same meaning as in panel (B).D Transition rates between the two rotational states obtained from the experiments (circles) plotted with the predictions from Kramers' theory, Eq. ( 9), with a global fit parameter Dθ = 0.05 s −1 (solid line) and Dθ fitted to the probability distribution p(θ) separately for each value ω 0 δt (squares).Error bars represent the standard error.

A B C
FIG. 3. Architecture of the reservoir computer.A Topology of the reservoir with a single physical node with a discrete-time delay δT = 1 in this example.The angle θ of the active particle at different timesteps T is considered as a virtual node.Each virtual node state θ(T ) is nonlinearly coupled to previous states θ(T − 1) and θ(T − δT − 1) through the physics of the system (Eq.4).B Sketch of the information processing in the single node RC with P in = 2 and P out = 4 in this example.The external signal X of each computation step is multiplexed by a weight matrix W in to P 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 θ states of previous P out timesteps using a weight matrix W out that is trained using ridge regression.Biases 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 out does not add more information to the RC, nor increases the reservoir dimensionality, which is determined by the number of independent variables [43] 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 [43], 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 W out .The states correlating more to the current output obtain higher weight magnitudes.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 in and P out .Without noise (Fig. 5A), the results indicate a larger error with low P in or P out and poor performances for P in < 2 and P out < 50.
The performance is improving for increasing P in/out and with an approximate linear correlation to obtain optimal results Fig. 5A.Best results for the deterministic system are obtained with For the noisy active particles nodes (Fig. 5B), the RC outputs are unstable for P out < 100.The freerunning prediction has a high probability to yield fast diverging outputs to yield a large NRMSE (bottom subplot).With P out ≥ 100 a stable output is achieved.A trend towards lower NRMSE with higher n 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 hist for the prediction of our RC improves the stability and the quality of the prediction even under extremely noisy conditions.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 δt.The tuning varies the memory of each node as indicated in Fig. 2D by the relaxation time τ θ 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 δt resulting from deterministic simulations.
The optimum around θ 0 ≈ 0.2 and the correlation with θ 0 = const.indicate that the dynamical properties of the physical node, in particular the fading memory as characterized by the relaxation time τ θ (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) [27].
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 [27,50,51] that creates the fading memory as a basic requirement for reservoir computing [8].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 [12,44].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 [37].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 [52].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 [27].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 [53,54].
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 [55], for instance in neurons [56], both with positive and negative effect [53,57].In physical RC approaches, noise is commonly a major limiting factor for the performances [43,47,58] although subtle noises are reported to be beneficial as well [8,19,46,59,60]).Yet, general strategies for the noise suppression are unclear except increasing the reservoir size [61], 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.

Fig. 3 |
Fig. 3 | Transition to a rotational dynamical state for a single active particle.A Graphical construction of condition (3) for a transition from an non-rotational state (red-shaded region) to a rotational state (green-shaded region).The red line (sin θ) and the black dashed line with slope 1/(ω0δt) intersect at several θ.The solution θ = θ+ in the green region and its chirally inverse image θ− in the third quadrant (not shown) correspond to co-and counter-clockwise rotation.B Experimentally measured propulsion angles (maxima of the histograms in Fig. 2C) as a function of ω0δt, exhibiting a bifurcation at ω0δt ≈ 0.76.The dashed line corresponds to the analytical prediction of the theoretical model (5), neglecting the inevitable instrumental delay Δt.The solid line shows the solution of the refined theoretical model, which includes the instrumental delay Δt = 64 ms of our setup in addition to the programmed delay δt.The colored dots indicate the control parameter values studied in Fig.2and the linked small color plots show the corresponding potentials of mean force, determined from the propulsion angle histograms in Fig.2C, together with a fit of the refined analytical model, including the instrumental delay Δt (see Sec. 2.2 and 3 of Supplementary Information).The only free parameter for fitting is the effective temperature of the system.C Relaxation time τ of a single active particle as determined experimentally from the autocorrelation of the propulsion angle fluctuations (Eq.(8), data points).The solid lines correspond to the refined version of the theoretical prediction (Eq.(7)), including the instrumental delay Δt (see Sec. 2.2 of Supplementary Information for details).The colored dots have the same meaning as in panel (B).D Transition rates between the two rotational states obtained from the experiments (circles) plotted with the predictions from Kramers' theory, Eq. (9), with a global fit parameter Dθ = 0.05 s −1 (solid line) and Dθ fitted to the probability distribution p(θ) separately for each value ω 0 δt (squares).Error bars represent the standard error.

Fig. 3 |
Fig.3| Transition to a rotational dynamical state for a single active particle.A Graphical construction of condition (3) for a transition from an non-rotational state (red-shaded region) to a rotational state (green-shaded region).The red line (sin θ) and the black dashed line with slope 1/(ω0δt) intersect at several θ.The solution θ = θ+ in the green region and its chirally inverse image θ− in the third quadrant (not shown) correspond to co-and counter-clockwise rotation.B Experimentally measured propulsion angles (maxima of the histograms in Fig.2C) as a function of ω0δt, exhibiting a bifurcation at ω0δt ≈ 0.76.The dashed line corresponds to the analytical prediction of the theoretical model (5), neglecting the inevitable instrumental delay Δt.The solid line shows the solution of the refined theoretical model, which includes the instrumental delay Δt = 64 ms of our setup in addition to the programmed delay δt.The colored dots indicate the control parameter values studied in Fig.2and the linked small color plots show the 2.2 of Supplementary Information for details).The colored dots have the same meaning as in panel (B).D Transition rates between the two rotational states obtained from the experiments (circles) plotted with the predictions from Kramers' theory, Eq. (9), with a global fit parameter Dθ = 0.05 s −1 (solid line) and Dθ fitted to the probability distribution p(θ) separately for each value ω 0 δt (squares).Error bars represent the standard error.

FIG. 2 .
FIG. 2.Impulse response of teh active particle node.A Respons of θ(t) to an impulsive input u(t) = δ(t) measured in the experiment with v 0 δt/R 0 of 0.05 and 0.16.The grey curves denote the measured θ(t) traces.They strongly fluctuate due to the Brownian motion of the active particles.The red curves denote the means of 500 θ(t) traces in each case.B Comparison of the mean impulse responses of θ(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 eff (θ) with different v 0 δt/R 0 .U eff (θ) 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 (θ + , θ − ).After perturbations ∆θ (solid arrows), θ relaxes to one of the minima (dashed arrows).D Relaxation time τ θ of θ as function of v 0 δt/R 0 evaluated in deterministic simulations with δt = 0.6 s varying v 0 from 0 µm s −1 to 2.74 µm s −1 .θ is perturbed with ∆θ from one of the states θ +,− at t = 0 s, then relaxes to θ(τ θ ) − θ +,− = ∆θ/10 .τ θ diverges to infinity at the transition point of U eff (θ) (dash-dotted line) at v 0 δt/R 0 = 0.46.

Fig. 3 |
Fig.3| Transition to a rotational dynamical state for a single active particle.A Graphical construction of condition (3) for a transition from an non-rotational state (red-shaded region) to a rotational state (green-shaded region).The red line (sin θ) and the black dashed line with slope 1/(ω0δt) intersect at several θ.The solution θ = θ+ in the green region and its chirally inverse image θ− in the third quadrant (not shown) correspond to co-and counter-clockwise rotation.B Experimentally measured propulsion angles (maxima of the histograms in Fig.2C) as a function of ω0δt, exhibiting a bifurcation at ω0δt ≈ 0.76.The dashed line corresponds to the analytical prediction of the theoretical model (5), neglecting the inevitable instrumental delay Δt.The solid line shows the solution of the refined theoretical model, which includes the instrumental delay Δt = 64 ms of our setup in addition to the programmed delay δt.The colored dots indicate the control parameter values studied in Fig.2and the linked small color plots show the

bFIG. 4 .
FIG. 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 node P in = 20 virtual nodes and n hist = 199 historical states (i.e.P out = 400).The simulations use a active particle speed of v 0 = 2.78 µm s −1 and a delay of δt = 0.05 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 out from 200 to 400 evaluated by 50 repeated predictions for each P 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 σ 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 out .G Means of the output weights W out trained in experiments.The elements of W out are in turn the weight for the output bias b out , the input signal X, and the historical θ states of the physical nodes (see Sec. S4.1 of Supplementary Information).Each N node P in = 20 weights for θ 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 τ = 17 (Eq.7).

Fig. 4G plots
Fig. 4G plots the W out trained in experiments with P out of 200 and 400.The first 2000 elements of W 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 out = 400) correspond to smaller weights in amplitude, indicating the weaker correlations to the current output.The highest magnitude of W out appears at around n = −17, which coincides with the delay parameter τ = 17 of the target MGS (Eq.7).The trained W out thereby partly reveal the property of the target signal.

FIG. 5 .
FIG. 5. Simulation results of the RC performance as function of the RC configuration.NRMSE of 200 steps MGS prediction versus P in and P out (A-B), v 0 /R 0 and δt (C-E).W in and b 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 hist = P out /P 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 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 W in are selected from {−1, 1} (C) and {−2, 2} (D) respectively.The white dash-dotted curves denote the contour lines of v 0 δt/R 0 .The yellow dashed curves denote the transition points of U eff (θ) from single well to double well form.An instrumental feedback latency δt F = 0.125 s is considered, thereby the transition curves do not coincide with the v 0 δt/R 0 contour lines.For more details see Sec.S1-S2 of Supplementary Information.The subplot on the right illustrates the corresponding U eff (θ), also θ perturbations and relaxations.E Comparison of NRMSE results from panel (C) and (D) as function of v 0 δt/R 0 with v 0 /R 0 = 1.01.The dashed line denotes the transition point of U eff .