|
Size: 22551
Comment:
|
← Revision 292 as of 2023-04-19 19:00:56 ⇥
Size: 34293
Comment:
|
| Deletions are marked like this. | Additions are marked like this. |
| Line 1: | Line 1: |
| ## page was renamed from \TutPac | |
| Line 3: | Line 2: |
| This tutorial introduces the concept of phase-amplitude coupling (PAC) and the metrics used in Brainstorm to estimate it. Those tools are illustrated on three types of data: simulated recordings, rat intra-cranial recordings and MEG signals. == Phase-amplitude coupling == Illustrated introduction and mathematical background. == Simulated recordings == Step-by-step instructions with as many screen captures as possible: generation and analysis of the signals. == Rat recordings == How to download the data. Step-by-step instructions to analyze the recordings. == MEG recordings (CURRENTLY BEING UPDATED - UNFINISHED) == Step-by-step instructions to analyze the wMNE source signals for Phase Amplitude Coupling. In order to do this part of the tutorial you will need to get the file sample_resting.zip from the [[http://neuroimage.usc.edu/bst/download.php|Download]] page. Preparation of the anatomy, basic pre-processing and source modeling will be only mentioned briefly and will be similar to the continuous recordings tutorials found here: [[http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawViewer|Continuous Recordings Tutorial]] === Step 1: Pre-processing === The first steps include importing the anatomy and the functional data and projecting the sources. If you you require more information than the brief overview provided here detailed description and guidance for all the steps can be found in the Continuous Recordings tutorial or within the tutorials for the '12 Easy steps for Brainstorm', all of which are available from this page: [[http://neuroimage.usc.edu/brainstorm/Tutorials|Tutorials]] Before doing Phase amplitude coupling analysis we need the pre-processed files to analyze. We need to: * Start a new protocol * Import the anatomical and functional data * Do some basic pre-processing (artifact correction) * Project sources ==== First steps ==== * Create a new protocol in Brainstorm by going to file/newProtocol and call it 'PACTutorial' * Create a new subject by going to file - new Subject and label it 'PACTutorialSubj1'<<BR>> * Select "Use individual anatomy" for default anatomy * Select "Use one channel file per condition" for default channel file ==== Anatomy ==== * Click on the anatomy tab at the top left of Brainstorm {{attachment:AnatMenu.gif}} * Right click on 'PACTutorialSubj1' and click on Import anatomy folder * You now need to select the 'Anatomy' folder inside the sample_resting folder you should have downloaded. Click 'open' when you have selected this folder.\ * You now need to define the fiducial points * The MRI coordinates should be (+/- a few millimeters): . NAS: x=128, y=225, z=135 . LPA: x=54, y=115, z=107 . RPA: x=204, y=115, z=99 . AC: x=133, y=137, z=152 . PC: x=132, y=108, z=150 . IH: x=133, y=163, z=196 (anywhere on the midsagittal plane) ==== Functional data ==== * The sample_resting download contains two 10 minute resting state runs. We are going to use the first one which is the one labelled 'subj002_spontaneous_20111102_02_AUX.ds'. * Click on the functional data tab, the middle button, directly to the left of the anatomy button * Right click on the 'PACTutorialSubj1' and click on 'review raw file' * This creates a link so Brainstorm can read the raw file * It will give you an option to 'refine using head points' - select 'Yes' ==== Pre-Processing ==== All data should be pre-processed and checked for artifacts prior to doing analyses such as PAC (including marking bad segments, and correcting for artifacts such as eye blinks and heartbeats with SSPs). For the purposes of this tutorial, we will artifact correct with SSPs but will not go through marking out bad sections. When using your own data reviewing the raw data for bad sections and using clean data is of the utmost importance. SSPs for cardiac and eye artifact; * Open the 'Link to raw file'. From the SSP menu select 'Detect eye blinks'. Put this as the Channel name and * The EOG channel in this data is labelled 'EEG058' in this data file. Put this as the 'Channel name:' and click run. <<BR>> {{attachment:blink.gif}} * Repeat same procedure for 'Detect heartbeats' using 'EEG057' which is the channel name of the ECG in this data set. . * Then use the ‘Compute SSP: Eyeblinks’ and ‘Compute SSP: Heartbeats’ to project away these artifacts from the data. Make sure to write 'MEG' in the 'Sensor Types or Names' option box if it is not already. For consistency with this tutorial use (only) the first component for each SSP. * For more information regarding dealing with artifacts and SSPs, view this tutorial: [[http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawSsp|Artifact Tutorial]] '''__Regarding sin removal:__''' PAC analysis involves examining a very wide band of frequencies, often the examining the entire range of 2Hz - 150Hz or more. This band contains the frequencies contaminated by line noise, of either 50 or 60 Hz and their harmonics. Here we will not do sin removal in part because it would too long for the sake of a tutorial but also because it is reasonable the not perform it. The PAC function looks for high frequencies occuring specifically certain phases of low signals such that the ubiquitous nature of line contamination effectively cancels it out for being identified as PAC. (Similarly, doing sin removal results in no 60 Hz anywhere, such that the function also identifies no PAC). To demonstrate that we can safely proceed without sin removal, consider the following PAC maps performed on the same time series with the only difference being lin noise removal (60 Hz and 120Hz). The one on the left is the raw signal and the one on the right had sin removal performed. {{attachment:SinCheck.gif}} Best practice regarding line noise removal and PAC estimates is: ??? ==== Importing ==== Once the data is pre-processed and ready for further analysis we will now import the data into Brainstorm, project the sources and do the PAC analysis * To import the data right click on the 'Link to raw file' and click on the first option which is 'Import in database' {{attachment:Import.gif}} * You should get the option box as above. We are going to leave all these options at default. When we click 'Import' brainstorm will now create a new file with the data imported with our SSPs applied and DC offset removed. * Note: importing this long recording will create a new large file (~3 gb) and may take a couple minutes. Due to the way Brainstorm reads raw files vs. imported files it is will not it is more computationally demanding to open the imported files (try opening the 'Link to raw file' vs. opening the imported file) which is why we did all the steps in which we need to scroll through the data with the link to the raw file. ==== Project Sources ==== * The imported file should have saved as a new condition in our tree in the brainstorm database. At this point we still have the sensor data and now want to projec the data into source space. We will need a head model and noise covariance matrix (as well as the imported anatomy) in order to do this. * To compute the head model, right click on the newly imported file (which should be labelled 'Raw(0.00s,600.00s) and click on the compute head model. Use the overlapping spheres model and keep all of the options at their default values. * In the original zip download folder there is an empty room recording from when this resting data was collected. Right click on the subject in the database, click 'review raw file' and select this file. Then right click on the 'Link to raw file' of the noise recording, go to the 'noise covariance' submenu and click on 'compute from recordings'. An option box will pop up, within which you should keep all the default values. Click okay to create the noise covariance matrix. This new file should now be available in the tree. Right click on the Noise covariance file and click 'copy to other conditions' to copy this file to all the other conditions were we need it. <<BR>> * Further information as well as the importance and relevance of noise covariance is described here: [[http://neuroimage.usc.edu/brainstorm/Tutorials/TutNoiseCov|Noise Covariance Tutorial]] * You should now have a condition with imported data, a head model and noise covariance that looks like this: PIC HERE * Right click on the raw file again and click 'compute sources'. Use the Minimum norm estimate (wMNE) and keep all the default settings. If you have the first 10 minute files from the resting_sample file projected onto the individual anatomy with one cardiac SSP and one Blink SSP active. We are now ready to run the PAC analysis. === Step 2: Using the PAC Function === Once you have the sources projected onto the anatomy proceed with the following instructions to use the PAC function on the source data. ==== The Function ==== * The function for Phase Amplitude Coupling analysis is found in the frequency menu in the process selection menu. {{attachment:PACMenu.gif}} * Drag and drop the sources file into the dropbox in the process 1 tab. Click on run, go to frequency and click on Phase Amplitude Coupling. ==== Process Options ==== Once you click on 'Phase-amplitude coupling' you should get a pop-up box with the following options. {{attachment:PACOptions65.gif}} . '''Time Window''': The time segment of the input file to be analyzed for PAC. . '''Nesting Frequency Band (low)''': The frequency band of interest for the frequency for phase (the low, nesting frequencies). * This can be a wide exploratory range (2 - 30 Hz) or a much smaller and specific range (ex. theta: 4-8Hz) '''Nested Frequency Band (high)''': The frequency band of interest for the frequencies for amplitude (the high, nested frequencies). * ''' '''The frequency band of interest for the frequencies for amplitude (the high, nested frequencies). This can be a wide exploratory range (ex. 40 - 250 Hz) or a smaller and specifc range (ex. low gamma: 40-80Hz) * Note: The nested frequency can only be as high as your sampling rate has the resolution to yield. '''Processing Options '''(These options should be left at default options unless you know how to use them) . '''Parallel processing toolbox: ''' PAC analysis of each time series (of a vertex or sensor) is independent of the PAC analysis of every other time series. This function is done with a loop and each iteration of the loop is independent of the one previously. The parallel processing toolbox uses a parfor loop in which multiple time series can be processed in parallel. . '''Use Mex files:''' Mex files are available for running this process and contribute to speeding up the computation time. . '''Number of signals to process at once: '''The file that is given to the function is processed in blocks, and this option signifies how many time series are in each block. '''Output Options''' . '''Save average PAC across trials: '''If this box is selected and multiple files have been given to the process, then the process will save the average PAC across all trials in a new file. . '''Save the full PAC maps: '''If this option is selected then the full PAC comodulogram will be saved for each time series (ie - the process saves the values for each and every frequency pairing for each and every time series). If this option is not selected then the process only save the values at the maxPAC - the frequency pairing with the most PAC coupling. * Saving the full PAC maps entails saving (for each time series) a matrix of [NumberOfTimeSeries x 1 x HighFreqs x LowFreqs]. This will very quickly create very large files (in source space it will save a matrix roughly [15000 x 1 x 39 x 12] which is a very large file. Only click this option if you do want to use / look at the maps. Note: The Phase Amplitude Coupling function can be a very computationally expensive and time consuming process. The longer the time segment is and the more time series there are(so more sensors or vertices) the longer this will take. Consider using scouts or downsampling if the computation time is too long. The files that are saved at the end of the process are not dependent on the length of the time series (so if harddrive is an issue using shorter time segments will not help). We will first test the process by computed the PAC for a single vertice. This will allow us to examine what the PAC process does and visualize the result. * In the PAC option box, if the options are not already filled out, fill in time window as the full length of the file (0 - 599.9996) and the frequency options as the wide bands of low: 2 - 30 and high: 40 - 150. * We are going to arbitrarily compute this on source #65. (Vertices in the source data are labelled with numbers). Write '65' in the source indices option. * Click on run. * PAC files that are computed will be saved under the file containing the time series from which they were computed under the name 'MaxPAC'.It should look like this in brainstorm. {{attachment:224DB.gif}} * Double click on the 'MaxPAC file to open the PAC map - the comodulogram. You should open a file which looks like this. {{attachment:comodulogram224.gif}} * We have this map available because we selected the 'Save Full PAC maps' and can therefore now visualize all the frequencies pairings and their PAC strenghs for each time series. * The small white circle indicates the circle indicates the PAC pairing with the strongest coupling and it is the results relevant to this pairing that are displayed on the top of the comodulogram. * You can click to any frequency pairing and the relevant results will be printed at the bottom of the comodulogram.This allows you to explore the values represented in the comodulogram. Here I clicked on PAC at around at 9 Hz nesting and can check that the coupling strength is similar to the maxPAC value of 11.47 Hz. {{attachment:comodulogram224sel.gif}} === Step 3: Verifying with Canolty Maps === Canolty maps are a type of Time Frequency decomposition that offer another way to visualize the data and serve as a complimentary tool to visualize and assess Phase-Amplitude Coupling. DESCRIPTION OF THE PROCESS ==== The Function ==== The Canolty Map's function is also found in the Frequency tab from the process functions. {{attachment:CanoltyMenu.gif}} There are two ways to use Canolty maps - you can manually input a low frequency of interest or you can give it the maxPAC file and it will take the low frequency at the maxPAC value. * Process 1 tab - Drop a file of time series into the process one tab and manually select the low-frequency of interest. * Process 2 tab - Drop a file of time series into File A and the maxPAC file (from the file in A) into File B. This process will make the Canolty maps by finding (for each vertex) the low frequency defined in the maxPAC file and use that to create the Canolty map. We will continue by doing the Process2 version to compliment are maxPAC results. Click on the Process2 tab. In the FileA box drop the original time series (the source data file). In the FileB box drop the maxPAC file that we just created for source. (IMAGE HERE) When you click on the Canolty Maps (process2) function you should get a an options box like this. {{attachment:Canolty2Options.gif}} ==== Process Options ==== . '''Time Window: ''' the time segment from the input file to be used to compute the Canolty map. . '''Epoch Time: '''''' '''The length of the epochs used. . '''Source Indices:''' which time series from the given file to compute the Canolty Maps for. . ''' ''' '''Number of signals to process at once: ''' This process is also done in blocks and this option allows for setting the block size (can be left at the default value). . '''Save averaged low frequency signals:''' In order to create the Canolty maps, Brainstorm filters the input time series at the low frequency of interest. This option saves that filtered signal, which can be useful for viisualization. The only difference in the Process1 version of Canolty Maps is the additional required field of Nesting Frequency. In this case you can enter in any low frequency of interest with which to compute the Canolty Map(s). {{attachment:Canolty1Options.gif}} * First we will run the Canolty map on the single vertice PAC we ran on vertice number 224. * To do this using the Process2 version with the maxPAC file, click on the Process2 button at the bottom of the dropbox and drop the source data into FileA and the MaxPAC file for vertice 224 into FileB {{attachment:224canoltyBox.gif}} * When in the Process2 format, clicking on run will only show the processes available. Canolty maps is still under the frequency tab, but the drop down menu will look a bit different. {{attachment:canolty2menu.gif}} * Click run and the option box for the Canolty process will pop up. In the case where the given maxPAC file does not include PAC values for all the time series in FileA (such as this case where the maxPAC only contains the PAC values for 1 of the vertices in FileA) Brainstorm does not automatically determine which PAC information it has and this information has to be given. * In the 'Source Indices' option write '224' so that it will use the maxPAC for this vertice and compute the Canolty map. {{attachment:canolty224Options.gif}} * Canolty maps are saved under the time series file from which they come - similar to the maxPAC. However, if the lowFreq signal is saved the link to this file is saved below. In Brainstorm your database should look something like this. {{attachment:canotlyTree224.gif}} * Double click on the 'Canolty map' file to open it. You should see the following image. At the top the image the vertice number and low freq are written. {{attachment:Canolty2-224.gif}} * Here we can see that the Canolty map corroborates what was represented in the maxPAC file. With the data plotted controlling the phase of the low freq (here - 11.47Hz) we can see that the amplitude of the gamma (indicated by the colour) is patterned such that it appears to * We can also verify the frequencies of the nested frequencies. Similar to our PAC comodulogram, we can see in the Canolty map that the nested frequency is predominantly around 80 Hz, but that we see (less) PAC occuring at other frequencies, such as around 120Hz. * In the Canolty map itself there is no representation of the lowFrequncy used. It can be useful to visualize the low frequency, especially if the Canolty maps look a bit 'messier'. The low frequency is accessble through the other link in the brainstorm database. Double click on the file named 'Canolty ERP'. This will open the time series filtered at the low frequency of interest that was used to compute the canolty map. {{attachment:CanoltyERPLabel.gif}} * TEXT {{attachment:C-ERP-224.gif}} * By arranging the Canolty map and Low Frequency file you can get a sense of how the low frequency oscillation and high frequency amplitude relate to each. The time selection on the two is synchronized, so try clicking at particular parts of the low frequncy file and examining the power in the Canolty map. You should notice gamma amplitude is low at the peak and trough of the oscillation and high between them. * There are other more quantitative ways of verifying the phase of the coupling (such as by using the value extracted by the maxPAC function. This is for visualization purposes, a sanity check and converging evidence of the coupling. {{attachment:Canolty+ERP224.gif|Canolty+ERP224.gif}} xx * You may remember that in the PAC comodulogram for vertice 224 the maxPAC value was at 11.47 Hz but that there was also other areas of high PAC, including an almost equal coupling intensity at 8.3 Hz. * Canolty maps only portray information relevant to the low frequency used to create the map - therefore we cannot make any conclusions about PAC at lowFreq = 8.3 with the Canolty map we have made with lowFreq = 11.47 Hz. * We will now examinethe PAC at lowFreq = 8.3 with a new Canolty map using the Process1 version. * Since 8.3 Hz is not the low frequency at the maxPAC pairing in the maxPAC pair we cannot examine this by giving the maxPAC file, we must manually specify it as a low frequency of interest. * Go back to Process1 and drop the source time series. {{attachment:Process1canolty.gif}} * Press run, go to frequency and click on the Canolty maps process * We need to specify the vertice of interest (224) and the low frequency of interest (8.3 Hz). Fill out the option box as follows. {{attachment:Canolty1-8.3Options.gif}} * Press run * Open the resultant file and you should see the following Canolty map * Again, we can see that when filtered based on the low frequency of 8.3 Hz the high frequency amplitude rises and falls in a consistent manner, supporting that there is indeed PAC. {{attachment:canolty-224-8.3.gif}} * We can also visualize the relation by using the low frequency filtered signal that we saved again. {{attachment:Canotly224-8.3+ERP.gif|Canotly224-8.3+ERP.gif}} An alternative check we can use Canolty maps is to verify that in the case where the PAC function indicates very low levels of Phase-amplitude coupling, that the Canolty map function also corroborates this. * Open the PAC map for vertice #224. Now we want to find a lowFreq where the PAC function did not indicate much coupling. IMAGE HERE: 224 - PACmap 3.57 {{attachment:224-3.57.gif}} * Take the lowFreq of 3.57. * Now use the lowFreq value of 3.57 to create a Canolty map, in the exact same way we did with the lowFreq of 8.30, changing only this value in the parameters. * You should get a Canolty map that looks like this IMAGE HERE: 224C - 3.57 {{attachment:224C-3.57.gif}} * Here we can see that the Canolty map displays nothing that looks like consistent coupling between the low frequency of 3.57 Hz and any high frequencies. This is what we expected based on the PAC map. |
''Authors: Soheila Samiee, Thomas Donoghue, Francois Tadel, [[https://www.neurospeed-bailletlab.org/sylvain-baillet|Sylvain Baillet]]'' This tutorial introduces the concept of phase-amplitude coupling (PAC) and the metrics used in Brainstorm to estimate it. These tools are illustrated first on simulated recordings. In separate tutorials, we illustrate how to use them on [[Tutorials/Resting|resting-state MEG recordings]]. <<TableOfContents(2,2)>> == Introduction == ==== Phase-amplitude coupling (PAC) ==== The oscillatory activity in multiple frequency bands is observed in different levels of organization from micro-scale to meso-scale and macro-scale. Studies have been shown that brain functions are achieved with simultaneous oscillations in different frequency bands [Schutter and Knyazev, 2012]. Classical studies in this field were only focused on rhythms in each of these frequency bands, and it has been reported that these rhythms are linked to perception and cognition [Cohen, 2008]. However, it is revealed that not only examining brain activity in each single frequency band, but also the relation and interaction between oscillations in different bands, can be informative in understanding brain function. Thus, this concept increasingly received interest especially in the field of cognitive neuroscience. This interaction between several oscillations is also known as cross-frequency coupling (CFC). Two forms of recognized CFC in brain rhythms are: phase amplitude coupling (PAC), and phase-phase coupling (PPC). In the first type, which is also called nested oscillations, the phase of the lower frequency oscillation (nesting) drives the power of the coupled higher frequency oscillation (nested), that results in synchronization of amplitude envelope of faster rhythms with the phase of slower rhythms. The second form is amplitude independent phase locking between ''n'' cycles of high frequency oscillation and ''m'' cycles of low frequency one. That's why it is also called ''n:m'' phase synchrony [Palva et al., 2005]. Among these two types, phase-amplitude coupling received more interests. It has been shown that behavioral tasks can modulate the phase amplitude coupling [Voytek et al., 2010], and also it is potentially involved in sensory integration, memory process, and attentional selection [Lisman and Idiart, 1995, Lisman, 2005, Schroeder and Lakatos, 2009]. This coupling is observed in several brain regions including hippocampus, basal ganglia, and neocortex; and these observations are reported in rats, mice, sheep, and monkeys, as well as humans [Tort et al., 2010]. Following figure shows a schematic of phase amplitude coupling. In the top signal, we have the sum of a fast and slow oscillations, where the power of fast oscillation's envelope changes with the phase of the slower oscillation, which is a sample of PAC. The bottom signal shows only the fast oscillation and the variation in its power. As it is obvious from comparison of two signals, the fast rhythm's power is always maximum, at a certain phase of slower oscillation, this phase is called coupling phase. {{attachment:pac_schematic2.jpg}} Measures of cross frequency phase amplitude coupling can monitor the relationship between the activities that modulate low frequency oscillations like sensory or motor inputs, and the local cortical activities such as local computations that are correlated to amplitude of higher frequency oscillation [Canolty and Knight, 2010]. All these interesting features of this coupling resulted in proposing several methods for its measuring. Each of these methods has certain limitations and also advantages over the others, and can be used for a particular purpose; that's why no preferred standard method has been chosen for this estimation yet [Tort et al., 2010]. One of these methods, which is implemented in Brainstorm, is called Mean Vector Length (MVL), proposed by Canolty et al. (2006). On the other hand, since phase-amplitude coupling changes over time with interinsic events and external stimuli, time-resolved estimation of this phenomena received increasing interest. Among suggested methods for estimating the dynamics of phase-amplitude coupling, tPAC [Samiee & Baillet, 2017], is implemented in Brainstorm. These methods (MVL and tPAC) and the step by step instruction of using them are described in this section of the tutorial. ==== Two measures of PAC: ==== ==== 1. Time-resolved PAC estimation: tPAC ==== This approach in principle searches for the <<latex($$f_P$$)>> oscillation with strongest PAC to <<latex($$f_A$$)>> bursts, over a time-window, which slides on the input electrophysiological data. Following figure, summarizes the steps of the algorithm (see [Samiee & Baillet, 2017] for more details). ---- {{attachment:Fig1_method_Final.jpg||width="650",height="252"}} ==== 2. Mean Vector Length (Modulation Index) ==== Canolty et al. (2006) pointed out that a time series defined in the complex plane by <<latex($$A_{f_A}. e^{i \phi_{f_p}}$$)>> could be used to extract a phase-amplitude coupling measure. In this formula <<latex($$A_{f_A}$$)>> is the envelope of fast oscillation, and <<latex($$\phi_{f_p}$$)>> is the phase of slow oscillation. Therefore, after filtering in fast and slow oscillation, and extracting the phase of slow, and the amplitude of fast rhythm; each instantaneous fast oscillation amplitude component in time is represented by the length of the complex vector, whereas the slow oscillation phase of the same time point is represented by the vector angle (see following figure). ---- {{attachment:mvl_step3.jpg||width="398",height="362"}} {{attachment:mvl_steps.jpg||width="189",height="177"}} ---- At the absence of phase-amplitude coupling, the plot of the <<latex($$A_{f_A} .e^{i \phi_{f_p}}$$)>> time series in the complex plane is characterized by a roughly uniform circular density of vector points, symmetric around zero, because the <<latex($$A_{f_A}$$)>> values (averaged over cycles of slow oscillation) are approximately the same for all phases. If there is modulation of the <<latex($$f_A$$)>> amplitude by the <<latex($$f_P$$)>> phase, the <<latex($$A_{f_A}$$)>> would be higher at certain phases than others. This higher amplitude for certain angles will lead to a “bump” in the polar plot of the <<latex($$A_{f_A}.e^{i \phi_{f_p}}$$)>>, leading to loss of symmetry around zero. This loss of symmetry can be inferred by measuring the length of the vector obtained from the mean over all points in the complex plane. It is thus assumed that a symmetric distribution as it occurs during lack of coupling leads to a small mean vector length (because the points in the different phases would cancel each other), whereas the existence of coupling leads to a larger mean vector length [Tort et al. 2010]. For more detail read [Canolty et al, 2006]. '''References''' Canolty RT, Edwards E, Dalal SS, Soltani M, Nagarajan SS, Kirsch HE, Berger MS, Barbaro NM, Knight RT (2006). [[http://www.ncbi.nlm.nih.gov/pubmed/16973878|High gamma power is phase-locked to theta oscillations in human neocortex]], Science, 313(5793), 1626-1628. Canolty RT, Knight RT (2010). [[http://www.ncbi.nlm.nih.gov/pubmed/20932795|The functional role of cross-frequency coupling]]. ''Trends in cognitive sciences'', ''14''(11), 506-515 Cohen MX (2008). [[http://www.ncbi.nlm.nih.gov/pubmed/18061683|Assessing transient cross-frequency coupling in EEG data]]. ''Journal of neuroscience methods'', ''168''(2), 494-499. Palva JM, Palva S, Kaila K (2005). [[http://www.ncbi.nlm.nih.gov/pubmed/15829648|Phase synchrony among neuronal oscillations in the human cortex]]. ''The Journal of Neuroscience'', ''25''(15), 3962-3972. Samiee, Soheila, and Sylvain Baillet (2017). [[https://doi.org/10.1016/j.neuroimage.2017.07.051|Time-resolved phase-amplitude coupling in neural oscillations]]. ''NeuroImage'', 59, 270-279.Schutter DJ, Knyazev GG(2012). [[http://www.ncbi.nlm.nih.gov/pubmed/22448078|Cross-frequency coupling of brain oscillations in studying motivation and emotion]]. ''Motivation and emotion'', ''36''(1), 46-54. Tort AB, Komorowski R, Eichenbaum H, Kopell N (2010). [[http://www.ncbi.nlm.nih.gov/pubmed/20463205|Measuring phase-amplitude coupling between neuronal oscillations of different frequencies]]. ''Journal of neurophysiology'', ''104''(2), 1195-1210. Voytek B, Canolty RT, Shestyuk A, Crone NE, Parvizi J, Knight RT (2010). [[http://www.ncbi.nlm.nih.gov/pubmed/21060716|Shifts in gamma phase–amplitude coupling frequency from theta to alpha over posterior cortex during visual tasks]].''Frontiers in human neuroscience'', ''4''. == Simulate signals == For PAC analysis you can use your own data, or generate a synthesized data containing cross-frequency phase-amplitude coupling, with your preferred parameters. Here we explain how to produce simulated data. The model used for data generation is a simple method introduced in [Tort et al. 2010]. In this section, we generate a dataset of synthesized signal and analyze it with available phase-amplitude coupling estimation tool in brainstorm. * Select the menu File > Create new protocol > "'''TutorialPac'''" and select the options: * "'''Yes, use protocol's default anatomy'''", * '''"No, use one channel file per condition'''". * Right-click on the TutorialPac folder > New subject > '''Subject01''' * In the Process1 tab, leave the file list empty and click on the button '''[Run]''' * Select the process "'''Simulate > Simulate PAC signals'''". * Set all the parameters as follows:<<BR>><<BR>> * {{attachment:process_simulate.gif||width="332",height="348"}} * After running the process, you'll get a synthesized data file: <<BR>><<BR>> * {{attachment:pac_file.gif||width="391",height="173"}} * Double-click on it to display the generated signal:<<BR>><<BR>> * {{attachment:pac_sample_synthesized.png||width="502",height="225"}} === Process options === * '''Subject name''': Specify in which subject the simulated file is saved. * '''Condition name''': Specify the folder in which the file is saved (default: "Simulation"). * '''Signal sampling frequency''': Sampling frequency of the simulated signal in Hz. * '''Nested and nesting frequency''': Definition of the phase-amplitude coupling. * '''Coupling intensity''': Value between 0 and 1 which determines how coupled the two oscillations are. * '''Coupling phase''': Phase of slow oscillation (phase driver) where we have the maximum envelope of fast oscillation (high frequency bursts). For example, if <<latex($$\phi$$)>>=90 degree, the maximum coupling happens in the peaks of slow oscillation. * '''Signal-to-noise ratio''': Determines the power of the noise that will be added to the pure coupled oscillations. This noise composed of two parts: a component with 1/f spectrum which simulate the background brains activity, and a white noise which represents the noise of data recording. == PAC estimation with MVL == === Input data === You can use your own data set, or the signal that is generated with brainstorm simulator. Here we used the simulated signal that is generated with the steps explained in the previous section. === PAC estimation === To extract the phase-amplitude coupling from this signal: * Drag the simulated file to the Process1 tab: <<BR>><<BR>> * {{attachment:pac_simul2.png||width="368",height="176"}} * Select the process '''"Frequency > Phase-amplitude coupling'''" * Set the parameters as follows:<<BR>><<BR>> * {{attachment:process_pac.gif}} * The PAC map calculated from your simulated signal is saved in a new file. Double click on it.<<BR>><<BR>> . {{attachment:pac_on_simulated_data.gif||width="500",height="215"}} * This coupling map, also known as comodulogram, shows the cross-frequency coupling in frequency domain. The pseudo color of the map indicates the level of coupling between each pair of low and high frequencies. The abscissa represents the frequencies analyzed as phase frequency, while the frequencies for amplitude are represented in the ordinate axis. * Here, in this example, there is a bump around 6-8 Hz for the phase and 70-85 Hz for the amplitude, which contains the pair that we initially used for synthesizing the data.The parameters of the point with maximum coupling intensity in the map will be shown in the title of the figure. These parameters are the maximum pac level, and the corresponding <<latex($$f_A$$)>> and <<latex($$f_p$$)>>, which are 0.3, 6.14 Hz, and 80.18, respectively Hz in this case. === Process options === __'''Input options'''__: * '''Time window''': The time segment of the input file to be analyzed for PAC. * '''Row names or indices''': The name of channels to be processed in the case of input files containing multiple signals. If empty, it will process all the signals in the file. __'''Estimator options'''__: * '''Nesting frequency band (low)''': Frequencies to consider for slow oscillation (frequency for phase). * Can be a wide exploratory range (2-30Hz) or a smaller specific range (theta: 4-8Hz) * '''Nested frequency band (high)''': Freq to consider for high freq bursts (frequency for amplitude) * Can be a wide exploratory range (40-250 Hz) or a smaller specific range (low gamma: 40-80Hz) * Note: Nested frequencies can only be as high as your sampling rate has the resolution to yield * '''Total number of frequency bins''': Set the total number of frequency bins estimated, for both the nesting and nested frequency bands. __'''Loop options'''__: Should be left at default options unless you know how to use them. * '''Use mex-files:''' Compiled mex-files are available for running this process and contribute to speeding up the computation time. * '''Number of signals to process at once: '''The file that is given to the function is processed in blocks of signals, and this option signifies how many signals are in each block. Increasing this value can make the process a bit faster but requires a large amounts of memory. __'''Output options'''__: * '''Save average PAC across trials''': If this box is selected and multiple files have been given to the process, then the process will save the average PAC across all trials in a new file.<<BR>>Be careful with this option: although it is offered, averaging over frequencies (such as is required to do this) can be problematic and results from this option should be interpreted carefully. * '''Save only the maximum PAC value''': If this option is selected then only the frequency pair with the maximum PAC coupling is saved for each signal in input (maxPAC), instead of the full PAC comodulogram. The default option is to save all the PAC values for each and every frequency pairing for each and every signal. * Saving the full PAC maps entails saving a matrix of [Nsignals x 1 x Nfa x Nfp], this will very quickly create very large files. Check this option if the files generated are too large to display, save or process on your computer. * The maxPAC is not a very stable estimator of the most significant frequency pairing in the comodulograms, so it can be usefull to save all the values for a visual inspection of the full PAC maps. === File contents === The files saved by this process have the same structure as the time-frequency files, with an additional "sPAC" field. To review its contents, right-click on the PAC file > File > View file contents. . {{attachment:file_contents.gif||width="426",height="469"}} Some of the relevent fields in the maxPAC files: * '''TF''': [Nsignals x 1 x 1], Coupling strength of the maximally coupled frequency pairing (maxPAC) * '''Options''': All the parameters that were given to the process * '''sPAC''': All the additional information related to the output of the PAC function * '''DirectPAC''': [Nsignals x Ntime x Nlow x Nhigh], PAC values for all the frequency pairings * '''NestingFreq''': [Nsignals x 1], Low frequency at the maxPAC * '''NestedFreq''': [Nsignals x 1], High frequency at the maxPAC * '''LowFreqs''': [Nlow x 1], List of low frequencies used for the DirectPAC matrix * '''HighFreqs''': [Nhigh x 1], List of high frequencies used for the DirectPAC matrix === Recommendations === ==== PAC estimation using the MVL algorithm ==== In order to have a reliable result from this method it is required to use a signal which length is at least''' ten cycles '''of the slowest oscillation in your low oscillation band. Considering this point is more important in analysis of real databases, where the noise level (and/or background brain activity) can be higher than synthesized data, and the coupling intensity can be low. Thus, if you want to examine the coupling for slow oscillations in [2, 14] Hz, it would be better to use a signal with minimum length of 10 cycles of the slowest oscillation, which would be 10 x 0.5 = 5 S. == Time-resolved PAC estimation with tPAC == === Input data === For tPAC estimation you can use your own data set or simulated data. For illustration here, we simulated three time series with phase-amplitude coupling using [[http://neuroimage.usc.edu/brainstorm/Tutorials/TutPac#Simulate_signals|brainstorm tool]] and combine them. The script to reproduce such signal can be found [[https://github.com/rcassani/brainstorm-scripts/blob/main/tpac_3modes.m|here]]. The resulting time series looks like: {{attachment:simulated_data.png||width="601",height="193"}} During the first half of the signal (10 s duration), the phase of slow oscillation at 9 Hz is coupled to the amplitude of a faster rhythm at 115 Hz. In the second half, the first coupling mode is terminated and the average of two other simultaneous modes appears. These modes are <<latex($$f_{P_2}$$)>> = 13 Hz, <<latex($$f_{A_2}$$)>> = 145 Hz, and <<latex($$f_{P_3}$$)>>=5 Hz, <<latex($$f_{A_3}$$)>> =87 Hz, respectively. The signal-to-noise ratio is set to 6 dB, and the preferred coupling phase in the three modes are 270, 0, and 180°, respectively. The coupling strength is kept constance and identical in all modes. === Time-resolved PAC === To extract the phase-amplitude coupling from this signal: * Drag the simulated file (or real data) to the Process1 tab: . <<BR>> {{attachment:data_in_box.png||width="363",height="158"}} <<BR>> * Select the process '''"Frequency > Time-resolved Phase-amplitude coupling''' '''> tPAC'''" * Set the parameters as follows (these parameters are explained later in the text):<<BR>> {{attachment:tPAC_set_parameters.png||width="450",height="506"}} <<BR>> * The PAC map calculated from your simulated signal is saved in a new file. Double click on it. <<BR>> {{attachment:tPAC_time_fa.png||width="416",height="260"}} <<BR>> * You can smooth the map with activating the smooth display button: . {{attachment:tPAC_time_fa_smooth.png||width="496",height="239"}} * This coupling map, shows the cross-frequency coupling in time-<<latex($$f_A$$)>> domain. The pseudo color of the map indicates the level of coupling in each time-interval for each <<latex($$f_A$$)>>. * Here in this example there is one mode around <<latex($$f_A = 120$$)>> Hz in the first 10 seconds, and two other modes appear in the second half of the signal on around 85 Hz, and 145 Hz, which reflect the oscillations that we initially used for synthesizing the data. * This map is a 2D projection of three dimensional phase-amplitude coupling (time, <<latex($$f_P$$)>>, and <<latex($$f_A$$)>>). For investigating other dimensions you can project it on other planes of <<latex($$f_P - f_A$$)>> and time-<<latex($$f_P$$)>> using available tools in brainstorm (see below). === Process options === __'''Input options'''__: * '''Time window''': The time segment of the input file to be analyzed for time-resolved PAC. * '''Row names or indices''': The name of channels to be processed in the case of input files containing multiple signals. If empty, it will process all the signals in the file. * __'''Buffer'''__: Is 2 seconds of extra data included in ... * '''No:''' zero padding will be used for edge effect reduction * '''yes''': 2 seconds of data (from both sides) will be used as buffer for reducing edge effect in filtering (recommended, especially if the input data is long the margins are not super important) __'''Estimator options'''__: * '''Frequency for phase band (low)''': Frequencies to consider for slow oscillation. It can be a wide exploratory range (2-30Hz) or a smaller specific range (theta: 4-8Hz) * '''Frequency for amplitude band (high)''': Frequencies to consider for high frequency bursts. It can be a wide exploratory range (40-250 Hz) or a smaller specific range (low gamma: 40-80Hz). * Note: The lower bound of the frequency-for-amplitude band must be at least twice the upper bound of the frequency-for-phase band. * '''Frequency resolution for frequency for amplitude''': This value indicates the bandwidth of the subbands in which the band <<latex($$f_A$$)>> will be divided. *'''Length of sliding time-window''': Should at least include one full cycle of minimum <<latex($$f_A$$)>> of interest (entered above). * Decreasing its length will improve temporal resolution, but decrease resolution for frequency for phase, and will also increase the processing time * If you have a event-related design for your study with more than 100 trials, and need very high temporal resolution (e.g. below 5 ms), you can use [[https://doi.org/10.1016/j.neuroimage.2012.09.023|event-related PAC (ERPAC) algorithm]] (not implemented in brainstorm). * '''Number of signals to process at once''': The file that is given to the function is processed in blocks of signals, and this option indicates how many signals will be processed simultaneously. Increasing this value can make the process a bit faster but requires a large amounts of memory. __'''Output options'''__: * '''Save average PAC across trials''': If this checkbox is selected and multiple files have been given to the process, then the process will save the average PAC across all trials in a new file. {{{#!wiki note Be careful with this option: although it is offered, averaging over frequencies (such as is required to do this) can be problematic and results from this option should be interpreted carefully. }}} === File contents === The files saved by this process have the same structure as the time-frequency files, with an additional "sPAC" field. To review its contents, right-click on the PAC file > File > View file contents. {{attachment:tPAC_output_file.png||width="549",height="706"}} '''Note:''' Preferred phase of coupling extracted for each pair of <<latex($$f_P$$)>>-<<latex($$ f_A$$)>> is saved under sPAC.DynamicPhase. It can be exported to matlab and visualized there. __'''Comodulogram of time-resolved PAC'''__: To extract the comodulogram from the time-resolved PAC: * Drag the time-resolved PAC file/files to the Process1 tab: * Note: if you have more than one time-resolved PAC file (e.g. from different trials), it is better to extract on comodulogram from all of them together to increase the SNR (considering the sparse nature of PAC in the three-dimensional time-<<latex($$f_P$$)>>-<<latex($$f_A$$)>> plane * Select the process '''"''''''Frequency > Time-resolved Phase-amplitude coupling''' '''> Extracting comodulogram from tPAC maps''''''"''' * Set the parameters as follows (these parameters are explained later in the text): . <<BR>> {{attachment:tPAC_comod_set_parameters.png||width="420",height="257"}} <<BR>> * The PAC map calculated from your simulated signal is saved in a new file. Double click on it. . <<BR>> {{attachment:tPAC_comod.png||width="320",height="251"}} <<BR>> * This coupling map, shows the cross-frequency coupling in <<latex($$f_P$$)>>-<<latex($$f_A$$)>> domain. The pseudo color of the map indicates the level of coupling for each pair of frequencies. * Here in this example there is one mode arouond <<latex($$f_A = 120$$)>> Hz in the first 10 seconds, and two other modes apear in the second half of the signal on around 85 Hz, and 145 Hz, which reflect the oscillations that we initially used for synthesizing the data. === Process Options === __'''Input options'''__: * '''Time window''': The time segment of the input file to be analyzed for extracting comodulogram. * '''Time-resolved Comodulogram''': * Input should be in second * It is recommended to leave it as default: 0 -- extract one comodulogram from the total file * Not recommended: If you would like to see the movie of time resolved comodulogram (in current set up of the brainstorm it can only be played in matlab and only once when it is computed), enter the window length to estimate several comodulograms each with that length from input tPAC map. -- You can zoom in/out in played movie in matlab. __'''Method'''__: * Method for estimating comodulogram. The method will be the same for both options if you have only one input (like this example with simulated data. * '''All sources/channels together''': select this option if all sources belong to the same scout, or you do not expect different responses from them. (This will increase the SNR for comodulogram estimation, and make it more accurate and less noisy) * '''Each source/channel separately: ''' . {{attachment:move_between_comod.png}} __'''Time-resolved PAC in '''____'''time - <<latex($$f_P$$)>>'''__'''__ domain__''' ... == Further analysis of PAC files == All functions for further analysis of PAC files (e.g., deriving the average of PAC maps, producing comodulograms, running a z-score transform with respect to baseline levels or surrogate data, etc. ) are available as processed located in Frequency > Time-resolved Phase-amplitude coupling. {{attachment:PAC_further_analysis.png||width="420",height="257"}} '''Note: '''Basic functions for averaging or running statistics are not yet avalaible with Brainstorm [working on it...]. Yet, PAC output files are stored as Matlab files and contents can be readily processed with custom scripts of your own. Yet, surrogate analysis is only available for single band tPAC maps, and can be accessed here: {{attachment:PAC_surrogate.png||width="420",height="257"}} Further, tPAC traces can be z-scored with respect to surrogate data using the Process2 tab. == Practical suggestions for PAC analysis == How to best get started with PAC analysis? How do you decide about setting parameters? Here, we provide basic recommendations that we hope can be helpful: === General suggestions === * First, decide which frequency bands of interest for the frequency for phase (fP) and the frequency for amplitude (fA). These bands are usually hypothesized from the neurophysiological mechanisms you wish to study. If you have no previous clue from previous literature etc. you may want to start with visualizing the power spectrum density (PSD) of the empirical signals. It is recommended to define the fP band around a visible peak of the PSD [Aru et al. 2015], which indicates strong rhythmic signalling. * If you are not sure about the fA band, which can be elusive in PSD plots, you may first define it as a relatively wide frequency band. The resulting comodulogram output (i.e. a map with <<latex($$f_P$$)>> and <<latex($$f_A$$)>> in horizontal and vertical axis, respectively, and a pseudo color illustrating the intensity of the coupling between each pair of frequencies) may indicate the restricted range of fA where PAC is more strongly expressed, for further single-band PAC analysis (see above). * The intensity of PAC coupling is often stronger for slower fP oscillations. This is due in part to greater signal strengths at lower frequency in electrophysiological data (think 1/f power decrease). We therefore recommend that you subdivide a possibly wide frequency band of interest for <<latex($$f_P$$)>> in to 2 or 3 narrower bands. This approach may also provide insight about possibly co-existing modes of PAC (i.e. having several [fP,fA] pairs coupled with PAC). It also prevents PAC modes at slower fP masking other modes. For example, if you want to investigate coupling between slow rhythms in from delta to beta bands (2-30 Hz) for <<latex($$f_P$$)>> and gamma oscillations for <<latex($$f_A$$)>>, we recommend you run PAC twice: first with delta to theta (2-8 Hz) for <<latex($$f_P$$)>>, then again with alpha to gamma (8-30 Hz). === Time-resolved analysis === For time-resolved analysis, if you are not sure about the exact frequencies that you want to focus on, one practical pipeline is to extract comodulogram with a relatively wide band (based on your initial hypothesis), and then find the main mode of coupling from the comodulogram. In the next step, focus on that particular mode of coupling to evaluate dynamical changes in coupling intensity (or preferred phase of coupling) with a single band analysis (single band option does not provide comodulogram but reflect a trace of dynamical changes in coupling parameters). With this strategy the temporal resolution will be optimized for the main mode of coupling, and also the computational load for intensive analysis (and checking significance of outcomes) will be much lower than looking at the whole interval. Following steps provide more details: 1. Start with comodulogram and find the main mode (or modes) of coupling based on comodulogram. 1. To decide about length of sliding time window in this step, you can start with two cycles of the slowest oscillation in the band of interest for <<latex($$f_P$$)>>. For example, if you are interested in investigating potential PAC between theta to alpha oscillations (4-12 Hz) for <<latex($$f_P$$)>> and gamma rhythms for <<latex($$f_A$$)>>, a conservative short sliding time window can be set to two cycles of 4 Hz (i.e. <<latex($$ 2 \times \frac{1}{4 Hz} = 2 \times 0.25 s = 0.5 s = 500 ms$$)>>. Certainly, if you have long inputs increasing the length of the window would improve the SNR. The only point that you should keep in mind is that, tPAC do sparse estimation of PAC in a 3D space of fP, fA and time; therefore, you need enough number of points in time dimension to be able to retrieve a reasonable projection of it on fp - fA plain (called comodulogram). Hence, the maximum length of the window can be selected based on the input length. For example, in the case where your data is around 5 minutes long, if you set the window length to 20 s, considering the 50% overlap of windows you will end up with 5 * 60s / 10s = 30 points along time dimension (which is used in averaging for comodulogram extraction). It is not recommended to go below a certain number (e.g 30). 1. For extracting comodulogram you can also use none-time-resolved option (e.g. MVL) which use the whole time window and does not divide it into smaller chunks. It should be noted that in one experiment we observed higher resolution from comodulogram extracted using time-resolved method compared to a traditional option (see Fig. 5 in [Samiee & Baillet 2017]: https://www.sciencedirect.com/science/article/pii/S1053811917306195#fig5). The reason behind this improvement could be similar to the improvement in power spectral density estimation using Welch approach compared to original Bartlett's method. 1. If you can find an obvious mode of coupling in the comodulogram, it would be better to investigate that mode of coupling using single band analysis for all regions of interest and time intevals (rather than maxPAC value extracted from comodulogram). 1. If there is no obvious mode of coupling in averaged comodulogram, it is not recommended to go further in analysis. It would be better to take a more close look at modes availble in each single recording, and consistancy of preferred phase of coupling in time and in different trials/subjects. 1. If you observe a blub in comodulogram which looks like a vertical line, you need to further investigate your data to make sure this observed coupling is not a pseudo coupling coming from sharp event artifacts (e.g. spikes in LFP) or asymmetric oscillations (e.g. oscillations with saw tooth waveform shape). See [Cole & Voytek 2017; Aru et al. 2015] for more details. === A quick summary of the most important practical points in PAC analysis === '''1. Temporal resolution and minimum <<latex($$f_P$$)>>''' * Temporal resolution, and frequency resolution are complementary to each other. The minimum <<latex($$f_P$$)>> that can be evaluated is equal 1/sliding time window, which means we need at least one cycle of minimum <<latex($$f_P$$)>> of interest in each sliding time window. For example if your <<latex($$f_P$$)>> band of interest is [4 - 8 Hz], minimum window length is 1/4 = 0.25 s. '''2. Suggestion for deciding about initial <<latex($$f_P$$)>> band''' * If you do not have any priori hypothesis about which frequency band could be involved in the task as <<latex($$f_P$$)>>, you can look at the power spectrum of the original signal. One of the requirements for having a PAC is a clear peak in that power spectrum around <<latex($$f_P^*$$)>>. For example if there is a peak around theta oscillations you can define your band of interest as 4 - 8 Hz. '''3. Coupling strength for the main coupled pair''' * After detection of the best coupled <<latex($$f_A$$)>> and <<latex($$f_P$$)>>, you can run tPAC with single <<latex($$f_A$$)>> band option only around that pair of interest. The window length can also get decreased to include at least two cycles of <<latex($$f_P^*$$)>> (May increase the temporal resolution compared to initial investigation because <<latex($$f_P^*$$)>> could be higher than the minimum <<latex($$f_P$$)>> in our band of interest). The output will be a time-series of PAC strength for that band of interest. <<EmbedContent(http://neuroimage.usc.edu/bst/get_feedback.php?Tutorials/TutPac)>> |
Phase-amplitude coupling
Authors: Soheila Samiee, Thomas Donoghue, Francois Tadel, Sylvain Baillet
This tutorial introduces the concept of phase-amplitude coupling (PAC) and the metrics used in Brainstorm to estimate it. These tools are illustrated first on simulated recordings. In separate tutorials, we illustrate how to use them on resting-state MEG recordings.
Contents
Introduction
Phase-amplitude coupling (PAC)
The oscillatory activity in multiple frequency bands is observed in different levels of organization from micro-scale to meso-scale and macro-scale. Studies have been shown that brain functions are achieved with simultaneous oscillations in different frequency bands [Schutter and Knyazev, 2012]. Classical studies in this field were only focused on rhythms in each of these frequency bands, and it has been reported that these rhythms are linked to perception and cognition [Cohen, 2008]. However, it is revealed that not only examining brain activity in each single frequency band, but also the relation and interaction between oscillations in different bands, can be informative in understanding brain function. Thus, this concept increasingly received interest especially in the field of cognitive neuroscience. This interaction between several oscillations is also known as cross-frequency coupling (CFC).
Two forms of recognized CFC in brain rhythms are: phase amplitude coupling (PAC), and phase-phase coupling (PPC). In the first type, which is also called nested oscillations, the phase of the lower frequency oscillation (nesting) drives the power of the coupled higher frequency oscillation (nested), that results in synchronization of amplitude envelope of faster rhythms with the phase of slower rhythms. The second form is amplitude independent phase locking between n cycles of high frequency oscillation and m cycles of low frequency one. That's why it is also called n:m phase synchrony [Palva et al., 2005].
Among these two types, phase-amplitude coupling received more interests. It has been shown that behavioral tasks can modulate the phase amplitude coupling [Voytek et al., 2010], and also it is potentially involved in sensory integration, memory process, and attentional selection [Lisman and Idiart, 1995, Lisman, 2005, Schroeder and Lakatos, 2009]. This coupling is observed in several brain regions including hippocampus, basal ganglia, and neocortex; and these observations are reported in rats, mice, sheep, and monkeys, as well as humans [Tort et al., 2010].
Following figure shows a schematic of phase amplitude coupling. In the top signal, we have the sum of a fast and slow oscillations, where the power of fast oscillation's envelope changes with the phase of the slower oscillation, which is a sample of PAC. The bottom signal shows only the fast oscillation and the variation in its power. As it is obvious from comparison of two signals, the fast rhythm's power is always maximum, at a certain phase of slower oscillation, this phase is called coupling phase.
Measures of cross frequency phase amplitude coupling can monitor the relationship between the activities that modulate low frequency oscillations like sensory or motor inputs, and the local cortical activities such as local computations that are correlated to amplitude of higher frequency oscillation [Canolty and Knight, 2010].
All these interesting features of this coupling resulted in proposing several methods for its measuring. Each of these methods has certain limitations and also advantages over the others, and can be used for a particular purpose; that's why no preferred standard method has been chosen for this estimation yet [Tort et al., 2010]. One of these methods, which is implemented in Brainstorm, is called Mean Vector Length (MVL), proposed by Canolty et al. (2006).
On the other hand, since phase-amplitude coupling changes over time with interinsic events and external stimuli, time-resolved estimation of this phenomena received increasing interest. Among suggested methods for estimating the dynamics of phase-amplitude coupling, tPAC [Samiee & Baillet, 2017], is implemented in Brainstorm.
These methods (MVL and tPAC) and the step by step instruction of using them are described in this section of the tutorial.
Two measures of PAC:
1. Time-resolved PAC estimation: tPAC
This approach in principle searches for the 
oscillation with strongest PAC to 
bursts, over a time-window, which slides on the input electrophysiological data. Following figure, summarizes the steps of the algorithm (see [Samiee & Baillet, 2017] for more details).
2. Mean Vector Length (Modulation Index)
Canolty et al. (2006) pointed out that a time series defined in the complex plane by 
could be used to extract a phase-amplitude coupling measure. In this formula 
is the envelope of fast oscillation, and 
is the phase of slow oscillation.
Therefore, after filtering in fast and slow oscillation, and extracting the phase of slow, and the amplitude of fast rhythm; each instantaneous fast oscillation amplitude component in time is represented by the length of the complex vector, whereas the slow oscillation phase of the same time point is represented by the vector angle (see following figure).
At the absence of phase-amplitude coupling, the plot of the 
time series in the complex plane is characterized by a roughly uniform circular density of vector points, symmetric around zero, because the values (averaged over cycles of slow oscillation) are approximately the same for all phases. If there is modulation of the amplitude by the phase, the would be higher at certain phases than others. This higher amplitude for certain angles will lead to a “bump” in the polar plot of the 
, leading to loss of symmetry around zero. This loss of symmetry can be inferred by measuring the length of the vector obtained from the mean over all points in the complex plane. It is thus assumed that a symmetric distribution as it occurs during lack of coupling leads to a small mean vector length (because the points in the different phases would cancel each other), whereas the existence of coupling leads to a larger mean vector length [Tort et al. 2010]. For more detail read [Canolty et al, 2006].
References
Canolty RT, Edwards E, Dalal SS, Soltani M, Nagarajan SS, Kirsch HE, Berger MS, Barbaro NM, Knight RT (2006). High gamma power is phase-locked to theta oscillations in human neocortex, Science, 313(5793), 1626-1628.
Canolty RT, Knight RT (2010). The functional role of cross-frequency coupling. Trends in cognitive sciences, 14(11), 506-515
Cohen MX (2008). Assessing transient cross-frequency coupling in EEG data. Journal of neuroscience methods, 168(2), 494-499.
Palva JM, Palva S, Kaila K (2005). Phase synchrony among neuronal oscillations in the human cortex. The Journal of Neuroscience, 25(15), 3962-3972.
Samiee, Soheila, and Sylvain Baillet (2017). Time-resolved phase-amplitude coupling in neural oscillations. NeuroImage, 59, 270-279.Schutter DJ, Knyazev GG(2012). Cross-frequency coupling of brain oscillations in studying motivation and emotion. Motivation and emotion, 36(1), 46-54.
Tort AB, Komorowski R, Eichenbaum H, Kopell N (2010). Measuring phase-amplitude coupling between neuronal oscillations of different frequencies. Journal of neurophysiology, 104(2), 1195-1210.
Voytek B, Canolty RT, Shestyuk A, Crone NE, Parvizi J, Knight RT (2010). Shifts in gamma phase–amplitude coupling frequency from theta to alpha over posterior cortex during visual tasks.Frontiers in human neuroscience, 4.
Simulate signals
For PAC analysis you can use your own data, or generate a synthesized data containing cross-frequency phase-amplitude coupling, with your preferred parameters. Here we explain how to produce simulated data. The model used for data generation is a simple method introduced in [Tort et al. 2010]. In this section, we generate a dataset of synthesized signal and analyze it with available phase-amplitude coupling estimation tool in brainstorm.
Select the menu File > Create new protocol > "TutorialPac" and select the options:
"Yes, use protocol's default anatomy",
"No, use one channel file per condition".
Right-click on the TutorialPac folder > New subject > Subject01
In the Process1 tab, leave the file list empty and click on the button [Run]
Select the process "Simulate > Simulate PAC signals".
Set all the parameters as follows:
After running the process, you'll get a synthesized data file:
Double-click on it to display the generated signal:
Process options
Subject name: Specify in which subject the simulated file is saved.
Condition name: Specify the folder in which the file is saved (default: "Simulation").
Signal sampling frequency: Sampling frequency of the simulated signal in Hz.
Nested and nesting frequency: Definition of the phase-amplitude coupling.
Coupling intensity: Value between 0 and 1 which determines how coupled the two oscillations are.
Coupling phase: Phase of slow oscillation (phase driver) where we have the maximum envelope of fast oscillation (high frequency bursts). For example, if
=90 degree, the maximum coupling happens in the peaks of slow oscillation. Signal-to-noise ratio: Determines the power of the noise that will be added to the pure coupled oscillations. This noise composed of two parts: a component with 1/f spectrum which simulate the background brains activity, and a white noise which represents the noise of data recording.
PAC estimation with MVL
Input data
You can use your own data set, or the signal that is generated with brainstorm simulator. Here we used the simulated signal that is generated with the steps explained in the previous section.
PAC estimation
To extract the phase-amplitude coupling from this signal:
Drag the simulated file to the Process1 tab:
Select the process "Frequency > Phase-amplitude coupling"
Set the parameters as follows:
The PAC map calculated from your simulated signal is saved in a new file. Double click on it.
- This coupling map, also known as comodulogram, shows the cross-frequency coupling in frequency domain. The pseudo color of the map indicates the level of coupling between each pair of low and high frequencies. The abscissa represents the frequencies analyzed as phase frequency, while the frequencies for amplitude are represented in the ordinate axis.
Here, in this example, there is a bump around 6-8 Hz for the phase and 70-85 Hz for the amplitude, which contains the pair that we initially used for synthesizing the data.The parameters of the point with maximum coupling intensity in the map will be shown in the title of the figure. These parameters are the maximum pac level, and the corresponding
and 
, which are 0.3, 6.14 Hz, and 80.18, respectively Hz in this case.
Process options
Input options:
Time window: The time segment of the input file to be analyzed for PAC.
Row names or indices: The name of channels to be processed in the case of input files containing multiple signals. If empty, it will process all the signals in the file.
Estimator options:
Nesting frequency band (low): Frequencies to consider for slow oscillation (frequency for phase).
- Can be a wide exploratory range (2-30Hz) or a smaller specific range (theta: 4-8Hz)
Nested frequency band (high): Freq to consider for high freq bursts (frequency for amplitude)
- Can be a wide exploratory range (40-250 Hz) or a smaller specific range (low gamma: 40-80Hz)
- Note: Nested frequencies can only be as high as your sampling rate has the resolution to yield
Total number of frequency bins: Set the total number of frequency bins estimated, for both the nesting and nested frequency bands.
Loop options: Should be left at default options unless you know how to use them.
Use mex-files: Compiled mex-files are available for running this process and contribute to speeding up the computation time.
Number of signals to process at once: The file that is given to the function is processed in blocks of signals, and this option signifies how many signals are in each block. Increasing this value can make the process a bit faster but requires a large amounts of memory.
Output options:
Save average PAC across trials: If this box is selected and multiple files have been given to the process, then the process will save the average PAC across all trials in a new file.
Be careful with this option: although it is offered, averaging over frequencies (such as is required to do this) can be problematic and results from this option should be interpreted carefully.Save only the maximum PAC value: If this option is selected then only the frequency pair with the maximum PAC coupling is saved for each signal in input (maxPAC), instead of the full PAC comodulogram. The default option is to save all the PAC values for each and every frequency pairing for each and every signal.
- Saving the full PAC maps entails saving a matrix of [Nsignals x 1 x Nfa x Nfp], this will very quickly create very large files. Check this option if the files generated are too large to display, save or process on your computer.
- The maxPAC is not a very stable estimator of the most significant frequency pairing in the comodulograms, so it can be usefull to save all the values for a visual inspection of the full PAC maps.
File contents
The files saved by this process have the same structure as the time-frequency files, with an additional "sPAC" field. To review its contents, right-click on the PAC file > File > View file contents.
Some of the relevent fields in the maxPAC files:
TF: [Nsignals x 1 x 1], Coupling strength of the maximally coupled frequency pairing (maxPAC)
Options: All the parameters that were given to the process
sPAC: All the additional information related to the output of the PAC function
DirectPAC: [Nsignals x Ntime x Nlow x Nhigh], PAC values for all the frequency pairings
NestingFreq: [Nsignals x 1], Low frequency at the maxPAC
NestedFreq: [Nsignals x 1], High frequency at the maxPAC
LowFreqs: [Nlow x 1], List of low frequencies used for the DirectPAC matrix
HighFreqs: [Nhigh x 1], List of high frequencies used for the DirectPAC matrix
Recommendations
PAC estimation using the MVL algorithm
In order to have a reliable result from this method it is required to use a signal which length is at least ten cycles of the slowest oscillation in your low oscillation band.
Considering this point is more important in analysis of real databases, where the noise level (and/or background brain activity) can be higher than synthesized data, and the coupling intensity can be low. Thus, if you want to examine the coupling for slow oscillations in [2, 14] Hz, it would be better to use a signal with minimum length of 10 cycles of the slowest oscillation, which would be 10 x 0.5 = 5 S.
Time-resolved PAC estimation with tPAC
Input data
For tPAC estimation you can use your own data set or simulated data. For illustration here, we simulated three time series with phase-amplitude coupling using brainstorm tool and combine them. The script to reproduce such signal can be found here. The resulting time series looks like:
During the first half of the signal (10 s duration), the phase of slow oscillation at 9 Hz is coupled to the amplitude of a faster rhythm at 115 Hz. In the second half, the first coupling mode is terminated and the average of two other simultaneous modes appears. These modes are 
= 13 Hz, 
= 145 Hz, and 
=5 Hz, 
=87 Hz, respectively. The signal-to-noise ratio is set to 6 dB, and the preferred coupling phase in the three modes are 270, 0, and 180°, respectively. The coupling strength is kept constance and identical in all modes.
Time-resolved PAC
To extract the phase-amplitude coupling from this signal:
- Drag the simulated file (or real data) to the Process1 tab:
Select the process "Frequency > Time-resolved Phase-amplitude coupling > tPAC"
Set the parameters as follows (these parameters are explained later in the text):
- The PAC map calculated from your simulated signal is saved in a new file. Double click on it.
- You can smooth the map with activating the smooth display button:
This coupling map, shows the cross-frequency coupling in time-
domain. The pseudo color of the map indicates the level of coupling in each time-interval for each . Here in this example there is one mode around
Hz in the first 10 seconds, and two other modes appear in the second half of the signal on around 85 Hz, and 145 Hz, which reflect the oscillations that we initially used for synthesizing the data. This map is a 2D projection of three dimensional phase-amplitude coupling (time,
, and ). For investigating other dimensions you can project it on other planes of 
and time- using available tools in brainstorm (see below).
Process options
Input options:
Time window: The time segment of the input file to be analyzed for time-resolved PAC.
Row names or indices: The name of channels to be processed in the case of input files containing multiple signals. If empty, it will process all the signals in the file.
Buffer: Is 2 seconds of extra data included in ...
No: zero padding will be used for edge effect reduction
yes: 2 seconds of data (from both sides) will be used as buffer for reducing edge effect in filtering (recommended, especially if the input data is long the margins are not super important)
Estimator options:
Frequency for phase band (low): Frequencies to consider for slow oscillation. It can be a wide exploratory range (2-30Hz) or a smaller specific range (theta: 4-8Hz)
Frequency for amplitude band (high): Frequencies to consider for high frequency bursts. It can be a wide exploratory range (40-250 Hz) or a smaller specific range (low gamma: 40-80Hz).
- Note: The lower bound of the frequency-for-amplitude band must be at least twice the upper bound of the frequency-for-phase band.
Frequency resolution for frequency for amplitude: This value indicates the bandwidth of the subbands in which the band
will be divided. Length of sliding time-window: Should at least include one full cycle of minimum
of interest (entered above). - Decreasing its length will improve temporal resolution, but decrease resolution for frequency for phase, and will also increase the processing time
If you have a event-related design for your study with more than 100 trials, and need very high temporal resolution (e.g. below 5 ms), you can use event-related PAC (ERPAC) algorithm (not implemented in brainstorm).
Number of signals to process at once: The file that is given to the function is processed in blocks of signals, and this option indicates how many signals will be processed simultaneously. Increasing this value can make the process a bit faster but requires a large amounts of memory.
Output options:
Save average PAC across trials: If this checkbox is selected and multiple files have been given to the process, then the process will save the average PAC across all trials in a new file.
Be careful with this option: although it is offered, averaging over frequencies (such as is required to do this) can be problematic and results from this option should be interpreted carefully.
File contents
The files saved by this process have the same structure as the time-frequency files, with an additional "sPAC" field. To review its contents, right-click on the PAC file > File > View file contents.
Note: Preferred phase of coupling extracted for each pair of -
is saved under sPAC.DynamicPhase. It can be exported to matlab and visualized there.
Comodulogram of time-resolved PAC:
To extract the comodulogram from the time-resolved PAC:
- Drag the time-resolved PAC file/files to the Process1 tab:
Note: if you have more than one time-resolved PAC file (e.g. from different trials), it is better to extract on comodulogram from all of them together to increase the SNR (considering the sparse nature of PAC in the three-dimensional time-
- plane
Select the process "Frequency > Time-resolved Phase-amplitude coupling > Extracting comodulogram from tPAC maps"
- Set the parameters as follows (these parameters are explained later in the text):
- The PAC map calculated from your simulated signal is saved in a new file. Double click on it.
This coupling map, shows the cross-frequency coupling in
- domain. The pseudo color of the map indicates the level of coupling for each pair of frequencies. Here in this example there is one mode arouond
Hz in the first 10 seconds, and two other modes apear in the second half of the signal on around 85 Hz, and 145 Hz, which reflect the oscillations that we initially used for synthesizing the data.
Process Options
Input options:
Time window: The time segment of the input file to be analyzed for extracting comodulogram.
Time-resolved Comodulogram:
- Input should be in second
- It is recommended to leave it as default: 0 -- extract one comodulogram from the total file
- Not recommended: If you would like to see the movie of time resolved comodulogram (in current set up of the brainstorm it can only be played in matlab and only once when it is computed), enter the window length to estimate several comodulograms each with that length from input tPAC map. -- You can zoom in/out in played movie in matlab.
Method:
- Method for estimating comodulogram. The method will be the same for both options if you have only one input (like this example with simulated data.
All sources/channels together: select this option if all sources belong to the same scout, or you do not expect different responses from them. (This will increase the SNR for comodulogram estimation, and make it more accurate and less noisy)
Each source/channel separately:
Time-resolved PAC in time - domain
...
Further analysis of PAC files
All functions for further analysis of PAC files (e.g., deriving the average of PAC maps, producing comodulograms, running a z-score transform with respect to baseline levels or surrogate data, etc. ) are available as processed located in Frequency > Time-resolved Phase-amplitude coupling.
Note: Basic functions for averaging or running statistics are not yet avalaible with Brainstorm [working on it...]. Yet, PAC output files are stored as Matlab files and contents can be readily processed with custom scripts of your own.
Yet, surrogate analysis is only available for single band tPAC maps, and can be accessed here: 
Further, tPAC traces can be z-scored with respect to surrogate data using the Process2 tab.
Practical suggestions for PAC analysis
How to best get started with PAC analysis? How do you decide about setting parameters? Here, we provide basic recommendations that we hope can be helpful:
General suggestions
- First, decide which frequency bands of interest for the frequency for phase (fP) and the frequency for amplitude (fA). These bands are usually hypothesized from the neurophysiological mechanisms you wish to study. If you have no previous clue from previous literature etc. you may want to start with visualizing the power spectrum density (PSD) of the empirical signals. It is recommended to define the fP band around a visible peak of the PSD [Aru et al. 2015], which indicates strong rhythmic signalling.
If you are not sure about the fA band, which can be elusive in PSD plots, you may first define it as a relatively wide frequency band. The resulting comodulogram output (i.e. a map with
and in horizontal and vertical axis, respectively, and a pseudo color illustrating the intensity of the coupling between each pair of frequencies) may indicate the restricted range of fA where PAC is more strongly expressed, for further single-band PAC analysis (see above). The intensity of PAC coupling is often stronger for slower fP oscillations. This is due in part to greater signal strengths at lower frequency in electrophysiological data (think 1/f power decrease). We therefore recommend that you subdivide a possibly wide frequency band of interest for
in to 2 or 3 narrower bands. This approach may also provide insight about possibly co-existing modes of PAC (i.e. having several [fP,fA] pairs coupled with PAC). It also prevents PAC modes at slower fP masking other modes. For example, if you want to investigate coupling between slow rhythms in from delta to beta bands (2-30 Hz) for and gamma oscillations for , we recommend you run PAC twice: first with delta to theta (2-8 Hz) for , then again with alpha to gamma (8-30 Hz).
Time-resolved analysis
For time-resolved analysis, if you are not sure about the exact frequencies that you want to focus on, one practical pipeline is to extract comodulogram with a relatively wide band (based on your initial hypothesis), and then find the main mode of coupling from the comodulogram. In the next step, focus on that particular mode of coupling to evaluate dynamical changes in coupling intensity (or preferred phase of coupling) with a single band analysis (single band option does not provide comodulogram but reflect a trace of dynamical changes in coupling parameters). With this strategy the temporal resolution will be optimized for the main mode of coupling, and also the computational load for intensive analysis (and checking significance of outcomes) will be much lower than looking at the whole interval.
Following steps provide more details:
- Start with comodulogram and find the main mode (or modes) of coupling based on comodulogram.
To decide about length of sliding time window in this step, you can start with two cycles of the slowest oscillation in the band of interest for
. For example, if you are interested in investigating potential PAC between theta to alpha oscillations (4-12 Hz) for and gamma rhythms for , a conservative short sliding time window can be set to two cycles of 4 Hz (i.e. 
. Certainly, if you have long inputs increasing the length of the window would improve the SNR.
The only point that you should keep in mind is that, tPAC do sparse estimation of PAC in a 3D space of fP, fA and time; therefore, you need enough number of points in time dimension to be able to retrieve a reasonable projection of it on fp - fA plain (called comodulogram). Hence, the maximum length of the window can be selected based on the input length. For example, in the case where your data is around 5 minutes long, if you set the window length to 20 s, considering the 50% overlap of windows you will end up with 5 * 60s / 10s = 30 points along time dimension (which is used in averaging for comodulogram extraction). It is not recommended to go below a certain number (e.g 30).
For extracting comodulogram you can also use none-time-resolved option (e.g. MVL) which use the whole time window and does not divide it into smaller chunks. It should be noted that in one experiment we observed higher resolution from comodulogram extracted using time-resolved method compared to a traditional option (see Fig. 5 in [Samiee & Baillet 2017]: https://www.sciencedirect.com/science/article/pii/S1053811917306195#fig5). The reason behind this improvement could be similar to the improvement in power spectral density estimation using Welch approach compared to original Bartlett's method.
- If you can find an obvious mode of coupling in the comodulogram, it would be better to investigate that mode of coupling using single band analysis for all regions of interest and time intevals (rather than maxPAC value extracted from comodulogram).
- If there is no obvious mode of coupling in averaged comodulogram, it is not recommended to go further in analysis. It would be better to take a more close look at modes availble in each single recording, and consistancy of preferred phase of coupling in time and in different trials/subjects.
If you observe a blub in comodulogram which looks like a vertical line, you need to further investigate your data to make sure this observed coupling is not a pseudo coupling coming from sharp event artifacts (e.g. spikes in LFP) or asymmetric oscillations (e.g. oscillations with saw tooth waveform shape). See [Cole & Voytek 2017; Aru et al. 2015] for more details.
A quick summary of the most important practical points in PAC analysis
1. Temporal resolution and minimum
Temporal resolution, and frequency resolution are complementary to each other. The minimum
that can be evaluated is equal 1/sliding time window, which means we need at least one cycle of minimum of interest in each sliding time window. For example if your band of interest is [4 - 8 Hz], minimum window length is 1/4 = 0.25 s.
2. Suggestion for deciding about initial band
If you do not have any priori hypothesis about which frequency band could be involved in the task as
, you can look at the power spectrum of the original signal. One of the requirements for having a PAC is a clear peak in that power spectrum around 
. For example if there is a peak around theta oscillations you can define your band of interest as 4 - 8 Hz.
3. Coupling strength for the main coupled pair
After detection of the best coupled
and , you can run tPAC with single band option only around that pair of interest. The window length can also get decreased to include at least two cycles of (May increase the temporal resolution compared to initial investigation because could be higher than the minimum in our band of interest). The output will be a time-series of PAC strength for that band of interest.