  The Neurophysiological Biomarker Toolbox (NBT) Phase Locking Value

Coordinated activation over different brain regions can be measured using the statistical method which goes under the name of phase synchronization.
Synchronization phenomena have been widely explore in chaotic systems and in the field of nonlinear dynamics.
Applied to biological time series (irregular, non-stationary, non-linear and noisy) of the brain electrical activities, phase synchronization detection can provide an indication of short-range synchronies, commonly interpreted as subserving 'perceptual binding' among adjacent or same brain regions, and more interesting it can describe long-range synchronization patterns between widely separated brain regions, thought to subserve cognitive mechanisms, such as memory, emotion, motor planning.

Related: Phase Lag Index

Phase synchronization differs from coherence measure: the latter does not separate the effects of amplitude and phase in the interrelations between two signals, while in the former the phase component is obtained separately from the amplitude component for a given frequency

This page is divided into two sections:

Extended Phase Synchrony tutorial

Two coupled frictionless harmonic oscillators, $a(t)$ and $b(t)$, are in phase synchrony when the PHASE LOCKING relationship is verified:

$φ_{n,m} = nΦ_{a}(t)-mΦ_{b}(t) = constant$

where n and m are small integers that define the frequency equality $nΦ_{a}(t) = mΦ_{b}(t)$ of the coupled slow and fast oscillations, $Φ_{a}(t)$ and $Φ_{b}(t)$ are the phases of the two oscillators and $φ_{n,m}$ is the relative phase.
In noisy and/or chaotic systems, the phase locking condition is replaced by the weaker PHASE ENTRAINMENT condition:

$|φ_{n,m}| = |nΦ_{a}(t)-mΦ_{b}(t)| < constant$

or by the even weaker condition of FREQUENCY LOCKING:

$<ω_{n,m}> = n<ω_{a}(t)>-m<ω_{b}(t)> = 0$

where <> denotes averaging over time, $<ω_{n,m}>$ is the relative frequency, $ω_{a}(t)$ and $ω_{b}(t)$ are the frequency of the oscillators expressed as time derivative of the phase of the two oscillators, respectively.

In biological signal, such as EEG time series, synchronization is measured computing the PHASE LOCKING VALUE, which is the mean phase coherence of an angular distribution:

$R = |1/N ∑_{j = 0} ^{N-1} e^{iφ_{n,m}}| = 1-CV$

where $φ_{n,m}$ is the relative phase, N is the number of samples of the data set, CV denotes the circular variance of an angular distribution obtained by transforming the relative phase angles onto the unit circle in the complex plane.

R has values in [0 1], R reaches the value 1 if and only if the condition of strict phase locking is obeyed, whereas for a uniform distribution of phases (for example in unsynchronized time series) R = 0.

Instantaneous Phase is computed using Hilbert transform (see link for definition). Since this requires integration over infinite time, which cannot we performed for a finite length data set, 10% of the calculated instantaneous phase values should be discarded on each side of the dataset. Instantaneous phase difference can be computed using either wavelet analysis or Hilbert transformation. Studies show that these two approaches are equivalent for the analysis of EEG signals.

Rossler Coupled Systems

For illustration, here, we consider a commonly used example with two coupled (1,2) non identical Rossler systems subjected to noisy perturbations:

$\dot{x}_{1,2} = -ω_{1,2}y_{1,2} - z_{1,2} + ξ_{1,2} + ε(x_{2,1}-x_{1,2})$

$\dot{y}_{1,2} = ω_{1,2}x_{1,2} + ay_{1,2}$

$\dot{z}_{1,2} = b + z_{1,2} (x_{1,2}-c)$

where a, b, c are constant parameters of the Rossler systems, $ω_{1,2}$ are the natural frequencies of the two systems, $ε$ governs the strength of the coupling, $ξ_{1,2}$ are two Gaussian delta-correlated noise.

The code in this link solves the coupled system of equations, extracts the instantaneous phase of the two oscillators computes the relative phase when the systems are uncoupled ($ε=0$, and not affected by noise), coupled ($ε=1$, and not affected by noise) and affected by noise.

Indeed, as you see in the figure below, increasing the coupling strength the systems are more and more synchronized resulting in small oscillations of relative phase around a constant value also in presence of noise.

