Functional connectivity: Coherence

Initialize NeuroAnalyzer 
using NeuroAnalyzer
using DSP
using Plots
using Statistics
eeg = load("files/eeg.hdf")
e10 = epoch(eeg; ep_len = 10)

Coherence measures the consistency of the phase relationship between two signals across trials or time windows - quantifying how reliably two channels oscillate together at a given frequency.

Where cross-correlation operates in the time domain, coherence operates in the frequency domain: it asks not just whether two signals are related, but at which frequencies their relationship is stable.

Cross-Power Spectral Density

Coherence is computed from the Cross-Power Spectral Density (CPSD) - a frequency-domain measure of the relationship between two signals, describing how the power of one signal relates to another at each frequency.

In practice, for discrete signals \(x[n]\) and \(y[n]\), the CPSD is estimated by multiplying the Fast Fourier Transform (FFT) of one signal by the complex conjugate of the FFT of the second signal:

\[ P_{xy}(f) = X(f) \times Y^*(f) \]

where:
\(X(f)\) is the Fourier transform of \(x\),
\(Y^*(f)\) is the complex conjugate of the Fourier transform of \(y\).

The result \(P_{xy}\) is complex-valued at each frequency, encoding two distinct pieces of information:

  • Magnitude \(\vert P_{xy} \vert\) - the strength of the relationship between \(x\) and \(y\) at that frequency
  • Phase \(angle P_{xy}\) - the phase difference between the two signals at that frequency

Coherence is then obtained by normalizing the CPSD by the auto-spectral densities of both signals:

\[ C_{xy}(f) = \frac{\vert P_{xy}(f) \vert^2}{P_{xx}(f) \times P_{yy}(f)} \]

which scales the result to the range \([0, 1]\), removing the influence of each signal’s individual power and isolating the consistency of their phase relationship.

Computing CPSD:

cp, f = cpsd(e10,
             e10,
             ch1 = "Fp1",
             ch2 = "Fp2",
             ep1 = 10,
             ep2 = 10,
             method = :mt)

Three estimation methods are available via the method parameter:

  1. :fft - Fast Fourier transform; fast and straightforward, suitable for stationary signals
  2. :mt - multi-taper; reduces spectral leakage by averaging across multiple orthogonal tapers, producing a more stable estimate
  3. :stft - short-time Fourier transform; computes CPSD in sliding time windows, suitable for non-stationary signals where the frequency content changes over time

Plotting magnitudes:

plot_psd(f,
         abs.(cp[1, :, 1]),
         xlabel = "Frequency [Hz]",
         ylabel = "Magnitude")

Magnitude peaks in the CPSD indicate strong connectivity between the two signals at those specific frequencies - the larger the magnitude, the more consistently the two signals co-vary at that frequency across the observation window.

Plotting phases:

plot_phase(DSP.angle.(cp[1, :, 1]),
           f)

Phase information reveals the temporal relationship between the two signals at each frequency - a consistent non-zero phase indicates that one region systematically leads or lags the other, which can imply a directional influence between them.

CPSD can also be compared across conditions at the same electrode location - for example, rest vs. task, or healthy vs. patient groups - to identify frequency-specific changes in inter-regional connectivity associated with cognitive state or pathology.

Coherence

Definition

Coherence is defined as the normalized cross-spectral density between two signals \(x\) and \(y\):

\[ C_{xy}(f) = \frac{\vert S_{xy}(f)\vert ^2}{S_{xx}(f) \times S_{yy}(f)} \]

where:
\(S_{xy}(f)\) is the cross-spectral density of \(x\) and \(y\),
\(S_{xx}(f)\) and \(S_{yy}(f)\) are the auto-spectral densities (power spectra) of \(x\) and \(y\) respectively.

Key properties:

  • Values range from 0 (no consistent relationship) to 1 (perfectly consistent phase relationship at that frequency)
  • Coherence of 1 at a given frequency means the two signals maintain a constant phase difference at that frequency across all observations
  • Coherence is symmetric - \(C_{xy}(f) = C_{yx}(f)\) - and therefore cannot reveal directionality
  • Like covariance, coherence requires averaging across trials or time segments to produce a stable estimate; a single trial will always yield coherence of 1

Magnitude-squared coherence (MSC) - the squared absolute value of coherence, as defined above, emphasizing strong relationships

Both coherence and magnitude-squared coherence (MSC) measure the linear relationship between two signals in the frequency domain - differing primarily in whether phase information is retained.

Coherence is complex-valued, encoding both the magnitude and the phase of the relationship between two signals at each frequency:

\[ C_{xy}(f) = \frac{P_{xy}(f)}{\sqrt{P_{xx}(f) \times P_{yy}(f)}} \]

where:
\(P_{xy}\) is the cross-power spectral density,
\(P_{xx}\), \(P_{yy}\) are the auto-power spectral densities of \(x\) and \(y\).

  • Magnitude \(\vert C_{xy} \vert\) - strength of the linear relationship (0 = none, 1 = perfect)
  • Phase \(\angle C_{xy}\) - phase difference between the signals at each frequency

