A Matlab library for solving quasi-static volume conduction problems using the boundary element method

A Matlab library for solving quasi-static volume conduction problems using the boundary element method

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 256–263 journal homepage: www.intl.elsevierhealth.com/j...

290KB Sizes 8 Downloads 35 Views

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 256–263

journal homepage: www.intl.elsevierhealth.com/journals/cmpb

A Matlab library for solving quasi-static volume conduction problems using the boundary element method a,b ¨ M. Stenroos a,∗ , V. Mantynen , J. Nenonen c a b c

Helsinki University of Technology, Laboratory of Biomedical Engineering, P.O. Box 2200, FI-02015 TKK, Finland Helsinki University Central Hospital, BioMag Laboratory, P.O. Box 340, FI-00029 HUS, Finland Elekta Neuromag Oy, P.O. Box 68, FI-00511 Helsinki, Finland

a r t i c l e

i n f o

a b s t r a c t

Article history:

The boundary element method (BEM) is commonly used in the modeling of bioelectro-

Received 7 March 2007

magnetic phenomena. The Matlab language is increasingly popular among students and

Received in revised form

researchers, but there is no free, easy-to-use Matlab library for boundary element com-

21 August 2007

putations. We present a hands-on, freely available Matlab BEM source code for solving

Accepted 18 September 2007

bioelectromagnetic volume conduction problems and any (quasi-)static potential problems that obey the Laplace equation. The basic principle of the BEM is presented and discretization of the surface integral

PACS: 02.70.Pt

equation for electric potential is worked through in detail. Contents and design of the library

41.20.Cv

are described, and results of example computations in spherical volume conductors are

41.20.Gz

validated against analytical solutions. Three application examples are also presented. Further information, source code for application examples, and information on obtaining

41.90.+e

the library are available in the WWW-page of the library: http://biomed.tkk.fi/BEM.

87.80.−y

© 2007 Elsevier Ireland Ltd. All rights reserved. Keywords: Bioelectricity Biomagnetism Boundary element method

1.

Introduction

The boundary element method (BEM) [1] is often used in the modeling of bioelectrical or biomagnetic phenomena, such as electro- and magnetocardiography (ECG, MCG) [2–4] and electro- and magnetoencephalography (EEG, MEG) [5]. Typically, BEM computer programs have been programmed with a low-level programming language, e.g., Fortran or C. For these languages, there are many code libraries, from which one can build one’s own BEM solver with some effort. To the authors’



knowledge, there is no complete, freely available C or Fortran library that is suitable for solving bioelectromagnetic volume conduction problems. Currently, the Matlab®1 language is becoming increasingly popular among students and researchers of physics and biomedical engineering. It is not uncommon that a young engineer is fluent in using Matlab, but hardly familiar with C, not even speaking of Fortran. The main advantage – and one reason behind the popularity – of Matlab is the ease of programming: memory allocation is automatic, routines for matrix and vector computations are readily implemented,

Corresponding author. Tel.: +358 94513171; fax: +358 94513182. E-mail addresses: [email protected]fi (M. Stenroos), ¨ [email protected]fi (V. Mantynen), 1 [email protected]fi (J. Nenonen). http://www.mathworks.com/. 0169-2607/$ – see front matter © 2007 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2007.09.004

257

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 256–263

there is a large variety of algorithms and toolboxes, and the source code is platform independent. There are some Matlab toolboxes for EEG/MEG analysis that contain BEM modeling. The BrainStorm2 package contains indepth functions for solving EEG and MEG forward problems with the BEM. Using these functions outside the BrainStorm toolbox or for other kinds of problems is, however, not straightforward: it requires prior understanding on the principles of the BEM, adaptation of the functions for the problem at hand, and thorough validation. The BrainStorm BEM kernel is used, for example, in the FieldTrip toolbox.3 In this paper, we describe the methodology and implementation of a Matlab BEM source code library, which is available from the corresponding author. The library is simple, robust and modular. It is suitable as well for hands-on learning of the BEM as for modeling of various bioelectromagnetic volume conduction phenomena. In addition, the library can also be applied to any scalar surface potential problem that obeys the Laplace or Poisson equation. It is also straightforward to use only one part of the library, e.g., the element integral kernel.

