Zeffiro User Interface for Electromagnetic Brain Imaging: a GPU Accelerated FEM Tool for Forward and Inverse Computations in Matlab

This article introduces the Zeffiro interface (ZI) version 2.2 for brain imaging. ZI aims to provide a simple, accessible and multimodal open source platform for finite element method (FEM) based and graphics processing unit (GPU) accelerated forward and inverse computations in the Matlab environment. It allows one to (1) generate a given multi-compartment head model, (2) to evaluate a lead field matrix as well as (3) to invert and analyze a given set of measurements. GPU acceleration is applied in each of the processing stages (1)–(3). In its current configuration, ZI includes forward solvers for electro-/magnetoencephalography (EEG) and linearized electrical impedance tomography (EIT) as well as a set of inverse solvers based on the hierarchical Bayesian model (HBM). We report the results of EEG and EIT inversion tests performed with real and synthetic data, respectively, and demonstrate numerically how the inversion parameters affect the EEG inversion outcome in HBM. The GPU acceleration was found to be essential in the generation of the FE mesh and the LF matrix in order to achieve a reasonable computing time. The code package can be extended in the future based on the directions given in this article.


Introduction
This article introduces the Zeffiro interface (ZI) [1] for electromagnetic brain imaging (Figure 1).The name Zeffiro is Italian for a 'gentle breeze'.ZI is aimed as an accessible and fast tool for utilizing the finite element method (FEM) in the electromagnetic (soft field) brain imaging applications in which the forward model follows from the Maxwell's equations [2].FEM is a volumetric modeling technique for solving equations in a bounded domain, such as the brain and the head [3,4].FEM allows discretizing the three-dimensional tissue parameter distributions, such as the conductivity, in an accurate way, taking into account anisotropic structures such as the fibrous white matter of the brain.In addition to the volumetric features, the finite element (FE) mesh can be accurately adapted to the complex internal surface layers of the head.The FEM can be applied to model an electromagnetic source within the grey matter (GM) of the brain in electro-/magnetoencephalography (EEG/MEG) [5,6], where the neural activity of the brain is being sensed remotely by measuring the electromagnetic field evoked by the activity.Consequently, via FEM one can compute an EEG/MEG lead field matrix which maps a given candidate solution to the expected measurement outcome.Furthermore, the same quasi-static set of Maxwell's equations that predicts the electric potential field of a neural source can be also approached via external magnetic field or current injections, where direct or alternating currents applied through electrodes act as the source of the electromagnetic field.Such applications include, for example, transcranial electric and magnetic stimulation (tES and tMS) [7,8] and electrical impedance tomography (EIT) [9].In tES, the brain activity is evoked through external stimuli, whereas in EIT the goal is to map the conductivity distribution within the computation domain.In contrast to EEG/MEG and tES/tMS, EIT constitutes a non-linear inverse problem where generating a lead field (LF) matrix requires linearizing the forward model with respect to a given current distribution.FEM is a powerful tool in EIT, since it allows generating forward simulation for an arbitrary volumetric conductivity distribution, in contrast to the Boundary Element Method (BEM) [10] which the dominating method in the brain applications, but assumes that the conductivity is a layer-wise constant parameter.However, until recently, FEM computations have been considered as computationally expensive regarding the practical usage.
ZI is an attempt to provide a simple and accessible and multimodal open source platform for FEM based forward and inverse computations within a given brain and head geometry.It has been implemented for the Matlab (The MathWorks Inc.) environment.In a typical workflow, one first imports a set of triangular tissue surface layers and after that (1) generates the FE mesh for a given multi-layer head model, (2) forms a basis for computing a LF matrix and inverting a given set of measurements, and (3) visualizes the reconstructions on the FE mesh.To tackle the issue of the high computational cost, ZI uses graphics processing unit (GPU) acceleration for each of the processing stages (1)-( 3).In its current configuration, it includes forward solvers for E/MEG and linearized EIT as well as a set of inverse solvers utilizing the hierarchical Bayesian model (HBM).The ZI platform and function library has been designed to be easily expandable and to allow implementing basically any FEM based linear or linearized forward simulation which can be formulated as a product between a LF matrix and a candidate solution vector.
In this paper, we describe the mathematics behind ZI, describe the principal operations and usage, concentrating especially on EEG/MEG, and introduce some central points for the developer perspective.We also show computed examples and analyze numerically how the inversion parameters should be selected, when inverting brain activity with HBM.