Lag Synchronization

Interestingly, it is possible to identify LAG SYNCHRONIZATION, meaning that phase synchronization between the system appears at specific time lag, τ. This time lag is obtained my computing the minimum of a Similarity Function defined as follow:

$Sˆ2(τ) = \frac {<[x_{2}(t+τ)-x_{1}(t)]ˆ2>}{[]ˆ1/2}$

$σ = min_{τ}(S(τ))$.

If $x_{1} = x_{2}$, $σ = 0$ for $τ=0$, and COMPLETE SYNCHRONIZATION is verified, and the systems are perfectly coupled.

If the similarity function has a minimum for nonzero time shift, it means that a time lag exists between the two processes.

If $x_{1}$ and $x_{2}$ are completely independent $S(τ)∼1$ for all time shifts $τ$.

In the following figure you can see how increasing the coupling strength not only the mean relative frequency reduces, but also the similarity function has minimum in correspondence of smaller time lag.

You can obtain this plot using the code in this link , where we vary the coupling factor from 0 to 0.4, for two coupled nonidentical Rossler systems.

Phase Locking Value in simulated EEG signals

Considering EEG time series correlated non linear oscillators, we simulate two real signals $u$ and $v$, linear combination of the components $x_{1}$ and $x_{2}$ of two uncoupled Rossler systems affected by noise ($ε=0$).

mu = 0.02;
u = mu*x2+(1-mu)*x1;
v = mu*x1+(1-mu)*x2;

We compute the instantaneous phases using Hilbert transform:

h1=hilbert(u);
h2=hilbert(v);
[phase1]=unwrap(angle(h1));
[phase2]=unwrap(angle(h2));
% since the calculation of the hilbert transform requires integration over
% infinite time; 10% of the calculated instantaneous values are discarded
% on each side of every window
perc10w =  floor(nit/10);
phase1 = phase1(perc10w:end-perc10w);
phase2 = phase2(perc10w:end-perc10w);

We compute the Phase Locking Value (PLV):

n = 1;
m = 1;
RP=n*phase1-m*phase2; % relative phase
P_L_V=abs(sum(exp(i*RP))/length(RP));

In the following figure, we show how increasing the coupling strength of the two systems the PLV approaches 1.

n:m Synchronization Indices

Until now we made the hypothesis of locking ratio equal to $n:m = 1:1$, which is plausible in EEG time series given that signals are obtained from the same physiological system (i.e. the brain). In some chases you may want to examine synchronization mechanisms measured at different districts (i.e. synchronization between MEG and Electromiographic signals ) or you want to be sure that the assumption n:m = 1:1 is correct for every couple of EEG time series. In order to determine the n:m ratio synchronization indices are introduced. n:m synchronization indices are used also to characterize the strength of synchronization.

Here we present three main indexes that are implemented in the function nbt_n_m_detection . These indices can be computed for different n:m ratio, the n:m ratio that gives the largest indices is, then, selected. You can also use the function nbt_n_m_surrogate_detection where surrogate data are used to derive significance levels for the indices.

1. Index based on Shannon Entropy

It is defined as:

$\tilde{ρ}_{n,m} = \frac{S_{max}-S}{S_{max}}$

where $S = - ∑_{k = 0} ^{N} p_{k} ln p_{k}$ is the entropy of the distribution of the cyclic relative phase $Ψ_{n,m} = φ_{n,m}mod 2π$, $S_{max} = ln N$, where N is the number of bins used for the distribution. The optimum number of bins (N) is set as $e ˆ{(0.626+0.4 ln (L-1))}$, where L is the number of data points.

$0≤ \tilde{ρ}_{n,m}≤1$, where $\tilde{ρ}_{n,m} = 0$ corresponds to a normal distribution (no synchronization), $\tilde{ρ}_{n,m} = 1$ corresponds to a Dirac-like distribution (perfect synchronization). Such distribution can be observed only in the ideal case of phase locking of noise-free quasilinear oscillators.

2. Index based on the intensity of the first Fourier mode of the distribution

It is defined as:

$γ^2_{n,m} = ^2 + ^2$

$0≤γ^2_{n,m}≤1$. The advantage of this index is that its computation involves no parameters: we do not need to choose the number of bins as we do not calculate the distribution itself.

