#  A Unifying Framework for Information Processing in Stochastically Driven
Dynamical Systems

Tomoyuki Kubota<EMAIL_ADDRESS>Graduate School of Information
Science and Technology, The University of Tokyo.  Hirokazu Takahashi  Graduate
School of Information Science and Technology, The University of Tokyo.  Kohei
Nakajima  Graduate School of Information Science and Technology, The
University of Tokyo.

###### Abstract

A dynamical system can be regarded as an information processing apparatus that
encodes input streams from the external environment to its state and processes
them through state transitions. The information processing capacity (IPC) is
an excellent tool that comprehensively evaluates these processed inputs,
providing details of unknown information processing in black box systems;
however, this measure can be applied to only time-invariant systems. This
paper extends the applicable range to time-variant systems and further reveals
that the IPC is equivalent to coefficients of polynomial chaos (PC) expansion
in more general dynamical systems. To achieve this objective, we tackle three
issues. First, we establish a connection between the IPC for time-invariant
systems and PC expansion, which is a type of polynomial expansion using
orthogonal functions of input history as bases. We prove that the IPC
corresponds to the squared norm of the coefficient vector of the basis in the
PC expansion. Second, we show that an input following an arbitrary
distribution can be used for the IPC, removing previous restrictions to
specific input distributions. Third, we extend the conventional orthogonal
bases to functions of both time and input history and propose the IPC for
time-variant systems. To show the significance of our approach, we demonstrate
that our measure can reveal information representations in not only machine
learning networks but also a real, cultured neural network. Our generalized
measure paves the way for unveiling the information processing capabilities of
a wide variety of physical dynamics which has been left behind in nature.

††preprint: APS/123-QED

##  I Introduction

Dynamical systems driven by external stimuli can be universally found in
nature, especially in biology. The dynamical aspects of information processing
found in biology have long been a source of inspiration for researchers who
wish to create a high-speed, energy efficient, and robust real-time
information processing device, which resolves a von Neumann bottleneck [1].
Reservoir computing (RC) [2, 3, 4] is a bioinspired information processing
paradigm that capitalizes on this dynamical perspective and has been widely
utilized in various fields in recent years [5, 6, 7, 8, 9, 10, 11, 12]. It
consists of a type of learning framework for recurrent neural networks (RNNs),
whose intermediate layer is referred to as the reservoir. In a reservoir
composed of $N$-nodes, the $i^{\rm th}~{}(i=1,2,\ldots,N)$ node state
$x_{i,t}$ at the $t^{\rm th}$ time step can be written as follows:

$\displaystyle x_{i,t+1}$ $\displaystyle=$ $\displaystyle
f\left(\sum_{j=1}^{N}w_{ij}x_{j,t}+w_{in,i}u_{t}\right),$ (1)

where $f$ is the activation function, and $w_{ij}$ and $w_{in,i}$ are the
internal and input weights, respectively. To emulate the target output
$z_{t}$, we use linear regression to obtain an estimate of $z_{t}$,
$\tilde{z}_{t}$, as follows:

$\displaystyle\tilde{z}_{t}$ $\displaystyle=$
$\displaystyle\tilde{\bm{w}}_{out}^{\top}\cdot{\bm{x}}_{t},$ (2)
$\displaystyle\tilde{\bm{w}}_{out}$ $\displaystyle=$
$\displaystyle\arg\min_{{\bm{w}}_{out}}\sum_{t=1}^{T}\left(z_{t}-{\bm{w}}_{out}^{\top}\cdot{\bm{x}}_{t}\right)^{2},$
(3)

where ${\bm{w}}_{out}$ and $\tilde{\bm{w}}_{out}\in\mathbb{R}^{N}$ are the
weight and solution vector for the target output, respectively. This learning
method leverages the dynamical resource through training without affecting the
state of the reservoir but places a constraint on ${\bm{x}}_{t}$. To perform
reproducible computation, RC requires ${\bm{x}}_{t}$ to be the same response
against identical input time-series (i.e., the state needs to be an echo
function, which is a function of only the past input time-series
${\bm{u}}_{t}=\\{u_{t-s}\\}_{s=1}^{t}$). This dynamical property is referred
to as the echo state property (ESP) [13, 14, 15] or the fading memory property
(FMP) [3, 16, 17], which are slightly different from each other (see the
Appendix for further details). According to these properties, various
activation functions can be used for the reservoir node. Furthermore, as the
reservoir is not limited to a computer-generated system, it can be replaced
with a real physical system. A reservoir using such a physical system is
called a physical reservoir [18]. Some of the above systems have a wide range
of dynamics that are not readily found in conventional neural networks.

In the literature, a measure called information processing capacity (IPC) [19]
has been proposed to quantify the information processing capability of
dynamical systems that have ESPs or FMPs. The IPC measures the type and
quantity of input history that is handled and held in the system by
decomposing the system state into an orthogonal basis [20, 21, 22]. Using the
$N$-dimensional state ${\bm{x}}_{t}\in\mathbb{R}^{N}$ and the one-dimensional
stochastic input $\zeta_{t}\in\mathbb{R}$ at the $t^{\rm
th}~{}(t\in\mathbb{Z})$ time step, the input-driven dynamical system (IDS)
determines the next state, as shown below:

$\displaystyle{\bm{x}}_{t+1}={\bm{f}}\left({\bm{x}}_{t},\zeta_{t}\right),$ (4)

where ${\bm{f}}$ maps
$\mathbb{R}^{N}\times\mathbb{R}\rightarrow\mathbb{R}^{N}$. The IPC evaluates
the emulation ability of the $i^{\rm th}$ target output
$z^{(i)}_{t}\in\mathbb{R}~{}(i=1,2,\ldots)$ from
${\bm{x}}_{t}~{}(t=1,2,\ldots,T)$. $z^{(i)}_{t}~{}(i=1,2,\ldots)$ is
represented by the product of the $n_{s}^{(i)}$-th order polynomial of the
random variable delayed by $s$ steps, $\zeta_{t-s}$,

$\displaystyle z^{(i)}_{t}$ $\displaystyle=$
$\displaystyle\prod_{s=1}^{\infty}\mathcal{F}_{n_{s}^{(i)}}(\zeta_{t-s}),$ (5)

where $\mathcal{F}_{n}(\zeta)$ represents the $n^{\rm th}$ order orthogonal
polynomial of $\zeta$. From Eqs. (2) and (3), we obtain an estimate of
$z^{(i)}_{t}$, $\tilde{z}^{(i)}_{t}$. When $z_{t}$ is an orthogonal function
of the independent variables $\\{\zeta_{t-s}\\}_{s=1}^{\infty}$, the IPC is
defined using a normalized emulation error of the reservoir, as follows:

$\displaystyle C\left({\bm{X}},{\bm{z}}^{(i)}\right)$ $\displaystyle=$
$\displaystyle
1-\frac{\min_{\bm{w}}\sum_{t=1}^{T}\left(z^{(i)}_{t}-{\bm{w}}^{\top}{\bm{x}}_{t}\right)^{2}}{\sum_{t=1}^{T}\left(z^{(i)}_{t}\right)^{2}}$
(6) $\displaystyle=$
$\displaystyle\frac{{\bm{z}}^{(i)\top}{\bm{X}}\left({\bm{X}}^{\top}{\bm{X}}\right)^{-1}{\bm{X}}^{\top}{\bm{z}}^{(i)}}{{\bm{z}}^{(i)\top}{\bm{z}}^{(i)}},$

where ${\bm{w}}\in\mathbb{R}^{N}$,
${\bm{X}}=\left[{\bm{x}}_{1}\cdots{\bm{x}}_{T}\right]^{\top}\in\mathbb{R}^{T\times
N}$, and ${\bm{z}}^{(i)}=\left[z^{(i)}_{1}\cdots
z^{(i)}_{T}\right]^{\top}\in\mathbb{R}^{T}$ are the weight vector, state, and
target output, respectively. In this case, the uniform random variable and the
Legendre polynomial are the stochastic variable $\zeta$ and orthogonal
polynomial $\mathcal{F}_{n}(\zeta)$, respectively, although the combination is
not restricted. For example, a Gaussian random variable and a Hermite
polynomial are also suitable [19]. Therefore, the IPC is a measure used to
evaluate the input information held by the state with the emulation ability of
the orthogonal basis.

In this connection, there is a theory about a deterministic dynamical system
with a stochastic input $\zeta_{t}$ in a different context. The system can be
described as an operator of $\\{\zeta_{t-s}\\}_{s=1}^{\infty}$ according to
the polynomial chaos expansion [23], which is a series expansion using the
target outputs of IPC as the bases (i.e., multivariate orthogonal polynomials
of the random variables, $z_{t}^{(i)}$, described in Eq. [5]). The polynomial
chaos expansion has been frequently utilized to determine the evolution of
uncertainty in a dynamical system when there is probabilistic uncertainty in
the system parameters [24, 25]. These multivariate polynomials are referred to
as polynomial chaoses (PCs), and the space spanned by PCs is called
homogeneous chaos [23], expressed as

$\displaystyle{\bm{x}}_{t}=\sum_{i=1}^{\infty}{\bm{c}}_{i}z_{t}^{(i)},$ (7)

where ${\bm{c}}_{i}\in\mathbb{R}^{N}$ is the $i^{\rm th}$ coefficient vector.
If the input $\zeta_{t}$ follows a certain distribution, the PC is determined
based on its orthogonality. Let ${\bm{\zeta}}$ be a vector notation of $\zeta$
sampled $T$ times ${\bm{\zeta}}=\\{\zeta_{1},\ldots,\zeta_{T}\\}$. If a
weighting function $w({\bm{\zeta}})$ specific to the PCs
$\\{z_{t}^{(i)}\\}_{i=1}^{\infty}$ exists, the following orthogonality
relations should be satisfied:

$\displaystyle\left<z_{t}^{(i)}z_{t}^{(j)}\right>$ $\displaystyle=$
$\displaystyle\left<(z_{t}^{(i)})^{2}\right>\delta_{ij},$ (8)
$\displaystyle\left<f({\bm{\zeta}})g({\bm{\zeta}})\right>$ $\displaystyle=$
$\displaystyle\sum_{\bm{\zeta}}w({\bm{\zeta}})f({\bm{\zeta}})g({\bm{\zeta}}),$
(9)

where $\delta_{ij}$ is the Kronecker delta function. The PC expansion can use
various combinations of input and polynomials. The generalized polynomial
chaos (gPC)[24] supplies the PC for specific combinations—e.g., the Hermite
polynomial for Gaussian distributions, the Legendre polynomial for uniform
distributions, and the Charlier polynomial for Poisson distributions are
available as $\mathcal{F}_{n}(\zeta)$. Furthermore, the arbitrary polynomial
chaos (aPC) [25] is the PC for random variables following an arbitrary
probability distribution by using the Gram-Schmidt orthogonalization procedure
(See the Appendix for both schemes).

As previously described, the IPC and PC expansion have a number of
similarities and differences. Both schemes use the IDS with stochastic inputs
and multivariate orthogonal polynomials, whereas the types of input
distributions and orthogonal polynomials are not as limited in the IPC as they
are in the PC. One objective of this paper is to establish a connection
between the IPC and the IDS expanded with PCs. Thus, we first aim to reveal
this relationship by deriving the IPC from the state expanded with PCs.
Second, to enlarge the applicable range of IPC, we extend this relationship
for time-variant systems. So far, the IPC assumed that the system is a
function of only the input time-series. However, a solution of the dynamical
system in Eq. (4) is a function of the input time-series and time with a given
initial state. Introducing time-dependent orthogonal bases to PCs, we aim to
derive the IPC for time-variant systems and illustrate that input information
processing is performed by coupling the terms of time and input time-series.
Finally, to demonstrate the potential of our approach, we apply our theory to
three cases—a model that is frequently used as a benchmark task in the context
of temporal machine learning, an artificial neural network, and a real,
cultured neural network—to reveal their information processing within the
systems.

##  II Methods

###  II.1 Singular value decomposition

Singular value decomposition (SVD) breaks down state
$\bm{X}\in\mathbb{R}^{T\times N}$ into

$\displaystyle{\bm{X}}={\bm{P}}{\bm{\Sigma}}{\bm{Q}}^{\top},$ (10)

where ${\bm{P}}\in\mathbb{R}^{T\times r}$ and ${\bm{Q}}\in\mathbb{R}^{N\times
r}$ are matrices whose column and row vectors, respectively, are singular
vectors, ${\bm{\Sigma}}={\rm
diag}\\{\sigma_{1},\ldots,\sigma_{r}\\}\in\mathbb{R}^{r\times r}$ is a
diagonal matrix containing the singular values $\sigma_{i}~{}(i=1,\ldots,r)$,
and $r(\leq N<T)$ is the rank of ${\bm{X}}$.

###  II.2 ESN

Let the $i^{\rm th}~{}(i=1,2,\ldots,N)$ state of the ESN at the $t^{\rm th}$
step be $x_{i,t}$. The state equation is given as follows:

$\displaystyle x_{i,t+1}$ $\displaystyle=$ $\displaystyle
f\left(\rho\sum_{j=1}^{N}w_{i,j}x_{j,t}+\iota\omega_{i}\zeta_{t}\right),$ (11)

where $w_{ij}$ was initialized with the uniform random number in the range of
$[-1,1]$ and multiplied by a constant so that the maximum eigenvalue of the
matrix ${\bm{W}}=[w_{ij}]$ was 1. The input weight $\omega_{i}$ was also set
to a uniform random number in the range of $[-1,1]$, and $\rho$,
$\iota(=0.1)$, and $N(=50)$ represent the spectral radius of $\rho{\bm{W}}$,
input intensity, and the number of nodes, respectively. $f$ denotes the
activation function, three types of which were used to solve the NARMA10 task
[26, 4]: a linear function

$\displaystyle f(y)=y,$ (12)

a hyperbolic tangent function

$\displaystyle f(y)=\tanh(y),$ (13)

and an analog integrator function

$\displaystyle
f(y)=\left(1-\frac{1}{\tau}\right)x_{i,t}+\frac{1}{\tau}\tanh(y),$ (14)

where $\tau$ was set to 1.25.

###  II.3 One-dimensional ESN

The one-dimensional ESN is described as

$\displaystyle x_{t+1}$ $\displaystyle=$ $\displaystyle\tanh(\rho
x_{t}+u_{t}),$ (15) $\displaystyle u_{t}$ $\displaystyle=$
$\displaystyle\mu+\sigma\zeta_{t},$ (16)

where $x_{t}$, $u_{t}$, and $\zeta_{t}$ are the state, input, and random
variable at the $t^{\rm th}$ step, respectively, and $\rho$ and $\sigma$ are
an inner weight and input intensity, respectively. To treat the bounded and
time-invariant state, we chose $\rho=0.95$ and $\mu=0$.

###  II.4 NARMA10 benchmark task

Non-linear autoregressive moving average (NARMA) tasks were introduced to test
the performance of RNNs [27]. Within these tasks, the target output $y_{t}$ is
generated with the NARMA model
$y_{t+1}=g(y_{t},y_{t-1},\ldots,\zeta_{t},\zeta_{t-1},\ldots)$ and input
series $\\{\zeta_{t}\\}_{t=1}^{T}$. The RNN receives the same input
$\zeta_{t}$ and emulates $y_{t}$ by modifying its weights. NARMA tasks,
including NARMA2, NARMA10 [27], and NARMA30 [28], evaluate the ability to
emulate the model. In particular, NARMA10 has been widely used as a benchmark
task for RC [26, 29, 4, 30, 5, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42,
43, 44, 45]. Some specific examples include testing the performance of ESNs
[26, 29, 4] as well as evaluating the computational capability of physical
systems, such as quantum systems [30, 36, 46], analog circuits [5], opto-
electronic architectures [31], carbon nanotubes [32], soft robotic systems
[33, 34, 47], and other dynamical systems [35, 36, 37, 38, 40, 39]. Thus,
NARMA10 is one of the most representative benchmark tasks and is broadly
utilized to compare and evaluate the computational capabilities of dynamical
systems.

