Home
Zeffiro Interface (ZI), © 2018- Sampsa Pursiainen & ZI Development Team, is an open source code package constituting an accessible tool for multidisciplinary finite element method (FEM) assisted forward and inverse simulations in complex geometries. ZI has been started as a project to enable the use of finite elements (FE) in Electro/Magnetoencephalography (EEG/MEG). Its parameter structure is modular to allow a wider applicability in multiphysics applications of inverse problems. Install ZI using zeffiro_downloader to allow automatic updates between the local and remote repositories. To obtain the downloader in matlab command line, use:
urlwrite('https://raw.githubusercontent.com/sampsapursiainen/zeffiro_interface/master/zeffiro_downloader.m','zeffiro_downloader.m');
The the repository can be downloaded to folder 'my_folder' by the command:
zeffiro_downloader('folder_name','my_folder');
With ZI, one can segment a realistic multilayer geometry of a human brain and generate a multi-compartment FE mesh, if triangular surface grids (for example in .STL, .DAT, or .ASC file format) are available. A suitable surface segmentation can be produced, for example, with the FreeSurfer software suite (Copyright © FreeSurfer, 2013). A folder containing multiple meshes FreeSurfer's .ASC format can be built as a single segmentation via import utility. An example of such a segmentation with an import file (.ZEF) has been included in the script folder. Also a parcellation created with FreeSurfer can be imported to enable distinguishing different brain regions and, thereby, analysing the connectivity of the brain function over a time series. Multiple compartments can be defined as active, allowing the analysis of the sub-cortical strucures. In each compartment, the orientation of the activity can be either normally constrained or unconstrained. The main routines of ZI can be accelerated significantly in a computer equipped with a graphics computing unit (GPU). Processes that will benefit from GPU parallelization include, in particular, the forward simulation phase which can involve high-resolution FE mesh and lead field matrix generation. Alternatively, these processes can be parallelized in a multicore CPU to obtain an enhanced performance if GPU is not available. In ZI, the forward simulation routines tie together with advanced inversion and optimization routines. The reconstructed distributions can be presented on the 3D FE mesh, e.g., as a time-lapse.
The ZI code utilizes the Matlab platform (The MathWorks, Inc.). A recent Matlab version (R2020a upwards) is recommendable. The name Zeffiro originates from an Italian high-end road bicycle. It is also Italian for 'a gentle breeze' referring to the aimed user-friendliness of ZI. ZI current has a plugin functionality which allows one to add more components easily. A brief introduction to the essential features of the interface can be found below. Quick video examples on how to use ZI can be found in YouTube.
In case any questions arise or there is an interest to develop the package, e.g., to write or test an inversion routine, please feel free to ask anything, e.g., by email sampsa.pursiainen(at)iki.fi or sampsa.pursiainen(at)tuni.fi. ZI is not designed to be used in clinical applications. The authors do not take the responsibility of the results obtained with ZI using clinical data.
Figure 1: A screenshot of the interface with figure, mesh, parcellation, and option tools opened.
A scientific article of ZI code itself has been published in Neuroinformatics. The essential mathematical techniques used in the interface have been reviewed and validated in (Miinalainen et al., 2019 and Pursiainen, 2012). The present FE approach is based on a current preserving H(div) source model in which the primary current distribution of the neural activity is modeled via vector valued and divergence conforming FE basis functions. The most well-known example of those is the linear Whitney (Raviart-Thomas) function. The earliest applications of these to EEG/MEG can be found in the dissertation works Tanzer, 2005 and Pursiainen, 2008. The divergence conforming basis functions can be interpreted as dipolar sources (Pursiainen et al., 2011) which form a piecewise continuous vector field in the brain. In comparison to the classical (monopolar) FE source modeling techniques, the current preserving approach provides a comparable and even superior performance for dipole approximation. The present technique combines linear (face-intersecting) and quadratic (edgewise) H(div) basis functions via the Position Based Optimization (PBO) method (Bauer et al., 2015) and the 10-source stencil in which 4 face and 6 edge functions are applied for each tetrahedron containing a source (Pursiainen et al., 2016). Alternatively, the linear Whitney functions, i.e., a 4-source stencil (4 face functions) can be applied for lower memory consumption. The hierarchical Bayesian reconstruction approach utilized in the inversion was introduced in (Calvetti et al., 2009) and applied for a realistic head model, e.g., in (Lucka et al., 2012).
ZI creates a single data structure (struct) zef in Matlab's base workspace. All the parameters and variables, such as the lead field matrix, measurement data and reconstruction, can be accessed via the zef structure. It (current project) can be also saved or loaded on a hard disk. This way, the lead field and other variables will be easily accessible for the user. The export functions enable saving separate variables, for example, the lead field matrix which, among other things, enables merging EEG and MEG lead fields for parallel EEG/MEG. The main window contains the essential functions needed for mesh and lead field generation and visualizing the results. Additional adjustments can be made via the Miscellaneous options dialog box (Edit menu) which allows controlling, e.g., smoothing, meshing and plotting parameters. The items of the Inverse tools menu enable inversion of the data. For external inverse procedure development, the list of most important variables is the following:
- the lead field matrix
zef.L, - the measurement data
zef.measurements(a matrix or a cell array with the number of rows and equal to that ofzef.Land the time steps in the dataset, respectively), - source locations
zef.source_positions(source positions corresponding to the columns ofzef.Lin the respective order), - source directions
zef.source_directions(an empty array, if Cartesian orientations are used, meaning that the source position/orientation pair for the columns ofzef.Lis given by the following pattern: position 1, xyz; position 2, xyz; position 3, xyz...), - the reconstruction
zef.reconstruction(candidate solution forzef.L*zef.reconstruction=zef.measurements).
The reconstruction can be visualized using the interface. For creating a new inversion tools menu item, one needs a .fig file and two .m files for initializing and updating the parameters. New items will be welcome and appreciated.
The current version also allows using a graphics card (GPU) to speed up the mesh segmentation as well as forward (lead field) and inversion computations. Since GPU can speed up these processes substantially, it is used as a default setting. GPU-based FE computations allow producing a higly accurate FE mesh and the corresponding lead field matrix in a few minutes. For generating and processing a high-definition 1 mm resolution FE mesh, it is recommendable to use a computer equipped with at least 16 / 32GB RAM for EEG / MEG, respectively, together with a GPU with at least 4GB RAM. The GPU is set by default to a moderate parallelization level. The level of parallelization and, thereby, the GPU performance can be controlled via the ParallelVectors parameter (default 50) which can be set in the zeffiro_interface.ini file (the larger the value the higher the memory consumption). The GPU device number can be set in the same file. If a GPU cannot be applied, then it is recommendable to use a considerably lower number of parallel vectors, for example five. While running the interface, the GPU device can be set in the Miscellaneous options dialog box of the Edit menu.
The following reference computation times:
- mesh generation 1329 s,
- EEG leadfield 2362 s (128 electrodes, relative residual 1E-06),
- MEG leadfield 4826 s (154 magnetometers, relative residual 1E-06),
- source Interpolation 212 S,
were obtained for 1 mm resolution 6-compartment mesh with 36M elements, 6M nodes and ~0.5M sources using Lenovo P910 ThinkStation equipped with 2 x Intel Xeon E5-2697Av4 (RAM 256 GB) and 2 x NVIDIA Quadro P6000 (RAM 24 GB).
An accurate tetrahedral FE mesh can be created based on a polygonal (triangular) surface segmentation of the head and its compartments. To create the FE mesh, proceed as follows (see also this video):
(1) Go to the home folder and start Zeffiro with the command zeffiro_interface.
(3) Start a project by defining the Sensor locations and (magnetometer) orientations together with the points and triangles of the surface meshes via the corresponding push-buttons in the main window. Each layer can be fine-tuned through rotation and scaling. Any changes to the segmentation, will be automatically committed when starting the mesh generation process and can be also obtained by clicking the Apply transform button.
(4) Check that each tissue layer will have the desired conductivity value (Sigma). The mutual segmentation order of the compartments can be controlled via the Priority option.
(5) Visualize the surfaces by pushing the Visualize surfaces button (Figure 2). The Visible check box for each tissue layer will determine whether or not the layer in question will be shown. Checking the Cutting plane check box will add a cutting plane. The drop-down menu Visualization should be set to Sigma, meaning that the surface structure of the volume conductor model is examined.
Figure 2: Triangular surface meshes cut by an axial plane. The compartments correspond to the color labels given in the bottom part of the image.
(6) After defining the surface segmentation, the volumetric (tetrahedral) FE mesh can be generated by pressing the Volume mesh button. The mesh resolution (tetrahedron size) can be defined via the Mesh resolution. If the Mesh smoothing box is checked, the surfaces can be smoothed by some amount, increasing the smoothness of the lead field. The Smoothing strength value controls the smoothing strength (the larger value the stronger smoothing). Once the FE mesh has been generated, the conductivity values can be updated, if needed, by pressing the Update sigma button.
(7) When the FE mesh is ready, the result can be visualized with the Visualize volume button with Sigma (conductivity structure) as the Visualization mode. Especially, if the structure and thickness of the grey and white matter layer is not of a satisfactory quality, then it is recommendable to adjust the Priority (mutual segmentation order) of these layers. An example of a volumetric mesh is given in Figure 2.
(8) After finishing, save the project.
Note: When segmenting models with thin layers (e.g., skull, scalp and grey matter), it is preferable to set Meshing accuracy to one in order to achieve the best possible perfomance. Using a value between 0 and 1 will speed up the segmentation process, but can also reduce the quality of the outcome for thin layers. The default value is 0.1.
Figure 3: A volumetric tetrahedral mesh with 1 mm resolution obtained based on the surface segmentation of Figure 2.
The next step is to compute the lead field matrix (see also the video).
(1) Start by choosing the lead field type (EEG or MEG) in the Imaging method drop-down list.
(2) Give the desired number of sources in the Source count box. This is a rough estimate for the eventual number of individual current preserving sources placed in the grey matter compartment. In order to achieve a sufficient source coverage for the whole grey matter compartment, the source counts needs to be large, for example 0.1M or above, especially, if a high resolution is used.
(3) Choose in the Source directions drop-down list whether the source orientations will be Cartesian or mesh-based. In the latter case, the virtually random orientations of the FE basis functions will be used. The mesh-based sources constitute very precise dipole approximations and can, therefore, help to improve the inversion accuracy, depending on the applied inverse approach (Miinalainen and Pursiainen, 2017).
(4) Press the Lead field button to start the computation.
(5) If the LF source interpolation box is checked, the resulting lead field will be interpolated to enambe visualization of inverse estimates (reconstructions). The interpolation can be started also separately by pressing the Source interpolation button when the lead field matrix is ready.
(6) Once the lead field matrix has been computed, it is a good idea to save the project.
The Make all button will run the FE mesh and lead field matrix generation routines together with the LF source interpolation procedure, allowing one to perform the whole forward simulation process with a single click. The simplest way to obtain a forward simulation is to run Make all after importing a segmentation as a set of ASCII files as described below.
ZI allows one to import a surface segmentation consisting of a set of ASCII files.
For importing, the import fole (.ZEF) file should be placed in the folder containing
the segmentation and opened via the menu item Import ASCII segmentation from folder.
The filescan be either DAT files containing either points or triangles, or ASC
files exported from the FreeSurfer Software Suite.
In the former case, the folder must contain two files per each triangular tissue surface mesh (filename_points.dat and filename_triangles.dat), whereas in the latter case a single file is needed (filename.asc) per a mesh. Each line in the list below corresponds to a single mesh. Each compartment in the segmentation is described by one or more meshes which will be automatically merged in the import process. The compartment identifiers are the following:
sensor_points, sensor_directions, white_matter, grey_matter, csf, skull, skin, detail_1, detail_2, ..., detail_22.
Of these, a mesh for each tissue compartment is specified by a single comma-separated line of the following form:
filename, compartment name, scaling, sigma, priority, activity, name, invert, extension
Here, the filename appears without any extensions; compartment_name is as in the list above, scaling, sigma and priority parameters are as in ZI's segmentation window with 0 corresponding to the default value; activity is a number describing the activity of the compartment (0 = inactive, 1 = constrained activity, 2 = unconstrained activity, or 3 = inner cortex [for white_matter only]); name is the compartment name as it appears in ZI; invert is for inverting an inward-pointing surface normal (0=not inverted, 1=inverted); and extension is either ASC (asc) or DAT (dat) for FreeSurfer and other meshes, respectively.
In the special case of sensor_points and sensor_directions, it is only possible to use the DAT format. Each one of these is imported via a single file (filename.dat), and the line for importing is of the form:
filename, compartment name, scaling, x-translation, y-translation, z-translation, xy-rotation, yz-rotation, zx-rotation
The scaling, translating and rotation parameters are as they appear in ZI's segmentation tool, and selecting 0 for them, means that the default value will be used.
Any compartment-specific parameter can be imported only once. In case there are multiple definitions in the list for a single compartment, only the first one counts. The parcellation tool will distinguish the meshes merged into a single compartment displaying them as part 1, 2, etc. in the respective order.
(1) Import data. This can be done either via the Ìnverse tools menu or by directly associating zef.measurements with a dataset. The number of rows in the measurement array should coincide with that of the sensors. The number of columns is considered as the number of the recorded time steps. If the dataset is a cell array, then the cell entries are considered as different segments.
(2) Open the ÌAS MAP Estimation dialog box.
(3) Set the parameters (Sampling frequency, Time interval start, Low-cut frequency, High-cut frequency, FFT time window and Time step) to match with the dataset and the investigated frequency range. Choose the desired Data segment.
(4) Set the Hyperprior (either Gamma or Inverse gamma). This means fixing the outlier model (hyperprior) for the prior variance (Calvetti et al., 2009). Choose Shape parameter and Scaling parameter for the hyperprior. The former one of these controls the shape of the hyperprior (the strength of the outliers) and the latter one sets the initial prior variance. Note: If the scaling parameter is 1.5 and the gamma hyperprior is used, the reconstructions will correspond to the classical Minimum Norm Estimate (MNE) and Minumum Current Estimate (MCE), when the number of the iteration steps is 1 and >1, respectively.
(5) Set the Likelihood STD (likelihood standard deviation) to match the estimated noise level. The value is relative to Data normalization.
(6) Choose the desired number of iterations (the more steps the more focal solution).
(7) Press start. The reconstruction will be computed for each time step in the dataset.
(8) Visualize the reconstruction either on the surface segmentation or in the volumetric mesh. The type of the visualization can be chosen in the Visualization drop-down menu. The Snapshot/Movie button lets you print the image/time-lapse of the reconstruction into a file. The video codec to be applied can be defined in the file zeffiro_interface.ini. Example visualizations can be found in Figures 3-5.
Figure 3: A reconstruction produced through two IAS MAP iteration steps. Here, the normal component of the activity is visualized on the pial surface of the original segmentation.
Figure 4: The reconstruction of Figure 3 visualized on the volumetric FE mesh (1 mm resolution) with an axial cut on the top.
Figure 5: The reconstruction of Figure 3 visualized on FE mesh with a deeper axial cut together with a coronal cut. The lateral ventricles are also visible.
In ZI, a mesh imported from FreeSurfer can be equipped with a parcellation. A parcellation can be imported using the Parcellation tool. A colortable and points (labels) produced by FreeSurfer can be imported via Colortable and Points buttons. Colortable refers to a .MAT file containing the vertices, label, and colortable variables generated by FreeSurfer's Matlab interface function read_annotation. The points refer to an .ASC file produced by FreeSurfer's mri_mergelabels function. Several parcellation segments (for the left and right hemisphere) can be merged. The segmentation including all of its active subcompartments can be included in the parcellation by pressing the Segmentation button.
After importing the segmentation, one can select some or all of the regions within the parcellation list (Figure 6) and interpolate them into the source basis of the lead field by pressing Interpolate. The intepolated regions can be visualized on the brain by selecting the option Parcellation for both the colormap (in Options) and the visualization type (in Mesh tool). To visualize a parcellation or to restrict an inverse estimate to the selected regions the status of the parcellation must be Active (Figure 6).
Time series will analyze any inverse estimate in the selected regions (a single one or a sequence) with respect to the selected reconstruction type. The results, e.g., correlation, can be plotted by pressing Plot (Figure 1). In the Options menu, one can select whether a point-wise value or a quantile taken over a region is considered.
Figure 6: Parcellation tool in the Active mode.
Figure 7: A parcellation interpolated on a 1 mm brain model and visualized with the Parcellation option (for colormap and visualization type).
Figure 8: Time series correlation calculated for the selected brain regions after reconstructing the brain activity.
ZI allows using electrodes, magnetometers and magnetic field gradiometers as the sensors. A cap of N electrodes can be defined with a file containing an N-by-3 array of points. Defining a set of N magnetometers will necessitate defining an additional N-by-3 array representing the sensor orientations. Further, a set of N gradiometers can be defined via N-by-6 array in which the the first three columns correspond to the measurement orientations and the following three to the directions in which the gradient is taken.
In addition to the classical point electrodes, ZI allows one to use also the complete electrode model (CEM). A file defining N CEM electrodes will need to include an N-by-6 array in which the first three columns define the points and the following two the outer and inner circle radii of a ring electrode in millimeters and the last one thecontact impedance in Ohms. That is, each line should be of the form [x_coord y_coord z_coord outer_radius inner_radius impedance] with (x_coord, y_coord, z_coord) specifying the center point of the electrode.
Figure 9: A head model with CEM ring electrodes attached on the skin.
Choosing EIT as the imaging method in the main panel, Zeffiro allows computing and inverting a lead field for electrical impedance tomography (EIT) which can be used, e.g., in stroke/hemorrhage detection and tracking. EIT necessitates using the complete electrode model (here a ring electrode model). Ring electrodes can be used also for EEG in a similar manner. Again, the reconstruction for EIT (Figures 10 and 11) can be computed using the same inversion routines as for EEG/MEG. Notice that visualizing an EIT reconstruction necessitates choosing Basis as source directions in the main window and Linear scale together with the Value as the visualization type in the Miscellaneous options window.
Figure 10: A conductivity distribution including a synthetic anomaly (15 mm hemorrhage) which is depicted by the grey sphere.
Figure 10: An EIT reconstruction of the synthetic anomaly (hemorrhage).