Inferring Causality in Synchronised Chaotic Systems using Convergent Cross Mapping

(1)

Inferring Causality in Synchronised Chaotic Systems using

Convergent Cross Mapping

Ildefonso Ferreira Pica

September 26, 2020

(2)

Abstract

When studying complex, multi-component dynamical systems we are often interested in under-standing the coupling scheme of the different subsystems. A general assumption is that a better understanding of the coupling relations i.e. their directionality, and strength amounts to causal understanding. The ultimate approach to characterize the causal relationships between system components is to identify the underlying physical mechanisms, and describe these using mathemat-ical models. In experimental settings, however, we typmathemat-ically have only partial knowledge about the physical mechanisms and must therefore rely on measurements to determine causal relationships. In this work we have focused on evaluating the effectiveness of CCM for causal inference in synchronised systems. Specifically, we set out to identify under which coupling conditions CCM produces false positive results, and whether or not the delay CCM method can successfully prevent false causal inference. To achieve this goal, we performed numerical simulations of discrete chaotic systems with different topologies, and for varying coupling strengths. More specifically, we examined a system of two identical H´enon maps, with a single driver and response system, without feedback. In addition, we examined a second system of identical H´enon maps, with a single driver providing common input to two identical, but isolated response systems.

We tested both the standard and delay CCM methods, and found that while the standard CCM method produced several false positive results (particularly for the system with common forcing), the delay CCM method was able to considerably reduce this quantity for the unidirec-tional coupling system, and only slightly for the system with common forcing. We also tested two measures for synchrony (mean phase coherence [1], and coarse-grained information rate [2]) with the idea that these might complement the CCM analysis, and help point out cases where the method fails due to synchrony.

Title: Inferring Causality in synchronised Chaotic Systems using Convergent Cross Mapping Author: Ildefonso Ferreira Pica, ildefonso.ferreirapica@student.uva.nl

Supervisors: Jo˜ao Afonso, PhD Alfonso Renart, PhD Second grader: Rick Quax, PhD Date: September 26, 2020

Graduate School of Informatics Faculty of Science

University of Amsterdam

(3)

Acknowledgements

I would like to thank my supervisors Jo˜ao Afonso and Alfonso Renart for giving me the oppor-tunity to work on this interesting project, and for being very supportive throughout throughout my project. I enjoyed working with them, and I am grateful for their advice and the interesting discussions we’ve had.

I would also like to thank the other members of the Renart lab including: Raphael Steinfeld, Mafalda Valente, Andr´e Monteiro and Davide Reato, for providing a welcoming and stimulating atmosphere. And I would particularly like to thank two other members of the lab, Tiago Costa, and Juan Casti˜neiras for the nice discussions we’ve had on a variety of interesting topics, and for helping me gain more insight into some of the theoretical aspects of my project.

(4)

1 Introduction

When studying complex, multi-component dynamical systems we are often interested in under-standing the coupling scheme of the different subsystems. A general assumption is that better understanding of the coupling relations i.e. their directionality and strength, amounts to causal understanding. The ultimate approach to characterize the causal relationships between system components is to identify the underlying physical mechanisms, and describe these using mathemat-ical models. In experimental settings, however, we typmathemat-ically have only partial knowledge about the physical mechanisms and must therefore rely on measurements to determine causal relationships. One popular framework for causal inference from observational data is that of (Wiener-)Granger causality [3, 4], which holds that a variable (the driver) is causally connected to another variable (the response) if the driver signal significantly assists in predicting the future of the response signal, beyond the degree to which the response predicts itself. Granger causality is typically implemented using autoregressive models [4, 5], but has been extended to methodologies that can deal with feedback (bidirectional causal relations) and nonlinear relations between variables. These methods have found widespread use in various fields ranging from economics [5] to neuro-science [6, 7]. An important property of the linear Granger methods are that they assume signals to be separable [8]. Separability implies that a response signal is a linear combination of the dynamics of the response system and its driver, meaning that the signal of candidate drivers can be excluded from the data set to determine their impact on the predictability of the response sys-tem’s behaviour. This separability property does not, in general, apply to deterministic nonlinear systems [9].

Recently, Sugihara et al. [9] introduced a statistical method called Convergent Cross Mapping (CCM), specifically aimed at causal inference for deterministic nonlinear systems. The method is based on Takens’ theorem [10], an important result from dynamical systems theory, which infor-mally states that the attractor of a dynamical system can be reconstructed from time dependent observations of the dynamics of the system, by using time lags as the axes of the phase space. The attractor is a very helpful mathematical object as it characterizes the long-term behaviour of a dynamical system. Importantly, Takens’ theorem states that we can define one-to-one map-pings between the attractors reconstructed from scalar observations of different subsystems that are coupled together. CCM is designed to test the existence of these one-to-one mappings (cross mapping), and thus infer causality. The method is related to the method of analogs [11] used by Lorenz for nonlinear predictions with atmospheric data [12], and the nonlinear interdependence measure used by Quiroga et al. [13] to measure synchronisation in EEG signals. What sets CCM apart from these earlier methods, is that it involves the explicit test of a property labeled by Sugihara et al. [9] as convergence. Convergence describes the increase in cross mapping ability, as a function of the number of data points used for the attractor reconstructions. The idea is that it a one-to-one mapping exists between different subsystems, then the addition of more data will fill in the attractors, thus making the cross mapping more accurate. If, however, no relation between the subsystems exists, then the addition of more data will not result in more accurate cross mapping. Convergence, thus becomes the primary condition for causal inference using CCM. The CCM method has been applied in a wide range of contexts ranging from the study of the causal relation between sea surface temperature and the size of fish populations [9, 14], to the predictability of synchronised bursting behaviour in groups of neurons using the signals of single neurons [15], and the examination of the causal effects of different global environmental factors in the occurrence of influenza [16]. Given that all these cases involve the dynamic behaviours of many interlinked complex systems, it is very encouraging and somewhat surprising that the CCM method works. Even more so, because it does not require any explicit knowledge of the physical mechanisms that govern the system dynamics. However, some authors have raised concerns about the reliability of the method in cases where different systems are (strongly) driven by common input [17]. Indeed, CCM has been shown to produce false positive results for strong coupling relations, where different systems become synchronised [9, 18].

Several attempts have been made to fix this shortcoming. One proposed solution, called extended CCM (referred to as delay CCM in this report) introduced by Ye et al. [18], adds the

(6)

condition that the cross mapping accuracy should be maximized when the time series are shifted to account for the delay in the coupling relation between the systems. With this, Ye et al. [18] have essentially encoded the condition that causes should precede effects. A second approach, introduced by Deyle et al. [16], involves the use of specially constructed surrogates to deal with entrainment to periodic behaviour by a common driver e.g. seasonality.

In this project we have focused on evaluating the effectiveness of CCM for causal inference in synchronised systems. Specifically, we set out to identify under which coupling conditions CCM produces false positive results, and whether or not the delay CCM method can successfully pre-vent false causal inference. To achieve this goal, we performed numerical simulations of discrete chaotic systems with different topologies, and for varying coupling strengths. More specifically, we examined a system of two identical H´enon maps, with a single driver and response system, with-out feedback. In addition, we examined a second system of identical H´enon maps, with a single driver providing common input to two identical, but isolated response systems. We characterized the synchronous regimes for these systems by estimating their maximal sub-Lyapunov exponents [19, 20], as a function of the coupling strength. We tested both the standard and delay CCM methods, and found that while the standard CCM method produced several false positive results (particularly for the system with common forcing), the delay CCM method was able to consider-ably reduce this quantity for both the unidirectional coupling system, the system with common forcing. We also tested two measures for synchrony (mean phase coherence [1], and coarse-grained information rate [2]) with the idea that these might complement the CCM analysis, and help point out cases where the method fails due to synchrony.

