11institutetext:  Artificial Intelligence and Cyber Futures Institute, Charles
Sturt University, Bathurst, NSW 2795, Australia
Department of Mathematics, Swinburne University of Technology, Hawthorn,
Victoria 3122, Australia

# Reservoir computing based on solitary-like waves dynamics of film flows: a
proof of concept

Ivan S. Maksymov 11 Andrey Pototsky 221122

###### Abstract

Several theoretical works have shown that solitons—waves that self-maintain
constant shape and velocity as they propagate—can be used as a physical
computational reservoir, a concept where machine learning algorithms designed
for digital computers are replaced by analog physical systems that exhibit
nonlinear dynamical behaviour. Here we propose and experimentally validate a
novel reservoir computing (RC) system that for the first time employs
solitary-like (SL) waves propagating on the surface of a liquid film flowing
over an inclined surface. We demonstrate the ability of the SL wave RC system
(SLRC) to forecast chaotic time series and to successfully pass essential
benchmark tests, including a memory capacity test and a Mackey-Glass model
test.

##  1 Introduction

Reservoir computing (RC) [1] and its foundational concepts of context of
reverberated input histories [2], Liquid State Machine (LSM) [3] and Echo
State Network (ESN) [4] underpin an important class of machine learning (ML)
algorithms that can perform certain functions of a biological brain [3] and
that can be trained to recognise patterns and to forecast the response of
nonlinear dynamical systems that exhibit chaotic behaviour [5]. This
functionality enables using RC systems in many areas of science and
technology, including physics, psychology and finance [6, 7].

Figure 1: (a) ESN-based RC algorithm.(b) Sketch of the proposed SLRC. (c)
Photograph of SL waves of natural origin seen on the bowl of a drinking
fountain that spouts water.

The ability of an RC system to perform complex forecasting tasks originates
from its structural organisation [Fig. 1(a)] that separates a fast linear
readout layer from a much slower artificial neural network called the
reservoir [4]. The reservoir consists of a network of randomly connected
artificial neurons, and it converts a time-varying input into a spatio-
temporal pattern of neural activations that is then transformed into practical
outputs using the readout layer. While the readout layer is task-specific and
memory-less, the universal reservoir is based on the principles of nonlinear
dynamical systems and has certain memory [8]. Such a structure enables a
reservoir consisting of a relatively small number of artificial neurons to be
more efficient in resolving certain practical problems than a computer program
run on a supercomputer. Thus, one can draw an analogy between a reservoir and
some biological brains [3], including the brains of fruit flies and mosquitos
that, on average, have just about 200,000 neurons (in contrast to mammalian
brains that have tens of billions neurons) but enable the insects to perform
complicated tasks while they navigate and find food [9].

Although in the framework of LSM the analogy between a reservoir and a liquid
was made rather at an abstract level [3], follow-up research works have
demonstrated that physical liquids inherently possess all characteristics
required for the creation of a reservoir [10, 11, 12, 13, 14, 15, 16]. There
have also been theoretical proposals [17, 18, 19, 20] of RC systems based on
optical solitons [21]. However, there have been no experimental demonstrations
of soliton-based RC system yet.

In this paper, we experimentally demonstrate a physical RC system that
exploits the dynamics of solitary-like (SL) waves [Fig. 1(b)] that originate
from the spatio-temporal evolution of liquid films flowing over an inclined
surface [22, 23, 24]. We confirm the plausibility of the SL waves based RC
(SLRC) system using several standard test protocols designed to benchmark
computer programs that implement the concept of RC. SL waves are ubiquitous in
nature and can be readily observed in many everyday-life situations [Fig.
1(c)]. Subsequently, the principles of SLRC demonstrated in this work should
be applicable to numerous physical system that support roll waves [25, 24],
including microfluidic devices [26].

##  2 Algorithmic reservoir computing

A standard ESN algorithm [Fig. 1(a)] that we employ as a reference uses the
following nonlinear update equation:

