Compute phase-amplitude coupling#

This example demonstrates how phase-amplitude coupling (PAC) can be computed with PyBispectra.

import numpy as np

from pybispectra import compute_fft, get_example_data_paths, PAC

Background#

PAC quantifies the relationship between the phases of a lower frequency \(f_1\) and the amplitude of a higher frequency \(f_2\) within a single signal, \(\textbf{x}\), or across different signals, \(\textbf{x}\) and \(\textbf{y}\).

The method available in PyBispectra is based on the bispectrum, \(\textbf{B}\). The bispectrum has the general form

\(\textbf{B}_{kmn}(f_1,f_2)=<\textbf{k}(f_1)\textbf{m}(f_2)\textbf{n}^* (f_2+f_1)>\) ,

where \(kmn\) is a combination of signals with Fourier coefficients \(\textbf{k}\), \(\textbf{m}\), and \(\textbf{n}\), respectively; and \(<>\) represents the average value over epochs. The computation of PAC follows from this [1]

\(\textbf{B}_{xyy}(f_1,f_2)=<\textbf{x}(f_1)\textbf{y}(f_2)\textbf{y}^* (f_2+f_1)>\) ,

\(\textrm{PAC}(\textbf{x}_{f_1},\textbf{y}_{f_2})=|\textbf{B}_{xyy}(f_1, f_2)|\) .

The bispectrum can be normalised to the bicoherence, \(\boldsymbol{\mathcal{B}}\), using the threenorm, \(\textbf{N}\), [2]

\(\textbf{N}_{xyy}(f_1,f_2)=(<|\textbf{x}(f_1)|^3><|\textbf{y}(f_2)|^3> <|\textbf{y}(f_2+f_1)|^3>)^{\frac{1}{3}}\) ,

\(\boldsymbol{\mathcal{B}}_{xyy}(f_1,f_2)=\Large\frac{\textbf{B}_{xyy} (f_1,f_2)}{\textbf{N}_{xyy}(f_1,f_2)}\) ,

\(\textrm{PAC}_{\textrm{norm}}(\textbf{x}_{f_1},\textbf{y}_{f_2})=| \boldsymbol{\mathcal{B}}_{xyy}(f_1,f_2)|\) .

where the resulting values lie in the range \([0, 1]\). Furthermore, PAC can be antisymmetrised by subtracting the results from those found using the transposed bispectrum, \(\textbf{B}_{xyx}\), [3]

\(\textrm{PAC}_{\textrm{antisym}}(\textbf{x}_{f_1},\textbf{y}_{f_2})=| \textbf{B}_{xyy}-\textbf{B}_{xyx}|\) .

In the context of analysing PAC between two signals, antisymmetrisation allows you to correct for spurious estimates of coupling arising from interactions within the signals themselves in instances of source mixing, providing a more robust connectivity metric [4]. The same principle applies for the antisymmetrisation of the bicoherence.

Loading data and computing Fourier coefficients#

We will start by loading some simulated data containing coupling between the 10 Hz phase of one signal and the 60 Hz amplitude of another. We will then compute the Fourier coefficients of the data, which will be used to compute PAC.

# load simulated data
data = np.load(get_example_data_paths("sim_data_pac_bivariate"))
sampling_freq = 200  # sampling frequency in Hz

# compute Fourier coeffs.
fft_coeffs, freqs = compute_fft(
    data=data,
    sampling_freq=sampling_freq,
    n_points=sampling_freq,
    verbose=False,
)

print(
    f"FFT coeffs.: [{fft_coeffs.shape[0]} epochs x {fft_coeffs.shape[1]} "
    f"channels x {fft_coeffs.shape[2]} frequencies]\nFreq. range: {freqs[0]} "
    f"- {freqs[-1]} Hz"
)
FFT coeffs.: [30 epochs x 2 channels x 101 frequencies]
Freq. range: 0.0 - 100.0 Hz

As you can see, we have Fourier coefficients for 2 channels across 30 epochs, with 101 frequencies ranging from 0 to 100 Hz with a frequency resolution of 1 Hz. We will use these coefficients to compute PAC.

Computing PAC#

To compute PAC, we start by initialising the PAC class object with the Fourier coefficients and the frequency information. To compute PAC, we call the compute() method. By default, PAC is computed between all channel and frequency combinations, however we can also specify particular combinations of interest.