Let the state and input at the $t^{\rm th}$ step be $y_{t}$ and $u_{t}$,
respectively. The NARMA10 model is given by

$\displaystyle y_{t+1}$ $\displaystyle=$ $\displaystyle\alpha y_{t}+\beta
y_{t}\sum_{s=0}^{9}y_{t-s}+\gamma u_{t}u_{t-9}+\delta,$ (17) $\displaystyle
u_{t}$ $\displaystyle=$ $\displaystyle\mu+\kappa\zeta_{t},$ (18)

where the default constant parameters $(\alpha,\beta,\gamma,\delta)$ are set
to $(0.3,\ 0.05,\ 1.5,\ 0.1)$, $\zeta_{t}$ is the random variable at the
$t^{\rm th}$ step and follows a uniform distribution in the $[-1,1]$ interval,
while $\mu$ and $\kappa$ are the average of $\zeta_{t}$ and the input
intensity parameter, respectively. This paper uses two ranges:
$u_{t}\in[0,\sigma]$ ($\mu=\kappa=\sigma/2$) and $u_{t}\in[-\sigma,\sigma]$
($\mu=0$, $\kappa=\sigma$). All the initial values of
$y_{t}~{}(t=0,1,\ldots,9)$ were set to zero, except when the basin of
attraction was examined.

###  II.5 The limit cycle system

The simple limit cycle system with radius $r$ and azimuth $\theta$ in polar
coordinates [48] is discretized as follows:

$\displaystyle r_{t+1}$ $\displaystyle=$ $\displaystyle(1+\tau)r_{t}-\tau
r_{t}^{3}+\tau u_{t},$ (19) $\displaystyle\theta_{t+1}$ $\displaystyle=$
$\displaystyle\theta_{t}+\tau\omega,$ (20) $\displaystyle u_{t}$
$\displaystyle=$ $\displaystyle\mu+\sigma\zeta_{t},$ (21)

where $\omega=2\pi/3$ and $\tau=0.1$ are the angular velocity and time step
width, respectively, and $\zeta_{t}$ is the uniform random number in the
$[-1,1]$ interval. Therefore, the input $u_{t}$ follows the uniform
distribution in the $[\mu-\sigma,\mu+\sigma]$ range ($\mu=0.2$ and
$\sigma=1.5$) and is applied in the radial direction. The Cartesian
coordinates are given by $x_{t}=r_{t}\cos\theta_{t}$ and
$y_{t}=r_{t}\sin\theta_{t}$.

##  III Results

###  III.1 The equivalence of the IPC and coefficient in PC expansion

To show the relationship between the IPC and PC expansion, we derive the IPC
from the state expanded in terms of PCs. First, we transform the IPC into a
simpler form using SVD, which reduces the $N$ state time-series to $r(\leq N)$
linearly independent time-series vectors
${\bm{p}}_{j}\in\mathbb{R}^{T}~{}(j=1,\ldots,r)$. Using the decomposed state,
we can rewrite the IPC relative to ${\bm{z}}^{(i)}$ as

$\displaystyle
C\left({\bm{X}},{\bm{z}}^{(i)}\right)=\sum_{j=1}^{r}\left({\bm{p}}_{j}^{\top}{\bm{\phi}}^{(i)}\right)^{2},$
(22)

where ${\bm{\phi}}^{(i)}={\bm{z}}^{(i)}/||{\bm{z}}^{(i)}||$ is the normalized
output. Next, assuming that ${\bm{x}}_{t}$ can be expanded with PCs, the state
${\bm{X}}$ is described as

$\displaystyle{\bm{X}}=\sum_{i=1}^{\infty}{\bm{\phi}}^{(i)}\cdot\hat{\bm{c}}_{i}^{\top}={\bm{\Phi}}\hat{\bm{C}}^{\top},$
(23)

where ${\bm{\Phi}}=\left[{\bm{\phi}}^{(1)}~{}{\bm{\phi}}^{(2)}\cdots\right]$
and $\hat{\bm{C}}=[\hat{\bm{c}}_{1}~{}\hat{\bm{c}}_{2}\cdots]^{\top}$ are the
basis matrix and coefficient matrix, respectively. Comparing the decomposed
state ${\bm{X}}={\bm{P}}{\bm{\Sigma}}{\bm{Q}}^{\top}$ (see Methods) and Eq.
(23), we obtain the matrix

$\displaystyle{\bm{P}}$ $\displaystyle=$
$\displaystyle{\bm{\Phi}}{\bm{\Lambda}}^{\top},$ (24)

where ${\bm{P}}=[{\bm{p}}_{1}~{}{\bm{p}}_{2}\cdots{\bm{p}}_{r}]^{\top}$ and
${\bm{\Lambda}}={\bm{\Sigma}}^{-1}{\bm{Q}}^{\top}\hat{\bm{C}}=[{\bm{\lambda}}_{1}~{}{\bm{\lambda}}_{2}\cdots]$
are a matrix form of linearly independent vectors and a constant matrix,
respectively, while ${\bm{\lambda}}_{i}\in\mathbb{R}^{r}$. From Eqs. (22) and
(24), the IPC becomes

$\displaystyle C({\bm{X}},{\bm{z}}^{(i)})=||{\bm{\lambda}}_{i}||^{2}.$ (25)

Because ${\bm{P}}$ is also written as
${\bm{P}}=\sum_{i=1}^{\infty}{\bm{\phi}}^{(i)}\cdot{\bm{\lambda}}_{i}^{\top}$,
Eq. (25) illustrates that the computation of the $i^{\rm th}$ IPC is
equivalent to expanding the state with the PCs and calculating the squared
norm of the coefficient of the $i^{\rm th}$ PC of the temporal bases
${\bm{p}}_{j}$ expanded by the PCs.

In addition, the IPC has an important property of summation. The total
capacity is described as

$\displaystyle
C_{tot}=\sum_{i=1}^{\infty}C\left({\bm{X}},{\bm{z}}^{(i)}\right)=\sum_{j=1}^{r}||{\bm{\Phi}}^{\top}{\bm{p}}_{j}||^{2}.$
(26)

Eq. (26) yields the sum of the squared norm of ${\bm{p}}_{j}~{}(j=1,\ldots,r)$
projected into an infinite-dimensional space that contains the orthogonal
vectors ${\bm{\phi}}^{(i)}~{}(i=1,2,\ldots)$; therefore, if all of the
${\bm{p}}_{j}$ are functions of only past input history,
$||{\bm{\Phi}}^{\top}{\bm{p}}_{j}||$ becomes one, resulting in

$\displaystyle C_{tot}=r~{}(\leq N).$ (27)

We call this condition an integrity property. Under the assumption that the
state is a function of only past input history, IPCs hold the integrity
property in information processing. These results provide a new perspective
that the coefficient of PC expansion represents the amount of processed input.

###  III.2 Demonstration of comprehensive computational capabilities using
general input distribution

To illustrate that various orthogonal polynomials can be used as target
outputs for the IPC, we examine the total capacities for eight target gPCs and
four target aPCs using one-dimensional ESNs, which is a time-invariant system.
FIG. 1(a) shows the IPC breakdowns with eight types of gPCs in the Askey
scheme. The utilized gPCs of random distributions are the Hermite-chaos of a
Gaussian distribution, Laguerre-chaos of a gamma distribution, Jacobi-chaos of
a beta distribution, Legendre-chaos of a uniform distribution, Charlier-chaos
of a Poisson distribution, Krawtchouk-chaos of a binomial distribution,
Meixner-chaos of a negative binomial distribution, and Hahn-chaos of a
hypergeometric distribution. The total IPCs are all one, suggesting that the
gPCs form a complete orthogonal system with any combination of distributions,
and orthogonal polynomials and can be used for the IPC.

Furthermore, to demonstrate that aPCs can also be used as the target outputs
for the IPC, we estimated the IPCs with the ESN given four types of inputs
following a mixed Gaussian distribution, Pareto distribution, Zipf
distribution, and Bernoulli distribution, which do not follow the Askey
scheme. To investigate the IPC of the ESN, we built Gram-Schmidt-chaoses. As
shown in FIG. 1(a), the total IPCs were one, indicating that the aPCs formed a
complete orthogonal system with the input distribution and Gram-Schmidt-
chaoses and can be used for the IPC.

These results suggest that both the gPC and aPC are suitable as the target
output of the IPC.

Figure 1:  Demonstration of theories. (a) IPC breakdowns of the one-
dimensional ESN with eight types of target generalized/arbitrary PCs of a
random variable $\zeta_{t}$: Hermite-chaos with Gaussian random variables,
Laguerre-chaos with gamma random variables, Jacobi-chaos with beta random
variables, Legendre-chaos with uniform random variables, Charlier-chaos with
Poisson random variables, Krawtchouk-chaos with binomial random variables,
Meixner-chaos with negative binomial random variables, Hahn-chaos with
hypergeometric random variables, and Gram-Schmidt-chaos with mixed-Gauss
random variables, Pareto random variables, Zipf random variables, and
Bernoulli random variables. The $n^{\rm th}$-order capacities $(n=1,\ldots,8)$
in the total capacity $C_{tot}$ are summarized. The variable $p$ represents
the probability distribution of $\zeta_{t}$. (b) The TIPC of the limit cycle
system in Eqs. (19)–(21). The phase plane of $x_{t}$ and $y_{t}$ in the system
(upper panel). The systems without input ($\mu=\sigma=0$; red line) and with
input ($\mu=0.2,\sigma=1.5$; black line) are shown. The first-order capacities
of $P_{1}(\zeta_{t-s_{1}})\cos(\omega\tau t+\alpha_{1,s_{1}})$ with the delay
step $s_{1}$ (middle panel). The second-order capacities relative to delay
steps $s_{1}$ and $s_{2}$ (lower panel). The diagonal ($s_{1}=s_{2}$) and the
upper left triangle ($s_{1}<s_{2}$) dots show the capacities of
$P_{2}(\zeta_{t-s_{1}})\cos(\omega\tau t+\alpha_{2,s_{1}})$ and
$P_{1}(\zeta_{t-s_{1}})P_{1}(\zeta_{t-s_{2}})\cos(\omega\tau
t+\alpha_{2,s_{1},s_{2}})$, respectively. In the white area, the TIPC
calculation is omitted.

###  III.3 Extending the integrity in information processing to a time-variant
domain

The information processing integrity shown so far holds only if the system is
a function of a past input series. Using the derived relationship between the
IPC and PC expansion, we extend the application range of IPC to time-varying
IDSs. First, we introduce a classification of IDSs based on the conditions
imposed on the solution of Eq. (4). From the connection between the IPC and
PC, the IPC is calculated by extracting the second-order process
${\bm{x}}_{t}$ in the $1\leq t\leq T$ range from the original sequence, which
is obtained according to Eq. (4). By recursively using Eq. (4) from $t=0$, the
solution is clearly determined from the time $t$, input sequence
${\bm{\zeta}}_{t}=\\{\zeta_{t-s}\\}_{s=1}^{t}$, and initial state
${\bm{x}}_{0}$, as follows:

$\displaystyle{\bm{x}}_{t}={\bm{\xi}}(t,{\bm{\zeta}}_{t};{\bm{x}}_{0}),$ (28)

where ${\bm{\xi}}$ is determined by $\bm{f}$ in Eq. (4). However, the state is
described as an operator of only $\\{\zeta_{t-s}\\}$, showing that the PC
expansion assumes that ${\bm{x}}_{t}$ is time-invariant. Therefore, two
conditions are imposed on the time-series. First, every $i^{\rm
th}~{}(i=1,\ldots,N)$ state time-series $x_{i,t}~{}(t=1,\ldots,T)$ must be a
second-order process. The specific condition is that the second moment of
$x_{i,t}$ must be finite:

$\displaystyle\left<x_{i,t}^{2}\right><\infty~{}(i=1,\ldots,N),$ (29)

where $x_{i,t}$ is the $i^{\rm th}$ element of ${\bm{x}}_{t}$,
$\left<x_{i,t}^{2}\right>=\sum_{t=1}^{T}(x_{i,t}-\bar{x}_{i})^{2}/T$, and
$\bar{x}_{i}=\sum_{t=1}^{T}x_{i,t}/T$. Finiteness is a prerequisite for
expansion convergence in the sense of $L_{2}$ due to the Cameron-Martin
theorem. Second, the solution ${\bm{\xi}}$ should be time-invariant. The
extracted time-series needs to be described only with the input time-series
${\bm{\zeta}}_{t}$. Thus,
${\bm{x}}_{t}={\bm{\xi}}(t,{\bm{\zeta}}_{t};{\bm{x}}_{0})$ is time-invariant
if

$\displaystyle{\bm{\xi}}(t-\tau,{\bm{\zeta}}_{t-\tau};{\bm{x}}_{0})$
$\displaystyle=$
$\displaystyle{\bm{\xi}}(t,{\bm{\zeta}}_{t-\tau};{\bm{x}}_{0})~{}{\rm
for}~{}\tau\in\mathbb{Z}.$ (30)

Therefore, the PC expansion assumes the non-divergence and time-invariance of
the system.

In contrast, if ${\bm{x}}_{t}$ is time-variant, the homogeneous chaos is no
longer a complete orthogonal system, and ${\bm{x}}_{t}$ is represented by the
time-dependent polynomial chaos (TDPC) as

$\displaystyle{\bm{x}}_{t}=\sum_{i=1}^{\infty}{\bm{\gamma}}_{i}({\bm{x}}_{0})\psi_{t}^{(i)}\phi_{t}^{(i)},$
(31)

where ${\bm{\gamma}}_{i}({\bm{x}}_{0})\in\mathbb{R}^{N}$, $\phi_{t}^{(i)}$,
and $\psi_{t}^{(i)}$ are the $i^{\rm th}$ coefficient vector depending on
${\bm{x}}_{0}$, polynomial chaos, and the time-dependent basis, respectively.
We can construct a complete orthogonal system from $t$ and
$\\{\zeta_{t-s}\\}_{s=1}^{\infty}$ by adding bases of $t$ that are orthogonal
to the PCs of $\\{\zeta_{t-s}\\}$. In the following discussion, we assume that
${\bm{x}}_{0}$ is given and fixed. By converting the summation over the input
time-series ${\bm{\zeta}}_{t}$ into a summation over time $t$, Eq. (8) can be
rewritten as follows:

$\displaystyle\left<f_{t}^{(i)}f_{t}^{(j)}\right>$ $\displaystyle=$
$\displaystyle\left<(f_{t}^{(i)})^{2}\right>\delta_{ij},$ (32)
$\displaystyle\left<f_{t}^{(i)}f_{t}^{(j)}\right>$ $\displaystyle=$
$\displaystyle\sum_{t=1}^{T}f_{t}^{(i)}f_{t}^{(j)},$ (33)

and can be satisfied by the PCs. Next, we define the time-dependent basis that
satisfies this inner product. For example, $\\{1,\cos\omega t,\sin\omega
t,\cos 2\omega t,\sin 2\omega t,\ldots\\}$ ($\omega=2\pi/T$) is used for the
Fourier series expansion and is a complete orthogonal basis of time. Since
$\phi_{t}^{(i)}$ and $\psi_{t}^{(j)}$ are uncorrelated with each other, we
expect that the following orthogonality relation will be satisfied for a
sufficiently long period $T$:

$\displaystyle\left<\phi_{t}^{(i)}\psi_{t}^{(i)}\phi_{t}^{(j)}\psi_{t}^{(j)}\right>$
$\displaystyle=$
$\displaystyle\left<\left(\phi_{t}^{(i)}\psi_{t}^{(i)}\right)^{2}\right>\delta_{ij}.$
(34)

