Lead/Lag directionality is not generally equivalent to causality in nonlinear systems: Comparison of phase slope index and conditional mutual information

Applications of causal techniques to neural time series have increased extensively over last decades, including a wide and diverse family of methods focusing on electroencephalogram (EEG) analysis. Besides connectivity inferred in defined frequency bands, there is a growing interest in the analysis of cross-frequency interactions, in particular phase and amplitude coupling and directionality. Some studies show contradicting results of coupling directionality from high frequency to low frequency signal components, in spite of generally considered modulation of a high-frequency amplitude by a low-frequency phase. We have compared two widely used methods to estimate the directionality in cross frequency coupling: conditional mutual information (CMI) and phase slope index (PSI) . The latter, applied to infer cross-frequency phase–amplitude directionality from animal intracranial recordings, gives opposite results when comparing to CMI. Both metrics were tested in a numerically simulated example of unidirectionally coupled Rössler systems, which helped to find the explanation of the contradictory results: PSI correctly estimates the lead/lag relationship which, however, is not generally equivalent to causality in the sense of directionality of coupling in nonlinear systems, correctly inferred by using CMI with surrogate data testing.


Introduction
The phenomenon of cross-frequency coupling (CFC) in brain electrophysiological signals has recently received much attention in human and animal studies related to various conditions of health and disease.Beyond the synchronization phenomena in particular temporal scales, i.e., within particular frequency bands, the cross-frequency coupling enriches the cooperative behavior of neuronal networks and apparently plays an important functional role in neuronal computation, communication, and learning (Canolty and Knight, 2010).Many computational methods for detecting CFC, understood as a statistical relationship between a combination of amplitude, phase, and frequency of two distinct frequency bands, have been proposed (Jirsa and Müller, 2013;Yakubov et al., 2022;Yeh et al., 2023) and critically evaluated (Aru et al., 2015;Yeh et al., 2023).Although CFC could refer to any interaction or coherence between frequencies (Osipova et al., 2008), the phaseamplitude coupling (PAC) -in which the phase of slower oscillations is coupled to the amplitude of neuronal activity in higher frequency bands * Correspondence to: Institute of Computer Science of the Czech Academy of Sciences, Pod Vodárenskou věží 2, 18200 Prague 8, Czech Republic.E-mail address: mp@cs.cas.cz(M.Paluš).
In order to understand coherence or synchronization phenomena, it is necessary to uncover how the synchronizing oscillatory processes are coupled, i.e., either mutually or directionally, when one process is influencing another in a causal relation (Paluš et al., 2001).However, studies that specifically measure cross-frequency directionality are just emerging (Besserve et al., 2010;Martínez-Cancino et al., 2020).In this study we compare two techniques proposed to infer directionality of CFC/PAC.Conditional mutual information (CMI) (Paluš et al., 2001;Paluš, 2014) also known as transfer entropy (TE) (Schreiber, 2000) has been proposed as information-theoretic generalization of Granger causality for nonlinear systems (Hlaváčková-Schindler et al., 2007).Nolte et al. (2008) proposed phase slope index (PSI) as a measure of directionality for oscillatory signals of close frequencies, applicable particularly in electroencephalogram (EEG) analysis.Jiang et al. (2015) . Arinyo-i-Prats et al. directionality by applying PSI to the phase of low-frequency oscillations and the phase extracted from the amplitude of high frequency oscillations, calling the whole method Cross-Frequency Directionality (CFD).Furthermore, in this study we use Mutual Information (MI, see Paluš (2007) and references therein) to find lead-lag relationship between two oscillatory signals (Chakraborty and van Leeuwen, 2022).
In this technical note we apply and compare these methods using the same simulated nonlinear signals and EEG recordings from rats in the case of phase-amplitude coupling and directionality.CFD/PSI gives opposite results when comparing to CMI.Our explanation of the observed contradiction is that PSI correctly estimates the lead/lag relationship which, however, is not generally equivalent to causality in the sense of directionality of coupling in nonlinear systems, correctly inferred by using CMI with surrogate data testing.In the following, we firstly describe briefly both CMI and CFD.Secondly, we present simulated data and EEG recordings.Lastly, we apply both methods to the nonlinear simulations and the EEG data and discuss differences in the results obtained.