2 Convergent Cross Mapping

In this chapter we discuss the convergent cross mapping (CCM), a statistical method introduced by Sugihara et al. [9] to infer causality in dynamical systems, from time series data. We start with a discussion of it’s theoretical underpinnings in the Takens’ embedding theorem. We continue with a high level summary of the method, followed by a small exposition on the expected outcomes of a CCM analysis for various system topologies. We also discuss an extension of the standard CCM method, here called the delay CCM method, designed by Ye et al. [18] to more accurately identify causal drivers. We conclude this chapter by considering the quantity of data required for CCM.

2.1 Takens’ embedding theorem

In dynamical systems theory, causality can be understood as the dynamical coupling between variables that belong to the same dynamic system. Together, the variables of a coupled dynamical system define a manifold in phase/state space, which describe the behaviour of the system at equilibrium. Intuitively, a manifold is a topological space that locally resembles Euclidean space, but need not resemble Euclidean space globally.

Convergent Cross Mapping tests for dynamical coupling between variables by quantifying the extent to which the historical records of one dynamic variable can estimate the historical records of another dynamic variable. This estimate is unidirectional, meaning that two estimates can be defined i.e. by swapping the predictor and predicted variable. An asymmetry in the accuracy of these estimates, is an indication of an asymmetry in the dynamical relation between the variables in terms of coupling strength. In turn, the differences in directional coupling strength and/or time delays in the coupling relations, are interpreted as causal relations within the CCM-framework.

The cross-prediction of dynamics using historical records of a different variable relies on two principles, firstly, that like antecedents produce like consequents, which can be understood both as a deterministic/causal principle (i.e. regularity approach [21]) and a dynamical systems principle (i.e. attractors), and secondly, that the full behaviour of a system can be estimated using (scalar) time-dependent observations from at least one system variable, under certain conditions. This second principle is formalised in Takens’ embedding theorem [10], which entails a procedure for attractor reconstruction.

(7)

In dynamical systems theory, an attractor is a set of n-dimensional vectors to which a dynam-ical system evolves for a set of initial conditions. Once on the attractor, the dynamdynam-ical behaviour of the system will be determined by the attractor, unless the system is forced off the attractor e.g. by noise or an external force. Chaotic systems are characterised by so-called strange attractors. Dynamical systems with strange attractors are sensitive to initial conditions, meaning that arbi-trarily small differences in the initial conditions can lead to significantly different future behaviour of the system. Many biological systems display this chaotic property, including (populations of) neurons.

Since the long-term behaviour of chaotic systems are characterised by their attractors, a lot can be gained by analysing these sets. It is therefore not surprising that a lot of work has been put into developing methods for reconstructing attractors from empirical data. Of particular interest here is the delay embedding method, introduced by Floris Takens [10]. A delay embedding can be constructed for a time series X = X(1), . . . , X(T ), by composing E-dimensional vectors through concatenating E time points that are multiples of τ apart, as follows:

x(t) = [X(t), X(t − τ ), X(t − 2τ ), . . . , X(t − (E − 1)τ )],

where E is a positive integer greater than one called the embedding dimension, and τ is a non-zero integer denoting the lag delay. This can be done for all T − (E − 1) points of the time series. If E (and in practice τ ) is chosen appropriately, the resulting time delay embedding yields a so-called shadow manifold MX that has a 1:1 mapping to the original manifold M of the dynamical system.

In addition, MX is a diffeomorphic reconstruction of the attractor.

A shadow manifold can be created from each system variable, and since there is a bijection between each shadow manifold and the original attractor manifold, Takens’ theorem implies that there must be a bijection between each shadow manifold as well, as long as the system variables are coupled. The CCM method operationalises this fact to infer the directionality of the coupling between variables, and by extension infer their causal relations.

2.2 Convergent cross mapping (CCM)

CCM is a method that quantifies the dynamic coupling between different variables of a dynamic system, via state space reconstruction. It does so by measuring the accuracy of mappings between the shadow manifolds constructed for the variables of interest. Let X(T ) and Y (T ) be two dy-namically coupled variables, M the attractor manifold of the total system with dimension d, the embedding dimension E = 2d + 1, and τ the time delay, then CCM is performed by constructing the shadow manifolds MX and MY through time delay embedding with parameters E and τ , and

subsequently making cross predictions for varying library sizes L = L1, . . . , Ln. The library points

are randomly sampled from the predictor shadow manifold e.g. MX, where the minimum number

of points is equal to one plus the neighbourhood size (K) selected for the cross mapping. A typical value for the neighbourhood size is E + 2, which produces the smallest possible bounding simplex. The predictor library is used to estimate the scalar values of Y , denoted as Y |MX. This

is achieved by iterating over the points in Y , then selecting the concurrent vector X in MX,

and selecting its K nearest neighbours. The K nearest neighbours are regressed to X with an exponentially decaying function, so as to give more importance to the most similar vectors. This yields regression weights W_XK. The vectors in MY that are concurrent to the K nearest neighbour

vectors, are summed using weights W_XKto produce the estimate Y |MX. The scalar estimate Y |MX

is obtained by selecting the first element of the predicted delay embedded vectors. The accuracy of the prediction is typically estimated using the Pearson correlation or the Mean Absolute Error (MAE) between the original and predicted signal.

According to Sugihara et al. [9], causality can be inferred by analysing the change in prediction accuracy of the cross mapping, as the library size increases (we will refer to this notion of causality as Sugihara causality, throughout this report). As more points are added to the reconstructed predictor manifold, the distance between the points decreases, and increasingly better analogs are added. If the cross mapped variables are dynamically coupled, the better analogs in the

(8)

predictor manifold should correspond to better concurrent analogs in the target manifold, thereby increasing the accuracy of the cross mapping. This has been labeled ”convergence” by Sugihara et al. [9]. On the other hand, if the variables are not coupled, then there is no bijective relation between the two manifolds, and thus the inclusion of more points in the predictor manifold does not, on average, result in better analogs in the target manifold. If the variables are deterministic and dynamically coupled, and their observations (i.e. the measured time series) are sufficiently lengthy, have sufficient fidelity and are noise-free, then the cross mapping should converge to a Pearson correlation of 1. The cross mapping convergence can be estimated in both directions, which, in principle, allows for the detection of the directionality of the causal relation. A more formal description of the CCM algorithm can be found in section 4.2.1.

2.3 Expected CCM results for different causal relations

2.3.1 Bidirectional coupling

In the current context, bidirectional Sugihara causality refers to existence of bidirectional coupling between the variables of a dynamical system. This type of coupling relation produces a situation where each variable of the coupled pair influences the dynamics of the other variable. Of course, the coupling relation need not be symmetrical in terms of coupling strength i.e. one of the variables can assert a greater influence on the dynamics of the other, and in terms of coupling lag i.e. the influence of one variable can temporally precede the influence of the other variable (feedback) and the feedback influence need not be perceived instantaneously, but after some time lag.

Figure 2A shows the cross mapping between the shadow manifolds MX and MY of two

bidi-rectionally coupled variables X and Y . Due to the bidirectional coupling relation, the historical records of X contain information about the dynamics of Y , and vice versa, which can be used to make mutual predictions via cross mapping. Put differently, there is a functional relation be-tween MX and MY, which is reflected in the preservation of the closeness of the points as they