$\displaystyle{\bf x}_{n}=(1-\alpha){\bf x}_{n-1}+\alpha\tanh({\bf W}^{in}{\bf
u}_{n}+{\bf W}{\bf x}_{n-1})\,,$ (1)

where $n$ is the index denoting equally-spaced discrete time instances
$t_{n}$, ${\bf u}_{n}$ is the vector of $N_{u}$ input values, ${\bf x}_{n}$ is
a vector of $N_{x}$ neural activations of the reservoir, the operator
$\tanh(\cdot)$ is a sigmoid activation function of a neuron, ${\bf W}^{in}$ is
the input matrix containing $N_{x}\times N_{u}$ elements, ${\bf W}$ is the
recurrent weight matrix consisting of $N_{x}\times N_{x}$ elements and
$\alpha\in(0,1]$ is a parameter that controls the reservoir’s temporal
dynamics.

One calculates the output weights ${\bf W}^{out}$ by solving a system of
linear equations ${\bf Y}^{target}={\bf W}^{out}{\bf X}$, where ${\bf X}$ and
${\bf Y}^{target}$ are the state matrix and the target matrix that are
constructed using, respectively, ${\bf x}_{n}$ and the vector of target
outputs ${\bf y}_{n}^{target}$ as columns for each discrete time instant
$t_{n}$. The solution is obtained in the form ${\bf W}^{out}={\bf
Y}^{target}{\bf X^{\top}}({\bf X}{\bf X^{\top}}+\beta{\bf I})^{-1}$, where
${\bf I}$ is the identity matrix, $\beta\geq 0$ is a regularisation
coefficient and ${\bf X^{\top}}$ is the transpose of ${\bf X}$ [4]. Then, one
solves Eq. (1) for new input data ${\bf u}_{n}$ and computes the output vector
${\bf y}_{n}={\bf W}^{out}[1;{\bf u}_{n};{\bf x}_{n}]$ using a constant bias
and the concatenation $[{\bf u}_{n};{\bf x}_{n}]$.

##  3 Solitary-like wave physical reservoir system

Figure 1(b) shows a sketch of the experimental setup, where SL waves
propagating over a metal plate inclined with respect to the ground by the
angle $\theta=3$ o [27] are used as a computational reservoir. The electronic
modulation of the pump flow rate is used to control the thickness of the
liquid film and to create SL waves on its surface. Adding an organic
fluorescent dye to the liquid (tap water) and irradiating it with UV light, we
generate green fluorescence light and measure variations of its intensity
caused by the SL waves. The optical intensity profiles are recorded using an
overhead high-speed (210 frames per second) digital camera and then processed
using customised software [27].

Figure 2: SL waves produced by a short training section of the randomly
generated binary input signal, where each ’1’ and ’0’ value is represented as
0.25-second-long square pulses and 0.25-second-long gaps between them,
respectively.

To validate SLRC, we adopt the approach used to test the ability of ESN
computer programs to predict chaotic time series (TS) [8, 5, 15]. First, the
response of the SL waves to a sufficiently long data set is recorded in one
single measurement. Then, the first half of the input data set is used for
training of the RC system but its second half is reserved for tests. Thus,
using the first half of the input data we calculate ${\bf W}^{out}$, but the
second half the data set, ${\bf u}_{n}$, produces the neural activations ${\bf
x}_{n}$ needed to compute the output as ${\bf y}_{n}={\bf W}^{out}[1;{\bf
u}_{n};{\bf x}_{n}]$.

As the first TS that is intended to test the memory capacity of SLRC, we use a
randomly generated vector of binary ’0’ or ’1’ values. This data vector is
converted into an electric signal using a digital function generator, where
0.25-second-long square pulses correspond to the ’1’ states but 0.25-second-
long spaces between the pulses represent the ’0’ states. This electric signal
modulates the pump flow [Fig. 1(b)]: the ’1’ states of the TS are implemented
as a brief change in the flow rate that results in the generation of one SL
wave but the ’0’ states correspond to no flow modulation. For example, an
input data subset ’10101’ is represented by three 0.25-second-long pulses
separated by two 0.25-second-long periods of no flow modulations, but
’111001111’ corresponds to one 0.75-second-long pulse separated by 0.5 seconds
of no flow modulation from the following 1-second-long pulse. The chosen
duration of both the pulses and the gaps between them is based on the optimal
SL wave generation parameters established in [27].