Forward modeling
ZI employs the divergence conforming source model in which the primary current density belongs to the space H(div) = { w|∇ • w ∈ L 2 (Ω)}.The model is current preserving, meaning that the source field does not have local sources or sinks.The total current density i.e., the sum of the primary current and volume current terms, has zero divergence which implies that The boundary condition, (σ∇u) • n = 0, given on ∂Ω, ensures that there is no current flow through the boundary of the domain.This equation multiplied by v ∈ H 1 (Ω) and integrated by parts results in the weak form which has a unique solution, assuming that the ground level has been fixed [11,12].

Discretization
The electrical potential field u and the primary current density J P are approximated as the finite sums u h = N i=1 z i ψ i and J P h = K j=1 x j w j , where ψ 1 , ψ 2 , . . ., ψ N ∈ H 1 (Ω) and w 1 , w 2 , . . ., w K ∈ H(div) [13].The discretized fields u h and J P h correspond to vectors z = (z 1 , z 2 , . . ., z N ) and x = (x 1 , x 2 , . . ., x K ) which satisfy the linear system in which the system matrix A is of the form A ∈ R (N ×N ) , A i,j = Ω ∇ψ j • (σ∇ψ i ) dV and the matrix G on the right hand side maps the primary current density on its divergence

LF matrix
The linear forward model is of the form where L is the LF matrix and n is a noise vector.In EEG, the LF matrix is of the form L = RA −1 G, where R is a restriction matrix.In the simplest case, R is a restriction to point values, i.e., the nodes in the FE mesh presenting the measurement point set.In addition to point electrodes, ZI allows one to use the lumped version of the complete electrode model (CEM) [14] in which the entries of R follow from integral averages over the electrode patches.
Computationally, the LF matrix can be obtained by evaluating a transfer matrix of the form T = A −1 R T .The LF and transfer matrix for MEG can be formed through the Biot-Savart integral formula for the total current field in a similar manner as for EEG as shown in [15].The relationship between the EEG forward simulation, the CEM and linearized EIT has been briefly explained in Appendix A.

Implementation in ZI
ZI utilizes the second-order version of the H(div) source model [13] in which both linear and quadratic basis functions are used for the primary current.In [16,13], this model has been shown to surpass the accuracy of the classical direct source modeling approaches based on the partical integration and St. Venant's priciple and to be especially advantageous for thin cortices and inversion of distributions.
In the Cartesian orientation mode, each source position includes three sources, one for each Cartesian coordinate direction.The Cartesian orientations are obtained from the mesh-based ones using the Position Based Optimization (PBO) method [17] with an adaptive [18] 10-source stencil in which 4 face and 6 edge functions are applied for each tetrahedron containing a source [13].Alternatively, the modeling can be restricted to the linear Whitney functions for which the 4-source stencil (4 face functions) can be applied to obtain the Cartesian set of orientations.Also, a plain set of Whitney functions can be used, that is, the LF matrix can be formed directly based on the mesh-based set of Whitney basis functions as suggested in [18].Furthermore, the sources can be oriented normally to the cortex surface according to the physiological knowledge of neural activity [19,20].To accomplish this, one can either fix the (pre-inversion) source basis to the normal orientation, or to find the positive or negative normal component of the neural activity after the inversion process via projection.
In ZI, the source positions are distributed randomly (uniformly) in the GM compartment and can be extended also in the white matter (WM), if necessary, e.g., if the deep brain structures are not included in GM.