are mapped from the predictor manifold to the target manifold (blue lines). Because a functional relation exists for both cross mapping directions, we expect to observe convergence for both cross maps X|MY and Y |MX.

Sugihara et al. [9] have proposed that the information the response variable records about the driver becomes more ”distinct” as the coupling strength increases. Effectively, the shape of the attractor manifold is transformed such that the proximity of points in the response manifold more closely matches the proximity of contemporaneous points in the driver manifold (figure 3). According to Sugihara et al. [9], a stronger effect of a driver variable X on its response variable Y will imply faster convergence for the cross map X|MY. Sugihara et al. [9] rely on this

property to claim that the relative speed of convergence can be used to identify which variable of a bidirectionally coupled pair, has the greatest causal effect on the other.

2.3.2 Unidirectional coupling

In the case of unidirectional coupling, a driver variable X influences the dynamics of a response variable Y without feedback. Figure 2B shows the cross mapping between the shadow manifolds MX and MY, where X unidirectionally affects the dynamics of Y . Since X does not ’observe’ Y ,

the shadow manifold of X contains no historical records of the dynamics of Y . This is reflected by the lack of closeness in the proximal points in MX, after cross mapping to MY. The lack of

a functional relation from X to Y (the mapping from X to Y is not injective), means that we expect no convergence will occur for the cross map Y |MX. However, convergence is expected

for the cross map X|MY, since the shadow manifold of Y has a record of the dynamics of X

and a functional relation therefore exists from Y to X. This asymmetry in the CCM scores and the convergence of the cross maps, is what allows one to distinguish unidirectional coupling from bidirectional coupling.

(9)

Figure 2: Cross mapping for cases of bidirectional and unidirectional coupling. A: cross mappings between the shadow manifolds of two bidirectionally coupled variables (X ↔ Y ). B: same as A but for two unidirectionally coupled variables (X → Y ). Images taken from Sugihara et al. [9].

2.3.3 Common forcing

In the simplest case of common forcing, two response systems Y 1 and Y 2 receive input from a single driver variable X. The only coupling relations are thus the unidirectional couplings between X and the response systems. The shadow manifolds MY 1 and MY 2 both contain records of the

dynamical behaviour of X, and unidirectional convergence is thus expected for cross mappings X|MY 1 and X|MY 2 only. Since there is no direct or indirect coupling relation between systems

Y 1 and Y 2, neither variable should have a record of the dynamics of the other. Therefore, no convergence should be found for the cross maps between the response systems Y 1 and Y 2. These expectations hold for systems with arbitrarily many response systems, as long as there is no feedback between any of the response systems and the driver variable. In the case of feedback, the transitivity of the system dynamics will ensure that information about the feedbacking response system is passed on to the other response systems, via the driver.

In general, the response systems do not need to have the same governing equations, nor do the couplings between the driver and the various response systems need to be equal strong. In section 3.2 we will discuss the consequences of having equal response systems versus response systems with different governing equations, for cases of strong coupling. In this report we will further study the simplest case of a system with common forcing, consisting of two equal response systems that are equally strongly coupled to a single driver system (section 4.1).

2.4 Delay CCM

An important weakness of CCM is that it is not able to reliably identify the causal driver in a system if the variables are in a synchronised regime due to strong coupling [9]. In these cases, the cross mapping ability can show bidirectional convergence, even though the driving variable is only unidirectionally coupled to the driven variable. This makes CCM an unreliable method for the inference of causal directionality in cases where there is no way to verify the strength of the dynamical coupling between system variables.

(10)

Figure 3: Transformation of the attractor manifold of two unidirectionally coupled chaotic systems. As the coupling strength increases, the shape of the manifold of the response variable (Lorenz system) more close matches the shape of the manifold of the driver variable (Roessler system). In addition, the proximity of points in the response manifold more closely matches the proximity of contemporaneous points in the driver manifold. Image taken from Quiroga et al. [13].

An extension to CCM has been proposed by Ye et al. [18] to solve the problem of unreliable causal inference in synchronised systems, and the detection of proper causal relations in systems with time-delayed interactions. The suggested common solution to these problems is to explicitly consider time lags during cross mapping. This amounts to analysing the cross mapping perfor-mance as a function of the time displacement of one time series relative to the other, an approach analogous to the cross-correlation. According to Ye et al. [18], optimal cross mapping is achieved for time shifts proportional to the delay time of the coupling interactions between the variables. Importantly, this means that a driven variable can be identified by observing a higher cross map-ping ability for negative time shifts with respect to the driver, as this indicates that the future dynamical behaviour of the driven variable is determined by the prior behaviour of the driver, at a (fixed) lag determined by the time-delayed interactions.

Ye et al. [18] found that in the case of bidirectional coupling, the optimal cross map lags were negative for both cross mapping directions (figure 4A-C). Greater delay times for the feedback relation, resulted in a shift of the optimal cross map lag in proportion with the delay time. The authors also found that the driver and response variables could be distinguished using the delay CCM method, even in the case of generalized synchrony due to strong unidirectional coupling (figure 4D).

2.5 Data quantity required for CCM

Cross mapping relies on the availability of enough good analogs (proximal on the attractor), which implies a requirement for a high enough density of points in the reconstructed attractor. By extension, this implies that the behaviour of the system must be sampled for sufficient length to ensure that there are analogs in the first place. This last point can become problematic, especially if the behaviours of interest occur infrequently.

The amount of data required for CCM is dependent on the dimensionality of the system under consideration, and the desired accuracy of the cross mapping. Let us consider a time series X = x1. . . , xN of points sampled at a time interval δt. We now denote the number of analogs

(11)

Figure 4: Model demonstration of causal lags and optimal cross mapping using a 2-species logistic system with bidirectional coupling and a system with strong unidirectional coupling. A-C: mean cross map skill as a function of the coupling strength for various causal lags τd. The results were

obtained for a delay embedding with parameters E = 2 and τ = 1. 100 random libraries consisting of 200 vectors were used for cross mapping, and 1000 points from a hold-out set were predicted. D: same as A-C but for a unidirectional system with strong coupling c = 0.8. Images taken from Ye et al. [18].

within a small radius around a point xi with K(). The average time between two such points

is then given by:

τk =

(N − 1)δt

K() . (1)

This gives us an estimate of the amount of time needed for the system to return a state that is within a distance removed from the system state at xi. We can now express the fraction of

-analogs of a point xi over the total number of data points N as:

Ck() =

K()

N − 1, (2)

which given equation (1) gives us:

Ck() =

δt τk

. (3)

Now in order to make a cross mapping prediction with an accuracy δ, we need at least N ≥ τN.

Given the scaling law due to Grassberger and Procaccia [22]:

(12)

where C() is the correlation sum defined as: C() = lim N →∞ 1 N2 N X i,j=1 Θ( − kxi− xjk), (5)

we get that the number of points needed to make a δ-accurate prediction scales as: N ∼ −D, where D corresponds with the correlation dimension which is an estimate of the dimension of the subspace of the phase space that is spanned by the system true, called the effective dimension DM. Importantly, this means that for large systems, prediction becomes intractable for reasonable

fidelity . An additional complication is that both D and are typically unknown, particularly when one only has partial observations of a dynamical system.

3 Synchrony in chaotic systems

In this chapter we give an overview of the phenomenon of synchrony in chaotic systems. We describe two different forms of synchrony which produce different dynamical relations between coupled chaotic systems, in a coupling-dependent way. More specifically, the strength of the coupling relation between the different chaotic systems or system components, is what primarily determines the type of synchronisation that can be observed. In the final section of this chapter (section 3.3) we discuss the consequences of synchronisation for the performance of CCM analyses.

