Generalized synchronization in the presence of dynamical noise and its detection via recurrent neural networks

Amigó, José M.; Dale, Roberto; King, Juan C.; Lehnertz, Klaus

doi:10.1063/5.0235802

Given two unidirectionally coupled nonlinear systems, we speak of generalized synchronization when the responder “follows” the driver. Mathematically, this situation is implemented by a map from the driver state space to the responder state space termed the synchronization map. In nonlinear times series analysis, the framework of the present work, the existence of the synchronization map amounts to the invertibility of the so-called cross map, which is a continuous map that exists in the reconstructed state spaces for typical time-delay embeddings. The cross map plays a central role in some techniques to detect functional dependencies between time series. In this paper, we study the changes in the “noiseless scenario” just described when noise is present in the driver, a more realistic situation that we call the “noisy scenario.” Noise will be modeled using a family of driving dynamics indexed by a finite number of parameters, which is sufficiently general for practical purposes. In this approach, it turns out that the cross and synchronization maps can be extended to the noisy scenario as families of maps that depend on the noise parameters, and only for “generic” driver states in the case of the cross map. To reveal generalized synchronization in both the noiseless and noisy scenarios, we check the existence of synchronization maps of higher periods (introduced in this paper) using recurrent neural networks and predictability. The results obtained with synthetic and real-world data demonstrate the capability of our method.

The first description of synchronization of two coupled dynamical systems (two pendulum clocks hanging from a beam) is attributed to Christiaan Huygens in 1665. In 1990, Pecora and Carroll demonstrated that chaotic systems that are unidirectionally coupled (i.e., drive-response systems) can also synchronize. By synchronization in both of the previous cases, we mean that the two systems evolve in finite time to a dynamic with a constant relationship between their states. In turn, chaotic synchronization gave rise to generalized synchronization, where now the relationship between states may be arbitrary. Precisely, our work deals with a mathematical formulation of generalized synchronization in the more realistic case of drivers perturbed by dynamical noise. In addition, we do not assume knowledge of the states but only scalar observables of them in the form of time series. We also discuss other practical issues, most importantly, a method to detect generalized synchronization for both noiseless and noisy drivers, based on recurrent neural networks. The capability of this method is successfully tested with numerical simulations and real-world data.

I. INTRODUCTION

The framework of this paper is nonlinear time series analysis, and the topic is synchronization of two unidirectionally coupled nonlinear systems and its generalization when noise is present in the driving system. By synchronization we mean generalized (or general) synchronization in the sense of Afraimovich et al.¹ and Rulkov et al.,² i.e., there is a map, called the synchronization map that transforms the states of the driving system (driver) into states of the driven system (responder), possibly with a time delay or after an initial transient time. Identical (complete, full, etc.) synchronization corresponds then to the synchronization map being the identity between two structurally equal systems; other, more interesting examples include lag, intermittent-lag, and phase synchronizations.³ Synchronization, whether identical or generalized, plays an important role in many fields of science and engineering, particularly in nonlinear dynamics,^4,5 telecommunications,^4,6,7 neuroscience,^8–10 and cryptography;^11–14 see, e.g., Pikovsky et al.,¹⁵ and Pecora and Carroll¹⁶ for overviews and historical notes.

In the noiseless or fully deterministic scenario, synchronization has been extensively studied using a number of techniques, including cross (or mutual) prediction,^17,18 conditional Lyapunov exponents,^4,19,20 replica synchronization,²¹ asymptotic stability of the responder,²² nonlinear interdependence measures,^23–25 cellular nonlinear networks,^26,27 complexity measures extracted from symbolic representations,^28–30 reservoir computing,^31,32 and more. We will use prediction because predictability is a fingerprint of determinism, i.e., functional dependence.³³

Furthermore, it is known¹⁸ that, in the case of two unidirectionally coupled nonlinear systems with a noiseless driver, there exists typically a continuous map defined from the reconstructed state space of the responder to the reconstructed state space of the driver, which was called the closeness mapping in Amigó and Hirata³⁴ and will be called the cross map here. As it turns out, the definition of synchronization amounts to the invertibility (i.e., bijectivity) of the cross map; in fact, the inverse of the cross map is the “translation” of the synchronization map (if any) from the original domains (driver and responder state spaces) to the reconstructed ones. The existence of the cross map has been used to study interdependence and causal relationships in nonlinear time series analysis.^34–37 In a nutshell, these methods harness some actual or hypothetical property of the cross map (continuity, smoothness, or local expansiveness) to reveal, given bivariate time series of a coupled dynamics, what the driving system is. We will generalize the cross and synchronization maps to multi-time versions that are well suited to the application of recurrent neural networks in time series analysis.

The main objective of this paper is the extension of the cross and synchronization maps from noiseless to noisy drivers. To model noise in the driver, we replace the dynamic of a noiseless driver with a family of driving dynamics indexed by a finite number of parameters whose values are randomly chosen, an idea called finitely parameterized stochasticity.³⁸ To implement this idea in our setting, we will use the stochastic forcing approach of Stark et al.³⁹ First, noise is formulated as an autonomous dynamical system called a shift system, whose states comprise all possible noise realizations in form of parametric sequences; the $n$ th component of a given sequence indicates which is the chosen driving dynamic at time $n$ ⁠. Second, the noisy driver is then formulated as a non-autonomous system, namely, a system forced by that shift system. As a result, our approach to synchronization in the presence of dynamical noise is based on state space reconstruction for unidirectionally coupled systems⁴⁰ and stochastic forcing³⁹ and is sufficiently general for practical purposes. We will show that the cross and synchronization maps can be extended from the noiseless to the noisy scenario by incorporating an additional dependence on the noise parameters and only for typical driver states in the case of the cross map.

In addition to discussing theoretical results in the noiseless and noisy scenarios, we also explore synchronization with synthetic and real world data. Prompted by multi-time expressions for the synchronization map, we not only use perceptrons but also long short-term memory (LSTM) nets.⁴¹ In this approach, synchronization is detected via estimation of a responder state by a contemporaneous driver state or, in case of LSTM nets, by a segment of contemporaneous and past driver states. As the benchmark in numerical simulations, we chose nearest neighbor cross prediction because it is based precisely on the existence of the cross map, so it fits very well in our approach. In fact, the continuity of the cross map entails that near neighbors of a point in the responder state space map to near neighbors of its image in the driver state space (and vice versa when synchronization sets in). Therefore, one expects that this correspondence stays if low-amplitude noise affects the driving dynamic. We also apply LSTM nets to detect coupling directionality and synchronization in electroencephalograms (EEGs) from a subject with epilepsy. The optimization of the parameters and metaparameters of our numerical tools is beyond the scope of the present work.

To address the points described above, the rest of this paper is divided into a first, theoretical, and a second, numerical, part. Thus, in Sec. II, we first review the basics of our approach to make this paper self-contained. For didactic reasons, we start with the Takens and Stark (or forced Takens) embedding theorems (Sec. II A), along with the concept of cross map (Sec. II B); then we introduce the concept of generalized synchronization (Sec. II C) and discuss its relationship with the cross map (Sec. II D). Novel concepts such as the cross and synchronization maps of higher periods are introduced for further applications in Secs. VI and VII. The presentation is rigorous from a mathematical point of view but unnecessary technicalities are avoided. Along the way, practical issues are also considered with a view to the second part of the paper. Once the traditional, noiseless scenario has been presented, the noisy scenario is set in two steps: in Sec. III, we revisit stochastic forcing and an embedding theorem that is used in the second step, Sec. IV, where a unidirectional coupling with a noisy driver is modeled as stochastic forcing. The generalization of the cross and synchronization maps to the noisy scenario is the subject of Sec. V. The theoretical concepts of Secs. II–V are illustrated and put into practice in the second part of the paper. For this purpose, in Sec. VI, we resort to two unidirectionally coupled Hénon maps, synchronization being detected with recurrent neural networks and compare the results with those obtained with nearest-neighbor cross predictability. In Sec. VII, we tackle the applicability of our tools to the analysis of real data in the form of EEGs, where noise and bidirectional coupling are the rules. In this case, the “driver” is identified by the direction with the strongest coupling. Our findings are discussed in light of results published in the specialized literature. Finally, the main contributions and conclusions of this paper are summarized in Sec. VIII.

II. THE NOISELESS SCENARIO

This section is a compact, mathematically oriented account of the cross map, synchronization, and their interplay in the absence of noise.

A. Embedding theorems

Following Stark,⁴⁰ let

Y

be a non-autonomous dynamical system (the responder) evolving under the influence of an autonomous dynamical system

X

(the driver).

X

is also called the driving or forcing system, and

Y

is the driven or forced system. In the case of discrete-time deterministic dynamical systems or flows observed at discrete times, this situation is described by the difference equations