Conditional mutual information
Consider random variables ,  , with probability distribution functions PDF (), () and the joint PDF (, ).Mutual information (MI) is a measure of general statistical dependence between  and  : . (1) Using the definition of the Shannon entropy () = − ∑  () log (), the mutual information can be rewritten as (2) Based on the conditional entropy, we can define the conditional mutual information (;  |) of variables ,  given the variable  as Then, considering time series {()}, {()} as realizations of stochastic processes {()}, { ()}, or, alternatively (see Paluš and Vejmelka (2007) for details) as projections of trajectories of multidimensional dynamical systems  and  , it is possible to define a measure of directionality as time-delayed conditional mutual information  , : where  is the reconstruction time lag and  is the forward time lag, or prediction horizon in the Granger causality concept (Paluš and Vejmelka, 2007).See Supplementary material for further details.

Cross-frequency directionality/phase slope index
Nolte et al. (2008) introduced PSI as a robust method to infer coupling directionality between two signals,  and , of close frequencies, considering that a receiver follows the driver by a fixed time lag.PSI estimates the slope of the phase difference as a function of frequency in a given frequency band.The sign of the slope reflects the direction of interaction.
First, the data -the two signals,  and , are segmented into  windows and the complex coherence is computed as where  is the cross-spectral matrix across windows and  is the frequency.Then, PSI is defined as: where  is the frequency resolution, determined by the window length, ℑ denotes taking the imaginary part. is the set of frequencies over which the slope is summed.
For the inference of CFD, Jiang et al. (2015) extended the application of PSI computing the coupling between a low-frequency signal () with the frequency  and the amplitude (envelope) of faster oscillations () with the frequency .We refer to   and   as the temporal signals for the segment or window , and   is the envelope of   filtered at frequency .In this case, the complex coherence is rewritten as: where and are the Fourier transforms of   and   after applying a Hanning window and     defines the frequency resolution.
Then, for the definition of CFD, the PSI can be rewritten as: where  denotes the bandwidth for which the phase slope is calculated.

Surrogate data testing procedure
Estimates of dependence/directionality measures from finite number of samples are always nonzero and therefore it is necessary to relate the values computed from studied data to ranges obtained from uncoupled processes sharing statistical properties of analyzed data.This is the base of the surrogate data tests in which we manipulate the original data in a randomization procedure preserving original frequency spectra or variance on all relevant time scales.Here we use the FFT-based surrogate data (Paluš, 2007;Lancaster et al., 2018) for PSI in the simulated data example.The circular time-shifted surrogates, shown effective for causality inference (Manshour et al., 2021), are used for evaluating CMI significance in the real data application.The results of the surrogate data tests can be represented as the -score, e.g., for CMI, marked as , it is where   is the CMI value estimated from the studied data,   is the mean for 100 realizations of the surrogate data and  2  is their variance.Typically, the results are considered statistically significant for  > 2. See, e.g., the results for CMI in Fig. 1.Alternatively, the surrogate ranges can be illustrated as gray curve and whiskers (e.g., in Fig. 2), showing the mean ±2 (standard deviation) of a set of 100 surrogate data realizations.

Coupled Rössler systems
The simulated data generated by the unidirectionally coupled nonlinear Rössler systems is used as a benchmark for assessing the coupling directionality methods (details in Paluš and Vejmelka (2007)).
The autonomous system is defined as: and the response system: The term  (  1 −  1 ) in Eq. ( 12) is the diffusive coupling through which the system  influences the system  , written as  ← ← →  .For  = 0.15,  = 0.2,  = 10.0, and frequencies  1 = 1.015 and  2 = 0.985, we call this system  1 .Then we swap the values  1 and  2 calling the second system  2 .In both systems  is autonomous, driving →  by its positive values and the opposite direction by its negative values; the last row LAG is the lag maximizing cross-mutual information of the phases of both time series.The top panel presents these results for the  1 system, the bottom panel for  2 .The gray curve and whiskers illustrate the surrogate means and mean ±2 (standard deviations) for the set of surrogate data realizations.
through the diffusive coupling term ( 1 −  1 ), characterized by the coupling strength .However, for  1 ,  is faster than  , and in  2  is slower.Time series are generated as solutions of the systems (11), ( 12) for different values of  ∈ [0 − 0.25] for  1 and  ∈ [0 − 0.1] for  2 .See Supplementary material for further details.

Code availability
FORTRAN codes for CMI estimation are available at http://www.cs.cas.cz/mp/projects/sw/.Codes for CFD estimation from EEG data are available at https://github.com/VictorLopezMadrona/Matlab_Neuroscience, PSI code from https://doc.ml.tu-berlin.de/causality/.FORTRAN codes for generating and processing the simulated data are available upon request from M. Paluš.

