Welcome to this introductory tutorial for the ConPy package. Together with MNE-Python, we can use it to perform functional connectivity estimation of MEG data. This tutorial was written to be used as an exercise during my lectures. In lieu of my lecture, you can read this paper to get the theoretical background you need to understand the concepts we will be dealing with in this exercise.
Ok, let's get started! I have similated some data for you. It's already stored on the virtual server you are talking to right now. In this simulation, a placed a couple of dipole sources on the cortex, sending out a signal in a narrow frequency band. During the first part of the recording, these sources are incoherent with each other. During the second part, some of the sources become coherent with each other. Your task is to find out:
We will use MNE-Python and ConPy to aswer the above questions.
In [ ]:
# Don't worry about warnings in this exercise, as they can be distracting.
import warnings
warnings.simplefilter('ignore')
# Import the required Python modules
import mne
import conpy
import surfer
# Import and configure the 3D graphics backend
from mayavi import mlab
mlab.init_notebook('png')
# Tell MNE-Python to be quiet. The normal barrage of information will only distract us. Only display errors.
mne.set_log_level('ERROR')
# Configure the plotting interface to display plots in their own windows.
# The % sign makes it a "magic" command: a command ment for the notebook environment, rather than a command for Python.
%matplotlib notebook
# Tell MNE-Python and PySurfer where to find the brain model
import os
os.environ['SUBJECTS_DIR'] = 'data/subjects'
# Let's test plotting a brain (this is handled by the PySurfer package)
surfer.Brain('sample', hemi='both', surf='pial')
If the code in the above cell ran without any errors, we're good to go. Let's load the stimulation data. It is stored as an MNE-Python Epochs file. To load it, you must use the mne.read_epochs function. To see how it works, you need to take a look at the documentation for this function. You can call up the documentation of any function by appending a ? to the function name, like so:
In [ ]:
mne.read_epochs?
The documentation shows us that mne.read_epochs takes one required parameter (fname) and three optional parameters (proj, preload and verbose). You can recognize optional parameters by the fact that they have a default value assigned to them. In this exercise, you can always leave the optional parameters as they are, unless explicitly instructed to change them.
So, the only parameter of interest right now is the fname parameter, which must be a string containing the path and filename of the simulated data, namely: 'data/simulated-data-epo.fif'.
Go ahead and call the mne.read_epochs function to load the stimulated data. Store the result in a variable called epochs:
In [ ]:
# Write your Python code here
In [ ]:
# If your code in the above cell is correct, executing this cell print some information about the data
print(epochs)
"Epochs" are snippets of MEG sensor data. In this simulation, all sensors are gradiometers. There are two epochs, appoximately 10 second in length: one epoch corresponding to the (simulated) subject "at rest" and one epoch corresponding to the subject performing some task.
Most objects we'll be working with today have a plot method. For example, the cell below will plot the epochs object:
In [ ]:
# The semicolon at the end prevents the image from being included in this notebook
epochs.plot();
In the epochs plot, you can use the scrolling function of your mouse/trackpad to browse through the channels. The vertical dashed line indicates where one epoch ends and the next one begins.
To find out, let's plot the power spectal density (PSD) of the signal. The PSD is computed by applying a Fourier transform to the data of each MEG sensor. We can use the plot_psd method of the epochs object to show it to us. By default, it will show us the average PSD across the sensors, which is good enough for our purposes. Check the documentation of the epochs.plot_psd method to see what parameters are required (remember: you are free to ignore the optional parameters).
In [ ]:
# Write here the code to plot the PSD of the MEG signal
If you were to name one frequency at which the sources are sending out a signal, what would that frequency be? Fill in the answer below. We'll use it in the upcoming tasks:
In [ ]:
# Fill in the source frequency, in Hertz
source_frequency = ###
Ok, so now we know the frequency at which to look for sources. How many dipole sources did I use in the simulation? To find out, we must look at which sensors have the most activity at the frequency of the sources. The plot_psd_topomap method of the epochs object can do that for us. If you call it with the default parameters, it will plot so called "topomaps" for the following frequency bands:
| Name | Frequency band |
|---|---|
| Delta | 0-4 Hz |
| Theta | 4-8 Hz |
| Alpha | 8-12 Hz |
| Beta | 12-30 Hz |
| Gamma | 30-45 Hz |
Try it now: take a look at the documentation for the epochs.plot_psd_topomap method and plot some topomaps:
In [ ]:
# Write here the code to plot some PSD topomaps
Take a look at the topomap corresponding to the frequency band that contains the frequency at which the sources are sending out their signal. How many sources do you think I simulated? Fill in your answer below:
In [ ]:
number_of_sources = ###
Looking at the topomaps will give you a rough location of the sources, but let's be more exact. We will now use a DICS beamformer to localize the sources on the cortex.
To construct a DICS beamformer, we must first estimate the cross-spectral density (CSD) between all sensors. You can use the mne.time_frequency.csd_morlet function to do so. Go check its documentation.
You will find that one of the parameters is a list of frequencies at which to compute the CSD. Use a list containing a single frequency: the answer to Question 1 that you stored earlier in the source_frequency variable. In Python code, the list can be written like this: [source_frequency]. Store the result of mne.time_frequency.csd_morlet in a variable called csd.
In [ ]:
# Write here the code to construct a CSD matrix
In [ ]:
# If the code in the cell above is correct, executing this cell will plot the CSD matrix
csd.plot()[0]
If you examine the CSD matrix closely, you can already spot which sources are coherent with each other. Sssshhh! we'll look at it in more detail later. For now, let's compute the DICS beamformer!
The next functions to call are mne.beamformer.make_dics and mne.beamformer.apply_dics_csd. Lets examine them more closely.
mne.beamformer.make_dicsThis function will create the DICS beamformer weights. These weights are spatial filters: each filter will only pass activity for one specific location on the cortex, at one specific frequency(-band). In order to do this, we'll need a leadfield: a model that simulates how signals on the cortex manifest as magnetic fields as measured by the sensors. MNE-Python calls them "forward solutions". Luckily we have one lying around: the 'data/simulated-data-fwd.fif' file contains one. You can load it with the mne.read_forward_solution function. Take a look at the documentation for that function and load the forward solution in the variable fwd:
In [ ]:
# Write your code to read the forward solution here
In [ ]:
# If the code in the above cell is correct, executing this cell will plot the source grid
fwd['src'].plot(trans='data/simulated-data-trans.fif')
For this exercise, we use a very sparse source grid (the yellow dots in the plot). This grid is enough for our purposes and our computations will run quickly. For real studies, I recommend a much denser grid.
Another thing you'll need for the DICS beamformer is an Info object. This object contains information about the location of the MEG sensors and so forth. The epochs object provides one as epochs.info. Try running print(epochs.info) to check it out.
Now you should have everything you need to create the DICS beamformer weights using the mne.beamformer.make_dics function. Store the result in the variable filters:
In [ ]:
# Write your code to compute the DICS filters here
In [ ]:
# If the code in the above cell is correct, executing this cell will print some information about the filters
print('Filters have been computed for %d points on the cortex at %d frequency.' %
(filters['weights'].shape[1], filters['weights'].shape[0]))
print('At each point, there are %d source dipoles (XYZ)' % filters['n_orient'])
mne.beamformer.apply_dics_csdWith the DICS filters computed, making a cortical power map is straightforward. The mne.beamformer.apply_dics_csd will do it for you. The only new thing here is that this function will return two things (up to now, all our functions only returned one thing!). Don't panick. The Python syntax for dealing with it is like this:
power_map, frequencies = mne.beamformer.apply_dics_csd(...)
See? It returns both the power_map that we'll visualize in a minute, and a list of frequencies for which the power map is defined. Go read the documentation for mne.beamformer.apply_dics_csd and make the powermap:
In [ ]:
# Write your code to compute the power map here
In [ ]:
# If the code in the above cell is correct, executing the cell will plot the power map
power_map.plot(hemi='both', smoothing_steps=20);
Use the mouse/trackpad to rotate the brain around. Can you find the sources on the cortex? Even though I've simulated them as dipole sources, they show more as "blobs" in the power map. This is called spatial leaking and is due to various inaccuracies and limitations of the DICS beamformer filters.
The simulated recording consists of two parts (=epochs): during the first epoch, our simulated subject is at rest and the sources are not coherent. During the second epoch, our simulated subject is performing a task that causes some of the sources to become coherent. It's finally time for some connectivity estimation!
We'll first tackle one-to-all connectivity, as it is much easier to visualize and the results are less messy. Afterward, we'll move on to all-to-all connectivity.
For this, we must first define a "seed region": one of the source points for which we will estimate coherence with all other source points. A common choice is to use the power map to find the source point with the most power. To find this point, you can use the .argmax() method of the power_map.data object. This is a method that all data arrays have. It will return the index of the maximum element in the array, which in the case of our power_map.data array will be the source point with the maximum power.
Go find your seed point and store it in the variable seed_point:
In [ ]:
# Write your code to find the seed point here
In [ ]:
# If the code in the above cell is correct, executing this cell will plot the seed point on the power map
brain = power_map.plot(hemi='both', smoothing_steps=20) # Plot power map
# We need to find out on which hemisphere the seed point lies
lh_verts, rh_verts = power_map.vertices
if seed_point < len(lh_verts):
# Seed point is on the left hemisphere
brain.add_foci(lh_verts[seed_point], coords_as_verts=True, hemi='lh')
else:
# Seed point is on the right hemisphere
brain.add_foci(rh_verts[seed_point - len(lh_verts)], coords_as_verts=True, hemi='rh')
You may need to rotate the brain around to find the seed point. It should be drawn as a white sphere.
Up to now, we've been using all data. However, we know our sources are only coherent during the second part. Executing the cell below will split the data into a "rest" and "task" part.
In [ ]:
# Splitting the data is not hard to do.
epochs_rest = epochs['rest']
epochs_task = epochs['task']
To estimate connectivity for just the epochs_task part, we need to compute the CSD matrix on only this data. You've computed a CSD matrix before, so rince and repeat: compute the CSD on just the epochs_task data and store it in the csd_task variable:
In [ ]:
# Write your code here to compute the CSD on the epochs_task data
In [ ]:
# If the code in the above cell is correct, executing this cell will plot the CSD matrix
csd_task.plot()[0]
Now you are ready to compute one-to-all connectivity using DICS. It will take two lines of Python code. First, you'll need to use the conpy.one_to_all_connectivity_pairs function to compute the list of connectivity pairs. Then, you can use the conpy.dics_connectivity function to perform the connectivity estimation. Check the documentation for both functions (remember: you can leave all optional parameters as they are) and store your connectivity result in the con_task variable:
In [ ]:
# Write your code here to compute one-to-all connectivity for the "task" data
In [ ]:
# If the code in the above cell is correct, executing this cell will print some information about the connectivity
print(con_task)
To visualize the connectivity result, we can create a cortical map, where the value at each source point is the coherence between the source point and the seed region. The con_task object defines a .make_stc() method that will do just that. Take a look at its documentation and store the map in the coherence_task variable:
In [ ]:
# Write your code here to compute the coherence map for the epochs_task data
In [ ]:
# If the code in the above cell is correct, executing this cell will plot the coherence map
brain = coherence_task.plot(hemi='both', smoothing_steps=20)
lh_verts, rh_verts = coherence_task.vertices
if seed_point < len(lh_verts):
# Seed point is on the left hemisphere
brain.add_foci(lh_verts[seed_point], coords_as_verts=True, hemi='lh')
else:
# Seed point is on the right hemisphere
brain.add_foci(rh_verts[seed_point - len(lh_verts)], coords_as_verts=True, hemi='rh')
Which source points seem to be in coherence with the seed point? Double-click on the text-cell below to edit it and write down your answer.
Double-click here to edit this text cell. Pressing CTRL+Enter will transform it back into formatted text.
If you examine the coherence map, you'll find that the regions surrounding the seed point are coherent with the seed point. This is not because there are active coherent sources there, but because of the spatial leakage. You will always find this coherence. Make a one-to-all coherence map like you did above, but this time for the epochs_rest data (in which none of the sources are coherent). Store the connectivity in the con_rest variable and the coherence map in the coherence_rest variable:
In [ ]:
# Write your code here to compute connectivity for the epochs_rest data and make a coherence map
In [ ]:
# If the code in the above cell is correct, executing this cell will plot the coherence map
brain = coherence_rest.plot(hemi='both', smoothing_steps=20)
lh_verts, rh_verts = coherence_rest.vertices
if seed_point < len(lh_verts):
# Seed point is on the left hemisphere
brain.add_foci(lh_verts[seed_point], coords_as_verts=True, hemi='lh')
else:
# Seed point is on the right hemisphere
brain.add_foci(rh_verts[seed_point - len(lh_verts)], coords_as_verts=True, hemi='rh')
See? You'll find that also when no coherent sources are active, there is an area of coherence surrounding the seed region. This will be a major problem when attempting to estimate all-to-all connectivity.
One way to deal with the spatial leakage problem is to make a contrast between the "task" and "rest" segments. Since the coherence due to spatial leakage is the same for both segments, it should cancel out.
Connectivity objects, like con_task and con_rest have support for common math operators like +, -, * and /. Creating a constract between the two object is therefore as simple as con_task - con_rest. Im the cell below, make a new coherence map of the contrast and store it in the coherence_contrast variable:
In [ ]:
# Write your code here to compute a contrast between the "task" and "rest" connectivity and make a coherence map
In [ ]:
# If the code in the above cell is correct, executing this cell will plot the coherence map
brain = coherence_contrast.plot(hemi='both', smoothing_steps=20)
lh_verts, rh_verts = coherence_contrast.vertices
if seed_point < len(lh_verts):
# Seed point is on the left hemisphere
brain.add_foci(lh_verts[seed_point], coords_as_verts=True, hemi='lh')
else:
# Seed point is on the right hemisphere
brain.add_foci(rh_verts[seed_point - len(lh_verts)], coords_as_verts=True, hemi='rh')
If all went well, you'll see that the coherence due to spatial leakage has disappeared from the coherence map.
Use the conpy.all_to_all_connectivity_pairs function to compute the connectivity pairs in an all-to-all manner. Then, use the conpy.dics_connectivity function like before to create a connectivity object for the "task" and a connectivity object for the "rest" data segments. Store them in the all_to_all_task and all_to_all_rest variables. You'll notice that computing all-to-all connectivity takes a while... Finally, create the contrast between them and store it in the all_to_all_contrast variable.
In [ ]:
# Write your code to produce all-to-all connectivity estimates for the "rest" and "task" segments
# and the contrast between them.
In [ ]:
# If the code in the above cell is correct, executing this cell will print some information about the connectivity
print(all_to_all_contrast)
How to visualize this all-to-all connectivity? This is a question worth pondering a bit. But for this exercise, we can get away with producing a coherence map like we did with the one-to-all connectivity. The value of the coherence map is, for each source point, the sum of the coherence of all connections from and to the source point.
Executing the cell below will plot this coherence map. Can you spot the connectivity between the sources?
In [ ]:
# This cell will plot the coherence map
all_to_all_coherence = all_to_all_contrast.make_stc()
all_to_all_coherence.plot(hemi='both', smoothing_steps=20);