See pages 1-last of LeDA-SWAn.pdf


# Learning noise-induced transitions by multi-scaling reservoir computing

Zequn Lin  Faculty of Arts and Sciences, Beijing Normal University, Zhuhai
519087, China  Zhaofan Lu  Faculty of Arts and Sciences, Beijing Normal
University, Zhuhai 519087, China  Zengru Di  International Academic Center of
Complex Systems, Beijing Normal University, Zhuhai 519087, China  Ying Tang
<EMAIL_ADDRESS>International Academic Center of Complex Systems,
Beijing Normal University, Zhuhai 519087, China

###### Abstract

Noise is usually regarded as adversarial to extract the effective dynamics
from time series, such that the conventional data-driven approaches usually
aim at learning the dynamics by mitigating the noisy effect. However, noise
can have a functional role of driving transitions between stable states
underlying many natural and engineered stochastic dynamics. To capture such
stochastic transitions from data, we find that leveraging a machine learning
model, reservoir computing as a type of recurrent neural network, can learn
noise-induced transitions. We develop a concise training protocol for tuning
hyperparameters, with a focus on a pivotal hyperparameter controlling the time
scale of the reservoir dynamics. The trained model generates accurate
statistics of transition time and the number of transitions. The approach is
applicable to a wide class of systems, including a bistable system under a
double-well potential, with either white noise or colored noise. It is also
aware of the asymmetry of the double-well potential, the rotational dynamics
caused by non-detailed balance, and transitions in multi-stable systems. For
the experimental data of protein folding, it learns the transition time
between folded states, providing a possibility of predicting transition
statistics from a small dataset. The results demonstrate the capability of
machine-learning methods in capturing noise-induced phenomena.

##  I Introduction

Noise-induced transitions are ubiquitous in nature and occur in diverse
systems with multi-stable states [1]. Examples include switches between
different voltage and current states in the circuit [2], noisy genetic
switches [3], noise-induced biological homochirality of early life self-
replicators [4], protein conformational transitions [5, 6], and chemical
reactions [7] with the multi-stable probability distribution [8]. Learning
noise-induced transitions is vital for understanding critical phenomena of
these systems. In many scenarios, only time series are available without
mathematical equations known in prior. To effectively learn and predict noise-
induced transitions from time series, there is also a challenge of discerning
dynamics with both slow and fast time scales: fast relaxation around distinct
stable states and slow transitions between them, where the fast time-scale
signals are often referred to noise [9, 10]. Consequently, it remains elusive
to learn stochastic transitions from time series in general.

Recently, many efforts have been made to learn the dynamics from data by
machine-learning methods [11, 12, 13, 14, 15]. One type of approach uses the
automatic differentiation for identifying nonlinear dynamics, denoising time-
series data, and parameterizing the noisy probability distribution from data
[16]. Due to the non-convexity of the optimization problem, the method may
struggle to robustly handle large function libraries for the regression.
Another type of approach employs physics-informed neural networks for data-
driven solutions and discoveries of partial differential equations [17, 18],
or Koopman eigenfunctions from data [19]. However, the method requires an
extensive quantity of data to train the deep neural network, alongside precise
adjustment and refinement of the network. Despite the broad application of the
aforementioned methods, to our knowledge, they have not been utilized in
studying noise-induced transitions.

To investigate whether machine-learning methods can capture and predict noise-
induced transitions, we start with one machine-learning architecture,
reservoir computing (RC) [20, 21]. The training of reservoir computer is a
simple linear regression, which is less computationally expensive than the
neural network that requires the back propagation. The input layer of the
reservoir transforms time series into the space of the reservoir network,
while the output layer transforms the variables of the reservoir back to time
series. The output layer is trained to minimize the difference between the
input and output, by tuning the hyperparameters. The reservoir computing is
particularly effective for learning dynamical systems [22, 23], including
chaotic systems [24, 25, 26, 27]. A recent research started to apply the
reservoir computing to stochastic resonance [28], however, the functional role
of noise in shifting dynamics between stable states has not been investigated.
There is one previous attempt on employing reservoir computing for noise-
induced transitions [29]. Nevertheless, it relies on an impractical assumption
on knowing the equation for the deterministic part of dynamics, and then
employs reservoir computer to learn the separated fast time-scale series. In
practice, one usually lacks a prior knowledge about the deterministic
dynamics, and in some cases, this part even cannot be directly described by an
equation [6]. Thus, the question remains that can we forecast noise-induced
transitions solely based on data without any prior knowledge of the underlying
equation?

In this paper, we develop a framework of multi-scaling reservoir computing for
learning noise-induced transitions in a model-free manner. Our method is
inspired by the capability of reservoir computer to model dynamical systems
[20, 22, 30, 31, 32, 21, 33, 34], and especially the hyperparameter $\alpha$
in the reservoir was found to determine the time scale of reservoir dynamics
[35]. Given a multi-scale time series, we can thus tune the hyperparameter
$\alpha$ to match the slowly time-scale dynamics. After the reservoir captures
the slowly time-scale dynamics by fitting the output layer matrix, we can
separate the fast time-scale series as a noise distribution. During the
predicting phase, we utilize the trained reservoir computer to simulate the
slowly time-scale dynamics, and then add back the noise sampled from the
separated noise distribution (for white noise) or learnt from the second
reservoir (for colored noise). The whole protocol is iterated over time points
as the rolling prediction. Notably, the present method is different from
previous work that regards noise merely as a disturbance [16, 36], and instead
focuses on capturing noise-induced transitions from the data.

Figure 1:  Framework of learning noise-induced transitions by multi-scaling
reservoir computing. (a) The training data is a time series $\mathbf{u}$ with
slow and fast time scales, and the fast time-scale part can be considered as
noise, causing noise-induced transitions between stable states. In the
training phase, at each time step $t$, the reservoir takes into
$\mathbf{u}_{t}$ through a matrix $W_{in}$ and has reservoir state
$\mathbf{r}_{t}$ with a connection matrix $A$. The output matrix $W_{out}$ is
trained to fit the output time series to the training data at the next time
point. Tuning the hyperparameter $\alpha$ alters the time scale of the output
$\widetilde{\mathbf{u}}$, and a properly chosen $\alpha$ leads to a match with
the slowly time-scale data. Then, $\mathbf{u}-\widetilde{\mathbf{u}}$
separates the fast time-scale signal $\mathbf{\eta}$ as a noise distribution.
(b) In the predicting phase, the $\widehat{\mathbf{u}}_{s}$ is put into the
trained reservoir to generate $\widetilde{\mathbf{u}}_{s+1}$. In the next time
step $s+1$, the input $\widehat{\mathbf{u}}_{s+1}$ is the
$\widetilde{\mathbf{u}}_{s+1}$ plus the noise $\mathbf{\eta}_{s+1}$ sampled
from the separated noise distribution. This process is iterated as a rolling
prediction [37] to generate the time series $\widehat{\mathbf{u}}$. (c) The
evaluation on the predicted transition statistics. In the middle, different
colored lines of $\widehat{\mathbf{u}}$ represent replicates of the
predictions. The accuracy is evaluated by the statistics of transition time
and the number of transitions. PDF: probability density function.

To demonstrate the effectiveness of the present method, we apply it to two
categories of scenarios. One type has the data generated from stochastic
differential equations (SDE), and the other has the experimental data [6]. For
the first category with white noise, it includes a one-dimensional (1D)
bistable gradient system, two-dimensional (2D) bistable gradient and non-
gradient systems [38], 1D and 2D gradient systems with a tilted potential, a
2D tilted non-gradient system, and a 2D tristable system [39]. The present
approach can capture statistics of the transition time and the number of
transitions. For the first category with colored noise, we study a 1D bistable
gradient system with Lorenz noise [29], and accurately predict the specific
transition time, without the assumption of knowing the deterministic part of
dynamics as required in [29]. For the second category, we apply the approach
to the protein folding data [6], and explore the least amount of required data
for accurate training, which can help reduce the demand for measurements in
experiments.

##  II Results

###  II.1 The framework of multi-scaling reservoir computing

The reservoir computing has the following scheme [22, 24]:

$\displaystyle\mathbf{r}_{t+1}$
$\displaystyle=(1-\alpha)\mathbf{r}_{t}+\alpha\tanh(A\mathbf{r}_{t}+W_{in}\mathbf{u}_{t}),$
(1) $\displaystyle\mathbf{\widetilde{u}}_{t+1}$
$\displaystyle=W_{out}\mathbf{r}_{t+1}.$ (2)

