Causality Analysis: Identifying the Leading Element in a Coupled Dynamical System

Physical systems with time-varying internal couplings are abundant in nature. While the full governing equations of these systems are typically unknown due to insufficient understanding of their internal mechanisms, there is often interest in determining the leading element. Here, the leading element is defined as the sub-system with the largest coupling coefficient averaged over a selected time span. Previously, the Convergent Cross Mapping (CCM) method has been employed to determine causality and dominant component in weakly coupled systems with constant coupling coefficients. In this study, CCM is applied to a pair of coupled Lorenz systems with time-varying coupling coefficients, exhibiting switching between dominant sub-systems in different periods. Four sets of numerical experiments are carried out. The first three cases consist of different coupling coefficient schemes: I) Periodic–constant, II) Normal, and III) Mixed Normal/Non-normal. In case IV, numerical experiment of cases II and III are repeated with imposed temporal uncertainties as well as additive normal noise. Our results show that, through detecting directional interactions, CCM identifies the leading sub-system in all cases except when the average coupling coefficients are approximately equal, i.e., when the dominant sub-system is not well defined.


Introduction
Identifying the leading [sub-]system from a pair of coupled dynamical systems using only time-series is challenging when nothing or little is known about the underlying dynamics. The definitive approach to detect causal relationships between components of a system is to fully identify the underlying physical mechanisms and governing equations. However, we only have partial knowledge about the internal physical mechanisms in most cases and must resort to observed data to establish the existing causal relationships in such cases.
In linear systems, lead time of the driver (or equivalently, latency of the response) may indicate a causal relationship, which hence could be identified by [linear] lag-correlation analysis. In a nonlinear system, however, there may be no persistent lead-lag between the two signals, for example, due to feedback loops between the variables. Therefore, a linear lag-correlation analysis may not reliably determine the correct causal relationship in a nonlinear system [1]. For example, Fig 1 shows lead-lag switching between two coupled variables. If one considers short intervals of the time-series and takes lead time and strong correlation as the indicators of causation, an incorrect conclusion about the causal relationship between the variables could be made. In such circumstances, one may analyze the signal phases to establish temporal precedence and exploit the phase slope index to estimate the flow direction of information flux [2].
The concept of information transfer has been used as an indication of causality. A practical measure of information transfer is the transfer entropy which distinguishes directional communications between variables of a system [3]. This probabilistic measure has been used in a variety of fields such as neuroscience [4,5], the study of chemical processes [6], and cellular automata [7]. A more general form of the transfer entropy that has an augmented time-delay parameter can detect the propagation time in addition to the asymmetric information transfer between observed signals [8].
Before establishing the general notion of information transfer (and transfer entropy), a practical definition of causality was proposed by Wiener [9], and later adopted and formalized in terms of linear autoregression by Granger [10,11]. In the context of information transfer, Granger causality is shown to be equivalent to transfer entropy for Gaussian variables [12]. Although the standard [linear] Wiener-Granger method was first introduced for economic systems, the method and its extended nonlinear versions gained widespread use in other fields such as neuroscience and finance [13][14][15][16][17][18]. According to the Wiener-Granger definition of causality, a candidate driver is among the drivers of a response signal if the response prediction error increases significantly by removing the candidate driver data from the universe of all drivers. More precisely, given sets of interdependent variables X and Y, it is said that "X causes Y" if, in an appropriate statistical sense, X assists in predicting the future of Y beyond the degree to which Y already predicts its own future [12]. The Granger method is a forward method as it uses the driver data to predict the response. An important characteristic of the Granger method is that it requires the signals to be separable and have non-zero entropy rates [19]. In a separable system, it is possible to separate a candidate driver's data from other factors, thus enabling prediction to be conducted using data sets including and excluding the candidate driver. Separability is a restrictive conditions since many observed signals of interest are from deterministic nonlinear systems with feedbacks between state variables, resulting in mixed information from different sources.
Although information transfer and Granger causality provide valuable statistical information about the observed signals and their asymmetric connectivity, "efficient" causal relationship, as defined by Lizier and Prokopenko, between variables of a system can ultimately be identified by interventional methods, i.e., perturbing the candidate cause variable to investigate its direct influence on the response variable [20,21]. Information flow was proposed as a quantitative measure of intervention. This Bayesian probabilistic measure quantifies the distribution of a response variable as a result of "imposing" conditions. However, there are some practical limitations for applying this measure and detecting causal relationships in realistic systems. For example, one may need to know about the structure of the causal links in a network or the underlying rules of the causal interactions. But understanding such structures, not known a priori, is the purpose of a causality study. The back-door adjustment criterion was proposed to solve this problem under certain conditions [20][21][22].
Sugihara et al. [23] investigated the problem of causality from a new perspective, proposing the Convergent Cross Mapping (CCM) method for deterministic nonlinear systems with smooth manifolds. The fundamental idea of the CCM method is that if Y is causally influenced by X, then Y has signatures of X such that the historical record of Y can reliably estimate the state of X. Therefore, a better estimate of the driver variable, X, shows stronger causal influence on the response variable, Y.
It is important to note that: (i) unlike Granger causality that uses prediction as its fundamental basis, the CCM method is not based on the prediction of a variable; (ii) the CCM method does not apply Bayesian probability, which is the core concept of the transfer entropy and information flow measure; and (iii) the CCM method is not an interventional method. Therefore, it does not perturb the system to identify the micro-level causal effects in a system. It measures the correlation between the reconstructed and recovered manifolds of the observed signals (see x1).
Sugihara et al. [23] showed that CCM can identify unidirectional and bidirectional causation, and dominant driver, in weakly coupled nonlinear systems with constant coupling coefficients. They also considered [in Supplementary Materials] examples of systems with asymmetrical couplings, external forcing, and time delays [lagged influences]. They presented successful applications of the CCM method in ecology, biology, and geoscience [23][24][25]. Our study is motivated by the need to identify the dominant constituent in systems that are speculated to have time-varying interconnections with switching between the dominant elements in different periods. A challenging application is to identify the dominant variables of the global climate system from geophysical records of greenhouse gases concentrations and temperature proxy [26][27][28]. For this purpose, we consider coupled systems that (i) have variable coupling parameters in sequential periods, and (ii) the larger coupling coefficient is not fixed for a  specific sub-system. Thus, we may observe switching between the dominant sub-systems during successive periods. We investigate whether CCM, through detecting directional interactions, can identify the leading sub-system from a pair of coupled sub-systems under different conditions. In this paper, the leading sub-system is defined as the system with a larger average coupling coefficient over specified intervals of the time-series assuming that the coupled systems have the same time scales. To do that, we perform numerical experiments with different [linear] coupling schemes of Lorenz system [29] with asymptotic synchronization [30]. We consider periodic-constant (case I), normally distributed (case II), and mixed normally/nonnormally distributed (case III) coupling coefficients. We also investigate the CCM results in the presence of temporal uncertainties and additive noise (case IV). Our results indicate that CCM can identify the leading sub-system, except when the average coupling coefficients are approximately equal, i.e., when the leading sub-system is not well-defined.
This article is organized as follows. Section x1 covers mathematical details of the CCM method. Section x2 describes the coupling schemes. In section x3, we show results of the numerical experiments. In section x4, we summarize and discuss our findings.