2.

Theory

Bioelectrical volume conduction phenomena are characterized by the Poisson equation: ∇ · (∇) = ∇ · Jp ,

(1)

where  is the electric potential, Jp is the primary current density, and  is the conductivity. Biomagnetic volume conduction problems obey the Biot-Savart law:  = 0 B 4

 V

J(r ) × (r − r ) dV  , |r − r |3

(2)

 is the magnetic induction field, r and are position where B vectors in field and source spaces, and J is the current density containing both primary currents Jp and ohmic volume currents Jv = −∇. Inside a volume conductor, the potential  and normal  are continuous component of the current density (−∇) · n across boundary surfaces. At the outermost boundary of a finite volume conductor, the normal component of the current density is zero. In mathematical treatment it is, in addition, assumed that there are no point-like sources (singularities) at the boundaries. In a piece-wise homogeneous volume conductor with K boundary surfaces separating regions of different conductivities, the Poisson equation can be converted to surface integral form using the Green formulas and boundary conditions. For a field point inside the volume conductor, not on a boundary surface, we get [6]: r

1  k k (− − + ) 4 K

(r)(r)=s ∞ −

k=1

2 3

 (r ) Sk

(r−r ) −→ ·dS , r ∈ / Sk . |r − r |3

http://neuroimage.usc.edu/brainstorm/. http://www.ru.nl/fcdonders/fieldtrip/.

For a field point on a smooth boundary surface, the corresponding integral equation is [7]:

(r)=

k − k 1  − +  − l + l l + l 2 − − + + K

2s

 (r )



k=1

Sk

(r − r ) −→ · dS , |r − r |3

r ∈ Sl , (4)

k and  k are conductivities inside and outside the surwhere − + face k, and ∞ is the potential generated by the sources of the field in infinite, homogeneous medium of conductivity s . In practice, s is interpreted as unit conductivity, used only to get correct units in the ∞ formula. For primary current sources inside a bounded volume conductor, ∞ is

∞ (r) =

1 4s

 V

Jp · (r − r ) dV  . |r − r |3

(5)

This is the most common source scenario in bioelectromagnetic volume conduction problems, covering, e.g., current dipole models and distributed source current models (see the first and second examples in Section 5). In some modeling scenarios obeying Eqs. (3) and (4), ∞ is known a priori, and sources are assumed to lie outside the region of interest. Problems of this kind are commonly referred to as scattering problems (see the third example in Section 5). For the magnetic induction field outside the volume conductor, the integral equation is derived from the Biot-Savart law. Writing the current density in terms of the primary and volume currents and applying the divergence theorem leads to [8]:

 k k =B  0 − 0 B (− − + ) 4 K

k=1



(r − r ) −→ (r ) dS  × , |r − r |3 Sk

(6)

where  0 = 0 B 4

 V

Jp (r ) × (r − r ) dV  |r − r |3

(7)

is the magnetic field in infinite vacuum-like medium due to the primary current density Jp . In order to calculate the effects  we need to know the electric of the volume conductor on B, potential  at each conductivity interface Sk . Thus, to calculate  we first need to solve  from Eq. (4). B,

2.1.

Method of weighted residuals

The boundary element method is based on the method of weighted residuals [1]. Consider the problem: Lf (r) = g(r)

(8)

(3) where L is a linear differential or integral operator. First, approximate f (r) as a sum of N basis functions j (r) with some constant coefficients ϕj , and insert the approximation to the

258

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 256–263

original equation:

f (r) ≈

N 

clk = 2

