<figure><img src="image_1.png" /></figure>

Technical Note

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

Andreu Arinyo-i-Prats<sup>a,c,d</sup>, Víctor J. López-Madrona<sup>b</sup>, Milan Paluš<sup>a,*</sup>

<figure><img src="image_2.png" /></figure>

ARTICLE INFO

Dataset link: https://doi.org/10.20350/digitalC
SIC/12537

Keywords:

Coupling directionality

Cross-frequency coupling

Conditional mutual information

Phase slope index

EEG

Nonlinear systems

A B S T R A C T

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.

# 1. 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 phase-amplitude coupling (PAC) –in which the phase of slower oscillations is coupled to the amplitude of neuronal activity in higher frequency bands

— has attracted increasing interest (Canolty et al., 2006; De Hemptinne et al., 2013).

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) incorporated the PSI into CFC analysis in order to infer PAC

<footnote><sup>a</sup> Department of Complex Systems, Institute of Computer Science of the Czech Academy of Sciences, Pod Vodárenskou věží 2, Prague, 18200, Czech Republic</footnote>

<footnote><sup>b</sup> Institut de Neurosciences des Systèmes, Aix Marseille Univ, INSERM, INS, Inst Neurosci Syst, 27 Bd Jean Moulin, Marseille, 13005, France</footnote>

<footnote><sup>c</sup> Department of Archaeology and Heritage Studies, School of Culture and Society, Aarhus University, Jens Chr. Skous Vej 7, Building 1467, Aarhus, 8000, Denmark</footnote>

<footnote><sup>d</sup> Department of Archaeology, Simon Fraser University, Education Building 9635, 8888 University Drive, Burnaby, V5A 1S6, B.C., Canada</footnote>

<footnote>* 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š).</footnote>
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.

## 2. Methods

### 2.1. Conditional mutual information

Consider random variables X, Y, with probability distribution functions PDF p(x), p(y) and the joint PDF p(x, y). Mutual information (MI) is a measure of general statistical dependence between X and Y:

$$I(X;Y) = \sum_{x} \sum_{y} p(x,y) \log \frac{p(x,y)}{p(x)p(y)}. \qquad (1)$$

Using the definition of the Shannon entropy $H(X) = -\sum_{x} p(x) \log p(x)$, the mutual information can be rewritten as

$$I(X,Y) = H(X) + H(Y) - H(X,Y). \qquad (2)$$

Based on the conditional entropy, we can define the conditional mutual information $I(X; Y|Z)$ of variables X, Y given the variable Z as

$$I(X;Y|Z) = H(X|Z) + H(Y|Z) - H(X,Y|Z). \qquad (3)$$

Then, considering time series {$x(t)$}, {$y(t)$} as realizations of stochastic processes {$X(t)$}, {$Y(t)$}, or, alternatively (see Paluš and Vejmelka (2007) for details) as projections of trajectories of multidimensional dynamical systems X and Y, it is possible to define a measure of directionality as time-delayed conditional mutual information $T^{\tau, \eta}$:

$$T_{\tau \to \eta}^{\tau, \eta}(X(t); Y(t+\tau)|Y(t)) = \qquad (4)$$

$$I(x(t); y(t+\tau)|y(t), y(t-\eta), \dots, y(t-(n-1)\eta)),$$

where $\eta$ is the reconstruction time lag and $\tau$ is the forward time lag, or prediction horizon in the Granger causality concept (Paluš and Vejmelka, 2007). See Supplementary material for further details.

### 2.2. Cross-frequency directionality/phase slope index

Nolte et al. (2008) introduced PSI as a robust method to infer coupling directionality between two signals, i and j, 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, i and j, are segmented into K windows and the complex coherence is computed as $C_{ij}(f) = S_{ij}(f)/\sqrt{|S_{ii}(f)||S_{jj}(f)|}$, where S is the cross-spectral matrix across windows and f is the frequency. Then, PSI is defined as:

$$\tilde{\Psi}_{ij} = \mathfrak{F} \left( \sum_{f \in F} C_{ij}^{*}(f) C_{ij}(f + \Delta f) \right), \qquad (5)$$

where $\Delta f$ is the frequency resolution, determined by the window length, $\mathfrak{F}$ denotes taking the imaginary part. F 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 (x) with the frequency f and the amplitude (envelope) of faster oscillations (y) with the frequency v. We refer to $x^k$ and $y^k$ as the temporal signals for the segment or window k, and $y^{vk}$ is the envelope of $y^k$ filtered at frequency v. In this case, the complex coherence is rewritten as:

$$C(v, f) = \frac{\sum_{k=1}^{K} \mathbf{X}^k (\mathbf{Y}^{vk})^*}{\sqrt{\sum_{k=1}^{K} |\mathbf{X}^k|^2 \sum_{k=1}^{K} |\mathbf{Y}^{vk}|^2}}, \qquad (6)$$

where

$$\mathbf{X}^k = \text{FFT}(\mathbf{h}^T \mathbf{x}^k, n_{\text{FFT}}) \qquad (7)$$

and

$$\mathbf{Y}^{vk} = \text{FFT}(\mathbf{h}^T \mathbf{y}^{vk}, n_{\text{FFT}}) \qquad (8)$$

are the Fourier transforms of $x^k$ and $y^{vk}$ after applying a Hannin window and $n_{FFT}$ defines the frequency resolution.

Then, for the definition of CFD, the PSI can be rewritten as:

$$\Psi(v, f) = \mathfrak{F} \left( \sum_{f-\frac{\beta}{2}}^{f+\frac{\beta}{2}} C^{*}(v, f) C(v, (f + \Delta f)) \right), \qquad (9)$$

where $\beta$ denotes the bandwidth for which the phase slope is calculated.

### 2.3. 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 Z-score, e.g., for CMI, marked as I, it is

$$Z = \frac{I_d - \bar{I}_s}{\sigma_s}, \qquad (10)$$

where $I_d$ is the CMI value estimated from the studied data, $\bar{I}_s$ is the mean for 100 realizations of the surrogate data and $\sigma_s^2$ is their variance. Typically, the results are considered statistically significant for $Z > 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 $\pm 2\sigma$ (standard deviation) of a set of 100 surrogate data realizations.

### 2.4. Coupled Rössler systems

The simulated data generated by the unidirectionally coupled non-linear 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:

$$\begin{aligned} \dot{x}_1 &= -\omega_1 x_2 - x_3 \\ \dot{x}_2 &= \omega_1 x_1 + ax_2 \\ \dot{x}_3 &= b + x_3 (x_1 - c), \end{aligned} \qquad (11)$$
<figure><img src="image_3.png" /><figcaption>Fig. 1. Comoludograms of the EEG data. Results for the three main components Sch-IC (A, D, G, J), Im-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 Z-score. Only values of |Z| > 2 are shown.</figcaption></figure>

and the response system:

$$ \dot{y}_1 = -\omega_2 y_2 - y_3 + \epsilon (x_1 - y_1) $$

$$ \dot{y}_2 = \omega_2 y_1 + ay_2 \qquad (12) $$

$$ \dot{y}_3 = b + y_3 (y_1 - c). $$

The term $\epsilon (x_1 - y_1)$ in Eq. (12) is the diffusive coupling through which the system X influences the system Y, written as $X \to Y$.

For $a = 0.15$, $b = 0.2$, $c = 10.0$, and frequencies $\omega_1 = 1.015$ and $\omega_2 = 0.985$, we call this system $R_1$. Then we swap the values $\omega_1$ and $\omega_2$ calling the second system $R_2$. In both systems X is autonomous, driving
<figure><img src="image_4.png" /><figcaption>**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 $X \to Y$ direction, third row CMI for $Y \to X$ direction, fourth row PSI, indicating the direction $X \to Y$ 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 $R_1$ system, the bottom panel for $R_2$. The gray curve and whiskers illustrate the surrogate means and mean $\pm 2\sigma$ (standard deviations) for the set of surrogate data realizations.</figcaption></figure>

$Y$ through the diffusive coupling term $\epsilon(x_1 - y_1)$, characterized by the coupling strength $\epsilon$. However, for $R_1$, $X$ is faster than $Y$, and in $R_2$ $X$ is slower. Time series are generated as solutions of the systems (11), (12) for different values of $\epsilon \in [0-0.25]$ for $R_1$ and $\epsilon \in [0-0.1]$ for $R_2$. See Supplementary material for further details.

## 2.5. 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š.

## 2.6. 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 (Im-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 65536 samples of the animal s10.

## 2.7. 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).

## 2.8. 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š.

# 3. Results

## 3.1. 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**, **Im-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 Z-score for each measure (CMI in Fig. 1A–F, violet color; CFD in Fig. 1G–L, green color) for $|Z| > 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 65536 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.

