Parallel implementation of the accelerated BEM approach for EMSI of the human brain

Abstract

Boundary element method (BEM) is one of the numerical methods which is commonly used to solve the forward problem (FP) of electro-magnetic source imaging with realistic head geometries. Application of BEM generates large systems of linear equations with dense matrices. Generation and solution of these matrix equations are time and memory consuming. This study presents a relatively cheap and effective solution for parallel implementation of the BEM to reduce the processing times to clinically acceptable values. This is achieved using a parallel cluster of personal computers on a local area network. We used eight workstations and implemented a parallel version of the accelerated BEM approach that distributes the computation and the BEM matrix efficiently to the processors. The performance of the solver is evaluated in terms of the CPU operations and memory usage for different number of processors. Once the transfer matrix is computed, for a 12,294 node mesh, a single FP solution takes 676 ms on a single processor and 72 ms on eight processors. It was observed that workstation clusters are cost effective tools for solving the complex BEM models in a clinically acceptable time.

Appendix: BEM formulation

The electric potential ϕ and the magnetic field B due to a current dipole source p in a piecewise homogeneous volume conductor model of the head, satisfy the following integral equations [23]:

$$ \bar{\sigma}\phi({{\mathbf{r}}})=g({{\mathbf{r}}})+\frac{1}{4\pi}\sum_{k= 1}^{L} (\sigma_k^{-}-\sigma_k^{+}) \int_{S_k}\phi({{\mathbf{r}}} '){\frac{{\mathbf{R}}} {R^3}}\cdot {\rm d} {{\mathbf{S}}}_k({{\mathbf{r}}} '), $$

(7)

$${{\mathbf{B}}} ({{\mathbf{r}}})= {{\mathbf{B}}}_0({{\mathbf{r}}})+{\frac{\mu_0}{4\pi}}\sum_{k=1}^{L} (\sigma_k^{-}-\sigma_k^{+})\int_{S_k}\phi({{\mathbf{r}}} '){\frac{{{\mathbf{R}}}} {R^3}}\times {\rm d} {{\mathbf{S}}}_k({{\mathbf{r}}} '). $$

(8)

Here, the surfaces between different conductivity regions are denoted by S _k , k = 1...L. The inner and outer conductivities of S _k are represented by σ ⁻ _k and σ ⁺ _k , respectively. R = r−r′ is the vector between the field point r and the source point r′, and $\bar{\sigma}$ is the mean conductivity at the field point. The contribution of the primary sources, g and B ₀, are defined as follows:

$$g({{\mathbf{r}}})=\frac{1}{4\pi\sigma_0} \frac{ {{\mathbf{p}}}\cdot{{\mathbf{R}}}} {R^3}, $$

(9)

$${{\mathbf{B}}}_0({{\mathbf{r}}})= {\frac{\mu_0}{4\pi}}{\frac{{{\mathbf{p}}}\times{{\mathbf{R}}}} {R^3}}, $$

(10)

where σ ₀ is the unit conductivity and μ ₀ is the permeability of the free space. Equations (7) and (8) can be solved numerically by dividing the surfaces into elements and computing the surface integrals over these elements [7–9, 23]. The surface elements used in this study are triangular, quadratic and isoparametric. Using (7) for all nodes and integrating over all elements, a set of equations with N unknowns can be obtained (here N denotes the number of nodes in the BEM mesh). In matrix notation, this can be expressed as:

$$\begin{aligned} &\Phi = {{\mathbf{g}}}+{{\mathbf{C}}}\Phi \\ &({{\mathbf{I}}}-{{\mathbf{C}}})\Phi = {{\mathbf{g}}} \\ &{{\mathbf{A}}}\Phi = {{\mathbf{g}}} \\ \end{aligned}$$

(11)

where Φ is an N × 1 vector of node potentials, C is an N × N matrix whose elements are determined by the geometry and electrical conductivity of the head and g is an N × 1 vector representing the contribution of the primary sources. Φ is then obtained as:

$$ \Phi = {{\mathbf{A}}}^{-1}{{\mathbf{g}}} $$

(12)

To eliminate the singularity in the coefficient matrix A, the method of matrix deflation is employed [33]. IPA is implemented to overcome numerical errors caused by high conductivity difference around the skull layer [29]. Once Φ is computed, B is calculated from the potential values using (8). In matrix notation this can be written as follows:

$${{\mathbf{B}}}={{\mathbf{B}}}_{{\mathbf{0}}}+{{\mathbf{H}}}\Phi $$

(13)

If the number of magnetic sensors is m, B is an m × 1 vector representing the magnetic fields at the sensor locations, and B ₀ denotes the m × 1 vector of magnetic fields at the same sensor locations for an unbounded homogeneous medium. Here, H is an m × N coefficient matrix determined by the geometry and electrical conductivity of the head.