HBM
The main inversion utility in ZI utilizes the HBM [21,22] to find a focal reconstrution of the activity.In HBM, a posterior probability for x is defined via choosing the standard deviation (STD) of a Gaussian likelihood density, the hypermodel, i.e., the gamma (G) or inverse gamma (IG) hyperprior determining the actual prior, and the shape and scaling parameters β, θ 0 of the hyperprior.
For a given measurement vector y, the Bayes formula [23] for the posterior probability density of x is of the form where p(x) is the prior density, and p(y | x) the likelihood function, following from the measurement noise model [24].It is assumed that the noise term is a Gaussian zero mean random vector n = y − Lx with independent entries.Consequently, the likelihood is of the form p(y | x) ∝ exp(−(2σ 2 ) −1 Lx − y 2 ), where σ is the noise standard deviation.In HBM, the prior can be expressed in the following hierarchical form where θ is the primary hyperparameter of the model.The conditional part p(x | θ) of the prior is a zero mean Gaussian density, whose diagonal covariance matrix is predicted by the hyperprior p(θ).The hyperprior is assumed to have a long-tailed density, implying that x is likely to be a sparse vector corresponding focal activity.
In ZI, it is either G or IG density [21], whose shape and scaling are controlled by the parameters β and θ 0 , respectively.ZI allows approximating both the posterior maximizer, i.e., maximum a posteriori (MAP), and conditional mean (CM) of the posterior density.The MAP estimate can be obtained via the iterative alternating sequential (IAS) method [21] for the whole domain at once, for a region of interest (ROI), or multiple different resolution levels.Of these, the third strategy can be applied to enhance the visibility of the deep activity.The basic version of the IAS method can be found in Appendix B. The CM estimate can be computed within a ROI using the Markov chain Monte Carlo (MCMC) sampling.

Hardware requirements
ZI is principally designed to be used with a workstation or a high-end desktop computer with tens of gigabytes of RAM, a multi-core CPU and one or more GPUs.A project after generating the FE mesh and LF matrix is supposed consume several gigabytes of RAM.Achieving a one-millimeter accuracy in forward modeling will presumably require ≥ 64 GB of RAM during computations with the total project size, when saved in a file, being around 0.5-1 GB.

GPU operations
ZI utilizes a GPU to accelerate the FE mesh generation process, forward and inverse computations, source interpolation and decompositions, as well as to speed up visualizing the reconstructions and the computation geometry in 3D.This is vital in order to achieve a feasibly short, around one-hour, computation time for a one-millimeter FE mesh resolution which has been shown to be essential in order to achieve physiologically accurate inverse estimates [25].A GPU is a parallel processing unit which is somewhat limited with respect to the RAM in comparison to the motherboard.It can handle computation intensive operations very effectively, while memory intensive operations need to be avoided.The operations related to forward and inverse computations can be accelerated due to the fast processing of matrix-vector products in a GPU, and the other operations are mainly based on the acceleration of find and sort operations and evaluating those as blocks instead of singular entries.A script example of a loop performed in blocks with an adjustable size for optimal GPU-specific performance is illustrated in Figure 2.

Forward simulation
In the Matlab environment, the most essential speed up gain is related to the sparse FE matrix-vector products which need to be evaluated iteratively in the forward simulation stage.The GPU parallelization of the forward simulation is especially important because Matlab currently handles the sparse matrix products in a single processor thread.Evaluating a single column  Here the for loop has been implemented in blocks for fast processing in a GPU.The block size, i.e., the number of vectors that are processed in parallel, can be controlled in order to obtain the best possible GPU-specific trade-off between parallel processing and memory consumption.
t of the transfer matrix T necessitates solving a linear system of the form At = r T , where r is a single row of the restriction matrix R. To find the solution, ZI uses the preconditioned conjugate gradient (PCG) [26] method with a lumped diagonal preconditioner LDP in which the each diagonal entry is obtained as a row sum of the absolute entry values.LDP is an advantageous preconditioner regarding the limited GPU memory.While LDP is not optimal with respect to minimizing the iteration steps needed for convergence, it enables establishing a fast forward solver due to a GPU's high computational capacity.

IAS MAP iteration
The IAS MAP estimation algorithm can be found in Appendix B. In the iteration, the most time consuming step is the third one (B.2) in which the size of the matrix to be inverted is determined by the length of the data vector.If a high number of time steps will need to be processed, the fastest inversion is obtained by evaluating the matrix-vector product of (B.2) in a GPU.