3.1 Weak coupling: phase synchronisation

When considering synchronisation in chaotic systems, we need to start by specifying a notion of periodicity, which is non-trivial because chaotic signals are not strictly periodic in a mathematical sense. Nevertheless, some chaotic signals can be considered ’almost periodic’ [20]. These signals then consist of cycles that appear similar, despite (small) differences in their amplitude and period [20]. When this happens, we can count the number of cycles that occur within a sufficiently large time interval to compute the mean frequency: ¯f = NT/T . Here T denotes the length of the time

interval, and NT denotes the number of cycles within that time interval. Mean frequencies can

be used to describe the behavior of a system composed of chaotic subsystems in the same way periodic systems can be described. If the coupling between chaotic oscillators is strong enough, their mean frequencies become equal, however this does not imply that the signals coincide (figure 5B). Weak interactions between chaotic systems do not interfere with their chaotic behavior; the amplitudes of coupled chaotic oscillators remain irregular and uncorrelated.

There are several ways to define the phase of chaotic signals. One way follows from the aforementioned approach of defining cycles as orbits that pass through a particular region in phase space. In principal, this region can be arbitrarily small, but when considering empirical measurements one must define a region of sufficient size to have cycles. A notion of phase can then be formulated by saying that each cycle corresponds to a phase gain of 2π [23, 20]. Since the cycle lengths are irregular, the phase rotates non-uniformly i.e. it accelerates and decelerates irregularly [20]. As the frequency of the system dynamics become entrained through weak coupling, we can consider the frequency adjustments as phase shifts between the chaotic signals. This form of weak synchronisation is called phase synchronisation [20].

Phase synchronisation can also be examined for chaotic systems, without defining cycles, by estimating the phase angle from the analytic signal obtained via the Hilbert transform [23, 1]. Using this procedure [1] have defined a measure for phase synchronisation called the mean phase coherence. We have used this measure to examine the occurrence of phase synchronisation in a chaotic system with common forcing. Further details on the phase synchronisation analysis can be found in section 4.3.2, while the results can be found in section 5.4.

(13)

Figure 5: An example of phase and complete synchronisation between two chaotic signals. The left plots show example traces, while the right plot shows the signal of X2as a function of X1. A: the

two chaotic systems are in an uncoupled state, the right plot shows that there is no distinguishable relationship between the two signals. B: the systems are in the phase synchronisation regime. A circular structure has now become visible in the right plot, indicating that the mean frequencies of the systems have become entrained. Note that the amplitudes of the signals remain uncorrelated while the signals remain approximately in phase. C: the systems are completely synchronised due to strong coupling. The two signals have become identical, as can be clearly seen in the right plot. Image taken from Pikovsky et al. [20].

(14)

3.2 Strong coupling: complete and generalized synchronisation

Complete synchronisation, occurs between identical systems (with identical governing equations) under strong coupling conditions. This forces the dynamical behaviour of the coupled chaotic systems to become identical, after a finite transient, regardless of their initial conditions [19, 20]. In this case, both the mean frequencies and the amplitudes become more similar between the two systems.

The difference between self-sustained and forced chaotic systems is crucial for phase synchro-nisation, but irrelevant for the effects of complete synchronisation [20]. Self-sustained chaotic systems are described by autonomous equations (where there is no direct time dependence), and consequently the dynamics are independent of time shifts (time symmetrical). However, in forced chaotic systems the dynamics are equivalent only in instants of time shifted by the period of the forcing. Pikovsky et al. [20] note that complete synchronisation refers to the suppression of the differences between coupled identical systems, which is different from the process of entrainment or locking seen in periodic synchronisation. Complete synchronisation can only occur if the coupling strength exceeds a threshold that is proportional to the size of the largest Lyapunov exponent of the individual systems.

Complete synchrony is a special case of generalized synchrony, which is a functional relation between the states of two chaotic systems: x2 = F (x1), where xi are the time-dependent state

vectors of the chaotic system i [20]. In the case of complete synchronisation F is the identity function [19, 20]. Generalized synchrony is typically observed for unidirectionally coupled systems, also called master-slave systems, due to the fact that the dynamics of the driven variable become fully determined by the driver. The determination of synchrony in noisy systems is ambiguous, since the border of synchronisation will be smeared [20]. Moreover, in the case of a passive experiment, where the relevant parameters governing coupling and frequency detuning are not controlled by the experimenter, Pikovsky et al. [20] states that analysis can only detect the presence of an interaction between the observed systems, and not synchronisation. According to Pikovsky et al. [20], this is because synchronisation is a dynamical process, and not a single state.

In the case of complete synchronisation, it is not possible to determine the direction of informa-tion flow or, by extension, causality between synchronised subsystems, from data alone [2]. This is because under complete synchronisation the dynamics of the subsystems, and therefore their asso-ciated time series, become identical to that of their driver (figure 5C). According to Paluˇs et al. [2], the same holds for the case of generalized synchronisation, where the subsystem associated time series become related through a one-to-one nonlinear mapping. Paluˇs et al. [2] have introduced an information theoretic method (coarse-grained information rate) for the identification of complete and generalized synchrony. We have used this measure to examine the occurence of complete synchronisation in a chaotic system with common forcing. Further details on the coarse-grained information rate can be found in section 4.3.3, while the results of the method for our simulation experiment can be found in section 5.4.

3.3 CCM and synchrony

In a direct reply to Baskerville and Cobey [17], Sugihara et al. [24] address the application of CCM in cases where a response variable (is expected to) synchronise to a driver variable, like for the seasonal effects of environmental variables on the incidence of influenza [16]. Sugihara et al. [9] have previously explained that CCM can give the false appearance of bidirectional causality, in unidirectional systems with complete and generalized synchrony.

Deyle et al. [16] have addressed this issue by introducing a bootstrap method dubbed seasonal surrogates, to help distinguish driving effects from mutual effects caused by a common seasonal driver. Seasonal surrogates were computed for driving variables D(t) by modelling the seasonal behaviour of D(t) through averaging. This produced the seasonal model bD(t). Next, Deyle et al. [16] computed the seasonal anomalies as the difference between the observed values {D(t)} and the seasonal model: { eD(t)} = { bD(t)} − {D(t)}. The seasonal anomalies were randomly reordered, and added to the seasonal model to produce seasonal surrogate time series D∗_{. The CCM scores}

(15)

computed for the seasonal surrogates, were then used to define a bootstrap distribution. Deyle et al. [16] interpreted a significantly better CCM prediction on the observed data, as proof for a causal effect of the driver variable in addition to any seasonal effects. A more thorough analysis of the applicability of seasonal surrogates in cases of synchrony is required. Deyle et al. [16] used seasonal surrogates in a situation where effects of a common driver (CD) were expected, and a model for these effects could be defined. However, this is unlikely to be the case in general. Seasonal surrogates are unlikely to provide a solution when no appropriate CD model or explicit measurements of the CD exist. And when a model is defined the method of seasonal surrogates may still fail to help, as there currently exists no general way to test the validity of a CD model. These issues need to be resolved before CCM can be reliably used to distinguish synchrony in (non)interacting variables, from true causal effects.

In cases where the effects of synchrony on CCM results cannot be resolved e.g. using seasonal surrogates, it is critical to detect whether data was collected from a system in the synchronous regime. This is precisely what we aimed to do. In section 5.4 we show the results of our synchro-nisation analysis for a chaotic system with common forcing, that undergoes synchrosynchro-nisation as the coupling strength between the driver and response variables are increased.

