The sequence of analyses we will go through shows the matlab code that answers the following questions

• Extracting instantaneous phases and amplitude envelopes of LFPs using the Hilbert transform.
• Computing the single trials Modulation Index (Tort, 2010), and the Weighted Phase Locking Factor (Maris, 2011), across a broad range of frequencies
• Calculating single LFP pair surrogate statistics to identify significant CFC.
• More detailed analysis of the phase at which high frequency activity bursts synchronize.

The code that we have developed for this tutorial has been provided online under a GNU general public license. To continue with this tutorial, please download the code, and add it to your matlab search path.

This tutorial also relies on the hilbert function in the matlab signal processing toolbox (or alternatively, you can implement the hilbert function for octave).

The CircStat toolbox from Berens is also necessary for analyses related to the phase.

Throughout this tutorial, you will see reference to “filtereddata”. This is the data that has been bandpass filtered for different frequencies, and is of the form [time x frequency]. You can simulate raw data using ft_freqsimulation, and perform bandpass filtering with ft_preprocessing (or the more low-level function, ft_preproc_bandpassfilter). It is suggested to use a variable band-pass that is +/- 1/3 Hz of the centre frequency (see also, Voloh et al 2015, PNAS).

The MI is based on the Kullback-Leibler distance between two distributions. For cross frequency analysis, the empirical distribution is the amplitude distribution. As defined by Tort 2010, the distribution is divided into phase bins, each bin is filled with the instantaneous amplitude of one signal wherever the phase in another signal corresponds the phase bin. The value of each bin is the average of all entries in the bin. This amplitude distribution becomes a quasi-probabilistic distribution by normalizing each bin entry by the sum over all bins. The MI is the the Kullbeck-Lieber distance between this empirical distribution, and a hypothesized uniform distribution, normalized such that values of MI are bound between 0 and 1.

To calculate the MI, two pieces of information are mandatory; the instantaneous phase of one signal and the instantaneous amplitude of another signal (of course, the phase and amplitude may be derived from the same signal). The instantaneous phase and amplitude information can be extracted by first converting the signal x to an analytic signal. This can be achieved with the short-time Fourier transform, wavelet transform, or Hilbert transform. We will use the Hilbert transform in this tutorial. Applying the hilbert function to the filtered signal x, we a complex-valued analytic signal Hx.

The instantaneous phase of Hx can be derived with the angle function of matlab, while the amplitude envelope is derived using the abs function.

% assign data
x=filtereddata;

% get the analytic signal
Hx=hilbert(x);

% extract phase and amplitude info
xphase=angle(Hx);
xamp=abs(Hx);

With these signals in hand, we can call the main function that calculates the MI, get_mi.m . This function also requires the specification of the number of phase bins (nbins) that will be used to segregate the data and build the amplitude distribution.

This function (and accessory functions that perform the calculations) have been optimized to run for 2 or 3 dimensional phase and amplitude signals. The first dimensions must correspond to the time ( thus, concatenating trials in order to calculate a multi-trial MI must be done along the first dimension). The second dimension corresponds to frequency. In other words, each column of the phase and amplitude matrices is derived from the signal from a single channel that has been filtered for different frequencies. The amplitude and phase matrices can also be 3D, where the third dimension corresponds to independent recording (for example, individual trials).

Information about the expected MI is useful for statistical tests that determine confidence that we have truly observed coupling. This information can be obtained by a surrogate analysis, whereby the data is shuffled in some way and the MI recalculated on the novel, surrogate set. If statistics about the expected level of MI are desired, two more arguments are required to the main function. The first is the number of randomizations (nsurrogates). The second is a number specifying the type of randomization. 3 randomizations have been built into the function (type “ help get_MI ” for more information).

The output of the function is a structure out with 5 fields

• MI: the calculated MI
• MIp: the p-value associated with each MI
• MIsurr: the average MI value of the surrogate sets
• mean_amps: the amplitude distribution for which the Kullbeck-Lieber distance is calculated
• ninds: the number of amplitude entries that were used to calculate the average amplitude for each phase bin