Here, the vector $\mathbf{u}$ is an $n$-dimensional state vector, and the
initial condition is $\mathbf{u}_{0}$ with the lower script denoting the time,
$W_{in}$ is the input matrix with the values uniformly sampled in
$\left[-K_{in},K_{in}\right]$, $\mathbf{r}$ is the $N$-dimensional reservoir
state vector, $A$ is the adjacency matrix of an Erdős-Rényi network with
average degree $D$ to describe the reservoir connection between $N$ nodes, and
$\rho$ is the spectral radius of $A$. The $\tanh$ represents our activation
function for this study. The $\mathbf{\widetilde{u}}$ is the output vector and
$W_{out}$ is the output matrix. The $\alpha$ is the leak hyperparameter,
representing the time scale [35]. When we reformulate Eq. 1 in its continuous
form, it provides a more intuitive interpretation of this relationship:

$\displaystyle\frac{1}{\alpha}\dot{\mathbf{r}}=-\mathbf{r}+\tanh(A\mathbf{r}+W_{in}\mathbf{u}).$
(3)

In the training phase, only $W_{out}$ is trained to minimize the difference
between the output time series and the training data [40]. With the
regularization term, the loss function is given by

$\displaystyle L$
$\displaystyle=\sum_{t=1}^{T}||\mathbf{u}_{t}-W_{out}\mathbf{r}_{t}||^{2}+\beta||W_{out}||^{2},$
(4)

where $\beta$ is the regression hyperparameter. We then regress the matrix
$W_{out}$ by minimizing the loss function (Methods). By stacking the vector of
different time points as vector, such as
$\mathbf{U}\doteq[\mathbf{u}_{1},\dots,\mathbf{u}_{T}]$ and
$\mathbf{R}\doteq[\mathbf{r}_{1},\dots,\mathbf{r}_{T}]$ with $t=1,\dots,T$, it
can be rewritten as a compact form:

$\displaystyle
W_{out}=(\mathbf{U}\mathbf{R}^{\intercal})\cdot(\mathbf{R}\mathbf{R}^{\intercal}+\beta)^{-1},$
(5)

where $\intercal$ denotes the transpose. The problem of determining $W_{out}$
is a simple linear regression, which is less computationally expensive than
the neural network that requires the back propagation.

The framework for learning noise-induced transitions using multi-scaling
reservoir computing is summarized in FIG. 1. A reservoir acquires a time
series $\mathbf{u}$ that contains signals with both fast and slow time scales.
Given that $\alpha$ characterizes the time scale of reservoir computer, we
search for an appropriate value of $\alpha$ to capture the slow time scale
initially. After identifying an appropriate $\alpha$ value, additional
searches are conducted to find suitable values for other hyperparameters. This
process aims to improve the accuracy of the results and obtain the trained
slow-scale model ($W_{out}$). We utilize the trained slow-scale model to
separate the noise distribution from the original series. Following this, we
use the noise sampled from the separated distribution and employ the trained
slow-scale model to do rolling prediction.

In details, we leverage the conventional approach [40] to search for
appropriate hyperparameters. Given the multiple stable states in the training
set, if slow-scale model effectively captures the slowly time-scale dynamics,
the trajectories from various initial points (e.g., ten chosen points) should
converge to the corresponding stable state. Thus, we tune hyperparameters to
reach this property such that the reservoir dynamics matches the slowly time-
scale dynamics. Specifically, we first tune the hyperparameter $\alpha$
because it represents the time scale [35], and then refine the remaining
hyperparameters. The effectiveness of these adjustments is assessed by
evaluating whether the resulting ten time series converge to their
corresponding stable states. In the case of non-convergence, the
hyperparameter is adjusted in the opposite direction. When the hyperparameter
adjustments do not further improve convergence, we turn to the next
hyperparameter [40], as illustrated in FIG. 1(a).

After finding the proper hyperparameters, we utilize the trained slow-scale
model to separate the noise distribution. Within the training phase, at time
step $t$, the reservoir accepts the input $\mathbf{u}_{t}$, resulting in an
output $\mathbf{\widetilde{u}}_{t+1}$. Then, the noise at time step $t+1$ can
be computed as

$\displaystyle\mathbf{\eta}_{t+1}=\mathbf{u}_{t+1}-\mathbf{\widetilde{u}}_{t+1}.$
(6)

Following this, we obtain the noisy time series and the distribution as
depicted in FIG. 1(a). Once the trained slow-scale model and the noise
distribution have been obtained, we can implement a rolling prediction (FIG.
1(b)). In the predicting phase, at time step $s$, the reservoir accepts
$\mathbf{\widehat{u}}_{s}$, yielding the output
$\mathbf{\widetilde{u}}_{s+1}$. By adding $\mathbf{\eta}_{s+1}$, sampled from
the noise distribution, to the output $\mathbf{\widetilde{u}}_{s+1}$ as

$\displaystyle\mathbf{\widehat{u}}_{s+1}=\mathbf{\widetilde{u}}_{s+1}+\mathbf{\eta}_{s+1},$
(7)

the $\mathbf{\widehat{u}}_{s+1}$ is used as the input for the time step $s+1$.
The vector $\mathbf{\widehat{u}}$ is our prediction.

As illustrated in FIG. 1(c), to validate that the present method accurately
captures noise-induced transitions, we compare the prediction to the test
data. For white noise, as Gaussian white noise is memoryless, we quantify the
accuracy of our prediction by the statistics of noise-induced transitions,
instead of aiming at predicting a single transition, so our focus is on
learning the statistics of transition time and the number of transitions from
a set of trajectories. For colored noise, we aim to accurately forecast the
occurrence of a specific noise-induced transition.

We next proceed with two categories of examples. One category is data
generated from stochastic differential equations, including a 1D bistable
gradient system and a 2D bistable non-gradient system with white noise, as
well as a 1D bistable gradient system with colored noise. More examples are
provided in Supplementary Note: a 1D tilted bistable gradient system
(Supplement FIG. 3), a 2D bistable gradient system (Supplement FIG. 4), 2D
tilted bistable gradient (Supplement FIG. 5) and non-gradient (Supplement FIG.
6) systems, and a 2D tristable system (Supplement FIG. 7). The second category
focuses on experimental data, where we apply the present method to protein
folding data [6]. We also assess the performance of using a small part of the
dataset recorded (FIG. 5 and Supplement FIG. 2).

###  II.2 Examples

Table 1: The list of hyperparameters used for the different examples of the main text. Model | $\delta t$ |  $T_{\text{train}}$ (Time steps) |  $T_{\text{predict}}$ (Time steps) | $N$ | $K_{\text{in}}$ | $D$ | $\rho$ | $\alpha$ | $\beta$
---|---|---|---|---|---|---|---|---|---
Example 1 | $0.01$ | $10000$ | $10000$ | $800$ | $4$ | $4$ | $1.2\times 10^{-3}$ | 0.2 | $1\times 10^{-8}$
Example 2 set 1 | $0.01$ | $8000$ | $8000$ | $1000$ | $1.5$ | $3.2$ | $1.3\times 10^{-3}$ | 0.25 | $1\times 10^{-7}$
Example 2 set 2 | $0.01$ | $580$ | $300$ | $800$ | $0.996$ | $0.996$ | $0.806$ | $0.065$ | $1\times 10^{-7}$
Example 3 | $0.002$ | $20000$ | $20000$ | $1200$ | $1$ | $2.2$ | $1.7\times 10^{-3}$ | 0.3 | $1\times 10^{-7}$
Example 4 | real data | $25000$ | $100000$ | $800$ | $0.04$ | $1.8$ | $0.021$ | $0.22$ | $1\times 10^{-6}$

####  II.2.1 A bistable gradient system with white noise

Figure 2: Capturing stochastic transitions in a bistable gradient system with
white noise. (a) Schematic of noise-induced transitions in the bistable
gradient system with Gaussian white noise. (b) Generated time series from Eq.
9 ($b=5,c=0,\varepsilon=0.3,u_{1}(0)=1.5,\delta t=0.01$) with $t=30$ as ground
truth. (c) The trained slow-scale model transforms ten different start points
into ten different slowly time-scale series (color lines). Separate the noise
distribution in the training phase. (d) Result of prediction using the slow-
scale model and the noise distribution in (c), $t\in[100,130]$. (e) The number
of transitions for the test and predicted data matches. Transition refers to
the shift from one state to another. The duration of the prediction is
$10000\delta t$. (f) Histograms of transition time for the test and predicted
data. Transition time refers to the interval between the completion of the
previous transition and the completion of the current transition.