Data structures
When started, 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.Consequently, the LF and other variables will be easily accessible for the user.The export functions enable saving separate variables, for example, the LF matrix which, among other things, enables merging EEG and MEG LF matrices for parallel EEG/MEG.The main window contains the essential functions needed for mesh and LF generation as well as for 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: 1.The lead field matrix zef.L, 2. The measurement data zef.measurements(a matrix or a cell array with the number of rows and equal to that of zef.L and the time steps in the dataset, respectively), 3. source positions zef.source positions (source positions corresponding to the columns of zef.L in the respective order), 4. Source directions zef.source directions (an empty array, if Cartesian orientations are used, meaning that the source location/position for the columns of zef.L is given by the following pattern: position 1, xyz; position 2, xyz; position 3, xyz, ...), 5.The reconstruction is stored in zef.reconstruction.
The reconstruction can be visualized using the interface.For creating a new inversion tools menu item, one needs a .figfile and two .mfiles for initializing and updating the parameters.The workflow to produce a reconstruction with ZI is depicted as a diagram in Figure 3.The user first (1) imports the triangular surface model of each tissue type by loading the .datfile of the points and triangles in the main window (Figure 4).In the current version, up to nine different tissue types can be included in a single model.The default conductivity values for WM, GM, cerebrospinal fluid (CSF), skull, and scalp are 0.14, 0.33, 0.0064, and 0.43 S/m, respectively, according to [27,28].The conductivity distribution produced by the mesh generator is currently isotropic.However, the function for computing the LF matrix is compatible also with an anisotropic conductivity tensor.

Workflow
Then, (2) a uniform tetrahedral mesh is generated based on the surface segmentation.The process is robust and allows the tissue boundaries to intersect each other which is a typical situation with a real segmentation based on magnetic resonance imaging (MRI) data.The user can set the meshing priority (number) of each tissue compartment.A tetrahedron with nodes in several compartments will be associated with the one with the highest priority (smallest number).The ability to change the priority is essential in order to avoid holes in the thin skin, skull and GM compartments which is important regarding the modeling and inversion accuracy.The FE mesh can be also smoothed using the bi-Laplacian smoothing flow [29,30].
After the FE mesh generation phase, (3) the LF matrix can be computed for a given number of source positions.Additionally, in order to enable visualizing the reconstructions, (4) an interpolation process connecting the source positions and the FE mesh nodes will need to be run.Finally, after importing data, (5) the actual reconstruction can be computed and visualized with the inversion tools and visualization functions included in ZI

IAS MAP estimation process
The window of the IAS MAP estimation tool is shown in Figure 5.A reconstruction can be obtained as follows: 1. Import data.This can be done either via the Import 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.for the prior variance.Choose the Shape 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.If the shape parameter β is 1.5 and the G hyperprior is used, the reconstructions will correspond to the classical minimum norm estimate [6] and minimum current estimate [31], when the number of the iteration steps is 1 and > 1, respectively.Similarly, with IG there is a correspondence to the minimum support estimate [32]. 5. Set the Likelihood STD to match the estimated noise level.The value is relative to the 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 of the main window (Figure 4).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.

Inversion of distributions
To demonstrate the forward and inversion modeling capability of ZI, a 0.85 mm resolution FE mesh was generated for the WM, GM, CSF, skull, and scalp surfaces of the population head model1 [33] using ZI's default conductivity value set.A total of 72 ring electrodes with 1 kOhm impedance and outer and inner diameter of 10 mm were attached on the scalp.
The tetrahedral FE mesh obtained is shown in Figure 6.The total number of the elements and nodes in the FE mesh generated was 30 M and 5.2 M,   respectively.Using the FE mesh, an LF matrix was produced for 10000 source locations and Cartesian source orientations, and the sources were connected to FE mesh nodes via interpolation.

Hypermodel and parameter selection
To find a guideline for selecting the hypermodel and scaling parameter for the IAS MAP estimation process in EEG/MEG, we compared the localization of a simultaneously active pair of deep and superficial 10 NAm sources.The reconstruction was found as a center of mass of the primary current distribution within two 30 mm ROIs centered at the source locations.The accuracy was measured via evaluating the positioning and orientation (angle) difference with respect to the exact sources.
As the computation domain we used a six-compartment (WM, GM, CSF, Compact Skull, Spongious Skull, Skin) head model of a 49-year old male subject with default conductivity values.For the spongious part of the skull 0.028 S/m was selected [28].The EEG LF matrix was formed for a cap of 72 electrodes and the one for MEG was based on a helmet of 271 radial magnetometers.As the scaling parameter, we tested both 1E-5 and 1E-9.The shape parameter β was set to be 1.5 in each case.The relative STD of the (white) noise term was either 2 or 5 %.The details of the reconstructions computed can be found in Table 3.4.Each reconstruction was evaluated for 50 different realizations of the noise vector.The resulting samples were analyzed via box-plots.

