Many devices used in the clinical or the cognitive science domain perform electromagnetic measurements, or stimulation, on the human body. Devices measuring electric fields include electroencephalograms (EEG), cardiograms (ECG), myograms (EMG), while magnetoencephalograms (MEG) or magnetocardiograms (MCG) measure magnetic fields. Among stimulating devices, transcranial magnetic stimulation (TMS) uses magnetic coils to stimulate brain regions, while functional electric stimulation (FES) and electrical impedance tomography (EIT) impose an electric current or potential through contact electrodes.
To interpret and control the bioelectromagnetic phenomena involved with these devices, realistic physiological modeling is required, in terms of geometry and conductivity [1]. Accurate numerical solutions of the governing equations must be computed: obtaining the best accuracy possible for a given computational model is one of the goals of OpenMEEG, the opensource software package introduced in this article.
Electromagnetic propagation is governed by the Maxwell equations, coupling the electrical and magnetic fields. This coupling simplifies when the relevant frequencies are low enough for the quasistatic regime to hold. The electric potential is then governed by the law of electrostatics
\nabla \cdot (\sigma \nabla V)=\nabla \cdot {J}^{\text{p}},
(1)
Where σ is the conductivity field, and J^{p} is a dipolar source distribution within the domain. When considering brain activations, it represents average postsynaptic currents within pyramidal cortical neurons. A boundary condition must be imposed, typically controlling the normal current on the domain boundary:
\sigma \nabla V\cdot n=j.
(2)
The electric potential can be computed independently from the magnetic field, by solving (1) with boundary condition (2). The magnetic field B depends both on the electric potential V and on the current source distribution J^{p}, through the Biot and Savart law.
B\text{(}r\text{)}=\frac{{\mu}_{0}}{4\pi}{\displaystyle \int ({J}^{\text{p}}(}{r}^{\prime})\sigma \nabla V({r}^{\prime}))\times \frac{r{r}^{\prime}}{\Vert r\text{\u2013}{r}^{\prime}{\Vert}^{3}}d{r}^{\prime},
(3)
written here in the case where j = 0 on the boundary.
A forward problem consists of simulating V and/or B when σ, J^{p}, and boundary current j are known. The forward electro magnetostatics problem is wellposed, in contrast to the far more difficult, illposed inverse problem of estimating σ, or J^{p}, from partial boundary data. Still, obtaining an accurate solution for the forward problem is far from trivial. This paper presents a software package, OpenMEEG, that makes available to the community recent research efforts to improve the accuracy of forward solvers.
Forward solutions provide the relationship between the quantities of interest and the measurements. To obtain a good description of this relationship, one must model the sources, the conductivity, and the sensors.
The choice of conductivity model is especially delicate because it strongly constrains the numerical solutions that can be used to solve the problem. In all generality, the conductivity field σ should be modeled as a tensor field, because composite tissues such as bone and fibrous tissues have an effective conductivity that is anisotropic. Realistic, anisotropic conductivity models can however be difficult to calibrate and handle: simpler, semirealistic, conductivity models assign a different constant conductivity to each tissue, as depicted in Figure 1. There are three main types of conductivity models, and associated numerical methods:

1.
if the conductivity field can be described using simple geometries (with axilinear, planar, cylindrical, spherical or ellipsoidal symmetry), analytical methods can be derived, and fast algorithms have been proposed that converge to the analytical solutions for EEG [2, 3] and MEG [4];

2.
for piecewise constant conductivity fields as in Figure 1, Boundary Element Methods (BEMs) can be applied, resulting in a simplified description of the geometry only on the boundaries [5];