## 3.2. Analysis of coupled Rössler systems

Fig. 2 presents different proxies of coupling and directionality for both Rössler systems $R_1$ and $R_2$, as functions of coupling strengths $\epsilon$. Examples of typical signals $x_1$ and $y_1$ obtained by integrating Rössler systems $R_1$ and $R_2$, and related instantaneous phases $\phi_X$ and $\phi_Y$ are presented in Fig. 3.

Following Osipova et al. (2008), in the upper panels of Fig. 2 we computed the coherence (COH) between the series $x_1$ and $y_1$ in a spectral band including both (close) system frequencies $\omega_1$ and $\omega_2$ (see Fig. S1 in Supplementary material). The coherence increases
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Rössler systems data and instantaneous phases. From top to bottom: A segment of time evolution of variables $x_1$ and $y_1$ and related phases $\phi_X$ and $\phi_Y$ for the Rössler systems $R_1$ for $\epsilon = 0.0625$; a segment of time evolution of variables $x_1$ and $y_1$ and related phases $\phi_X$ and $\phi_Y$ for the Rössler systems $R_2$ for $\epsilon = 0.045$.</figcaption></figure>

until the synchronization threshold, ($\epsilon \approx 0.12$ for $R_1$, $\epsilon \approx 0.055$ for $R_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.

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 $X \to Y$, followed by the opposite direction $Y \to X$. CMI points to the system forcing $X \to Y$ until the level of full coherence, where also the opposite direction $Y \to X$ 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 $\epsilon$ (under the synchronization threshold) in which the PSI values for the $R_1$ system are significantly negative indicating the directionality $Y \to X$, i.e., opposite to the one simulated ($X \to Y$). Only for the coherent systems (large $\epsilon$), the significantly positive PSI shows consistently the correct directionality $X \to Y$.

In the case of the coupled $R_2$ system, the PSI becomes significantly positive from $\epsilon \approx 0.03$ till $\epsilon \approx 0.045$. This correct directionality detection is followed by PSI fluctuations around the transition to the synchronized state where, for large $\epsilon$, the significantly positive PSI again shows the correct directionality $X \to Y$.

Finally, in the bottom graphs (LAG) we show the lag which maximizes the lagged Mutual Information $I(\phi_X(t); \phi_Y(t+\tau))$ of the phases of the coupled Rössler oscillators. Thus, a positive lag $\tau$ means that $Y$ is delayed and follows $X$, or, in other words, $X$ is leading $Y$. A negative lag $\tau$ means that $Y$ is leading $X$. 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 $X \to Y$ for both $R_1$ and $R_2$ systems. It seems that the PSI correctly infers the lead/lag relationship which, however, is not generally equivalent to the direction of coupling.

### 3.3. 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össler systems $R_1$, for $\epsilon_1 = 0.0625$, and $R_2$, for $\epsilon_2 = 0.045$. A segment of analyzed signals $x_1$ and $y_1$ and related instantaneous phases $\phi_X$ and $\phi_Y$ are illustrated in Fig. 3. It can hardly

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Causality and lead/lag relationship in Rössler systems. CMI [nats] (A, B) and MI [nats] (C, D) as functions of time lag $\tau$, for the Rössler systems $R_1$ (A, C) with $\omega_1 = 1.015$, $\omega_2 = 0.985$, $\epsilon_1 = 0.0625$; and for the Rössler systems $R_2$ (B, D) with $\omega_1 = 0.985$, $\omega_2 = 1.015$, $\epsilon_1 = 0.045$. For CMI the red dashed curve shows the direction $X \to Y$, the solid blue curve $Y \to X$. For MI solid green curve indicates $Y$ leading $X$, dashed purple curve indicates $X$ leading $Y$. Vertical lines indicate lags maximizing MI.</figcaption></figure>

be determined which signal is leading which just by looking at the waveforms or phases, although these values of $\epsilon$ 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.

#### 3.3.1. CMI for Rössler systems

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

#### 3.3.2. MI for Rössler systems

We computed the lagged cross-MI of the phases $\phi_X(t)$ and $\phi_Y(t)$ for both cases of the coupled Rössler oscillators. In Fig. 4C, D, the dashed purple curve represents $I(\phi_X(t); \phi_Y(t+\tau))$ for the positive lag $\tau$; the solid dark green curve represents the negative lag $\tau$ flipped into positive values or, equivalently, the lag for $I(\phi_X(t+\tau); \phi_Y(t))$. 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 $R_2$ (Fig. 4D) the maximum is found in $I(\phi_X(t); \phi_Y(t+\tau))$ for $\tau = 18$, thus we can see that $Y$ lags $X$, or $X$ leads $Y$ by 18 samples which agrees with the direction of coupling, or causality relation given by Eqs. (11) and (12). However, for $R_1$ (Fig. 4C) we can see that the maximum is find in $I(\phi_X(t+\tau); \phi_Y(t))$ for $\tau = 7$, which is equivalent to $I(\phi_X(t); \phi_Y(t+\tau))$ for $\tau = -7$, therefore $X$ lags $Y$ or $Y$ leads $X$ 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.
<figure><img src="image_7.png" /><figcaption>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 shows the direction LFP → HFA, the solid blue curve HFA → LFP. For MI solid green curve indicates HFA leading LFP, dashed purple curve indicates LFP leading HFA. Vertical lines indicate lags maximizing MI.</figcaption></figure>

3.3.3. 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), while 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).