To demonstrate the present method, we first consider data generated from SDE.
The continuous-state and continuous-time Markovian stochastic dynamics can be
given as

$\displaystyle\dot{\mathbf{u}}=f(\mathbf{u})+\sigma\xi(t),$ (8)

where the vector $\dot{\mathbf{u}}$ is the time derivative, the deterministic
part of the dynamics is $f(\mathbf{u})$, and $\sigma$ corresponds to the noise
strength. The $\xi(t)$ is a $k$-dimensional Gaussian white noise with
$\langle\xi(t)\rangle=0,\langle\xi(t)\xi^{\intercal}(t^{{}^{\prime}})\rangle=\delta(t-t^{{}^{\prime}})I_{k}$,
where $I_{k}$ is the $k$-dimensional identity matrix,
$\delta(t-t^{{}^{\prime}})$ is the Dirac $\delta$ function, and
$\langle\cdots\rangle$ represents the noise average.

As a first example, we consider a 1D bistable gradient system with white noise
[9],

$\displaystyle\dot{u}_{1}=-b(-u_{1}+u_{1}^{3}+c)+\sqrt{2\varepsilon
b}\xi_{1}(t),\enspace t\geq 0,$ (9)

with a Gaussian white noise $\xi_{1}(t)$. The parameter $b$ denotes the
strength of the diffusion coefficient, $\varepsilon$ is the noise strength,
and $c$ controls the tilt of the two potential wells. The system has noise-
induced transitions between the two potential wells as illustrated in FIG.
2(a). We generated a time series with a duration of $20000\delta t$, with the
training set $t\in[0,100]$ and the predicting set $t\in[100,200]$. FIG. 2(b)
shows the first $3,000\delta t$ of the training set.

In the training phase, the tuning of hyperparameters for slow-scale model is
performed as outlined in subsection II.1, utilizing ten initial points to
search for proper hyperparameters. FIG. 2(c) illustrates the process, after
finding the proper hyperparameters listed in Table 1 (Example 1), we obtain
the trained slow-scale model. We use the generated slowly time-scale series to
separate the noise distribution, as demonstrated in FIG. 1(a). We employ the
trained slow-scale model and the separated noise distribution for rolling
prediction. The first $3000\delta t$ of the prediction is illustrated in FIG.
2(d). The prediction series demonstrates similar dynamics to the test data,
including noise-induced transitions.

In the evaluation, we generated $100$ replicates of time series from Eq. 9,
trained the model, and produced $100$ time series separately. Then, we compare
the statistics of the noise-induced transitions for these two sets of time
series, e.g., the number of transitions over $10000\delta t$ (FIG. 2(e)) and
the transition time (FIG. 2(f)). The match between the test and predicted data
demonstrates the effectiveness of our approach in capturing noise-induced
transition dynamics.

Figure 3: Predicting the accurate transition time series for a bistable
gradient system with colored noise. The system is the same as Eq. 10 [29]. (a)
The flowchart of predicting stochastic transitions with colored noise. The
process for obtaining noise $\zeta_{t}$ follows that in FIG. 1(a), and a
second reservoir takes into $\zeta_{t}$ through matrix $W_{in}^{*}$ and has
reservoir states $\mathbf{r}^{*}_{t}$ with a connection matrix $A^{*}$. The
output matrix $W_{out}^{*}$ is trained to learn noise. Using the trained slow-
scale model to predict noise-induced transitions. (b) Generated time series
($x(0)=y(0)=z(0)=1,b=1,c=0,\psi=0.08,\varepsilon=0.5,u_{1}(0)=-1.5,\delta
t=0.01$) with $8000\delta t$, where a noise-induced transition occurs in
$t\in[22,25]$ marked by the green dashed line. The noise data from $580\delta
t$ before the stochastic transition at $t=22$ is applied to predict the noisy
time series in $t\in[22,25]$. (c) The trained slow-scale model transforms ten
different start points into ten different slowly time-scale series (color
lines). Separate the noisy time series in the training phase. (d) By repeating
the process in (a) $50$ times with the same hyperparameters, $50$ predicted
$u_{1}(t)$ are obtained (fainter lines). The averaged predicted time series
(thick green) matches the test data (coral). (e) Absolute error of the
predicted $50$ time series and its mean value (thick green).

####  II.2.2 A bistable gradient system with colored noise

The prediction on a single stochastic transition may become possible when the
system has colored noise. To demonstrate that the present method is applicable
to such cases, we consider a system Eqs. 10, 11, 12 and 13 studied in [29],
where their method relies on the assumption of knowing the deterministic part
of the equation in prior. In contrast to their method, we do not assume any
prior knowledge of the deterministic part of the dynamical system and directly
learn the deterministic part and noise (FIG. 3(a)), enabling a prediction in a
model-free manner.

The system is a 1D bistable gradient system, as illustrated in FIG. 3(b):

$\displaystyle\dot{u}_{1}$
$\displaystyle=-b(-u_{1}+u_{1}^{3}+c)+\frac{\psi}{\varepsilon}y,$ (10)
$\displaystyle\dot{x}$ $\displaystyle=\frac{10}{\varepsilon^{2}}(y-x),$ (11)
$\displaystyle\dot{y}$ $\displaystyle=\frac{1}{\varepsilon^{2}}(28x-xz-y),$
(12) $\displaystyle\dot{z}$
$\displaystyle=\frac{1}{\varepsilon^{2}}(xy-\frac{8}{3}z).$ (13)

The parameter $b$ denotes the strength of the diffusion coefficient,
$\varepsilon$ corresponds to the noise strength, $\psi$ controls the influence
of the noise on the slow-scale dynamics, and $c$ controls the tilt of the two
potential wells. The noise $(x,y,z)$ is modeled by a dynamical system, the
Lorenz-63 model [41]. The system has stochastic transitions between the two
potential wells under the Lorenz noise.

To demonstrate the present method, we consider the time series for $8000\delta
t$ (FIG. 3(b)), where a stochastic transition is indicated prior to the green
dashed line. In the training phase, we obtain slow-scale model to learn the
deterministic part (FIG. 3(c)), and to separate noise. Hyperparameters are
listed in Table 1 (Example 2 set 1). In the predicting phase, accurately
forecasting the stochastic transitions requires to predict the noisy time
series. Thus, we utilize a second reservoir (FIG. 3(a)) to learn the
previously separated noise during the training phase. The hyperparameters for
the noisy time series are listed in Table 1 (Example 2 set 2). With the
deterministic component slow-scale model, we execute a rolling prediction to
predict transition.

To evaluate the accuracy of prediction, as in [29], we applied the same
hyperparameters to conduct $50$ predictions. These predictions were then used
alongside trained slow-scale model for $50$ times rolling prediction. The
average of the $50$ predictions outcomes closely approximates the actual time
series (FIG. 3(d)). Furthermore, FIG. 3(e) shows a near-zero average absolute
error between the $50$ predictions and the actual time series, indicating a
high accuracy. These results demonstrate that the present approach requires no
assumptions on knowing the deterministic part, underscoring its potential in
predicting a single stochastic transition under colored noise.

####  II.2.3 A bistable non-gradient system

Figure 4: Learning noise-induced transitions in a bistable non-gradient
system. (a) Schematic of noise-induced transitions in the 2D bistable non-
gradient system with Gaussian white noise. (b) Generated time series from Eqs.
14 and 15
($a=b=5,c=0,\varepsilon_{1}=\varepsilon_{2}=0.3,u_{1}(0)=0,u_{2}(0)=2,\delta
t=0.002$) with $t=40$ as ground truth. (c) The trained slow-scale model
transforms ten different start points into ten different slowly time-scale
series (color lines), $t\in[40,80]$. Separate the noise distribution. (d)
Result of prediction using the slow-scale model and the noise distribution in
(c). (e) The number of transitions for the 100 replicates simulated in
$t\in[40,80]$ and the 100 generated matches. The transition refers to the time
series in the $u_{1}$-direction that crosses the zero point and remains either
non-negative or non-positive for $50\delta t$. (f) Histograms of transition
time for the test and predicted data. When a transition occurs within the
system, the transition time is defined as the interval between two consecutive
zero crossings in the $u_{1}$ direction.

We next focus on investigating whether the noise-induced transitions can be
predicted by the present method for the 2D non-gradient systems. We consider a
bistable non-gradient system [38]:

$\displaystyle\dot{u}_{1}$
$\displaystyle=-b(-u_{1}+u_{1}^{3}+c)-au_{2}+\sqrt{2\varepsilon_{1}b}\xi_{1}(t),\enspace
t\geq 0,$ (14) $\displaystyle\dot{u}_{2}$
$\displaystyle=a(-u_{1}+u_{1}^{3}+c)-bu_{2}+\sqrt{2\varepsilon_{2}b}\xi_{2}(t),\enspace
t\geq 0,$ (15)

with Gaussian white noise $\xi_{1}(t)$ and $\xi_{2}(t)$. In this system, $b$
is the diffusion coefficient, $a$ represents the strength of the non-detailed
balance part, $\varepsilon_{1}$ and $\varepsilon_{2}$ are the noise strength,
and $c$ controls the tilt of the potential. The system has noise-induced
transitions between the two potential wells under the noise as illustrated in
FIG. 4(a). The presence of a non-detailed balance introduces a rotational
component to the time series, which adds difficulty to the prediction.

In the training phase, we generated a series consisting of $40000\delta t$.
The training set is $t\in[0,40]$, and the predicting set is $t\in[40,80]$.
FIG. 4(b) displays the training set. As the method in subsection II.1, the
deterministic part is plotted in FIG. 4(c). A proper set of hyperparameters is
listed in Table 1 (Example 3). We observe that the generated time series
starting from the ten initial points converge to two potential wells, where
the time series has rotational dynamics. In the predicting phase, we perform
rolling prediction within $t\in[40,80]$ (FIG. 4(d)).

In the evaluation, we predict $100$ replicates of the time series, and compare
them to $100$ replicates simulated from Eqs. 14 and 15. FIG. 4(e) presents
histograms of the number of transitions for the $100$ predicted time series
and $100$ test time series, in $t\in[40,80]$. FIG. 4(f) presents histograms of
transition time for the $100$ predicted time series and the $100$ test series,
in $t\in[40,80]$. The results demonstrate that, for the 2D bistable non-
gradient system, the present method accurately learns the dynamics and yields
precise estimations on the number of transitions and transition time.

Figure 5: Learning the stochastic transitions from the experimental data of
protein folding. The $u_{1}$ represents the end-to-end length of the protein.
We refer to transitions from around $u_{1}=15$ to around $u_{1}=30$ as upward
transitions, and vice versa as downward transitions. The right-pointing arrow:
reduction of training data. (a) Time series of the training set ($0-25000$
time steps). (b) The trained slow-scale model transforms ten different start
points into ten different slowly time-scale series (color lines). Separate
noise distribution in the training phase. (c) Result of prediction using the
trained slow-scale model and separated noise in (b) during time steps
$25000-50000$. (d) Histograms of upward transition time for the prediction and
the true data, where the length of the training set $(T_{train})$ is $25000$
time steps. Transition time refers to the interval between the completion of
the previous transition and the completion of the current transition. (e)
Histograms of downward transition time for the prediction and the true data,
with $T_{train}=25000$ time steps. (f-i) Similar histograms of upward and
downward transition time as in (d) and (e), and with different lengths of the
training sets, $T_{train}=7500$ time steps for (f) and (g), and
$T_{train}=6000$ time steps for (h) and (i). This match demonstrates that the
present method can still make predictions even when reducing the training
length to $T_{train}=7500$.

####  II.2.4 Experimental data of protein folding

To demonstrate that the present method can learn the functional role of noise
in stochastic transitions of experimental data, we apply the present method to
the protein folding data [6]. The protein can escape the free energy barrier
from one stable state to another stable state, variations in end-to-end length
correspond to different free energy, which can result in the protein being
either in an unfolded or native conformation. This process can be regarded as
a noise-induced transition. FIG. 5(a) shows the training data, where there are
transitions between two stable states.

In the training phase, with $T_{train}=25000$ time steps, we obtain the
trained slow-scale model and ten different slowly time-scale series with the
separated noise distribution (FIG. 5(b)). The proper hyperparameters are
listed in Table 1 (Example 4). In the predicting phase, we employ the trained
slow-scale model and the separated noise distribution to do rolling prediction
for $100000$ time steps. The first $25000$ time steps of prediction are
plotted in FIG. 5(c), showing the transitions between stable states and the
asymmetric dynamics.

In experiments, the available data is often limited, and it is essential to
determine the minimum amount of data required. Thus, we reduce the amount of
training data to $7500$ and $6000$ time steps separately. We generate
prediction for $100000$ time steps and then compare this prediction to the
true data to evaluate the impact of data length on prediction accuracy. The
results in FIG. 5(d-e) and FIG. 5(f-g) demonstrate that the present method can
learn the dynamics of protein folding from the data with around $7500$ time
steps. When $T_{train}=6000$ time steps, FIG. 5(h-i) show a larger error
between the predicted and true transition time. This suggests that $7500$ time
steps approximate the minimum data requirement for the present method in this
system, such that the behavior of protein folding can be effectively learnt
and simulated for more time steps. Therefore, the present approach holds
promising potential for streamlining the workload of experimentalist by
allowing to learn protein folding dynamics from a small dataset.

##  III Discussion

We have provided a general framework for noise-induced transitions that is
solely based on data. The present method exhibits considerable improvements
compared with the previous work [29]: A crucial distinction is that we can
obtain the deterministic part of dynamics without any prior knowledge about
underlying deterministic equation. We have applied the method to examples from
stochastic differential equations and the experimental data, where the method
can accurately learn noise-induced transitions and estimate transition
statistics in various systems.

The choice of hyperparameters affects the training. For example, the
hyperparameter $\alpha$ plays an important role: the larger $\alpha$
corresponds to the time series with fast time scale, while the smaller
$\alpha$ leads to the slow time scale [35]. We utilize this characteristic to
search for $\alpha$ to match with the slow dynamics and separate noise. For
the example with colored noise (subsubsection II.2.2), the $\alpha$ for noise
is smaller than that for the deterministic part (Table 1). We find that using
a smaller $\alpha$ to learn the noisy time series can better predict
transitions than using a larger $\alpha$. As the fast and chaotic noisy time
series is often difficult to predict, using a small $\alpha$ leads to smoother
noisy time series. While this selection may sacrifice some details of noise,
it enhances the ability to capture the major trend of noise and thus improves
transition prediction. In asymmetric systems, such as tilted bistable gradient
systems, the two distinct potential wells of the system exhibit different time
scales. Consequently, we need to employ two different sets of hyperparameters
(Supplement FIGs. 3, 5, and 6). If a time series is generated from an
asymmetric system, our multi-scaling reservoir computing approach can identify
the tilted dynamics and accurately simulate two types of transitions, even
without any prior information about its tilt and noise strength.

The effectiveness of learning slowly time-scale series can also be influenced
by other hyperparameters [40]. In our protocol, we utilize ten initial points
to search for appropriate hyperparameters. Given the inherent complexity of
neural network and random matrices in reservoir computer, developing a
universal and systematic strategy for selecting predefined hyperparameters is
a significant challenge [42, 43]. As potential ways of improvements, the
application of Bayesian optimization [44] and simulated annealing [45] can be
used to facilitate the search for hyperparameters. Additionally, the
convergence speed of the ten slowly time-scale series may have some
discrepancies compared with the actual dynamics. As a result, the noise
distribution separated by the slow-scale model during the training phase may
exhibit either lower or higher intensity compared with the actual noise
distribution. In this case, we can employ a noise factor to magnify or reduce
separated noise, aiming to predict transitions more accurately. To verify this
idea, we amplify the sampled noise by a factor of $1.1$ based on subsubsection
II.2.1, which improves the accuracy of the predictions (Supplement FIG. 1).

For experimental data, we have shown the feasibility of predicting the
experimental time series from a small dataset, as exemplified in the analysis
of the protein folding data [6]. We can capture the transitions and the tilted
potential dynamics. In this example, over short timescales, the protein
samples a local equilibrium involving the native and unfolded conformations.
However, if the measurement time is significantly extended, previously
inaccessible regions separated by high-energy barriers can be explored.
Consequently, in order to capture a wider variety of protein folding
transitions, it may be necessary to use longer training set, which would
result in higher computational costs. Additionally, we observed tilted
dynamics in the time series of protein folding. Even so, we can learn both the
upward and downward transitions by using only one set of hyperparameters. This
suggests that in this system the upward and downward transition time scales
might not differ significantly. The noise strength drives the system to cross
one energy barrier more frequently while encountering difficulty in crossing
another. Therefore, when dealing with a tilted dynamics time series, we can
utilize two distinct sets of hyperparameters to learn the different time
scales.