3.
general nonhomogeneous and anisotropic conductivity fields are handled by volumic methods; Finite Element Methods and Finite Difference Methods belong to this category.
This paper deals with Boundary Element Methods, whose advantage over volumic methods is to use an economic representation of the conductivity field (a few conductivity parameters  one per tissue, and a few triangular meshes to represent the interfaces). Until recently, all Boundary Element Methods used in bioelectromagnetics were based on a Green representation theorem, involving operators called single and doublelayer potentials. In OpenMEEG the approach is to consider an extended version of the Green representation theorem [6], involving more operators, and leading to a new BEM formulation, called the Symmetric Boundary Element Method [7].
The structure of the paper is as follows: the applications targeted by OpenMEEG are first presented, then some mathematical aspects of the Symmetric BEM are explained (details are presented in an appendix). Then, to motivate the use of OpenMEEG, a comparison study with four other solvers in the context of EEG forward modeling is presented. The accuracy of the solvers is tested on multiple random sphere models. The accuracy of OpenMEEG is also assessed by numerical experiments for MEG forward modeling. Finally a section provides technical details on the OpenMEEG software package while giving sample code for using OpenMEEG via Python or via the Fieldtrip Matlab toolbox. A more complete presentation of the usage of OpenMEEG via a command line interface is available in Additional file 1.
Target applications of OpenMEEG
The main purpose of OpenMEEG is to solve the forward problems arising in magneto and electroencephalography. OpenMEEG is therefore primarily developed in the field of brain research, but its development could directly benefit to other problems dealing with electromagnetic biosignals (see for example [8] in the context of electrocardiography). A geometrical model of nested meshes representing tissue interfaces must be provided to OpenMEEG, which does not perform any segmentation or meshing. OpenMEEG handles nested regions, and could in principle be extended to more general, disjoint regions [9]. The main output of a forward problem is a leadfield, i.e., the linear application relating sources at specific locations to sensor measurements. Although the principal target of OpenMEEG is magneto and electroencephalography, other types of bioelectromagnetic problems have also been handled with OpenMEEG: we hereforth describe its current scope.
Electroencephalography (EEG) is concerned with the variations of electric potential on the scalp, due to sources within the brain. At frequencies of interest, the quasistatic regime is valid, resulting in the electrostatics relation (1). The air surrounding the scalp is supposed nonconductive, hence the normal current vanishes on the scalp: j = 0 in boundary condition (2). The sources within the brain are represented by dipoles: in equation (1), the sources are {J}^{\text{p}}=\overrightarrow{q}{\delta}_{p} where p is a dipole position and \overrightarrow{q} the associated dipolar moment. OpenMEEG computes the electric potential and the normal current on each interface between two homogeneous tissues, due to electric sources within the brain. EEG sensors are electrodes, modeled in OpenMEEG as discrete positions on the scalp at which the potential can be measured (infinite impedance assumption). On these sensors, OpenMEEG computes the EEG leadfield, representing the linear relationship between source amplitude (for fixed position and orientation) and sensor values.
Magnetoencephalography (MEG) is concerned with magnetic fields produced by sources within the brain, which, like for EEG, are modeled as dipoles [10]. OpenMEEG makes use of the Biot and Savart relation (3) to compute the magnetic field, and hence requires the electric potential to be computed beforehand on all the interfaces [11]. Magnetometers or gradiometers can be modeled in OpenMEEG. Magnetometers are defined by their position and the direction of the field they measure. Gradiometers and more general sensors are handled by providing to OpenMEEG the positions, the orientations and the weights of integration points. For example with axial gradiometers present in CTF MEG systems, the forward field for a sensor is obtained by subtracting the two leadfields computed at the locations of two nearby magnetometers.
Electrical Impedance Tomography (EIT) infers characteristics of a conductive domain, by analyzing the potential resulting from the application of a current on the boundary. This method has been applied to calibrate the conductivity of EEG head models [12–14]. OpenMEEG computes the forward problem associated to EIT: given a conductive domain Ω defined by the interfaces between homogeneous regions, and their conductivity, and for a prescribed normal current j, OpenMEEG computes the potential V and the normal current σ∂_{
n
}
V on each interface by solving (1) for J^{p} = 0, and with boundary condition (2). The electrodes are modeled with a P0 approximation over the triangles of the scalp; the injected current is assumed constant over a triangle.
Intracranial electric potentials are measured in certain clinical settings, either on the surface (ElectroCorticography) or within the brain (intracranial EEG, or stereotaxic EEG). OpenMEEG is able to compute leadfields for such intracranial electrodes, in realistic head models.
Functional Electrical Stimulation (FES) provides a way to restore movement of paralyzed body regions by activating the efferent somatic axons. For this a current is applied to a nerve sheath, using specially designed electrodes. The precise location and time course of the applied current must be optimized in order to achieve the best selectivity, and to minimize the current intensity for a desired outcome. Optimizing the stimulation parameter in a realistic nerve model requires a forward model for FES. OpenMEEG provides such a tool, by combination of the concepts of EIT and of internal potential simulation (see section on the applications targeted by OpenMEEG) [15, 16].
Cortical Mapping is an inverse problem that aims to recover the potential and the normal current occurring on the surface of the cortex (i.e., under the skull bone), given EEG measurements on electrodes [17]. A particularly elegant solution to this problem has been proposed with the symmetric BEM [18], making it possible e.g. to solve further source localization problems.
Methods: implementation
OpenMEEG uses a Galerkin Boundary Element formulation, that jointly considers the electric potential and the normal current on each interface as unknowns of the problem. A P1 (piecewise linear) approximation is used for the electric potential, whereas the normal current \sigma \frac{\partial V}{\partial n} is approximated with P0 (piecewise constant) boundary elements. The most intricate part of the implementation concerns the assembly of the system matrix, requiring singular kernel integration over triangles. Those are double integrals. The inner integrals are computed with analytical schemes [19], whereas the outer integrals are computed with a 7point Gauss quadrature scheme [20]. An adaptive integration scheme recursively subdivides the triangles until the required precision is achieved. This adaptive integration has an influence on the accuracy, as will be exposed in the benchmark results further on.
Since the electric potential can only be computed up to a constant, the system matrix is deflated to make it invertible. In practice, it consists of imposing the constraint that the integral of the potential over the external layer is 0.
The magnetic field is computed by using the Biot and Savart equation with the Galerkin piecewise linear approximation for the potential [11].
As stated above, given a set of dipole positions and orientations, a set of sensors and a head model defining homogeneous conductivity regions, the M/EEG forward problem produces as output the leadfield. OpenMEEG computes such lead fields with the following procedure (see the tutorial in Additional file 1 for a detailed graphical representation of this flowchart):