%inputs
nbins=18;
nsurrogates=200;
randtype=2; %timesplice

out = get_mi(xphase,xamp,nbins,nsurrogates,randtype);

The calculation of the wPLF (Maris et al 2011) also depends on the complex-valued signals derived via the hilbert or wavelet transforms (after appropropriate pre-processing and filtering). However, unlike the MI, the wPLF is a complex valued index. It is obtained by taking the inner product of a complex signal with the amplitude envelope of another signal. In assessing cross frequency coupling, the entries of the complex signal carry information about the instantaneous phase (ie in the low frequency) while the amplitude envelope carries information about the instantaneous amplitude (ie in the high frequency signal).

Moreover, each signal is mean centred (ie the mean of the signal is subtracted from each entry) and normalized (ie divided by the signals norm). By applying these operations, the wPLF gives “the explained variance in the mean- centered amplitude envelope […] by the phase information in the normalized” complex signal (Maris et al 2011). The wPLF as defined by Maris is the average inner product taken across individual realizations (ie trials) of the experiment.

% 
% assign data
x=filtereddata;

% get the analytic signal
Hx=hilbert(x)

%get the complex- and real- valued signals
xphase=Hx;
xamp=abs(Hx);

The main function that calculates the wPLF and performs statistical testing is get_wplf.m. This functions requires at least two inputs; the complex valued, low frequency signal xphase, and the real-valued, high frequency amplitude envelope xamp. The inputs to this functions must be in a similarly format as inputs to get_mi.m; they must be 2D or 3D arrays, where the first dimension is the time, second dimension is the frequency, and the third (optional) dimension is independent realizations (eg trials) of the signals in the same channel. The bulk of the calculation (ie mean-centering, normalization, computation of inner-product) is performed by calc_wplf.m.

As before, statistical testing randomizes a signal and recomputes the wPLF on the new surrogate data set. By doing this many times, we can compare the observed and surrogate data and thus obtain a p-value. To test for statistical significance, two additional input arguments must be supplied; nsurrogates, which is the number of surrogates data sets that should be generated, and randType, which specifies the type of randomization.

The output of this function is a structure with 3 fields;

• wplf: the calculated wPLF
• wplfp: the p-value associated with individual wPLF entries
• surrwplf: the average wPLF across all surrogate wPLF values

%inputs to the function
nsurrogates=200;
randtype=2; %timesplice

out = get_wplf(xphase,xamp,nsurrogates,randtype);

… in progress

In assessing the functional significance of cross frequency coupling, we can use the structure of our task to look for changes in coupling; for example, by quantifying the change in coupling after cue-onset. Normalizing data to baseline activity also ameliorates confounds related to the different sites of recording, and thus allows for comparison between different LFP locations.

The function cfcdiff_wilx tests for a change in coupling between different temporal epochs. It requires two n x m matrices A and B, corresponding to (for example) the pre- and post-cue MI, where n is the phase dimension, and m is the amplitude dimension. This function performs a Wilcoxon signed rank test on the normalized difference between A and B (see the function normdiff) across the population of LFPs to determine frequency pairs that couple. Results are corrected for the false discovery rate (FDR) using the Bonferonni-Holms method. The output is an array stats (the results of the Wilcoxon test for each frequency pair) and p (the FDR corrected p-value).

The results of surrogate statistics can also be applied by supplying logical matrices Ap and Bp that will be sued to mask non-significant coupling to zero.

%difference in MI
d=normdiff(a,b);

%statistical testing
[stat,p]=cfcdiff_wilx(a,b);

Frequency-specific functional effects often localize to specific areas. In this way, oscillations in one area may influence or modulate oscillatory activity in another area. Thus, it is important to determine if a) cross frequency coupling effects localize to particular areas and b) to gauge the degree of connectedness.

To this end, we construct two matrices: a matrix A that details the distribution of all possible LFP pairs, and a matrix B detailing the distribution of just those channels that showed significant coupling. Both A and B must be n x n matrices, where n is the number of areas, the area for phase is along the rows, and the area for amplitude is along the columns. The diagonal of this matrix corresponds to within-area LFP pairs. These matrices are then fed into the function stats_anatomy, performs Pearson’s Chi-square test and a Z-score test on the frequency of coupling between area pairs to determine if there is area-specific coupling.