Figure 2 shows the SL waves produced by the training portion of the random TS.
While a periodic square pulse signal corresponding to a regular ’101010…’
binary sequence produces a pattern of equally spaced SL waves (not shown here,
see [27] for a relevant discussion), the SL wave profile in Fig. 2 follows the
random pattern of the input TS. The SL waves do not precisely reproduce the
training pulses, which is a typical picture seen at the training stage of
physical RC systems [28, 15]. As also established in algorithmic RC systems
[8], the response of the reservoir must not be a copy of the training signal
since in that case the trained RC system would not be able to produce an
output that resembles the target TS.

The SL wave signal in Fig. 2 corresponds to a single neural activation of the
reservoir. To generate more neural activations, we sample this signal with a
small discrete time step following the procedure described in [28], thus
producing 40 neural activations. We arrange the so-obtained activations into a
state matrix ${\bf X}$ and use it to calculate ${\bf W}^{out}$ following the
standard ESN algorithm [Fig. 1(a)].

As the second test, using the same training approach as above, we task the
SLRC to predict a Mackey-Glass TS (MGTS). MGTS is often employed to test RC
systems due to its nonlinear chaotic behaviour [29, 30, 15] and it is produced
by the delay differential equation [31]

$\displaystyle\dot{x}_{{}_{MG}}(t)$ $\displaystyle=$
$\displaystyle\beta_{{}_{MG}}\frac{x_{{}_{MG}}(\tau_{{}_{MG}}-t)}{1+x_{{}_{MG}}^{q}(\tau_{{}_{MG}}-t)}-\gamma_{{}_{MG}}x_{{}_{MG}}(t)\,,$
(2)

where we choose $\tau_{{}_{MG}}=17$ and set $q=10$, $\beta_{{}_{MG}}=0.2$ and
$\gamma_{{}_{MG}}=0.1$. The resulting TS modulates liquid flow in the
experiment. However, in contract to the square pulses, the MGTS modulation is
continuous. We note that the chosen set of MGTS parameters is often used to
test the performance of standard ML algorithms, including autoregressive
moving average (ARMA) [29]. Algorithmic RC systems are known to outperform
ARMA in forecasting chaotic TS [32]. Thus, we also use the MGTS test to verify
whether SLRC could outperform ARMA.

##  4 Benchmark test

Memory capacity. A reservoir suitable for applications in RC should have an
internal memory [8]. We test the memory capacity of SLRC using the $k$-delay
task applied to the random binary TS, where we calculate the correlation
coefficient $r({\bf y}^{target}_{n-k},{\bf y}_{k,n})$ between a delayed target
signal ${\bf y}^{target}_{n-k}$ and the predicted output ${\bf y}_{k,n}$
produced by the RC system trained using the training data ${\bf u}_{k,n}$ with
a discrete time delay $k$ [8]. The function $r^{2}({\bf y}^{target}_{n-k},{\bf
y}_{k,n})$ can accept any value from 0 to 1, being 0 an indicator of the full
loss of correlation and 1 corresponding a 100% correlation between the correct
target data and the prediction made by the RC system. Calculating $r^{2}$ as a
function of the discrete delay time $k$, we determine the memory capacity (MC)
as

$\displaystyle MC=\sum_{k=1}^{k_{max}}r^{2}({\bf y}^{target}_{n-k},{\bf
y}_{k,n})\,,$ (3)

being $k_{max}$ the maximum delay considered in the analysis.