assemble the system matrix involving boundary integral operators on the discretized surfaces;

for a specified set of dipole positions and orientations, assemble a discretized, sourcerelated term;

solve the linear system relating the two above matrices (providing V and \sigma \frac{\partial V}{\partial n} on each interface)

interpolate the scalp potential (for EEG) or apply the discretized Biot and Savart relation (for MEG).
Mathematical details, as well as the practical usage of OpenMEEG can be found in the Appendix.
Methods: benchmark comparison study
In order to motivate the use of the OpenMEEG software, we have conducted a set of numerical experiments that compare the accuracy and the robustness of OpenMEEG with alternative M/EEG forward solvers. The comparison in the context of EEG forward modeling is run with the 4 alternative BEM solvers (see [21] for a comparison with FEM). The precision of the MEG forward solver is demonstrated using known analytical properties of the magnetic field when considering sphere models and with two other solvers.
Publicly available M/EEG forward solvers
Several software projects have the ability to solve the M/EEG forward problems. MNE, BrainStorm, EEGLAB (via the NFT Toolbox), Fieldtrip, Simbio, OpenMEEG and SPM, which shares with Fieldtrip the same M/EEG forward solvers.
Fieldtrip and SPM offer two implementations of the BEM. The first one, called Dipoli, was written by Oostendorp [22] and is not open source (only binary files for UNIX systems are available), while the second one, called BEMCP, is opensource and was written by Christophe Phillips during his PhD [23]. The Dipoli solver implements a linear collocation method with Isolated Skull Approach (ISA) [24], whereas BEMCP implements a simple linear collocation method. The Dipoli implementation details can be found in [22]. Another implementation with linear collocation (with and without ISA) can be found in the Matlab toolbox called Helsinki BEM [25]. The MNE package written in pure C code by Hämäläinen also offers a linear collocation (with ISA) implementation of the BEM. The Brainstorm toolbox in its latest version uses only sphere models for both EEG and MEG. Recently the EEGLAB toolbox also has provided a package called NFT that is based on a BEM called METU [26]. Finally, the Simbio forward solver also consists of a BEM with linear collocation and ISA.
Commercial software packages are not listed here. However, to our knowledge, commercial products that provide a BEM solver for forward modeling implement a linear collocation method with ISA. This is for example the case of ASA [27]. Note also that beyond the field of brain research, alternative BEM implementation exist. The SCIRun/BioPSE project contains for example a BEM implementation based on [19] that can solve the forward problem of electrocardiography [28].
Accuracy measures
The accuracy of forward solvers can be assessed for simple geometries such as nested spheres, by comparison with analytical results. The precision of a forward solution is tested with two measures. the RDM (Relative Difference Measure) and the MAG (Magnitude ratio) [29].
The RDM between the forward field given by a numerical solver g
_{
n
} and the analytical solution g
_{
a
} is defined as:
RDM({g}_{\text{n}},{g}_{\text{a}})=\Vert \frac{{g}_{n}}{\Vert {g}_{n}\Vert}\frac{{g}_{a}}{\Vert {g}_{a}\Vert}\Vert \in \left[0,2\right],
while the MAG between the two forward fields is defined as:
MAG({g}_{n},{g}_{a})=\frac{\Vert {g}_{n}\Vert}{\Vert {g}_{a}\Vert}.
In both of these expression, the norm is the discrete ℓ^{2} norm over the set of sensor measurements. The closer to 0 (resp. to 1) the RDM (resp. the MAG), the better it is.
Geometrical models
The comparisons were made both on classic regular sphere meshes as in Figure 2, and on random meshes. A random sphere mesh with unit radius and N vertices is obtained by randomly sampling 10N 3D points, normalizing them, meshing their convex hull and decimating the obtained triangular mesh from 10N to N vertices. This process guarantees an irregular meshing while avoiding flat triangles. The BEM solvers are tested with three nested sphere shells which model the inner and outer skull, and the skin. The radii of the 3 layers are set to 88, 92 and 100, while the conductivities of the 3 homogeneous volumes are set to 1, 1/80 (skull) and 1. In this benchmark, the units are arbitrary, but in practice, units should be expressed with the International System of Units (SI). For each randomly generated head model, it was tested that they were no intersection between each mesh. For every head model, solvers are tested with the same 5 dipoles positioned on the zaxis with orientation (1,0,1) and various distances to the inner layer (cf. Figure 2). As expected, the accuracy of the solvers decreases as this distance gets small.