LeadField matrix

đź’¬ Support & Discussions
,
1

Hi everyone, I am trying to obtain a leadfield matrix (channels x brain regions) which in my case would be 62x200 (I am using a brain atlas with 200 parcels). Any idea on how to move from the “canonical” leadfield matrix (channels x sources) to the one I’ve just described. Since mne has a function to extract the parcels time series, there should be an easy way to extract the leadfield matrix without scripting it from scratch.

Thanks in advance

2

what would be for you the forward field of a parcel? putting the full parcel with a constant activation of 1?

Alex

3

Hi @agramfort thanks for reply to this. It would be the same as the “canonical” leadfield array but instead having the vertices on the second dim, it will have the average of the values for the vertices that falls in the parcels. Anyway, I’ve just solved the issue myself…just some extra-work with mne.decimate

Thanks

4

ok that’s equivalent to what I suggest above. It’s as if you were to put a constant activity over each parcel.

you would just want the leadfield as a numpy array?

Alex

5

@agramfort please find the piece of code I wrote below for a sanity check and for other people who might encounter the same issue. Could you please take a look at it?

stcs = mne.read_source_estimate(op.join(eeg_path + 'sub_0014-stc'))

lh_regmap = nib.freesurfer.io.read_annot(op.join(subjects_dir + \
                                    subject +  '/label/lh.Schaefer2018_200Parcels_7Networks_order.annot'))[0]

rh_regmap = nib.freesurfer.io.read_annot(op.join(subjects_dir + \
                                    subject +  '/label/rh.Schaefer2018_200Parcels_7Networks_order.annot'))[0]

regmap = np.concatenate([lh_regmap, rh_regmap+np.max(lh_regmap)])

vtx2use = np.concatenate([stcs.vertices[0], stcs.vertices[1]+np.max(stcs.vertices[0])])
regions_idx = regmap[vtx2use]


fwd = mne.read_forward_solution(op.join(eeg_path + 'sub_0014-fwd'))
leadfield = fwd['sol']['data']
print("Leadfield size : %d sensors x %d dipoles" % leadfield.shape)



fwd_fixed = mne.convert_forward_solution(fwd, surf_ori=True, force_fixed=True,
                                         use_cps=True)
leadfield = fwd_fixed['sol']['data']
print("Leadfield size : %d sensors x %d dipoles" % leadfield.shape)



new_leadfield = []
for xx in range(1,np.unique(regions_idx).shape[0]):
    new_leadfield.append(np.mean(np.squeeze(leadfield[:,np.where(regions_idx==xx)]),axis=1))
    
new_leadfield = np.array(new_leadfield).T

print("Leadfield size : %d sensors x %d regions" % new_leadfield.shape)```


Output is: 

Leadfield size : 62 sensors x 23121 dipoles
Leadfield size : 62 sensors x 7707 dipoles
Leadfield size : 62 sensors x 200 regions


Could you please confirm that is correct? 

Thanks
6

it’s hard to tell without running the code. Can you share snippet using the mne-sample-data so it works also on my machine?

thanks
Alex

7

@agramfort I apologize but I don’t have the mne-sample-data here. Anyway it should be easy to run. Please set up:
stcs = the source estimate
lh_regmap = path to the annotation file of the .annot file (in your case it would be the lh.aparc.annot)
rh_regmap = same as lh_regmap but for the right hemisphere.

Here is the code.

stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, inv_method, pick_ori=None, return_generator=True)
lh_regmap = nib.freesurfer.io.read_annot(op.join(subjects_dir + \
                                    subject +  '/label/lh.aparc.annot'))[0]

rh_regmap = nib.freesurfer.io.read_annot(op.join(subjects_dir + \
                                    subject +  '/label/rh.aparc.annot'))[0]

Hope this helps

8

hi Davide,

yes it’s easy for me to do but it’s still minutes for me. You can easily do this yourself too.
Let’s balance the efforts here so we can collectively progress.

Alex

1 Like
9

You can download the data by running mne.datasets.sample.data_path()

10

Greetings,

Apologies for the delay, @richard and @agramfort . Kindly find the code intended for verifying accuracy at the following link:

Feel free to execute it directly on Colab.

Your feedback is greatly appreciated.

Please let me know

Thanks

11

ok I get now what you do here. You average columns of the forward within a parcel.
Note that this is not what we consider MNE standard practice to get one time course per label.
We rather use so called “label time course” see eg mne.extract_label_time_course — MNE 1.4.2 documentation2023-08-12T16:28:59Z and mne.SourceEstimate — MNE 1.4.2 documentation

These typically do not do plan averages in the parcel. One key reason is that sources have a directions so by averaging the raw time courses obtained with fixed orientation or pick normal you have cancellation effects. You will have the same problem with your approach as you have columns of the forward solution and these can be basically of opposite polarity for dipoles eg located on opposite walls of a gyrus.

Alex