Figure 3(a) plots $r^{2}$ for $k_{max}=20$ and produces $MC=13.4$ bits. We
observe that the last non-zero value of $r^{2}$ corresponds to $k=40$ thus
showing that SLRC, which has $N=40$ neurons, retains some memory of the past
inputs up to the theoretical fundamental limit $MC=N$ of algorithmic RC
systems [8]. Yet, in practice, in algorithmic RC systems $MC\approx 0.3N$ due
to an accumulated effect of rounding errors in the network update equations
[8]. Obtaining $N/MC\approx 0.35$ for SLRC, we conclude that imperfections of
the experimental setup and of the image processing technique used to detect SL
waves play the role of rounding error in an algorithmic RC system.

The value of MC reported in Fig. 3(a) is an order of magnitude higher than
that reported in [18], where an optical soliton RC system was theoretically
suggested. We hypothesise that the application of SL waves in the field of RC
since may be advantageous since they possess unique nonlinear properties that
are not available with solitary wave of other physical nature. Indeed, SL
waves can merge [22] instead of passing one through another, with the latter
being the case of well-known optical [21] and matter [33] soliton waves, which
we speculate may increase the memory of SLRC to the past inputs.

To support our arguments, we first refer to the literature on the memory
capacity of RC systems and we confirm that the behaviour of $r^{2}$ observed
in Fig. 3(a) closely resembles that of reservoirs consisting of artificial
neurons with highly nonlinear activation functions [34, 35]. We also task the
SLRC with a parity check test used to evaluate the performance of recurrent
neural network that operate in a nearly chaotic dynamic regime [36]. In this
test, the RC system predicts a TS defined as $PARITY({\bf u}_{n-k},{\bf
u}_{n-k-1},{\bf u}_{n-k-2})$ for increasing delays $k$. The so-generated TS is
a binary sequence itself and its forecast is a challenging task for any RC
system since the parity function is not linearly separable and its prediction
requires the reservoir to have memory. Similarly to Eq. (3), we calculate PC
memory capacity by summation of the parity check values up to the maximum
delay $k_{max}$.

Figure 3: Memory capacity of the SLRC determined using (a) $k$-delay and (b)
parity check tasks.

Figure 3(b) shows the outcome of the parity check test that produces $PC=12.4$
bits. In the literature on physical RC systems, the value of $PC$ is usually
slightly smaller than the value of $MC$ and, typically, it reaches 2 bits [37,
28, 38, 16]. The value of $PC$ for algorithmic RC systems is larger than that
of physical RC systems [36]. While for the SLRC the value of $PC$ is also
slightly smaller than that of $MC$, it is six times larger than the typical
value for physical RC systems used in this work as a reference.

A plausible explanation of this difference may be the fact that the phase
space of a standard algorithmic reservoir is chosen so that is does not
contain multiple periodic or chaotic attractors or fixed points, the presence
of which may disrupt the echo state property [30]. Physical RC systems are
also often designed according to this rule. However, in practice, many
physical RC systems with complex dynamics not only can reproduce the behaviour
of algorithmic RC systems having an echo state property but they can also
resolve complex nonlinear problems with a higher efficiency [36, 15]. Since SL
waves exhibit a complex nonlinear behaviour [22], similarly to the RC systems
discussed in [36, 15] SLRC should be capable of resolving some complex
nonlinear problems with high efficiency.

We also note that in Fig. 3(b) at the zero delay the parity check produces the
memory of approximately $0.96$ bit instead of the theoretically expected 1 bit
corresponding to 100% correlation. The exact origin of this discrepancy in our
proof-of-principle experiment is unknown. However, the following explanations
may be plausible. Firstly, it may be an artefact since it is not observed in
the MC test. Secondly, in algorithmic RC systems, this situation is known to
originate from long-term reverberations in the reservoir that deteriorate the
recall of immediately past inputs [2, 8]. In the SLRC, reverberations may be
caused by technical imperfections of the inlet fixture and by small
reflections of SL from the downstream edge of the inclined plate, which both
are well-established facts [27].

