I want to use eLORETA to estimate sources in a lesioned brain (stroke patients). The pial surface reconstruction with FS looks reasonable, showing a chunk of the cortex missing. As the forward model tries to fit sources in the cortex by placing sources between the pial and white matter surface, we end up having sources placed on the surface of the lesioned area. However, we know that the surface of the lesion is not cortical. So, I want to modify the source model by removing any source within some threshold distance of the wall of the lesion. I'll do that step in Python. So, the process would look like that:
1 - Import an FS subject
2 - Export the source model from BS
3 - Process the source model as described and save this modified version
4 - Import the modified version of the source model in BS to replace the original
5 - Run OpenMEEG
6 - Continue analyses as usual
First, does this pipeline make sense with BS pipeline? Second, what are the files in the BS file system database that corresponds to the head model that I would need to modify? I first thought about modifying tess_cortex_pial_low.mat but that does not seem right as this is the pial surface, not the source model... unless sources are placed systematically at the vertices of this surface in which case removing the corresponding vertices could make sense?
In Brainstorm, when using a surface source model, all the vertices of the cortex surface are used as possible sources of activity. If you want to exclude some brain regions from the minimum norm source estimation, removing the vertices from the cortex surface (tess_cortex_pial_low.mat) is the easiest approach. You need to do this before computing the forward model.
To remove vertices from the cortex surface, the simplest approach is to create a scout with the vertices to remove (either from the interface or programmatically), and then from the Scout tab, use the menu Scout > Edit surface > Remove selected scouts.
An alternative is to use a volume source model: https://neuroimage.usc.edu/brainstorm/Tutorials/TutVolSource
In this case, the source space is saved in the GridLoc field of the forward model, and there is not interactive menu to removed grid points from it. It is however possible to edit the GridLoc field from the headmodel file and simply remove some points from it, before computing the MNE solution.
I'm following up on the discussion for source reconstruction of EEG data in stroke subjects. Following your recommendation, Christian and I created a modified surface and removed lesion scout to eliminate potential sources of neural activities originating from the peri-lesion surface. Apparently, this seems to be working fine for removing the lesion; however, when used for creating head models, the openMEG did not show any differences between reconstructed source fields before and after lesion removal. We think this might be due to 2 potential reasons: 1) openMEG still continues to use the intact brain surface for modeling, or 2) it somehow automatically recreates vertices for damaged brain areas and does the modeling similar to the original surface reconstruction. We would appreciate any input you may have to help us do source modeling correctly on the stroke data. I've uploaded two sample figures for the lead fields calculated using openMEG before and after lesion surface removal.
What you show is actually quite reassuring: you obtain similar source maps if you change a bit the source space.
The surface that is used to display the MNE results when you double-click on them is really the source space that is used for estimating the sources. You should not fear that Brainstorm uses some cached version of the cortex surface.
If you look at inverse results computed from these forward models: unless all the signal observed in EEG is identified as originated from the area you removed from the cortex surface, you should not expect major differences in the two versions, especially if you look at the flattened unconstrained source estimation in a 3D figure.
However, the results cannot be strictly equal. If you want to observe this: create a scout with one vertex only in the two versions of the cortex (the same one), set the scout function to "All", extract the scout time series for both source files, concatenate the two files and display the results, or plot the signals from Matlab: they should have the same shape but if you zoom in, they should not overlap perfectly.
Hi @Francois. I am not convinced that the approach you proposed is working as expected. Computing the largest relative difference between the lead field gains obtained with the original mesh (x) and the mesh missing a large lesion (\hat{x}) only for the vertices present in both mesh using the formula
I get a maximal error of 5e-12, which seems to me to be only numerical accuracy errors. A similar approach in MNE-Python gives some hot spots with relative differences up to about 1.5. Although, I would expect some stability in the estimation algorithm, I would also expect that a big chunk of cortex missing would change the gains in the rest of the cortex. Any advices?
I get a maximal error of 5e-12, which seems to me to be only numerical accuracy errors. A similar approach in MNE-Python gives some hot spots with relative differences up to about 1.5.
Does this refer to the Gain matrix or the leadfield sensitivity, as shown in the screen captures above?
Can it be that some of the leadfield post-processing (depth weighting, dipoles orientations, whitening) is done at the time of the inverse computation in Brainstorm, while it is done before your tests in MNE-Python?
I would expect some stability in the estimation algorithm, I would also expect that a big chunk of cortex missing would change the gains in the rest of the cortex.
The Gain matrix indicates how each source dipole independently projects on the sensors. I'm not sure I understand why you would expect it to be different for a dipole A depending on whether there is an adjacent dipole B or not.
In my understanding, the gains in the Brainstorm "headmodel" file should be identical (minus numerical approximations) for a dipole that is present in the two files, even if most of the rest of the cortex is missing in one of the two.
The Gain matrix indicates how each source dipole independently projects on the sensors. I'm not sure I understand why you would expect it to be different for a dipole A depending on whether there is an adjacent dipole B or not.
Hmm... right... I think you are right with respect to the lead field. Sorry, I did not take enough time to think this through and I had more in mind the inverse modeling (i.e., how the sources for a single dipole will be different if you add another dipole because the second dipole will now be responsible for some of the EEG activity that was attributed to the first dipole originally... but you are right, the lead field is independent of all that). This kind of validation of the effect of a lesion would need to be more on sources distribution than the lead field, I suppose. Arguably, though, a big chunk of missing cortex would also impact the head model and, hence, the lead field, but it is probably something that would be better captured by a FEM than a BEM model (much as the effect of ventricles are)... I'll do some further checks. Thanks @Francois for the input.
Jumping in the discussion a bit late, sorry: a lesion would alter the head model/gain vectors of sources if its conductivity parameters are different from surrounding tissues ; for instance if this is a resected portion of cortex, then its conductivity is that of CSF. A FEM/BEM model can explicitly take these aspects into account in the head model, again, if the region in question is properly delineated. See this recent paper: https://pubmed.ncbi.nlm.nih.gov/35781078/
From the scout tab, you can load a MRI binary mask or a labelled volume parcellation as surface or volume scouts using one of the file formats "Volume mask or atlas" in the scout import file selection window.
These are research considerations for new possible tools to be developed in the future for MEG/EEG forward modeling. At the present time, I don't have any solution to suggest for handling this with the current OpenMEEG BEM workflow in Brainstorm.
It would be easier to think about with the DUNEuro FEM approach, but I wouldn't know how to add easily the lesion masks in the tetrahedral FEM mesh - hopefully we'll get there at some point.
Something I would change is the way brainstorm deals with tissue names when imported from SimNIBS4. In in_tess_simnibs.m, line 60 and below, the function assigns some default names to the tissues, but this won't work with exceptions such as including an extra tissue, or for instance having the lesion instead of "electrodes".
case 'SIMNIBS4'
% Relabel the electrodes and gel
MeshMat.Tissue(MeshMat.Tissue==100) = 11;
MeshMat.Tissue(MeshMat.Tissue==500) = 12;
% Default tissue labels (from file: final_tissues_LUT.txt)
MeshMat.TissueLabels = {'white', 'gray', 'csf', 'skull', 'scalp', 'eyes', 'compact', 'spongy', 'blood', 'muscle', 'electrode', 'gel'}; % I WOULD CHANGE THIS
this won't work with exceptions such as including an extra tissue, or for instance having the lesion instead of "electrodes".
I made this choice of using a hard-coded list after getting the confirmation from the SimNIBS developers that the list of tissues was always the same in the CHARM output.
How would you obtain a list of tissues different from this one?
Note that this code is dedicated specifically to the import of SimNIBS4/CHARM outputs.
Different pipelines should probably be handled in a different use case.
If this exchange on the SimNIBS GitHub issue is not covering all the possible use cases of SimNIBS4/CHARM, please let us know, either here or on GitHub.
A solution would be to always read the file final_tissues_LUT.txt instead of relying on a hardcoded list.
But I am reluctant to do this if we don't have a clear reason (= a use case) for doing it, as it requires the presence of an additional file (and development time).
As it is now, it works if the user only needs the "standard" tissues. In case of additional tissues (e.g. lesion), and, for example, absence of electrode or gel, 1) the new tissue should not be labeled >12 or the pipeline will crash and 2) if labeled between 10 and 12, it will be assigned as electrode or gel whatsoever.
Before checking the code in brainstorm, I hoped that it would rely on the final_tissues_LUT.txt. That would make the pipeline more flexible.
In the meantime, as a temporary solution I just changed the code myself by editing the label names in the brainstorm function.
Can you please describe how you obtain a different list of tissues than the one used in the code?
Is this coming from some personal scripts, or standardized and publicly shared processing pipelines (such as CHARM)?