r) j (

ϕj

(9)

N 

ϕj L j (r) − g(r) = RN (r),

(10)

j=1

where RN is the residual of the approximated solution. Next, the residual is minimized by taking the inner product of the residual and weight functions wi (r), i = 1, . . . , N, over the solution domain ˝, and setting this integral to zero:

 RN (r)wi (r) d˝ = 0



wi (r)L j (r) d˝ =

˝

j=1

g(r ) Sk

(r − r ) −→ · dS , |r − r |3

r ∈ Sl ,

(18)

where in case of two indices in the superscript, the first one labels the field coordinate and the second one the (primed) source coordinate. With this notation, Eq. (4) gets the form

l = bl −

K 

clk Dlk [k ].

(19)

k=1

(12)

l (r) ≈

l

wi (r)g(r) d˝.

N 

˝

Lϕ = g,

(13)

where ϕ and g are N × 1-vector with elements ϕi , gi , and L is an N × N matrix with elements:

 wi (r)L j (r) d˝.

ϕil

l  i ( r ),

(20)

i=1

Now the discretized problem can be written in matrix form:

Lij =



Next, discretize the potential  at each surface into Nl basis functions:



ϕj

1 D [g](r) = 4

(11)

˝

 N 

(17)

l + l − +

lk

j=1



k − k − +

(14)

where il (r) is the i:th basis function at surface Sl and ϕil is the corresponding constant coefficient. Then multiply with the weight function wlj , and integrate in the field space over Sl to get

 Nl  ϕil

i=1

 wlj

Sl

l i

dS = Sl



k

wlj bl

dS −

K N  

c

k=1

lk

m=1

k ϕm

Sl

wlj Dlk [

k m ] dS.

(21)

˝

The coefficients ϕj can, in principle, be solved from Eq. (13) by inverting the matrix L. The most simple choice for the weight function is the Dirac ı function. With this choice, the integrals in previous equations are simplified to evaluations of the integrands at the definition points of the ı functions. This method is called the point collocation method, and the definition points of the ı functions are called the collocation points. Another popular way is the Galerkin method, where same functions serve as basis and weight functions so that the residual is minimized in the whole solution domain [1].

Now we are ready to specify the basis and weight functions. In this work, we use only point collocation weighting with constant or linear basis functions [1]. The most simple boundary element transfer matrix is obtained by choosing the constant basis functions. The basis and weight functions are formed after discretizing each surface Sl to nl nodes and tl triangular elements Til . Constant basis functions and corresponding weight functions are defined as

2.2. Boundary element solution of the integral equations

wlj = ı(r − cjl ),

In the boundary element method, the basis and weight functions are defined on the boundary surfaces. In the BEM implementation chosen for this work, these surfaces are tessellated into triangular elements. Now we discretize Eq. (4) along guidelines presented in the previous section. The pioneer work of this kind of discretization was done for electric ˇ potential by Barnard et al. [9] and for magnetic field by Hora´ cek [10]. First, we introduce the following notations: l

 = (r), l

b =

r ∈ Sl

2s l −

l + +

∞ (r),

(15) l

r ∈ S

(16)

 l  i (r)

=

1, r ∈ Til

(22)

0, r ∈ / Til

(23)

where cjl is the centroid of triangle Tjl . Applying these functions to Eq. (21), we get the equation in a more simple form: k

ϕjl

=

blj



K t  

c

k=1

lk

k lk ϕm ˝jm ,

(24)

m=1

lk is the unit-normalized solid angle spanned by T k where ˝jm m

at cjl :

lk ˝jm

1 = 4

 k Tm

(cjl − r ) |cjl − r |3

−→ · dS  .

(25)

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 256–263

towards zero at nodes j and k [13]:

Now write Eq. (24) in matrix form:

l

l

ϕ =b −

K 

lk

lk k

c  ϕ ,

(26)

k=1

where ϕl and bl have dimension [tl × 1], ϕk has [tk × 1], and lk has [tl × tk ]. Eq. (26) is, however, not enough for solving the potential on Sl : from the right side we see that we have to know the potential on all surfaces. This is done by collecting corresponding equations for all the surfaces together into one matrix-vector equation:

⎞ ⎛ c11 11 c12 12 · · · c1K 1K ⎞ b ⎟ .. ⎜ ϕ2 ⎟ ⎜ b2 ⎟ ⎜ .. ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ . c21 21 . ⎜ ⎟ ⎜ . ⎟ = ⎜ ⎟−⎜ ⎟ ⎜ . ⎟ ⎜ .. ⎟ ⎜ . . ⎟ . .. .. ⎝ . ⎠ ⎝ . ⎠ ⎝ .. ⎠ ⎛

ϕ1





⎛ ×

ϕ1



cK1 K1

⎜ ϕ2 ⎟ ⎜ ⎟ ⎜ . ⎟. ⎜ . ⎟ ⎝ . ⎠

···

···



cKK KK

(27)

(28) (29)

where † labels the inversion of the possibly preconditioned ˜ cannot be directly matrix. In a finite volume conductor, (I + ) inverted, because one of the eigenvalues is zero. Physically this means that the zero level of the potential is undefined. The inversion can be done with standard techniques after defining the zero level. This procedure is called the matrix deflation. It was first presented by Lynn and Timlake [11] and later thoroughly explained and discussed by Fischer et al. [12].

2.3.

=



hlijk (r)

(30)

Tl i{jk}

wlj = ı(r − rjl ),

ri · (rj × rk )

(32)

l . r ∈ / Tijk

0,

With linear basis functions, the discretization procedure is essentially the same as in Eqs. (24)–(28), but the indexing (subscripts, vector and matrix sizes) is based on nodes instead of triangles. The largest change occurs in the element integrals: a linear basis function extends over all the triangles that belong to the node of interest. Thus, Eq. (25) gets the form:

2.4.

 k Tmno

hkmno (r )

(rjl − r ) |rjl − r |3

−→ · dS  .

(33)

(31)

where the sum goes over all the triangles that have node i as a l so vertex, and the shape function hlijk is defined on triangle Tijk that the shape function has value 1 at node i and falls linearly

Calculation of element integrals

All surface integrals over the triangular elements are in this work calculated analytically. The solid angle integrals in Eq. (25) are calculated with the formula derived by van Oosterom and Strackee (Eq. (8) in [15]), and the shape function integrals of Eq. (33) are calculated as presented by de Munck (Eq. (19) in [13]). The element integrals for the magnetic fields are calculated with constant potential approach as specified by de Munck (Eq. (13) in [13]) and with linear potential approach as presented by Ferguson et al. (Eq. (12) in [14]). When the field point r is in the integration domain, the integrals in Eqs. (25) and (33) contain a singularity. In this work, this so called auto-solid angle problem is solved with the method suggested by de Munck [13]: the value of the auto-solid angle is set to such a value that the solid angle spanned by the whole surface is −2. The auto-solid angle is discussed by Ferguson and Stroink [7].

3.

Linear basis functions

In the BEM with constant basis functions, the number of basis functions for a surface of N nodes is the same as the number of triangles, approximately 2N. A more compact transfer matrix is obtained by choosing nl linear basis functions with collocation points in the nodes of the mesh: l  i (r)



Discretization of Eqs. (3) and (6) follows the guidelines presented in this section; the conductivity coefficients and surface integrals change according to the differences in the original integral equations. Discretization of Eq. (6) is thoroughly presented in [14].

Writing all this in brief:

˜ b, ⇒ ϕ = (I + )

⎧   ⎨ r · (rj × rk ) , r ∈ Tl ijk

{no}

ϕK

˜ ⇒ (I + )ϕ ˜ =b ϕ = b − ϕ

hlijk (r) =

lk ˝jm =

1

bK

ϕK

259

Library description

The Matlab library contains the source code of all computational tools needed for boundary element modeling of a (quasi-)static volume conductor problem that obeys any of the integral equations presented in Section 2. In addition, the library can be applied to any static potential problem, which can be characterized with equations similar to Eqs. (3) and (4). The user has to supply the triangle meshes for all boundaries, the volume conductivities, and the discretized primary current density or infinite medium potential. The meshes are described with an array containing the node coordinates and another array containing the triangle descriptions as index pointers to the node array; the format of the triangle description array is the same one that is used by the Matlab mesh processing routines. The library contains functions in following categories:

260

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 256–263

• Mesh preprocessing: Tools for calculating normal vectors, node neighborhoods, etc., for a triangle mesh. • Model building: Tools for discretization of all integral equations described in this paper, including functions from analytical calculation of the element integrals to building of the whole BEM transfer matrix. • Basic examples: Thorough examples of the following problems: ◦ Potential on the surface and inside a homogeneous volume conductor. ◦ Potential on the surface and inside a two-compartment volume conductor. ◦ Magnetic induction field outside a homogeneous volume conductor. ◦ Scripts for performing validation calculations described in Section 4. For validation purposes, these examples are carried out in a spherical mesh with a current dipole source. ◦ Application examples: See more in Section 5. ◦ Visualization: Functions for displaying triangle meshes, scalar data on the meshes, point sets, and vector sets. These tools use Matlab visualization functions; they do not add new features, but are more straightforward to use with the BEM library than the original Matlab routines. ◦ Analytical: Functions for analytical calculation of potential or magnetic induction field generated by a current dipole source: ◦ Potential inside and on the surface of a homogeneous, spherical volume conductor. ◦ Magnetic induction field outside a homogeneous, spherical volume conductor. ◦ Potential inside and on the surface of a concentric, twolayer, spherical volume conductor, when the source is in the centroid of the model. ◦ Potential and magnetic induction field of a dipole in infinite medium. ◦ Miscellaneous: A set of functions for, e.g., error evaluation.

3.1.

Library implementation

In the design and implementation of the library, the key concepts are simplicity, independence on computational platform, and ease of use. Hereby no (pre)compiled code, executables, or linked libraries are used. The whole library is provided as Matlab source code. Matlab built-in functions are used, when feasible. Care has been taken to keep the code readable. Most of the computational cost comes from the calculation of the element integrals. These core functions are vectorized and optimized with aid of the Matlab profiler tool. The example codes are written so that they can be understood, when read along with this article. The example codes thus serve as starting point for the user of the library. Comments contain references to equations in this paper. The core routines are not commented or explained; an interested reader should refer to the original publications by de Munck [13] and Ferguson et al. [14]. The library is developed using Matlab 6.5 (R13) in Intel® Pentium IV PCs equipped with Microsoft® Windows 2000 and XP Professional operating systems. Testing has been

carried out with Matlab versions 6.5 and above in various Windows and Unix systems (Linux, Mac OS 10.8.4). So far only compatibility problems have been related to those visualization functions, which utilize OpenGL rendering. These problems are typically due to the configuration of graphics drivers or operating systems.

3.2.

Distribution and licensing of the library

The Matlab BEM library is available from the corresponding author via the WWW-page http://biomed.tkk.fi/BEM. The source code is open, but the library is not free software.4 For example, the user of the library is not allowed to distribute the library or use the library in a commercial purpose. More information on the library and licensing, including the license text, is provided in the WWW-page http://biomed.tkk.fi/BEM of the library. When the library is used, this article is to be cited.

4.

Validation

The boundary element solver was validated by comparing the results produced with the solver to the results obtained with analytical formulas. The calculations were done in a homogeneous spherical volume conductor. For this geometry, there are analytical solutions in closed form for both the electric  [18] propotential  [16,17] and magnetic induction field B duced by a dipolar source. The sphere was discretized with 162 nodes/320 triangular elements with a Matlab tool written by F. Lindgren.5 The computations were performed with varying source positions at different eccentricities following the ideas presented in [7]: first, the distance of the sources from the centroid of the volume conductor was chosen. Then, 100 random points on this surface were generated. The surface potential generated by radial and tangential unit dipoles in each position was computed both analytically and with the BEM solver. The magnetic field was calculated in 162 random points on an eccentric spherical shell surrounding the volume conductor. In these calculations, the radius of the sphere was the same in both numerical and analytical calculations; no radius scalings like those presented in [7] and [19] were used. The error of the BEM solution was assessed in collocation points with relative error (RE) and correlation (CC) measures: RE =

|ϕa − ϕn | |ϕa |

(34)

CC =

(ϕa − ϕ¯ a ) · (ϕn − ϕ¯ n ) , |ϕa − ϕ¯ a ||ϕn − ϕ¯ n |

(35)

where ϕa is a vector of analytically calculated values of the potential at the collocation points, ϕ¯ a is the mean of the aforementioned potential, and ϕn and ϕ¯ n are their numerical counterparts. In the analytical solution, the zero level was set according to the definition done in the deflation of the BEM transfer matrix: the sum of the potential over all the collocation points was set to zero.

4 5

http://www.fsf.org/licensing/essays/free-sw.html. http://www.maths.lth.se/matstat/staff/finn/prog/dnl/trisphere.m.

261

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 256–263

Fig. 1 – Mean relative error of the potential for (a) radial and (b) tangential dipoles at different source depths for constant and linear basis functions.

Mean relative errors for the surface potential are displayed in Fig. 1, and for the magnetic field in Fig 2. The shapes of the error curves in Fig. 1 resemble closely those published by Ferguson and Stroink in Fig. 3 of [7]. Straight comparison of the error plots is not possible, as Ferguson and Stroink used in their analysis a measure, which scales and offsets the solutions so that the error is minimized. Fig. 11 in [7] shows, however, an example with different error statistics. The “No Fit”-case is equivalent to our RE measure; the error curve has the same shape as our RE curve of the constant collocation method. The absolute value of the error is about twice as large as in our results. This is likely explained by the accuracy of the element mesh: while Ferguson and Stroink used 180 triangles, we used 320 triangles. Scripts for generating the plots shown in Figs. 1 and 2 and corresponding plots for correlation coefficient CC are included in the library. Also included is a code for evaluating the error in case of a dipole in the origin of a spherical, two-layer volume conductor. Examples on computation times for matrix building operations are shown in Table 1. The computations were done

Table 1 – Computation times (in s) for matrix building operations with different mesh sizes and basis functions Number of nodes Basis functions ˜ † (I + ) Constant Linear / S) ϕ (r ∈ Constant Linear B Constant Linear

642

2562

0.3 1.5

5.1 20.1

198.9 327.7

0.1 0.7

0.2 2.5

1.0 9.9

0.5 1.3

2.1 5.0

8.2 19.7



˜ refers to Eq. (29), and matrices ϕ (r ∈ The matrix (I + ) / S) and B to the discretized integrals of Eqs. (3) and (6), respectively. The number / S) and of fieldpoints (equal to the number of rows) in matrices ϕ (r ∈ B is 100.

in a standard desktop computer with Pentium 4 Processor (2.8 GHz), 1 GB memory, and Windows 2000 operating system. All computations were repeated 10 times; the reported numbers are mean values. In construction of the transfer matrix ˜ † , the constant potential approach is significantly faster (I + ) with sparse meshes than the linear potential. The relative time difference between the constant and linear approaches diminishes as the number of nodes in the mesh increases: with dense meshes, the relative time needed for the matrix inversion grows, and the constant potential transfer matrix has about four times as many matrix elements as the linear transfer matrix has.

5.

Fig. 2 – Mean relative error for the magnetic induction field of tangential dipoles at different source depths for constant and linear basis functions.

162

Application examples

In addition to the validation codes, the library contains three application examples. The commented source codes and visualized results of these examples are presented in the WWW-page of the library. In the first example, an application of electrocardiographic inverse problem is presented. The data used in this exam-

262

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 256–263

ple was created with the ECGSIM program6 [20]. A region of transmural myocardial ischemia in the region of left anterior descending coronary artery was simulated by reducing the amplitude of the action potential locally. This source configuration produced a body surface potential map with elevation of the ST-segment potential on the left anterior thorax (see, e.g. [21]). The thorax was modeled as a homogeneous volume conductor using the ECGSIM thorax mesh. Cardiac sources were approximated with an equivalent model containing normal current density on endo- and epicardial surfaces of the ECGSIM heart mesh. A lead field matrix between the current sources and the body surface potential was set up and inverted with aid of truncated singular value decomposition regularization (the Matlab function pinv). The spatial distribution of the estimated current densities was similar to the reference potential on the heart surfaces; the region of elevated epicardial potential was localized correctly. The second example deals with magnetoencephalographic forward problem. An accurate mesh (1995 nodes) describing the outer surface of the brain was created with the SOVITA software [22]. The mesh was co-registered with the field reconstruction surface of the Neuromag-122 MEG helmet [23]. A linear boundary element model for the magnetic induction field on the reconstruction surface was set up, and fields generated by various current dipoles in the cortex were simulated. According to our visual inspection, the results were very similar to those obtained with the commercial Neuromag source modeling software. In the third example, we present a simple electrostatic scattering problem. A dielectric sphere is brought into an initially homogenous electric field, where it causes dipolar scattering field due to the electric polarization. The problem obeys Eqs. (3) and (4), where conductivities are only replaced with permittivities, and the infinite medium potential is replaced by the electric potential of the initial electric field. A linear boundary element model was set up, and potential was calculated both inside, outside, and on the surface of the scatterer. The numerical results were compared to the analytically calculated ones (see, e.g. [24]); relative error was everywhere below 1%.

6.

Discussion

We have described the basic theory behind boundary element modeling of the (quasi-)static volume conduction problems, and contents of our Matlab BEM source code library. The basic example computations have been validated against analytical calculations. The validation results of the electric potential are similar to those published by Ferguson and Stroink [7]. The computation times on a standard PC are relatively short. If there is no need for very dense meshes (say, over 4000 nodes), the performance of Matlab and our library on a standard PC should be sufficient for most applications in the field of bioelectromagnetism. With more dense meshes, the memory capacity needed in the matrix building and inversion is the limiting

factor (Matlab performs all computations in double precision). Our library uses only point collocation weighting. In some studies, however, the Galerkin method is preferred to the collocation method [5,19], especially in the modeling of MEG. The Galerkin method is straightforward to implement on basis of this BEM library. Essentially one has to evaluate all the integral operations over the whole surfaces instead of the collocation points. These outer integrals (see Eq. (21)) are typically evaluated numerically. This will increase the computational costs, especially in case of the linear potential approach: in a mesh with N nodes and 2N triangles and an M-point quadrature, one needs to evaluate the linear shape functions h (Eq. (32)) in 2NM points instead of N points needed with the collocation method. In calculation of the electric potential, the error over the surface is reported to be slightly smaller in the Galerkin solution than in the collocation solution (Figs. 1 and 4 in [19]). The relative performance of different weighting methods depends, however, on the chosen error evaluation technique: error evaluation in the mesh nodes favors the collocation method [19], while error evaluation over the whole surface favors the Galerkin method. And, using the mesh nodes as error evaluation points might be a more natural approach, if the electrodes are located at the nodes of the mesh (see discussion in [5]). Herewith one may ask, whether benefits gained with the Galerkin method are worth the computational cost (see Fig. 9 in [19]). The library contains all the essential computational routines needed for BEM modeling of the bioelectromagnetic volume conduction phenomena formulated in terms of primary currents. Applications defined without modeling of the primary currents or the infinite medium potential, e.g., the epicardial potential problem [25], cannot be computed with the current version of the library. If the reader would like to perform epicardial potential modeling with the library, he may contact the corresponding author for further information and possible auxiliary functions. Tools for modeling of the anatomy are not included in the library. Segmentation of anatomical images, triangulation of surface meshes, and image registration are reviewed by ¨ onen ¨ Lotj [22], who also provides tools for constructing surface models from magnetic resonance images. Other free and commercial mesh generation tools are listed by Schneiders.7 This paper and the library provide both a hands-on tutorial for learning of the BEM and a flexible toolkit for solving bioelectromagnetic volume conduction problems. The source code is freely available for academic use from the corresponding author via the WWW-page http://biomed.tkk.fi/BEM.

Acknowledgements ¨ onen ¨ The authors thank Jyrki Lotj (VTT, Tampere, Finland) for the surface mesh of the brain. M. Stenroos thanks The ¨ Foundation of Technology in Finland and V. Mantynen thanks

7

6

http://www.ecgsim.org.

http://www-users.informatik.rwth-aachen.de/∼roberts/ software.html.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 256–263

Instrumentarium Science Foundation and Alfred Kordelin Foundation for financial support. [13]

references [14] [1] C. Brebbia, J. Dominguez, Boundary elements: An Introductory Course, McGraw-Hill, 1989. [2] J. Nenonen, Solving the inverse problem in magnetocardiography, IEEE Eng. Med. Biol. 13 (1994) 487–496. [3] R. MacLeod, D. Brooks, Recent progress in inverse problems in electrocardiology, IEEE Eng. Med. Biol. 17 (1998) 73–83. [4] R. Gulrajani, The forward and inverse problems of electrocardiography, IEEE Eng. Med. Biol. 17 (1998) 84–101. [5] J. Mosher, R. Leahy, P. Lewis, EEG and MEG: forward solutions for inverse methods, IEEE Trans. Biomed. Eng. 46 (1999) 245–259. [6] D. Geselowitz, On bioelectric potentials in an inhomogeneous volume conductor, Biophys. J. 7 (1967) 1–11. [7] A. Ferguson, G. Stroink, Factors affecting the accuracy of the boundary element method in the forward problem. I. Calculating the surface potentials, IEEE Trans. Biomed. Eng. 44 (1997) 1139–1155. [8] D. Geselowitz, On the magnetic field generated outside an inhomogeneous volume conductor by internal current sources, IEEE Trans. Magn. MAG-6 (1970) 346–347. [9] A. Barnard, I. Duck, M. Lynn, W. Timlake, The application of electromagnetic theory in electrocardiology. II. Numerical solution of the integral equations, Biophys. J. 7 (1967) 463–490. ˇ [10] B. Hora´ cek, Digital model for studies in magnetocardiography, IEEE Trans. Magn. MAG-9 (1973) 440–444. [11] M. Lynn, W. Timlake, The use of multiple deflations in the numerical solution of singular systems of equations, with applications to potential theory, SIAM J. Numer. Anal. 5 (1968) 303–322. [12] G. Fischer, B. Tilg, F. Hanser, B. Messnarz, P. Wach, On modeling the Wilson terminal in the boundary and finite

[15] [16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24] [25]

263

element method, IEEE Trans. Biomed. Eng. 49 (2002) 217–224. J. de Munck, A linear discretization of the volume conductor boundary integral equation using analytically integrated elements, IEEE Trans. Biomed. Eng. 39 (1992) 986–989. A. Ferguson, X. Chang, G. Stroink, A complete linear discretization for calculating the magnetic field using the boundary element method, IEEE Trans. Biomed. Eng. 41 (1994) 455–460. A. van Oosterom, J. Strackee, The solid angle of a plane triangle, IEEE Trans. Biomed. Eng. 30 (1983) 125–126. F. Wilson, R. Bayley, The electric field of an eccentric dipole in a homogeneous spherical conducting medium, Circulation 1 (1950) 84–92. D. Yao, Electric potential produced by a dipole in a homogeneous conducting sphere, IEEE Trans. Biomed. Eng. 47 (2000) 964–966. J. Sarvas, Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem, Phys. Med. Biol. 32 (1987) 11–22. S. Tissari, J. Rahola, Error analysis of a Galerkin method to solve the forward problem in MEG using the boundary element method, Comput. Methods Prog. Biomed. 72 (2003) 209–222. A. van Oosterom, T. Oostendorp, ECGSIM: an interactive tool for studying the genesis of QRST waveforms, Heart 90 (2004) 165–168. ˇ B.M. Hora´ cek, G.S. Wagner, Electrocardiographic ST-segment changes during acute myocardial ischemia, Cardiac Electrophysiol. Rev. 6 (2002) 196–203. ¨ onen, ¨ J. Lotj Construction of patient-specific surface models from MR images: application to bioelectromagnetism, Comput. Methods Prog. Biomed. 72 (2003) 167–178. ¨ al ¨ ainen, ¨ J. Knuutila, A. Ahonen, M. Ham M. Kajola, P. Laine, O. Lounasmaa, L. Parkkonen, J. Simola, C. Tesche, A 122-channel whole-cortex SQUID system for measuring the brain’s magnetic fields, IEEE Trans. Magn. 29 (1993) 3315–3320. J. Vanderlinde, Classical Electromagnetic Theory, Kluwer Academic, 2004, pp. 185–187. ˇ B. Hora´ cek, J. Clements, The inverse problem of electrocardiography: A solution in terms of single- and double-layer sources on the epicardial surface, Math. Biosci. 144 (1997) 119–154.