In this paper, we call the time and input time-series-dependent basis in Eq.
(31), $\psi_{t}^{(i)}\phi_{t}^{(i)}$, the TDPC and define the space spanned by
the time-dependent polynomial chaoses as time-dependent homogeneous chaos
(TDHC).

The IPCs of ${\bm{x}}_{t}$ in the TDHC can be estimated by replacing the
target output ${\phi}_{t}^{(i)}$ in Eq. (22) with the TDPC
$\psi_{t}^{(i)}\phi_{t}^{(i)}$. We define the IPC with the target TDPCs as
temporal information processing capacity (TIPC). As with the IPC of a time-
invariant system, the TIPC of a time-variant system is equivalent to expanding
the state with TDPCs and calculating the squared norm of the coefficient of
the $i^{\rm th}$ TDPC of the temporal bases ${\bm{p}}_{j}~{}(j=1,\ldots,r)$
expanded by the TDPCs. Note that to remove the terms that do not include the
input, the time-average of ${\bm{x}}_{t}$ and the time-varying terms—for
example, $A_{n}\cos(2\pi f_{n}t+\theta_{n})~{}(n=1,2,\ldots)$, where the
amplitude $A_{n}$, the frequency $f_{n}$, and the phase $\theta_{n}$ are
estimated by the Fourier transform—are subtracted. Since TDPCs constitute a
complete orthogonal system, the norm of each of the $r$-time-series vectors in
the space is one, and the total TIPC retains $r$. This extension can reveal
aspects of the information processing performed by time-variant systems.

###  III.4 Demonstration of complete computational capabilities in a time-
variant system

To demonstrate the extension of the IPC, we show the TIPC of a two-dimensional
limit cycle system. As shown in FIG. 1(b), input that follows a uniform
distribution forces the system to fluctuate around the non-input state. Using
this time-variant system, we calculated the TIPCs, whose TDPCs were
constructed using the PC $\\{\phi_{t}\\}$ for $n_{s}~{}(\sum_{s}n_{s}<5)$,
$s<10$, and temporal basis $\psi_{t}\in\\{\cos\Omega t,\sin\Omega t,\cos
2\Omega t,\sin 2\Omega t,\ldots,\cos\frac{T}{2}\Omega t,\sin\frac{T}{2}\Omega
t\\}$, where $\Omega=2\pi/T$. From the estimated TIPCs, the state on the two-
dimensional plane ${\bm{x}}_{t}=[x_{t}~{}y_{t}]^{\top}$ can be expressed as
follows:

$\displaystyle{\bm{x}}_{t}$ $\displaystyle=$
$\displaystyle\sum_{s_{1}=1}^{18}{\bm{p}}_{s_{1}}P_{1}(\zeta_{t-s_{1}})\cos\left(\omega\tau
t+\alpha_{1,s_{1}}\right)$ (35) $\displaystyle+$
$\displaystyle\sum_{s_{1},s_{2}}{\bm{q}}_{s_{1},s_{2}}P_{1}(\zeta_{t-s_{1}})P_{1}(\zeta_{t-s_{2}})\cos\left(\omega\tau
t+\alpha_{2,s_{1},s_{2}}\right)$ $\displaystyle+$
$\displaystyle\sum_{s_{1}=2}^{8}{\bm{q}}_{s_{1}}P_{2}(\zeta_{t-s_{1}})\cos\left(\omega\tau
t+\alpha_{2,s_{1}}\right),$

where $P_{n}(\zeta)$ represents the $n^{\rm th}$-order Legendre polynomial,
and the phases $\alpha_{1}$ and $\alpha_{2}$ depend on the initial values.
Note that although TIPCs depend on the initial values of the system in
general, the final outcome of TIPCs is the same in this case because the
choice of the initial values does not affect the coefficient vectors of Eq.
(35). The coefficient vectors for the first and second-order terms are
${\bm{p}}$ and ${\bm{q}}\in\mathbb{R}^{2}$, respectively, indicating that the
TIPCs of the system are composed of capacities of the product of the time-
varying basis vectors $\psi_{t}\in\\{\cos\omega\tau t,~{}\sin\omega\tau t\\}$
and PCs $\\{P_{n}(\zeta_{t-s})\\}$.

In addition, in the case of using the conventional IPC, the total IPC
saturates at the rank of the state matrix, $r$, only in the system with the
negative maximum Lyapunov exponent [19]. In the limit cycle system we adopted,
the Lyapunov exponents of Eqs. (19)–(21) are zero in the azimuthal and
negative in the radial direction; therefore, the maximum Lyapunov exponent is
zero. Furthermore, the system does not satisfy the ESP because it is time-
variant, and the phases $\alpha_{1}$ and $\alpha_{2}$ depend on the initial
values. Although the system, whose rank is two, does not satisfy these
conventional conditions, the total TIPC saturates at $C_{tot}=1.987$ ($1^{\rm
st}$-order, $C_{tot}=1.952$; $2^{\rm nd}$-order, $C_{tot}=3.45\times 10^{-2}$)
as shown in FIG. 1(b).

Therefore, these results suggest that the information processing that was
lacking with conventional IPC can be measured by adding time-varying bases,
and the total TIPC can reach the rank even if the maximum Lyapunov exponent is
not negative or the system state depends on the initial values.

###  III.5 Application #1: the benchmark task

To illustrate the usefulness of our theory, we demonstrate the information
processing capabilities of three systems. First, we analyze the computational
capabilities required to emulate a simple model for a time-series benchmark
test, which is a well-known NARMA10 model. The NARMA10 task is widely utilized
to evaluate the computational capabilities of dynamical systems, but the
meaning of predicting this model is unknown. We classified the model with
certain parameter regions to be time-invariant (see the Appendix). In FIG. 2,
$p$ expresses the probability of not diverging as a function of $\sigma$ for
different random series $\\{\zeta_{t}\\}$, and the model is stable with
certain parameters. Using the non-divergent model, we estimated its IPCs for
the target Legendre-chaoses with delayed time step $s(<16)$ and degree
$n_{s}~{}\left(\sum_{s}n_{s}<9\right)$.

Since the NARMA10 model is a one-dimensional system, the total capacity is
one, but the breakdown of the IPC changes with some parameters. FIG. 2(a) and
(b) show that the capacities using the uniform random variable in an
asymmetric range $u_{t}\in[0,\sigma]$ differ significantly from those using
the input in a symmetric range $u_{t}\in[-\sigma,\sigma]$. The capacity
breakdown with the symmetric input includes only the second-order capacities
($\sum_{s}n_{s}=2$). In contrast, the capacities with the asymmetric input
contain the first-order ones because the input term emerges as
$u_{t-9}u_{t}=\sigma^{2}/4(P_{1}(\zeta_{t})P_{1}(\zeta_{t-9})+P_{1}(\zeta_{t})+P_{1}(\zeta_{t-9})+1)$
in Eq. (17), including the first-order terms of $\zeta$. Hence, in the case of
using an asymmetric input, we can regard the model as the system receiving the
first-order inputs, which are retained for a few steps. These results suggest
that the input should be changed according to the dynamical system when one
uses the NARMA10 task; for example, as the nodes of an ESN are represented by
an odd function, such as a hyperbolic tangent, the ESN has only odd capacities
[19]. From this property and our results, the ESN with $u_{t}\in[0,\sigma]$
emulates the NARMA10 model, whereas the one with $u_{t}\in[-\sigma,\sigma]$
does not predict it at all.

Figure 2:  The IPC breakdown relative to $\sigma$ and the probability that
$y_{t}$ does not diverge, $p$. The labels indicate representative combinations
of $\\{\\{n_{s},s\\}\\}$, where $n_{s}$ is the degree of the polynomial, $s$
is the delayed time step of the input, and the desired output is
$\prod_{s}P_{n_{s}}(\zeta_{t-s})$. The labels for other combinations are
omitted. The probability $p$ represents the proportion of $y_{t}$ values that
do not diverge to infinity and is calculated using 100 input time-series
$\zeta_{t}$ generated from 100 random seeds. In all the figures, the
capacities are not stacked if the output diverged or $\sigma=0$. In (a) and
(b), the two ranges of $u_{t}$ are as follows: (a) $u_{t}$ follows the uniform
distribution in the range $[-\sigma,\sigma]$ and (b) in the range
$[0,\sigma]$.

###  III.6 Application #2: the machine learning network

Second, we analyzed the performance of a machine learning network that solved
the benchmark task using the ESN as an example. As shown in FIG. 3(a), we also
emulated target NARMA10 model using 50-node ESNs with three activation
functions—linear, hyperbolic tangent, and analog integral functions—and
compared the output of the ESN $\hat{y}_{t}$ and the target $y_{t}$ with the
normalized root-mean-square errors (NRMSEs). For all the functions, the NRMSE
decreased as the spectral radius $\rho$ increased and increased when $\rho\geq
1$. To analyze the outputs of the ESNs training the NARMA10 model, the IPCs of
the output $\hat{y}_{t}$ were estimated. FIG. 3(b)–(d) shows the change in the
IPC breakdown of the output $y_{t}$ with the increase in the spectral radius
$\rho$ of the ESN with the linear, hyperbolic tangent, and analog integral
functions, respectively. To emulate the NARMA10 model, whose IPC breakdown is
shown in FIG. 3(e), the nine types of Legendre-chaoses
$\\{P_{1}(\zeta_{t-s})\\}_{s=1,2,3,10,11,12}$ and
$\\{P_{2}(\zeta_{t-s})\\}_{s=1,2,3}$ are mainly required in a certain ratio.
According to FIG. 3(b)–(d), for any activation function, as $\rho(<1)$
increases, $\\{P_{1}(\zeta_{t-s})\\}_{s=1,2,3,10,11,12}$ approaches the
required rate. However, the three types of second-order IPCs
$\\{P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9})\\}_{s=1,2,3}$ are almost zero in
ESNs with an activation function. Therefore, in the NARMA10 task with ESNs,
performance is compared using only the first-order IPCs.

To investigate why the second-order IPCs did not appear in the breakdown, we
estimated the IPCs from ESN states. FIG. 3(f)–(h) shows the change in the IPC
breakdown with the increase in the spectral radius $\rho$ of the ESN with a
linear, hyperbolic tangent and an analog integral function, respectively. The
rank $r$ of the ESN state increases with $\rho(\leq 1)$ and corresponds to
$C_{tot}$. As these three breakdowns have low-target second-order IPCs (each
capacity is less than 0.2), the ESNs cannot emulate second-order Legendre-
chaoses.

The above results demonstrate that our method is capable of decomposing the
computational capability of the machine learning network before and after
training, as well as clarifying whether the computational components required
for the task have been extracted through learning and exist in the original
network.

Figure 3:  The performance breakdown of NARMA10-trained ESNs with the IPC.
Frame (a) shows the NRMSE of the NARMA10 task relative to the spectral radius
$\rho$ of ESNs composed of linear, hyperbolic tangent, or analog integrator
functions. Panels (b), (c), and (d) display the relationship between $\rho$
and the IPC breakdown of the output of the NARMA10-trained ESN with linear,
hyperbolic tangent, and analog integrator functions, respectively. The labels
indicate representative combinations of $\\{\\{n_{s},s\\}\\}$, where $n_{s}$
is the degree of the polynomial, $s$ is the delayed time step of the input,
and the desired output is $\prod_{s}P_{n_{s}}(\zeta_{t-s})$. In the case of
first-order IPCs, the representative combinations—{{1,1}}, {{1,2}}, {{1,3}},
{{1,10}}, {{1,11}}, and {{1,12}}—are shown, and other combinations are
summarized as “rest of 1st order”. The representative combinations of second-
order IPCs—{{1,1},{1,10}}, {{1,2},{1,11}}, and {{1,3},{1,12}}—are barely held
by the ESN, and thus, other combinations are summarized as “rest of 2nd
order”. Frame (e) shows the IPC breakdown of the target NARMA10 model with
input in the $[0,0.45]$ range. Panels (f), (g), and (h) represent the
relationship between $\rho$ and the IPC breakdown of the ESN state with
linear, hyperbolic tangent, and analog integrator functions, respectively. The
red line denotes the rank of state $r$, which corresponds to the total IPC in
the time-invariant domain.

###  III.7 Application #3: the real neural network

Finally, to show the broad applicability of our theory, we prepared a
dissociated culture of neurons for a physical reservoir, which is an open non-
equilibrium system that fluctuates due to external inputs and has parameters
that can be considered time-dependent. Real neurons extracted from the
cortices of rat embryos were pharmacologically and mechanically isolated and
then seeded on an electrode array. After the culture matured (FIG. 4[a]), we
constructed a physical reservoir using electrodes with active neurons (FIG.
4[b]; see the Appendix). We repeatedly applied bipolar waves with a 10, 20, or
30 ms interpulse interval (IPI) to 29 stimulus electrodes, whose amplitude
$\zeta_{t}$ follows a Gaussian distribution with mean $\mu=$ 200, 300, or 400
mV and standard deviation $\sigma=50$ mV (FIG. 4[c]). Furthermore, we computed
the number of spikes in an IPI-width bin from $N=792$ measurement electrodes
as the reservoir states (FIG. 4[d]). As a result, we obtained a long single
trajectory of the activation of the electrodes according to the input stream,
which consisted of 20,000 time steps and was used for our TIPC analysis.

Using these data, we computed the TIPCs of the physical system. FIG. 4(a)
shows the first-order TIPCs of the delay step $s$, which contain the memory
function (MF) [49, 50] and the four temporal memory functions (TMFs), which
had time-varying target $z_{t}=P_{1}(\zeta_{t-s})\cos(n\omega t)$ or
$z_{t}=P_{1}(\zeta_{t-s})\sin(n\omega t)$ $(n=1,2)$. Note that $n\omega=2\pi
n/T~{}(n=1,\ldots,T/2)$ denotes the frequency. The TMFs monotonically decay
with an increase in $s$, as well as the MF. Next, to investigate the frequency
characteristics of the TIPC, we plotted the TIPC with $s=1$. FIG. 4(b) and (c)
shows the relationship between the frequency and TIPCs with the time-varying
cosine and sine targets, respectively. We term such a graph the TIPC spectrum.
Both spectra have larger TIPCs at lower frequencies. Therefore, the input was
processed by the coupling terms of the low frequency sinusoidal wave and the
past input, suggesting that the superposition of these waves represents a
gradual trend—which may be caused, for example, by synaptic plasticity and
neural adaptation—and thus, information processing was also embedded in the
non-stationary changes.

FIG. 4(d)–(f) illustrates the total capacities $C_{tot}$ with different $\mu$
and IPI. Every $C_{tot}$ contained time-varying IPCs, and as the IPI
decreased, the ratio of time-varying IPC to $C_{tot}$ increased. As with the
time-invariant case, the time-variant IPC increased as the degree decreased.
Therefore, the total capacity of the dissociated culture of neurons contained
time-varying IPCs in all cases. This type of information processing could not
be elucidated by the conventional IPC, suggesting that our proposed measure is
effective.