3.3.4. MI for EEG

The lagged cross-MI was again used to establish the lead/lag rela-
tionship, 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 I (*HF amplitude*(*t*); *LF phase*(*t* + *τ*))
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.

4. 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 high-
frequency 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 in-
ferred 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 max-
imizes 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 com-
putational neuroscience, results inferred by CMI are relevant, not the
lead/lag relations, uncovered by CFD/PSI.

The conditional mutual information (transfer entropy) is a promis-
ing directionality measure for analyzing the cross-frequency interac-
tions, 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.

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š.

CRediT authorship contribution statement

**Andreu Arinyo-i-Prats:** Writing – original draft, Visualization, Software, Investigation, Data curation, Conceptualization. **Víctor J. López-Madrona:** Writing – review & editing, Software, Resources, Data curation. **Milan Paluš:** Software, Methodology, Funding acquisition, Formal analysis, Conceptualization, Supervision, Visualization, Writing – review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial
interests or relations that could have appeared to influence the work
reported in this paper.
## Data availability

The EEG data are available at DIGITAL.CSIC: Canals-Gamoneda, S., Álvarez-Salvado, E., 2020. Electrophysiology of rat hippocampus Novelty and TMaze.

https://doi.org/10.20350/digitalCSIC/12537.

The simulated data were obtained by numerical integration of the coupled Rössler systems using the Bulirsch-Stoer method which uses adaptive integration steps, however, the final sampling time 0.314 was prescribed. Details or data are available upon request from M. Paluš.

## Acknowledgments

This study is supported by the Czech Science Foundation (grant No. GF21-14727K); and by the Czech Academy of Sciences, Praemium Academiae awarded to M. Paluš.

## Appendix A. Supplementary information about causality and directionality of coupling, used methods and simulated data.

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.neuroimage.2024.120610.

## References

Aru, J., Aru, J., Prieseemann, V., Wibral, M., Lana, L., Pipa, G., Singer, W., Vicente, R., 2015. Untangling cross-frequency coupling in neuroscience. Curr. Opin. Neurobiol. 31, 51–61. http://dx.doi.org/10.1016/j.conb.2014.08.002.

Besserve, M., Schölkopf, B., Logothetis, N.K., Panzeri, S., 2010. Causal relationships between frequency bands of extracellular signals in visual cortex revealed by an information theoretic analysis. J. Comput. Neurosci. 29, 547–566.

Canals-Gamoneda, S., Álvarez-Salvado, E., 2020. Electrophysiology of rat hippocampus novelty and TMaze. http://dx.doi.org/10.20350/digitalCSIC/12537.

Canolty, R.T., Edwards, E., Dalal, S.S., Soltani, M., Nagarajan, S.S., Kirsch, H.E., Berger, M.S., Barbaro, N.M., Knight, R.T., 2006. High gamma power is phase-locked to theta oscillations in human neocortex. Science 313, 1626–1628.

Canolty, R.T., Knight, R.T., 2010. The functional role of cross-frequency coupling. Trends Cogn. Sci. 14, 506–515.

Chakraborty, N., van Leeuwen, P.J., 2022. Using mutual information to measure time lags from nonlinear processes in astronomy. Phys. Rev. Res. 4, 013036. http://dx.doi.org/10.1103/PhysRevResearch.4.013036.

