Note

Go to the end to download the full example code

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

## 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)
)
```

## References#

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