3. Index based on conditional probability

If the oscillators are strongly nonlinear then the distribution of $Ψ_{n,m}(t)$ is non-uniform even in the absence of noise. We observe the phase of the second oscillator ($Φ_{2}$) at the instants of time when the phase of the ﬁrst one attains a ﬁxed value $θ$(phase stroboscope):

$η = Φ_{2} mod 2πn|_{Φ_{1} mod 2πm = θ}$

To account for the n:m locking, the phases are wrapped into intervals $[0 2πm]$ and $[0 2πn]$, respectively. Repeating this procedure for all $0≤θ≤2π$ and averaging, we get a statistically signiﬁcant synchronization index. Practically, if we deal with the time series, we can introduce binning for the phase of the ﬁrst oscillator, i.e. divide the interval $[0 2πm]$ into N bins. Next, we denote the values of $Φ_{1} mod 2πm$ falling into the l-th bin as $θ_{l}$ and the number of points within this bin as $M_{l}$. Then, we compute $M_{l}$ corresponding values $η_{i,l}$, where $i = 1,...,M_{l}$

If the oscillators are not synchronized, then we expect $η_{i,l}$ to be uniformly distributed on the interval $[0 2πn]$, otherwise these quantities group around some value and their distribution is unimodal. To quantify it, we compute:

$Λ_{l} = 1/M_{l} ∑_{i = 1}^{M_{l}} e^{i(η_{i,l}/n)}$

The case of complete dependence between both phases corresponds to $|Λ_{l}| = 1$, whereas $|Λ_{l}|$ vanishes if there is no dependence at all. To improve the statistics, we average over all N bins and get the synchronization index.

$λ_{n,m} = 1/N ∑_{i = 1}^{N}|Λ_{l}|$

According to the deﬁnition above $λ_{n,m}$ measures the conditional probability for $Φ_{2}$ to have a certain value provided $Φ_{1}$ is in a certain bin.

Computing Phase Locking Value using NBT

In NBT Phase Locking Value (PLV) is computed by the function PhaseLockingObject = nbt_doPhaseLocking(Signal,SignalInfo,FrequencyBand,interval,filterorder) which updates the Phase Locking Value Biomarker Object. This function computes the PLV for each couple of EEG signals, for a specified frequency band (FrequencyBand) in HZ and a specified time interval in second (interval). You can also specify the order of the fir filter used to extract the desired frequency components (default value is filterorder = 4).

When you set the filter order remember that the signal length must be at least 3 times the filter order. When you set the frequency band remember that the highest of the frequencies must be minor that the half of the sampling frequency

Since nbt_doPhaseLocking uses Hilbert Transform to compute the phases of the signals, for signal longer than 4096/fs seconds (fs is the sampling frequency), the PLV is averaged over a windowed signal (hanning windows of width 4096 samples with 20% overlapping), resulting in a compromise between statistical accuracy and stationarity. This procedure might take more computational time

The command here after computes the PLV biomarker object for the dataset, Signal, SignalInfo, in the alpha frequency range (8-13 Hz), for the first 5 seconds (interval = [0 5]) of the signal. Remember to indicate a directory where to save your phase locking biomarker using the function nbt_SaveClearObject.

PhaseLockingValue8_13Hz = nbt_doPhaseLocking(Signal,SignalInfo,[8 13],[0 5]);
nbt_SaveClearObject('PhaseLockingValue8_13Hz',SignalInfo,SaveDir);

To plot the PLV matrix, type:

PhaseLockingValue8_13Hz.plotPLV The Phase Locking Value Biomarker Object contains information also on the synchronization indices (i.e. PhaseLockingValue8_13Hz.indexE based on Shannon Entropy,PhaseLockingValue8_13Hz.indexF based on intensity of the first Fourier mode, PhaseLockingValue8_13Hz.indexCF based on conditional probability).

You can also create and run the function nbt_runPhaseLocking(Signal, SignalInfo, SaveDir) to process and save more than one PLV biomarker at time, as follow in this example code

tutorial/phase_locking_value.txt · Last modified: 2014/11/20 23:04 by Simon-Shlomo Poil
The NBTwiki platform - version 2.8 - 9 May 2013