The recent Python library for machine learning dynamical models from time
series, Deeptime [46], has not directly handled the stochastic transitions.
There are three categories of methods in Deeptime related to our work. The
first category is the deep dimension reduction, such as deep Koopman networks.
These methods can find stable states from time series and identify proper time
scales to separate slow and fast dynamics [47]. However, they have not
achieved to predict stochastic transitions between stable states. The second
category method is the sparse identification of nonlinear dynamics (SINDy).
Although it can identify nonlinear dynamics from noisy data [16], it mainly
regards noise as a disturbance and is not designed for learning noise-induced
phenomena. The third category comprises Markov state models and hidden Markov
models. These methods are applicable to stochastic processes with discrete
states, whereas the present study focuses on stochastic processes with
continuous states.

There are more potential applications of the present method in various
domains. In physics, such as in Nagumo’s tunnel diode neuron model, the
inclusion of a nonlinear resistor [2] can lead to noise-induced transitions.
We can apply the present method to learn these transitions between different
voltage and current stable states. The present approach may also be extended
to analyze transition of trajectories between different dynamical phases of
spins [48, 49], and to study the switches of different meander states in
climate change [50]. Furthermore, we can further extend the approach for
examples involving a wider variety of noise types, such as band-limited
thermal noise, power law noise, shot noise, and impulsive noise, where the
conditional generative adversarial network [51] may be employed to model these
noise.

##  IV Methods

In this section, we reformulate the loss function to derive the expression for
the output matrix $W_{out}$ and discuss the hyperparameters in the present
method. The loss function is given as Eq. 4. In detail, we should write the
loss function as a sum from all the parameters to do linear regression. Then,
the regression becomes simply a sum of vectors:

$\displaystyle L$
$\displaystyle=\sum_{t=1}^{T}[||\mathbf{u}_{t}-W_{out}\mathbf{r}_{t}||^{2}+\beta||W_{out}||^{2}]$
$\displaystyle=\sum_{t=1}^{T}[(\mathbf{u}_{t}-W_{out}\mathbf{r}_{t})^{\intercal}(\mathbf{u}_{t}-W_{out}\mathbf{r}_{t})+\beta||W_{out}||^{2}]$
$\displaystyle=\sum_{t=1}^{T}[(\mathbf{u}_{t})^{\intercal}\mathbf{u}_{t}-(W_{out}\mathbf{r}_{t})^{\intercal}\mathbf{u}_{t}-(\mathbf{u}_{t})^{\intercal}W_{out}\mathbf{r}_{t}+(W_{out}\mathbf{r}_{t})^{\intercal}W_{out}\mathbf{r}_{t}+\beta||W_{out}||^{2}].$
(16)

As the loss function is convex (prove that the zero gradient is indeed the
local minimum, one needs to differentiate once more to obtain the Hessian
matrix and show that it is positive definite. This is provided by the Gauss-
Markov theorem), the optimum solution lies at the zero gradient by

$\displaystyle\partial_{W_{out}}L$
$\displaystyle=\sum_{t=1}^{T}[-2(\mathbf{r}_{t})^{\intercal}\mathbf{u}_{t}+2(\mathbf{r}_{t})^{\intercal}W_{out}\mathbf{r}_{t}+2\beta
W_{out}]=0,$ (17)

which leads to the regression:

$\displaystyle
W_{out}=\sum_{t=1}^{T}\left[(\mathbf{u}_{t})\cdot(\mathbf{r}_{t})^{\intercal}\right]\cdot[(\mathbf{r}_{t})\cdot(\mathbf{r}_{t})^{\intercal}+\beta]^{-1},$
(18)

where we have neglected the notation of identity matrix and identify vector.
By stacking the vector of different time points as vector, it can be rewritten
as a compact form Eq. 5.

There are six hyperparameters in the present method. The variable $N$
represents the number of reservoir nodes, which determines reservoir size. In
most instances, performance improves with larger reservoir [52]. However,
using large reservoir might lead to overfitting, requiring the application of
suitable regularization techniques [40]. The hyperparameter $K_{in}$
represents the scaling factor for the input matrix $W_{in}$. The average
degree of the reservoir connection network is denoted by $D$, and we choose
the connection matrix $A$ to be sparse [20]. This approach stems from the
intuition that decoupling the state variables can result in a richer encoding
of the input signal [52]. The spectral radius of the reservoir connection
network, denoted as $\rho$, represents a critical characteristic of the
dynamics of the reservoir state. Notably, it affects both the nonlinearity of
the reservoir and its capacity to encode past inputs in its state [52, 36].
The $\alpha$ represents the leak parameter, it represents the time scale [35].
The $\beta$ is the regularization parameter [40].

Data availability: The authors declare that the data supporting this study are
available within the paper.

Code availability: A pytorch implementation of the present algorithm will be
publicly available upon the acceptance of the manuscript.

## Acknowledgments

We acknowledge Xingang Wang for helpful discussions. We thank Rafael Tapia-
Rojo for sharing the data of protein folding. This work is supported by
Project 12322501, 12105014 of National Natural Science Foundation of China.
Y.T. acknowledges the start-up research funding (28705-310432101) by Beijing
Normal University. The HPC is supported by Interdisciplinary Intelligence
SuperComputer Center of Beijing Normal University, Zhuhai.

## Author contributions

Y.T. had the original idea for this work, Z.Q.L. and Z.F.L. performed the
study, and all authors contributed to the preparation of the manuscript.