Figure 4: MGTS test of SLRC. (a) SL waves produced by a short training section
of the MGTS input signal. (b) One-timestep ahead prediction made by SLRC
(dots) compared with the expected target signal (line). Note the different
timescales use to plot the training and prediction stage data.

Mackey-Glass time series. We test the ability of SLRC to predict MGTS using a
one-step ahead prediction tasks that is often employed to test the accuracy of
neural network models [6, 39, 40]. In this task, the value of a TS at a time
step $t_{n}$ is mapped to a future time step $t_{n+1}$ and the resulting
prediction is compared with the expected target signal to calculate the
accuracy of the reservoir.

Using the approach outlined in [15], we resample the MGTS signal so that the
frequency of the principal peak in its Fourier spectrum equals 2 Hz, which is
the frequency of SL waves investigated in the relevant experimental work [27].
In Fig. 4(a) we plot a section of the training MGTS signal and the traces of
the SL waves created by it. As with the square pulses in Fig. 2, we observe an
SL wave pattern that follows but not precisely fit the training MGTS, which
indicates that the trained RC system will be able to generalise the learnt
input signal patterns [8].

In Fig. 4(b), we can see that the SLRC correctly predicts the target MGTS. We
calculate the normalised mean square error (NMSE), which is a standard measure
of the accuracy of algorithmic RC systems [4], as

$\displaystyle NMSE({\bf y},{\bf y}^{target})=\frac{\Bigl{\langle}\left\|{\bf
y}_{n}-{\bf
y}_{n}^{target}\right\|^{2}\Bigl{\rangle}}{\Bigl{\langle}\left\|{\bf
y}_{n}-\langle{\bf y}_{n}^{target}\rangle\right\|^{2}\Bigl{\rangle}}\,,$ (4)

where $\left\|\cdot\right\|$ is the norm and $\langle\cdot\rangle$ denotes the
averaging over time. An $NMSE=0$ means that the RC system produces a 100%
correct prediction but $NMSE=1$ means that the RC system has fully failed to
predict the target data. We obtain $NMSE_{SLRC}\approx 6.6\times 10^{-4}$.
Then, we task software that implements the ESN model Eq. (1) [30] to forecast
MGTS and obtain $NMSE_{ESN}\approx 1.75\times 10^{-7}$. A three orders of
magnitude difference between NMSE of SLRC and ESN is inconsequential for our
analysis and is attributed to a noise present in the SLRC measurements. We
also use software that implements an ARMA model [41] and predicts MGTS in a
one-step-ahead regime. We obtain $NMSE_{ARMA}\approx 1.5\times 10^{-3}$, which
shows that the accuracy of ARMA is lower than that of RC systems, which is an
established fact [32].

##  5 Theoretical analysis

It has been shown that considerable memory capacity of physical RC systems can
be achieved under a condition, when the system’s response to one pulse of a
chaotic TS does not fully decay before it is presented with the following
pulses, which enables the system to retain some traces of its response to
several consecutive pulses of the input TS [28]. SLRC has a similar memory
mechanism that is based on the property of SL waves of higher amplitude to
overtake and absorb SL waves of lower amplitude [22]. We theoretically
investigate this mechanism using a model, where, for the sake of illustration,
we excite SL waves using three square pulses.

We employ a simplified hydrodynamic Shkadov model that is based on the
boundary layer approximation of the Navier-Stokes equation under the
assumption of the parabolic longitudinal velocity profile [42, 43] and that is
formulated in terms of two coupled equations for the local film thickness
$h(x,t)$ and the fluid flux across the layer $q(x,t)$

$\displaystyle\rho\left[\partial_{t}q+\frac{6}{5}\partial_{x}\left(\frac{q^{2}}{h}\right)\right]$
$\displaystyle=$ $\displaystyle-\frac{3\mu q}{h^{2}}+\sigma
h\partial_{x}^{3}h$ $\displaystyle-$ $\displaystyle\rho
g\cos(\theta)h\partial_{x}h+\rho g\sin(\theta)h,$
$\displaystyle\partial_{t}h+\partial_{x}q$ $\displaystyle=$ $\displaystyle
f(x,t),$ (5)

