Plugins
This page contains information on the plugins currently available in Zeffiro Interface.
Electroencephalography (EEG) and magnetoencephalography (MEG) are two non-invasive neuroimaging techniques focused on localizing neural activity within the brain by measuring the electromagnetic field on the surface of the skull. To address the inverse problem associated with MEG and EEG, various approaches have been employed, including regularization, machine learning, and probabilistic methods.
Among probabilistic methods, one notable approach is SESAME (SEquential SemiAnalytic Monte Carlo Estimation) by Sommariva and Sorrentino, which is based on Sequential Monte Carlo (SMC) samplers. SMC samplers are iterative algorithms designed to approximate complex probability distributions. SESAME, as a specific SMC sampler, provides an approximation of the posterior distribution concerning the number, position, and intensity of active neural sources. This outcome enables the estimation of point estimates and the quantification of associated uncertainty. In particular, SESAME is based on the assumption that the positions of sources remain constant during the analysed time window. For more information, please refer to [1], [2], [3].
In 2022, it was implemented as an inverse tool for the Zeffiro Interface. In this tutorial, we will provide some examples of how to use SESAME within the Zeffiro interface: initially, two examples will be presented in which a single time instant will be analysed, then the same two examples will be proposed but analysing a time window, finally, we will proceed to visualize the dynamics of a reconstruction.
The first step is to open a new project.
Then in “example_projects” let us select “multicompartment_head_project”, which is a compartmental model of the head. At this point, all windows of Zeffiro will close automatically, don't worry: after a few seconds, they will reopen. Now we can visualize the model's surface: in “Mesh visualization tool” we must select "Visualize surfaces", then in “Figure tool” we will see the model’s surface.
Now we need to create a mesh: in “Mesh tool” we select “Create FEM mesh”; now a window like the one shown in the figure will open, displaying the time required to create the mesh. The required time varies based on your computer's processing power. If you do not have a particularly powerful computer and cannot afford to wait for an extended period, you can increase the mesh resolution. However, please keep in mind that this might lead to a loss of accuracy.
Now we can visualize the volume: like the case of the surface, in “Mesh visualization tool” we must select “Visualize volume”, then in “Figure tool” we will see the model’s volume. Now we can also see the mesh by pressing “Toggle edges” in “Figure tool”.
As can be seen in the figure, the mesh is of adaptive type; this is because we have selected “Refinement” and “Mash smoothing” in “Mesh tool” before creating a mesh.
Before proceeding with the forward problem, we need to select a script from the “Mesh tool”. At the top right of the mesh tool, we can select a script that best represents the problem we are interested in; in this case, we select the first one and then press “Run script”.
Now, a loading window will open, like when the mesh was created. Once the loading is complete, we can proceed with the forward problem.
Now we can proceed with the forward problem.
What we want to do is insert one or more electrical currents inside the model. To do this, we will not use the predefined synthetic sources on Zeffiro, but we'll load different ones by entering the following commands in the MATLAB command window:
load scripts/data/SEP_synth_source_data.mat
zef.synth_source_data = SEP_synth_source_data Now in Menu tool we select “Forward tools” and then “Find synthetic source”.
Now the window “Find synthetic source” will open.
On the left, we will find a series of synthetic sources. For this first example, we will select a single source, for instance the first one, Deep(14).
Be sure to change the selection from “All sources” to “Selected sources”; otherwise, all the sources listed on the left side of the window will be considered.
Before generating synthetic data, it is possible to visualize the position of the sources within the model by using the bottom “Plot source(s)”. Now the source can be displayed in the “Figure tool”. To view the source, especially for the deeper ones, it is necessary to increase the transparency of the surface (“Transp. Sur.”).
At this point, to generate synthetic data, it is necessary to press Generate time sequence and Create synthetic data in this order in the window “Find synthetic source”.
Now we are ready to plot the butterfly plot: in the Menu Tool we select Forwards tools and then Butterfly plot; at this point, window “butterfly_plot” will open. We need to make sure that the “Sampling frequency (Hz)” is the same in both windows “Find synthetic source” and window “butterfly_plot”.
Now it’s important to open another “Figure tool” window, otherwise Zeffiro will plot a two-dimensional graph on a three-dimensional one, consequently the butterfly plot won’t be displayed correctly. So, in “Menu tool” we select “Multi-tools” and then “figure tool”, we can also press Ctrl+5; so, a new “Figure tool” window will open, and now we can plot the butterfly plot by clicking “Plot” in window “butterfly_plot”.
Now a graph like the one in the figure should be displayed.
It’s possible that the “Time interval start” needs to be changed in window “butterfly_plot”: if in this case we had set the start time to 0, a graph like the one below would have been obtained. It’s clear, therefore, that no neural activity was detected in this time window. Thus, it’s necessary to change the “Time interval start”.
After generating the synthetic data, we can proceed to the inverse problem: reconstructing neural currents from measurements of the electromagnetic field present on the surface of the skull.
In “Menu tool” we select “Inverse tools” and then “SESAME” which is the last item we will find.
Now the window “SESAME App” will open.
In this window, we need to consider a few things:
- The first item we find is “Signal to noise ratio”; it refers to the signal-to-noise ratio estimation given to SESAME in decibels. Determining the ideal signal to noise ratio (SNR) for a specific dataset is not immediately evident due to the potential variety of noise sources.
- The second item is “Sample size”; the ideal sample size is around 100.
- “Sampling frequency (Hz)”; also in this case, this value should match what is found in the window “Find Synthetic Source”.
- “Low-cut frequency (Hz)”; the default setting for the low-cut frequency is 30Hz. However, it should be adjusted to a value that is appropriate for the specific data processing in the task question. The choice of 30Hz comes from a reference in a textbook suggesting the use of a 30-3000Hz bandwidth for localizing peaks of somatosensory evoked responses (Aminoff’s Electrodiagnosis in Clinical Neurology), available at trepo.tuni.fi/handle/10024/136311. However, setting the value to 0 for simulated data is the best choice.
- “High-cut frequency (Hz)”; it sets the upper bound for the applied signal.
- “Time interval start (s)”; it marks the beginning of the time interval from which the algorithm starts; it’s important to choose a start time when the signal isn’t just noise, otherwise, the algorithm won’t detect neural currents. To find the optimal start time, we can look at the butterfly plot when a field variation is recorded (simply click near the peak) and use the time indicated as “Time interval start”.
- “Time window (s)”; it is the length of the time interval that is considered as a single time point; the data is averaged over the time window before a reconstruction is found. The default length is 0, i.e., that there is no averaging.
- “Number of time steps”; it refers to the number of time points for which the brain activity is reconstructed (one reconstruction per time value).
- “Time step (s)”; it is the time interval that the inverse routine jumps forward between the reconstructions.
After these premises, we modify the algorithm parameters in a way that they are optimal, or nearly so, for our case.
At this point, we can press “Start” in the window “SESAME App” to initiate the algorithm.
Once the algorithm has finished, we can view the result on the model by clicking Plot source (s). In this case as well, we need to pay attention to which window “Figure tool” we are viewing the reconstructed sources. If you haven't interacted with any windows, the last selected one will be the butterfly plot window. Therefore, you need to select the one with the brain model – for example, by slightly rotating the model or performing a similar action – and then proceed to plot the reconstructed sources. If you're unable to view them, the advice is to close the brain model window and open a new one, following the same steps as you did for the butterfly plot (Ctrl+5). Then, press “Visualize surfaces” in the window “Mesh visualization tool” to view the model's surface, click Plot source (s) in the window “Find Synthetic Source” to view the previously created synthetic source, and finally plot the reconstructed source.
Now we can see the synthetic source in green and the reconstructed source in black. In this case, it can be observed that the performed reconstruction is a good reconstruction, and this is due to the input of the correct parameters in the window “SESAME App” before executing the SESAME algorithm.
To view the details of the reconstruction, we need to type the following command in the Command Window:
zef.SESAME Below is a brief description of some attributes that can be found:
- data are the sensor measurements
- sourspace are the head discretization coordinates
- pmap is the posterior probability map
- mod_sel is the posterior distribution on the number of dipoles
- estimated_dipoles are the source space grid points indices in which a source is estimated
- Q_estimated are the intensity of the sources estimated in the final iteration
- QV_estimated are the moments of the sources estimated in the final iteration
- dipole_positions are the coordinates of the estimated dipoles
We will now see an example in which two synthetic sources are present.
We need to repeat all the steps of the example with a single source until the window “Find synthetic source” is opened, alternatively, if you have previously tried the example with a single source or encountered other scenarios using the multicompartment head model, you can directly open the window “Find synthetic source”.
In the window “Find synthetic source” we will now select two synthetic sources by first choosing one and then, while holding down the Ctrl key, selecting the other. For example, let us choose the sources cortex(20) and deep(20).
It's important to always remember to switch from “All sources” to “Selected sources”. Then we need to click “Generate time sequence” and subsequently “Create synthetic data”. Just like in the previous case, to visualize the selected sources on the model, we need to click Plot source(s).
So, let us open a new window “Figure tool” to visualize the butterfly plot, always making sure that the “Sampling frequency (Hz)” is the same in both windows “Find synthetic source” and “butterfly_plot”.
In this case, we can see that the peak in the butterfly plot is around 0.125s.
Consequently, we will set the parameters for the SESAME algorithm as follow:
The reconstructed sources appear as those in the figure.
In this case as well, we can visualize the details of the reconstructed sources using the command
zef.SESAMEIf we try running several times the algorithm and plot the solution each time, we can observe that the estimated sources are different each time, especially for the deep source. To attempt to solve this problem, we can run the algorithm over a time window (in the previous example, we were looking at a single instant).
Let's consider a “Time window (s)” parameter of 0.002 in the window “SESAME App”, so that the SESAME algorithm takes measurements from 0.124 seconds when the peak begins to rise until 0.126 seconds when the peak starts to decrease.
Since “Sampling frequency (Hz)” is set to 5000, we will fix the value of “Time step (s)”, in window “SESAME App”, at 0.0002, which is 1/5000, to consider all the samples.
Once the algorithm has been run with the parameters just described, to visualize the estimated sources, it is necessary press “Plot source(s)”.
In this way, by running the algorithm several times, you can observe that the estimated sources do not differ too much from one run to another; however, it can still happen that the approximation for the deep source is not good. It is important to note that estimating deep sources is more challenging than estimating cortical sources.
For completeness, let's revisit the first example but with a nonzero time window.
Observing the butterfly plot, we can see that in this case as well, the peak rises and falls within 0.002 seconds, specifically it rises around 0.118 seconds and starts decreasing around 0.12 seconds.
So, in this case as well, we will set, in the window “SESAME App”, the value 0.002 for “Time window (s)” and 0.0002 for “Time step (s)”, for the same reason discussed earlier.
The visualization of the estimated sources is as follows:
If we look at the reconstruction obtained previously, we can notice that the latter is better than the first one, even though the first one was also good reconstruction.
An interesting analysis could be to visualize the dynamics of the dipoles.
To do this, let us go back to the window “Find synthetic source” and this time select “All sources”. Now, press “Generate time sequence” and then “Create synthetic data”. In this way, a synthetic data containing all the sources on the left side of the window “Find synthetic source” has been generated.
Now we can visualize the butterfly plot, always remembering to change in the window “Butterfly_plot” the “Samling frequency (Hz)” and set it to 5000, in other words, set it to the same value that is present in the window “Find synthetic source”.
After pressing “Plot” in the window “Butterfly_plot”, we can visualize the butterfly plot graph, and we can notice that there are multiple sources at different time points.
Now we can proceed with the reconstruction using the SESAME algorithm. So, in the window “SESAME App”, we set the parameters as shown in the figure.
Notice that all the parameters are the same as in the first two examples except for the “Number of time steps”, which is set to 16. This adjustment accounts for the fact that starting from 0.118 (“Time interval start (s)”) seconds and taking steps of 0.001 (“Time step (s)”) seconds, it allows us to include the last peak by the 16th step, Because the last step starts at 0.134 seconds, which is where the last peak is rising. Now in the window “Menu tool” let us select “Multi-tools” and than “Dynamical plot queue”. Now the window “Dynamical plot queue” will open, and at the bottom right, scroll until you find the “zef_plot_SESAME_dipoles” option. Right-click on this and then click “Add to script pipeline” and now we will see it in the top left corner. Under “Type” we need to click on “static” and change it to “dynamical”, as shown in the figure.
In the window “Mesh visualization tool”, change "Visualization type" to “Distribution (surface)” and set both “Transparency rec./surf.” to 0.9, otherwise, the estimated sources wouldn't be visible.
Now open a new window “Figure tool” (Ctrl+6), and then in the window “Mesh visualization tool” press “Visualize surface” to see the dynamics of the active sources.
Dynamic.dipole.plot.mp4
We can also visualize a single frame of the reconstruction by entering the frame number we are interested in under “Frame start, stop”, in the window “Mesh visualization tool”, and then by pressing “Visualize surface”. In the example in the figure, the twelfth frame is displayed.
To see the posterior, we can modify the “Plot scale” option in the window “Mesh visualization tool”, set it to “Logarithmic”, and set the “Plot threshold” to 1e-12.
By pressing “Visualize surfaces” in the window “Mesh visualization tool”, the result will be as follows
Dynamic.dipole.plot.log.mp4
To save the work you've done, simply select “Project” in the window “Menu tool” and then choose “Save as...”.
[1] Sorrentino A., Luria G., and Aramini R. "Bayesian multi-dipole modelling of a single topography in MEG by adaptive sequential Monte Carlo samplers." Inverse Problems 30.4 (2014): 045010.
[2] Sommariva S., and Sorrentino A. "Sequential Monte Carlo samplers for semi-linear inverse problems and application to magnetoencephalography." Inverse Problems 30.11 (2014): 114020.
[3] Viani A., Luria G., Bornfleth H., and Sorrentino A. "Where Bayes tweaks Gauss: Conditionally Gaussian priors for stable multi-dipole estimation." Inverse Problems and Imaging 15.5 (2021): 1099–1119.