EEG data
Electrophysiological recordings from the hippocampus of one rat were acquired with a 32 channels silicon probe (Neuronexus Technologies, Michigan, USA).The sampling frequency was 5 kHz and the data was filtered between 0.5 Hz and 300 Hz, and at 50 Hz and 100 Hz to remove the net noise.Then, independent component analysis (Makarov et al., 2010) was applied on the whole dataset, obtaining three main components related to the CA3 Schaffer collateral (Sch-IC) and the projections from the entorhinal cortex to the CA1 stratum lacunosum-moleculare (lm-IC) and to the dentate gyrus through the perforant pathway (PP-IC).For details see López-Madrona et al. (2020).In this work we used the first 65 536 samples of the animal s10.

Ethics statement
An adult male Long-Evans rat, with a weight of 250-300 g, was recorded while the animal freely explored a known 50 × 50 cm open field.All animal experiments were approved by the Animal Care and Use Committee of the Instituto de Neurociencias de Alicante, Alicante, Spain, and comply with the Spanish (law 32/2007) and European regulations (EU directive 86/609, EU decree 2001-486, and EU recommendation 2007/526/EC).

Data availability
The EEG data are available at DIGITAL.CSIC (Canals-Gamoneda and Álvarez-Salvado, 2020).The simulated data were obtained by numerical integration of the coupled Rössler systems using the Bulirsch-Stoer method (Press et al., 1992) which uses adaptive integration steps, however, the final sampling time 0.314 was prescribed.Details or data are available upon request from M. Paluš.

EEG cross-frequency directionality from CMI and CFD
In Fig. 1 we compare the CMI and CFD results, returning considerable differences in the phase-amplitude directionality for the three EEG signals from the rat hippocampus, Sch-IC, lm-IC and PP-IC, in the frequency bands from 5 to 15 Hz for the LF phase and from 40 to 250 Hz for the HF amplitude.In the color-scale we show the -score for each measure (CMI in Fig. 1A-F, violet color; CFD in Fig. 1G-L, green color) for || > 2. In the case of the CFD, we separated positive values, detecting the coupling from the low frequency phase to the high frequency amplitude (PAC), from negative values, signifying the opposite coupling from the high frequency amplitude to the low frequency phase (APC).Note that CFD only decides for one specific direction of coupling at each pair of frequency bands, as it cannot simultaneously return both positive and negative values.Thus, unlike CMI which tests the two directions separately, there are no overlaps of directions detected by CFD.
For the comparison of the methods, we only used a sequence of 65 536 samples from one animal (s10).Still, for this reduced dataset, the comoludogram produced using CFD shows the APC first described by López-Madrona et al. (2020).
In the case of CMI (A-F, violet color) we can see that the dominant direction of coupling is that of phase ← ← → amplitude, i.e., from low frequency to high frequency oscillations, while in the case of CFD (G-L, green color) the dominant coupling is directed from high-frequency to low-frequency oscillations.Following Osipova et al. (2008), in the upper panels of Fig. 2 we computed the coherence (COH) between the series  1 and  1 in a spectral band including both (close) system frequencies  1 and  2 (see Fig. S1 in Supplementary material).The coherence increases until the synchronization threshold, ( ≈ 0.12 for  1 ,  ≈ 0.055 for  2 ), when the systems are fully coherent.The gray line with whiskers represent means ± two standard deviations of the surrogate data distribution, representing the null hypothesis of independence, i.e., uncoupled systems.

Analysis of coupled Rössler systems
In the middle panels of Fig. 2 we show the CMI averaged for a range of lags (cf Fig. 4A, B) for both directions, first for  ← ← →  , followed by the opposite direction  ← ← → .CMI points to the system forcing  ← ← →  until the level of full coherence, where also the opposite direction  ← ← →  becomes significant.Recall that, in nonlinear systems, coupling direction can be inferred when the systems are coupled but not yet fully synchronized (Paluš and Vejmelka, 2007).
In the lower-middle panels of Fig. 2 we show the results of the PSI applied to the signals of both oscillators.There are many values of  (under the synchronization threshold) in which the PSI values for the  1 system are significantly negative indicating the directionality  ← ← → , i.e., opposite to the one simulated ( ← ← →  ).Only for the coherent systems (large ), the significantly positive PSI shows consistently the correct directionality  ← ← →  .In the case of the coupled  2 system, the PSI becomes significantly positive from  ≈ 0.03 till  ≈ 0.045.This correct directionality detection is followed by PSI fluctuations around the transition to the synchronized state where, for large , the significantly positive PSI again shows the correct directionality  ← ← →  .Finally, in the bottom graphs (LAG) we show the lag which maximizes the lagged Mutual Information  (   ();   ( + ) ) of the phases of the coupled Rössler oscillators.Thus, a positive lag  means that  is delayed and follows , or, in other words,  is leading  .A negative lag  means that  is leading .Comparing the sign of the lag maximizing the MI and the significant values of PSI, we can see an agreement in most cases.Recall that the directionality of coupling is in the direction  ← ← →  for both  1 and  2 systems.It seems that the PSI correctly infers the lead/lag relationship which, however, is not generally equivalent to the direction of coupling.