## References

  * Horsthemke and Lefever [2006] W. Horsthemke and R. Lefever,  _Noise-Induced Transitions: Theory and Applications in Physics, Chemistry, and Biology_ , 2nd ed. (Springer-Verlag, Berlin, 2006).
  * Semenov  _et al._ [2016] V. V. Semenov, A. B. Neiman, T. E. Vadivasova, and V. S. Anishchenko, Noise-induced transitions in a double-well oscillator with nonlinear dissipation, Phys. Rev. E 93, 052210 (2016).
  * Assaf  _et al._ [2011] M. Assaf, E. Roberts, and Z. Luthey-Schulten, Determining the stability of genetic switches: Explicitly accounting for mrna noise, Phys. Rev. Lett. 106, 248102 (2011).
  * Jafarpour  _et al._ [2015] F. Jafarpour, T. Biancalani, and N. Goldenfeld, Noise-induced mechanism for biological homochirality of early life self-replicators, Phys. Rev. Lett. 115, 158101 (2015).
  * Qian [2002] H. Qian, From discrete protein kinetics to continuous brownian dynamics: A new perspective, Protein Sci. 11, 1 (2002).
  * Tapia-Rojo  _et al._ [2023] R. Tapia-Rojo, M. Mora, S. Board, J. Walker, R. Boujemaa-Paterski, O. Medalia, and S. Garcia-Manyes, Enhanced statistical sampling reveals microscopic complexity in the talin mechanosensor folding energy landscape, Nat. Phys. 19, 52 (2023).
  * Hänggi  _et al._ [1990] P. Hänggi, P. Talkner, and M. Borkovec, Reaction-rate theory: fifty years after kramers, Rev. Mod. Phys. 62, 251 (1990).
  * Tang  _et al._ [2023] Y. Tang, J. Weng, and P. Zhang, Neural-network solutions to stochastic reaction networks, Nat. Mach. Intell. 5, 376 (2023).
  * Forgoston and Moore [2018] E. Forgoston and R. O. Moore, A primer on noise-induced transitions in applied dynamical systems, SIAM Rev. 60, 969 (2018).
  * Hartmann  _et al._ [2013] C. Hartmann, R. Banisch, M. Sarich, T. Badowski, and C. Schütte, Characterization of rare events in molecular dynamics, Entropy 16, 350 (2013).
  * Tanaka  _et al._ [2019] G. Tanaka, T. Yamane, J. B. Héroux, R. Nakane, N. Kanazawa, S. Takeda, H. Numata, D. Nakano, and A. Hirose, Recent advances in physical reservoir computing: A review, Neural Netw. 115, 100 (2019).
  * Xiong and Zhao [2019] Y. Xiong and H. Zhao, Chaotic time series prediction based on long short-term memory neural networks, Sci. China Phy. Mech. Astron. 49, 120501 (2019).
  * Karniadakis  _et al._ [2021] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Physics-informed machine learning, Nat. Rev. Phys. 3, 422 (2021).
  * Zhao [2021] H. Zhao, Inferring the dynamics of “black-box” systems using a learning machine, Sci. China Phy. Mech. Astron. 64, 270511 (2021).
  * Li  _et al._ [2023a] Q. Li, T. Wang, V. Roychowdhury, and M. K. Jawed, Meta-learning generalizable dynamics from trajectories, arXiv preprint arXiv:2301.00957 (2023a).
  * Kaheman  _et al._ [2022] K. Kaheman, S. L. Brunton, and J. N. Kutz, Automatic differentiation to simultaneously identify nonlinear dynamics and extract noise probability distributions from data, Mach. Learn.: Sci. Technol. 3, 015031 (2022).
  * Raissi  _et al._ [2019] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys. 378, 686 (2019).
  * Cuomo  _et al._ [2022] S. Cuomo, V. S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, and F. Piccialli, Scientific machine learning through physics–informed neural networks: Where we are and what’s next, J. Sci. Comput. 92, 88 (2022).
  * Lusch  _et al._ [2018] B. Lusch, J. N. Kutz, and S. L. Brunton, Deep learning for universal linear embeddings of nonlinear dynamics, Nat. Commun. 9, 4950 (2018).
  * Jaeger [2001] H. Jaeger, The “echo state” approach to analysing and training recurrent neural networks-with an erratum note, Bonn, Germany: German National Research Center for Information Technology GMD Technical Report 148, 13 (2001).
  * 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 comput. 14, 2531 (2002).
  * Pathak  _et al._ [2018] J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach, Phys. Rev. Lett. 120, 024102 (2018).
  * Kim and Bassett [2023] J. Z. Kim and D. S. Bassett, A neural machine code and programming framework for the reservoir computer, Nat. Mach. Intell. 5, 622 (2023).
  * Jaeger and Haas [2004] H. Jaeger and H. Haas, Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication, Science 304, 78 (2004).
  * Zimmermann and Parlitz [2018] R. S. Zimmermann and U. Parlitz, Observing spatio-temporal dynamics of excitable media using reservoir computing, Chaos 28, 043118 (2018).
  * Fan  _et al._ [2020] H. Fan, J. Jiang, C. Zhang, X. Wang, and Y.-C. Lai, Long-term prediction of chaotic systems with machine learning, Phys. Rev. Res. 2, 012080 (2020).
  * Kim  _et al._ [2021] J. Z. Kim, Z. Lu, E. Nozari, G. J. Pappas, and D. S. Bassett, Teaching recurrent neural networks to infer global temporal structure from local examples, Nat. Mach. Intell. 3, 316 (2021).
  * Zhai  _et al._ [2023] Z.-M. Zhai, L.-W. Kong, and Y.-C. Lai, Emergence of a resonance in machine learning, Phys. Rev. Res. 5, 033127 (2023).
  * Lim  _et al._ [2020] S. H. Lim, L. Theo Giorgini, W. Moon, and J. S. Wettlaufer, Predicting critical transitions in multiscale dynamical systems using reservoir computing, Chaos 30, 123126 (2020).
  * Carroll [2018] T. L. Carroll, Using reservoir computers to distinguish chaotic signals, Phys. Rev. E 98, 052209 (2018).
  * Nakai and Saiki [2018] K. Nakai and Y. Saiki, Machine-learning inference of fluid variables from data using reservoir computing, Phys. Rev. E 98, 023111 (2018).
  * Weng  _et al._ [2019] T. Weng, H. Yang, C. Gu, J. Zhang, and M. Small, Synchronization of chaotic systems and their machine-learning models, Phys. Rev. E 99, 042203 (2019).
  * Grigoryeva and Ortega [2018] L. Grigoryeva and J.-P. Ortega, Echo state networks are universal, Neural Netw. 108, 495 (2018).
  * Zhang  _et al._ [2021] H. Zhang, H. Fan, L. Wang, and X. Wang, Learning hamiltonian dynamics with reservoir computing, Phys. Rev. E 104, 024205 (2021).
  * Tanaka  _et al._ [2022] G. Tanaka, T. Matsumori, H. Yoshida, and K. Aihara, Reservoir computing with diverse timescales for prediction of multiscale dynamics, Phys. Rev. Res. 4, L032014 (2022).
  * Du  _et al._ [2023] Y. Du, Q. Li, H. Fan, M. Zhan, J. Xiao, and X. Wang, Inferring attracting basins of power system with machine learning, arXiv preprint arXiv:2305.14374 (2023).
  * Fang  _et al._ [2023] C. Fang, Y. Lu, T. Gao, and J. Duan, Reservoir computing with error correction: Long-term behaviors of stochastic dynamical systems, arXiv preprint arXiv:2305.00669 (2023).
  * Tang  _et al._ [2018] Y. Tang, S. Xu, and P. Ao, Escape rate for nonequilibrium processes dominated by strong non-detailed balance force, J. Chem. Phys. 148, 064102 (2018).
  * Belkacemi  _et al._ [2021] Z. Belkacemi, P. Gkeka, T. Lelièvre, and G. Stoltz, Chasing collective variables using autoencoders and biased trajectories, J. Chem. Theory Comput. 18, 59 (2021).
  * Lukoševičius [2012] M. Lukoševičius, A practical guide to applying echo state networks, in  _Neural Networks: Tricks of the Trade: Second Edition_ (Springer, 2012) pp. 659–686.
  * Lorenz [1963] E. N. Lorenz, Deterministic nonperiodic flow, J. Atmos. Sci. 20, 130 (1963).
  * Jiang and Lai [2019] J. Jiang and Y.-C. Lai, Model-free prediction of spatiotemporal dynamical systems with recurrent neural networks: Role of network spectral radius, Phys. Rev. Research 1, 033056 (2019).
  * Gauthier  _et al._ [2021] D. J. Gauthier, E. Bollt, A. Griffith, and W. A. Barbosa, Next generation reservoir computing, Nat. Commun. 12, 5564 (2021).
  * Yperman and Becker [2016] J. Yperman and T. Becker, Bayesian optimization of hyper-parameters in reservoir computing, arXiv preprint arXiv:1611.05193 (2016).
  * Ren and Ma [2022] B. Ren and H. Ma, Global optimization of hyper-parameters in reservoir computing, Electron. Res. Arch. 30, 2719 (2022).
  * Hoffmann  _et al._ [2021] M. Hoffmann, M. Scherer, T. Hempel, A. Mardt, B. de Silva, B. E. Husic, S. Klus, H. Wu, N. Kutz, S. L. Brunton, _et al._ , Deeptime: a Python library for machine learning dynamical models from time series data, Mach. Learn.: Sci. Technol. 3 015009 3, 015009 (2021).
  * Li  _et al._ [2023b] R. Li, H. Wang, and Y. Li, Learning slow and fast system dynamics via automatic separation of time scales, in  _Proc. ACM SIGKDD Int. Conf. Knowl. Discov. Data Min._, KDD ’23 (Association for Computing Machinery, New York, NY, USA, 2023) p. 4380–4390.
  * Casert  _et al._ [2022] C. Casert, I. Tamblyn, and S. Whitelam, Learning stochastic dynamics and predicting emergent behavior using transformers, arXiv preprint arXiv:2202.08708 (2022).
  * Tang  _et al._ [2022] Y. Tang, J. Liu, J. Zhang, and P. Zhang, Solving nonequilibrium statistical mechanics by evolving autoregressive neural networks, arXiv preprint arXiv:2208.08266 (2022).
  * Vanden-Eijnden and Weare [2013] E. Vanden-Eijnden and J. Weare, Data assimilation in the low noise regime with application to the kuroshio, Mon. Weather Rev. 141, 1822 (2013).
  * Wunderlich and Sklar [2022] A. Wunderlich and J. Sklar, Data-driven modeling of noise time series with convolutional generative adversarial networks, Mach. Learn.: Sci. Technol. 4, 035023 (2022).
  * Verzelli [2022] P. Verzelli,  _Learning dynamical systems using dynamical systems: the reservoir computing approach_ , Ph.D. thesis, Università della Svizzera italiana (2022).

##  V Supplementary Information

##  VI Supplementary note

In the Supplementary Note, we first decrease the length of the training set to
observe its effect on training. We also apply the present method to more
systems, including a 1D tilted gradient system, a 2D gradient system, 2D
tilted gradient and non-gradient systems, and a 2D tristable system. Figures
of the results are in Supplementary Figures.