This function performs the test in two ways: first, on the global distribution of all area pairs; and second, on just the inter-area pairs (i.e. ignoring within-area coupling). We borrow lightly from graph theory in the conceptualization of the inter-area coupling. We essentially quantify the out-degree and in-degree and compare them to an expected value.

%inputs to the function
A=sum(allareas,3);
B=sum(tgcareas,3);

%call function. results displayed in command window.
stats_anatomy(A,B)

Cross frequency coupling assumes that high frequency activity couples to a specific phase of the low frequency oscillation. Knowledge of the preferred phase of coupling can inform of of possible circuits underlying coupling.

The preferred phase can be readily ascertained from the MI and wPLF calculations. In the case of the former, the preferred phase is simply the circular average amplitude distribution (i.e. the circular average of the phase bins, weighted by the average high frequency amplitude). In this calculation, it is also important to consider the influence of differences in the number of data points used in calculating the average amplitude over each phase bin. The output of get_mi feeds directly into the function get_miphase:

miphase = get_miphase(out.mean_amps,out.ninds)

The MI value itself does not differentiate between a unimodal or multimodal amplitude distribution. Thus, it is possible that the circular average is actually a trough of the underlying amplitude distribution. To control for possible multimodal effects, one solution is to find all local maxima in the amplitude distribution and calculate the preferred phase over these peaks, and then compare it to the preferred phase as determined above. The function get_miphase_peaks does exactly that, which is a wrapper for the maltab function findpeaks. Because the amplitude distribution represents circular data, this function first doubles up the amplitude distribution. We must also supply the relative threshold value (i.e. data must be higher than the threshold), and optionally also supply the minimum peaks distance (i.e. don’t find local maxima that are too close together).

threshold=0.7; 
mpd=3;

pksphases = get_miphase_peaks(out.mean_amps, threshold, mpd);

Using the wPLF the preferred phase is the phase of the (complex) wPLF value. The get_wplf output feed directly into the get_wplfphase function:

wplfphase = get_wplfphase(out.wplf)

Comparing distributions of preferred coupling phases can be done using the CircStat toolbox from Berens. Assessing divergence from uniformity is done via either Rayleigh’s test (circ_rtest, for data sampled from von Mises distributions) or the non-parametric Hodjes-Ajne test (circ_otest, for data sampled from a NON von Mises distributions). Comparison of distributions (eg between correct and error trials) can be done via the Kuiper test (circ_kuipertest).

One general proposal for how CFC may be instantiated is when oscillatory activity becomes phase aligned to a specific (endogenous or exogenous) event. If this was consistent across trials, then we should see a substantial rise in phase consistency time-locked to the event.

To quantify this possibility, we can calculate the phase consistency via a Rayleigh’s Z statistic across trials. This is performed via the get_rayz_time function. This function accepts as input a matrix X, which must have at least two dimensions corresponding to the time and number of trials. By default, it is assumed trials are organized along dimension 2, but we can specify the trial dimension along which to calculate the Rayleigh Z with a second input dim.

We can then determine the time at which most LFPs were phase consistent. To do this, we find the peak phase consistency (maximum of Z in the code below) and its corresponding time bin for each LFP. This allows us to perform Pearson’s Chi-square test to determine if phase consistency of LFPs tends to cluster in time. The function rayz2zpks performs this analysis.

Note that get_rayz_time must be embedded in a script to calculate the Rayleigh Z for individual LFPs, and it’s outputs concatenated into a matrix of the form N x t, where N is the number of LFPs and t is the time.

X=lfpphases; % this has the dimensions LFP x time x trial
dim=3;

% get time-resolved Rayeligh Z score
[P, Z] = get_rayz_time(X,dim); %Z and P have the appropriate form to input to rayz2zpks

% assess significance of temporal clustering of phase consistent LFPs
time=-500:1:500; %ms
bins=10;

[stats, edges, bincentre] = ray2zpks(Z,time,bins);