Here, we specify the indices to compute PAC on. indices is expected to be a tuple containing two lists for the indices of the seed and target channels, respectively. The indices specified below mean that PAC will only be computed across frequencies between the channels (i.e. 0 -> 1). By leaving the frequency arguments f1s and f2s blank, we will look at all possible frequency combinations.

pac = PAC(
    data=fft_coeffs, freqs=freqs, sampling_freq=sampling_freq, verbose=False
)  # initialise object
pac.compute(indices=((0,), (1,)))  # compute PAC

pac_results = pac.results.get_results()  # return results as array

print(
    f"PAC results: [{pac_results.shape[0]} connection x "
    f"{pac_results.shape[1]} f1s x {pac_results.shape[2]} f2s]"
)
PAC results: [1 connection x 101 f1s x 101 f2s]

We can see that PAC has been computed for one connection (0 -> 1), and all possible frequency combinations, averaged across our 30 epochs. Whilst there are > 10,000 such frequency combinations in our [101 x 101] matrix, PAC for those entries where \(f_1\) would be higher than \(f_2\), as well as where \(f_2 + f_1\) exceeds the frequency bounds of our data, cannot be computed. In such cases, the values are numpy.nan.

Plotting PAC#

Let us now inspect the results. Here, we specify a subset of frequencies to inspect around the simulated interaction. If we wished, we could also plot all frequencies. Note that the Figure and Axes objects can also be returned for any desired manual adjustments of the plots. In this simulated data example, we can see that the bispectrum indeed identifies the occurrence of 10-60 Hz PAC between our seed and target channel.

fig, axes = pac.results.plot(f1s=(5, 15), f2s=(55, 65))
PAC | Bispectrum, Seed: 0 | Target: 1

Antisymmetrisation for across-signal PAC#

In this simulated data example, interactions are only present between the signals, and not within the signals themselves. This, however, is not always the case, and estimates of across-site PAC can be corrupted by coupling interactions within each signal in the presence of source mixing. To combat this, we can employ antisymmetrisation [3]. The example below shows some such simulated data consisting of two independent sources, with 10-60 Hz PAC within each source (top two plots), as well as a mixing of the underlying sources to produce 10-60 Hz PAC between the two signals (bottom left plot). When appyling antisymmetrisation, however, we see that the spurious across-signal PAC arising from the source mixing is suppressed (bottom right plot). Antisymmetrisation is therefore a useful technique to differentiate genuine across-site coupling from spurious coupling arising from the within-site interactions of source-mixed signals.

# load real data
data = np.load(get_example_data_paths("sim_data_pac_univariate"))
sampling_freq = 200

# compute Fourier coeffs.
fft_coeffs, freqs = compute_fft(
    data=data,
    sampling_freq=sampling_freq,
    n_points=sampling_freq,
    verbose=False,
)

# compute PAC
pac = PAC(
    data=fft_coeffs, freqs=freqs, sampling_freq=sampling_freq, verbose=False
)
pac.compute(
    indices=((0, 1, 0), (0, 1, 1)),
    f1s=(5, 15),
    f2s=(55, 65),
    antisym=(False, True),
)
pac_standard, pac_antisym = pac.results

vmin = np.min(
    (
        np.nanmin(pac_standard.get_results()),
        np.nanmin(pac_antisym.get_results()),
    )
)
vmax = np.max(
    (
        np.nanmax(pac_standard.get_results()),
        np.nanmax(pac_antisym.get_results()),
    )
)

# plot unsymmetrised PAC within & between signals
fig_standard, axes_standard = pac_standard.plot(
    f1s=(5, 15), f2s=(55, 65), cbar_range=(vmin, vmax)
)

# plot antisymmetrised PAC between signals
fig_antisym, axes_antisym = pac_antisym.plot(
    nodes=(2,), f1s=(5, 15), f2s=(55, 65), cbar_range=(vmin, vmax)
)
  • PAC | Bispectrum, Seed: 0 | Target: 0
  • PAC | Bispectrum, Seed: 1 | Target: 1
  • PAC | Bispectrum, Seed: 0 | Target: 1
  • PAC (antisymmetrised) | Bispectrum, Seed: 0 | Target: 1

References#

Total running time of the script: (0 minutes 1.844 seconds)

Gallery generated by Sphinx-Gallery