where $\theta$ is the inclination angle, $\mu$ the dynamic viscosity and
$\sigma$ the liquid-air surface tension. The term $f(x,t)$ in the second
equation in Eqs. (5) is introduced to mimic the time-dependent influx of the
fluid induced by the inlet. We use a Gaussian
$f(x,t)=\exp{(-a(x-x_{0})^{2})}s(t)$ centred around the location of the inlet
$x_{0}$ with parameter $a$ chosen sufficiently large to ensure that the width
of the Gaussian is significantly smaller than the typical length of the
surface waves. The binary function $s(t)=f_{a}\\{0,1\\}$ is introduced to
describe the time modulation of the fluid influx and $f_{a}$ stands for the
amplitude of the influx. We numerically integrate Eqs. (5) in the $(-L\leq
x\leq L)$ domain of total length of $2L=400$ cm using periodic boundaries. The
inlet is located at $x_{0}=100$ cm. The relatively large domain length enables
us to observe the dynamics of the waves over a time period of up to 15
seconds, before the waves travel over the distance of $2L$.

To demonstrate the interaction and collision of waves with different
amplitudes, we excite three SL waves on a 0.7-mm-thick water film using three
flow modulation phases with a gradually increasing amplitude. All three phases
have a duration of 0.5 seconds and they are separated one from another by 0.5
seconds of no flow modulation [Fig. 5(a)]. The space-time plot of the excited
waves is shown in Fig. 5(b), where the arrows labelled as $(1,2,3)$ highlight
the propagation direction of the three SL waves. As anticipated, the first
smaller wave propagates slower than the second larger wave. These two waves
eventually collide approximately 100 cm further downstream and form a new wave
with a complex profile. Such an interaction scenario is known as an inelastic
wave collision [22]. All three waves $(1,2,3)$ eventually form one compact
structure that consists of many smaller waves around 150 cm further downstream
from the inlet location.

Figure 5: (a) Modulation function $s(t)$ used to excite three SL waves with
different amplitudes. (b) Space-time plot of the temporal evolution of the
three SL waves, labelled by the arrows $(1,2,3)$, excited by the modulation
function shown in panel (a). The first smaller wave $(1)$ propagates slower
than the second and the third larger waves $(2,3)$. The second wave
inelastically collides with the first wave about 100 cm further downstream.

Thus, we can see that, when the SL waves are detected very close to the inlet,
they do not have enough time to develop and interact, which, in turn, means
that the memory capacity of SLRC is low. The memory capacity will be low also
when SL waves are detected far downstream, where a significant transformation
of the waveform implies that the ability of the reservoir to memorise the
input signal is compromised. Subsequently, it is plausible that the optimal
memory capacity will be achieved when SL waves are detected midway the inlet
and the far downstream region where the SL waves merge, which is the approach
adopted in the experiments conducted in this work.

##  6 Discussion

Before we conclude, we discuss several potential improvements of SLRC that can
be done in the future. Firstly, the experimental setup used in this work holds
the potential to be extended to enable measurements, where the RC system could
generate a new TS using its own predictions as the input data. Such an
operating regime is called the generative mode and its implementation is a
challenging task for both physical and algorithmic RC systems [8, 30, 15]. In
particular, the operation in a generative mode requires introducing a feedback
circuit that converts the output of the readout layer into new input data at
every discrete time step. However, transfer of data from the readout layer,
which, typically, is a computer program itself, to the input layer of the RC
system introduces a time delay that can be several orders of magnitude larger
than the timescale of the reservoir dynamics. For instance, if the reservoir’s
natural frequency is of order of several GHz [28], all readout calculations
must be done at the GHz scale or above. In the case of SLRC, where the
frequencies of SL waves are of order of several Hz, this problem can be
resolved using a laser vibrometry technique to detect the SL waves [44] and
then processing the optical signal using ultrafast opto-electronic and
photonics systems that have already been adopted in the field of ML [45, 46].
The time delay introduced by such a processing scheme will be negligibly small
compared with the temporal dynamics of SL waves.