Inversion of distributions
The total time of the mesh generation, LF matrix evaluation and interpolation processes were 21, 39 and 3.5 min, respectively, with NVIDIA Quadro P6000 [34] GPU.Using GPU was also found to be necessary to obtain a feasible computation time as GPU accelerated these routines by more than a factor of ten.An example a multiresolution IAS reconstruction obtained using synthetic data with 3 % relative (white) noise STD corresponding to a simultanously active dipolar thalamic (e.g.N14) and somatosensory (e.g.N20/P20) component is illustrated in Figure 7. Since the thalamus is not included in the GM compartment of PHM, the sources positions cover both GM and WM.A reconstruction for somatosensory activity obtained when restricting the sources into the GM can be found in Figure 8 which includes, additionally, also a IAS MAP reconstruction of a conductivity anomaly produced with ZI using a linearized EIT model and synthetic data (Appendix A) within a spherical geometry.

Hypermodel and parameter selection
Figures 9 and 10 show comparisons between the exact and reconstructed sources for EEG and MEG.The box-plot analysis of the reconstruction accuracy can be found in Figures 11 and 12, respectively.Table 3.4 lists the hypermodel and scaling parameter corresponding to the smallest obtained median position and angle difference for the deep and superficial position.The obvious tendency is that the superficial source can be localized more accurately than the deep one.In EEG, the best results were obtained with the value 1E-5 and 1E-9 for hyperpriors G and IG, respectively.Regarding the superficial source, G was observed to yield the superior results for the position, whereas IG was found to be advantageous with respect to the angle.For the deep part, IG with θ 0 = 1E-9 was observed to be preferable regarding both quantities.In MEG, G yielded the most accurate position estimate in each case.IG performed appropriately for angle in the superficial part, whereas for deep part, the best orientation accuracy was obtained with IG and G for the noise levels 2 % and 5 %, respectively.The general observation for both EEG and MEG was that the low scaling parameter value θ 0 = 1E-9 was preferable for localizing the deep source, especially, with the lower 2 % noise level.