Method
In the CCM method, it is assumed that a response signal has signatures from its driver so that the approximate behavior of the driver can be estimated (or recovered) from the response signal. Thus, a better estimate of the driver shows stronger causal influence on the response variable. To implement the CCM method, we use the concept of one-to-one mapping between the original smooth manifold of the [full] system and the compact reconstructed phase spaces (shadowing manifolds) of the observed signals [31,32]. If the two observed signals belong to the same dynamical system, a one-to-one mapping between the two reconstructed phase spaces could be established by considering the original manifold as an interface. If the dimension of the reconstructed phase spaces are selected based on the criteria of the false neighborhood method [33], arbitrary nearby points on the driver's reconstructed phase space map to the nearby points on the original attractor. Because of the causal influence of the driver, the nearby points on the original manifold stay close on the reconstructed phase space of the response signal. The stronger causal influence of the driver on the response, the closer the mapped points of the response signal on its shadowing manifold [23].
We assume, without loss of generality, that x(t) is the driver and y(t) is the response. The time-series x(t) and y(t) are sampled at equally spaced time intervals Δt, ( where n = 1, 2, Á Á Á, N, and N is the total number of data points in each time-series. We set t 0 = 0 to simplify the notation. We use [x(t), y(t)] and [x n , y n ] notations interchangeably. CCM uses the time-delay coordinates to reconstruct the phase spaces of the L-point windows from the time-series x(t) and y(t) [34]. The L-point windows are defined as ( for i = 1 to i = N + 1 − L and L min < = L < = L max . We define L min and L max below. The L-point windows sweep the entire time-series x(t) and y(t) (see Fig 2a).
We use the average mutual information measure, I, a nonlinear generalization of the correlation function, to select the time lags τ as an integer multiple of Δt. This measure shows the average amount of information about a signal at t + t 0 , i.e., s(t + t 0 ), when s(t) is observed where P(s(t), s(t + t 0 )) is the joint probability of measuring s(t) and s(t + t 0 ). We choose τ where the first minimum of I(t 0 ) occurs [31,35,36]. We employ the false neighborhood method to determine the embedding dimension, E. The correct embedding dimension prevents self overlapping of the projected manifold. Therefore, "if E is qualified as an embedding dimension by the embedding theorem [31,32], then any two points which stay close in the E-dimensional reconstructed space will be still close in the (E + 1)-dimensional reconstructed space. Such a pair of points are called true neighbors, otherwise, they are called false neighbors" [33].
By these two parameters, i.e., τ and E, we generate the time-delay coordinates as well as the reconstructed phase space resembling the original manifold of the attractor. The E-dimensional time-delay coordinate vectors corresponding to W i;L x and W i;L y are shown as m i;L x ðtÞ ¼ ½xðt À ðE À 1ÞtÞ; Á Á Á ; xðt À tÞ; xðtÞ m i;L y ðtÞ ¼ ½yðt À ðE À 1ÞtÞ; Á Á Á ; yðt À tÞ; yðtÞ Δt. Limits of t are chosen such that the first/last component of m i, L vectors corresponds to the first/last point of W i, L windows. The reconstructed phase spaces are After reconstructing the phase spaces, sufficient numbers of nearest neighbor points are selected for each E-dimensional vector in M i;L y , referred to as a Y-central point, Y c (see Fig 2b). The neighbor points are denoted by Y j , ordered by their distance, d j , to Y c , from nearest to farthest. This means Y 1 is the nearest neighbor to Y c . The number of selected neighbor points should be equal or larger than E + 1 so the simplex method can represent the dynamics on an E-dimensional manifold [37]. This requirement specifies L min = (E + 2) + (E − 1)(τ/Δt). We use the corresponding contemporaneous points of the Y j in M x , denoted by X j , to calculate the position of a single point that represents them,X . We define the exponentially decaying weights based on the distances between the Y c and its neighbor points, d j , as where k Á k is the Euclidean norm. The position of the representative point in M x is calculated asX The set of these recovered points,X È É , is the recovered phase space. If x(t) and y(t) are dynamically coupled,X È É should be strongly correlated to M x = {X c }, where X c is the point corresponding to Y c (see the star and black filled circle in Fig 2b). In addition, as L increases, the density of the points in the reconstructed phase space increases, hence the distances between the nearest neighbors shrink. Thus,X converges to a vicinity of X c that becomes smaller as the causal effect of the driver increases. Therefore, the CCM coefficient, ρ(L), as a measure of causal relationship between the [candidate] driver and the response signal is defined as the correlation coefficient betweenX and X c , averaged over all possible L-point windows, If there is a causal relationship between the signals, ρ converges to a constant equal or less than one. Note that h Á i L indicates averaging over all L-point windows, andX j M y emphasizes thatX is estimated given M y . If we study the opposite roles, i.e., y(t) as the driver and x(t) as the response, the CCM coefficient would be shown by r Y c ;Ŷ jM x . We simplify this notation by setting It is important to note that: (i) the CCM method, through detecting directional dynamical interactions, can investigate causal relationship over any window of the time-series that contains sufficient information about the attractor of the dynamical system. In the case of nonstationary signals, the CCM results during the considered time interval are presumably valid if the original manifold and the reconstructed phase spaces do not have significant changes, i.e., if the near-stationary condition holds. We do not discuss the application of CCM to nonstationary signals in this paper. (ii) CCM only considers the past data and does not attempt to predict the signals. Therefore, CCM could be robustly applied to chaotic systems, where predictions may rapidly diverge from the true/observed states. (iii) The identified causal relationship by the CCM method is not exclusive, i.e., the identified driver might be one of the many drivers of the system. Therefore, in systems with multiple variables, CCM can be applied and repeated for different combinations of the candidate driver and response signals. (iv) Our results (not shown here) show that the CCM method is not sensitive to increases in the embedding dimension if it is chosen by the false neighborhood method. By increasing the embedding dimension beyond what is prescribed by the false neighborhood method, computation cost increases because the search process for the nearest neighbors should be performed in a larger space dimension. However, the CCM results remain almost unchanged because the extra dimensions do not add new information. In addition, the CCM results are not sensitive to small changes of τ if the embedding parameter remains in the vicinity of the first minimum of the average mutual information measure.

Experiment design with coupled Lorenz systems
In order to study the applicability of the CCM method to [strongly] coupled systems with variable coupling coefficients, we apply it to a pair of identical Lorenz systems [29], L X and L Y , coupled with different coupling schemes. We choose to study synchronized dynamical systems, in which the coupling is canonically between two identical [sub-]systems coupled simultaneously with a linear coupling term where the dominant element is clearly identified by its larger coupling coefficient. Study of synchronized dynamical systems has expanded beyond the canonical case to included "generalized synchronization" [38], in which the two systems differ. Other work in synchronized dynamical systems has studied "lagged synchronization" [39]. Since this application of the CCM method is a new, we choose to focus on the canonical representation of the problem and allow for similar extensions in future work to address these more complex variations of the general problem.
The general form of L X and L Y is given by Eq (9).
The parameters of the Lorenz-63 systems are (σ, ρ, β) = (16, 45.92, 4) throughout this paper. μ and η are the coupling coefficients, and k μ and k η are the binary on-off switches. Variations of the coupling coefficients, activation of the on-off switches, and duration of successive periods are discussed in section x3.
We investigate whether CCM can distinguish the leading system by detecting the difference between the CCM coefficients rX jY and rŶ jX . We probe one signal from each system, x 1 (t) and y 1 (t) in the following experiments. When the two systems are strongly correlated, for example, due to a strong feedback loop, rX jY and rŶ jX are close to each other for large L's. (This is expected, unless the two systems have different amplitudes and time scales, as in the case of the ocean-atmosphere models [40].) Thus, one must investigate small differences between rX jY and rŶ jX .
A brief description of the four experiments is as follows. In case I, L X and L Y are coupled by a periodic on-off switching mechanism and constant coupling coefficients during all periods. We cover 100 different combinations of coupling coefficients. In case II, both coupling coefficients in all periods are normal random variables with specified means and standard deviations.
In case III, one of the coupling coefficients is a normal random variable and the other one is from a Weibull (skewed non-normal) distribution [41]. In case IV, we consider the original pairs of signals from cases II and III and impose relative temporal shifts (representing temporal uncertainties) as well as additive Gaussian noise. In all experiments, we observe a strong correlation between the probed signals x 1 (t) and y 1 (t) due to bidirectional couplings in system Eq (9).
As discussed in x1, ρ(L) converges to a constant by increasing the length of the L-point windows, provided that the two signals are dynamically connected. We test convergence of ρ(L) over the range L = 25 to 500. We observe negligible variations in the CCM coefficients for L ! 300. We choose L = 500 to report all the CCM coefficients. This choice of L is equivalent to 0.5 time unit since we sample the signals 1000 times per non-dimensional time unit.

Experiment results
Below we present the results of four sets of numerical experiments with the system Eq (9).
3.1 Case I: Coupled Lorenz systems, periodic-constant coupling coefficients. In this case, coupling coefficients μ and η assume constant integers 1 to 10, resulting in 100 combinations, during 10 consecutive periods of the time-series, each spanning 10 non-dimensional time units. The switching mechanism is controlled by k μ and k η , such that in odd periods k μ = 0 and k η = 1 and in even periods k μ = 1 and k η = 0 as described in Eq (10). We observe similar patterns in rŶ jX and rX jY (not shown here) such as vertical and horizontal bands. For example, we see a blue horizontal band near μ = 1 in Fig 3a, showing that the influence of L Y on L x is small. Therefore,Ŷ È É (the recovered phase space of L Y from L X ), is poorly correlated to M y (the reconstructed phase space of L Y ), hence the small values for rŶ jX on the band. Another observed pattern is the high values of ρ in the regions with the strongest coupling between the two sub-systems (high μ and η), showing the large mutual influence between the two systems.
To address the main question of this study we investigate rŶ jX À rX jY . We expect that if L Y is the leading sub-system, then rŶ jX is larger than rX jY and similarly, rX jY is larger than rŶ jX if L X is the dominant sub-system. Referring to Fig 3b, we observe that rŶ jX > rX jY when μ > η (upper diagonal of Fig 3b) and rX jY > rŶ jX when η > μ (lower diagonal of Fig 3b).
This numerical experiment with 100 combinations of coupling coefficients shows that even when ρ values are close to each other (rŶ jX % rX jY ), CCM can capture the leading system through small differences between the CCM coefficients, except for μ % η where neither system is leading.

Case II: Coupled Lorenz systems, normally distributed coupling coefficients
In section x3.1, we studied CCM's capability to identify the leading system in coupled systems with periodic constant coupling coefficients. A more general and realistic case is a system with randomly distributed coupling coefficients and with random lengths of time periods. In this section, we choose the coupling coefficients from normal distributions. Note that these coefficients are constant during each period. In Eq (9), we set η = η N and μ = μ N , where subscript N stands for normal random variables. We set the mean and standard deviation of the coupling coefficients such that the ratio of the averaged coupling coefficients Z N = m N covers the approximate range (0.5,5) and we have enough data points both above and below 1. The switching coefficients, k η and k μ , are kept on for all periods in order to have bidirectional communication between the two sub-systems. We also choose the duration of each period, T i , from a normal distribution as described in Eq (11). We also observe a close-to-linear relationship between Z N = m N and Δρ in the considered range of the coupling coefficients.
We extend this analysis using the Monte Carlo method. We generate 1000 independent realizations of the coupled system for each α 2 int [1:10] (thus, we have 10,000 individual data points) and calculate the difference between the CCM coefficients at L = 500. As we observe in Fig 4b, Δρ is negative for almost all realizations when ð Z N = m N Þ < 1 and positive when ð Z N = m N Þ > 1. Note that the distribution of the points on the two sides of ð Z N = m N Þ ¼ 1 depends on the choices of η N and μ N in Eq (11). Fig 4a and 4b shows that CCM can determine the leading system when the coupling coefficients are normally distributed except when Z N % m N , i.e., when the leading sub-system is not well-defined.
Because ρ(L) is obtained by averaging over all possible L-point windows (see Eq (8)), one might be concerned whether the difference between ρ values, i.e., Dr ¼ rX jY À rŶ jX , is statistically significant and meaningful. To address this concern, we perform a non-parametric onesample Kolmogorov-Smirnov test for all data points of Fig 4a. Results of this test, for rX jY and rŶ jX , fail to reject the null hypothesis that the data in the ρ vectors come from a standard normal distribution at 0.001 significance level. Next, a group t-test shows that for all the data points the obtained t-values exceed the required minimum corresponding to a critical p-value of 0.01. Therefore, we can conclude that the differences, Δρ, are significant at 0.01 level. For example, for the smallest absolute value of Δρ in Fig 4a, the calculated t-value equals 3.16 which is larger than 2.58 corresponding to 0.01 critical p-value. These numerical values correspond to the sample size of the L-point windows at L = 500. Therefore, we can reject the null hypothesis that there is no significant difference between the mean values of ρ.

Case III: Coupled Lorenz systems, normally and non-normally distributed coupling coefficients
For numerical experiments x3.1 and x3.2, we considered constant and normally distributed coupling coefficients, respectively. In case II, the difference between the two coupling coefficients is also a normal variable because both coefficients are normal variables. In case III, to further generalize our study, we break the symmetry of the coupling coefficients by choosing one of them from a normal and the other from a non-normal random distribution. By this selection, the difference between the coupling coefficients is not a normal variable. We set η = η N as in x3.2 but choose μ = μ Wb from Weibull distributions [41]. As before, each time-series consists of twenty periods with normally distributed lengths. (During each period, the coupling coefficients remain constant.) The current setting is summarized in Eq (12).
2Þ j a 2 int½1 : 11 8 > > > > < > > > > : Parameters of the Weibull distribution are chosen in order to have sufficient values of Z N = m Wb both below and above one. The same reasoning also applies to setting Eq (11) of experiment case II. Fig 5 shows the asymmetric Weibull probability density function for the parameters of Eq (12).
Similar to cases I and II, we investigate Dr ¼ rX jY À rŶ jX as a function of Z N = m Wb at L = 500 (see Fig 6a). Similar to Fig 4a, we observe that Δρ is negative for ð Z N = m Wb Þ < 1 and is positive for ð Z N = m Wb Þ > 1. This result supports our expectation about the ability of the CCM method to distinguish the leading sub-system. We employ the Monte Carlo method to generate 1000 independent realizations of the coupled Lorenz systems with setting Eq (12) which results in a total of 11,000 individual data points. Then we calculate Δρ for each realization at L = 500, shown in Fig 6b. We again observe that Δρ is negative for ð Z N = m Wb Þ < 1 and positive for ð Z N = m Wb Þ > 1, except for a small fraction of points around Z N = m Wb % 1 where the leading system is not well-defined. Therefore, CCM can distinguish the leading sub-system in a coupled system with normal/non-normal coefficients, except when the two sub-systems have almost equal coupling coefficients. We also note that in both cases II and III, detection of the leading sub-system is more reliable for larger values of ð Z= mÞ.
In contrast to case II, we observe a nonlinear relationship between Δρ and Z N = m Wb , i.e., the rate of change of Δρ decreases rapidly as Z N = m Wb increases. Thus, we see a faster saturation in Fig 6 compared to Fig 4. Also, non-symmetric distribution of the data points around the vertical line Z N = m Wb ¼ 1 in Fig 6b is due to the range of Z N = m Wb according to Eq (12), which has a longer tail on the right hand side ( Z N = m Wb > 1). If we plot ÀDr ¼ rŶ jX À rX jY vs. m Wb = Z N , we would see a similar trend as in Fig 6b, although some ranges are non-overlapping due to the asymmetric distribution of the coupling coefficients. Similar to case II, we repeat the non-parametric one-sample Kolmogorov-Smirnov test and t-test for Δρ values of Fig 6a. In this case, except one data point that lies closest to the vertical axis Z N = m Wb ¼ 1, all other data points have t-values larger than the required minimum corresponding to a critical p-value of 0.01. Therefore, we can reject the null hypothesis and conclude that Δρ is significant at 0.01 level.

Case IV: Coupled Lorenz systems, temporally shifted and noisy signals
In this section, we investigate the ability of the CCM method for detecting the leading system in the presence of chronological uncertainties (temporal uncertainty between the two timeseries) as well as additive Gaussian white noise. We simulate the chronological uncertainties by relative shifting of the two time-series. We take the signals from cases II and III and apply relative shifts of 2.5 and 5% of the full length of the time-series.
For the additive Gaussian white noise, we consider a normal random variable with zero mean and standard deviation set at 5 and 10% level of the standard deviation of the original

Summary
There is often interest in determining causal relationship between variables of a physical system. Establishing a method for causality analysis is hence of paramount importance [22]. Linear lead-lag analysis is commonly applied to underpin causal relationship when there are some evidences about the internal mechanisms of the system. But quite frequently, two signals show inconsistent lead-lag behavior in different time windows due to nonlinearities of the system. Therefore, lead-lag analysis cannot conclusively determine causality over the full time span of observed signals. Introducing Granger causality [10] was a turning point in the field of causality analysis. Based on Wiener's definition of causality, Granger proposed a practical method that is founded on the idea of predicting a variable both with and without a candidate driver. If forecast is significantly improved when the information of the candidate driver is included in the set of predictors, the Granger method concludes that the candidate signal is a driver. Transfer entropy was introduced as a measure of directional communication between elements of a system [3]. Later, it was shown that Granger causality and transfer entropy are equivalent for Gaussian variables [12]. Although information transfer measures provide important understanding about the directional interconnections in a system, they do not identify efficient causal relationships. It is shown that interventional methods and information flow can identify micro-level causal relationships [21]. Sugihara et al. [23] proposed a new idea for causality analysis, applicable to deterministic nonlinear systems, based on cross convergent mapping between shadowing manifolds. They successfully applied the CCM method to analyze causality in weakly coupled systems with constant coupling coefficients. This study was inspired by identifying the leading element in systems that are speculated to have time-varying internal connections with probable change of dominant sub-system in different periods, for example, the earth system with its many interconnected sub-systems. For the first time, we addressed applicability of the CCM method to coupled systems with timevarying coupling coefficients and switching between dominant elements in different periods. We conducted numerical experiments with I) periodic-constant, II) normal, and III) mixed normal and non-normal coupling coefficients. In experiment IV, we imposed temporal uncertainties and additive noise to the observed time-series of cases II and III. We investigated whether the CCM method can identify the leading sub-system that has a larger average coupling coefficient over the entire span of the time-series.
Our main conclusions are: 1. If the averaged coupling coefficients are not approximately equal, i.e., Z≉ m, and a leading system exists, then the CCM coefficient of the leading system is significantly larger than the CCM coefficient of the system with a smaller average coupling coefficient.
2. If the ratio of the average coupling coefficients is close to one ð Z= m % 1Þ, the leading system is not well-defined and the CCM method is not applicable.
3. The CCM results are quite robust to temporal uncertainties and moderate levels of additive Gaussian noise.
4. For normally distributed coupling coefficients (case II), a close-to-linear relationship between Δρ and the ratio of the average coupling coefficients is observed (in the range of the selected coefficients).
5. For mixed normally and non-normally distributed coupling coefficients (case III), a [saturating] nonlinear relationship between Δρ and the ratio of the average coupling coefficients is observed (in the range of the selected coefficients).
According to these observations, we conclude that when the ratio of the average coupling coefficients is not close to one, the CCM method can detect the leading sub-system in a set of two coupled systems with time-varying coupling coefficients-even in the presence of chronological uncertainties or additive Gaussian noise.
There are still questions regarding applicability of CCM to systems with time delays [lagged influences], different time scales, different embedding dimensions, non-identical sub-systems, nonlinear influences, and non-smooth manifolds. Application of CCM to high-dimensional systems and non-stationary signals demands future studies too. Sufficiency of the number of observed data points for a reliable CCM analysis is another question that remains open. All of these open questions call for potential extensions of our study.
We anticipate the CCM method can be employed to study causal relationships between variables of systems such as those in atmospheric science, biology, ecology, epidemiology, and sociology.