4 Materials and Methods

All analyses were performed in MATLAB R2019a [25] using custom-written scrips, except for the CIR analyses, which were performed in R v3.5.3 [26] using custom-written scrips and the infotheo package [27].

4.1 Model system: coupled H´

enon maps

To test the effect of synchronisation on the performance of the CCM method, we studied a chaotic system composed of three Hénon maps. The Hénon map [28] is a discrete-time dynamical system that captures the most essential features of the famous Lorenz system [29], while introducing a tuneable dissipation parameter [30]. The Hénon map has several nice properties, including that it is time invertible, meaning that each point has a unique past; it is dissipative, meaning that it contracts areas in phase space, and it has a trapping region in which the attractor lies [30]. The fact that the this chaotic system is discrete, aided in the analysis of its dynamics using state-space reconstruction (e.g. time delay embedding) based methods, as the reconstruction process is more straightforward for discrete systems than for continuous systems. Specifically, the estimation of proper time delay embedding parameters is more involved for continuous systems, where special care must be taken to define an appropriate time delay τ . In contrast, for discrete chaotic systems τ = 1 is a natural choice.

We defined a system with a single driver subsystem X that provided unidirectional input to two identical response subsystems Y1 and Y2. This particular system topology has been labelled

by Pearl et al. [31] as a ”fork”. The driver subsystem X was defined as:

x(t + 1) = a − x(t)2+ bu(t), (6)

u(t + 1) = x(t) (7)

while the response systems, Y1 and Y2, were defined as:

y(t + 1) = a −cx(t) + (1 − c)y(t)y(t) + bv(t), (8)

v(t + 1) = y(t) (9)

where c ∈ [0, 1] is the coupling parameter that determines the input strength from the driver sub-system X to subsub-systems Y1and Y2. We chose parameter values a = 1.4 and b = 0.3 corresponding

with the classical H´enon map, which displays chaotic behaviour. Figure 6A shows an example signal for the driver and one of the response subsystems, for moderate coupling (c = 0.4), while figures 6(B-C) show the corresponding strange attractors of the respective subsystems.

(16)

Figure 6: Example signals and the attractors of the coupled H´enon maps, for a moderate coupling strength (c = 0.4). A: example signals consisting of 100 map iterations of the driver system (blue) and the response system (red). B-D: the attractors of the driver (blue) and response systems (red) shown in 3D lagged coordinate space. The attractors are overlayed in D to show their relative positions.

(17)

4.2 Convergent Cross Mapping

4.2.1 The CCM algorithm

The CCM method combines two principles to perform causal inference: cross mapping, and con-vergence. After constructing the shadow manifolds for the variables of interest, cross mapping is performed by projecting the exponentially weighted average of the nearest neighbours of all points in the predictor manifold to the target manifold. This yields a prediction that improves as more data is added to predictor manifold, if the predictor variable is dynamically driven by the target variable (or the variables are synchronised). The improvement of the cross prediction as a function of the data length, is called convergence, and it is the critical condition for Sugihara causality. I will now give a more technical description of the CCM algorithm proposed by Sugihara et al. [9]. For further details we refer the reader to the supplementals of their paper.

We define two time series of interest X(t) and Y (t), their corresponding attractor manifold M of dimension d, the reconstruction dimension E = 2d + 1, and the time delay τ . The algorithm consists of the following steps:

1. Choose a projection direction (i.e. select the predictor variable) e.g. X.

2. Sample (randomly) L library points from X, where K ≤ L ≤ N , with K the number of nearest neighbours (the simplex projection requires at least E + 1 neighbours), and N is the total data length. To determine convergence, several libraries must be sampled of increasing size.

3. Create the shadow manifold MX from L, through time-delay embedding, by defining the

E-dimensional vectors:

X(t) = [X(t), X(t − τ ), . . . , X(t − τ (E − 1)], (10)

for each t ≥ E − 1.

4. For each point Y (i), choose the corresponding vector X(i), and selects its K nearest neigh-bours.

5. Denoting the time index of the kth_{nearest neighbour as t}

k, calculate weights wk, using the

Euclidean distances between X(i) and X(tk). The weights function proposed by Sugihara

et al. [9] is: wk= 1 W max exp −kX(tk) − X(i)k dmin , wmin , (11) with W = P

kwk, dmin = mink(kX(tk) − X(i)k), and wmin equal to a small value (e.g.

machine epsilon), to improve the computational stability.

6. Compute bY (t) as the weighted sum:

b

Y (t) = X

k≤E+1

wkY (tk). (12)

7. Finally, evaluate the accuracy of the predictions bY (t), for example by calculating the Pearson correlation coefficient (suggested by Sugihara et al. [9]). Repeat steps 3-7 for all libraries L, and determine whether there is convergence i.e. whether the prediction accuracy of the cross mapping improves as L increases.

Given the outcome of the CCM analysis, Sugihara causality is inferred if convergence can be observed. Convergence can be defined for an arbitrary similarity measure between the true and cross-map predicted time series. Selecting the Pearson correlation coefficient as our similarity

(18)

measure, we can formally define convergence as follows, let L be the length of the time series used to cross-map, if the assumptions of Takens’ theorem are satisfied, then:

lim

L→∞ρL( bY , Y ) = 1, (13)

where ρ is the Pearson correlation coefficient. Even in the case of a causal relation, cross-mapping may not always converge to ρ = 1. In particular, the presence of noise, in the dynamical system itself or due to measurement errors, can lower the convergence point. This fact is predicted by Takens’ theorem, as Takens’ theorem requires the temporal dynamics of the system to be C2_,

which is not generally true for Stochastic Differential Equations. However, if the noise is sufficiently small, the CCM approach will still function, but a perfect cross-mapping prediction will not be achieved i.e. the cross-mapping will converge to some ρ < 1.

4.2.2 CCM analysis

We performed a CCM analysis for all 101 generated time series, for the different coupling strengths. The three time series consisted of the signals corresponding to variables x(t) and y1,2(t). All time

series consisted of 5 × 103time steps. The systems were iterated for a total of 6 × 103 time steps, but the first 103 _{time steps were considered a transient period, and were not included in any}

further analyses.

We performed all time delay embeddings using an embedding dimension E = 4, and time delay τ = 1. We used K = E + 1 = 5 nearest neighbours for the cross mapping, which is the minimal number of nearest neighbours required. To counteract spurious cross mapping due to the proximity of subsequent time points on the attractor, we prevented the selection of nearest neighbors from a window of 14 time points centered at each target point (exclusion radius = 7). We obtained the value for the exclusion radius by computing the partial autocorrelation function for the driver and response systems, and selecting the largest lag with a non-zero correlation score. We varied the library size from 10 to 103 _{points, in steps of 10. The library was constructed by}

randomly sampling data points, with replacement, from the first 2 × 103 points of the time series. Predictions were made for the final 3 × 103 data points of the time series. The CCM scores were obtained by computing the Pearson correlation coefficient between the target time series, and the cross map prediction. We quantified the convergence (∆CCM; delta CCM in figures) by taking the difference between the CCM scores for the minimal and maximal library sizes. The entire CCM procedure was repeated 10 times per coupling strength. All CCM scores shown here were obtained by averaging over the repetitions. Convergence was deemed significant if the mean ∆CCM score for the true cross mapping was larger than the mean ∆CCM score plus twice the standard deviation of the ∆CCM score obtained for random cross mapping, where random neighbours were used for the cross mapping instead of nearest neighbours.

4.2.3 The delay CCM algorithm

As mentioned in section 4.2.1, the CCM method relies on two principles for causal inference: cross mapping, and convergence. However, these two principles are not sufficient conditions for causal inference, as they do not prevent false positives in cases of strong coupling between dynamical systems or coupled variables of a single dynamical system. The delay CCM method introduced by Ye et al. [18] adds an additional requirement based on the principle that effects cannot precede causes. This principle is operationalised by extending the normal CCM algorithm as follows:

1. Define a maximal cross map lag lmax.

2. Shift the time series by l ∈ [−lmax, lmax] time steps. Note that the set of cross map lags does

not need to contain all consecutive lags between negative and positive lmax.

(19)

The delay CCM method shows how the cross mapping ability depends on time shifts between the time series. As mentioned in section 2.4, Ye et al. [18] have shown that the value of the cross map lag which produces the highest cross map score is proportional to the delay time of the coupling interactions between coupled dynamical systems. More importantly, the delay CCM method yields a more trustworthy procedure for causal inference. If the optimal cross map lag for X|MY is negative and convergence occurs, then this is an indication that system X affected

the future dynamical behaviour of Y . This is because the historical records of the dynamics of X, contained in the shadow manifold of Y , correspond with the dynamics of X that temporally preceded the dynamics of Y .

There are no theoretical restrictions on the size of lmax, however larger portions of the data

become unavailable for cross mapping as l increases. This is because for a shift of size l, 2 × l data points lack a contemporaneous point in the other time series. Thus the data size effectively limits the size of the maximal cross map lag. In addition, it is undesirable to pick a large lmax,

as the cross mapping part of the CCM algorithm needs to be repeated for each lag. Ideally, an appropriate maximal cross map lag is selected on the basis of some prior knowledge about the characteristic delay times of the dynamical system under consideration. Recently an optimization method (chaotic ant swarm) has been proposed for the identification of delay times in chaotic systems. We refer the interested reader to [32] for further details.

4.2.4 Delay CCM analysis

We performed a delay CCM analysis to evaluate the ability to differentiate synchronisation between directly coupled variables, and indirect synchronisation between mutually uncoupled variables that received common input. The time delay embedding and cross mapping were performed as described in section 4.2.3. For the delay CCM analysis we selected a maximal cross map lag lmax= 10, as this

provided a good trade-off between the computational work and the determination of the optimal cross map lag.

4.2.5 CCM causal inference conditions

Both the original and delay CCM analyses have testable conditions on the basis of which causality is inferred. For the original CCM analysis we defined the following condition to infer causality:

(C1) Causality is inferred when δCCM > δCCMRN+ 2ˆσRN.

Here the l.h.s. denotes the mean δCCM score, the first term of the r.h.s.denotes the mean δCCM score obtained for cross mappings with random neighbours instead of the usual nearest neighbours, and ˆσRN denotes the corresponding sample standard deviation. All measures were obtained for

10 CCM repetitions, where the predictor library (section 4.2.2) was randomly resampled with replacement for each repeat.

The directionality of a causal relationship between two variables can be inferred using condition (C1): if the condition is met for a single cross map direction, then the causal relation is assumed to be unidirectional and the predictor of the converging cross map is assumed to be the response variable, while the predictee is assumed to be the driver variable. However, the normal CCM analysis is prone to false causal inferences in cases where different system variables are weakly synchronized.

The delay CCM analysis uses time shifts to infer the coupling delay of the causal relations between variables. Using the delay CCM analysis condition (C1) can be supplemented with the following conditions for the inference of the causal direction:

(C2) The predictor of the converging cross map is a driver variable if the optimal cross map lag is positive.

(C3) The predictor of the converging cross map is a response variable if the optimal cross map lag is negative.

(20)

These conditions were introduced by Ye et al. [18]. In the case of unidirectional coupling, conditions (C2) and (C3) are met, whereas with bidirectional coupling condition (C2) is met for both cross maps (figure 4). Ye et al. [18] showed that these conditions yield a more reliable analysis, that can better cope with weak synchronization due to strong coupling. In section 5.2 we show results of a delay CCM analysis for a unidirectional system, which are in agreement with Ye et al. [18].

In the case of a system with common coupling like the system tested here (section 4.1), normal cross mapping between the mutually uncoupled response systems can yield false causal inferences when the common driver is strongly coupled to both response systems (section 5.2). We found that the following two conditions can be used to detect certain cases where convergence leads to false causal inference:

(F1) Causality is wrongly inferred if convergence is observed for one cross mapping direction only, and the corresponding optimal cross map lag is positive.

(F2) Causality is wrongly inferred if convergence is observed for both cross mapping directions, and the corresponding optimal cross map lags are positive.

We have used conditions (F1) and (F2) with the delay CCM analysis to improve the accuracy of the causal analysis of a system with common coupling (section 5.2).

4.3 Detecting synchrony

4.3.1 Maximal sub-Lyapunov exponent

To determine whether a particular coupling strength induced synchronisation, we estimated the maximal sub-Lyapunov exponent (SLE) (also called the horizontal or transversal Lyapunov expo-nent [33, 20]) for the response systems. This quantity measures the degree to which perturbations along the chosen transversal axis (in our case this will be the axis defined by the driver system X) grow or decay asymptotically, as a function of time [20]. A negative maximal SLE indicates that the synchronous state, between the driver X and response system Y , is stable [33, 34, 20]. What this means is that the response system Y will synchronise to its driver X after a finite transient, irrespective of their initial conditions [19]. The same holds for two response systems that are sufficiently strongly coupled to a common driver. If the response systems are identical, complete synchronisation will occur after a finite transient, and the dynamics of the of the response systems will be indistinguishable, irrespective of their initial conditions.

We used the method outlined by Schiff et al. [19] to compute the largest SLEs of the response system Y , as a function of the coupling strength. The algorithm consists of the following steps:

1. Randomly select an initial point (x0, y0) on the attractor of the response system, and a

random n dimensional unit vector ˆe0, where n is the dimension of the subspace formed by

the response system (2 in our case).

2. Iteratively compute the following quantities for 1 ≤ i ≤ imax (where imax is the maximal

number of iterations):

ai= |δYY (yi−1, xi−1, c)ˆei−1| (14)

ˆ ei=

δYY (yi−1, xi−1, c)ˆei−1

ai

(15)

here c is the coupling strength, and δYY (yi−1, xi−1, c) is the Jacobian matrix of the response

system Y , with respect to the response system Y only (i.e. without the driver), evaluated at the i − 1th _{iterate of the orbit starting from (x}

0, y0) [19]. The magnitude of the matrix

vector product between the Jacobian matrix of g and the random unit vector ˆei−1is denoted

(21)

3. Finally, we approximated the largest sub-Lyapunov exponent as: λ ≈ 1 t t X i=1 log(ai). (16)

The Jacobian of the response system, with respect to the response system only, can be derived from the governing equations (8) and (9), and is given by:

δYY = " −cx(t) + (1 − c)2y(t) 0.3 1 0 # (17)

We varied the coupling strength c from 0 (no coupling) to 1 (maximal coupling) in 101 equally spaced steps, and estimated the maximal SLE using 104 _{map iterates. We repeated this process}

100 times, with different initial conditions [x0, u0, y0, v0], and took the average of these repetitions

as our final estimate of the maximal SLEs.

4.3.2 Phase synchronisation analysis

To study phase synchronisation in coupled chaotic systems, a typical approach is to estimate the instantaneous phase angle of the (oscillating) systems, and check whether a phase locking condition |nφ1− mφ2| < c is satisfied, where c is some constant [23, 1]. Rosenblum et al. [23] have noted

that the phase angle is not uniquely defined for chaotic oscillators, and that while there exist some cases where specific observables yield more intuitive definitions, in principle the phase angle of chaotic oscillators can be computed from any system observable.

A typical approach for phase angle estimation is through the Hilbert transform of the observed signal s(t): ˜ s(t) = π−1P.V. Z ∞ −∞ s(τ ) t − τdτ (18)

where P.V. denotes that the integral is evaluated in the sense of the Cauchy principal value [23]. Using ˜s(t), the complex valued analytic signal ψ(t) is then given by ψ(t) = s(t) + i˜s(t). Finally, the phase angle φ(t) of the observed signal s(t) is obtained, as with any complex function, by taking the inverse tangent of the ratio of the parts of ψ(t):

φ(t) = arctans(t)˜

s(t). (19)

We can now define the relative phase by taking the difference between the phases for the signals of the different systems: ϕ(t) = φX(t) − φY(t). From this phase difference we can now compute

the mean phase coherence R, defined by Mormann et al. [1], as follows:

R = 1 N N −1 X j=0 eiϕ(j∆t) (20)

where 1/∆t is the sampling rate of the time series. Using Euler’s formula we can turn equation (19) into the following equivalent expression:

R =      1 N N −1 X j=0 sin[ϕ(j∆t)]   2 +   1 N N −1 X j=0 cos[ϕ(j∆t)]   2   (21)

We computed the mean phase coherence, using equation (21), between the unidirectionally coupled driver and response variable (X,Y ), and the mutually uncoupled response variables (Y 1,Y 2) for 101 time series with coupling strengths ranging from 0 to 1. All time series consisted of 5 × 103

system iterations. We used the MATLAB built-in functions hilbert and angle to compute the Hilbert transform and compute the phase angle, respectively.

(22)

4.3.3 Coarse-grained information rates

Paluˇs et al. [2] have introduced the coarse-grained information rate (CIR), which measures the complexity of a system’s dynamics over time, and can be used to detect complete and generalized synchronisation from time series. The CIR is defined as:

i(X) = ∆τ τmax− τmin+ ∆τ τmax X τ =τmin I(X; Xτ), (22)

and the mutual CIR is defined as:

i(X, Y ) = 1 2τmax τmax;τ 6=0 X τ =−τmax I(X; Yτ), (23)

where τmin / max denote the minimal and maximal shift between the time series, respectively, and

∆τ denotes the time shift step size. The usual choice for τmin = ∆τ = 1 sample (i.e. step) [2].

Lastly, I(X; Y ) denotes the mutual information between variables X and Y , which is defined as:

I(X; Y ) = X x∈X X y∈Y pX,Y(x, y) log pX,Y(x, y) pX(x)pY(y) (24)

For the detection of synchrony in chaotic systems, Paluˇs et al. [2] have introduced the following relations:

i(X) = i(X, Y ) = i(Y ), (25) min(i(X), i(Y )) ≤ i(X, Y ) ≤ max(i(X), i(Y )) (26)

to determine whether a chaotic system is in either a complete or generalized synchronised state, respectively [2].

We computed the CIRs from the time series of all system variables X, Y 1 and Y 2. In addition, we computed the mutual CIRs for system pairs (X, Y ) and (Y 1, Y 2). All CIRs were computed for 101 time series with coupling strengths ranging from 0 to 1. All time series consisted of 5 × 103 system iterations. We computed the CIRs using a maximal shift τmax = 10 time steps, however

different values of τmaxproduced qualitatively identical results. We discretized the data using the

discretize function from the infotheo R package [27] using the equal frequencies algorithm. After discretization, we performed the mutual information computations using the function mutinfor-mation from the infotheo package.

(23)

5 Results

5.1 Analysis of the maximal sub-Lyapunov exponents (SLE)

Figure 7 shows the maximal SLE as a function of the coupling strength. The grey shaded areas (S1: c ∈ [0.45, 0.51], and S2: c ∈ [0.7, 1]) mark regions where the maximal SLE is negative, indicating that the response system has synchronised with the driver system [20].

Figure 7: Maximal sub-Lyapunov exponents (SLEs) of the response systems Y1 and Y2, as a

function of the coupling strength from X to Y . The grey sections S1 ( c ∈ [0.45, 0.51]) and S2 (c ∈ [0.7, 1]) show regions where the maximal SLEs are negative. This is an indication of complete synchronisation of the response systems to the dynamics of the driver system.

5.2 CCM analysis

Figures 8 and 9 show two examples of CCM results for systems with weak (c = 0.03) and moderate coupling (c = 0.4), respectively. Figures 8A and 9A show the CCM results for the true coupling relation. The direction of the coupling can be correctly inferred in both cases, as convergence is observed for the cross map from the response to the driver variable (red line), but not for the cross map in the opposite direction (blue line). The convergence is authentic, as no convergence is observed for the cross mapping using random neighbours (figures 8C and 9C). Figures 8B and 9B show the results of the cross mapping between the mutually uncoupled response variables. A CCM analysis should (ideally) not show any convergence for these cases, which it doesn’t. The CCM score is non-zero for the moderate coupling case (figure 9B), but this is in itself not a sufficient condition for Sugihara causality, as no convergence is observed. The scores obtained for random cross mapping between Y 1 and Y 2 (figure 9D) are higher than for the random cross mapping between X and Y (figure 9C). This is due to the increased similarity between signals Y 1 and Y 2, compared to X and Y . It is important to note here that, while both response systems are differently initialized, and their corresponding time series are therefore different, both systems have the same underlying attractor, which means that, in the limit of infinite data, CCM would show convergence if we are allowed to shift the time series by an arbitrary amount (this point will be further addressed in the discussion).

(24)

Figure 8: Example CCM results for a system of coupled H´enon maps with weak coupling (coupling strength = 0.03). A: CCM scores for the cross maps between the driver X and response system Y , as a function of the library size. The red curve corresponds with the cross map X|MY, where the

shadow manifold of Y is used to make predictions about X. This cross mapping corresponds with the true causal relation X → Y , which is reflected in the upward trend of the curve as the library size increases (convergence). The blue curve corresponds with the cross mapping in the other direction, and shows no convergence, as expected. B: same as A but for the response systems Y 1 and Y 2, which received the same input from the driver system X. Since there is no direct causal relation between these two systems, there is no convergence, as expected. C: same as A but the cross mappings were performed using random neighbours instead of nearest neighbours. D: same as C but for the response systems Y 1 and Y 2.

(25)

Figure 9: Example CCM results for a system of coupled H´enon maps with moderate coupling (coupling strength = 0.4). A: CCM scores for the cross maps between the driver X and response system Y , as a function of the library size. The red curve corresponds with the cross map X|MY,

where the shadow manifold of Y is used to make predictions about X. This cross mapping corresponds with the true causal relation X → Y , which is reflected in the upward trend of the curve as the library size increases (convergence). The blue curve corresponds with the cross mapping in the other direction, and shows no convergence, as expected. B: same as A but for the response systems Y 1 and Y 2, which received the same input from the driver system X. Since there is no direct causal relation between these two systems, there is no convergence, as expected. C: same as A but the cross mappings were performed using random neighbours instead of nearest neighbours. D: same as C but for the response systems Y 1 and Y 2.

Figure 10 shows the results of the CCM analysis, for all coupling strengths. The results are superimposed over the maximal SLEs, to allow for an easy comparison between the CCM results and the degree of synchronisation between the driver and response systems. Figures 10A and 10B show the maximal CCM scores (obtained for the largest library size of 103) for the cross mapping between X and Y , and Y 1 and Y 2, respectively. The maximal CCM scores for X|MY are non-zero

for all non-zero coupling strengths, and become 1 in region S2 (figure 10A). The CCM analysis provided the correct directionality of the driver-response relation for 80% of coupling strengths that did not lead to complete synchronisation (see figure 7): convergence was observed for X|MY

and not for Y |MX (figure 10C). The CCM analysis consistently failed to produce the correct

results when the coupling strength exceeded 0.7, where the ∆CCM values were approximately equal for the cross mappings in both directions (figure 10C). However, using the outcomes of the delay CCM analysis, the accuracy of the causal inference increased to 93.8%. Out of a total of 65 inferences for coupling strengths outside regions S1 and S2, only 4 led to the false detection of a bidirectional coupling relation. For coupling strengths 0.1, 0.2 and 0.69 the best cross map lag of Y |MX was negative. For a coupling strength of 0 (no coupling), the delay CCM analysis

(26)

The CCM results between the response variables Y 1 and Y 2 (figure 10B and 10D) within region S2, are similar to the results obtained for the cross mapping between X and Y (figure 10A and 10C). This is expected, since the time series of X, Y 1 and Y 2 are (almost) identical for the coupling strengths in region S2. High maximal CCM scores for the response systems were also observed in region S1, and between region S1 and S2 (figure 10B). In total, 33.9% of ∆CCM scores were significant (figure 10) for coupling strengths that did not produce complete synchrony, meaning that we would incorrectly infer a uni- or bidirectional coupling relation between systems Y 1 and Y 2. Using the delay CCM analysis and false causality conditions (F1) and (F2) (section 4.2.5), we find a final false positive rate of 21.5%. This signifies an improvement of 36.6%, over the 33.9% false positive scores obtained without the delay CCM method.

Figure 10: CCM scores as a function of the coupling strength from X to Y , overlayed on figure 7. The maximal CCM scores were obtained for a library size of 2 × 103, and the ∆CCM scores were obtained by computing the difference between the maximal CCM scores and the CCM scores obtained for the smallest library size of 10 points. The coupling strength was varied from 0 to 1 in 101 steps of 0.1. The CCM scores were obtained by averaging over 10 repetitions, where the libraries were randomly sampled with replacement. A: maximal CCM scores for the cross maps between the driver X and response system Y . The red curve corresponds with the cross maps X|MY, which quantifies the true causal relation X → Y . The blue curve corresponds with the

cross map in the opposite direction, and should not produce positive results (i.e. convergence). B: same as A but for the response systems Y 1 and Y 2, which received the same input from the driver system X. C: same as A but for the ∆CCM scores. D: same as C but for the response systems Y 1 and Y 2.

(27)

5.3 Delay CCM analysis

Figure 11 shows the results for the delay CCM analysis between X and Y . Figures 11A and 11C correspond with the maximal and ∆CCM scores for X|MY, respectively. While figures 11B and

11D correspond with the same scores for Y |MX. Figure 12 shows similar results for the cross

mappings between the response systems Y 1 and Y 2. The lighter colors indicate a higher CCM score. These results are similar to the results for the normal CCM analysis, in that high CCM scores and convergence are observed for the true coupling direction, starting from weak coupling, until complete synchronisation results in the false inference of a bidirectional relation between X and Y (figure 11). In the synchronous region S2, very high maximal CCM scores are obtained for cross-map lags ranging from -5 to +5 (figure 11A). Since the X and Y signals are identical to each other in S2, due to complete synchronisation, cross mapping here amounts to self prediction. The temporal extent of high maximal cross mapping scores, thus indicate the backward and forward predictability of the H´enon map. The low ∆CCM scores in region S2 for cross-map lags in the interval [-3, 1] are a result of high CCM scores obtained for the minimal library size, rather than a lack of convergence. This is expected, since these temporally close points have similar delay embedding vectors (i.e. at least 1 and at most E-1 of the vector elements are identical), which means that the geometrical relation between the putative nearest neighbour points (sampled for the simplex projection) and the target point, is very similar for the reconstructed attractors under these small temporal shifts (cross-map lags). However, as the cross-map lags increase, the geometrical relation is destroyed (but the topological relation is preserved) and the cross mapping becomes less accurate for small library sizes.

Figure 11: CCM scores between X and Y , as a function of the coupling strength from X to Y , and the cross-map lag. The maximal CCM scores were obtained for a library size of 2 × 103_{, and}

the ∆CCM scores were obtained by computing the difference between the maximal CCM scores and the CCM scores obtained for the smallest library size of 10 points. The coupling strength was varied from 0 to 1 in 101 steps of 0.1, while the cross map lag was varied from -10 to +10 in steps of 1. A: maximal CCM scores for the cross maps X|MY, which quantify the true causal relation

X → Y . B: same as A but for the cross maps Y |MX. C: same as A but for the ∆CCM scores.

(28)

Figure 12: CCM scores between Y 1 and Y 2, as a function of the coupling strength from X to Y , and the cross-map lag. The maximal CCM scores were obtained for a library size of 2 × 103, and the ∆CCM scores were obtained by computing the difference between the maximal CCM scores and the CCM scores obtained for the smallest library size of 10 points. The coupling strength was varied from 0 to 1 in 101 steps of 0.1, while the cross map lag was varied from -10 to +10 in steps of 1. A: maximal CCM scores for the cross maps Y 1|MY 2. Since there is no causal relation

between Y 1 and Y 2, the CCM analysis should not give positive results (i.e. convergence). B: same as A but for the cross maps Y 2|MY 1. C: same as A but for the ∆CCM scores. D: same as

C but for the cross maps Y 2|MY 1.

Figure 13 shows the cross map lags that produce the highest CCM scores (i.e. best cross map lag), per coupling strength. The black circles mark cross maps with unidirectionally significant convergence, whereas the white circles indicate significant convergence in both cross mapping directions. Figure 13A shows the results for the cross maps between the driver and response systems. The X|MY cross maps are all maximized for negative cross map lags, in accordance with

the claim by Ye et al. [18] that the best CCM scores are obtained for shifts that are proportional to the true coupling lag. The negative cross map lags for X|MY indicate that the dynamics of

time series X precedes the dynamics of Y , which is correct since the value of X at time t, is used to compute the value of Y at time t + 1. The majority of cross maps Y |MX, achieve the best scores

for positive cross map lags. This is again in accordance with the claim by Ye et al. [18], however the estimation for these cross maps is less accurate, as indicated by the increased variance of the optimal cross map lags. Figure 13B shows the results for the cross maps between the response systems. Since the results for the cross maps between X and Y , and Y 1 and Y 2 are the same in region S2, we cannot tell the difference between the two in the case of (very) strong coupling. However, in region S1 and the region between S1 and S2, there are a relatively large number of coupling strengths for which we can distinguish the true coupling relation between X and Y from the false coupling relation between Y 1 and Y 2. More specifically, the negative best cross map lags in figure 13B indicate that Y 1 and Y 2 are synchronised for almost all coupling strengths c ≥ 0.45. Thus, the delay CCM method is potentially helpful to distinguish a direct coupling relation from false coupling relations as a result of common input, as long as the coupling strength between the driver and response systems is not sufficient to cause complete synchronisation. In the following

Inferring Causality in Synchronised Chaotic Systems using Convergent Cross Mapping