Discussion
This paper has introduced the Zeffiro interface (ZI) [1] for brain imaging.ZI provides a GPU accelerated approach to using the FEM [3,35] in applications involving forward modeling and inversion of electromagnetic fields [2].The central aim of ZI is to provide an accessible and fast modeling tool for the FEM.When aided by a GPU, ZI allows one to invert a given set of EEG data for a physiologically accurate [25] sub-one-millimeter volumetric multi-compartment headmodel within a reasonable one hour's time.
GPU acceleration was found to be necessary, specifically, in the generation of the FE mesh and the evaluation of the LF matrix in order to achieve a feasibly short computaion time.This is especially important, since Matlab does not currently parallelize the sparse matrix operations in a CPU which is why the performance difference between CPU and GPU computations, both applicable in ZI, is especially pronounced.A standard laptop GPU was observed to be generally not sufficient for high-resolution forward simulation regarding its speed and memory, whereas a workstation GPU performed decently.The workstation GPUs utilized in this study (Nvidia Quadro P6000, P4000 and M40000 [34]) had at least 8 GB of RAM.The minimum requirement is 2 GB based on the maximal memory consumption during the computations.
Consequently, for generating the FE mesh and the LF matrix, the office is suggested to be equipped with at least one workstation or powerful desktop as mentioned above.A laptop was, however, found to be sufficient for the inversion process.An important property, in this regard, is the ability to save the zef struct, including all the variables, the LF matrix, and FE mesh of ZI, as a single .matfile, so that the project can be shared and utilized for further calculation easily.Thereby, multiple neuroscientists will be able to use and make further calculation or analyze on the single file on different computers for various purposes.Handling the variables as a single struct data type within Matlab was also found to be advantageous from the viewpoint of having access to the all the variables not only via the ZI's front-end, but also using the command prompt which allows the scientists and developers to use external functions easily together with ZI.
The central advantage of ZI's tetrahedral mesh generator is that, unlike many widely used software, e.g., TetGen [36] and Netgen [37], it allows the tissue boundaries to intersect each other without collapsing.Whereas an ordinary FE mesh generator attempts to optimize the elements with respect to both volumetric and surface geometry, and will, therefore, easily crash when meshing a set of thin or intersecting surface layers, ZI utilizes a robust uniform mesh structure and priority based compartment labeling which enables, via control of the priority parameter, optimizing the layer structures also after the actual meshing process.This is essential in practice, since the tissue segmentation routines utilizing MRI images often do not distinguish the different compartments smoothly.
The present computed examples show that ZI enables robust inversion of the primary current distribution, and that the resulting reconstruction can be visualized as a volumetric field, restricting the activity according the physiological normal constraint [19,20] if needed.Since the current is inverted and visualized using the FE mesh, the volumetric structure of the domain has a strong influence on the final reconstruction.For example, if the deep structures of the head model are absent which was the the case with the PHM [33] used this study.In ZI, the solution is to extend the sources to white matter which works appropriately based on the present results.Differences between source spaces extended over GM and WM, especially, regarding localization of deep sources will be an important future topic.
As shown by the results, the hypermodel and scaling parameter θ 0 (initial prior variance) selection for the IAS MAP estimation process will affect on the source localization accuracy in EEG/MEG.An advantageous choice for depth localization seems to be IG coupled with the value θ 0 = 1E-9.For the superficial areas, G with a slightly larger value, e.g., θ 0 = 1E-5 appears beneficial.The guidelines found in this study are in line with previous findings [21] and might be optimized later on in an application specific context.
The preliminary inversion results obtained with EIT, suggest that ZI can also be extended to non-linear problems and inversion of scalar-valued fields.Since the CEM and the ability to inject current patterns into the domain, which is necessary in EIT [14,9], have already been implemented in ZI, it would be relatively easy to add in solvers for other applications involving current patterns, such as tES.
Currently, Zeffiro is the only open source Matlab toolbox to use GPU assisted FEM for EEG/MEG, but there are many other software tools aimed at EEG/MEG neuroimaging.For instance, the Duneuro [38,39] and SimBio [40] are an open source libraries with similar functions as Zeffiro but utilizing C++ language.Brainstorm [41] and Fieldtrip [42] are both implemented for Matlab as well.However, the core forward modeling approach of Brainstorm is the BEM method, and the Fieldtrip does not have an advanced forward and inverse modeling functions.Moreover, none of these are currently capable of advanced FEM or GPU computations.The MNE-Python toolbox [43] is the leading option for Python which also allows utilizing a GPU.However, MNE Python is also limited to the BEM forward simulation.We emphasize that the ability to utilize the FEM in ZI is essential, to allow extending it to other inverse imaging modalities.
In the future work, we will investigate further the forward and inversion accuracy that is obtainable with ZI.An interesting direction will be to examine improved techniques to find a reconstruction with the HBM approach, e.g., the MCMC method restricted into a ROI.Implementing more inverse modalities for ZI and developing machine learning based analysis methods for it will also be potential directions.The data structure of the current Zeffiro can also be enhanced so as to make the computation faster and occupying less memory.
That is, the optimization is done alternatingly for the unknown vector x and the hyperparameter θ.The algorithm can be, further, written as 1.Set k = 0 and θ (0) = (θ 0 , θ 0 , . . ., θ 0 ). 2. Set L (k) = LD 5. Set k = k + 1 and go back to 2., if k is less than the total number of iterations defined by the user.