De Hemptinne, C., Ryapolova-Webb, E.S., Air, E.L., Garcia, P.A., Miller, K.J., Ojemann, J.G., Ostrem, J.L., Galilianakis, N.B., Starr, P.A., 2013. Exaggerated phase-amplitude coupling in the primary motor cortex in parkinson disease. Proc. Natl. Acad. Sci. 110, 4780–4785.

Gupta, K., Paluš, M., 2021. Cross-scale causality and information transfer in simulated epileptic seizures. Entropy 23 (526), http://dx.doi.org/10.3390/e23050526.

Hlaváčková-Schindler, K., Paluš, M., Vejmelka, M., Bhattacharya, J., 2007. Causality detection based on information-theoretic approaches in time series analysis. Phys. Rep. 441, 1–46.

Jiang, H., Bahramisharif, A., van Gerven, M.A., Jensen, O., 2015. Measuring directionality between neuronal oscillations of different frequencies. Neuroimage 118, 359–367.

Jirsa, V., Müller, V., 2013. Cross-frequency coupling in real and virtual brain networks. Front. Comput. Neurosci. 7, http://dx.doi.org/10.3389/fncom.2013.00078.

Lancaster, G., Iatsenko, D., Pidde, A., Ticcinelli, V., Stefanovska, A., 2018. Surrogate data for hypothesis testing of physical systems. Phys. Rep. 748, 1–60. http://dx.doi.org/10.1016/j.physrep.2018.06.001.

López-Madrona, V.J., Pérez-Montoyo, E., Álvarez-Salvado, E., Moratal, D., Herreras, O., Pereda, E., Mirasso, C.R., Canals, S., 2020. Different theta frameworks coexist in the rat hippocampus and are coordinated during memory-guided and novelty tasks. Elife 9, e57313. http://dx.doi.org/10.7554/eLife.57313.

Makarov, V.A., Makarova, J., Herreras, O., 2010. Disentanglement of local field potential sources by independent component analysis. J. Comput. Neurosci. 29, 445–457.

Manshour, P., Balasis, G., Consolini, G., Papadimitriou, C., Paluš, M., 2021. Causality and information transfer between the solar wind and the magnetosphere-ionosphere system. Entropy 23 (390), http://dx.doi.org/10.3390/e23040390.

Martínez-Cancino, R., Delorme, A., Wagner, J., Kreutz-Delgado, K., Sotero, R.C., Makeig, S., 2020. What can local transfer entropy tell us about phase-amplitude coupling in electrophysiological signals? Entropy 22, 1262.

Nolte, G., Ziehe, A., Nikulin, V.V., Schlögl, A., Krämer, N., Brismar, T., Müller, K.R., 2008. Robustly estimating the flow direction of information in complex physical systems. Phys. Rev. Lett. 100, 234101.

Osipova, D., Hermes, D., Jensen, O., 2008. Gamma power is phase-locked to posterior alpha activity. PLoS One 3, e3990.

Paluš, M., 2007. From nonlinearity to causality: statistical testing and inference of physical mechanisms underlying complex dynamics. Contemp. Phys. 48, 307–348. http://dx.doi.org/10.1080/00107510801959206.

Paluš, M., 2014. Cross-scale interactions and information transfer. Entropy 16, 5263–5289. http://dx.doi.org/10.3390/e16105263.

Paluš, M., Komárek, V., Hrnčíř, Z., štěrbová, K., 2001. Synchronization as adjustment of information rates: Detection from bivariate time series. Phys. Rev. E 63, 046211.

Paluš, M., Vejmelka, M., 2007. Directionality of coupling from bivariate time series: How to avoid false causalities and missed connections. Phys. Rev. E 75, 056211.

Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1992. Numerical Recipes in FORTRAN 77: The Art of Scientific Computing. Cambridge University Press.

Schreiber, T., 2000. Measuring information transfer. Phys. Rev. Lett. 85 (461).

Yakubov, B., Das, S., Zomorrodi, R., Blumberger, D.M., Enticott, P.G., Kirkovski, M., Rajji, T.K., Desarkar, P., 2022. Cross-frequency coupling in psychiatric disorders: A systematic review. Neurosci. Biobehav. Rev. 138, 104690. http://dx.doi.org/10.1016/j.neubiorev.2022.104690.

Yeh, C.H., Zhang, C., Shi, W., Lo, M.T., Tinkhauser, G., Oswal, A., 2023. Cross-frequency coupling and intelligent neuromodulation. Cyborg Bionic Syst. 4 (0034), http://dx.doi.org/10.34133/cbsystems.0034.