{\begin{cases} x_{t + 1} = f (x_{t}), \\ y_{t + 1} = g (x_{t}, y_{t}), \end{cases}

(1)

where (i)

x_{t} \in M_{X}

is the state of system

X

at time

t

⁠; (ii)

y_{t} \in M_{Y}

is the state of system

Y

at time

t

⁠; (iii)

M_{X}

and

M_{Y}

are compact manifolds of dimensions

\dim_{X} \geq 0

and

\dim_{Y} \geq 1

⁠, respectively; (iv)

f : M_{X} \to M_{X}

is a

C^{1}

diffeomorphism (i.e., a

C^{1}

invertible map such that

f^{- 1}

is also

C^{1}

⁠, where

C^{1}

is shorthand for continuously differentiable); and (v)

g : M_{X} \times M_{Y} \to M_{Y}

is a

C^{1}

map such that

g (x, \cdot)

is a diffeomorphism of

M_{Y}

for every

x \in M_{X}

⁠. Alternatively, we say that there is a (unidirectional) coupling from

X

to

Y

and use the shorthand

X ⇝ Y

⁠. Since we assume that

f

and

g (x, \cdot)

are invertible, we may take

t \in Z

⁠, although in applications time series have a beginning that we will set at

t = 0

⁠.

Remark 1

By defining the map

[f, g] (x, y) = (f (x), g (x, y)),

(2)

the forced system becomes an autonomous dynamical system on the full state space

M_{X} \times M_{Y}

⁠, called the skew product of

f

and

g

⁠. Due to the properties (ii)–(v) above,

[f, g]

is a

C^{1}

diffeomorphism of the compact manifold

M_{X} \times M_{Y}

⁠.

As stated in the Introduction, we are mainly interested in nonlinear time series analysis. So, suppose further that the only information available about the systems $X$ and $Y$ are scalar observations $φ_{X} (x_{t})$ of the states $x_{t}$ and $φ_{Y} (y_{t})$ of the states $y_{t}$ ⁠, where the observation functions $φ_{X} : M_{X} \to R$ and $φ_{Y} : M_{Y} \to R$ are assumed to be $C^{1}$ ⁠. To reconstruct the state spaces of the driver $X$ and the responder $Y$ from the corresponding observations $φ_{X} (x_{t})$ and $φ_{Y} (y_{t})$ ⁠, we use the Takens and Stark theorems, respectively, which we remind below for further reference.

[Takens theorem42]

Theorem 2

[Takens theorem⁴²]

If

d \geq 2 \dim_{X} + 1

⁠, then the map

E_{f, φ_{X}} : M_{X} \to R^{d}

defined as

E_{f, φ_{X}} (x) = (φ_{X} (f^{0} (x)), φ_{X} (f^{1} (x)), \dots, φ_{X} (f^{d - 1} (x)))

(3)

is an embedding for generic

f

and

φ_{X}

⁠.

As usual, $f^{0}$ is the identity, $f^{1} = f$ ⁠, and $f^{n}$ is the $n$ th iterate of $f$ ⁠. Here, “generic $f$ and $φ_{X}$ ” formally means that the set ${f, φ_{X}}$ for which $E_{f, φ_{X}}$ is an embedding (i.e., a $C^{1}$ diffeomorphism onto its image) is open and dense in the $C^{1}$ -topology (uniform convergence of a map and its derivative) of the respective function spaces, namely, diffeomorphisms of $M_{X}$ ⁠, and $C^{1}$ maps from $M_{X}$ to $R$ ⁠. In general, a property is generic in a topological space $T$ if it holds on a residual subset $S \subset T$ ⁠, i.e., on a subset that contains a countable intersection of open sets. It turns out that an open and dense set of maps $f$ for which $E_{f, φ_{X}}$ is an embedding for generic $φ_{X}$ is built by those $C^{1}$ diffeomorphisms of $M_{X}$ that have only a finite number of periodic orbits of period less than $d$ ⁠, and the eigenvalues of each such periodic orbits are distinct (Stark,⁴⁰ Theorem 2.2).

Remark 3

Theorem 2 was generalized by Sauer et al.⁴³ in two ways. First, by replacing “generic” with “probability-one” (in the sense of prevalence). Second, by replacing the manifold $M_{X}$ by a compact invariant set $A$ that may have fractal box-counting dimension, and the restriction $d \geq 2 \dim_{X} + 1$ (which comes from Whitney’s embedding theorem⁴⁴) by $d \geq 2$ boxdim $(A) + 1$ ⁠, where boxdim $(A)$ is the box-counting dimension of $A$ ⁠.

For our purposes, we need to generalize Theorem 2 to the forced dynamic

(x_{t + 1}, y_{t + 1}) = [f, g] (x_{t}, y_{t})

defined by the diffeomorphism in the full state space

M_{X} \times M_{Y}

⁠. As before, set

{[f, g]}^{0} =

identity,

{[f, g]}^{1} = [f, g]

and

{[f, g]}^{t + 1} = [f, g] ° [f, g]^{t}

for the iterates of

[f, g]

⁠, so that

[f, g]^{t} (x, y) = (f^{t} (x), g^{(t)} (x, y)),

(4)

where

g^{(t)} (x, y) : M_{X} \times M_{Y} \to M_{Y}

is recursively defined by

g^{(0)} (x, y) = y

and

g^{(t)} (x, y) = g (f^{t - 1} (x), g^{(t - 1)} (x, y))

(5)

for

t \geq 1

⁠. Application of the Takens theorem to the skew product

[f, g]

would provide a map

E_{[f, g], φ_{X, Y}} : M_{X} \times M_{Y} \to R^{D}

⁠, with

D \geq 2 (\dim_{X} + \dim_{Y}) + 1

⁠, which would be an embedding for open dense sets of diffeomorphisms of

M_{X} \times M_{Y}

and observation

C^{1}

maps

φ_{X, Y} (x, y) : M_{X} \times M_{Y} \to R

⁠, in the

C^{1}

topology of the respective function spaces. However, what we need for applications to nonlinear time series analysis is an embedding for generic maps

f

⁠,

g

⁠, and observation maps

φ_{Y}

on

M_{Y}

⁠, and this is not guaranteed by this approach.

The generalization of the Takens theorem to forced dynamics that we need is the following due to Stark.

[forced Takens theorem40]

Theorem 4

[forced Takens theorem⁴⁰]

If

D \geq 2 (\dim_{X} + \dim_{Y}) + 1

⁠, then the map

E_{f, g, φ_{Y}} : M_{X} \times M_{Y} \to R^{D}

defined as

\begin{aligned} E_{f, g, φ_{Y}} (x, y) = (φ_{Y} (g^{(0)} (x, y)), φ_{Y} (g^{(1)} (x, y)), \dots, φ_{Y} (g^{(D - 1)} (x, y))) \end{aligned}

(6)

is an embedding for generic

f

⁠,

g

⁠, and

φ_{Y}

⁠.

Specifically, generic $g$ means that $E_{f, g, φ_{Y}}$ is an embedding for an open and dense set of diffeomorphisms $g (x, y)$ (such that $g (x, \cdot)$ is a diffeomorphism of $M_{Y}$ for every $x \in M_{X}$ ⁠) in the $C^{1}$ -topology of $M_{X} \times M_{Y}$ ⁠. In this case, an open and dense set of maps $f$ for which $E_{f, g, φ_{Y}}$ is an embedding for generic $g$ and $φ_{Y}$ is built by those $C^{1}$ diffeomorphisms of $M_{X}$ whose periodic orbits of period less than $2 d$ are isolated and have distinct eigenvalues (Stark,⁴⁰ Theorem 3.1).

B. The cross map

Hereinafter, we tacitly assume that $f$ ⁠, $g$ ⁠, $φ_{X}$ ⁠, and $φ_{Y}$ are generic in the sense of Theorems 2 and 4. Also, “smooth” stands for $C^{1}$ smoothness in the following.

Given the scalar observations

{(φ_{X} (x_{t}))}_{t \in Z}

and

{(φ_{Y} (y_{t}))}_{t \in Z}

⁠, Theorems 2 and 4 allow to “reconstruct” the (possibly unknown) dynamics of the underlying systems

X

and

Y

in the manifolds

N_{X} = E_{f, φ_{X}} (M_{X}) \subset R^{d}

and

N_{Y} = E_{f, g, φ_{Y}} (M_{X} \times M_{Y}) \subset R^{D}

called the reconstructed driver and responder state spaces, respectively, by means of the time-delay vectors

x_{t} = E_{f, φ_{X}} (x_{t}) = (φ_{X} (x_{t}), φ_{X} (x_{t + 1}), \dots, φ_{X} (x_{t + d - 1})) \in R^{d}

(7)

and

y_{t} = E_{f, g, φ_{Y}} (x_{t}, y_{t}) = (φ_{Y} (y_{t}), φ_{Y} (y_{t + 1}), \dots, φ_{Y} (y_{t + D - 1})) \in R^{D} .

(8)

In turn, the dynamics

x_{t + 1} = f (x_{t})

in

M_{X}

translates into the reconstructed driving dynamics

x_{t + 1} = (E_{f, φ_{X}} ° f ° E_{f, φ_{X}}^{- 1}) (x_{t}) =: \tilde{f} (x_{t})

(9)

in

N_{X}

⁠, while the dynamics

(x_{t + 1}, y_{t + 1}) = [f, g] (x_{t}, y_{t})

in

M_{X} \times M_{Y}

translates into the reconstructed coupled dynamics

y_{t + 1} = (E_{f, g, φ_{Y}} ° [f, g] ° E_{f, g, φ_{Y}}^{- 1}) (y_{t}) =: \tilde{[f, g]} (y_{t})

(10)

in

N_{Y}

⁠, the manifolds

N_{X}

and

N_{Y}

being diffeomorphic copies of

M_{X}

and

M_{X} \times M_{Y}

⁠, respectively. Therefore, all coordinate-independent properties of

f

and

[f, g]

can be determined in

N_{X}

and

N_{Y}

⁠.

Remark 5

Without loss of generality, it can be assumed that $d = D$ ⁠. In nonlinear time series analysis, where the underlying dynamical system is unknown, the embedding dimension of a time series is usually chosen by the method of false nearest neighbors.⁴⁵

Let

Π_{X} : M_{X} \times M_{Y} \to M_{X}

be the projection onto

M_{X}

⁠, i.e.,

Π_{X} (x, y) = x

⁠. From the diagram

\begin{array}{rcccl} M_{X} \times M_{Y} ∋ & (x_{t}, y_{t}) & \overset{Π_{X}}{⟶} & x_{t} & \in M_{X} \\ E_{f, g, φ_{Y}}^{- 1} & ↑ & ↓ & E_{f, φ_{X}} \\ N_{Y} ∋ & y_{t} & x_{t} & \in N_{X} \end{array}

(11)

along with the smoothness of the embeddings

E_{f, g, φ_{Y}}^{- 1}, E_{f, φ_{X}}

⁠, and the projection

Π_{X}

⁠, we conclude the following proposition.

Proposition 6

A unidirectional coupling

X

⇝ Y

necessarily implies the existence of a smooth map

Φ := E_{f, φ_{X}} ° Π_{X} ° E_{f, g, φ_{Y}}^{- 1} : N_{Y} \to N_{X}

(12)

called the cross map of the coupling

X

⇝ Y

⁠, which sends

y_{t}

to

x_{t}

⁠, i.e.,

x_{t} = Φ (y_{t}) .

(13)

Intuitively, Eq. (13) spells out that the responder signal carries information about the dynamics of the driver because of the time evolution law $y_{t + 1} = g (x_{t}, y_{t})$ ⁠.

Remark 7

Equation is equivalent to the existence of a map

x_{t} = Φ^{(k)} (y_{t - k}),

(14)

for any

k \in Z

⁠, where

Φ^{(k)} : N_{Y} \to N_{X}

⁠. Indeed, from

x_{t} = {\tilde{f}}^{k} (x_{t - k})

(15)

[see Eq. ] and

x_{t - k} = Φ (y_{t - k})

⁠, it follows

Φ^{(k)} = {\tilde{f}}^{k} ° Φ = E_{f, φ_{X}} ° f^{k} ° E_{f, φ_{X}}^{- 1} ° Φ

(16)

and, hence,

x_{t} = \frac{1}{K} \sum_{k = 0}^{K - 1} Φ^{(k)} (y_{t - k}) =: Φ_{K} (y_{t}, y_{t - 1}, \dots, y_{t - K + 1})

(17)

for all

K \geq 1

⁠. Note that

Φ^{(k)}

and

Φ_{K}

are continuous, and

Φ^{(0)} = Φ_{1} = Φ

⁠.

By changing the summation limits in (17), one can construct other similar multi-time expressions. For definiteness, we will use only definition (17).

Definition 8

We call the continuous map $x_{t} = Φ^{(k)} (y_{t - k})$ the cross map of order $k \in Z$ ⁠, and the continuous map $x_{t} = Φ_{K} (y_{t}, y_{t - 1}, \dots, y_{t - K + 1})$ the cross map of period $K \geq 1$ ⁠.

The continuity of the cross map $Φ$ has been used in nonlinear time series analysis to discriminate functional (deterministic, causal, etc.) relationships between observations due to coupled dynamics from statistical correlation. In its simplest version, the continuity of the cross map $x_{t} = Φ (y_{t})$ belonging to the coupling $X ⇝ Y$ implies that, given an open ball $B_{ε} (x_{t}) \subset R^{d}$ with center $x_{t}$ and arbitrary radius $ε > 0$ ⁠, there exists an open ball $B_{δ} (y_{t}) \subset R^{D}$ with center $y_{t}$ and radius $δ = δ (ε) > 0$ such that $Φ (B_{δ} (y_{t})) \subset B_{ε} (x_{t})$ ⁠. Therefore, the $k$ nearest neighbors $y_{t_{1}}, \dots, y_{t_{k}}$ of a time-delay vector $y_{t} \in N_{Y}$ in a time series ${(y_{t})}_{0 \leq t \leq T}$ of the responder are mapped by $Φ$ to close neighbors $x_{t_{1}}, \dots, x_{t_{k}}$ of the contemporaneous vector $x_{t} \in N_{Y}$ in the corresponding time series ${(x_{t})}_{0 \leq t \leq T}$ of the driver. Methods that take advantage of the continuity of $Φ$ in this way to unveil the coupling $X ⇝ Y$ include cross prediction,¹⁸^, convergent cross mapping,³⁵ and continuity scaling.³⁷

C. Generalized synchronization

According to Rulkov et al.² and Pikovsky et al.,¹⁵ the systems

X ⇝ Y

are in generalized (or general) synchronization if there exists a continuous map

h :

M_{X} \to M_{Y}

such that

y_{t} = h (x_{t}),

(18)

for all

t \in Z

⁠. That is, the responder follows the driver but in a weaker form than in identical synchronization, which corresponds to

h

being the identity (i.e.,

X

and

Y

are structurally the same and

y_{t} = x_{t}

⁠). We will also say that

Y

is synchronized to

X

if Eq. holds and call

h (x)

the synchronization map.

Therefore, in case of synchronization the full state space

M_{X} \times M_{Y}

shrinks into the subspace

{(x, y) \in M_{X} \times M_{Y} : y = h (x)}

⁠, which is the graph of the synchronization map

x \mapsto h (x)

⁠. This subspace is usually called the synchronization manifold, even when

h

is not smooth. It follows then that the projection map from

M_{X} \times M_{Y}

onto

M_{X}

⁠,

Π_{X} (x, y) = x

⁠, is invertible,

Π_{X}^{- 1} (x) = (x, h (x)),

(19)

and the range of

Π_{X}^{- 1} : M_{X} \to M_{X} \times M_{Y}

is the synchronization manifold.

Remark 9

Plug the driver dynamic

x_{t} = f (x_{t - 1})

into Eq. to derive

\begin{aligned} y_{t} = h (x_{t}) = (h ° f) (x_{t - 1}) = \dots = (h ° f^{k}) (x_{t - k}) =: h^{(k)} (x_{t - k}), \end{aligned}

(20)

where

h^{(k)} = h ° f^{k} :

M_{X} \to M_{Y}

is continuous for all

k \geq 0

⁠, and

h^{(0)} = h

⁠. Equation with

k \geq 1

corresponds to generalized synchronization for responders with an internal delay loop. When

f

is invertible (as in our case), the generalized synchronization of the systems

X

and

Y

can indistinctly be defined by

y_{t} = h (x_{t})

or, more generally, by

y_{t} = h^{(k)} (x_{t - k})

for

k \in Z

⁠; in the latter case,

h = h^{(k)} ° f^{- k}

⁠. The maps

h^{(k)} :

M_{X} \to M_{Y}

are called synchronization maps of order

k

⁠, orders

0

and

1

being the usual choices in applications.

From Eq. , it trivially follows that

y_{t} = \frac{1}{K} \sum_{k = 0}^{K - 1} h^{(k)} (x_{t - k}) =: h_{K} (x_{t}, x_{t - 1}, \dots, x_{t - K + 1}),

(21)

for all

K \geq 1

⁠, in case

Y

is synchronized with

X

⁠. By changing the summation limits in Eq. , one can construct other similar expressions. For definiteness, we will use only Eq. in this paper, so that

h_{1} = h

⁠.

Therefore, to detect synchronization of a time series ${y_{t}}_{t \geq 0}$ with another time series ${x_{t}}_{t \geq 0}$ ⁠, we can look for functional dependencies of the form (21) with $K > 1$ rather than $y_{t} = h (x_{t})$ ⁠. If ${x_{t}}_{t \geq 0}$ is a deterministic time series (i.e., $x_{t + 1} = f (x_{t})$ ⁠) and ${y_{t}}_{t \geq 0}$ is synchronized with it, then Eq. (21) holds with a continuous map $h_{K} : M_{X}^{K} \to M_{Y}$ such that $h_{K} (x_{t}, f^{- 1} (x_{t}), \dots, f^{- K + 1} (x_{t})) = h (x_{t})$ ⁠. The point is that, in time series analysis, multi-time dependencies like (21) can be efficiently detected by recurrent neural nets, as we discuss in Sec. VI.

Definition 10

We call the continuous map $y_{t} = h_{K} (x_{t}, \dots, x_{t - K + 1})$ in Eq. (21), the synchronization map of period $K \geq 1$ ⁠.

The synchronization maps of order

k

⁠,

h^{(k)} = h ° f^{k}

⁠, satisfy a number of straightforward relations involving also the function

g (x, y)

⁠. Indeed, in case of synchronization, the dynamic of

X ⇝ Y

simplifies to

{\begin{cases} x_{t + 1} = f (x_{t}), \\ y_{t + 1} = g (x_{t}, h (x_{t})) . \end{cases}

(22)

Comparing with

y_{t + 1} = h (f (x_{t}))

shows that

h (x)

fulfills the functional relation

h (f (x)) = g (x, h (x)) .

(23)

Replace

x

with

f^{k} (x)

in Eq. to obtain

h^{(k + 1)} (x) = h^{(k)} (f (x)) = g (f^{k} (x), h^{(k)} (x)) .

(24)

Recursion of Eq. leads to alternative formulas for synchronization maps of arbitrary periods involving function

g

⁠.

Contingent upon the structure of $g (x, y)$ ⁠, the synchronization map $h (x)$ can sometimes be written in closed form; see Pikovsky et al.¹⁵ and Parlitz⁴⁶ for an example with a baker map. Interestingly, the parameters of that example can be fine-tuned so that the cross sections $x^{(2)} = const$ of $h (x^{(1)}, x^{(2)})$ are Weierstrass functions, i.e., continuous functions that are nowhere differentiable.

The definition of synchronization can be weakened by requiring condition only asymptotically. In more formal terms, we say that the responder

Y

is asymptotically synchronized to the driver

X

if there exists a continuous map

h :

M_{X} \to M_{Y}

such that

lim_{t \to \infty} ‖ y_{t} - h (x_{t}) ‖ = 0,

(25)

where

‖ \cdot ‖

is a distance in

M_{Y}

⁠. In this case, the synchronization manifold becomes an attracting set in

M_{X} \times M_{Y}

⁠.

A direct consequence of asymptotic synchronization is the asymptotic stability of the responder. We say that the responder

Y

is asymptotically stable if all orbits converge to the same orbit regardless of the initial condition, that is, if given two responses

{(y_{t})}_{t \geq 0}

and

{({\tilde{y}}_{t})}_{t \geq 0}

to a signal

{(x_{t})}_{t \geq 0}

from the driver with different initial conditions

y_{0} \neq {\tilde{y}}_{0}

⁠, then

lim_{t \to \infty} ‖ y_{t} - {\tilde{y}}_{t} ‖ = 0.

(26)

Asymptotic stability of the responder is weaker than asymptotic synchronization because the existence of a hypothetical synchronization map does not follow from Eq. (26). This is the case, for example, when a periodic driver has gone through a period-doubling bifurcation⁴⁷ or there is a multistability in the responder, i.e., a driver signal ${(x_{t})}_{t \geq 0}$ can elicit two or more stable responses.¹⁵

On the other hand, the asymptotic stability of the responder provides a simple method to test synchronization called the auxiliary system method.²¹ This method boils down to check Eq. (26) for two initial conditions $y_{0} \neq {\tilde{y}}_{0}$ ⁠; if (26) does not hold, then $Y$ is not synchronized to $X$ ⁠.

D. Relationship between generalized synchronization and the cross map

According to Proposition 6, a coupling $X ⇝ Y$ implies the existence of the cross map $x_{t} = Φ (y_{t})$ ⁠, whereas the synchronization map $y_{t} = h (x_{t})$ exists in seemingly exceptional cases (unless the coupling is strong enough). Despite this notable difference, both maps are closely related, as we will now see.

Let

Π_{Y} (x, y) = y

be the projection map from

M_{X} \times M_{Y}

onto

M_{Y}

⁠. On the one hand, from the diagram

\begin{array}{rcccl} M_{X} ∋ & x_{t} & \overset{h}{⟶} & y_{t} & \in M_{Y} \\ E_{f, φ_{X}} & ↓ & ↑ & Π_{Y} ° E_{f, g, φ_{Y}}^{- 1} \\ N_{X} ∋ & x_{t} & \overset{Φ}{⟵} & y_{t} & \in N_{Y} \end{array},

(27)

we have that if

Φ

is invertible, then

h

exists and

h = Π_{Y} ° E_{f, g, φ_{Y}}^{- 1} ° Φ^{- 1} ° E_{f, φ_{X}} .

(28)

On the other hand, Eq. spells out that, if

h

exists, then

Π_{X}

is invertible. From Eq. , it follows then that

Φ^{- 1} = E_{f, g, φ_{Y}} ° Π_{X}^{- 1} ° E_{f, φ_{X}}^{- 1} .

(29)

The bottom line of Eqs. and is the following.

Proposition 11

The systems $X ⇝ Y$ are synchronized if and only if the cross map $x_{t} = Φ (y_{t})$ is invertible and bicontinuous (i.e., $Φ^{- 1}$ is continuous). If this case, the synchronization map $y_{t} = h (x_{t})$ and $y_{t} = Φ^{- 1} (x_{t})$ are related through the expressions (28) and (29), where $Π_{X}^{- 1} (x) = (x, h (x))$ ⁠.

In other words, $Φ^{- 1} : N_{X} \to N_{Y}$ is the synchronization map of the systems $X ⇝ Y$ in the reconstructed state spaces (if it exists and is continuous); see Pecora et al.⁴⁸ for numerical methods to test whether two time series are related by a map with properties such as continuity, invertibility, smoothness, and more. As a rule, the relationship $x \mapsto (x, y)$ is multivalued owing to folds in the manifold $M_{X} \times M_{Y}$ ⁠, so generalized synchronization is rather an exception. Multivalued synchronization maps, corresponding to noninvertible cross maps, have been considered, e.g., in Rulkov et al.⁴⁹ and Parlitz.⁴⁶

To lift

y_{t} = h_{K} (x_{t}, x_{t - 1}, \dots, x_{t - K + 1})

⁠, the synchronization map of period

K

, to the reconstructed state spaces

N_{X}

and

N_{Y}

⁠, use the reconstructed driver dynamic

\tilde{f} : N_{X} \to N_{X}

defined in Eq. to derive

y_{t} = Φ^{- 1} (x_{t}) = (Φ^{- 1} ° \tilde{f}) (x_{t - 1}) = \dots = (Φ^{- 1} ° {\tilde{f}}^{k}) (x_{t - k}),

(30)

where

k \in Z

and

Φ^{- 1} ° {\tilde{f}}^{k} = Φ^{- 1} ° E_{f, φ_{X}} ° f^{k} ° E_{f, φ_{X}}^{- 1} = {(Φ^{(- k)})}^{- 1},

(31)

by the definition of the cross map of order

k

⁠, Eq. . Therefore, the synchronization map of period

K \geq 1

becomes

y_{t} = \frac{1}{K} \sum_{k = 0}^{K - 1} {(Φ^{(- k)})}^{- 1} (x_{t - k}) =: H_{K} (x_{t}, x_{t - 1}, \dots, x_{t - K + 1})

(32)

in the reconstructed spaces. Note that

H_{1} = Φ^{- 1}

⁠.

Definition 12

The continuous map $y_{t} = H_{K} (x_{t}, x_{t - 1}, \dots, x_{t - K + 1})$ defined in Eq. (32) will be called the reconstructed synchronization map of period $K \geq 1$ ⁠.

We will harness $H_{K}$ with $K > 1$ in the applications with synthetic data (Sec. VI) and real world data (Sec. VII) via recurrent neural networks. We remark already at this point that, unlike the synthetic data of Sec. VI, real world data are generally bidirectionally coupled, as happens with the EEGs of Sec. VII. Following the standard approach, we will measure the coupling strength between pairs of EEGs in both directions, the “driver” being identified by the direction with the strongest coupling. In case of equal strengths, the systems are assumed to be synchronized.

As mentioned in Sec. II B, in case of unidirectional coupling (the framework of this paper) the relationship

x_{t} = Φ (y_{t})

due to the coupled dynamic

X ⇝ Y

can be unveiled numerically from the time series

{(x_{t})}_{0 \leq t \leq T}

and

{(y_{t})}_{0 \leq t \leq T}

⁠.^35,37,48 Thus, in the method of nearest-neighbor cross prediction, one estimates

x_{t}

or

x_{t + 1}

based on the nearest neighbors of

y_{t}

for any

0 \leq t \leq T

to test for the existence of

Φ

⁠. Likewise, if

Φ^{- 1}

exists and is continuous (i.e.,

X

and

Y

are synchronized), then one can also discern the inverse relationship

y_{t} = Φ^{- 1} (x_{t})

by the same techniques. As a matter of fact, in the case of a bijective and bicontinuous

Φ

⁠, there is a one-to-one relation between the neighborhoods of nearest neighbors of

y_{t}

and

x_{t}

⁠, so, if

Q_{T} (Φ)

and

Q_{T} (Φ^{- 1})

are fidelity metrics of the respective estimations, then

Q_{T} (Φ) - Q_{T} (Φ^{- 1}) ≃ 0,

(33)

where it applies that the longer the time series, the better the predictions and, hence, the smaller

Q_{T} (Φ) - Q_{T} (Φ^{- 1})

⁠. In other words, the continuity of the cross map and its inverse can be exploited via to test two time series

{(x_{t})}_{0 \leq t \leq T}

and

{(y_{t})}_{0 \leq t \leq T}

for general synchronization.

Remark 13

The nonexistence of the cross map can uncover common drivers. Indeed, if $Z ⇝ X$ and $Z ⇝ Y$ ⁠, then $z_{t} = Φ (x_{t})$ and $z_{t} = \tilde{Φ} (y_{t})$ ⁠, so that $Φ (x_{t}) = \tilde{Φ} (y_{t})$ ⁠. Here, $Φ : N_{X} \to N_{Z}$ and $\tilde{Φ} : N_{Y} \to N_{Z}$ are the cross maps associated to the coupled dynamics $Z ⇝ X$ and $Z ⇝ Y$ ⁠, respectively. However, there is no cross map between $N_{X}$ and $N_{Y}$ ⁠, unless $Φ$ or $\tilde{Φ}$ is invertible and continuous [so that $x_{t} = (Φ^{- 1} ° \tilde{Φ}) (y_{t})$ or $y_{t} = ({\tilde{Φ}}^{- 1} ° Φ) (x_{t})$ ], in which case $X$ or $Y$ is synchronized with $Z$ ⁠.

To wrap up this section, let us point out that the diagram is a particularization of the diagram

\begin{array}{rcccl} M_{X} ∋ & x_{t} & \overset{h^{(k)}}{⟶} & y_{t + k} & \in M_{Y} \\ E_{f, φ_{X}} & ↓ & ↑ & Π_{Y} ° E_{f, g, φ_{Y}}^{- 1} \\ N_{X} ∋ & x_{t} & \overset{Φ^{(- k)}}{⟵} & y_{t + k} & \in N_{Y} \end{array}

(34)

to

k = 0

⁠, where

h^{(k)}

is the synchronization map of order

k

[Eq. ,

h^{(0)} = h

⁠) and

Φ^{(- k)}

is the cross map of order

- k

[Eq. ,

Φ^{(0)} = Φ

]. Thus, Eq. is the special case

k = 0

of the relationship

\begin{aligned} h^{(k)} & = Π_{Y} ° E_{f, g, φ_{Y}}^{- 1} ° {(Φ^{(- k)})}^{- 1} ° E_{f, φ_{X}} \\ = Π_{Y} ° E_{f, g, φ_{Y}}^{- 1} ° Φ^{- 1} ° E_{f, φ_{X}} ° f^{k}, \end{aligned}

(35)

where we used and in the second line. In this regard, note that

h^{(k)}

is invertible if and only if

h

is invertible [since

h^{(k)} = h ° f^{k}

by Eq. ] and, likewise,

Φ^{(- k)}

is invertible if and only if

Φ

is invertible [since

Φ^{(- k)} = E_{f, φ_{X}} ° f^{- k} ° E_{f, φ_{X}}^{- 1} ° Φ

by Eq. ].

III. DYNAMICAL NOISE AS STOCHASTIC FORCING

Random or “noisy” dynamical systems can be modeled in different ways, from the perhaps simplest ones, such as switching systems^50,51 and iterated function systems⁵² to nonautonomous dynamical systems⁵³ and full-fledged random dynamical systems, described by random differential equations.⁵⁴ In our setting, a natural way to turn a noiseless dynamic, say, $x_{t + 1} = f (x_{t})$ on a compact manifold $M_{X}$ ⁠, into a noisy one is to replace the map $f$ with a family of maps ${f_{ω_{t}}}_{t \in Z}$ ⁠, where the index $ω_{t}$ is a (possibly multi-component) parameter belonging to a suitable space that is randomly chosen at each discrete time $t$ ⁠. For example, $f_{ω_{t}} (x_{t}) = f (x_{t}) + ω_{t}$ ⁠, $ω_{t} \in M_{X}$ ⁠, models additive dynamical noise, while $f_{ω_{t}} (x_{t}) = ω_{t} f (x_{t})$ models multiplicative dynamical noise. This approach, sometimes called “finitely parameterized stochasticity”³⁸ is sufficient for most practical applications.⁵⁵ So, the term noisy dynamical system will refer hereafter to such implementation of dynamical noise via stochastic processes in the parameter space; our parameter spaces will be compact topological sets.

At this point, we recall that each stationary (discrete-time) stochastic process corresponds in a canonical way to a so-called shift system, which is a dynamical system whose states are the realizations of the stochastic process considered.⁵⁶ In other words, stationary stochastic processes can be modeled as autonomous dynamical systems. It is, therefore, not surprising that shift systems allow to formulate noisy dynamical systems as forced systems in a manner formally similar to the noiseless case. However, before getting to that point, we need to introduce the concepts and notation.

Let

Ω

be a compact topological space of parameters and let

Ω^{Z}

be the set of all two-sided sequences of elements of

Ω

⁠,

ω = (\dots, ω_{- t}, \dots, ω_{- 1}, ω_{0}, ω_{1}, \dots, ω_{t}, \dots)

endowed with the product topology. As a result,

Ω^{Z}

is a compact topological space too. Furthermore, let

σ : Ω^{Z} \to Ω^{Z}

be the (left) shift map,

σ (ω) = (\dots, ω_{- t + 1}, \dots, ω_{0}, \overset{*}{ω_{1}}, ω_{2}, \dots, ω_{t + 1}, \dots),

where the asterisk marks the zeroth component; component-wise,

{[σ (ω)]}_{t} = ω_{t + 1}

for all

t \in Z

⁠. The shift map is a homeomorphism of

Ω^{Z}

⁠.

In addition, the continuous (or topological) dynamical system

(Ω^{Z}, σ)

can be further promoted to a measure-preserving dynamical system by introducing a

σ

-invariant (probability) measure

μ

on the Borel sigma-algebra

B^{\times}

of

Ω^{Z}

via the finite-dimensional probability distributions of the given or desired

Ω

-valued stochastic process.⁵⁶ For instance, the

n

-dimensional marginal probabilities

P (ω_{i_{1}} \in B_{i_{1}}, \dots, ω_{i_{n}} \in B_{i_{n}})

on

(Ω^{Z}, B^{\times})

are given by

μ_{n} (B_{i_{1}} \times \dots \times B_{i_{n}}) := μ (B_{i_{1}} \times \dots \times B_{i_{n}} \times \prod^{i \neq i_{1}, \dots, i_{n}} Ω),

(36)

where

B_{i_{1}}, \dots, B_{i_{n}}

are Borel sets (e.g., open sets) of

Ω

⁠. The resulting dynamical system

Σ = (Ω^{Z}, B^{\times}, μ, σ)

is the shift system mentioned above. When the measure

μ

consists of finitely many atoms, then

{f_{ω_{t}}}

is an iterated function system.⁵² Product measures

μ_{n} (B_{i_{1}} \times B_{i_{2}} \times \dots \times B_{i_{n}}) = μ_{1} (B_{i_{1}}) μ_{1} (B_{i_{2}}) \dots μ_{1} (B_{i_{n}})

correspond to independent (memoryless) processes such as coin tossing and white noise.

Following Stark et al.,³⁹ a noisy dynamical system

X

is then modeled by the skew product

{\begin{cases} ω_{t + 1} = {[σ (ω)]}_{t}, \\ x_{t + 1} = f (ω_{t}, x_{t}) = f_{ω_{t}} (x_{t}), \end{cases}

(37)

where we suppose that

f_{ω_{t}} = f (ω_{t}, \cdot) : M_{X} \to M_{X}

is a diffeomorphism for all

ω_{t} \in Ω

⁠. Alternatively,

{\begin{cases} ω_{t + 1} = {[σ^{t} (ω)]}_{0}, \\ x_{t + 1} = f ({[σ^{t} (ω)]}_{0}, x_{t}) = f_{{[σ^{t} (ω)]}_{0}} (x_{t}), \end{cases}

(38)

i.e., the parameter of the dynamic at time

t

is the

0

-component of the shifted sequence

σ^{t} (ω)

⁠.

Due to the formal similarity of Eq. (37) with Eq. (1) for the forced dynamic $X ⇝ Y$ ⁠, the modeling (37) of a noisy dynamical system is called stochastic forcing.³⁹ Indeed, here we have $Σ ⇝ X$ ⁠, where the shift system $Σ = (Ω^{Z}, B^{\times}, μ, σ)$ is also an autonomous dynamical system and $X$ is randomly forced by $Σ$ since $x_{t + 1} = f (ω_{t}, x_{t})$ ⁠. This parallelism also carries over to the embedding $E_{f, φ_{X}} : M_{X} \to R^{d}$ ⁠, Eq. (3), as follows.

Let

ω = {(ω_{t})}_{t \in Z}

be a two-sided sequence of points in

Ω

and set

f_{ω_{t}, \dots, ω_{0}} (x) = (f_{ω_{t}} ° \dots ° f_{ω_{0}}) (x),

(39)

for all

x \in M_{X}

⁠, so

f_{ω_{t}, \dots, ω_{0}} : M_{X} \to M_{X}

for all

t \geq 0

⁠. For every

ω

⁠, define the map

E_{f, φ_{X}, ω} : M_{X} \to R^{d}

as

E_{f, φ_{X}, ω} (x) = (φ_{X} (x), φ_{X} (f_{ω_{0}} (x)), \dots, φ_{X} (f_{ω_{d - 2}, . ., ω_{0}} (x))) .

(40)

Note that

E_{f, φ_{X}, ω}

actually depends on the

d - 1

parameters

ω_{0}, ω_{1}, \dots, ω_{d - 2}

⁠.

[Stark et al.39]

Theorem 14

[Stark et al.³⁹]

If $d \geq 2 \dim_{X} + 1$ ⁠, then there exists a residual set of $(f, φ_{X})$ such that for any $(f, φ_{X})$ in this set there is an open dense set of sequences $ω \in Ω^{Z}$ such that the map $E_{f, φ_{X}, ω}$ is an embedding.

Finally, let us point out that Theorem 14 generalizes readily to the case of noisy observations. For example, if the observation function

φ_{X} (x)

is replaced by the noisy observation function

φ_{X, η} (x) = φ_{X} (x, η)

⁠, where

η \in {(Ω^{'})}^{Z}

⁠,

Ω^{'}

is a compact set, and

Σ^{'}

is the corresponding shift space, then the map

\begin{aligned} E_{f, φ_{X}, ω, η} (x) = (φ_{X, η_{0}} (x), φ_{X, η_{1}} (f_{ω_{0}} (x)), \dots, φ_{X, η_{d - 1}} (f_{ω_{d - 2}, . ., ω_{0}} (x))) \end{aligned}

(41)

is an embedding for generic

(ω, η) \in Σ \times Σ^{'}

⁠. See Stark et al.³⁹ for more detail and other possibilities. Therefore, we may assume hereafter that the observations are noiseless for notational simplicity.

IV. COUPLED DYNAMICS AND NOISE

Next, we show that the skew product includes the case of two unidirectionally coupled systems

X ⇝ Y

⁠, where the driver is a noisy dynamical system, namely,

{\begin{cases} x_{t + 1} = f_{ω_{t}} (x_{t}), \\ y_{t + 1} = g (x_{t}, y_{t}) . \end{cases}

(42)

Equivalently,

x_{t + 1} = f_{ω_{t}, ω_{t - 1}, \dots, ω_{0}} (x_{0}),

(43)

see Eq. , and

y_{t + 1} = (g_{x_{t}} ° g_{x_{t - 1}} ° \dots ° g_{x_{0}}) (y_{0}) = g_{x_{t}, x_{t - 1}, \dots, x_{0}} (y_{0})

(44)

for

t \geq 0

⁠, where

g_{x} (y) := g (x, y)

⁠.

Remark 15

Of course, if

ω_{t} = ω_{0}

for all times

t

⁠, then we recover the noiseless case with

f := f_{ω_{0}}

⁠,

y_{0} = g^{(0)} (x_{0}, y_{0})

and

g_{x_{t - 1}, \dots, x_{0}} (y_{0}) = g_{x_{t - 1}, \dots, x_{1}} (g (x_{0}, y_{0})) = g^{(t)} (x_{0}, y_{0})

(45)

for

t \geq 1

⁠; see Eqs. and .

Some basic facts about the noisy driving dynamic $x_{t + 1} = f (ω_{t}, x_{t})$ follow.

Fact 1. Since the parametric sequence

ω = {(ω_{t})}_{t \in Z}

is a trajectory of an

Ω

-valued random process modeled by the shift space

Σ = (Ω^{Z}, B^{\times}, μ, σ)

⁠, the noisy orbit

\begin{aligned} ξ (x, ω) & = (\dots, x, f_{ω_{0}} (x), f_{ω_{1}, ω_{0}} (x), \dots, f_{ω_{t - 1}, \dots, ω_{0}} (x), \dots .) \\ = (\dots, x_{0}, x_{1}, x_{2}, \dots, x_{t}, \dots) \in M_{X}^{Z} \end{aligned}

(46)

is a trajectory of an

M_{X}

-valued random process. In general, i.i.d. parametric sequences

ω

(commonly used in applications) do not generate i.i.d. noisy orbits

ξ = ξ (x, ω)

⁠. Additive noise

x_{t + 1} = f (x_{t}) + ω_{t}

is a plain example when the invariant measure of

f : M_{X} \to M_{X}

is not uniform.

Fact 2. Since the $ω_{t}$ ’s are the outcomes of a stationary process, the $x_{t}$ ’s are also the outcomes of a stationary process. Indeed, the definition $x_{t + 1} = f (ω_{t}, x_{t})$ is time-invariant due to the stationarity in the generation of the $ω_{t}$ ’s.

Fact 3. Under additional assumptions, $ξ (x, ω) = {(x_{t})}_{t \in Z}$ can match any arbitrary stationary sequence in $M_{X}$ (i.e., a trajectory of a stationary $M_{X}$ -valued random process) by fine tuning the sequence $ω$ ⁠. For example, assume the following mild proviso.

Condition 16

$f (\cdot, x) = f_{x} : M_{P} \to M_{X}$ is an embedding for each $x \in M_{X}$ ⁠.

Under this condition, given $x_{t}$ ⁠, the relationship between $ω_{t}$ and $x_{t + 1}$ is one-to-one for each $t$ ⁠, which implies that the equation $f_{x_{t}} (ω_{t}) := f (ω_{t}, x_{t}) = x_{t + 1}$ can be solved for $ω_{t}$ in a unique way. Therefore, the noisy orbit $ξ (x, ω)$ can be recursively transformed into any stationary sequence $η \in M_{X}^{Z}$ by choosing $x_{0} = η_{0}$ and $ω_{t}$ as the unique solution of $f (ω_{t}, x_{t}) = η_{t + 1}$ for $t = 0, 1, 2, \dots$ and $t = - 1, - 2, \dots$ ⁠. Henceforth, we assume that Condition 16 is met.

Fact 4. In particular, by Condition 16, the relationship between $ξ = {x_{t}}_{t \in Z}$ and ${x_{0}, ω}$ is one-to-one, i.e., the function $(x, ω) \mapsto ξ$ is invertible, where $ξ_{0} = x$ ⁠. Therefore, we may indistinctly talk of $x$ and $ω$ ⁠, or $ξ = {(x_{t})}_{t \in Z}$ ⁠. In practice, one chooses $ω$ so that the noisy orbit $ξ$ of $x$ deviates from the noiseless orbit by small perturbations.

By Fact 2, we can view the noisy dynamic (42) as stochastic forcing, the compact manifold $M_{X}$ being the parameter set and the orbits $ξ = {(x_{t})}_{t \in Z}$ of the noisy driver playing the role of the parameter sequences $ω = {(ω_{t})}_{t \in Z}$ ⁠. This being the case, replace in Theorem 14(i) the sequence $ω \in Ω^{Z}$ with the noisy orbit $ξ = {(x_{t})}_{t \in Z} \in M_{X}^{Z}$ ⁠, (ii) the map $f : Ω \times M_{X} \to M_{X}$ with $g : M_{X} \times$ $M_{Y} \to M_{Y}$ ⁠, (iii) $f_{ω_{t}} (x_{t}) = f (ω_{t}, x_{t})$ with $g_{x_{t}} (y_{t}) = g (x_{t}, y_{t})$ ⁠, (iv) $f_{ω_{t}, ω_{t - 1}, . ., ω_{0}}$ with $g_{x_{t}, x_{t - 1}, \dots, x_{0}}$ ⁠, and (v) $φ_{X}$ with $φ_{Y}$ ⁠, to derive the following result.

Theorem 17

If

δ \geq 2 \dim_{Y} + 1

⁠, then there exists a residual set of

(g, φ_{Y})

such that for any

(g, φ_{Y})

in this set there is an open dense set of sequences

ξ = {(x_{t})}_{t \in Z} \in M_{X}^{Z}

such that the map

E_{g, φ_{Y}, ξ} : M_{Y} \to R^{δ}

defined by

E_{g, φ_{Y}, ξ} (y) = (φ_{Y} (y), φ_{Y} (g_{x_{0}} (y)), \dots, φ_{Y} (g_{x_{δ - 2}, . ., x_{0}} (y)))

(47)

is an embedding.

According to Eq. (47), $E_{g, φ_{Y}, ξ} : M_{Y} \to R^{δ}$ depends actually on the $δ - 1$ parameters $x_{0}, x_{1}, \dots, x_{δ - 2}$ ⁠. The points $(x_{0}, \dots, x_{δ - 2})$ are dense in the finite-dimensional manifold $M_{X}^{δ - 1}$ if and only if the points $x_{k}$ are dense in $M_{X}$ for each $0 \leq k \leq δ - 2$ ⁠. It follows that $E_{g, φ_{Y}, ξ}$ is an embedding for a residual set of $(g, φ_{Y})$ and dense sets of points ${x_{0}, \dots, x_{δ - 2}}$ in $M_{X}$ ⁠.

Remark 18

One can extend the map

E_{g, φ_{Y}, ξ}

from

M_{Y}

to

M_{X} \times M_{Y}

by defining

E_{g, φ_{Y}, x_{1}, \dots, x_{δ - 2}} : {x_{0}} \times M_{Y} \to R^{δ}

as

E_{g, φ_{Y}, x_{1}, \dots, x_{δ - 2}} (x_{0}, y) := E_{g, φ_{Y}, ξ} (y) .

(48)

Yet,

E_{g, φ_{Y}, x_{1}, \dots, x_{δ - 2}}

does not allow to reconstruct the full state space

M_{X} \times M_{Y}

because, according to Theorem 17, in general, this map is an embedding only for a dense set of points

(x_{0}, y) \in M_{X} \times M_{Y}

⁠. Nevertheless, this result can be useful in applications where one can assume

x_{0}

to be fixed and generic, like in time series analysis.

V. THE NOISY SCENARIO

In this section, we discuss some changes and limitations introduced by noise in the conventional framework of Sec. II. Since the driver dynamic now explicitly depends on time through the noise, so do the main concepts like state reconstruction, cross map, and synchronization map.

A. State reconstruction

Let

ω \in Ω^{Z}

be a parametric sequence and suppose

d \geq 2 \dim_{X} + 1

⁠. Then, according to Theorem 14, the map

E_{f, φ_{X}, ω} : M_{X} \to R^{d}

defined in Eq. is generically an embedding. Similar to the noiseless case, we define the manifolds

N_{X, ω} = E_{f, φ_{X}, ω} (M_{X}) \subset R^{d}

(49)

(each one diffeomorphic to

M_{X}

⁠), and the noisy time-delay vectors

\begin{aligned} x_{t} & = E_{f, φ_{X}, σ^{t} (ω)} (x_{t}) \\ = (φ_{X} (x_{t}), φ_{X} (x_{t + 1}), \dots, φ_{X} (x_{t + d - 1})) \in N_{X, σ^{t} (ω)}, \end{aligned}

(50)

with

x_{t + k} = f_{ω_{t + k - 1}, \dots, ω_{t}} (x_{t})

for

k \geq 1

⁠. Then, the driver dynamics

x_{t + 1} = f_{ω_{t}} (x_{t})

translate to

x_{t + 1} = F_{σ^{t} (ω)} (x_{t})

(51)

in the reconstructed state spaces, where the map

F_{ω} : N_{X, ω} \to N_{X, σ (ω)}

defined as

F_{ω} = E_{f, φ_{X}, σ (ω)} ° f_{ω_{0}} ° E_{f, φ_{X}, ω}^{- 1}

(52)

is a diffeomorphism, provided that

E_{f, φ_{X}, ω}

and

E_{f, φ_{X}, σ (ω)}

are embeddings. At variance with the noiseless case, the reconstructed dynamic

x_{t} \mapsto x_{t + 1}

hops from a diffeomorphic copy

N_{X, σ^{t} (ω)}

of

M_{X}

to another diffeomorphic copy

N_{X, σ^{t + 1} (ω)}

⁠.

Likewise, let

ξ = ξ (x, ω) \in M_{X}^{Z}

be a noisy orbit of

x = x_{0}

[Eq. ] and suppose

δ \geq 2 \dim_{Y} + 1

⁠. Then, according to Theorem 17, the map

E_{g, φ_{Y}, ξ} : M_{Y} \to R^{δ}

defined in Eq. is generically an embedding. Define the manifolds

N_{Y, ξ} = E_{g, φ_{Y}, ξ} (M_{Y}) \subset R^{δ}

(53)

(each one diffeomorphic to

M_{Y}

⁠) and the noisy time-delay vectors

\begin{aligned} y_{t} & = E_{g, φ_{Y}, σ^{t} (ξ)} (y_{t}) \\ = (φ_{Y} (y_{t}), φ_{Y} (y_{t + 1}), \dots, φ_{Y} (y_{t + δ - 1})) \in N_{Y, σ^{t} (ξ)}, \end{aligned}

(54)

with

y_{t + 1} = g (x_{t}, y_{t}) =: g_{x_{t}} (y_{t})

and

y_{t + k} = g (x_{t + k - 1}, g_{x_{t + k - 2}, \dots, x_{t + 1}, x_{t}} (y_{t})) =: g_{x_{t + k - 1}, \dots, x_{t + 1}, x_{t}} (y_{t}),

(55)

for

k \geq 2

⁠. Then, similar to , the map

G_{ξ} : N_{Y, ξ} \to N_{Y, σ (ξ)}

defined as

G_{ξ} = E_{g, φ_{Y}, σ (ξ)} ° g_{ξ_{0}} ° E_{g, φ_{Y}, ξ}^{- 1}

(56)

is a diffeomorphism, provided that

E_{g, φ_{y}, ξ}

and

E_{g, φ_{y}, σ (ξ)}

are embeddings, and it holds

y_{t + 1} = G_{σ^{t} (ξ)} (y_{t}) .

(57)

Again, the range of

G_{σ^{t} (ξ)}

depends on

t

through

σ^{t} (ξ)

⁠, but all of them are diffeomorphic copies of

M_{Y} .

B. Cross map

The definition of the cross map of the systems $X ⇝ Y$ in (12) hinges on the reconstruction of both the driver state space $M_{X}$ and the full state space $M_{X} \times M_{Y}$ ⁠. However, according to Remark 18, the latter reconstruction is generally only possible in the noisy scenario for a dense set of driver states.

This being the case, we are going to define the cross map

x_{t} = Φ_{σ^{t} (ω)} (y_{t})

for two time series

{(x_{t})}_{t \geq 0}

and

{(y_{t})}_{t \geq 0}

of time-delay vectors obtained from a noisy orbit

ξ = ξ (x_{0}, ω)

of the driver

X

[Eq. ] and the corresponding response from the system

Y

⁠, respectively. This limited approach suffices for the needs of time series analysis, where the focus in practice is on (finite segments of) single orbits rather than on manifolds, and points and parameters can be considered generic. For ease of notation, we will write

E_{g, φ_{Y}, σ^{t} (ω)} (x_{t}, y_{t}) = E_{g, φ_{Y}, x_{t + 1}, \dots, x_{t + δ - 2}} (x_{t}, y_{t}) = y_{t},

(58)

see Eq. , since the relation

σ^{t} (ξ) \leftrightarrow (x_{t}, σ^{t} (ω))

is one-to-one by Condition 16 in Sec. IV.

To define the cross map in the presence of dynamical noise,

x_{t} = Φ_{σ^{t} (ω)} (y_{t})

⁠, we mimic the definition of the cross map in Eq. in the form

\begin{array}{rcccl} M_{X} \times M_{Y} ∋ & (x_{t}, y_{t}) & \overset{Π_{X}}{⟶} & x_{t} & \in M_{X} \\ E_{g, φ_{Y}, σ^{t} (ω)}^{- 1} & ↑ & ↓ & E_{f, φ_{X}, σ^{t} (ω)} \\ N_{Y, σ^{t} (ξ)} ∋ & y_{t} & x_{t} & \in N_{X, σ^{t} (ω)}, \end{array}

(59)

under the assumption that

E_{g, φ_{Y}, σ^{t} (ω)} (x_{t}, y_{t})

is an embedding for the considered states

x_{t} \in M_{X}

⁠. Hence,

x_{t} = Φ_{σ^{t} (ω)} (y_{t}) := (E_{f, φ_{X}, σ^{t} (ω)} ° Π_{X} ° E_{g, φ_{Y}, σ^{t} (ω)}^{- 1}) (y_{t}) .

(60)

Let us check that $Φ_{σ^{t} (ω)} (y_{t})$ becomes $Φ (y_{t})$ ⁠, Eq. (12), when the noise is switched off in Eq. (60), i.e., when $ω = \bar{ω}$ with ${\bar{ω}}_{t} = ω_{0}$ for all $t \in Z$ ⁠. We suppose that the maps $E_{f, φ_{X}, σ^{t} (ω)}$ and $E_{g, φ_{Y}, σ^{t} (ω)}$ are embeddings for $ω = \bar{ω}$ ⁠.

In that case,

E_{f, φ_{X}, σ^{t} (\bar{ω})} (x_{t}) = E_{f, φ_{X}} (x_{t})

with

f := f_{ω_{0}}

⁠; see Eqs. and . Similarly, by Eqs. and with

D = δ - 1

⁠, and setting

\bar{ξ} = (x_{0}, \bar{ω}) = {(f^{t} (x_{0}))}_{t \in Z}

⁠,

\begin{aligned} E_{g, φ_{Y}, σ^{t} (\bar{ω})} (x_{t}, y_{t}) \\ = E_{g, φ_{Y}, σ^{t} (\bar{ξ})} (y_{t}) \\ = (φ_{Y} (y_{t}), φ_{Y} (g_{x_{t}} (y_{t})), \dots, φ_{Y} (g_{x_{t + δ - 2}, . ., x_{t}} (y_{t}))) \\ = (φ_{Y} (g^{(0)} (x_{t}, y_{t})), φ_{Y} (g^{(1)} (x_{t}, y_{t})), \dots, φ_{Y} (g^{(δ - 1)} (x_{t}, y_{t}))) \\ = E_{f, g, φ_{Y}} (x_{t}, y_{t}) . \end{aligned}

(61)

Comparison with Eq. shows that

Φ_{σ^{t} (\bar{ω})} (y_{t}) = Φ (y_{t})

⁠, as it should.

C. Synchronization map

Generalized synchronization can be extended from the noiseless dynamic

[f, g]

to the noisy dynamic

[f_{ω}, g]

⁠,

ω \in Ω^{Z}

⁠, where the dynamic changes at every time step, by requiring

y_{t} = h_{σ^{t} (ω)} (x_{t}) .

(62)

Definition 19

We say that the responder $Y$ is synchronized to a driver $X$ perturbed by the noise $ω \in Ω^{Z}$ ⁠, if there is a sequence of continuous maps $h_{σ^{t} (ω)} : M_{X} \to M_{Y}$ such that Eq. (62) holds for all $t \in Z$ ⁠.

More generally, the synchronization map of order

k \geq 0

⁠,

h^{(k)} = h ° f^{k}

[Eq. ], generalizes to

h_{σ^{t} (ω)}^{(0)} := h_{σ^{t} (ω)}

and

y_{t} = h_{σ^{t} (ω)}^{(k)} (x_{t - k}) := h_{σ^{t} (ω)} ° f_{ω_{t - 1}, \dots, ω_{t - k}} (x_{t - k}),

(63)

for

k \geq 1

in the noisy case, while the synchronization map of period

K \geq 1

generalizes to

y_{t} = \frac{1}{K} \sum_{k = 0}^{K - 1} h_{σ^{t} (ω)}^{(k)} (x_{t - k}) =: h_{K, σ^{t} (ω)} (x_{t}, \dots, x_{t - K + 1}) .

(64)

To define the synchronization map in the reconstructed spaces, we replace

y_{t} = Φ^{- 1} (x_{t})

with

y_{t} = H_{σ^{t} (ω)} (x_{t})

in the “noisy” version of diagram ,

\begin{array}{rcccl} M_{X} ∋ & x_{t} & \overset{h_{σ^{t} (ω)}}{⟶} & y_{t} & \in M_{Y} \\ E_{f, φ_{X}, σ^{t} (ω)}^{- 1} & ↑ & ↑ & E_{g, φ_{Y}, σ^{t} (ξ)}^{- 1} \\ N_{X, σ^{t} (ω)} ∋ & x_{t} & \overset{H_{σ^{t} (ω)}}{⟶} & y_{t} & \in N_{Y, σ^{t} (ξ)} \end{array} .

(65)

Here, we used Eqs. and and assume that the maps

E_{f, φ_{X}, σ^{t} (ω)}

and

E_{g, φ_{Y}, σ^{t} (ξ)}

are embeddings. Then, it follows from

\begin{aligned} y_{t} & = (E_{g, φ_{Y}, σ^{t} (ξ)} ° h_{σ^{t} (ω)} ° E_{f, φ_{X}, σ^{t} (ω)}^{- 1}) (x_{t}) \\ =: H_{σ^{t} (ω)} (x_{t}) . \end{aligned}

(66)

Furthermore, by Eq. ,

\begin{aligned} x_{t} & = F_{σ^{t - 1} (ω)} (x_{t - 1}) \\ = \dots = (F_{σ^{t - 1} (ω)} ° \dots ° F_{σ^{t - k} (ω)}) (x_{t - k}) = F_{σ^{t - 1} (ω)}^{(k)} (x_{t - k}), \end{aligned}

(67)

so that

\begin{aligned} y_{t} & = \frac{1}{K} \sum_{k = 0}^{K - 1} (H_{σ^{t} (ω)} ° F_{σ^{t - 1} (ω)}^{(k)}) (x_{t - k}) \\ =: H_{K, σ^{t} (ω)} (x_{t}, \dots, x_{t - K + 1}) \end{aligned}

(68)

generalizes the reconstructed synchronization map of period

K

⁠, Eq. , to the noisy case.

To check that

y_{t} = H_{σ^{t} (ω)} (x_{t})

becomes

y_{t} = Φ^{- 1} (x_{t})

when the noise is switched off in Eq. , replace

E_{g, φ_{Y}, σ^{t} (ξ)}^{- 1} : y_{t} ⟶ y_{t}

(69)

on the right column of diagram with

Π_{Y} ° E_{g, φ_{Y}, σ^{t} (ω)}^{- 1} : y_{t} \overset{E_{g, φ_{Y}, σ^{t} (ω)}^{- 1}}{⟶} (x_{t}, y_{t}) \overset{Π_{Y}}{⟶} y_{t}

(70)

so that, according to Eq. ,

Π_{Y} ° E_{g, φ_{Y}, σ^{t} (ω)}^{- 1} (y_{t})

becomes

Π_{Y} ° E_{f, g, φ_{Y}}^{- 1} (y_{t})

in the noiseless case

ω = \bar{ω}

⁠, i.e.,

{\bar{ω}}_{t} = ω_{0}

for all

t \in Z

⁠. Finally, set

h (x_{t}) = h_{σ^{t} (\bar{ω})}

to convert diagram to diagram , thus identifying

H_{σ^{t} (\bar{ω})} (x_{t})

with

Φ^{- 1} (x_{t})

⁠, as it should be.

The numerical simulations of Sec. VI show that synchronization is robust against dynamical noise for strong enough couplings and, hence, can occur in the presence of dynamical noise. On the other hand, if synchronization occurs in the presence of dynamical noise but disappears when noise is switched off, then one speaks of noise-induced synchronization.⁵⁷

Finally, we can generalize the concepts of asymptotic synchronization and stability of the responder in the presence of noise in the driver as follows. We say that

Y

is asymptotically synchronized to the noisy driver

X

if the definition of synchronization, Eq. , holds only asymptotically, i.e.,

lim_{t \to \infty} ‖ y_{t} - H_{σ^{t} (ω)} (x_{t}) ‖ = 0,

(71)

where

‖ \cdot ‖

is a norm in

R^{\dim_{Y}}

⁠. It follows then that

Y

is asymptotically stable, i.e., the orbits of

Y

converge to

H_{σ^{t} (ω)} (x_{t})

⁠, regardless of their initial conditions. Asymptotic stability can easily be checked in practice. As in the noiseless case, it is a handy method to rule out synchronization.

VI. NUMERICAL SIMULATIONS

Unlike identical synchronization, which can be easily visualized, generalized synchronization is more difficult to detect. As mentioned in the Introduction, there exists an extensive literature on methods to detect functional dependency (and generalized synchronization for that matter) between two time series. The functional dependency targeted in this section is the synchronization map of a certain period $K > 1$ given in Eqs. (32) and (68) for noiseless and noisy drivers, respectively. For this reason, we use recurrent neural networks of the type long short-term memory (LSTM), which excel at predicting data from time series and are robust to noise. In fact, the LSTM nets outperformed the perceptrons (⁠ $K = 1$ ⁠) in the numerical simulations below, so we will only report the results obtained with the former. As a benchmark we use nearest-neighbor cross prediction (Sec. II B) because it is based on the continuity of the cross map (and its inverse in case of synchronization). In addition, nearest-neighbor cross prediction is robust against noise, particularly if the neighborhoods are well populated.

A. Models

For the numerical simulations we chose two unidirectionally coupled Hénon maps with several structural parameters and varying coupling strength. This testbed, first proposed by Schiff et al.¹⁷ and studied with the normalized mutual error, has been revisited several times in the literature, e.g., in Quian Quiroga et al.,²⁴ where the authors use the conditional Lyapunov exponent and the so-called nonlinear interdependencies.²³

Thus, the equations of the driver

X

⁠, with states

x = (x^{(1)}, x^{(2)})

in a trapping region of the attractor, are

{\begin{cases} x_{t + 1}^{(1)} = 1.4 - {(x_{t}^{(1)})}^{2} + (b_{1} + ω_{t}) x_{t}^{(2)}, \\ x_{t + 1}^{(2)} = x_{t}^{(1)}, \end{cases}

(72)

where

b_{1}

is a constant and

ω_{t}

are i.i.d. random numbers in the interval

[- A, A]

⁠, the noiseless scenario corresponding to

A = 0

⁠. The observation function is

φ_{X} (x_{t}) = x_{t}^{(1)}

⁠, i.e., the projection on the first component.

The equations of the responder

Y

⁠, with states

y = (y^{(1)}, y^{(2)})

⁠, are

{\begin{cases} y_{t + 1}^{(1)} = 1.4 - [C x_{t}^{(1)} y_{t}^{(1)} + (1 - C) {(y_{t}^{(1)})}^{2}] + b_{2} y_{t}^{(2)}, \\ y_{t + 1}^{(2)} = y_{t}^{(1)}, \end{cases}

(73)

where

b_{2}

is a constant and

C

is the coupling strength. For

C = 0

⁠, systems

X

and

Y

are uncoupled. The observation function is again the projection on the first component,

φ_{Y} (y_{t}) = y_{t}^{(1)}

⁠.

The parameter settings are as follows.

The settings for the constants $b_{1}$ and $b_{2}$ are the same as in Schiff et al.¹⁷ and Quian Quiroga et al.²⁴ So, we first set $b_{1} = b_{2} = 0.3$ ⁠, the standard values of the Hénon map, to study the coupling of identical systems (Model Hénon 0.3–0.3), which allows identical synchronization (i.e., $y_{t} = x_{t}$ ⁠) for $C = 1$ ⁠. Then, to study the coupling of non-identical systems, we set $b_{1} = 0.3$ ⁠, $b_{2} = 0.1$ (Model Hénon 0.3–0.1) and $b_{1} = 0.1$ ⁠, $b_{2} = 0.3$ (Model Hénon 0.1–0.3).
For the previous choices of $b_{1}$ and $b_{2}$ ⁠, we found that the driver orbits can diverge for noise amplitudes $A > 0.013$ ⁠, so we restrict them to the interval $0 \leq A \leq 0.013$ ⁠. The amplitudes used in the figures below are $A = 0$ (noiseless driver), $0.005$ and $0.013$ ⁠.
The range of the coupling strength $C$ is $0 \leq C \leq 1.2$ ⁠; the increment of $C$ in the figures below is $Δ C = 0.05$ ⁠.
For each case described above (identical/non-identical systems, noiseless/noisy driver), one series ${(x_{t})}_{0 \leq t \leq T - 1}$ and one series ${(y_{t})}_{0 \leq t \leq T - 1}$ were generated with seeds $x_{0} = (0, 0.9)$ and $y_{0} = (0.75, 0)$ ⁠, and length $T = 10^{5}$ (after discarding the first $1000$ points). Since we are only interested in synchronization, one series per case suffices because of asymptotic stability (Sec. II C).
The embedding dimension in the noiseless and noisy scenarios is $d = 5$ ⁠, i.e., $x_{t} = (x_{t}^{(1)}, \dots, x_{t + 4}^{(1)})$ and $y_{t} = (y_{t}^{(1)}, \dots, y_{t + 4}^{(1)})$ ⁠, $0 \leq t \leq T - 5$ ⁠. A posteriori justification for this choice are the excellent results obtained in the benchmark below.

The methods to test for synchronization in the noiseless case (⁠ $A = 0$ ⁠) and noisy cases (⁠ $A = 0.005$ ⁠, $0.013$ ⁠) are the following.

Our first method unveils synchronization by detecting functional dependencies, namely, the existence of the synchronization map of period $K = 10$ for time-delay vectors, i.e.,

$y_{t} = H_{10} (x_{t}, x_{t - 1}, \dots, x_{t - 9})$

(74)
[see Eq. (32)]. To do this, we used a three-layer neural network to predict $y_{t}$ based on $x_{t}, \dots, x_{t - 9}$ ⁠. Specifically, (i) the input layer consisted of an LSTM net with $5$ units, hidden states of dimension $10$ (corresponding to the inputs $x_{t}, \dots, x_{t - 9}$ ⁠) and the activation function $ReLU (x) = max {0, x}$ ⁠; (ii) the intermediate layer had 25 neurons and the activation function $Sigmoid (x) = 1 / (1 + e^{- x})$ ⁠; and (iii) the output layer had 5 neurons. Hence, the output layer returns $5$ states, corresponding to the 5 components of ${\hat{y}}_{t}$ ⁠, the prediction of $y_{t}$ ⁠. The network was trained with an $80 %$ of the data (the first $80 000$ time-delay vectors) and stochastic gradient descend, while the remaining $20 %$ of the data was used for testing. The accuracy of the predictions ${\hat{y}}_{t}$ output by the neural network based on the testing data $x_{t}$ ⁠,…, $x_{t - 9}$ (i.e., for each $t = 80 000$ ⁠,…, $99 990$ ⁠) was measured by their mean squared error (MSE), mean absolute error (MAE), and mean absolute scaled error (MASE). The matching of these three metrics in both the training and testing phases discards overfitting. Furthermore, since predictions based on data patterns are robust against low levels of noise, we expect this method to work well in both the noiseless and noisy cases.
As a benchmark we used nearest-neighbor cross prediction which, in the noiseless case, estimates $x_{t}$ based on the continuity of the cross map $x_{t} = Φ (y_{t})$ and corresponding nearest neighbors.¹⁸ Following the convergent cross mapping (CCM) method, we measured the accuracy of those estimations by $r (x, \hat{x})$ ⁠, the Pearson correlation coefficient of the estimates ${\hat{x}}_{t}$ obtained with the $d + 1 = 6$ nearest neighbors of $y_{t}$ ⁠. Since $T = 10^{5}$ ⁠, the time series are sufficiently long to obtain good estimates (actually only the first $10 000$ points were used), so $r (x, \hat{x}) ≃ 1$ ⁠. On the contrary, if $r (y, \hat{y})$ is the Pearson correlation coefficient of the estimates ${\hat{y}}_{t}$ obtained via the $d + 1 = 6$ nearest-neighbors of $x_{t}$ ⁠, then we expect $r (y, \hat{y}) ≃ 0$ ⁠, unless $Y$ synchronizes with $X$ ⁠, in which case $y = Φ^{- 1} (x)$ and $r (y, \hat{y}) ≃ 1$ (this time due to the continuity of $Φ^{- 1}$ ⁠). We conclude that if

$Δ r = r (x, \hat{x}) - r (y, \hat{y})$

(75)
and $X ⇝ Y$ ⁠, then
- $0 \leq$ $Δ r \leq 1$ ⁠, and
- $Δ r ≃ 0$ signalizes synchronization, except when $r (x, \hat{x}) = 0 = r (y, \hat{y})$ ⁠, i.e., when $X$ and $Y$ are uncoupled.
In the noisy cases, the situation is qualitatively the same, thanks to the robustness of nearest-neighbor cross prediction against low levels of noise. See, e.g., Sugihara et al.,³⁵ Mønster et al.⁵⁸ and the book by Datseris and Parlitz⁵⁹ for CCM algorithms to compute (75).

B. Results

Out of the accuracy results obtained with the LSTM network and testing data, we are going to discuss only the MSE vs $C$ curves since the other two curves, MAE and MASE vs $C$ ⁠, are similar for all models.

The results of the numerical simulations are depicted in Figs. 1–3 for Method 1 [panels (a)] and Method 2 [panels (b)] and the three models Hénon 0.3–0.3 (Fig. 1), 0.3–0.1 (Fig. 2), and 0.1–0.3 (Fig. 3). Comparison of both panels for each case and $C > 0$ shows an excellent agreement of both methods on the synchronization states, i.e., MSE $(C) = 0$ in panels (a) and $Δ r (C) = 0$ in panels (b). As noted above, $Δ r (0) = 0$ in all cases owing to the fact that $X$ and $Y$ are uncoupled for $C = 0$ ⁠; such numerical artifacts can be easily filtered out by checking whether $r (x, \hat{x}) ≃ 0$ and $r (y, \hat{y}) ≃ 0.$

FIG. 1.

View large Download slide

Numerical results for the model Hénon 0.3–0.3, i.e., $b_{1} = 0.3$ in (72) and $b_{2} = 0.3$ in (73). (a) MSE vs the coupling strength $C$ for a noiseless driver (noise amplitude $A = 0$ ⁠) and a noisy driver (⁠ $A = 0.005$ ⁠, $0.013$ ⁠) obtained with an LSTM net. (b) $Δ r$ vs $C$ for a noiseless driver (⁠ $A = 0$ ⁠) and a noisy driver (⁠ $A = 0.005$ ⁠, $0.013$ ⁠) obtained via $6$ -nearest-neighbor cross prediction. See text for more detail.

FIG. 2.

View large Download slide

Numerical results for the model Hénon 0.3–0.1, i.e., $b_{1} = 0.3$ in (72) and $b_{2} = 0.1$ in (73). The information displayed in the panels (a) and (b) is the same as in Fig. 1.

FIG. 3.

View large Download slide

Numerical results for the model Hénon 0.1–0.3, i.e., $b_{1} = 0.1$ in (72) and $b_{2} = 0.3$ in (73). The information displayed in the panels (a) and (b) is the same as in Fig. 1.

The main conclusions from the numerical results can be summarized as follows.

Small-amplitude noise does not destroy all the states of “strong” synchronization (i.e., due to strong enough couplings), it only shifts the synchronization threshold to higher values. So, synchronization can also occur in the presence of dynamical noise.
Synchronization due to strong enough couplings is robust against small-amplitude dynamical noise, while synchronization states with a weak coupling strength can be unstable whatever the amplitude of the noise. This fact is illustrated in the Model Hénon 0.1–0.3 (Fig. 3), where synchronization is detected in the interval $0.5 ≲ C ≲ 0.6$ for $A = 0$ ⁠.
Weakly coupled systems can be asymptotically synchronized, which can be detected via the auxiliary systems method both in the noiseless and noisy cases. Indeed, Table I shows the intervals of coupling strengths for which the responder is asymptotically stable, obtained with the auxiliary system method. Therefore, synchronization can occur only for couplings in the corresponding interval (as it does). Note that Table I excludes the spurious synchronization $Δ r (0) = 0$ ⁠.
In general, when the noise amplitude increases, the threshold of stable synchronization moves towards stronger couplings. However, the Model Hénon 0.3–0.1 (Fig. 2) shows that there can be parameter settings for which that threshold is virtually the same for the noise amplitudes considered here.

TABLE I.

Coupling strengths in the range 0 ≤ C ≤ 1.2 for which the responder is asymptotically stable.

	Hénon 0.3–0.3	Hénon 0.3–0.1	Hénon 0.1–0.3
A = 0	0.40 ≤ C ≤ 1.20	0.15 ≤ C ≤ 1.20	0.40 ≤ C ≤ 1.20
A = 0.005	0.40 ≤ C ≤ 1.20	0.20 ≤ C ≤ 1.20	0.40 ≤ C ≤ 1.20
A = 0.013	0.40 ≤ C ≤ 1.20	0.20 ≤ C ≤ 1.20	0.40 ≤ C ≤ 1.20

	Hénon 0.3–0.3	Hénon 0.3–0.1	Hénon 0.1–0.3
A = 0	0.40 ≤ C ≤ 1.20	0.15 ≤ C ≤ 1.20	0.40 ≤ C ≤ 1.20
A = 0.005	0.40 ≤ C ≤ 1.20	0.20 ≤ C ≤ 1.20	0.40 ≤ C ≤ 1.20
A = 0.013	0.40 ≤ C ≤ 1.20	0.20 ≤ C ≤ 1.20	0.40 ≤ C ≤ 1.20

Finally, let us point out that we also performed numerical simulations with perceptrons to detect the possible existence of the conventional synchronization map (period 1). The results were similar but not as sharp regarding weak synchronization states as the results obtained with LSTM nets to detect synchronization maps of period $K > 1$ ⁠. However, no performance analysis of Method 1 with respect to the period $K$ was carried out and, hence, no attempt was made to optimize the parameter $K$ (which, anyway, depends on the data at hand).

VII. APPLICATION TO REAL-WORLD DATA: EEGs

The purpose of this section is to illustrate the application of Method 1 to real world data, specifically, intracranial EEG recordings from a subject with epilepsy. Therefore, we will not scrutinize here the complexity of such signals but rather check whether our findings align with results obtained in previous studies.

First of all, we notice that real observations ${(φ_{X} (x_{t}))}_{1 \leq t \leq T}$ and ${(φ_{Y} (y_{t}))}_{1 \leq t \leq T}$ of coupled systems $X$ and $Y$ ⁠, respectively, can deviate from our assumptions in Secs. II–VI in two important issues: nonstationarity or/and bidirectionality of the coupling, as it actually occurs with the time series in this section. To meet these challenges, this time we will apply Method 1 (Sec. VI) in both directions $X ⇝ Y$ and $Y ⇝ X$ ⁠, on sufficiently short data segments to ensure approximate stationarity.

There is a subtlety, though. In the unidirectionally coupling $X ⇝ Y$ studied in the previous sections, we detected synchronization by detecting a functional dependency between $y_{t}$ and $x_{t}, \dots, x_{t - K + 1}$ ⁠, namely, $y_{t} = H_{K, X ⇝ Y} (x_{t}, \dots, x_{t - K + 1})$ ⁠, where $H_{K, X ⇝ Y}$ is the reconstructed synchronization map of period $K$ of the coupling $X ⇝ Y$ ⁠, defined in Eq. (32). The robustness to noise of $H_{K, X ⇝ Y}$ allowed us then to extend our conclusions to signals contaminated with low-amplitude noise. If, for the sake of this argument, we think of a bidirectional coupling $X ↭ Y$ as the joint action of two separate unidirectional couplings $X ⇝ Y$ and $Y ⇝ X$ ⁠, then $y_{t}$ will depend on $x_{t}, \dots, x_{t - K + 1}$ (whether $X$ and $Y$ are synchronized or not) through the cross map of period $K$ of the coupling $Y ⇝ X$ ⁠, i.e., $y_{t} = Φ_{K, Y ⇝ X} (x_{t}, \dots, x_{t - K + 1})$ ⁠; see Eq. (17) with $y$ and $x$ swapped. Therefore, here we expect $y_{t}$ to depend on $x_{t}, \dots, x_{t - K + 1}$ in general.

The bottom line is that, by using Method 1 in the directions $X ⇝ Y$ and $Y ⇝ X$ ⁠, we will be able to detect the “dominant driver” or the “coupling directionality” of the bidirectional coupling $X ↭ Y$ ⁠. To this end, we are going to measure the strength of the coupling in both directions via the accuracy of the predictions of $x_{t}$ and $y_{t}$ made by LSTM nets in short non-overlapping segments over the entire EEGs, the dominant driver being given by the direction with the strongest coupling. In case of equal strengths, the systems are assumed to be synchronized.

A. Data description

The data that we are going to analyze is the following; see Lehnertz and Dickten⁶⁰ for more detail.

The signals are EEGs recorded intracranially from a subject with epilepsy during $86 090$ s (23 h, 54 m, 50.4 s) with 48 electrode contacts at a sampling frequency of 200 Hz (sampling time\break = 5 ms). The subject had signed informed consent that her/his clinical data might be used and published for research purposes, and the study protocol had previously been approved by the ethics committee of the University of Bonn. The recording started at 7:00 am, corresponding to the initial sampling interval $t = 1$ ⁠, and ended at the final sampling time $t_{final} = 17.217 984 \times 10^{6}$ ⁠. The epileptic convulsions occurred at the following sampling times.
- Average time of a first group of subclinical seizures: ${\bar{t}}_{C 1} = 3. 4082 \times 10^{6}$ (⁠ $17 041$ s). By subclinical seizures we mean localized seizure activity on the EEG with no obvious clinical activity.
- Average time of a second group of subclinical seizures: ${\bar{t}}_{C 2} = 4. 6082 \times 10^{6}$ (⁠ $23, 041$ s).
- Average time of a third group of subclinical seizures: ${\bar{t}}_{C 3} = 6. 8762 \times 10^{6}$ (⁠ $34 381$ s).
- Onset time of a clinical seizure (the only one in the whole series): $t_{C 4} = 17.1842 \times 10^{6}$ (⁠ $85 921$ s).
A schematic of the implanted electrodes can be found in Fig. 2 of Lehnertz and Dickten.⁶⁰ The electrodes contacts are divided into the following three categories:
- focal (F), which comprises electrode contacts located within the seizure-onset zone;
- neighbor (N), which groups electrode contacts not more than two contacts distant to those of category $F$ ⁠;
- other (O), gathering all remaining electrode contacts.
To designate the electrode contacts and their (approximately) 24 h recordings, we use the same labels as in Ref. 60. For example, $X = TR 01$ means that the system $X$ is the source of the EEG ${(φ_{X} (x_{t}))}_{1 \leq t \leq t_{final}}$ recorded at the electrode contact TR01.
For the sake of our analysis, we will consider the following five pairs $(X, Y)$ of electrode contacts.
- Case 1: $(X, Y) = (TR 05$ – $TR 06)$ in the categories (F–F).
- Case 2: $(X, Y) = (TR 07$ – $TBPR 1)$ in the categories (F–N).
- Case 3: $(X, Y) = (TR 05$ – $TL 05)$ in the categories (F–O).
- Case 4: $(X, Y) = (TBAR 1$ – $TLL 04)$ in the categories (N–O).
- Case 5: $(X, Y) = (TLR 04$ – $TLL 04)$ in the categories (O–O).
As in the previous numerical simulations, the embedding dimension of the systems $X$ and $Y$ is 5. Thus,

$x_{t} = (φ_{X} (x_{t}), φ_{X} (x_{t + 1}), \dots, φ_{X} (x_{t + 4}))$

(76)
$1 \leq t \leq t_{final} - 4$ ⁠, are the time-delay vectors corresponding to system $X$ ⁠, and analogously with the EEG ${(φ_{Y} (y_{t}))}_{1 \leq t \leq t_{final}}$ generated by the system $Y$ ⁠.
For approximate stationarity,^64,65 we partitioned the time series $(x_{t})$ and $(y_{t})$ ⁠, $1 \leq t \leq t_{final} - 4$ ⁠, into 1434 non-overlapping segments

$S_{X, n} = {(x_{t})}_{12 000 (n - 1) + 1 \leq t \leq 12 000 n - 4}$

(77)
and

$S_{Y, n} = {(y_{t})}_{12 000 (n - 1) + 1 \leq t \leq 12 000 n - 4}$

(78)
of 11 996 points (⁠ $≃ 60$ s) each, $n = 1, 2, \dots, 1434$ ⁠, and a last pair of segments

$S_{X, 1435} = {(x_{t})}_{1 72 08 001 \leq t \leq 1 72 17 980}$

(79)
and

$S_{Y, 1435} = {(y_{t})}_{1 72 08 001 \leq t \leq 1 72 17 980},$

(80)
comprising only 9980 points (⁠ $≃ 50$ s). The segments $1 \leq n \leq 720$ ⁠, correspond to the daylight hours (7 am–7 pm), while the segments $721 \leq n \leq 1435$ correspond to the night hours. The clinical seizure occurs in the segment $n = 1433$ ⁠, i.e., the third to last segment of the series, and it initiates just one second after the beginning of that segment (⁠ $t_{C 4} = 85 921$ s).
As in Sec. VI, we use the first $80 %$ of the data of each $n$ th segment $S_{X, n}$ and $S_{Y, n}$ as training data, and the remaining $20 %$ as testing data. So, this time we obtain two accuracy measures: (i) ${MSE}_{X ⇝ Y} (n)$ ⁠, for the predictions of $y_{t}$ output by the LSTM net, based on $x_{t}, \dots, x_{t - K + 1}$ with testing data of the segments $S_{Y, n}$ and $S_{X, n}$ ⁠, and (ii) ${MSE}_{Y ⇝ X} (n)$ ⁠, for the predictions of $x_{t}$ output by the LSTM net, based on $y_{t}, \dots, y_{t - K + 1}$ with testing data of the segments $S_{X, n}$ and $S_{Y, n}$ ⁠. As in Sec. VI, we set $K = 10$ here. Of course, the parameter $K$ can be fine-tuned for optimal results, but this is an issue not contemplated in the present work.

B. Results

Since, at variance with the numerical simulations in Sec. VI, we have here bidirectionally coupled signals and two prediction accuracy measures

{MSE}_{X ⇝ Y} (n)

and

{MSE}_{Y ⇝ X} (n)

⁠, we are going to use the coupling directionality index

Δ MSE (n) = \frac{{MSE}_{X ⇝ Y} (n) - {MSE}_{Y ⇝ X} (n)}{{MSE}_{X ⇝ Y} (n) + {MSE}_{Y ⇝ X} (n)},

(81)

for each pair of data segments

S_{X, n}

and

S_{Y, n}

⁠,

1 \leq n \leq 1435

⁠, so that

$- 1 \leq Δ MSE (n) \leq + 1$ and
$Δ MSE (n) \geq 0$ if and only if ${MSE}_{Y ⇝ X} (n) \leq {MSE}_{X ⇝ Y} (n)$ ⁠, i.e., knowledge of $x_{t}$ in the $n$ th segment leads to better predictions of $y_{t}$ than the other way around.

Therefore, if $Δ MSE (n) > 0$ (respectively, $Δ MSE (n) < 0)$ ⁠, then we conclude that $X$ is the dominant driver (respectively, $Y$ is the dominant driver). This interpretation agrees with other approaches based on the cross map,^34,35 transfer entropy,^66,67 phase dynamics,⁶⁸ etc. Otherwise, if $Δ MSE (n) = 0$ ⁠, then $X$ and $Y$ are assumed to be synchronized in segment $n$ (although it might be difficult to discern this situation from “no-coupling”).

Figure 4 plots $Δ MSE (n)$ vs the segment number $n$ ⁠, $1 \leq n \leq 1435$ ⁠. The numerical results are summarized in Table II for the 24 h EEGs, and in Table III for 12 h EEGs corresponding to daylight hours (⁠ $1 \leq n \leq 720$ ⁠) and night hours (⁠ $721 \leq n \leq 1435$ ⁠). It was not possible to highlight the clinical seizure in Fig. 4 because it occurs in the segment $n = 1433$ ⁠, so any visual marks at that point are indistinguishable from the right margin of the corresponding panel.

FIG. 4.

Top to bottom: plots of Δ MSE ( n ), the directionality indicator (81), obtained using the segments S X , n and S Y , n, 1 ≤ n ≤ 1435, given in Eqs. (77)–(79), for cases 1–5. The clinical seizure occurs in the segment n = 1433, too close to the right margin to be marked. See Sec. VII A for detail.

View large Download slide

Top to bottom: plots of $Δ MSE (n)$ ⁠, the directionality indicator (81), obtained using the segments $S_{X, n}$ and $S_{Y, n}$ ⁠, $1 \leq n \leq 1435$ ⁠, given in Eqs. (77)–(79), for cases 1–5. The clinical seizure occurs in the segment $n = 1433$ ⁠, too close to the right margin to be marked. See Sec. VII A for detail.

TABLE II.

Results of cases 1–5 with 24 h EEGs.

Case	ΔMSE(n) > 0	Dominant electrode
1 (F–F)	$16 %$ of segments	Y = TR06 (F) dominates X = TR05 (F)
2 (F–N)	$64 %$ of segments	X = TR07 (F) dominates Y = TBPR1 (N)
3 (F–O)	$56 %$ of segments	X = TR05 (F) dominates Y = TL05 (O)
4 (N–O)	$5 %$ of segments	Y = TLL04 (O) dominates X = TBAR1 (N)
5 (O–O)	$7 %$ of segments	Y = TLL04 (O) dominates X = TLR04 (O)

Case	ΔMSE(n) > 0	Dominant electrode
1 (F–F)	$16 %$ of segments	Y = TR06 (F) dominates X = TR05 (F)
2 (F–N)	$64 %$ of segments	X = TR07 (F) dominates Y = TBPR1 (N)
3 (F–O)	$56 %$ of segments	X = TR05 (F) dominates Y = TL05 (O)
4 (N–O)	$5 %$ of segments	Y = TLL04 (O) dominates X = TBAR1 (N)
5 (O–O)	$7 %$ of segments	Y = TLL04 (O) dominates X = TLR04 (O)

In view of Fig. 4 and Tables II and III, we can draw the following general conclusions.

The coupling directionality, as measured by $Δ MSE (n)$ ⁠, depends on the segment $n$ ⁠. The overall dominance of the signals is stable with respect to day and night, although the dominance degrees, as measured by the percentages of segments contributing to the dominant coupling direction in the first or second 12 h, respectively, are different in all cases.
Except in case 2, the sign of $Δ MSE (n)$ during the epileptic convulsions (⁠ $n = 1433$ ⁠) coincides with the overall coupling direction. In fact,

Case 1 2 3 4 5

ΔMSE(1433) −0.58 −0.87 0.60 −0.68 −0.40

Case 1 2 3 4 5

ΔMSE(1433) −0.58 −0.87 0.60 −0.68 −0.40

View Large
According to Osterhage et al.,^62,63 an important question in epileptology is whether the pathological interaction between the seizure-onset zone (label F) and other brain areas (labels N and O), a phenomenon called focal driving, can also be identified during seizure-free periods. Cases 2 and 3 in Table II answer this question affirmatively, that is, our method detects focal driving in the analyzed EEG.
In addition, cases 2 and 3 in Table III indicate that focal driving is not diminished during sleep.
More generally, Table III shows that focal driving does not appear to be influenced by other (possibly “stronger”) synchronization phenomena such as sleep.
The dominance degree is rather high in cases 1, 4, and 5, with $Δ MSE (n) < 0$ over 80% of the segments both in the 24 h and 12 h EEGs. Note that the interaction in those cases is local (cases 1 and 5) or it does not involve the seizure generating area (case 4).

The above findings are in line with the results of more comprehensive studies by Lehnertz and Dickten,⁶⁰ Dickten et al.,⁶¹ and Osterhage et al.,^62,63 which empirically demonstrates the capability of our LSTM net-based method.

TABLE III.

Results of cases 1–5 with 12 h EEGs (day and night).

Case	ΔMSE(n) > 0 day	ΔMSE(n) > 0 night
1 (F–F)	$18 %$ (Y dominant)	$15 %$ (Y dominant)
2 (F–N)	$70 %$ (X dominant)	$58 %$ (X dominant)
3 (F–O)	$51 %$ (X dominant)	$60 %$ (X dominant)
4 (N–O)	$6 %$ (Y dominant)	$4 %$ (Y dominant)
5 (O–O)	$9 %$ (Y dominant)	$6 %$ (Y dominant)

Case	ΔMSE(n) > 0 day	ΔMSE(n) > 0 night
1 (F–F)	$18 %$ (Y dominant)	$15 %$ (Y dominant)
2 (F–N)	$70 %$ (X dominant)	$58 %$ (X dominant)
3 (F–O)	$51 %$ (X dominant)	$60 %$ (X dominant)
4 (N–O)	$6 %$ (Y dominant)	$4 %$ (Y dominant)
5 (O–O)	$9 %$ (Y dominant)	$6 %$ (Y dominant)

VIII. CONCLUSION

Synchronization of two unidirectionally coupled dynamical systems $X ⇝ Y$ is a classical topic in nonlinear dynamics. It is defined by the existence of a continuous function $y_{t} = h (x_{t})$ between the states $x_{t}$ of the driver $X$ and the states $y_{t}$ of the responder, called the synchronization map. While $h$ ⁠, when it exists, points from the state space of the driver (domain) to the state space of the responder (range), the cross map $x_{t} = Φ (y_{t})$ always exists in that framework, is continuous, and points in the opposite direction between the corresponding reconstructed state spaces. In the standard, noiseless scenario, the existence of the synchronization map (i.e., synchronization between $X$ and $Y$ ⁠) amounts to $Φ$ being invertible and bicontinuous. These and other fundamentals of generalized synchronization in the absence of dynamical noise were presented in a self-contained and unified way in Sec. II, with emphasis on the relationship between the cross map and the synchronization map.

In this context, the main contributions of the present paper are the following.

Introduction of higher-period versions of the cross and synchronization maps in Eqs. (17) and (21), the period-1 versions corresponding to the conventional concepts. They are based on the corresponding maps of order $k$ ⁠, defined in Eqs. (14) and (20), and, actually, they may be defined in many different ways. A synchronization map of period 10 was used in the numerical simulation of Sec. VI because it gave better results than the conventional map in the detection of synchronization. Higher-period cross maps were invoked in Sec. VII to understand the sign of the directionality index (81) when the coupling is bidirectional. Optimization of the period was not discussed because it is beyond the scope of this paper.
Generalizations of the synchronization map and the cross map when the driver is noisy in Secs. V B and V C, respectively. To this end, the dynamical noise was modeled as stochastic forcing. The generalizations consist of families of maps that depend on noise parameters and coincide with their conventional counterparts when the noise is switched off. As usual, those generalizations have the wished properties under some formal provisos, e.g., generacy of the maps and parameters involved, as well as the driver states in the case of the cross map. However, this does not mean that they are only useful in theory; they can be also useful in practice, e.g., in laboratory or numerical experiments, where typical properties are taken for granted and even noise parameters may be known.
Application of LSTM nets to detect synchronization in synthetic data. This method harnesses the existence of synchronization maps of higher periods in both the noiseless and noisy scenarios. To be more precise, in the numerical simulations of Sec. VI, synchronization was revealed in the reconstructed state spaces by detecting a period-10 synchronization map [Eq. (32)] using an LSTM net and predictability. The dynamical systems were two coupled Hénon maps with several parameter settings and noise amplitudes. As a benchmark, we used nearest-neighbor cross prediction based on the continuity of the cross map and its inverse (in case of synchronization). The agreement of the results obtained with the two methods was excellent. The results also showed the robustness of both methods against noise.
Application of LSTM nets to detect the dominant driver in real-world data in Sec. VII. The real-world data consisted of a 24 h EEG from a subject with epilepsy. To cope with the nonstationarity of the signal and the bidirectionality of the couplings between different brain areas, we partitioned the EEG in non-overlapping 60 s segments and measured the coupling dominance by the directionality index $Δ MSE (n)$ defined in Eq. (81); this index is based on the mean square error of predictions made by LSTM nets in the $n$ th segment. The results agreed with results published in the literature, in particular, the existence of focal driving and its robustness to the day/night cycle.

ACKNOWLEDGMENTS

J. M. Amigó and R. Dale were financially supported by Agencia Estatal de Investigación, Spain, Grant No. PID2019-108654GB-I00/AEI/10.13039/501100011033. J. M. Amigó was also supported by Generalitat Valenciana, Spain, Grant No. PROMETEO/2021/063.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Ethics Approval

The studies involving human participants were reviewed and approved by the ethics committee of the University of Bonn. The subject provided written informed consent to participate in this study.

Ethics approval for experiments reported in the submitted manuscript on animal or human subjects was granted. The studies involving human participants were reviewed and approved by the ethics committee of the University of Bonn. The subject provided written informed consent to participate in this study.

Author Contributions

José M. Amigó: Conceptualization (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Validation (equal); Writing – original draft (lead); Writing – review & editing (equal). Roberto Dale: Data curation (equal); Software (lead); Validation (equal); Visualization (lead); Writing – review & editing (equal). Juan C. King: Data curation (equal); Software (lead); Validation (equal); Visualization (equal); Writing – review & editing (equal). Klaus Lehnertz: Data curation (lead); Methodology (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal).

DATA AVAILABILITY

The EEG dataset presented in this article is not publicly available due [state restrictions such as privacy or ethical restrictions]. Requests to access the dataset should be directed to [email protected].

REFERENCES

1.

V. S.

Afraimovich

,

N. N.

Verichev

, and

M. I.

Rabinovich

, “

Stochastic synchronization of oscillation in dissipative systems

,”

Radiophys. Quantum Electron.

29

,

795

–

803

(

1986

).

https://doi.org/10.1007/BF01034476

Google Scholar

Crossref

2.

N. F.

Rulkov

,

M. M.

Sushchik

,

L. S.

Tsimring

, and

H. D. I.

Abarbanel

, “

Generalized synchronization of chaos in directionally coupled chaotic systems

,”

Phys. Rev. E

51

,

980

–

994

(

1995

).

https://doi.org/10.1103/PhysRevE.51.980

Google Scholar

Crossref

3.

M. G.

Rosenblum

,

A. S.

Pikovsky

, and

J.

Kurths

, “

From phase to lag synchronization in coupled chaotic oscillators

,”

Phys. Rev. Lett.

78

,

4193

–

4196

(

1997

).

https://doi.org/10.1103/PhysRevLett.78.4193

Google Scholar

Crossref

4.

L. M.

Pecora

and

T. L.

Carroll

, “

Synchronization in chaotic systems

,”

Phys. Rev. Lett.

64

,

821

–

824

(

1990

).

https://doi.org/10.1103/PhysRevLett.64.821

Google Scholar

Crossref

PubMed

5.

S.

Boccaletti

,

J.

Kurths

,

G.

Osipov

,

D. L.

Valladare

, and

C. S.

Zhou

, “

The synchronization of chaotic systems

,”

Phys. Rep.

366

,

1

–

101

(

2002

).

https://doi.org/10.1016/S0370-1573(02)00137-0

Google Scholar

Crossref

6.

K. M.

Cuomo

,

A. V.

Oppenheim

, and

S. H.

Strogatz

, “

Synchronization of Lorenz-based chaotic circuits with applications to communications

,”

IEEE Trans. Circuits Syst. II

40

,

626

–

633

(

1993

).

https://doi.org/10.1109/82.246163

Google Scholar

Crossref

7.

L.

Kocarev

and

U.

Parlitz

, “

General approach for chaotic synchronization with applications to communication

,”

Phys. Rev. Lett.

74

,

5028

–

5031

(

1995

).

https://doi.org/10.1103/PhysRevLett.74.5028

Google Scholar

Crossref

PubMed

8.

T.

Nowotny

,

R.

Huerta

, and

M. I.

Rabinovich

, “

Neuronal synchrony: Peculiarity and generality

,”

Chaos

18

,

037119

(

2008

).

https://doi.org/10.1063/1.2949925

Google Scholar

Crossref

PubMed

9.

J. M.

Amigó

,

T. S.

Mosqueiro

, and

R.

Huerta

, “

Predicting synchronization of three mutually inhibiting groups of oscillators with strong resetting

,”

Appl. Math. Inf. Sci.

9

,

2245

–

2256

(

2015

).

https://doi.org/10.12785/amis/090505

Google Scholar

10.

P.

Hövel

,

A.

Viol

,

P.

Loske. L. Merfort

, and

V.

Vuksanović

, “

Synchronization in functional networks of the human brain

,”

J. Nonlinear Sci.

30

,

2259

–

2282

(

2020

).

https://doi.org/10.1007/s00332-018-9505-7

Google Scholar

Crossref

11.

L. J.

Kocarev

,

K. S.

Halle

,

K.

Eckert

,

L. O.

Chua

, and

U.

Parlitz

, “

Experimental demonstration of secure communications via chaotic synchronization

,”

Int. J. Bifurcation Chaos

2

,

709

–

713

(

1992

).

https://doi.org/10.1142/S0218127492000823

Google Scholar

Crossref

12.

J. M.

Amigó

,

L.

Kocarev

, and

J.

Szczepanski

, “

Discrete Lyapunov exponent and resistance to differential cryptanalysis

,”

IEEE Trans. Circuits Syst. II

54

,

882

–

886

(

2007

).

https://doi.org/10.1109/TCSII.2007.901576

Google Scholar

Crossref

13.

W.

Kinzel

,

A.

Englert

, and

I.

Kanter

, “

On chaos synchronization and secure communication

,”

Philos. Trans. R. Soc. A

368

,

379

–

389

(

2010

).

https://doi.org/10.1098/rsta.2009.0230

Google Scholar

Crossref

14.

Chaos Synchronization and Cryptography for Secure Communications: Applications for Encryption, edited by S. Banerjee (Information Science Reference, Hershey, PA, 2011).

15.

A.

Pikovsky

,

M.

Rosenblum

, and

J.

Kurths

,

Synchronization: A Universal Concept in Nonlinear Sciences

(

Cambridge University Press

,

Cambridge

,

2001

).

Google Scholar

Crossref

16.

L. M.

Pecora

and

T. L.

Carroll

, “

Synchronization of chaotic systems

,”

Chaos

25

,

097611

(

2015

).

https://doi.org/10.1063/1.4917383

Google Scholar

Crossref

PubMed

17.

S. J.

Schiff

,

P.

So

,

T.

Chang

,

R. E.

Burke

, and

T.

Sauer

, “

Detecting dynamical interdependence and generalized synchrony through mutual prediction in a neural ensemble

,”

Phys. Rev. E

54

,

6708

–

6724

(

1996

).

https://doi.org/10.1103/PhysRevE.54.6708

Google Scholar

Crossref

18.

M.

Le Van Quyen

,

J.

Martinerie

,

C.

Adam

, and

F. J.

Varela

, “

Nonlinear analyses of interictal EEG map the brain interdependences in human focal epilepsy

,”

Physica D

127

,

250

–

266

(

1999

).

https://doi.org/10.1016/S0167-2789(98)00258-9

Google Scholar

Crossref

19.

B. R.

Hunt

,

E.

Ott

, and

J. A.

Yorke

, “

Differentiable generalized synchronization of chaos

,”

Phys. Rev. E

55

,

4029

–

4034

(

1997

).

https://doi.org/10.1103/PhysRevE.55.4029

Google Scholar

Crossref

20.

K.

Pyragas

, “

Conditional Lyapunov exponents from time series

,”

Phys. Rev. E

56

,

5183

–

5188

(

1997

).

https://doi.org/10.1103/PhysRevE.56.5183

Google Scholar

Crossref

21.

H. D. I.

Abarbanel

,

N. F.

Rulkov

, and

M. M.

Sushchick

, “

Generalized synchronization of chaos: The auxiliary system approach

,”

Phys. Rev. E

53

,

4528

–

4535

(

1996

).

https://doi.org/10.1103/PhysRevE.53.4528

Google Scholar

Crossref

22.

L.

Kocarev

and

U.

Parlitz

, “

Generalized synchronization, predictability and equivalence of unidirectionally coupled dynamical systems

,”

Phys. Rev. Lett.

76

,

1816

–

1819

(

1996

).

https://doi.org/10.1103/PhysRevLett.76.1816

Google Scholar

Crossref

PubMed

23.

J.

Arnhold

,

P.

Grassberger

,

K.

Lehnertz

, and

C. E.

Elger

, “

A robust method for detecting interdependences: Application to intracranially recorded EEG

,”

Physica D

134

,

419

–

430

(

1999

).

https://doi.org/10.1016/S0167-2789(99)00140-2

Google Scholar

Crossref

24.

R.

Quian Quiroga

,

J.

Arnhold

, and

P.

Grassberger

, “

Learning driver-response relationships from synchronization patterns

,”

Phys. Rev. E

61

,

5142

–

5148

(

2000

).

https://doi.org/10.1103/PhysRevE.61.5142

Google Scholar

Crossref

25.

R.

Quian Quiroga

,

A.

Kraskov

,

T.

Kreuz

, and

P.

Grassberger

, “

Performance of different synchronization measures in real data: A case study on electroencephalographic signals

,”

Phys. Rev. E

65

,

041903

(

2002

).

https://doi.org/10.1103/PhysRevE.65.041903

Google Scholar

Crossref

26.

R.

Sowa

,

A.

Chernihovskyi

,

F.

Mormann

, and

K.

Lehnertz

, “

Estimating phase synchronization in dynamical systems using cellular nonlinear networks

,”

Phys. Rev. E

71

,

061926

(

2005

).

https://doi.org/10.1103/PhysRevE.71.061926

Google Scholar

Crossref

27.

D.

Krug

,

H.

Osterhage

,

C. E.

Elger

, and

K.

Lehnertz

, “

Estimating nonlinear interdependencies in dynamical systems using cellular nonlinear networks

,”

Phys. Rev. E

76

,

041916

(

2007

).

https://doi.org/10.1103/PhysRevE.76.041916

Google Scholar

Crossref

28.

R.

Monetti

,

W.

Bunk

,

T.

Aschenbrenner

, and

F.

Jamitzky

, “

Characterizing synchronization in time series using information measures extracted from symbolic representations

,”

Phys. Rev. E

79

,

046207

(

2009

).

https://doi.org/10.1103/PhysRevE.79.046207

Google Scholar

Crossref

29.

J. M.

Amigó

,

R.

Monetti

,

T.

Aschenbrenner

, and

W.

Bunk

, “

Transcripts: An algebraic approach to coupled time series

,”

Chaos

22

,

013105

(

2012

).

https://doi.org/10.1063/1.3673238

Google Scholar

Crossref

PubMed

30.

J. M.

Amigó

and

K.

Keller

, “

Permutation entropy: One concept, two approaches

,”

Eur. Phys. J. Spec. Top.

222

,

263

–

273

(

2013

).

https://doi.org/10.1140/epjst/e2013-01840-1

Google Scholar

Crossref

31.

D.

Ibáñez-Soria

,

J.

García-Ojalvo

,

A.

Soria-Frisch

, and

G.

Ruffini

, “

Detection of generalized synchronization using echo state networks

,”

Chaos

28

,

033118

(

2018

).

https://doi.org/10.1063/1.5010285

Google Scholar

Crossref

PubMed

32.

T.

Lymburn

,

D. M.

Walker

,

M.

Small

, and

T.

Jüngling

, “

The reservoir’s perspective on generalized synchronization

,”

Chaos

29

,

093133

(

2019

).

https://doi.org/10.1063/1.5120733

Google Scholar

Crossref

PubMed

33.

H.

Kantz

and

T.

Schreiber

,

Nonlinear Time Series Analysis

, 2nd ed. (

Cambridge University Press

,

Cambridge

,

2004

),

Google Scholar

34.

J. M.

Amigó

and

Y.

Hirata

, “

Detecting directional couplings from multivariate flows by the joint distance distribution

,”

Chaos

28

,

075302

(

2018

).

https://doi.org/10.1063/1.5010779

Google Scholar

Crossref

PubMed

35.

G.

Sugihara

,

R.

May

,

H.

Ye

,

C.-H.

Hsieh

,

E.

Deyle

,

M.

Fogarty

, and

S.

Munch

, “

Detecting causality in complex ecosystems

,”

Science

338

,

496-500

(

2012

).

https://doi.org/10.1126/science.1227079

Google Scholar

Crossref

PubMed

36.

D.

Harnack

,

E.

Laminski

,

M.

Schünemann

, and

K. R.

Pawelzik

, “

Topological causality in dynamical systems

,”

Phys. Rev. Lett.

119

,

098301

(

2017

).

https://doi.org/10.1103/PhysRevLett.119.098301

Google Scholar

Crossref

PubMed

37.

X.

Ying

,

S.-Y.

Leng

,

H.-F.

Ma

,

Q.

Nie

,

Y.-C.

Lai

, and

W.

Lin

, “

Continuity scaling: A rigorous framework for detecting and quantifying causality accurately

,”

Research

2022

,

9870149

(

2022

).

https://doi.org/10.34133/2022/9870149

Google Scholar

Crossref

PubMed

38.

M. R.

Muldoon

,

D. S.

Broomhead

,

J. P.

Huke

, and

R.

Hegger

, “

Delay embedding in the presence of dynamical noise

,”

Dyn. Stab. Syst.

13

,

175

–

186

(

1998

).

https://doi.org/10.1080/02681119808806259

Google Scholar

Crossref

39.

J.

Stark

,

D. S.

Broomhead

,

M. E.

Davies

, and

J.

Huke

, “

Delay embeddings for forced systems. II: Stochastic forcing

,”

J. Nonlinear Sci.

13

,

519

–

577

(

2003

).

https://doi.org/10.1007/s00332-003-0534-4

Google Scholar

Crossref

40.

J.

Stark

, “

Delay embeddings for forced systems. I: Deterministic forcing

,”

J. Nonlinear Sci.

9

,

255

–

332

(

1999

).

https://doi.org/10.1007/s003329900072

Google Scholar

Crossref

41.

I.

Goodfellow

,

Y.

Bengio

, and

A.

Courville

,

Deep Learning

(

The MIT Press

,

Cambridge, MA

,

2016

).

Google Scholar

42.

F.

Takens

, “Detecting strange attractors in turbulence,” in Dynamical Systems and Turbulence, Lecture Notes in Mathematics Vol. 898, edited by D. A. Rand and L. S. Young (Springer Verlag, 1981), pp. 366–381.

43.

T.

Sauer

,

J. A.

Yorke

, and

M.

Casdagli

, “

Embedology

,”

J. Stat. Phys.

65

,

579

–

616

(

1991

).

https://doi.org/10.1007/BF01053745

Google Scholar

Crossref

44.

H.

Whitney

, “

Differentiable manifolds

,”

Ann. Math.

37

,

645

–

680

(

1936

).

https://doi.org/10.2307/1968482

Google Scholar

Crossref

45.

M. B.

Kennel

,

R.

Brown

, and

H. D. I.

Abarbanel

, “

Determining embedding dimension for phase-space reconstruction using a geometrical construction

,”

Phys. Rev. A

45

,

3403

–

3411

(

1992

).

https://doi.org/10.1103/PhysRevA.45.3403

Google Scholar

Crossref

PubMed

46.

U.

Parlitz

, “

Detecting generalized synchronization

,”

Nonlinear Theory Appl. IEICE

3

,

113

–

127

(

2012

).

https://doi.org/10.1587/nolta.3.113

Google Scholar

Crossref

47.

U.

Parlitz

,

L.

Junge

, and

L.

Kocarev

, “

Subharmonic entrainment of unstable period orbits and generalized synchronization

,”

Phys. Rev. Lett.

79

,

3158

–

3161

(

1997

).

https://doi.org/10.1103/PhysRevLett.79.3158

Google Scholar

Crossref

48.

L. M.

Pecora

,

T. L.

Carroll

, and

J. F.

Heagy

, “

Statistics for mathematical properties of maps between time series embeddings

,”

Phys. Rev. E

52

,

3420

–

3439

(

1995

).

https://doi.org/10.1103/PhysRevE.52.3420

Google Scholar

Crossref

49.

N. F.

Rulkov

,

V. S.

Afraimovich

,

C. T.

Lewis

,

J. R.

Chazottes

, and

A.

Cordonet

, “

Multivalued mappings in generalized chaos synchronization

,”

Phys. Rev. E

64

,

016217

(

2001

).

https://doi.org/10.1103/PhysRevE.64.016217

Google Scholar

Crossref

50.

J. M.

Amigó

,

P.

Kloeden

, and

A.

Giménez

, “

Switching systems and entropy

,”

J. Differ. Equ. Appl.

19

,

1872

–

1888

(

2013

).

https://doi.org/10.1080/10236198.2013.788166

Google Scholar

Crossref

51.

J. M.

Amigó

,

P.

Kloeden

, and

A.

Giménez

, “

Entropy increase in switching systems

,”

Entropy

15

,

2363

–

2383

(

2013

).

https://doi.org/10.3390/e15062363

Google Scholar

Crossref

52.

M. F.

Barnsley

and

S.

Demko

, “

Iterated function systems and the global construction of fractals

,”

Proc. R. Soc. London, Ser. A

339

,

243

–

375

(

1985

) See https://www.jstor.org/stable/2397690.

Google Scholar

53.

P. E.

Kloeden

and

M.

Rasmusen

, “

Nonautonomous Dynamical Systems

(

American Mathematical Society

,

Providence, RI

,

2011

).

Google Scholar

Crossref

54.

L.

Arnold

,

Random Dynamical Systems

(

Springer Verlar

,

Berlin

,

2003

).

Google Scholar

55.

Y.

Hirata

,

J. M.

Amigó

,

S.

Horai

,

K.

Ogimoto

, and

K.

Aihara

, “

Forecasting wind power ramps with prediction coordinates

,”

Chaos

31

,

103105

(

2021

).

https://doi.org/10.1063/5.0061705

Google Scholar

Crossref

PubMed

56.

K.

Petersen

,

Ergodic Theory

(

Cambridge University Press

,

Cambridge

,

2000

).

Google Scholar

57.

O. I.

Moskalenko

,

A. E.

Hramov

,

A. A.

Koronovskii

, and

A. A.

Ovchinnikov

, “

Effect of noise on generalized synchronization of chaos: Theory and experiment

,”

Eur. Phys. J. B

82

,

69

–

82

(

2011

).

https://doi.org/10.1140/epjb/e2011-11019-1

Google Scholar

Crossref

58.

D.

Mønster

,

R.

Fusaroli

,

K.

Tylén

,

A.

Roepstorff

, and

J. F.

Sherson

, “

Causal inference from noisy time-series data—Testing the convergent cross-mapping algorithm in the presence of noise and external influence

,”

Future Gener. Comput. Syst.

73

,

52

–

62

(

2017

).

Google Scholar

Crossref

59.

G.

Datseris

and

U.

Parlitz

,

Nonlinear Dynamics

(

Springer Nature

,

2022

).

Google Scholar

Crossref

60.

K.

Lehnertz

and

H.

Dickten

, “

Assesing directionality and strength of coupling through symbolic analysis: An application to epilepsy patients

,”

Philos. Trans. R. Soc. A

373

,

20140094

(

2015

).

https://doi.org/10.1098/rsta.2014.0094

Google Scholar

Crossref

61.

H.

Dickten

,

S.

Porz

,

C. E.

Elger

, and

K.

Lehnertz

, “

Weighted and directed interactions in evolving large-scale epileptic brain networks

,”

Sci. Rep.

6

,

34824

(

2016

).

https://doi.org/10.1038/srep34824

Google Scholar

Crossref

PubMed

62.

H.

Osterhage

,

F.

Mormann

,

T.

Wagner

, and

K.

Lehnertz

, “

Measuring the directionality of coupling: Phase versus state space dynamics and applications to EEG time series

,”

Int. J. Neural Syst.

17

,

139

–

148

(

2007

).

https://doi.org/10.1142/S0129065707001019

Google Scholar

Crossref

PubMed

63.

H.

Osterhage

,

F.

Mormann

,

T.

Wagner

, and

K.

Lehnertz

, “

Detecting directional coupling in the human epileptic brain: Limitations and potential pitfalls

,”

Phys. Rev. E

77

,

011914

(

2008

).

https://doi.org/10.1103/PhysRevE.77.011914

Google Scholar

Crossref

64.

C.

Rieke

,

K.

Sternickel

,

R. G.

Andrzejak

,

C. E.

Elger

,

P.

David

, and

K.

Lehnertz

, “

Measuring nonstationarity by analyzing the loss of recurrence in dynamical systems

,”

Phys. Rev. Lett.

88

,

244102

(

2002

).

https://doi.org/10.1103/PhysRevLett.88.244102

Google Scholar

Crossref

PubMed

65.

C.

Rieke

,

F.

Mormann

,

R. G.

Andrzejak

,

T.

Kreuz

,

P.

David

,

C. E.

Elger

, and

K.

Lehnertz

, “

Discerning nonstationarity from nonlinearity in seizure-free and preseizure EEG recordings from epilepsy patients

,”

IEEE Trans. Biomed. Eng.

50

,

634

–

639

(

2003

).

https://doi.org/10.1109/TBME.2003.810684

Google Scholar

Crossref

PubMed

66.

T.

Schreiber

, “

Measuring information transfer

,”

Phys. Rev. Lett.

85

,

461

–

464

(

2000

).

https://doi.org/10.1103/PhysRevLett.85.461

Google Scholar

Crossref

PubMed

67.

M.

Staniek

and

K.

Lehnertz

, “

Symbolic transfer entropy

,”

Phys. Rev. Lett.

100

,

158101

(

2008

).

https://doi.org/10.1103/PhysRevLett.100.158101

Google Scholar

Crossref

PubMed

68.

M. G.

Rosenblum

and

A.

Pikovsky

, “

Detecting direction of coupling in interacting oscillators

,”

Phys. Rev. E

64

,

045202

(

2001

).

https://doi.org/10.1103/PhysRevE.64.045202

Google Scholar

Crossref

© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International (CC BY-NC-ND) license (https://creativecommons.org/licenses/by-nc-nd/4.0/).

This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License .

Generalized synchronization in the presence of dynamical noise and its detection via recurrent neural networks

I. INTRODUCTION

II. THE NOISELESS SCENARIO

A. Embedding theorems

[Takens theorem42]

[forced Takens theorem40]

B. The cross map

C. Generalized synchronization

D. Relationship between generalized synchronization and the cross map

III. DYNAMICAL NOISE AS STOCHASTIC FORCING

[Stark et al.39]

IV. COUPLED DYNAMICS AND NOISE

V. THE NOISY SCENARIO

A. State reconstruction

B. Cross map

C. Synchronization map

VI. NUMERICAL SIMULATIONS

A. Models

B. Results

VII. APPLICATION TO REAL-WORLD DATA: EEGs

A. Data description

B. Results

VIII. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Ethics Approval

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Generalized synchronization in the presence of dynamical noise and its detection via recurrent neural networks

I. INTRODUCTION

II. THE NOISELESS SCENARIO

A. Embedding theorems

[Takens theorem42]

[forced Takens theorem40]

B. The cross map

C. Generalized synchronization

D. Relationship between generalized synchronization and the cross map

III. DYNAMICAL NOISE AS STOCHASTIC FORCING

[Stark et al.39]

IV. COUPLED DYNAMICS AND NOISE

V. THE NOISY SCENARIO

A. State reconstruction

B. Cross map

C. Synchronization map

VI. NUMERICAL SIMULATIONS

A. Models

B. Results

VII. APPLICATION TO REAL-WORLD DATA: EEGs

A. Data description

B. Results

VIII. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Ethics Approval

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only