Figure 4:  TIPCs of the dissociated culture of neurons. Frame (a) shows the
dissociated culture of neurons on the electrode array. Panel (b) shows the
spontaneous firing rate of the culture on the electrode array, where $(x,y)$
denotes the position of electrode. The rate is filtered by the Gaussian
kernel. The red dots represent the 29 stimulation electrodes. Frames (c) and
(d) represent the amplitude of stimulus $\zeta_{t}$, which follows a Gaussian
distribution ($\mu=200$ mV), and the spike count of each electrode
${\bm{x}}_{t}$ ($\mu=200$ mV and IPI = 10 ms), respectively. Frame (e)
illustrates the temporal memory functions (TMFs) of delay step $s$. (e)–(g)
show the TMFs when IPI = 30 ms and $\mu$ = 400 mV. TMFs with a time-invariant
target $z_{t}=P_{1}(\zeta_{t-1})$ and four time-variant targets—i.e.,
$z_{t}=P_{1}(\zeta_{t-1})\cos(n\omega t)$ and $P_{1}(\zeta_{t-1})\sin(n\omega
t)~{}(n=1,2)$—are plotted. Panels (f) and (g) show the first-order TIPC
spectrum with the cosine $\zeta_{t}=P_{1}(\zeta_{t-1})\cos(n\omega t)$ and
sine target $\zeta_{t}=P_{1}(\zeta_{t-1})\sin(n\omega t)$, respectively.
Frames (h), (i), and (j) depict the total capacity of the culture with IPI =
10, 20, and 30 ms, respectively, and the horizontal axis is the mean of
amplitude $\mu$. The hatched bar represents the total capacity with time-
varying $n^{\rm th}$ order polynomials.

##  IV Discussion

###  IV.1 The relationship between an attractor and information processing

In the present paper, random variables were given as input to the time-
invariant systems, whose states were represented by echo functions, which
depend not on the time $t$ but only on the past input time-series. In the case
where the system with no input converges to a fixed-point (e.g., the
intersection point of the NARMA10 system), this function illustrates that the
system stays at a fixed-point attractor and fluctuates around the fixed-point
due to the input. The conventional IPC targets time-invariant systems and
quantifies the input processing of a state that depends only on the input
time-series. Therefore, the conventional IPC sometimes represents the input
information processing performed around a certain fixed-point attractor.

In addition, the IPC is the squared norm of the coefficient vector of the
temporal basis vectors expanded with PCs. Since these coefficient vectors
obviously change depending on the fixed-point, different types of information
processing are performed at different fixed points. However, the IPC extended
for time-variant systems was applied to the limit cycle system, which does not
satisfy the prerequisites for RC. From its TIPC estimates, the expanded
solution contains the coupling terms of the input and the time-varying terms,
$P_{1}(\zeta_{t-s_{1}})\cos(\omega\tau t+\alpha_{1,s_{1}})$,
$P_{2}(\zeta_{t-s_{1}})\cos(\omega\tau t+\alpha_{2,s_{1}})$, and
$P_{1}(\zeta_{t-s_{1}})P_{1}(\zeta_{t-s_{2}})\cos(\omega\tau
t+\alpha_{2,s_{1},s_{2}})$, showing that the processed input represents the
amplitude scale of the sinusoidal. Therefore, the fluctuation of a periodic
attractor determined by input can be processed around the limit cycle.

Based on these findings, we conclude that the conventional IPC can evaluate
computational capabilities around a fixed point, while the TIPC can also
evaluate capabilities around a periodic attractor. Since the recent RC
framework exploits not only fixed-points or periodic attractors but also
chaotic ones [51, 52, 53, 54], future work should examine the relationship
between various attractors and information processing.

###  IV.2 Methods to utilize a time-variant system as a computational resource

In demonstrating the TIPCs for a limit cycle system, we showed that the state
of the system can include coupling terms of the past input time-series. Since
the coupling term is represented by the product of the time-varying basis
$\psi_{t}$ and PC $\phi_{t}$, the conventional IPC for time-invariant terms
could not quantify the amount of information. The TIPC indicates that even in
a time-variant system, information can be processed by the PC $\phi_{t}$ in
the coupling term. In addition, to date, RC has trained a static readout
weight by linear regression, assuming that the ESP or FMP is satisfied; the
state in the reservoir is a function of the input time-series. Since the
target output is described by the input history, performance decreases when a
time-dependent reservoir (e.g., an ESN with a spectral radius of $\rho\geq 1$)
is used. As the number of input history terms in the expanded state decreases,
the performance drops, while the number of coupling terms increases. The input
information in a coupling term can be used for the task and can be extracted
by giving a readout weight that cancels out $\phi_{t}$ in the coupling term
(e.g., a time-varying weight). Therefore, even in a reservoir where the ESP or
FMP does not hold, the input time-series may be processed by TDPCs, and the
task can be successfully solved by using new types of readout.

In this connection, a method already exists for exploiting periodic systems as
computational resources. As shown in FIG. 5, the time-multiplexing technique
[5] switches the input $u_{t}$ with the time width $\tau$ (FIG. 5[a]) and
extracts ${\bm{x}}(t)=[x(t+\tau),x(t+2\tau),\ldots,x(t+(N-1)\tau)]^{\top}$ to
virtually increase the number of nodes in the reservoir and improve the
computational capabilities. If $x(t)$ is a periodic function and oscillates
with a period specific to $u(t)$ (FIG. 5[b]), applying an input with the same
width $\tau$ as the period can extract time-invariant virtual nodes because
the $i^{\rm th}$ virtual node always corresponds to $x(t)$ at a certain phase
and is not affected by the periodic fluctuation. However, if the period does
not match $\tau$ (FIG. 5[c]), the phase shifts, and the scheme cannot exploit
the computational capabilities of the periodic system. Therefore, time-
multiplexing can be interpreted as a method capable of extracting a
computational resource by transforming a time-varying system into a time-
invariant system. Thus, a time-variant system can process input information
through coupling terms, and rich input information can be extracted from time-
variant systems by designing readout for the systems. In the future,
information processing using a time-varying reservoir and novel design methods
for readouts will be investigated.

Figure 5:  The time-multipleixing method transforms a periodic system into a
time-invariant one. (a) The input time-series $\\{u_{t}\\}$ is switched with
time width $\tau$. Frames (b) and (c) show a periodic state time-series whose
frequency is and is not consistent with that of the input, respectively. In
(b), the circles represent virtual nodes
${\bm{x}}_{1}(t)=\left[x_{1}\left(t\right),x_{1}\left(t+\frac{\tau}{N}\right),\ldots,x_{1}\left(t+\frac{(N-1)\tau}{N}\right)\right]^{\top}$.

###  IV.3 Extension to TIPC with Multiple Input Variables

In the present paper, TIPC was limited to one type of input, but it could be
easily extended to a multiple input version. Let $M$ independent stochastic
inputs that follow multiple arbitrary distributions be
$\zeta_{t}^{(1)},\ldots,\zeta_{t}^{(M)}\in\mathbb{R}$, and the state equation
in Eq. (4) is modified as
${\bm{x}}_{t+1}={\bm{f}}({\bm{x}}_{t},\zeta_{t}^{(1)},\ldots,\zeta_{t}^{(M)})$,
whose state for a time-invariant system can be expanded by aPCs with multiple
input variables [55]. As all of these aPCs are orthogonal, and their
orthogonality is defined by the same inner product as Eqs. (32) and (33),
which are also common to TDPCs, it is clear that we can expand the state for a
time-variant system using the TDPCs and define its TIPCs. Physical systems
often receive various inputs from the external environment, resulting in such
multiple input-driven systems. For example, in the dissociated culture of
neurons, the ranks of the state are 724–792, but the total capacities are less
than 3.5 (FIG. 4[c]), which are much smaller than their ranks. One possible
speculation regarding this issue is that in our scheme, the state may be
expressed as a function of multiple stochastic inputs, including the one we
gave (e.g., electrical noise, synaptic noise, thermal noise, and shot noise
[56]). Of all the inputs, the ones we could observe is limited, and the total
capacity computed only from the inputs did not reach the rank. Therefore,
physical systems can receive unobservable inputs, which disturb the
examination of all capacities.

##  V Conclusion

This paper attempted to clarify the unknown relationship between the PC
expansion and IPC by deriving the IPC from the PC-expanded system. To
illustrate this relation, we showed that the IPC can be measured using gPCs
and aPCs. In addition, using the NARMA10 model, we concretely described the
relationship and showed the usefulness of our theory. Next, taking into
account the characteristics of the general solution of the input-driven
dynamical system, we proposed the IPC for time-variant systems—called the
TIPC. To demonstrate that the time-variant system has such TIPCs, we
investigated the TIPC breakdown of the limit cycle system. The primary results
are summarized as follows:

  * •

Using SVD, we can obtain the orthogonal temporal basis vectors from the state
time-series. These vectors can be expanded with PCs to obtain the coefficient
vector of each basis. The IPC is equivalent to the squared norm of the
coefficient of the PC used as the target output. Therefore, the expanded basis
coefficients represent the amount of input processing information.

  * •

Using eight types of polynomials in the Askey scheme and Gram-Schmidt PCs, we
estimated the IPCs of a one-dimensional ESN, whose total IPCs were equal to
one. These results indicate that various types of PCs within the Askey scheme
and Gram-Schmidt orthogonalization are suitable for the target output of the
IPC.

  * •

We calculated the TIPC of a time-variant system using the simple limit cycle
and revealed that the input information processing is performed by time- and
input time-series-dependent terms.

  * •

IPC analysis revealed that the NARMA10 model is mainly composed of
$P_{1}(\zeta_{t-s})~{}(s=1,2,3,10,11,12)$ and
$P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9})~{}(s=1,2,3)$. The NARMA10 benchmark
task can be solved by holding the nine types of input information in a
reservoir. Consequently, combining the IPC and PC expansion provides a clear
and concise picture of information processing.

  * •

The dissociated culture of neurons had not only the time-invariant IPC but
also the time-variant IPC, suggesting that the TIPC reveals the information
processing in a trend—for example, synaptic plasticity and neural
adaptation—that has been left behind to date.

The above results suggest that the connection between the IPC and PC expansion
allows for a simpler description of information processing in dynamical
systems. In the future, information processing using time-variant systems
should also be elucidated. This scheme can be applied to non-stationary
systems and thus may be suitable for elucidating information processing in
neural circuits which has been overlooked so far. In addition, it can be
applied not only to the neural systems but also to other physical systems that
can be time-variant—for example, fluid, quantum, spintronics, and optical
systems. For example, recently it is reported that some types of vortex
generated when a fluid flows past a bluff body can be used as an information
processing device [57]. In their analysis, they found that near the critical
Reynolds number, where the flow exhibits a twin vortex before the onset of the
Karman vortex shedding associated with the Hopf bifurcation, the information
processing capability was maximized. This was also characterized by the
breakdown of ESP. According to our results, it may be possible to evaluate the
type and amount of information processing even in the Karman vortex shedding,
which would be a direction for future work.

###### Acknowledgements.

This paper is based on results obtained from the Exploration of Neuromorphic
Dynamics towards Future Symbiotic Society project commissioned by NEDO,
KAKENHI grant (17K20090), AMED (JP18dm0307009) and Asahi Glass Foundation. H.
T. was partly supported by JSPS KAKENHI (20H04252), AMED (JP20dm0307009), the
Naito Foundation, and the Asahi Glass Foundation. K. N. was supported by the
JSPS KAKENHI Grant Number JP18H05472 and MEXT Quantum Leap Flagship Program
(MEXT Q-LEAP) Grant Number JPMXS0118067394.

##  Appendix A An ESN and LSM

In an ESN composed of $N$-nodes, the $i^{\rm th}~{}(i=1,2,\ldots,N)$ node
state $x_{i,t}$ at the $t^{\rm th}$ time step can be written as follows:

$\displaystyle x_{i,t+1}$ $\displaystyle=$ $\displaystyle
f\left(\sum_{j=1}^{N}w_{ij}x_{j,t}+w_{in,i}u_{t}\right),$ (A1)

where $f$ is the activation function, and $w_{ij}$ and $w_{in,i}$ are the
internal and input weights, respectively.

To emulate the target output $z_{t}$, we use linear regression to obtain an
estimate of $z_{t}$, $\tilde{z}_{t}$, as follows:

$\displaystyle\tilde{z}_{t}$ $\displaystyle=$
$\displaystyle\tilde{\bm{w}}_{out}^{\top}\cdot{\bm{x}}_{t},$ (A2)
$\displaystyle\tilde{\bm{w}}_{out}$ $\displaystyle=$
$\displaystyle\arg\min_{{\bm{w}}_{out}}\sum_{t=1}^{T}\left(z_{t}-{\bm{w}}_{out}^{\top}\cdot{\bm{x}}_{t}\right)^{2},$
(A3)

where ${\bm{w}}_{out}$ and $\tilde{\bm{w}}_{out}\in\mathbb{R}^{L}$ are the
weight and solution vector for the target output, respectively. This learning
method does not affect the state of the reservoir but places a constraint on
${\bm{x}}_{t}$. Jaeger [2] and Maass [3] independently developed RC by
integrating ESNs and LSMs, respectively. The prerequisites for ${\bm{x}}_{t}$
differ in ESNs and LSMs.

An ESN requires ${\bm{x}}_{t}$ to be an echo function, which is a function of
only the past input time-series ${\bm{u}}_{t}=\\{u_{t-s}\\}_{s=1}^{t}$. This
dynamical property is referred to as the ESP [13, 14, 15]. We examine this
feature by giving the same input time-series to two systems with different
initial values and checking whether the two states coincide after a
sufficiently long period of time. If the input is noise, this phenomenon is
called noise-induced synchronization [58, 59]. Furthermore, when the input is
generated from a deterministic system, the phenomenon wherein the state is
synchronized with the input is called generalized synchronization [60].
Therefore, ESP is related to the synchronization phenomenon of nonlinear
dynamical systems.

In addition, LSMs impose a prerequisite on the power series expansion of
states. If the system is time-invariant—i.e., the state does not depend on
time—and retains exponentially decaying inputs, its state can be approximated
by the Volterra series [61], which is a series expansion with non-orthogonal
bases. Accordingly, the Volterra series operator [62] can be expanded in a
non-orthogonal power series, as follows:

$\displaystyle{\bm{x}}_{t}=\sum_{n=1}^{\infty}\prod_{s_{1}=1}^{t}\prod_{s_{2}=s_{1}}^{t}\cdots\prod_{s_{n}=s_{n-1}}^{t}{\bm{g}}_{s_{1}\cdots
s_{n}}u_{t-s_{1}}\cdots u_{t-s_{n}},$

where ${\bm{g}}_{s_{1}\cdots s_{n}}\in\mathbb{R}^{L}$ is the $n^{\rm th}$
Volterra kernel. Such a memory decay feature is called the FMP [3, 16, 17] and
is considered a prerequisite for LSMs.

##  Appendix B Polynomial chaos expansion

A deterministic dynamical system with a stochastic input $\zeta_{t}$ can be
described as an operator of $\\{\zeta_{t-s}\\}_{s=1}^{\infty}$ according to
polynomial chaos expansion [23], which is a series expansion using the target
outputs of IPC as the bases—i.e., multivariate orthogonal polynomials of the
random variables, $z_{t}^{(i)}$, described in Eq. (5). These multivariate
polynomials are referred to as polynomial chaoses (PCs), and the space spanned
by the PCs is called homogeneous chaos [23], expressed as

$\displaystyle{\bm{x}}_{t}=\sum_{i=1}^{\infty}{\bm{c}}_{i}z_{t}^{(i)},$ (B1)

where ${\bm{c}}_{i}\in\mathbb{R}^{L}$ is the $i^{\rm th}$ coefficient vector.

Figure B1:  The relationship between polynomials in the Askey scheme. Each polynomial is represented with the hypergeometric series ${}_{r}F_{s}(p)$, where $p$ represents the number of parameters substituted in $a_{i}~{}(i=1,2,\ldots,r)$, $b_{i}~{}(i=1,2,\ldots,s)$, or $z$. An upper polynomial with the limit of a certain parameter or parameters becomes a lower polynomial connected with a line.  Table B1:  Sets of random variables and polynomial chaoses  Support of random variable | Random variable $\zeta_{t}$ | Polynomial chaos $z_{t}^{(i)}$
---|---|---
Continuous | $(-\infty,\infty)$ | Gaussian | Hermite-chaos
| $[0,\infty)$ | Gamma | Laguerre-chaos
| $[a,b]$ | Beta | Jacobi-chaos
| $[a,b]$ | Uniform | Legendre-chaos
Discrete | $\\{0,1,\ldots\\}$ | Poisson | Charlier-chaos
| $\\{0,1,\ldots,N\\}$ | Binomial | Krawtchouk-chaos
| $\\{0,1,\ldots\\}$ | Negative binomial | Meixner-chaos
| $\\{0,1,\ldots,N\\}$ | Hypergeometric | Hahn-chaos

