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.
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 some central package features can be found below. Script examples can be found in the scripts folder.
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 its default windows open: segmentation, mesh, visualization, figure, and menu tool.
The original publication of the ZI code package is He et al. 2019. Package development was originally started as a project to combine high-resolution FE head models with hierarchical Bayesian inverse approaches. The essential mathematical techniques used in the interface have been reviewed and validated in various papers, for example, (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).
Figure 2: Examples of adapted high-resolution FE mesh structures with details smaller than 0.5 mm: a refined cortical structure for an improved source modelling accuracy in EEG (1st from left), refined sub-thalamic arteries (2nd from left), refined thalamus providing an enhanced forward simulation accuracy in the sub-cortical part of the brain (3rd and 4th from left).
When opened as a graphical user interface ZI creates a single project data structure (struct) zef in Matlab's workspace. All the parameters and variables, such as the lead field matrix, measurement data and reconstruction, can be accessed via this structure. The project can be also saved or loaded on a hard disk as a single file. This way, the access to lead field and other variables is straightforward and fast. The default set of windows - segmentation, mesh, visualization, figure, and menu tool - contains functions needed for mesh and lead field generation and visualizing the results. The items of the Inverse tools menu enable inversion of the data. For 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).
A reconstruction can be visualized using the figure and mesh visualization tool. 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.
Figure 3: Simulated EEG reconstructions depicting brain activity: the first one of the reconstructions includes a cone field and the second one an inflated set of cortical and sub-cortical surfaces.
The current version 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 0.5 - 1 mm resolution FE mesh, it is recommendable to use a computer equipped with at least 16 / 32GB RAM for EEG / MEG, respectively. A high-memory GPU is recommended but not necessary. 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). Similarly the parameter ParallelProcesses controls the number of processes run simultaneously in CPU. The GPU device number can be set in the same file.
The following reference computation times:
- mesh labeling 1329 s,
- EEG leadfield 2362 s (128 electrodes, relative residual 1E-06),
- MEG leadfield 4826 s (154 magnetometers, relative residual 1E-06),
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). In this test, a uniform multi-compartment mesh was generated, while generating an adapted FE 1 mm mesh can be expected to take 2 to 4 times as long.
The next step is to compute the lead field matrix (see also the video).
(1) A lead field generation routine can be run using the Run script button of the mesh tool. The list of different routines is given in the box above.
(2) Give the desired number of sources in the Source/Field count box. This is a rough estimate for the eventual number of individual current preserving sources placed in the grey matter and other active compartments. 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 Directions drop-down list whether the source orientation type will be Cartesian, Normal or Basis. In the last 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, might improve inversion accuracy, depending on the applied inverse approach (Miinalainen and Pursiainen, 2017). However, the Cartesian and normal sources can be assumed to give the best performance in most cases.
(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 once it is finished to enable visualization of inverse estimates (reconstructions). The interpolation is a quick process that can be run also separately by pressing the Source interpolation button when the lead field matrix is ready.
(6) Once the lead field matrix has been computed, the heavy forward modelling phase is done and it is a good idea to save the project. A project with a lead field matrix can be easily operated in a lower performance computer, for example, a laptop, whereas high-resolution volumetric simulations are the easiest to be run using a workstation.
The following gives an example on how to invert data using IAS Inversion, which is one of the inverse tools.
(1) Import measurement 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.
(2) Open the ÌAS Inversion dialog box.
(3) Set the parameters (Sampling frequency, Time interval start, Low-cut frequency, High-cut frequency, Time window and Time step) to match with the dataset and the investigated frequency range. Choose the desired Data segment.
(4) The hyperprior can be tuned in the Gaussian prior options menu, where either Gamma or Inverse gamma hypermodel can be chosen. This means fixing the outlier model (hyperprior) for the prior variance (Calvetti et al., 2009). The SNR and PM-SNR (in Gaussian prior options) will tune the Shape parameter and Scaling parameter of 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. For more about the role of these two parameters, see Rezaei et al., 2020
(5) Set the SNR (likelihood standard deviation) to match the estimated noise level. The value is relative to Data normalization.
(6) Choose the desired number of IAS MAP 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 type and Reconstruction type drop-down menus. The Snapshot/Movie button lets you print the image/time-lapse of the reconstruction into a file.
Figure 6: A reconstruction produced through two IAS MAP iteration steps. Here, the normal component of the activity is visualized as a contour plot on the pial surface of the original segmentation.
Figure 7: The reconstruction of Figure 6 visualized on the volumetric FE mesh (0.65 mm resolution on cortex) with an axial cut on the top.
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.
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 modelling can be done 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).