Magnitude-squared coherence (MSC) discards the phase and retains only the strength of the relationship. It is always a real number in \([0, 1]\):

\[ MSC_{xy}(f) = \frac{\vert P_{xy}(f) \vert^2}{P_{xx}(f) \times P_{yy}(f)} = \vert C_{xy}(f)\vert^2 \]

Interpretation is the same as coherence magnitude: 0 indicates no linear relationship at that frequency, 1 indicates a perfect linear relationship.

Imaginary Coherence (IC) uses only the imaginary component of the cross-spectrum, measures coherence while ignoring zero-lag synchronization (which may be due to volume conduction), making it insensitive to volume-conducted signals that appear with zero phase lag.

Formula:

\[ \text{IC} = \vert \Im(C_{xy}(f)) \vert \]

Interpretation:

  • IC = 0: No synchronization.
  • IC > 0: Significant synchronization at frequency \(f\).

Advantages:

  • Reduces volume conduction effects.
  • Frequency-specific: Can be computed for specific frequency bands.

Which to use:

  • Use coherence when both strength and phase lag/lead are of interest - e.g. for directional connectivity analysis.
  • Use MSC when only the strength of frequency-specific coupling is needed - e.g. for detecting synchronization between channels.

Computing coherence, imaginary part of coherence and magnitude-squared coherence using multi-taper method:

coh_data = NeuroAnalyzer.coherence(
               e10,
               e10,
               ch1 = "Fp1",
               ch2 = "Fp2",
               method = :mt
            )

Three estimation methods are available via the method parameter:

  1. :fft - Fast Fourier transform; fast and straightforward, suitable for stationary signals
  2. :mt - multi-taper; reduces spectral leakage by averaging across multiple orthogonal tapers, producing a more stable estimate
  3. :stft - short-time Fourier transform; computes CPSD in sliding time windows, suitable for non-stationary signals where the frequency content changes over time

Getting phases:

coh_ph = DSP.angle.(coh_data.coh)

Getting magnitudes:

coh = abs.(coh_data.coh)

Plotting coherence between O1 and O2 in the alpha band:

coh_data = NeuroAnalyzer.coherence(
               e10,
               e10,
               ch1 = "O1",
               ch2 = "O2",
               ep1 = 1,
               ep2 = 1,
               flim = band_frq(eeg, band = :alpha)
           )

# magnitude
coh_mag = abs.(coh_data.coh)
plot_psd(coh_data.f,
         coh_mag[1, :, 1],
         ylabel = "",
         title = "Coherence (magnitude)")
plot_psd(coh_data.f,
         coh_data.msc[1, :, 1],
         ylabel = "",
         title = "MSC")
plot_psd(coh_data.f,
         coh_data.imcoh[1, :, 1],
         ylabel = "",
         title = "Coherence (imaginary part)")

Also, a dedicated plot_coherence() function is available:

coh_data = NeuroAnalyzer.coherence(
               eeg,
               eeg,
               ch1 = get_channel(eeg, type = "eeg")[1:4],
               ch2 = get_channel(eeg, type = "eeg")[5:8],
               ep1 = 1,
               ep2 = 1,
               flim = band_frq(eeg, band = :alpha)
           )

# get magnitude and simplify values
coh, _ = areduce(abs.(coh_data.coh), coh_data.f, n = 0.1)
mscoh, f = areduce(coh_data.msc, coh_data.f, n = 0.1)

# generate label pairs
l = paired_labels(get_channel(eeg, type = "eeg")[1:4],
                  get_channel(eeg, type = "eeg")[5:8])

# plot coherence values along frequencies
plot_coherence(coh[1, :, 1],
               f,
               title = "Coherence $(l[1])")

For multi-channel data:

plot_coherence(mscoh[:, :, 1],
               f,
               clabels = l,
               title = "Magnitude-squared coherence")
plot_coherence(mscoh[:, :, 1],
               f,
               clabels = l,
               avg = true,
               title = "Averaged magnitude-squared coherence")
plot_coherence(mscoh[:, :, 1],
               f,
               clabels = l,
               ci95 = true,
               title = "Magnitude-squared coherence")

Plotting averaged coherence values as connections:

c = zeros(19, 19)
coh_avg = round.(mean(coh[:, :, 1], dims = 2), digits = 2)
c[1, 5] = coh_avg[1]
c[2, 6] = coh_avg[2]
c[3, 7] = coh_avg[3]
c[4, 8] = coh_avg[4]
plot_locs(eeg,
          ch = "eeg",
          connections = c,
          gui = false)

Tip: Plotting connections requires channel locations.

Tip: Line thickness represents connectivity value. If mono=true, positive values are represented with gray lines, negative - with dotted gray lines. If mono=false, positive values are represented with red lines, negative - with blue lines.

To present values numerically, set the option weights to false.

plot_locs(eeg,
          ch = "eeg",
          ch_labels = false,
          connections = c,
          weights = false,
          gui = false)

Tip: If mono=false, positive values are in red and negative values are in blue.