Apart from the operation in the generative mode, the introduction of a time
delay between the output and input of the reservoir may help increase the
memory capacity of the RC system [47, 28, 48], also enabling the reservoir to
operate using just a single artificial neuron [49]. Yet, in the field of fluid
dynamics, the introduction of a feedback control has been used to stabilise
liquid film flows [50], which may also help improve the performance of SLRC.

##  7 Conclusions

We have experimentally demonstrated a physical computational reservoir system
enabled by solitary-like waves propagating on the surface of a liquid thin
film flowing over an inclined surface. While computational reservoirs based on
optical and matter-wave solitons have been theoretically predicted but not
experimentally demonstrated yet, solitary wave dynamics of liquid films has
unique features that we have employed in the field of machine learning for the
first time. We trained and exploited the solitary-like wave reservoir to
perform several essential benchmark tests used to evaluate the performance of
machine learning algorithms designed for digital computer, showing that the
solitary-like wave reservoir has considerable memory capacity and can predict
chaotic time series more efficiently than certain standard algorithmic machine
learning models.

## References

  * [1] Nakajima K. Fisher I. Reservoir Computing (Springer, Berlin) 2021\.
  * [2] Kirby K. G. Proc. Florida AI Research Symposium (FLAIRS)199166.
  * [3] Maass W., Natschläger T. Markram H. Neural Comput.1420022531.
  * [4] Lukoševičius M. Jaeger H. Comput. Sci. Rev.32009127.
  * [5] Gauthier D. J., Bollt E., Griffith A. Barbosa W. A. S. Nat. Commun.1220215564.
  * [6] Marcellino M., Stock J. H. Watson M. W. J. Econom.1352006499.
  * [7] Small M. Applied Nonlinear Time Series Analysis: Applications in Physics, Physiology and Finance (World Scientific, Singapore) 2005.
  * [8] Jaeger H. A tutorial on training recurrent neural networks, covering BPPT, RTRL, EKF and the “echo state network” approach (GMD Report 159, German National Research Center for Information Technology) 2005.
  * [9] Raji J. I. Potter C. J. PLOS ONE1620211.
  * [10] Fernando C. Sojakka S. Pattern recognition in a bucket in proc. of Advances in Artificial Life, edited by Banzhaf W., Ziegler J., Christaller T., Dittrich P. Kim J. T. (Springer, Berlin) 2003 pp. 588–597.
  * [11] Jones B., Stekel D., Rowe J. Fernando C. Is there a liquid state machine in the bacterium escherichia coli? in proc. of 2007 IEEE Symposium on Artificial Life 2007 pp. 187–191.
  * [12] Nakajima K. Aoyagi T. IEICE Tech. Rep.1152015109.
  * [13] Nakajima K. Jpn. J. Appl. Phys.592020060501.
  * [14] Goto K., Nakajima K. Notsu H. N. J. Phys.232021063051.
  * [15] Maksymov I. S., Pototsky A. Suslov S. A. Phys. Rev. E1052021044206.
  * [16] Matsuo T., Sato D., Koh S.-G., Shima H., Naitoh Y., Akinaga H., Itoh T., Nokami T., Kobayashi M. Kinoshita K. ACS Appl. Mater. Interfaces14202236890.
  * [17] Jiang J. Lai Y.-C. Phys. Rev. Research12019033056.
  * [18] Zeng Q., Wu Z., Yue D., Tan X., Tao J. Xia G. Appl. Opt.5920206932.
  * [19] Marcucci G., Pierangeli D. Conti C. Phys. Rev. Lett.1252020093901.
  * [20] Silva N. A., Ferreira T. D. Guerreiro A. New J. Phys.232021023013.
  * [21] Kivshar Y. S. Agrawal G. P. Optical Solitons: From Fibers to Photonic Crystals (Academic Press, New York) 2003.
  * [22] Liu J. Gollub J. P. Phys. Fluids619941702.
  * [23] Chang H.-C. Annu. Rev. Fluid Mech.261994103.
  * [24] Kalliadasis S., Ruyer-Quil C., Scheid B. Velarde M. G. Falling Liquid Films (Springer-Verlag, London) 2012.
  * [25] Balmforth N. J. Mandre S. J. Fluid Mech.51420041.
  * [26] Hu X. Cubaud T. Phys. Rev. Lett.212018044502.
  * [27] Maksymov I. S. Pototsky A. Appl. Sci.1320221888.
  * [28] Watt S. Kostylev M. Phys. Rev. Appl.132020034057.
  * [29] Thissen U., van Brakel R., de Weijer A. P., Melssen W. J. Buydens L. M. C. Chemom. Intell. Lab. Syst.69200335.
  * [30] Lukoševičius M. A practical guide to applying echo state networks in Neural Networks: Tricks of the Trade, Reloaded, edited by Montavon G., Orr G. B. Müller K.-R. (Springer, Berlin) 2012 pp. 659–686.
  * [31] Mackey M. C. Glass L. Science1971977287.
  * [32] Maciel L., Gomide F., Santos D. Ballini R. Exchange rate forecasting using echo state networks for trading strategies in proc. of 2014 IEEE Conference on Computational Intelligence for Financial Engineering and Economics (CIFEr) 2014 pp. 40–47.
  * [33] Remoissenet M. Waves Called Solitons: Concepts and Experiments (Springer) 1994.
  * [34] Ganguli S., Huh D. Sompolinsky H. PNAS105200818970.
  * [35] Verzelli P., Alippi C. Livi L. Sci. Reps.9201913887.
  * [36] Bertschinger N. Natschläger T. Neural Comput.1620041413.
  * [37] Furuta T., Fujii K., Nakajima K., Tsunegi S., Kubota H., Suzuki Y. Miwa S. Phys. Rev. Appl.102018034063.
  * [38] Liao Z., Wang Z., Yamahara H. Tabata H. Chaos Solitons Fract.1532021111503.
  * [39] Rodan A. Tino P. IEEE Trans. Neural Netw.222011131.
  * [40] Morales G. B., Mirasso C. R. Soriano M. C. Neurocomputing4612021705.
  * [41] Khan S. Mackey Glass Time Series Prediction Using Least Mean Square (Matlab Central File Exchange) 2023.
  * [42] Shkadov V. Y. Izv. Akad. Nauk SSSR, Mekh. Zhidk. Gaza1196743.
  * [43] Shkadov V. Y. Izv. Akad. Nauk SSSR, Mekh. Zhidk. Gaza2196820.
  * [44] Maksymov I. Pototsky A. Phys. Rev. E1002019.
  * [45] Sorokina M. J. Phys. Photonics22020044006.
  * [46] Chembo Y. K. Chaos302020013111.
  * [47] Riou M., Torrejon J., Garitaine B., Araujo F. A., Bortolotti P., Cros V., Tsunegi S., Yakushiji K., Fukushima A., Kubota H., Yuasa S., Querlioz D., Stiles M. D. Grollier J. Phys. Rep. Appl.122019024049.
  * [48] Kondrashov A. V., Nikitin A. A., Nikitin A. A., Kostylev M. Ustinov A. B. J. Magn. Magn. Mater.5632022169968.
  * [49] Appeltant L., Soriano M. C., der Sande G. V., Danckaert J., Massar S., Dambre J., Schrauwen B., Mirasso C. R. Fischer I. Nat. Commun.22011468.
  * [50] Thompson A. B., Gomes S. N., Pavliotis G. A. Papageorgiou D. T. Phys. Fluids282016012107.