a u x v e c 1 =
( 1 / 3 ) * ( boundary p ( b o u n d a r y t ( : , 1 ) , : ) + boundary p ( b o u n d a r y t ( : , 2 ) , : ) + boundary p ( b o u n d a r y t ( : , 3 ) , : ) ) ' ; a u x v e c 2 = boundary p ( b o u n d a r y t ( : , 2 ) , : ) '− boundary p ( b o u n d a r y t ( : , 1 ) , : ) ' ; a u x v e c 3 = boundary p ( b o u n d a r y t ( : , 3 ) , : ) '− boundary p ( b o u n d a r y t ( : , 1 ) , : ) ' ; a u x v e c 4 = c r o s s ( a u x v e c 2 , a u x v e c 3 ) / 2 ; i n d v e c = z e r o s ( s i z e ( nodes , 1 ) , 1 ) ; o n e s v e c = o n e s ( l e n g t h ( a u x v e c 1 ) , 1 ) ; i n d v e c a u x = z e r o s ( l e n g t h I , 1 ) ; n o d e s a u x = nodes ( I , : ) ' ; par num = e v a l i n ( ' b a s e ' , ' z e f .p a r a l l e l v e c t o r s ' ) ; b a r i n d = c e i l ( l e n g t h I / ( 5 0 * par num ) ) ; i i n d = 0 ; f o r i = 1 : par num : l e n g t h I i i n d = i i n d + 1 ; b l o c k i n d = [ i : min ( i+par num −1 , l e n g t h I ) ] ; a u x v e c = n o d e s a u x ( : , b l o c k i n d ) ; a u x v e c = r e s h a p e ( a ux vec , 3 , 1 , l e n g t h ( b l o c k i n d ) ) ; a u x v e c 5 = a u x v e c 1 ( : , : , o n e s ( 1 , l e n g t h ( b l o c k i n d ) ) ) − a u x v e c ( : , o n e s v e c , : ) ; a u x v e c 2 = sum ( a u x v e c 5 .* a u x v e c 4 ( : , : , o n e s ( 1 , l e n g t h ( b l o c k i n d ) ) ) ) ; a u x v e c 3 = s q r t ( sum ( a u x v e c 5 .* a u x v e c 5 ) ) ; a u x v e c 3 = ( a u x v e c 3 .* a u x v e c 3 ) .* a u x v e c 3 ; a u x v e c 6 = sum ( a u x v e c 2 ./ a u x v e c 3 ) / ( 4 * p i ) ; i n d v e c a u x ( b l o c k i n d ) = a u x v e c 6 ( : ) ; end

Figure 2 :
Figure2:A script example of finding elements inside a given closed triangulated tissue layer surface based on determining the solid angle for a given point.The solid angle is zero for an exterior point and 4π for a point inside.Here the for loop has been implemented in blocks for fast processing in a GPU.The block size, i.e., the number of vectors that are processed in parallel, can be controlled in order to obtain the best possible GPU-specific trade-off between parallel processing and memory consumption.

Figure 3 :
Figure 3: A diagram for the workflow in ZI.

Figure 4 :
Figure 4: The main window of ZI.Illustrated is the 5.2 node and 30 M element FE mesh generated using the population head model (PHM).

Figure 5 :
Figure 5: The IAS MAP estimation tool.

2 .
Open the IAS 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 G or IG), i.e., the outlier model (hyperprior)

Figure 6 :
Figure 6: A comparison of between the surface model (left) of the brain and volumetric FE mesh (right) generated for the PHM using the meshing resolution of 0.85mm.

Figure 7 :EITFigure 8 :
Figure 7: Six examples of how a (multiresolution IAS) reconstruction can be visualized in ZI.The distribution can be visualized interpolated on the surface of the cortex or in the volume.It can be constrained to the positive or negative normal direction, or restricted to the tangential one.Here, normal and tangential refer to the orientation with respect to the visualized surface (either cut or not).In this case, two synthetic sources are simultaneously active: one (P20/N20) in area 3b in the posterior wall of the central sulcus and a deep (P14/N14) upward-pointed source.The first one of these is visible in the sulcus wall in normally constrained distribution and a deep component corresponding to the second one can be observed in the tangentially restricted volumetric cut.

Figure 9 :
Figure 9: Examples of the center of mass (red) found for the exact (cyan) synthetic deep and superficial source in the case of the synthetic EEG data.

Figure 10 :
Figure 10: Examples of the center of mass (red) found for the exact (cyan) synthetic deep and superficial source in the case of the synthetic MEG data.

Figure 11 :Figure 12 :
Figure 11: The distributions for position and angle differences in the source localization experiments in the case of synthetic EEG data.In each image the box-plot for the deep source is on the left and for the superficial one on the right.

Table 1 :
The parameters for the reconstructions computed in the EEG/MEG source localization experiment.

Table 2 :
The hypermodel and scaling parameter corresponding to the best median accuracy obtained in the EEG/MEG localization of a superficial and deep source.