Lead/lag relationships versus direction of coupling/causality
In order to better understand the above findings, we present detailed analysis of CMI and lagged cross-MI obtained from the instantaneous phases of the two Rösler systems  1 , for  1 = 0.0625, and  2 , for  2 = 0.045.A segment of analyzed signals  1 and  1 and related instantaneous phases   and   are illustrated in Fig. 3.It can hardly be determined which signal is leading which just by looking at the waveforms or phases, although these values of  were chosen so that both CMI and PSI give statistically significant results (Fig. 2).
We also present CMI and lagged cross-MI for EEG data to see the relations between the phase of 6.5 Hz (8 Hz) oscillations and the amplitude of 167 Hz (100 Hz) oscillations.These frequency pairs were chosen to fit into the areas of significant values of either PSI or CMI in Fig. 1.

CMI for Rössler systems
In Fig. 4A, B, CMI as a function of time lag  is presented for the two Rössler systems  1 (A), and  2 (B).In order to distinguish the two coupling directions, the red dashed curve represents CMI for the direction  ← ← →  ; and the solid blue line represents the opposite direction  ← ← → .The red dashed curve is always distinctively positive, while the blue curve stays near the zero value, thus the true directionality  ← ← →  was inferred for both cases,  1 and  2 , in accordance with the coupling defined in Eqs. ( 11) and (12).

MI for Rössler systems
We computed the lagged cross-MI of the phases   () and   () for both cases of the coupled Rössler oscillators.In Fig. 4C, D ) .The time lag for which MI is maximized (marked by the vertical lines) represents the best alignment of the two phase series and thus uncover the lead/lag relation of the two systems.For  2 (Fig. 4D) the maximum is found in  (   ();   (+) ) for  = 18, thus we can see that  lags , or  leads  by 18 samples which agrees with the direction of coupling, or causality relation given by Eqs. ( 11) and (12).However, for  1 (Fig. 4C) we can see that the maximum is find in  (   ( + );   () ) for  = 7, which is equivalent to  (   ();   ( + ) ) for  = −7, therefore  lags  or  leads  by 7 samples, which contradicts the linear understanding of causality.We can see that in nonlinear oscillatory systems the lead/lag relationship is not equivalent to causality or the direction of coupling.See also Supplementary material for another example of nonlinear oscillatory systems.

CMI for EEG
The lm-IC channel has been selected for the example of the CMI and MI analysis of EEG signals.For the relation of the 8 Hz phase and the 100 Hz amplitude (Fig. 5B) CMI returns high values for the directionality LF phase ← → HF amplitude (dashed red curve), for the direction HF amplitude ← ← → LF phase (solid blue curve) the CMI stays near the zero value.This result corresponds to the violet spot of significant coupling directionality in Fig. 1B.On the other hand, in the case of the 6.5 Hz phase and 167 Hz amplitude (Fig. 5A) CMI has similarly low values in both directions, i.e., no statistically significant causality was detected by CMI for this frequency pair, while a significant area in PSI comodulogram can be seen (Fig. 1K).

MI for EEG
The lagged cross-MI was again used to establish the lead/lag relationship, but in the EEG case between the phase of the low frequency oscillations and the amplitude of the high frequency oscillations.In both cases the maxima are located in the solid green curve (Fig. 5 C, D) depicting the mutual information  (  ();  ℎ( + ) ) as the function of the time lag .Thus, the HF amplitude is leading LF phase by 102 and 26 samples for the relation 6.5 Hz phase-167 Hz amplitude and 8 Hz phase-100 Hz amplitude, respectively.