###  VI.1 Assessing the impact of training set length on performance

To assess the performance of the present method with the limited training
data, we shorten the length of the training set based on Example 1 of the main
text. Additionally, we keep the length of the predicting set same as the
training set. Supplement FIG. 2(a-b) display the same results in FIG. 2(e-f)
of the main text, over the same hyperparameters and length ($10000\delta t$).
Next, we reduce the length of the training and predicting sets by half
($5000\delta t$), as illustrated in Supplement FIG. 2(c-d). Without any
significant changes in the results, we continue to reduce the length by half,
this means that the training and predicting sets are only one-fourth of their
original size ($2500\delta t$). However, the presented results (Supplement
FIG. 2(e-f)) exhibit larger errors compared with the results in Supplement
FIG. 2(a-d). Furthermore, we proceed to halve the length ($1250\delta t$) of
both the training and predicting sets once again. The results depicted in
Supplement FIG. 2(g-h) demonstrate significant errors. This example shows the
training set length impact on our training.

###  VI.2 More examples

####  VI.2.1 A one-dimensional tilted bistable gradient system

To demonstrate whether the noise-induced transitions can be predicted by the
present method for systems with the tilted double-well potential, we apply the
present method to the system depicted in Supplement FIG. 3. This system is the
same as Example 1 of the main text, and the parameter $c=0.25$ makes the
system to be tilted. The distinct time scales of upward and downward
transitions necessitate two sets of hyperparameters for effective learning.
The “upward” refers to the transitions from down state to up state, and the
“downward” refers to the transitions from up state to down state.

Supplement FIG. 3(a) shows that we focus on the upward transitions, and
Supplement FIG. 3(b-f) are the results. As the same method in Framework, we
find proper hyperparameters set $1$ (for upward transitions) listed in
Supplementary Table 1 to obtain the trained slow-scale model (for upward
transitions). The ten different slowly time-scale series and separated noise
distribution as illustrated in Supplement FIG. 3(b). In the predicting phase,
we perform rolling prediction with the trained slow-scale model (for upward
transitions) and the noise distribution. Then, we evaluate the transitions
within the test and predicted data.

We note a significant difference in the downward transition time between test
and predicted data (Supplement FIG. 3(c)), and with differences in the number
of downward transitions per $10000\delta t$ (Supplement FIG. 3(d)). Meanwhile,
our evaluation of the upward transition time (Supplement FIG. 3(e)) and the
number of upward transitions in $10000\delta t$ (Supplement FIG. 3(f)) match.
It indicates the effectiveness of hyperparameters set $1$ for upward
transitions.

Using the same training set, we switch to downward transitions, with results
shown in Supplement FIG. 3(g-l). As we shift focus, we observe the results
similar to the upward case. With the proper hyperparameters set $2$ (in
Supplementary Table 1) for downward transitions, we determine the
corresponding slow-scale model (for downward transitions). The results,
including downward transition time (Supplement FIG. 3(i)) and the number of
downward transitions in $10000\delta t$ (Supplement FIG. 3(j)) match.
Conversely, for upward transitions, a substantial error is observable in
Supplement FIG. 3(k), while Supplement FIG. 3(l) shows the result with an
insignificant error. It displays the effectiveness of hyperparameters set $2$
for downward transitions.

####  VI.2.2 A two-dimensional bistable gradient system

Before the work in Example 3 of the main text, we design a pre-experiment, to
ensure the present method can apply in a 2D bistable system. The system has
noise-induced transitions between the two potential wells under the noise as
illustrated in Supplement FIG. 4(a), and the scatter plot of the training set
(consisting of $20000$ data points) is depicted in Supplement FIG. 4(b).
According to the method in Framework, we generate ten different slowly time-
scale series in the training phase (Supplement FIG. 4(c)) and prediction
(consisting of 20000 data points) in the predicting phase (Supplement FIG.
4(d)). In Supplement FIG. 4(e-f), we compare the number of transitions and the
transition time (over $20000\delta t$) of generated data and predictions. The
proper hyperparameters for this system are listed in Supplementary Table 2.

####  VI.2.3 A two-dimensional tilted bistable gradient system

To demonstrate the present method can predict noise-induced transitions in a
2D tilted system, we first apply it to a 2D bistable gradient system with
white noise as illustrated in Supplement FIG. 5.

Similar to the 1D tilted double-well potential system, the upward and the
downward transitions necessitate two sets of hyperparameters. The “upward”
refers to the transitions from down state to up state, and the “downward”
refers the transitions from up state to down state. As shown in Supplement
FIG. 5, the same graphical scheme used in Supplement FIG. 3 is applied to
represent the results of the two-dimensional tilted bistable gradient system
with Gaussian white noise. Supplement FIG. 5(a-f) show the results for upward
transitions, and the hyperparameters set $1$ (for upward transitions) are
listed in Supplementary Table 3. Supplement FIG. 5(g-l) show the results for
downward transitions, and the hyperparameters set $2$ (for downward
transitions) are listed in Supplementary Table 3.

####  VI.2.4 A two-dimensional tilted bistable non-gradient system

To demonstrate the present method can predict noise-induced transitions in a
both tilted and non-gradient system, we add a rotation element to make this
system without detailed balance as illustrated in Supplement FIG. 6 to apply
our approach.

As shown in Supplement FIG. 6, the same graphical scheme used in Supplement
FIG. 5 is applied to represent the results of a two-dimensional non-gradient
system with double-well potential and Gaussian white noise. The “upward”
refers to the transitions from down state to up state, and the “downward”
refers the transitions from up state to down state. Supplement FIG. 6(a-f)
show the results for upward transitions, and the hyperparameters set $3$ (for
upward transitions) are listed in Supplementary Table 3. Supplement FIG.
6(g-l) show the results for downward transitions, and the hyperparameters set
$4$ (for downward transitions) are listed in Supplementary Table 3.

####  VI.2.5 A two-dimensional tristable system

To assess the performance of the present method in multi-stable systems, we
study a 2D tristable system [39] as illustrated in Supplement FIG. 7(a):

$\displaystyle\begin{split}\dot{u}_{1}&=10u_{1}e^{-u_{1}^{2}}\left[e^{-(u_{2}-1/3)^{2}}-e^{-(u_{2}-5/3)^{2}}\right]-8e^{-u_{2}^{2}}\left[(u_{1}-1)e^{-(u_{1}-1)^{2}}+(u_{1}+1)e^{-(u_{1}+1)^{2}}\right]\\\
&\ -0.8u_{1}^{3}+\sqrt{\frac{2}{\gamma}}\xi_{1}(t),\enspace t\geq
0,\end{split}$ (S1)
$\displaystyle\begin{split}\dot{u}_{2}&=-8u_{2}e^{-u_{2}^{2}}\left[e^{-(u_{1}-1)^{2}}+e^{-(u_{1}+1)^{2}}\right]+10e^{-u_{1}^{2}}\left[(u_{2}-\tfrac{1}{3})e^{-(u_{2}-1/3)^{2}}-(u_{2}-\tfrac{5}{3})e^{-(u_{2}-5/3)^{2}}\right]\\\
&\quad-0.8(u_{2}-\tfrac{1}{3})^{3}+\sqrt{\frac{2}{\gamma}}\xi_{2}(t),\enspace
t\geq 0.\end{split}$ (S2)

The system has noise-induced transitions under the Gaussian white noise
($\xi_{1}(t),\xi_{2}(t)$) as in Eq. (8) of the main text, and $\gamma$
corresponds to the noise strength.

To investigate transitions among three stable states, we specifically consider
the $u_{1}$ direction due to the overlapping of two wells in the $u_{2}$
direction. Therefore, our study focuses on a one-dimensional time series in
the $u_{1}$ direction. We generated a time series from Eqs. S1 and S2, in
$t\in[0,100]$. The time series is plotted for the first $30000\delta t$ in
Supplement FIG. 7(b).

In the training phase, with the training set $t\in[0,50]$, we obtain the
appropriate hyperparameters listed in Supplementary Table 4. The time series
from ten distinct initial points converge to the stable states
$(u_{1}=-1,u_{1}=0,u_{1}=1)$ as illustrated in Supplement FIG. 7(c).
Subsequently, we perform rolling prediction for $t\in[50,80]$ using the
separated noise and the trained slow-scale model. The result demonstrates that
the method learns the multi-state dynamics of the generated time series,
including the transitions among these stable states.

##  VII Supplement figures