###  B.1 Generalized polynomial chaos

Generalized polynomial chaos (gPC) [24] is PC composed of the univariate
polynomial $\mathcal{F}_{n}(\zeta)$ included in the Askey scheme [63] tree
(FIG. B1). The Askey scheme represents various orthogonal polynomials (e.g.,
Hermite, Jacobi, Laguerre, and Charlier) using the hypergeometric series
${}_{r}F_{s}$ of $x$, along with parameters $a_{1},\ldots,a_{r}$ and
$b_{1},\ldots,b_{s}$:

${}_{r}F_{s}\left(x;\begin{matrix}a_{1},\ldots,a_{r}\\\ b_{1},\ldots,b_{s}\\\
\end{matrix}\right)=\sum_{k=0}^{\infty}\frac{(a_{1})_{k}\cdots(a_{r})_{k}}{(b_{1})_{k}\cdots(b_{s})_{k}}\frac{x^{k}}{k!},$
(B2) $\displaystyle(a)_{n}=\begin{cases}1&(n=0)\\\
a(a+1)\cdots(a+n-1)&(n=1,2,\ldots)\end{cases},$

where $(a)_{n}$ is the Pochhammer symbol. For example, the Laguerre polynomial
with a parameter $\alpha$, $L^{(\alpha)}_{n}(\zeta)$ can be written as

$\displaystyle
L^{(\alpha)}_{n}(\zeta)=\frac{(\alpha+1)_{n}}{n!}{{}_{1}F_{1}}\left(\zeta;\begin{matrix}-n\\\
\alpha+1\end{matrix}\right).$

Note that in FIG. B1, an upper polynomial with the limit of a certain
parameter or parameters becomes a lower polynomial connected with a line; for
example, the Laguerre polynomial becomes the Hermite polynomial $H_{n}(\zeta)$
by taking the following limit of $\alpha$:

$\displaystyle\lim_{\alpha\rightarrow\infty}\left(\frac{2}{\alpha}\right)^{n/2}L_{n}^{(\alpha)}\left((2\alpha)^{1/2}\zeta+\alpha\right)=\frac{(-1)^{n}}{n!}H_{n}(\zeta).$

The $i^{\rm th}$ target $z^{(i)}_{t}$ is represented by the product of the
$n^{(i)}_{s}$-th order polynomial of $\zeta_{t}$ delayed by $k^{(i)}_{s}$
steps, $\mathcal{F}_{n^{(i)}_{s}}(\zeta_{t-k^{(i)}_{s}})$. When the sets of
degree $n$ and delay step $k$ are given by the $i^{\rm th}$ family of sets
$\mathcal{N}_{i}=\\{(n_{1}^{(i)},k_{1}^{(i)}),(n_{2}^{(i)},k_{2}^{(i)}),\ldots\\}$,
the target output is represented as

$\displaystyle z_{t}^{(i)}$ $\displaystyle=$
$\displaystyle\prod_{s}\mathcal{F}_{n_{s}^{(i)}}\left(\zeta_{t-k_{s}^{(i)}}\right).$
(B4)

Using the hypergeometrical series ${}_{r}F_{s}(\zeta)$,
$\mathcal{F}_{n}(\zeta_{t-k})$ can be applied to the following eight types of
orthogonal polynomials.

#### The Hermite polynomial and a Gaussian distribution

The $n^{\rm th}~{}(n=1,2,\ldots)$ order Hermite polynomial $H_{n}(\zeta)$ is
given by

$\displaystyle H_{n+1}=\zeta H_{n}-nH_{n-1},$ (B5)

where $H_{0}=1$, and $H_{-1}=0$. $\zeta$ follows a standard normal
distribution

$\displaystyle
f(\zeta)=\frac{1}{\sqrt{2\pi}}\exp\left(\frac{\zeta^{2}}{2}\right).$ (B6)

#### The Laguerre polynomial and a gamma distribution

The $n^{\rm th}~{}(n=1,2,\ldots)$ order Laguerre polynomial
$L^{(\alpha)}_{n}(\zeta)$ is given by

$\displaystyle L^{(\alpha)}_{n}(\zeta)$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{n}(-1)^{i}\left(\begin{matrix}n+\alpha\\\
n-i\end{matrix}\right)\frac{\zeta^{i}}{i!},$ (B7)

where $\zeta$ follows a gamma distribution

$\displaystyle
f(\zeta)=\frac{1}{\Gamma(\alpha+1)\beta^{\alpha+1}}\zeta^{\alpha}e^{-\zeta/\beta}.$
(B8)

The parameter $\alpha>-1$, $\beta=1$, and $\Gamma(\cdot)$ is the gamma
function. In FIG. 1(b), the parameter $\alpha$ was set to $1$.

#### The Jacobi polynomial and a beta distribution

The $n^{\rm th}~{}(n=1,2,\ldots)$ order Jacobi polynomial
$P^{(\alpha,\beta)}_{n}(\zeta):=P_{n}$ is given by

$\displaystyle P_{n+1}$ $\displaystyle=$
$\displaystyle\frac{(\gamma_{n}+1)(\gamma_{n}(\gamma_{n}+2)\zeta+\alpha^{2}-\beta^{2})P_{n}}{2(n+1)(\gamma_{n}-n+1)\gamma_{n}}$
(B9) $\displaystyle-$
$\displaystyle\frac{2(n+\alpha)(n+\beta)(\gamma_{n}+2)}{2(n+1)(\gamma_{n}-n+1)\gamma_{n}}P_{n-1},$

where $\gamma=2n+\alpha+\beta$, and $\zeta$ follows a beta distribution in the
range of $[-1,1]$:

$\displaystyle f(\zeta)$ $\displaystyle=$
$\displaystyle\frac{1}{B(\alpha,\beta)}\zeta^{\alpha-1}(1-\zeta)^{\beta-1},$
(B10) $\displaystyle B(\alpha,\beta)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}t^{\alpha-1}(1-t)^{\beta-1}dt.$ (B11)

In FIG. 1(c), the parameters $(\alpha,\beta)$ were set to $(-0.25,-0.25)$.

#### The Legendre polynomial and a uniform distribution

The $n^{\rm th}~{}(n=1,2,\ldots)$ order Legendre polynomial $P_{n}(\zeta)$ is
given by

$\displaystyle P_{n}(\zeta)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\lfloor
n/2\rfloor}\frac{(-1)^{-k}}{2^{n}}\left(\begin{matrix}n\\\
k\end{matrix}\right)\left(\begin{matrix}2n-2k\\\
n\end{matrix}\right)\zeta^{n-2k},$

where $\lfloor\cdot\rfloor$ represents the floor function, and $\zeta$ follows
a uniform distribution in the range of $[-1,1]$:

$\displaystyle f(\zeta)$ $\displaystyle=$ $\displaystyle\frac{1}{2}.$ (B13)

#### The Charlier polynomial and a Poisson distribution

The $n^{\rm th}~{}(n=1,2,\ldots)$ order Charlier polynomial $C_{n}(\zeta;a)$
is given by

$\displaystyle C_{n}(\zeta;a)=\sum_{i=0}^{n}\left(\begin{matrix}n\\\
i\end{matrix}\right)\left(\begin{matrix}\zeta\\\
i\end{matrix}\right)i!(-a)^{n-i},$ (B14)

where $\zeta$ follows a Poisson distribution

$\displaystyle f(\zeta)$ $\displaystyle=$
$\displaystyle\frac{e^{-a}a^{\zeta}}{\zeta!}.$ (B15)

In FIG. 1(e), the parameter $\alpha$ was set to $6$.

#### The Krawtchouk polynomial and a binomial distribution

The $n^{\rm th}~{}(n=1,2,\ldots)$ order Krawtchouk polynomial
$K_{n}(\zeta;p,N)$ is given by

$\displaystyle
K_{n}(\zeta;p,N)=\sum_{i=0}^{n}(-1)^{n-i}\left(\begin{matrix}N-\zeta\\\
n-i\end{matrix}\right)\left(\begin{matrix}\zeta\\\
i\end{matrix}\right)p^{n-i}(1-p)^{i},$ (B16)

where $\zeta$ follows a binomial distribution

$\displaystyle f(\zeta)$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}N\\\
\zeta\end{matrix}\right)p^{\zeta}(1-p)^{N-\zeta}.$ (B17)

In FIG. 1(f), the parameters $(p,N)$ were set to $(0.5,10)$.

#### The Meixner polynomial and a negative binomial distribution

The $n^{\rm th}~{}(n=1,2,\ldots)$ order Meixner polynomial
$M_{n}(\zeta;\beta,c):=M_{n}$ is given by

$\displaystyle
M_{n+1}=\frac{\\{(c-1)\zeta+n+c(n+\beta)\\}M_{n}-nM_{n-1}}{c(n+\beta)},$ (B18)

where $M_{0}=1$ and $M_{-1}=0$. $\zeta$ follows a negative binomial
distribution

$\displaystyle f(\zeta;\beta,c)$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}\zeta-1\\\
\beta-1\end{matrix}\right)(1-c)^{\beta}c^{\zeta-\beta}.$ (B19)

In FIG. 1(g), the parameters $(c,\beta)$ were set to $(0.2,10)$.

#### The Hahn polynomial and a hypergeometric distribution

The $n^{\rm th}$ order Hahn polynomial $Q_{n}(\zeta;\alpha,\beta,N):=Q_{n}$ is
given by

$\displaystyle Q_{n+1}$ $\displaystyle=$
$\displaystyle\frac{A_{n}+C_{n}-\zeta}{A_{n}}Q_{n}-\frac{C_{n}}{A_{n}}Q_{n-1},$
(B20) $\displaystyle A_{n}$ $\displaystyle=$
$\displaystyle\frac{(n+\alpha+\beta+1)(n+\alpha+1)(N-n)}{(2n+\alpha+\beta+1)(2n+\alpha+\beta+2)},$
(B21) $\displaystyle C_{n}$ $\displaystyle=$
$\displaystyle\frac{n(n+\alpha+\beta+N+1)(n+\beta)}{(2n+\alpha+\beta)(2n+\alpha+\beta+1)},$
(B22)

where $Q_{0}=1$ and $Q_{-1}=0$. $\zeta$ follows a hypergeometric distribution
with $m=-\alpha-1$ and $n=-\beta-1$:

$\displaystyle f(\zeta;m,n,N)$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}m\\\
\zeta\end{matrix}\right)\left(\begin{matrix}n\\\
N-\zeta\end{matrix}\right)p^{\zeta}(1-p)^{N-\zeta}.$

In FIG. 1(h), the parameters $(m,n,N)$ were set to $(100,50,20)$.

###  B.2 Arbitrary polynomial chaos

In addition, arbitrary polynomial chaos (aPC) [25] is the PC for random
variables following an arbitrary probability distribution, and we can compute
univariate polynomials by applying the Gram-Schmidt orthogonalization
procedure. The $n^{\rm th}$ Gram-Schmidt polynomial
$\mathcal{F}_{n}(\zeta_{t})=\psi_{t}^{(n)}$ is obtained from the following
equations:

$\displaystyle\psi_{t}^{(n)}$ $\displaystyle=$
$\displaystyle\zeta_{t}^{n}-\sum_{i=0}^{n-1}c_{i}^{(n)}\psi_{t}^{(i)},$ (B24)
$\displaystyle c_{i}^{(n)}$ $\displaystyle=$
$\displaystyle\frac{({\bm{\psi}}^{(i)})^{\top}\cdot{\bm{\zeta}}^{n}}{({\bm{\psi}}^{(i)})^{\top}\cdot{\bm{\psi}}^{(i)}},$
(B25)

where ${\bm{\zeta}}^{n}=[\zeta_{1}^{n}\cdots\zeta_{T}^{n}]^{\top}$ and
${\bm{\psi}}^{(i)}=[\psi_{1}^{(i)}\cdots\psi_{T}^{(i)}]^{\top}$. Note that
$\psi_{t}^{(0)}=1$. From the univariate polynomial $\mathcal{F}_{n}(\zeta)$
and Eq. (5), we can calculate the polynomial chaos for an arbitrary input
distribution, Gram-Schmidt polynomial chaos.

If the input $\zeta_{t}$ follows a certain distribution, PC, including gPC and
aPC, is determined based on its orthogonality. If a weighting function
$w({\bm{\zeta}})$ specific to PCs $\\{z_{t}^{(i)}\\}_{i=1}^{\infty}$ exists,
the following orthogonality relations should be satisfied:

$\displaystyle\left<z_{t}^{(i)}z_{t}^{(j)}\right>$ $\displaystyle=$
$\displaystyle\left<(z_{t}^{(i)})^{2}\right>\delta_{ij},$ (B26)
$\displaystyle\left<f({\bm{\zeta}})g({\bm{\zeta}})\right>$ $\displaystyle=$
$\displaystyle\sum_{\bm{\zeta}}w({\bm{\zeta}})f({\bm{\zeta}})g({\bm{\zeta}}),$
(B27)

where $\delta_{ij}$ is the Kronecker delta function. Several sets of random
variables and gPCs were proposed [24], as shown in TABLE B1. Note that
expansion, in terms of the Hermite-chaos of Gaussian variables and Charlier-
chaos of Poisson variables, converges in the sense of $L_{2}$ according to the
Cameron-Martin theorem [64] and Ogura [65], respectively.

##  Appendix C Threshold of IPC

To remove the estimation error, $C({\bm{X}},{\bm{z}})$ was set to zero if it
was smaller than the threshold $\epsilon$:

$\displaystyle
C_{\epsilon}({\bm{X}},{\bm{z}})=\theta_{\epsilon}\left(C({\bm{X}},{\bm{z}})\right)C({\bm{X}},{\bm{z}}),$
(C1)

where $\theta_{\epsilon}(\cdot)$ is the Heaviside step function. The threshold
is determined using random shuffle surrogates. We prepared $N(=200)$
surrogates that are $z_{t}$ shuffled in the time direction and calculated IPCs
using the surrogates to obtain $N$ capacities. Furthermore, we let the
significance level be $\alpha(=1)\%$ and chose the original IPC, which exceeds
1.2–3 times the value in the top $\alpha/2\%$ of $N$ capacities. The above
operation was performed for each $z_{t}$, and significant IPCs were obtained.

Figure C1:  Properties of the NARMA10 dynamical system with
$(\alpha,\beta,\gamma,\delta)=(0.3,0.05,1.5,0.1)$. Panel (a) shows the time-
series of Eq. (17) with $u_{t}\in[0,0.01]$ (black) and $u_{t}\in[-0.01,0.01]$
(red) and (b) the stability diagram of $w^{(1)}_{t}=z_{t}^{(1)}$ and
$w_{t}^{(2)}=\sum_{i=1}^{10}z_{t}^{(i)}$ as $z_{t}^{(10)}>\delta$. The two
nullclines, $w_{t}^{(2)}=\frac{1-\alpha}{\beta}-\frac{\delta}{\beta
w_{t}^{(1)}}$ (solid black) and
$w_{t}^{(2)}=-\frac{\alpha}{\beta}+\frac{z_{t}^{(10)}-\delta}{\beta
w_{t}^{(1)}}$ (solid red), are displayed. The signs of $\Delta w_{t}^{(1)}$
and $\Delta w_{t}^{(2)}$ are shown in TABLE D1. Note that the point
$(z_{t}^{(10)},(1-\alpha)/\beta-\delta/\beta z_{t}^{(10)})$ is a saddle point.
Panels (c) and (d) show the basins of attraction of Eq. (17) relative to
$w_{0}^{(1)}$ and $w_{0}^{(2)}$ with $u_{t}=0$ and $u_{t}\in[0,0.4]$,
respectively, as well as the two nullclines from (b). In (c), (d), and (f),
each dot shows the time step at which $y_{t}$ diverges to infinity, and the
white dot represents the initial values for which $y_{t}$ does not do so.
Frame (e) illustrates the bifurcation diagram of Eq. (17), while (f) shows the
basin of attraction of Eq. (17) relative to $\sigma$ and $\psi$, where the
initial values were $y_{i}=\psi~{}(i=0,1,\ldots,9)$. The three largest
Lyapunov spectra $\lambda_{i}~{}(i=1,2,3)$ relative to $\sigma$ are shown in
(g).