Discussion
In the present study we have shown disagreements of two coupling directionality measures: cross-frequency directionality (CFD) based on the phase slope index (PSI) and the conditional mutual information (CMI, also known as transfer entropy).First, in an animal study using EEG from rats (Fig. 1), CFD suggested that in the studied EEG the highfrequency amplitude drove the low-frequency phase.This finding was just a reproduction of the results of López-Madrona et al. (2020).On the contrary, the CMI results in Fig. 1 were dominated by the opposite causality or coupling directionality, meaning that the low-frequency phase drives the high-frequency amplitude.In order to understand this discrepancy, we applied CMI and PSI, the essence of CFD, on phases of unidirectionally coupled Rössler oscillators.The Rössler systems had close, but different natural frequencies.We studied both combinations of coupling: faster-to-slower and slower-to-faster system, and PSI inferred the correct direction of coupling and agreed with CMI only in one of these cases.Later we performed the same analysis with van der Pol oscillators (see Supplementary material) and found that correct and incorrect results of PSI were obtained apparently randomly, irrespectively of natural frequency ratios.
The conditional mutual information, however, in all cases correctly inferred the coupling direction given by the mathematical model used to generate the data.
Then we applied the lagged mutual information (MI) on the phases of both types of coupled Rössler systems.Identifying the lag which maximizes the lagged MI we were able to distinguish which system leads the other one.When PSI was statistically significant, its ''direction'' agreed with the lead/lag relationship given by the lagged MI.Thus we are able to claim that PSI correctly identifies the lead/lag relationship which is supposed to agree with the direction of causal influences in linear systems.In nonlinear system, however, the lead/lag relationship is not generally equivalent to the causality/coupling direction.Then we studied the lagged MI, as well as CMI for selected combinations of low-frequency phase and high-frequency amplitude of the rat EEG.We could see that for some frequency combinations CMI clearly inferred the coupling direction LF phase ← ← → HF amplitude, while MI showed that the high-frequency amplitude led the low-frequency phase.This analysis suggests that the cross-frequency coupling in the studied EEG is a nonlinear phenomenon for which the lead/lag relation, uncovered by both PSI and lagged MI, is not equivalent to the causal relation or coupling direction as it is understood in coupled nonlinear dynamical systems.Therefore, for constructions of mathematical models in computational neuroscience, results inferred by CMI are relevant, not the lead/lag relations, uncovered by CFD/PSI.The conditional mutual information (transfer entropy) is a promising directionality measure for analyzing the cross-frequency interactions, as already proposed by Besserve et al. (2010) or Martínez-Cancino et al. (2020) and demonstrated by Gupta and Paluš (2021) who applied CMI to data from simulated epileptic seizures and successfully uncovered all directional interactions between different time-scales included in the Epileptor model.We believe that CMI can contribute to better understanding of nonlinear cross-frequency interactions in the brain dynamics.

Fig. 1 .
Fig. 1.Comoludograms of the EEG data.Results for the three main components Sch-IC (A, D, G, J), lm-IC (B, E, H, K) and PP-IC (C, F, I, L).Areas of statistically significant directionality of cross-frequency phase-amplitude coupling obtained from CMI for the influence of low-frequency phase on high-frequency amplitude (LFP ← → HFA: A, B, C), and for HFA ← → LFP (D, E, F); and CFD for LFP ← → HFA (G, H, I), and for HFA ← → LFP (J, K, L); renormalized with the -score.Only values of || > 2 are shown.

Fig. 2 .
Fig. 2. Comparison of dependence and directionality measures for the coupled Rössler systems with varying coupling strength.Top row coherence COH, second row CMI for  ← →  direction, third row CMI for  ← →  direction, fourth row PSI, indicating the direction  ←→  by its positive values and the opposite direction by its negative values; the last row LAG is the lag maximizing cross-mutual information of the phases of both time series.The top panel presents these results for the  1 system, the bottom panel for  2 .The gray curve and whiskers illustrate the surrogate means and mean ±2 (standard deviations) for the set of surrogate data realizations.

Fig. 2
Fig. 2 presents different proxies of coupling and directionality for both Rössler systems  1 and  2 , as functions of coupling strengths .Examples of typical signals  1 and  1 obtained by integrating Rössler systems  1 and  2 , and related instantaneous phases   and   are presented in Fig. 3.FollowingOsipova et al. (2008), in the upper panels of Fig.2we computed the coherence (COH) between the series  1 and  1 in a spectral band including both (close) system frequencies  1 and  2 (see Fig.S1in Supplementary material).The coherence increases

Fig. 3 .
Fig.3.Rössler systems data and instantaneous phases.From top to bottom: A segment of time evolution of variables  1 and  1 and related phases   and   for the Rössler systems  1 for  = 0.0625; a segment of time evolution of variables  1 and  1 and related phases   and   for the Rössler systems  2 for  = 0.045.

Fig. 5 .
Fig. 5. Causality and lead/lag relationship in EEG.CMI [nats] (A, B) and MI [nats] (C, D) as functions of time lag , for EEG signal in the lm-IC channel frequency pair 6.5 Hz-167 Hz (A, C) and EEG signal in the lm-IC channel frequency pair 8 Hz-100 Hz (B, D).For CMI the red dashed curve the direction LFP → HFA, the solid blue curve HFA → LFP.For MI solid green curve indicates HFA leading LFP, dashed purple indicates LFP leading HFA.Vertical lines indicate lags maximizing MI.