Supplement FIG. 1: Results of increasing the intensity of sampled noise in
rolling prediction based on Example 1 of the main text. PDF: probability
density function. (a) The number of transitions of test and predicted data
over $10000\delta t$. (b) Histograms of transition time of test and predicted
data within $10000\delta t$. (c-d) Same as (a-b), and sampled noise is
amplified by a factor of 1.1 in rolling prediction.

Supplement FIG. 2: The influence of decreasing the length of training and
predicting sets based on Example 1 of the main text. PDF: probability density
function. (a) The number of transitions over the 10000 test and predicted data
points. (b) Histograms of transition time of the 10000 test and predicted data
points. (c-d) Results over the 5000 test and predicted data points. (e-f)
Results over the 2500 test and predicted data points. (g-h) Results over the
1250 test and predicted data points.

Supplement FIG. 3:  Capturing stochastic transitions in a 1D tilted bistable
gradient system with white noise. The system is the same as Example 1 of the
main text, and the parameter $c=0.25$, making the potential tilted. Generated
time series from Eq. (9) of the main text ($b=5,c=0.25,\varepsilon=0.3,\delta
t=0.01$), spanning a duration of $20000\delta t$, with the training set
$t\in[0,100]$ and the predicting set $t\in[100,200]$. (a) Schematic of noise-
induced upward transitions in the 1D tilted bistable gradient system with
Gaussian white noise. For figures (b-f), we focus on the upward transitions.
(b) The trained slow-scale model transforms ten different start points into
ten different slowly time-scale series (color lines). Separate the noise
distribution. (c) Histograms of downward transition time for the test and
predicted data. (d) The number of downward transitions for the test and
predicted data. The duration of the predicting set is $10000\delta t$ (10000
data points). (e) Histograms of upward transition time for the test and
predicted data. (f) The number of upward transitions for the test and
predicted data. The duration of the predicting set is $10000\delta t$ (10000
data points). (g-l) Same as (a-f), and focusing on the downward transitions.

Supplement FIG. 4:  Capturing stochastic transitions in a 2D bistable gradient
system with white noise. The system is the same as Example 3 of the main text,
and the parameter $a=0$. (a) Schematic of noise-induced transitions in the 2D
bistable gradient system with Gaussian white noise. (b) Generated time series
from Eqs. (14) and (15) of the main text
($a=0,b=5,c=0,\varepsilon_{1}=\varepsilon_{2}=0.3,\delta t=0.002$) with $t=80$
as ground truth. (c) The trained slow-scale model transforms ten different
start points into ten different slowly time-scale series (color lines),
$t\in[40,80]$. Separate the noise distribution in the training phase. (d)
Result of prediction in $t\in[40,80]$ using the slow-scale model and noise in
(c). (e) The number of transitions for the 100 replicates simulated in
$t\in[40,80]$ and the 100 generated matches. The transition refers to the time
series in the $u_{1}$-direction that crosses the zero point and remains either
non-negative or non-positive for $50\delta t$. (f) Histograms of transition
time for the test and predicted data. When a transition occurs within the
system, the transition time is defined as the interval between two consecutive
zero crossings in the $u_{1}$ direction.

Supplement FIG. 5: Capturing stochastic transitions in a 2D tilted bistable
gradient system with white noise. The system is the same as Example 3 of the
main text, and the parameter $c=0.25$ making the potential tilted, $a=0$,
$\delta t=0.002$. (a) Schematic of noise-induced upward transitions in the 2D
tilted bistable gradient system. For figures (b-f), we focus on the upward
transitions. (b) The trained slow-scale model transforms ten different start
points into ten different slowly time-scale series (color lines),
$t\in[40,80]$. Separate the noise distribution in the training phase. (c)
Histograms of downward transition time for the test and predicted data. (d)
The number of downward transitions for the test and predicted data. The
duration of the predicting set is $20000\delta t$ (20000 data points). (e)
Histograms of upward transition time for the test and predicted data. (f) The
number of upward transitions for the test and predicted data. (g-l) Same as
(a-f), and focusing on the downward transitions.

Supplement FIG. 6:  Learning noise-induced transitions in a 2D tilted bistable
non-gradient system. The system is the same as Example 3 of the main text, and
the parameter $c=0.25$ making the potential tilted, $a=2$ represents the
strength of the non-detailed balance part, $\delta t=0.001$. (a) Schematic of
noise-induced upward transitions in the 2D tilted bistable non-gradient system
with Gaussian white noise. For figures (b-f), we focus on the upward
transitions. (b) The trained slow-scale model transforms ten different start
points into ten different slowly time-scale series (color lines),
$t\in[25,50]$. Separate the noise distribution in the training phase. (c)
Histograms of downward transition time for the test and predicted data. (d)
The number of downward transitions for the test and predicted data. The
duration of the predicting set is $25000\delta t$ (25000 data points). (e)
Histograms of upward transition time for the test and predicted data. (f) The
number of upward transitions for the test and predicted data. (g-l) Results
focusing on the downward transitions.

Supplement FIG. 7: Capturing stochastic transitions in a 2D tristable gradient
system with white noise. (a) Schematic of the 2D tristable system. (b)
Generated time series from the Eqs. S1 and S2 ($\gamma=1,\delta t=0.001$) of
[39] with $t=30$ as ground truth. In the scenario where $\gamma=1$, the three
potential wells exhibit comparable depths. (c) The trained slow-scale model
transforms ten different start points into ten different slowly time-scale
series (color lines), only the $u_{1}$ direction is depicted. (d) Result of
prediction using the slow-scale model and separated noise in $t\in[50,80]$.

##  VIII Supplementary tables

Supplementary Table 1: Hyperparameters for subsubsection VI.2.1.  Set | Transitions | $\delta t$ |  $T_{\text{train}}$ (Time steps) |  $T_{\text{predict}}$ (Time steps) | $N$ | $K_{\text{in}}$ | $D$ | $\rho$ | $\alpha$ | $\beta$
---|---|---|---|---|---|---|---|---|---|---
1 | Upward | 0.01 | 10000 | 10000 | 800 | 4.2 | 4 | $1.1\times 10^{-3}$ | 0.18 | $1\times 10^{-8}$
2 | Downward | 0.01 | 10000 | 10000 | 800 | 3.8 | 4 | $1.2\times 10^{-3}$ | 0.25 | $1\times 10^{-8}$
Supplementary Table 2: Hyperparameters for subsubsection VI.2.2. $\delta t$ |  $T_{\text{train}}$ (Time steps) |  $T_{\text{predict}}$ (Time steps) | $N$ | $K_{\text{in}}$ | $D$ | $\rho$ | $\alpha$ | $\beta$
---|---|---|---|---|---|---|---|---
0.002 | 20000 | 20000 | 1200 | 1.5 | 2 | $1.6\times 10^{-3}$ | 0.36 | $1\times 10^{-7}$
Supplementary Table 3: Hyperparameters for subsubsection VI.2.3 and subsubsection VI.2.4. Set | Transitions | $\delta t$ |  $T_{\text{train}}$ (Time steps) |  $T_{\text{predict}}$ (Time steps) | $N$ | $K_{\text{in}}$ | $D$ | $\rho$ | $\alpha$ | $\beta$
---|---|---|---|---|---|---|---|---|---|---
1 | Upward | 0.002 | 20000 | 20000 | 1200 | 1.6 | 3.6 | $1.4\times 10^{-3}$ | 0.58 | $1\times 10^{-7}$
2 | Downward | 0.002 | 20000 | 20000 | 1200 | 1.5 | 3.5 | $1.1\times 10^{-3}$ | 0.64 | $1\times 10^{-7}$
3 | Upward | 0.001 | 25000 | 25000 | 1200 | 1.1 | 3 | $9.1\times 10^{-4}$ | 0.52 | $1\times 10^{-7}$
4 | Downward | 0.001 | 25000 | 25000 | 1200 | 1.3 | 3 | $8.9\times 10^{-4}$ | 0.6 | $1\times 10^{-7}$
Supplementary Table 4: Hyperparameters for subsubsection VI.2.5. $\delta t$ |  $T_{\text{train}}$ (Time steps) |  $T_{\text{predict}}$ (Time steps) | $N$ | $K_{\text{in}}$ | $D$ | $\rho$ | $\alpha$ | $\beta$
---|---|---|---|---|---|---|---|---
0.001 | 50000 | 50000 | 1200 | 1 | 2.8 | $1.7\times 10^{-3}$ | 0.45 | $1\times 10^{-6}$