##  Appendix D Classification of the NARMA10 model

###  D.1 Attractor analysis

To study the divergence and time-dependence of the NARMA10 model without
input, we analyzed its attractor. We defined new variables for a time-delay
system as $z_{t}^{(s)}\equiv y_{t+1-s}~{}(s=1,2,\ldots,10)$, and the model
without input ($\mu=\kappa=0$) was rewritten as follows:

$\displaystyle z_{t+1}^{(1)}$ $\displaystyle=$ $\displaystyle\alpha
z_{t}^{(1)}+\beta z_{t}^{(1)}\sum_{i=1}^{10}z_{t}^{(i)}+\delta,$ (D1)
$\displaystyle z_{t+1}^{(s)}$ $\displaystyle=$ $\displaystyle
z_{t}^{(s-1)}~{}(s=2,3,\ldots,10).$ (D2)

Eq. (D1) makes use of $z_{t}^{(1)}$ and $\sum_{i=1}^{10}z_{t}^{(i)}$; thus, we
defined $w_{t}^{(1)}=z_{t}^{(1)}$ and
$w_{t}^{(2)}=\sum_{i=1}^{10}z_{t}^{(i)}$, and the discrete derivatives $\Delta
w_{t}^{(1)}(=w_{t+1}^{(1)}-w_{t}^{(1)})$ and $\Delta
w_{t}^{(2)}(=w_{t+1}^{(2)}-w_{t}^{(2)})$ were derived from Eqs. (D1) and (D2)
as follows:

$\displaystyle\Delta w_{t}^{(1)}$ $\displaystyle=$
$\displaystyle(\alpha-1)w_{t}^{(1)}+\beta w_{t}^{(1)}w_{t}^{(2)}+\delta,$ (D3)
$\displaystyle\Delta w_{t}^{(2)}$ $\displaystyle=$ $\displaystyle\alpha
w_{t}^{(1)}+\beta w_{t}^{(1)}w_{t}^{(2)}+\delta-z_{t}^{(10)}.$ (D4)

From Eqs. (D3) and (D4), the nullclines can be obtained:

$\displaystyle w_{t}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{1-\alpha}{\beta}-\frac{\delta}{\beta w_{t}^{(1)}},$ (D5)
$\displaystyle w_{t}^{(2)}$ $\displaystyle=$
$\displaystyle-\frac{\alpha}{\beta}+\frac{z^{(10)}_{t}-\delta}{\beta
w_{t}^{(1)}}.$ (D6)

Note that the intersection point of the nullclines is a saddle point
$(w_{t}^{(1)},w_{t}^{(2)})=\bigl{(}z_{t}^{(10)},(1-\alpha)/\beta-\delta/\beta
z_{t}^{(10)}\bigr{)}$. Furthermore, we plotted the nullclines on the
$w_{t}^{(1)}$–$w_{t}^{(2)}$ plane to examine the increase and decrease in
$w_{t}^{(1)}$ and $w_{t}^{(2)}$ (FIG. C1[b] and TABLE D1). Since
$z^{(10)}_{t}$ depends on the time step $t$, the nullcline changes over time
and three regimes can be distinguished: $z^{(10)}_{t}<\delta$,
$z^{(10)}_{t}=\delta$, and $z^{(10)}_{t}>\delta$ (FIGs. C1[b] and D1). The
signs of $\Delta w^{(1)}$ and $\Delta w^{(2)}$ are shown in TABLEs D1–D3. From
the above results, we found that the model with no input has a fixed point
attractor at the saddle point.

Table D1:  The derivative table of $\Delta w_{t}^{(1)}$ and $\Delta w_{t}^{(2)}$ in the stability diagram ($z_{t}^{(10)}>\delta$; FIG. C1[b])  # area | $\Delta w_{t}^{(1)}$ | $\Delta w_{t}^{(2)}$ | # area | $\Delta w_{t}^{(1)}$ | $\Delta w_{t}^{(2)}$ | # area | $\Delta w_{t}^{(1)}$ | $\Delta w_{t}^{(2)}$ | # area | $\Delta w_{t}^{(1)}$ | $\Delta w_{t}^{(2)}$
---|---|---|---|---|---|---|---|---|---|---|---
(i) | $-$ | $-$ | (iii) | + | + | (v) | $-$ | + | (vii) | + | +
(ii) | + | $-$ | (iv) | $-$ | $-$ | (vi) | + | $-$ |  |  |
Table D2:  The derivative table of $\Delta w_{t}^{(1)}$ and $\Delta w_{t}^{(2)}$ in the stability diagram ($z_{t}^{(10)}=\delta$; FIG. D1[a])  # area | $\Delta w_{t}^{(1)}$ | $\Delta w_{t}^{(2)}$ | # area | $\Delta w_{t}^{(1)}$ | $\Delta w_{t}^{(2)}$ | # area | $\Delta w_{t}^{(1)}$ | $\Delta w_{t}^{(2)}$ | # area | $\Delta w_{t}^{(1)}$ | $\Delta w_{t}^{(2)}$
---|---|---|---|---|---|---|---|---|---|---|---
(i) | $-$ | $-$ | (iii) | + | + | (v) | $-$ | + | (vii) | + | +
(ii) | + | $-$ | (iv) | $-$ | $-$ | (vi) | + | $-$ |  |  |
Table D3:  The derivative table of $\Delta w_{t}^{(1)}$ and $\Delta w_{t}^{(2)}$ in the stability diagram ($z_{t}^{(10)}<\delta$; FIG. D1[b])  # area | $\Delta w_{t}^{(1)}$ | $\Delta w_{t}^{(2)}$ | # area | $\Delta w_{t}^{(1)}$ | $\Delta w_{t}^{(2)}$ | # area | $\Delta w_{t}^{(1)}$ | $\Delta w_{t}^{(2)}$ | # area | $\Delta w_{t}^{(1)}$ | $\Delta w_{t}^{(2)}$
---|---|---|---|---|---|---|---|---|---|---|---
(i) | $-$ | $-$ | (iii) | $-$ | + | (v) | $-$ | $-$ | (vii) | + | +
(ii) | + | $-$ | (iv) | + | + | (vi) | + | $-$ |  |  |
Figure D1:  Stability diagrams. The solid black and solid red lines represent
$\Delta w_{t}^{(1)}=0$ and $\Delta w_{t}^{(2)}=0$, respectively. The signs of
$\Delta w_{t}^{(1)}$ and $\Delta w_{t}^{(2)}$ in (a) and (b) are shown in
TABLEs D2 and D3, respectively. In (a) and (b), the stability diagrams of
$w_{t}^{(1)}$ and $w_{t}^{(2)}$ when $z_{t}^{(10)}=\delta$ and
$z_{t}^{(10)}>\delta$, respectively, are shown. Note that the intersection
point of the nullclines is
$(w_{t}^{(1)},w_{t}^{(2)})=(z_{t}^{(10)},(1-\alpha)/\beta-\delta/\beta
z_{t}^{(10)})$.

###  D.2 Divergence

Next, we examined the conditions under which the NARMA10 model diverges. Since
the nullclines are time-varying due to $z_{t}^{(10)}$, the increase or
decrease of $w_{t}^{(1)}$ and $w_{t}^{(2)}$ in each area indicated by FIG.
C1(b) and TABLE D1 can change. Thus, in the case of no input ($\mu=\kappa=0$),
we examined the time step at which $y_{t}$ diverges with the initial values of
$w_{0}^{(1)}$ and $w_{0}^{(2)}$. The initial values of $y_{t}$ were set as
$y_{0}=w_{0}^{(1)}$ and
$y_{1}=y_{2}=\cdots=y_{9}=(w_{0}^{(2)}-w_{0}^{(1)})/9$. As shown in FIG.
C1(c), the combination of initial values $(w_{0}^{(1)},w_{0}^{(2)})$ where
$y_{t}$ does not diverge is distributed in a complex manner. We attribute the
difference between the theoretical and numerical basins of attraction to the
change in $z_{t}^{(10)}$, which produced a complicated distribution of the
nullclines (see FIGs. C1[b] and D1, and TABLEs D1–D3).

In the presence of input, $y_{t}$ diverges due to the initial condition. As
shown in FIG. C1(d), the model with input $u_{t}\in[0,0.4]$ ($\mu=\kappa=0.2$)
diverges at an initial value similar to the one in the case of no input. FIG.
C1(e) shows the bifurcation diagram of $y_{t}$ given an initial value at which
$y_{t}$ converges to the fixed-point. Since the fixed-point around which
$y_{t}$ fluctuates is the saddle point, $y_{t}$ can diverge with the given
input. To investigate the divergence conditions caused by the input, we ran
the model over $10^{6}$ time steps and examined the time step at which $y_{t}$
diverges by altering the initial values $y_{0}=\cdots=y_{9}=\psi$ and input
intensity $\sigma$. FIG. C1(f) shows that $y_{t}$ diverges to infinity when
$\psi$ or $\sigma$ exceeds three thresholds: (i) $\psi<-5.15$, (ii)
$1.239<\psi$, and (iii) $\sigma>0.45$. In the (i) and (ii) cases, $y_{t}$
diverges at a shorter time step ($t<100$) than in (iii) because the divergence
is caused by the initial value $\psi$. However, in case of (iii), $y_{t}$
diverges according to the input. After converging to the fixed-point, $y_{t}$
can diverge successively when receiving large positive inputs. According to
the time steps in FIG. C1(f), when $y_{t}$ diverges as the runtime becomes
longer, $y_{t}$ diverges with smaller $\sigma$. Consequently, the threshold
(iii) is the boundary that depends on the input time-series. These results
suggest that $y_{t}$ in the vicinity of the fixed-point can diverge depending
on the input.

In FIG. 2, $p$ expresses the probability of not diverging as a function of
$\sigma$ for different random series $\\{\zeta_{t}\\}$. As previously shown,
once $y_{t}$ converges to the fixed-point, the model stochastically diverges
due to the input, and the probability $p$ depends on $\sigma$; thus, even
though $p$ is high, $y_{t}$ can potentially diverge. For example, two typical
ranges of input, $u_{t}\in[0,0.5]$ and $[0,1]$, have been used for the NARMA10
task; however, $y_{t}$ diverges in both cases (FIG. 2[b], $\sigma=0.5,1.0$)
because $\sigma$ is relevant to the average time until divergence, and the
runtime, $10^{6}$ time steps, is much longer than the time used for the
benchmark task. Consequently, the divergence probability of $y_{t}$, $p$,
changes depending on the parameter $\sigma$.

Therefore, although the NARMA10 model produces the fixed-point attractor, it
can potentially diverge depending on the initial values, input time-series,
and parameter settings.

###  D.3 Time-variance analysis

Finally, we investigated the time-dependence of the NARMA10 model. The state
of a dynamical system receiving noise input can transit from non-chaos to
chaos, an effect referred to as noise-induced chaos [66]. As chaotic behavior
is exhibited by a time-variant system, we investigated whether the system is
chaotic or ordered by calculating the maximum Lyapunov exponent $\lambda_{1}$.
Thus, we derived the Lyapunov spectrum of the model
$\lambda_{i}~{}(i=1,2,\ldots,10)$ based on the 10-dimensional time-delay
system in Eqs. (D1) and (D2). We expressed $z_{t}^{(s)}~{}(s=1,\ldots,10)$ as
a vector $\bm{z}_{t}=[z^{(1)}_{t}\cdots z^{(10)}_{t}]^{\top}$, and the
Jacobian matrix of Eq. (17) ${\bm{J}}_{t}\in\mathbb{R}^{10\times 10}$ with
respect to ${\bm{z}}_{t}$ can be written as follows:

$\displaystyle{\bm{J}}_{t}=\frac{\partial{\bm{z}}_{t+1}}{\partial{\bm{z}}_{t}}=\begin{bmatrix}X_{t}&Y_{t}&\cdots&Y_{t}&Y_{t}\\\
1&&&&\\\ &1&&&\\\ &&\ddots&&\\\ &&&1&\\\ \end{bmatrix},$ (D7)

where $X_{t}=\alpha+2\beta z_{t}^{(1)}+\beta\sum_{i=2}^{10}z_{t}^{(i)}$ and
$Y_{t}=\beta z_{t}^{(1)}$. Using the Jacobian matrices, the Lyapunov spectrum
was computed as follows:

$\displaystyle\lambda_{i}=\frac{1}{T}\sum_{t=1}^{T}\ln{\rho_{i}({\bm{J}}_{t+M-1}{\bm{J}}_{t+M-2}\cdots{\bm{J}}_{t})}~{}(i=1,2,\ldots,10),$
(D8)

where $\rho_{i}({\bm{J}}_{t+M-1}{\bm{J}}_{t+M-2}\cdots{\bm{J}}_{t})$ is the
$i^{\rm th}$ singular value of matrix
${\bm{J}}_{t+M-1}{\bm{J}}_{t+M-2}\cdots{\bm{J}}_{t}$, while $T$ and $M$ were
set to 6000 and 40, respectively. FIG. C1(g) shows the three largest Lyapunov
spectra, $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$, all of which are
negative relative to $\sigma$, indicating that the system does not demonstrate
chaos. Therefore, the NARMA10 model is not a chaotic and time-variant system.

The divergence and time-invariance analysis results revealed that the NARMA10
model converges to the fixed-point and varies in the vicinity of the point. We
considered the $y_{t}$ fluctuating around the fixed-point to be time-
invariant.

##  Appendix E Interconvertibility of the PC expansion and IPC

To clearly demonstrate that the PC expansion and IPC are interconvertible, we
derived an approximate model that has nearly the same breakdown of the IPC as
the original breakdown using Legendre-chaoses. From the above capacity
analysis, we narrowed the polynomial terms to
$P_{1}(\zeta_{t-s})~{}(s=1,2,\ldots)$ and
$P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9})~{}(s=1,2,\ldots)$, which yielded
significantly greater capacities. The expanded state is expressed as follows:

$\displaystyle
y_{t}=p+\sum_{s\in\mathcal{N}_{1}}q_{s}P_{1}(\zeta_{t-s})+\sum_{s\in\mathcal{N}_{2}}r_{s}P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9}),$

where $p$, $q_{s}$, and $r_{s}$ are coefficients for the Legendre-chaoses
$P_{0}=1$, $P_{1}(\zeta_{t-s})$, and $P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9})$,
respectively, while $P_{1}(\zeta)=\zeta$. $\zeta_{t}$ follows a uniform
distribution in $[-1,1]$, and $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$
represent the sets of delayed time steps $s$ for $P_{1}(\zeta_{t-s})$ and
$P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9})$, respectively. Let the normalized
Legendre-chaoses be $\phi_{0}=\frac{1}{\sqrt{T}}$,
$\phi_{1,t}^{(s)}=\frac{P_{1}(\zeta_{t-s})}{\sqrt{\sum_{t=1}^{T}P(\zeta_{t-s})^{2}}}$,
and
$\phi_{2,t}^{(s)}=\frac{P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9})}{\sqrt{\sum_{t=1}^{T}\\{P(\zeta_{t-s})P(\zeta_{t-s-9})\\}^{2}}}$,
and the state is represented as follows:

$\displaystyle
y_{t}=\hat{p}\phi_{0}+\sum_{s\in\mathcal{N}_{1}}\hat{q}_{s}\phi_{1,t}^{(s)}+\sum_{s\in\mathcal{N}_{2}}\hat{r}_{s}\phi_{2,t}^{(s)},$
(E2)

where $\hat{p}=p\sqrt{T}$,
$\hat{q}_{s}=q_{s}\sqrt{\sum_{t=1}^{T}P(\zeta_{t-s})^{2}}$, and
$\hat{r}_{s}=r_{s}\sqrt{\sum_{t=1}^{T}\\{P(\zeta_{t-s})P(\zeta_{t-s-9})\\}^{2}}$
are the modified coefficients for $p$, $q_{s}$, and $r_{s}$, respectively.

Detrending the state and using Eq. (6), the IPCs for $\phi_{1,t}^{(s)}$ and
$\phi_{2,t}^{(s)}$ become

$\displaystyle
C=\frac{q_{s}^{2}}{\sum_{s\in\mathcal{N}_{1}}q_{s}^{2}+\sum_{s\in\mathcal{N}_{2}}r_{s}^{2}}$
(E3)

and

$\displaystyle
C=\frac{r_{s}^{2}}{\sum_{s\in\mathcal{N}_{1}}q_{s}^{2}+\sum_{s\in\mathcal{N}_{2}}r_{s}^{2}},$
(E4)

respectively. Therefore, Eqs. (E3) and (E4) show that each IPC is the
normalized squared coefficient in the polynomial chaos expansion. To
demonstrate this model, we employ $\mathcal{N}_{1}=\\{1,2,3,10,11,12\\}$ and
$\mathcal{N}_{2}=\\{1,2,3\\}$. As shown in FIG. E1(a) and (b), for
$u_{t}\in[-\sigma,\sigma]$ and $[0,\sigma]$, respectively, the approximate
model successfully reproduced the original NARMA10 model. Furthermore, as in
FIG. E1(c) and (d), we confirmed that the IPC breakdown of the approximate
model reproduced the original breakdown.

Figure E1:  The proposed model for the benchmark task. Frames (a) and (b) show
the time-series data with $u_{t}\in[-\sigma,\sigma]$ and $[0,\sigma]$,
respectively. The solid black and solid red lines represent the output of the
NARMA10 and our models, respectively. Panels (c) and (d) show the IPC with
$u_{t}\in[-\sigma,\sigma]$ and $[0,\sigma]$, respectively. The labels
represent combinations of $\\{\\{n_{s},s\\}\\}$, where $n_{s}$ is the degree
of the polynomial, and $s$ is the delayed time step of the input. Here, the
desired output is $\prod_{s}P_{n_{s}}(\zeta_{t-s})$.

###  E.1 Derivation

Here, we derive the following equations:

$\displaystyle y_{t}$ $\displaystyle=$ $\displaystyle
p+\sum_{s\in\mathcal{N}_{1}}q_{s}P_{1}(\zeta_{t-s})$ (E5) $\displaystyle\ \
+\sum_{s\in\mathcal{N}_{2}}r_{s}P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9}),$
$\displaystyle p$ $\displaystyle=$
$\displaystyle\frac{1-\alpha}{20\beta}-\sqrt{\left(\frac{1-\alpha}{20\beta}\right)^{2}-\frac{\gamma\mu^{2}+\delta}{10\beta}},$
(E6) $\displaystyle q_{s}$ $\displaystyle=$
$\displaystyle\begin{cases}\gamma\mu\kappa~{}(s=1)\\\ (\alpha+10\beta
p)q_{s-1}+\sum_{j=0}^{s-2}\beta pq_{s-j-1}\\\ (s=2,3,\ldots,9)\\\
\gamma\mu\kappa+(\alpha+10\beta p)q_{s-1}+\sum_{j=0}^{s-2}\beta pq_{s-j-1}\\\
(s=10)\\\ (\alpha+10\beta p)q_{s-1}+\sum_{j=0}^{9}\beta pq_{s-j-1}\\\
(s=11,12,\ldots)\end{cases},$ $\displaystyle r_{s}$ $\displaystyle=$
$\displaystyle\begin{cases}\gamma\kappa^{2}~{}(s=1)\\\ (\alpha+10\beta
p)r_{s-1}+\beta q_{s-1}\sum_{j=0}^{9}q_{s+8-j}\\\
+\beta\sum_{j=0}^{s-2}(pr_{s-j-1}+q_{s+8}q_{s-j-1})\\\ (s=2,3,\ldots,10)\\\
(\alpha+10\beta p)r_{s-1}+\beta q_{s-1}\sum_{j=0}^{9}q_{s+8-j}\\\
+\beta\sum_{j=0}^{9}(pr_{s-j-1}+q_{s+8}q_{s-j-1})\\\
(s=11,12,\ldots)\end{cases},$ (E8)

where $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$ represent the sets of delayed
time steps $s\in\mathbb{N}$ for $P_{1}(\zeta_{t-s})$ and
$P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9})$, respectively.

First, $y_{t}$ in Eq. (17) is expanded using the Legendre-chaoses of input
time-series $\zeta_{t-s}~{}(s=1,2,\ldots,t)$ with time-varying coefficients as
follows:

$\displaystyle y_{t}$ $\displaystyle=$ $\displaystyle
p_{t}+\sum_{s=1}^{t}q_{t,s}P_{1}(\zeta_{t-s})$ (E9) $\displaystyle+$
$\displaystyle\sum_{s=1}^{t}r_{t,s}P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9})+\cdots,$

where $p_{t}$ denotes a time-varying term independent of
$\zeta_{t-s}~{}(s=1,2,\ldots,t)$, and $q_{t,s}$ and $r_{t,s}$ are the $s^{\rm
th}$ coefficients of $P_{1}(\zeta_{t-s})$ and
$P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9})$, respectively. The NARMA10 model with
$\zeta_{t}$ can be expressed as

$\displaystyle y_{t+1}$ $\displaystyle=$ $\displaystyle\alpha y_{t}+\beta
y_{t}\sum_{j=0}^{9}y_{t-j}$ (E10) $\displaystyle+$
$\displaystyle\gamma(\mu+\sigma\zeta_{t})(\mu+\sigma\zeta_{t-9})+\delta.$

According to Eqs. (E9) and (E10), $y_{t+1}$ is rewritten as

$\displaystyle y_{t+1}$ $\displaystyle=$ $\displaystyle\left(\alpha
p_{t}+\beta p_{t}\sum_{j=0}^{9}p_{t-j}+\gamma\mu^{2}+\delta\right)$ (E11)
$\displaystyle+$
$\displaystyle\gamma\mu\kappa\left(P_{1}(\zeta_{t})+P_{1}(\zeta_{t-9})\right)$
$\displaystyle+$ $\displaystyle\alpha\sum_{s=1}^{t}q_{t,s}P_{1}(\zeta_{t-s})$
$\displaystyle+$ $\displaystyle\beta
p_{t}\sum_{j=0}^{9}\sum_{s=1}^{t-j}q_{t-j,s}P_{1}(\zeta_{t-s-j})$
$\displaystyle+$
$\displaystyle\beta\sum_{j=0}^{9}p_{t-j}\sum_{s=1}^{t}q_{t,s}P_{1}(\zeta_{t-s})$
$\displaystyle+$ $\displaystyle\gamma\kappa^{2}\zeta_{t}\zeta_{t-9}$
$\displaystyle+$
$\displaystyle\alpha\sum_{s=1}^{t}r_{t,s}P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9})$
$\displaystyle+$ $\displaystyle\beta
p_{t}\sum_{j=0}^{9}\sum_{s=1}^{t-j}r_{t-j,s}P_{1}(\zeta_{t-s-j})P_{1}(\zeta_{t-s-j-9})$
$\displaystyle+$
$\displaystyle\beta\sum_{j=0}^{9}p_{t-j}\sum_{s=1}^{t}r_{t,s}P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9})$
$\displaystyle+$
$\displaystyle\sum_{s=1}^{t}q_{t,s}P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s})\sum_{j=0}^{9}\sum_{s=1}^{t-j}q_{t-j,s}P_{1}(\zeta_{t-s-j})$
$\displaystyle+$ $\displaystyle\cdots,$

where $P_{1}(\zeta)=\zeta$. When increasing $t$ by one in Eq. (E9), the
following equation is obtained:

$\displaystyle y_{t+1}$ $\displaystyle=$ $\displaystyle
p_{t+1}+\sum_{s=1}^{t+1}q_{t+1,s}P_{1}(\zeta_{t+1-s})$ $\displaystyle+$
$\displaystyle\sum_{s=1}^{t+1}r_{t+1,s}P_{1}(\zeta_{t+1-s})P_{1}(\zeta_{t+1-s-9})+\cdots.$

Equating the coefficients in Eqs. (E11) and (LABEL:eqA:xt+1_2) yields:

$\displaystyle p_{t+1}$ $\displaystyle=$ $\displaystyle\alpha
p_{t}+\sum_{j=0}^{9}\beta p_{t}p_{t-j}+\gamma\mu^{2}+\delta,$ (E13)
$\displaystyle q_{t+1,s}$ $\displaystyle=$
$\displaystyle\begin{cases}\gamma\mu\kappa~{}(s=1)\\\
\left(\alpha+\beta\sum_{j=0}^{9}p_{t-j}\right)q_{t,s-1}\\\
+\sum_{j=0}^{s-2}\beta p_{t}q_{t-j,s-j-1}~{}(s=2,3,\ldots,9)\\\
\gamma\mu\kappa+\left(\alpha+\beta\sum_{j=0}^{9}p_{t-j}\right)q_{t,s-1}\\\
+\sum_{j=0}^{s-2}\beta p_{t}q_{t-j,s-j-1}~{}(s=10)\\\
\left(\alpha+\beta\sum_{j=0}^{9}p_{t-j}\right)q_{t,s-1}\\\
+\sum_{j=0}^{9}\beta p_{t}q_{t-j,s-j-1}~{}(s=11,12,\ldots)\end{cases},$
$\displaystyle r_{t+1,s}$ $\displaystyle=$
$\displaystyle\begin{cases}\gamma\kappa^{2}~{}(s=1)\\\
\left(\alpha+\beta\sum_{j=0}^{9}p_{t-j}\right)r_{t,s-1}\\\ +\beta
q_{t,s-1}\sum_{j=0}^{9}q_{t-j,s+8-j}\\\
+\beta\sum_{j=0}^{s-2}(p_{t}r_{t-j,s-j-1}+q_{t,s+8}q_{t-j,s-j-1})\\\
(s=2,3,\ldots,10)\\\
\left(\alpha+\beta\sum_{j=0}^{9}p_{t-j}\right)r_{t,s-1}\\\ +\beta
q_{t,s-1}\sum_{j=0}^{9}q_{t-j,s+8-j}\\\
+\beta\sum_{j=0}^{9}(p_{t}r_{t-j,s-j-1}+q_{t,s+8}q_{t-j,s-j-1})\\\
(s=11,12,\ldots)\end{cases}.$

According to Eq. (E13), $p_{t}$ has a stable and an unstable equilibrium
point. If
$p_{t}<(1-\alpha)/20\beta+\sqrt{((1-\alpha)/20\beta)^{2}-(\gamma\mu^{2}+\delta)/10\beta}$,
it converges to the stable point. When $t$ is large enough, $p_{t}$ converges
to

$\displaystyle
p=\lim_{t\rightarrow\infty}p_{t}=\frac{1-\alpha}{20\beta}-\sqrt{\left(\frac{1-\alpha}{20\beta}\right)^{2}-\frac{\gamma\mu^{2}+\delta}{10\beta}}.$
(E16)

According to Eq. (LABEL:eqA:qt+1), $q_{t,s}$ also converges to

$\displaystyle q_{s}$ $\displaystyle=$
$\displaystyle\lim_{t\rightarrow\infty}q_{t+1,s}$ $\displaystyle=$
$\displaystyle\begin{cases}\gamma\mu\kappa~{}(s=1)\\\ (\alpha+10\beta
p)q_{s-1}+\sum_{j=0}^{s-2}\beta pq_{s-j-1}\\\ (s=2,3,\ldots,9)\\\
\gamma\mu\kappa+(\alpha+10\beta p)q_{s-1}+\sum_{j=0}^{s-2}\beta pq_{s-j-1}\\\
(s=10)\\\ (\alpha+10\beta p)q_{s-1}+\sum_{j=0}^{9}\beta pq_{s-j-1}\\\
(s=11,12,\ldots)\end{cases}.$

In the same manner, $r_{t,s}$ converges to

$\displaystyle r_{s}$ $\displaystyle=$
$\displaystyle\lim_{t\rightarrow\infty}r_{t+1,s}$ $\displaystyle=$
$\displaystyle\begin{cases}\gamma\kappa^{2}~{}(s=1)\\\ (\alpha+10\beta
p)r_{s-1}+\beta q_{s-1}\sum_{j=0}^{9}q_{s+8-j}\\\
+\beta\sum_{j=0}^{s-2}(pr_{s-j-1}+q_{s+8}q_{s-j-1})\\\ (s=2,3,\ldots,10)\\\
(\alpha+10\beta p)r_{s-1}+\beta q_{s-1}\sum_{j=0}^{9}q_{s+8-j}\\\
+\beta\sum_{j=0}^{9}(pr_{s-j-1}+q_{s+8}q_{s-j-1})\\\
(s=11,12,\ldots)\end{cases}.$

Therefore, when $t$ is large enough, and Eq. (E9) is approximated with the
constant term $p$ and the Legendre-chaoses of $P_{1}(\zeta_{t-s})$ and
$P_{1}(\zeta_{t-s})P_{1}(\zeta_{t-s-9})$, whose delayed time steps $s$ are
limited to sets $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$, respectively, Eqs.
(E5)–(E8) are obtained.

##  Appendix F How to compose the dissociated culture reservoir

All experiments were approved by the ethical committee of the University of
Tokyo and followed the “Guiding Principles for the Care and Use of Animals in
the Field of Physiological Science” estalished by the Physiological Society of
Japan. Embryonic rat cortices were dissected from E18 rats and used for
cortical cell cultures. The cortices were dissociated in 2 mL of 0.25%
trypsin-ethylenediaminetetraacetic acid (Trypsin-EDTA, Life Technologies),
from which cells were isolated by trituration, and 38,000 cells were seeded on
each microelectrode array (MEA; MaxWell Biosystems). For cell adhesion, 5 mL
of 0.05% Polyethileneimine (PEI; Sigma-Aldrich) and 5 ${\rm\mu}$l of 0.02
mg/ml Laminin (Sigma-Aldrich) were used before plating the cells. Then, after
24 hours, the plating media[67] were changed to growth media[68]. The plating
media were composed of Neurobasal $850\ {\rm\mu l}$ (Life Technologies), 10%
horse serum (HyClone), 0.5 mM GlutaMAX (Life Technologies), and 2% B27 (Life
Technologies). The growth media were composed of DMEM $850{\rm\mu l}$ (Life
Technologies), 10%horse serum (HyClone), 0.5 mM GlutaMAX (Life Technologies),
and 1 mM sodium pyruvate (Life Technologies). All experiments were conducted
in an incubator at 37∘C and 5% CO2. The MEAs were sealed with a lid to prevent
water evaporation and invasion of bacteria and fungus.

The MEA had 26,400 electrodes, which were placed 17.5 $\rm\mu$m apart and
arranged in a 120 $\times$ 220 grid. The MEA can simultaneously use up to
1,024 of 26,400 electrodes. We selected electrodes with a high firing rate as
measurement electrodes and electrodes on which the axon places, as stimulation
electrodes. We applied bipolar pulse stimuli with an amplitude of $\zeta_{t}$,
which followed a normal distribution with mean $\mu$ and standard deviation
$\sigma$, and of an interpulse interval (IPI) of 10, 20, and 30 ms to the
stimulation electrodes. Furthermore, a 6th-order Butterworth bandpass filter
and zero-phase IIR filter were applied to the voltage traces observed from the
measurement electrodes to extract 300–3000 Hz components. At all electrodes,
stimulus-induced artifacts were removed by eliminating traces $\pm$2 ms from
the stimulus times. The standard deviation of extracted signals was calculated
as follows [69]:

$\displaystyle\sigma={\rm
median}\biggl{\\{}\frac{|{\bm{x}}|}{0.6745}\biggr{\\}}.$ (F1)

If the amplitude of an extracted signal exceeded $4\sigma$, the value of the
spike train was set to one; otherwise, it was set to zero. As the measurement
frequency was 20 kHz, the above spike train was separated by a 1-ms time bin,
and if one or more spikes appeared in one bin, the modified spike train was
set to one; otherwise it was set to zero. The train was divided into bins by
IPI-width, and the number of spikes in the bin was used for the state
${\bm{x}}_{t}$.

## References

  * Stieg  _et al._ [2012] A. Z. Stieg, A. V. Avizienis, H. O. Sillin, C. Martin-Olmos, M. Aono, and J. K. Gimzewski, Emergent criticality in complex turing b-type atomic switch networks, Advanced Materials 24, 286 (2012).
  * Jaeger and Haas [2004] H. Jaeger and H. Haas, Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication, Science 304, 78 (2004).
  * Maass  _et al._ [2002] W. Maass, T. Natschläger, and H. Markram, Real-time computing without stable states: A new framework for neural computation based on perturbations, Neural computation 14, 2531 (2002).
  * Verstraeten  _et al._ [2007] D. Verstraeten, B. Schrauwen, M. d’Haene, and D. Stroobandt, An experimental unification of reservoir computing methods, Neural networks 20, 391 (2007).
  * Appeltant  _et al._ [2011] L. Appeltant, M. C. Soriano, G. Van der Sande, J. Danckaert, S. Massar, J. Dambre, B. Schrauwen, C. R. Mirasso, and I. Fischer, Information processing using a single dynamical node as complex system, Nature communications 2, 468 (2011).
  * Brunner  _et al._ [2013] D. Brunner, M. C. Soriano, C. R. Mirasso, and I. Fischer, Parallel photonic information processing at gigabyte per second data rates using transient states, Nature communications 4, 1 (2013).
  * Vandoorne  _et al._ [2014] K. Vandoorne, P. Mechet, T. Van Vaerenbergh, M. Fiers, G. Morthier, D. Verstraeten, B. Schrauwen, J. Dambre, and P. Bienstman, Experimental demonstration of reservoir computing on a silicon photonics chip, Nature communications 5, 1 (2014).
  * Du  _et al._ [2017] C. Du, F. Cai, M. A. Zidan, W. Ma, S. H. Lee, and W. D. Lu, Reservoir computing using dynamic memristors for temporal information processing, Nature communications 8, 2204 (2017).
  * Torrejon  _et al._ [2017] J. Torrejon, M. Riou, F. A. Araujo, S. Tsunegi, G. Khalsa, D. Querlioz, P. Bortolotti, V. Cros, K. Yakushiji, A. Fukushima, _et al._ , Neuromorphic computing with nanoscale spintronic oscillators, Nature 547, 428 (2017).
  * Moon  _et al._ [2019] J. Moon, W. Ma, J. H. Shin, F. Cai, C. Du, S. H. Lee, and W. D. Lu, Temporal data classification and forecasting using a memristor-based reservoir computing system, Nature Electronics 2, 480 (2019).
  * Marinella and Agarwal [2019] M. J. Marinella and S. Agarwal, Efficient reservoir computing with memristors, Nature Electronics 2, 437 (2019).
  * Lüdge and Röhm [2019] K. Lüdge and A. Röhm, Computing with a camera, Nature Machine Intelligence 1, 551 (2019).
  * Jaeger [2002] H. Jaeger,  _Tutorial on training recurrent neural networks, covering BPPT, RTRL, EKF and the ”echo state network” approach_ , Vol. 5 (GMD-Forschungszentrum Informationstechnik Bonn, 2002).
  * Yildiz  _et al._ [2012] I. B. Yildiz, H. Jaeger, and S. J. Kiebel, Re-visiting the echo state property, Neural networks 35, 1 (2012).
  * Manjunath and Jaeger [2013] G. Manjunath and H. Jaeger, Echo state property linked to an input: Exploring a fundamental characteristic of recurrent neural networks, Neural computation 25, 671 (2013).
  * Maass and Markram [2004] W. Maass and H. Markram, On the computational power of circuits of spiking neurons, Journal of computer and system sciences 69, 593 (2004).
  * Maass [2011] W. Maass, Liquid state machines: motivation, theory, and applications, in  _Computability in context: computation and logic in the real world_ (World Scientific, 2011) pp. 275–296.
  * Nakajima [2020] K. Nakajima, Physical reservoir computing—–an introductory perspective, Japanese Journal of Applied Physics 59, 060501 (2020).
  * Dambre  _et al._ [2012] J. Dambre, D. Verstraeten, B. Schrauwen, and S. Massar, Information processing capacity of dynamical systems, Scientific reports 2, 514 (2012).
  * Martínez-Peña  _et al._ [2020] R. Martínez-Peña, J. Nokkala, G. Giorgi, R. Zambrini, and M. Soriano, Information processing capacity of spin-based quantum reservoir computing systems, Cognitive Computation , 1 (2020).
  * Nokkala  _et al._ [2020] J. Nokkala, R. Martínez-Peña, G. L. Giorgi, V. Parigi, M. C. Soriano, and R. Zambrini, Gaussian states provide universal and versatile quantum reservoir computing, arXiv preprint arXiv:2006.04821 (2020).
  * Akashi  _et al._ [2020] N. Akashi, T. Yamaguchi, S. Tsunegi, T. Taniguchi, M. Nishida, R. Sakurai, Y. Wakao, and K. Nakajima, Input-driven bifurcations and information processing capacity in spintronics reservoirs, Phys. Rev. Research 2, 043303 (2020).
  * Wiener [1938] N. Wiener, The homogeneous chaos, American Journal of Mathematics 60, 897 (1938).
  * Xiu and Karniadakis [2002] D. Xiu and G. E. Karniadakis, The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM journal on scientific computing 24, 619 (2002).
  * Oladyshkin and Nowak [2012] S. Oladyshkin and W. Nowak, Data-driven uncertainty quantification using the arbitrary polynomial chaos expansion, Reliability Engineering & System Safety 106, 179 (2012).
  * Jaeger [2003] H. Jaeger, Adaptive nonlinear system identification with echo state networks, in  _Advances in neural information processing systems_ (2003) pp. 609–616.
  * Atiya and Parlos [2000] A. F. Atiya and A. G. Parlos, New results on recurrent network training: unifying the algorithms and accelerating convergence, IEEE transactions on neural networks 11, 697 (2000).
  * Verstraeten and Schrauwen [2009] D. Verstraeten and B. Schrauwen, On the quantification of dynamics in reservoir computing, in  _International Conference on Artificial Neural Networks_ (Springer, 2009) pp. 985–994.
  * Rodan and Tino [2010] A. Rodan and P. Tino, Minimum complexity echo state network, IEEE transactions on neural networks 22, 131 (2010).
  * Nakajima  _et al._ [2019] K. Nakajima, K. Fujii, M. Negoro, K. Mitarai, and M. Kitagawa, Boosting computational power through spatial multiplexing in quantum reservoir computing, Physical Review Applied 11, 034021 (2019).
  * Paquot  _et al._ [2012] Y. Paquot, F. Duport, A. Smerieri, J. Dambre, B. Schrauwen, M. Haelterman, and S. Massar, Optoelectronic reservoir computing, Scientific reports 2, 287 (2012).
  * Dale  _et al._ [2016] M. Dale, J. F. Miller, S. Stepney, and M. A. Trefzer, Evolving carbon nanotube reservoir computers, in  _International Conference on Unconventional Computation and Natural Computation_ (Springer, 2016) pp. 49–61.
  * Nakajima  _et al._ [2013] K. Nakajima, H. Hauser, R. Kang, E. Guglielmino, D. G. Caldwell, and R. Pfeifer, Computing with a muscular-hydrostat system, in  _2013 IEEE International Conference on Robotics and Automation_ (IEEE, 2013) pp. 1504–1511.
  * Nakajima  _et al._ [2018] K. Nakajima, H. Hauser, T. Li, and R. Pfeifer, Exploiting the dynamics of soft materials for machine learning, Soft robotics 5, 339 (2018).
  * Bürger  _et al._ [2015] J. Bürger, A. Goudarzi, D. Stefanovic, and C. Teuscher, Hierarchical composition of memristive networks for real-time computing, in  _Proceedings of the 2015 IEEE/ACM International Symposium on Nanoscale Architectures (NANOARCH´ 15)_ (IEEE, 2015) pp. 33–38.
  * Fujii and Nakajima [2017] K. Fujii and K. Nakajima, Harnessing disordered-ensemble quantum dynamics for machine learning, Physical Review Applied 8, 024030 (2017).
  * Hermans  _et al._ [2016] M. Hermans, P. Antonik, M. Haelterman, and S. Massar, Embodiment of learning in electro-optical signal processors, Physical review letters 117, 128301 (2016).
  * Yin  _et al._ [2012] J. Yin, Y. Meng, and Y. Jin, A developmental approach to structural self-organization in reservoir computing, IEEE transactions on autonomous mental development 4, 273 (2012).
  * Okumura  _et al._ [2019] T. Okumura, M. Tai, and M. Ando, Experimental study on parallel and analog optical reservoir computing with delayed feedback system for physical implementation, Nonlinear Theory and Its Applications, IEICE 10, 236 (2019).
  * Barazani  _et al._ [2020] B. Barazani, G. Dion, J.-F. Morissette, L. Beaudoin, and J. Sylvestre, Microfabricated neuroaccelerometer: Integrating sensing and reservoir computing in mems, Journal of Microelectromechanical Systems (2020).
  * Bianchi  _et al._ [2017] F. M. Bianchi, L. Livi, C. Alippi, and R. Jenssen, Multiplex visibility graphs to investigate recurrent neural network dynamics, Scientific reports 7, 44037 (2017).
  * Duport  _et al._ [2016] F. Duport, A. Smerieri, A. Akrout, M. Haelterman, and S. Massar, Fully analogue photonic reservoir computer, Scientific reports 6, 22381 (2016).
  * Hermans and Schrauwen [2012] M. Hermans and B. Schrauwen, Recurrent kernel machines: Computing with infinite echo state networks, Neural Computation 24, 104 (2012).
  * Inubushi and Yoshimura [2017] M. Inubushi and K. Yoshimura, Reservoir computing beyond memory-nonlinearity trade-off, Scientific reports 7, 10199 (2017).
  * Kan  _et al._ [2021] S. Kan, K. Nakajima, Y. Takeshima, T. Asai, Y. Kuwahara, and M. Akai-Kasaya, Simple reservoir computing capitalizing on the nonlinear response of materials: Theory and physical implementations, Physical Review Applied 15, 024030 (2021).
  * Tran and Nakajima [2020] Q. H. Tran and K. Nakajima, Higher-order quantum reservoir computing, arXiv preprint arXiv:2006.08999 (2020).
  * Torres  _et al._ [2019] E. Torres, K. Nakajima, and I. S. Godage, Information processing capability of soft continuum arms, in  _2019 2nd IEEE International Conference on Soft Robotics (RoboSoft)_ (IEEE, 2019) pp. 441–447.
  * Strogatz [2001] S. Strogatz, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering (studies in nonlinearity), (2001).
  * Jaeger [2001] H. Jaeger,  _Short term memory in echo state networks_ , Vol. 5 (GMD-Forschungszentrum Informationstechnik, 2001).
  * White  _et al._ [2004] O. L. White, D. D. Lee, and H. Sompolinsky, Short-term memory in orthogonal neural networks, Physical review letters 92, 148102 (2004).
  * Sussillo and Abbott [2009] D. Sussillo and L. F. Abbott, Generating coherent patterns of activity from chaotic neural networks, Neuron 63, 544 (2009).
  * Laje and Buonomano [2013] R. Laje and D. V. Buonomano, Robust timing and motor patterns by taming chaos in recurrent neural networks, Nature neuroscience 16, 925 (2013).
  * Nicola and Clopath [2017] W. Nicola and C. Clopath, Supervised learning in spiking neural networks with force training, Nature communications 8, 1 (2017).
  * Inoue  _et al._ [2020] K. Inoue, K. Nakajima, and Y. Kuniyoshi, Designing spontaneous behavioral switching via chaotic itinerancy, Science Advances 6 (2020).
  * Ahlfeld  _et al._ [2016] R. Ahlfeld, B. Belkouchi, and F. Montomoli, Samba: sparse approximation of moment-based arbitrary polynomial chaos, Journal of Computational Physics 320, 1 (2016).
  * Faisal  _et al._ [2008] A. A. Faisal, L. P. Selen, and D. M. Wolpert, Noise in the nervous system, Nature reviews neuroscience 9, 292 (2008).
  * Goto  _et al._ [2020] K. Goto, K. Nakajima, and H. Notsu, Computing with vortices: Bridging fluid dynamics and its information-processing capability, arXiv preprint arXiv:2001.08502 (2020).
  * Maritan and Banavar [1994] A. Maritan and J. R. Banavar, Chaos, noise, and synchronization, Physical Review Letters 72, 1451 (1994).
  * Toral  _et al._ [2001] R. Toral, C. R. Mirasso, E. Hernández-Garcıa, and O. Piro, Analytical and numerical studies of noise-induced synchronization of chaotic systems, Chaos: An Interdisciplinary Journal of Nonlinear Science 11, 665 (2001).
  * Lu  _et al._ [2018] Z. Lu, B. R. Hunt, and E. Ott, Attractor reconstruction by machine learning, Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 061104 (2018).
  * Boyd and Chua [1985] S. Boyd and L. Chua, Fading memory and the problem of approximating nonlinear operators with volterra series, IEEE Transactions on circuits and systems 32, 1150 (1985).
  * Volterra [1959] V. Volterra, Theory of functionals and of integral and integro-differential equations, (1959).
  * Askey and Wilson [1985] R. Askey and J. A. Wilson,  _Some basic hypergeometric orthogonal polynomials that generalize Jacobi polynomials_ , Vol. 319 (American Mathematical Soc., 1985).
  * Cameron and Martin [1947] R. H. Cameron and W. T. Martin, The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals, Annals of Mathematics , 385 (1947).
  * Ogura [1972] H. Ogura, Orthogonal functionals of the poisson process, IEEE Transactions on Information Theory 18, 473 (1972).
  * Crutchfield  _et al._ [1982] J. P. Crutchfield, J. D. Farmer, and B. A. Huberman, Fluctuations and simple chaotic dynamics, Physics Reports 92, 45 (1982).
  * Brewer  _et al._ [1993] G. J. Brewer, J. Torricelli, E. Evege, and P. Price, Optimized survival of hippocampal neurons in b27-supplemented neurobasal™, a new serum-free medium combination, Journal of neuroscience research 35, 567 (1993).
  * Potter and DeMarse [2001] S. M. Potter and T. B. DeMarse, A new approach to neural cell culture for long-term studies, Journal of neuroscience methods 110, 17 (2001).
  * Quiroga  _et al._ [2004] R. Q. Quiroga, Z. Nadasdy, and Y. Ben-Shaul, Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering, Neural computation 16, 1661 (2004).