IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. 8, AUGUST 1992

805

Equivalent Source Estimation Based on the Calculation of the Electric Field from Depth EEG Data Louis Lemieux and Andrd Leduc

Abstract-This paper describes a method for the estimation of equivalent source position and strength based on the estimation of the electric field from depth electroencephalographic (EEG) data. The calculation method for the electric field is based on a tetrahedral geometry. The proposed approach for source parametrization is twofold. Firstly, the distribution of electrical energy by the squared norm of the electric field vector can give an estimate of the source position, without having to assume a dipole source. Secondly, the average electric field can be related to the dipole magnitude and orientation of the equivalent source. Simulation results demonstrate the potential usefulness of the method. The effects of noise and sampling, and the geometry of the measurement system (i.e., implanted electrodes) relative to the source are also investigated through simulations.

I. INTRODUCTION TEREOTAXIC implantation of depth electrodes is done on a regular basis as part of the clinical investigation of patients who are potential candidates for surgical removal of the area of the epileptic focus [l]. Many of the quantitative methods currently applied on the data to help in the localization of the source(s) are related to signal analysis: spectral analysis, pattern recognition, and coherence analysis [ 11-[3]. Although the parametrization of dipole equivalent sources of electrical activity in the brain based on application of physical principles has been the topic of several publications over the years [4]-[6], those based on the physical analysis of depth EEG are scarce. Church [7], [8] emphasized the differences with the problem of surface EEG and proposed a method for some localization which is an adaptation to depth measurements of the single moving dipole approach (SMD). It has been demonstrated that the single point dipole is a valuable model for dipolar disk and anular sources, at least from the point of view of the localization based on scalp measurements [9]. However, it can be regarded as a limitation on the general applicability to neuronal generators of arbitrary nature, particularly in the case of epilepsy where the source of activity can be diffused and/or can spread rapidly. Furthermore, in the case

S

Manuscript received June 9, 1989; revised November 29, 1990. The authors are with the Department of Neuroelectronics, Montreal Neurological Institute, Montreal, P.Q., Canada H3A 2B4. IEEE Log Number 9201483.

of depth measurements, the relative importance of the higher multipolar moments is increased due to the proximity of the measuring electrodes to the source. A different approach of field characterization, related to the one we propose, is current source density (CSD) analysis. It is based on the calculation of the Laplacian of the potential to estimate the source current density [lo], [ll]. Estimation of the Laplacian requires the potential measurements to be done on a regular spatial array; therefore three-dimensional applications of CSD analysis necessitate the use of special electrodes. Up to now, applications of the CSD method are limited to studies at the microscopic level. The purpose of this paper is to present a novel approach to the problem of parametrizing the equivalent source of focal epileptic activity and neuronal generators in general. The theoretical basis is outlined, and a simulation study is presented. The proposed method is based on the calculation of the electric field in the brain from referential depth EEG, as acquired in ordinary clinical practice. It is possible to estimate the dipole moment of a source by calculating the average electric field over a certain volume. Furthermore, the calculation of the distribution of energy based on the estimated electric field can be used to estimate the position of a source. The initial motivation for this work was the search for a method of calculating source parameters with a minimum of assumptions, based on the application of physical principles; the obtained values would then be used as initial guesses for Church’s SMD method. The proposed method of source parametrization is independent of a dipole hypothesis. The only assumptions made are that the source be of relatively limited extent and located in the vicinity of the implanted electrodes. The latter is often guaranteed in clinical practice where noninvasive diagnostic techniques are used to help obtain an initial estimate of the position of the epileptic focus, and therefore to direct the implantation of depth electrodes in the target area. In order to make a clear demonstration of the validity of the approach we restrict ourselves to the problem of parametrizing a unique dominant source located in the immediate vicinity of the implanted electrodes. Through simulations we demonstrate the relevance of the calcula-

0018-9294/92$03.00 0 1992 lEEE

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. 8, AUGUST 1992

806

tions for ideal and realistic electrode configurations, and measure the effects of noise and spatial sampling. 11. PRINCIPLES The proposed method is based on the estimation of the electric field in the brain from sets of four voltage values measured with depth electrodes. The calculation of the electric field in the brain is useful in two ways for the parametrization of generators located near the implanted electrodes. Firstly, the average electric field over a spherical volume is proportional to the dipole moment of the sources located inside the volume. We suggest, therefore, that this principle be used to estimate the dipole moment of an equivalent source. Secondly, the squared magnitude of the electric field vector is proportional to the electrical energy dissipated per unit volume and unit time. The energy distribution reflects the position of the sources. We propose, therefore, that the position of a dominant source be estimated based on the distribution of energy in the brain.

Fig. 1 . Partitioning of the volumes of exploration. Example for only three implanted strands with six electrodes on each. The most superficial tetrahedron is represented by the wide solid line. Two other tetrahedra are represented by dashed lioes. In this case the total number of tetrahedra would be 15.

-/-----

A. Calculation of the Electric Field in the Brain Depth EEG data is obtained from stereotaxically implanted strands, each one comprising a number of evenly 1 mm). We define the spaced electrodes (dimensions volume of exploration to be the imaginary volume defined by the positions of the N electrodes,' i.e., the polyhedron enclosed within the imaginary surface consisting of triangular faces obtained by drawing lines between the outer electrodes (Fig. 1). Coordinates are defined according to Fig. 2. The proposed method for the estimation of source parameters is described in detail in Section 11-B and relies on the calculation of integrals of the electric field over the volume of exploration. In order to integrate properly one must partition the volume of integ,ration. This can be accomplished by dividing the volume into tetrahedra, the vertices of which are four electrodes. Fig. 3 illustrates such a partitioning for the case of only three strands. First, a preliminary partitioning of the sagittal (XoZ)projection plane is performed, on which the strands appear as points forming the vertices of triangles. For each triangle there is a corresponding volume which projects into the brain, stretching as far as the deepest electrode on each strand. These volumes can in turn be thought of as a pile of tetrahedra, the vertices of which are two electrodes along one of the strands and one on each of the other two. The s u m of these tetrahedra provides the desired partition. If strands are implanted on both sides of the head, the procedure must be performed for each side, and a proper partitioning must also be done for the (medial) volume which lies between the two sets of tetrahedra. Note that any partition done this way is not unique; it depends on the arbitrary choice of triangles (for cases with more than three

-

'In this paper the term electrode is used to designate an active, or measuring electrode.

X

J h

Y Fig. 2. Definition of head coordinates and fhe dipole vector, p', coordinates: 0' s 8 s 180", 0" s Q s 360": X i s the position vector of the dipole. The origin of the X n coordinate system is at the center of the spherical head model, the profile of which is outlined, showing the nose at the left. The hatted values represent the unit coordinate vectors.

strands), and on the (also arbitrary) choice of the first tetrahedron for each "pile." Therefore, for a given set of electrodes, there are a number of possible partitions. The number of tetrahedra in a given partition depends on the number of electrodes, and the number and relative position of the strands. In order to be able to apply the proposed method one must have at least three implanted strands each of which must not all lie in the same plane. The partitioning of the volume of exploration into tetrahedra has many practical advantages. It is flexible, as it adapts to any configuration of stereotaxically implanted electrodes; it is adaptable to subdural electrodes; it is simple since the tetrahedra are directly defined by the positions of the electrodes, and their volume may be calcu-

LEMIEUX AND LEDUC: SOURCE ESTIMATION BASED ON ELECTRIC FIELD CALCULATION

K

F LJ

I

J

Fig. 3. Example of possible electric field edge vectors in a tetrahedron for a given origin. The tetrahedron is defined by the points I , J , K, and L , which corresponds to electrode positions. Here the electric field points towards the inside of the tetrahedron and none of the components are null.

lated directly using the following formula: vTET,

=

i[AI

(1)

where [A] is the determinant of the matrix formed by the coordinates of the vertices of the tetrahedron:

In addition, the tetrahedron provides a simple device for calculating the electric field, Average Tetrahedron Field: The tetrahedron is the building block of the calculation since it is from the voltages at the four vertices that we estimate the electric field in the vicinity. Starting with the four voltages measured at the vertices of a tetrahedron, , +4 } , it+is possible to determine the average tetrahedron j e l d , ETET,. For a given vertex, Z, taken as an “origin,” the three aL- + I , give differences of potential, aJ - +f, ipK an estimate of the average electric field projected along three nonorthogonal axes (see Fig. 3)*:

---

+,,

+I

-

+j

=

+ -ETETm

(SZ, - q)

j = J,

K,L. (3)

Zi

The vectors in (3) represent the position ve5ors of the vertices. Solving this system of equations for & E T U , we obtain the orthogonal components of the average tetrahedron field. *Since the field is static and therefore conservative, the selection of any vertex yields the same result.

807

The average tetrahedron field taken as a set of vectors in space can be regarded as a 3-D generalization of bipolar EEG. It takes explicitly into account the relative positions of the electrodes, which is reflected in the vectorial nature of the data. The average tetrahedron field vector is assigned to the geometric center of the tetrahedron. It is an estimate of the average electric field inside the tetrahedron. The accuracy of this estimate depends on the following factors [12]: the degree of spatial fluctuation (or uniformity) of the field in the vicinity of the tetrahedron; the dimensions of the tetrahedron; and the geometry (for example: compact versus elongated) of the tetrahedron. The last two factors correspond to the efective uniformity of the field across the tetrahedron. Intuitively, we expect the error on the estimated field to decrease with tetrahedron size as the measurement becomes more ‘‘punctual.’’ As demonstrated in Fig. 4 the tetrahedron has the effect of smoothing, or averaging the field, therefore causing the estimated field to be attenuated for nearby sources. Furthermore, a more compact tetrahedron gives a better estimate of the local field since the field will be effectively more uniform inside the tetrahedron. Although we will refer to these observations throughout this paper, we will not discuss this topic further. We are interested only in “global” quantities which are derived from these local estimates such as the center of energy and the average field.

B. Estimation of Source Parameters Dipole Moment-Average Electric Field: In order to estimate the dipole moment of a source we propose the use of the following result: the average electric field over a spherical volume in an infinite homogeneous conducting medium due to an arbitrary charge distribution enclosed in the volume is proportional to the total dipole moment of the charge di~tribution.~ This is given by * (E ) = --p‘ (4) 3av where U is the medium conductivity; V, the volume of the sphere; p’, the dipole moment, in Am, defined relative to a coordinate system who’zorigin is located at the center of the sphere [13], and ( E ) is the average electric field inside volume V given by the integral:

Fig. 2 shows the angles and vector components which define the orientation of p’. Note that the dipole moment is independent of the choice of origin for sources with zero total current (i.e., the total source current is equal to the total sink current), such as the ones considered here. In order to use (4) to calculate the dipole moment of a ’This result, derived in electrostatics, can be applied to the case of quasistationary fields in a conducting medium by substituting U for E and the source current density for the density of free charge [ 111, [ 141.

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. 8, AUGUST 1992

808

The estimated average electric field, (2 ) c, is therefore given by the finite version of (5):

- 1 porallel

1021

3-\

‘A

( E ) c = -V1 cm Emurn

perpendiculor

where V is the volume of exploration, Emis the average tetrahedron field in tetrahedron m and U,,, is the volume of the tetrahedron. The dipole moment, p’c, is obtained by inverting (4) :

---

jic = - 3 a V ( Z ) c

-i

lo 10 - 3 ‘ 10

-’

(7)

Source Position-Center of Energy: For a linear conducting medium, such as brain tissue in the macroscopic approximation, the current density is proportional to the electric field (Ohm’s law) [14]): +

--t

J

=

aE.

(8)

I

1 DISTANCE (cm)

10

Fig. 4. Illustration of the effeet of the tetrahedron on the estimated field for a pure dipole. The average tetrahedronfield (symbols) and the electric field calculated at the center of the tetrahedron (lines: 1 / r 3 )are plotted as a function of the distance between the source and the center of the tetrahedron for two dipole orientations: parallel (continuous line): dipole pointing towards the tetrahedron; perpendicular (dashed line): dipole pointing perpendicular to the line joining the dipole and the tetrahedron. We notice the effect of averaging, or smoothing, of the field which becomes more important when the source is near and consequently the field more inhomogeneous. The tetrahedron is a perfect pyramid with edges of 1 cm (volume = 0.118 em3). The values at 0.18 cm are for dipoles located inside the tetrahedron.

dominant source, we must assume it to be located inside the volume of exploration. Furthermore, by dominant source, we mean that the field in V must be dominated by the activity of that source. This is required since the average field in Vis the sum of fields due to all the sources in the medium, both inside and outside V, while (4) only accounts for sources located in V. This condition will be discussed further in the section on simulations. In the case of the neuronal generators, the medium is the head, which is neither infinite nor homogeneous. It has nonetheless been shown that for depth EEG measurements, the brain can be approximated as homogeneous [15]. Moreover, the finite dimensions of the brain when modeled as a sphere has relatively little influence on the field in the vicinity of the source. These factors will influence the actual fields in the brain, and in general may cause deviations from (4). Furthermore, in the application to depth EEG measurements, the volume of integration (i.e., the volume of exploration) will not be spherical: the calculation of the average electric field by ( 5 ) will be done on voltage data measured with implanted electrodes, and the configuration of their arrangement, i.e., the volume they enclosed will not be spherical. However, the results of the simulations presented in a later section support the validity of the application of the above principle for finite conducting media and nonspherical volumes of integration.

The volume density of dissipated energy (per unit time), w, is given by [ 161

w = a1E12.

(9)

We call this quantity the density of electrical energy, or simply the energy density. As will be demonstrated in this paper, the energy density is a valuable tool for the estimation of the position of a unique source. The electroneutrality of the brain means that sources are really source-sink pairs or more complex arrangements of sources and sinks [ 111. For simple sources such as dipoles, quadrupoles, etc., the electric field intensity and therefore the energy density are maximum at the source and decrease with distance. For a point dipole the energy density falls off as the distance to the sixth power: for example, from one to two centimeters the energy density decreases by a factor of 64. For a dipole layer, the fall off rate of the energy density is slower than for point dipoles in the direction perpendicularto the layer near the source [ l l ] , [17]. For example, for a disk-shaped layer of thickness 3 mm and radius 1 cm, the energy density decreases by a factor of 16 between 1 and 2 ~ m We . ~ therefore expect the energy distribution to have a “sharp” maximum in either case. Furthermore, in all the above cases, the energy distribution has spatial symmetry about the position of the source. For the point dipole, the distribution has axial symmetry about the axis of the dipole and reflective symmetry about the plane perpendicular to the axis. The above arguments suggest the use of the center of the energy distribution as a mea+sure of the position of the source. The center ofenergy, XCE,is given by

wdV 4This rate of fall gradually increases with distance so as to gradually become equal to the case of a point dipole.

I -

LEMIEUX AND LEDUC: SOURCE ESTIMATION BASED ON ELECTRIC FIELD CALCULATION

809

where V is the volume of exploration. The center of en- havior of the method in terms of sampling, noise level, ergy is the energetic analog of the center of mass: it is the and geometrical factors. We are therefore interested in first moment of the energy distribution, the weight of a calculating the systematic errors on the source parameters given position being the estimated energy density at that resulting from: the finiteness of the measuring device (in position. Taking the case of a dipole source in an infinite terms of the number of electrodes), the nonsphericity of medium, it is easy to demonstrate that the position of the the volume of exploration, the inherent physiological center of energy is equal to the position of the source. For noise, and the relative configuration of the source/immore complex sources such as the ones described above planted electrodedhead. the center of energy will be pulled towards the regions of For most of the simulations the source is located inside greater energy density. the volume of exploration. As mentioned in the Section It is important to note that this method of calculating 11, the method we propose for the estimation of the source the position of a source is accurate only if the source is parameters is based on this assumption. The important located within the volume of exploration (or very close to problem of determining whether a source is effectively loit) because the possible results of (10) are confined to V . cated inside the volume of exploration is not addressed in As will be demonstrated in the simulations, the center of this paper. In order to achieve the above goals, two types of simenergy can nevertheless be useful in localizing sources ulations were performed: located outside the volume of exploration. As mentioned in Section 11-A, the average tetrahedron Distributed Dipoles: Here we consider the average refield gives an estimate of the average electric field in a sults of a number of simulations, N (generally loo), under tetrahedron. Therefore, the squared magnitude of the av- identical conditions, but for dipoles randomly distributed erage tetrahedron field is an estimate of the average en- in an imaginary volume, called a volume of distribution. ergy density inside a tetrahedron w. A finite version of Dipole position is randomly generated within a given vol(10) is given by ume of distribution. We then repeat this process for N 1 other dipoles with the same dipole moment. The volume WmZm m of distributions can be right parallelepipeds or spherical xc, = (1 1) shells (which include spheres). Source parameters are calWln culated for each case. For most simulations the volume of distributions are totally enclosed within the volume of exwhere the index m runs from 1 to nT, the total number of ploration, in which case they correspond approximately tetrahedra, and where Zm is the position of the geometric to the effective volumes of integration. The choice of the center of tetrahedron m. Note that in (10) we are not tak- dipole moment depends on the volume symmetry; for exing into account the volumes of the tetrahedra. Therefore, ample, for a cubic volume of exploration, it is not necthe weight of a given position is the energy density at that essary to simulate dipoles with dipole moments (1, 0, 0) position only. Adding the volume factor would merely and (0, 1, 0) because the fields will be the same relative give more weight to the larger tetrahedra, for which the to the volume (unless the volume of distribution is not error on the field is largar, as explained in Section 11-A. symmetric with regard to the volume of exploration, for Note that (11) can take into account a nonuniform con- example an off centered sphere). The statistics of the reductivity by using an appropriate (T in each tetrahedron. sults give us an estimate of the performance of the method Effective Volume of Exploration: It is easy to see that for dipoles located inside the volume of distribution, not (11) cannot give a result which is outside the volume de- just one arbitrarily located dipole. The shell volumes alfined by the positions of the centers of the tetrahedra, low us to study the effect of eccentricity relative to the which is smaller than the volume of exploration. We will volume of exploration. call this volume the effective volume of exploration. In Series: For a given dipole in a noisy environment, a practical terms this means that, for example, in a situation series of 100 simulations are performed, in which the where only three stereotaxically implanted strands are varying factor is the randomly generated (“instrumenavailable on one side of the head the spatial resolution in tal”) noise, and in some cases the parameters of the dithe plane perpendicular to the strands is null (and the efpole, such as the magnitude of the dipole moment. Each fective volume is null) since the centers of the tetrahedra such individual simulation is called an “instant” in a sethen all have approximately the same coordinates in that ries. The results are presented as a graphical series of the plane. One could nonetheless still have spatial resolution estimated values of the parameters of interest. This type along the axis parallel to the strands, as will be demonof simulation is used to illustrate the relative amplitudes strated in the simulations. of the fluctuations of estimated parameters for a specific source and configuration. 111. SIMULATIONS The purposes of the simulation are: firstly, to demonstrate the applicability of the principles; secondly, to A. Model of Head, Source and Noise measure the performance of the numerical methods under We use Wilson’s expression for the electrical potential realistic conditions; thirdly, to formally evaluate the be- due to an ideal dipole in a spherical homogeneous me-t

~

c c

-

810

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. 8, AUGUST 1992

dium, with radius = 7.5 cm [ l l ] . The values represent Difference in Distance from Center of Volume of Expoint electrode voltage measurements. ploration (DDV, in cm): The difference between the esInstrumental Noise: Noise is simulated by adding ran- timated distance from the center of the volume of explodomly generated numbers, U , directly to the data (uncor- ration and the real value. A negative value means that the related local noise). This noise is used to represent both estimated source is closer to the center of the volume of physiological noise, that is the electrical potentials due to exploration than the simulated source. the normal activity of the brain, and instrumental noise, Difference in Distance from Center of the Head (DDH, which includes electrode polarization, digitizing error, the in cm): The differences between the estimated radius relerror on the estimation of the position of the electrodes ative to the center of the head and the real value. This (i.e., electrode separation) as reflected in the calculated measure is used in the evaluation of the ability of the field^.^ method to locate superficial (i.e., near the scull) sources. The signal to noise ratio is defined as the ratio of the Error on Dipole Orientation (EDO, in Degrees): The maximum voltage measured to the amplitude of the noise, angle between the estimated dipole vector and the simuA [ - A IY IA ] , lated dipole vector:

and the relative noise level A,:

ED0 must be thought of as defining a cone around the simulated dipole vector gs.The range of E D 0 is (O”, 180” }: values larger than 90” mean that the cone points in a direction opposite to Zs. Relative Error on Dipole Magnitude (EDM, Relative Value):

The data for simulations with noise were generated by the following method: first, voltage values are calculated using Wilson’s formula, then the maximum voltage on the measuring array (without noise), is found and the EDM = I F S 1 - I p ’ c l I $SI amplitude of the pseudorandom noise to be added to the voltages, A , is set according to (12), for the desired sig- where I $ 1 is the magnitude of vector p’. The range of nal-to-noise ratio. In this way we obtain distributed di- values of EDM is { -00, 1.0}, corresponding to p c >> poles with constant SINmax.In the case of the series sim- ps and p s >> p c , respectively. EDM can also be repreulations, the amplitude A is set constant and for each sented as a percentage. It is a valid measure of the error simulation, or “instant,” randomly generated noise is on dipole magnitude only, given that E D 0 is small. added to the voltage values (which are constant, unless The above measures of error are the quantities of which the parameters of the source of the noise level are varied). averages and standard deviations are calculated in the disPhysiological Noise: Physiological noise can be simu- tributed dipole simulations. Note that in the case of EDO, lated by adding the fields due to a number of “noisy” arbitrarily oriented estimated dipoles would give an avdipoles N,, to the source field [ 6 ] .The noisy dipoles have erage of 90”, with a rather large standard deviation; this equal amplitude, p n , but random orientation, and are ran- must not be confused with the situation where the estidomly distributed in the whole spherical head model. The mated dipole is systematically at right angle with the real signal-to-noise ratio is then defined as the ratio of the dipole, for which the average error would still be 90” but source magnitude to the total magnitude of the noisy di- with a small standard deviation. poles, NnPn,

B. Measures of Error on Source Parameters Five measures of the error on the estimated source parameters will be used. These are: Absolute Error on the Position (DX, in cm):

where gcE is the estimated position of the source and is the position of the simulated source.

gs

’The error on the position of the tetrahedra due to the uncertainty on the positions of the electrodes is not taken into consideration. Note that in clinical practice it is possible to measure electrode positions within 1 mm.

IV. RESULTS A. Infuence of the Geometry of the Volume of Exploration In light of the fact that (7) is only exactly true for a spherical volume of integration, the effect of the geometry of the volume of exploration was measured by applying the method to three volumes: (pseudo-)sphere6: concentric with the spherical medium, V = 327 cm3; (pseudo-)half-sphere: upper half of the pseudosphere, V = 167 cm3; cube: V = 216 cm3. 6Pseudo-sphere in the sense that is represents a sphere to the extent possible with a finite number of points.

81 1

LEMIEUX AND LEDUC: SOURCE ESTIMATION BASED ON ELECTRIC FIELD CALCULATION

In all three cases 35 electrodes were used. In order to cancel any sampling effect, the dimensions as well as the position of the strands of the three volumes were chosen such that the average tetrahedral volume (the total volume of exploration divided by the number of tetrahedra in the partition), (U), and the geometry of the individual tetrahedra were nearly identical for all three cases. The volume of distributions were identical, that is, 100 simulations were done for dipoles randomly distributed inside a sphere of 1.8 cm radius, the center of which was 2 cm above the center of the head. The results are summarized in Table I. First we note that the results for the dipole moment for the spherical volume of exploration show good the calculated and the simulated dipole vectors. These results are in effect a validation of the calculation. Coming back to the effect of the geometry of the volume of exploration, we note that the relative error on dipole magnitude is larger for the cube and the half-sphere, which is to be expected since our estimation of the dipole moment is based on a spherical volume of integration. We note that in the case of dipole orientation, the error is still relatively small even for the case of the half-sphere. As to the error on the position, it is nearly the same for all three volumes. The relationship between the geometry of the field and that of the volume of exploration is evident in the case of the half-sphere for the calculation of the dipole orientation: for p’ = (0, 0, l ) , the field’s geometry “matches” that of the volume of exploration (which has rotation symmetry in the XoY plane) in such a way that the error on dipole orientation is small. Also, the dipole magnitude is underestimated, which reflects the fact that part of the field is “missing.” For the cubic volume, the dipole magnitude is overestimated. Overall, we can say that even for nonspherical volumes of integration such as the halfsphere, (7) gives a fairly accurate estimate of the dipole moment of a source located within the volume of exploration. This is due to the fact that the magnitude of the dipole field decreases rapidly with distance, 1 / r 3 , and therefore, the effect of the extra (or missing) contributions for a nonspherical volume is small.

B. Injuence of Sampling Sampling refers to the total number of electrodes and the tetrahedral volume, which is equivalent to the reciprocal of the volume density of electrodes. In the following simulations, we study the influence of the above factors for various volumes of integration on the estimation of the source parameters. For most simulations the total volume defined by the outer electrodes is kept constant for a given volume of exploration, while the number of electrodes inside the volume is varied, thus obtaining different tetrahedral volumes. Volumes of Integration and Distribution We have defined the volume of exploration, which corresponds to the volume of integration in (7), to be the

I

TABLE I EFFECT OF THE GEOMETRY OF THE VOLUME OF INTEGRATION ON THE ERRORS ON DIPOLE PARAMETERS Error on Dipole Parameters

Dipole Vector PY

Volume

p,

(A cm)

pr

DX (cm)

ED0 (deg)

EDM

2.8 3.1 4.0

-0.04 -0.03 -0.01

sphere sphere sphere

1

0

0

0

1

1

0 1

1

0.96 1.06 0.87

1/2-sphere 1/2-sphere Ij2-sphere

1 0 1

0 0 1

0 1 1

0.89 0.85 0.77

12.5 2.5 24

0.3 0.5

0 0

0 1

0.84 0.92

1

1

0.17

4.9 4.1 4.7

-0.2 -0.3 -0.2

cube cube cube

1 0 1

0.0

No simulated noise. Number of electrodes: N = 35. The average tetrahedron volumes are equal for the three volumes (( U ) = 2.8 cm3). Average results for 1 0 0 simulations. rel.: relative units.

volume corresponding to the configuration of the electrodes. In practical situations the geometry of this volume will vary from case to case. Four types of volumes of exploration were used in the simulations: I. Cubic (512 cm3): 25 strands arranged in a square array, with maximum of 5, evenly spaced electrodes on each probe for a maximum total of 125. The cube is concentric with the spherical head. This configuration allows us to measure the sampling properties of the system. We also used a scaled down version ( 1 / loththe linear dimensions: 0.512 cm3). ZI. Irregular (167 cm3): Six strands on one side of the head and two on the other: maximum total of 58 electrodes. The geometric center of the volume has coordinates (-0.2, 2.7, 1.4) cm. This volume is an attempt to reproduce a possible clinical situation (although somewhat generous in terms of implanted strands). See Fig. 5 for the case of 32 electrodes. III. Irregular with Subdural Electrodes (120 cm3): Same six strands on the left side of the brain as in volume 11, plus two subdural electrodes on the same side: total of 30 electrodes. This volume is used to study the effect of adding a subdural electrode on the determination of superficial sources (see Fig. 5). IV. Three Strands (12 cm3): Subset of three strands from I1 with eight electrodes on each probe: total of 24 electrodes. This volume is used to determine the usefulness of the method in a situation where only three strands are used, in particular the ability to differentiate superficial and deep sources, even ones located outside the volume of exploration. The volume of distribution for the large cube (I) is cubic, centered at the center of the cube, and has a volume of 262 cm3. For the irregular volume (11), 167 cm3, it is also a right parallelpiped of 83 cm3. For the case with subdural electrodes (111), we used ten dipoles distributed along the outer surface (near the skull) of the volume of exploration (at an average radius of 6.4 cm, compared to

_.

812

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. 8, AUGUST 1992

TABLE I1 ERROR ON POSITION FOR DIFFERENT SAMPLING CONFIGURATIONS Volume

V (cm3)

irregular

167

(U)

N

(cm3)

(#)

DX (cm)

SD (cm)

DDV (cm)

SD (cm)

~

2.0 2.5 4.8

58 32 26 16

0.96 0.98 1.00 1.37

0.4 0.4 0.4 0.4

-0.5 -0.5 -0.5 -0.8

0.5 0.5 0.5 0.7

1.3 1.8 2.7 5.3 21

125 100 75 50 15

0.61 0.80 0.95 1.33 2.02

0.2 0.3 0.3 0.4 0.5

-0.3 -0.4 -0.7 -1.2 -1.6

0.3 0.4 0.4 0.4 0.6

1.0 (11)

cubic (1)

512

~~

No simulated noise. p s = (1, 1 , 1 ) (relative units). Average results of 100 simulations for distributed dipoles. SD: standard deviation.

Fig. 5 . Example of volume of exploration (case (11) in the simulations for 32 electrodes). The projections of the strands are represented to scale on the lateral view: the solid round markers represent electrodes implanted in the left-hand side of the brain; empty markers, the ones on the left-hand side; and triangular markers, subdural electrodes (on the left side, for volume 111). The volume of exploration (round markers, volume 11) is represented by the dotted area in the back lateral and back views. Note the large tetrahedron in the medical region.

7.5 for the head), six of which are inside the volume. For the case of three strands (IV), the volume of distributions are: first, an elongated parallelepiped similar to the volume of exploration running parallel to it at a distance of approximately 3 cm; another parallelepiped which includes the volume of exploration, i.e., sources located near or close to the volume of exploration. In the case of the simulations with spherical shell volume of distributions, the thickness of the shells is 20% of their radius, R. This radius can be larger than the equivalent radius of the volume of exploration, which means that dipoles can be located outside the volume of exploration for those simulations.

Error on Position As shown in Table 11, the error on position estimation decreases as the average tetrahedron volume, ( U ) , decreases, although the effect is less accentuated for the case of the irregular volume. By comparing the results for the

r

two volumes of integration one could deduce that the geometry of the volume of exploration is an important factor in the relationship between DX and ( U ) since the error is almost double for the irregular volume (0.98 cm versus 0.61 cm). However, it is the geometry of the individual tetrahedra which is responsible for the difference in the results: in the case of the cubic volume (I) the strands are implanted in irregular array and therefore the individual tetrahedra have more regular (compact) geometries than for the irregular volume (11). The results for measure DDV, in particular the negative values, indicate that the method has a slight tendency to find positions which are too close to the center of the volume ofexploration: for the noiseless case the average displacement is 0.5 cm. This tendency is a reflection of the nature of the center of energy measure which incorporates contributions from the energy distribution in the whole volume of exploration, and only in that volume. Therefore, for sources located near the edge of the volume the center of energy will be pulled towards the center of the volume. Supeflcial Sources: A common clinical situation is the determination of the depth of the epileptic focus, in particular for superficial sources. Using the irregular volume of exploration (11) for dipoles located at about 1 cm from the scull (near the edge of the volume of exploration) we find an average value for the difference in the distance from center of the head, DDH, to be equal to - 1.4 cm for (N =)lo dipoles. For the volume of exploration with subdural electrodes (111), the average value of DDH reduces to -0.7 cm. Case of Only Three Strands: The configuration (IV) is used to determine the resolution of the system along the axis running parallel to the strands ( y axis), that is the ability to differentiate deep and superficial sources. In the case of only three implanted strands, the resolution in the sagittal plane ( X o Z ) is null since the centers of the tetrahedra all have the same coordinates in that plane. For ten simulations with dipoles distributed outside the volume of exploration at a distance of 3 cm, but of various depths, the average error on the estimated position along y , Dy, was found to be 0.75 cm (kO.4 cm). Again for ten dis-

~

813

LEMIEUX AND LEDUC: SOURCE ESTIMATION BASED ON ELECTRIC FIELD CALCULATION

tributed dipoles, this time located inside or close to the volume of exploration, this error was found to be 0.68 cm ( f 0 . 6 cm). We therefore conclude that the method is able to reliably determine the depth of a source which is relatively close to the implanted strands.

Error on Dipole Orientation The major factor in the estimation of the dipole orientation is the geometry of the volume of exploration, as demonstrated by the results contained in Tables I and 111, although the effect is not dramatic. As noted before, for the error on source position, for volume I1 the tetrahedra have more irregular geometries, which affects the estimate of the tetrahedron field. This is reflected in the results of Table 111. Moreover, the calculation of the dipole orientation, in the noiseless case, depends more on the total number of electrodes than on the density of electrodes, as demonstrated by the results in Table 111: as the density increases by a factor of 1000 (when we use the scaled down version of the same volume of exploration), the error less than doubles. This can be understood by considering the fact that our calculation of the dipole orientation is really a measure of the field configuration. In the case of a small (mathematical) dipole, this configuration is the same at all scales; therefore whether one uses a large or small array of electrodes does not significantly affect the results. Of course, this is only absolutely true for mathematical dipoles in an infinite conducting medium. The increase in the error is due to the finite dimensions of the conducting medium. This was verified by performing the same comparison for an infinite medium: we have found that the effect of the head boundary is to reduce the error on the dipole orientation compared to the infinite medium case, at least for a cubic volume of exploration. This is the reason why the error is larger for the reduced volume of exploration, for which the field effectively resembles more the infinite medium case. For the case of only three strands (IV) the average error on dipole orientation for the sources located close to the strands was 39" (+19") (averaged over three orientations, ten dipoles each: (1, l , l), (1, 0, 0) and (0, l , 0)). This shows that it is possible to obtain significant estimates of dipole orientation even for this minimal set of measurement points and very nonspherical volume of explogration. Error on Dipole Magnitude The results of Table I11 demonstrate that the method gives a fairly accurate estimate of the dipole magnitude for the irregular volume: EDM = 10% ( 30 %). The major factor in the estimation of the dipole magnitude is the geometry of the volume of exploration. As the results in Tables I and I11 demonstrate, the total number of electrodes, N, and the average tetrahedron volume, (U), for a given volume of exploration have relatively little influence on the estimated dipole magnitude (expect for the

TABLE I11 ON DIPOLE ORIENTATION AND MAGNITUDE FOR THREE VOLUMES OF ERROR INTEGRATION AND DIFFERENT SAMPLINGS Volume

V (cm')

irr. (11)

167

cubic

512

(1)

cubic (1' )

0.5

ED0 (deg)

SD (deg)

58 32 16

16 16 17

12 12 12

125 15 15

7 9 10

7

100 75 15

12 13 20

(U)

N

(cm')

(#)

1.o 2.0 4.8

1.3 2.1 21 0.002 0.003 0.021

5

6

8 10 12

EDM (rel.)

0.09 0.1 0.1

SD (rel.) 0.3 0.3 0.3

-0.6 -0.6 -0.2

0.2 0.2 0.3

-0.03 -0.01 0.3

0.2 0.3 0.3

No simulated noise. The third set of results (0.5 cm') is for the scaleddown version of the cube. Average results of 100 distributed dipole simulations. Averages for 4 dipole orientations: (1, 1, I), (1, 0, 0), (0,1 , 0), (0, 0, l), except for V = 512 cm', 0.5 cm': 2 orientations: (1, 1, l), (1, 0, 0). SD: standard deviation; rel.: relative units.

reduced cube where EDM becomes positive for small sampling). Another important factor is the proximity of the medium boundary as indicated by the results for the reduced cube as compared with the cube: the estimation of the dipole magnitude via (6) and (7) is degraded by the proximity of the head boundary, as opposed to the effect on the estimated dipole orientation (mentioned in the previous section). Note that the relatively large standard deviations reflect the wide range of values of EDM; in fact for some of the distributed dipoles the error is zero, while for others it is 100% (EDM = 1). For the case of three strands (IV) the error on dipole magnitude is found to be large: 0.66 (f0.3). Therefore the dipole magnitude is under-estimated.

C. Influence of Dipole Position on Errors To assess the effect of the position of the dipoles relative to the volume of exploration, we performed distributed dipoles simulations in which the volume of distributions where spherical shells of various radii, concentric with the center of the volume of exploration. As we mentioned in the section on the error on position, the estimated dipole positions tend to be too centric relative to the volume of exploration. Therefore, we expect the error on dipole position to be larger for dipoles located far from the center of the volume. As Fig. 6 demonstrates, both the error on position and on dipole orientation increase as the distance relative to the volume center increases. We note that even for dipoles outside the volume of exploration, R > 3.5 cm, it is still possible to get an estimation of dipole orientation in some cases: this reflects the fact that the measure of dipole orientation is not meaningless for dipoles located outside the volume, as long as the orientation of the dipole relative to the volume of exploration is such that the field's geometry inside the volume of exploration is relatively uniform, and parallel to the dipole vector.

814

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. 8, AUGUST 1992

I

I

+

I

I

I

l+p‘ ’

0 1 ‘ 0

’

’

0.25 0.5 0.75 1.0

N

”

I

lo

10

Ar

(a)

DISTANCE FROM CENTER (cm) 1

(a)

0

-5

I

n W

-10

I 0

,

0.25 0.5 0.75 1.0

N ’I

I

10

Ar

0 0

2

I

I

I

4

6

8

I

I

1

0

1

2

DISTANCE FROM CENTER (cm) (b) Fig. 6. Errors on dipole position and orientation as a function of dipole position relative to volume center. No simulated noise. N = 32, V = 167 cm3. Average results for 50 distributed dipoles in spherical shells of radius R . p = (1, 1, 1). The boundary of the integration volumes is at about R = 3.5 cm.

D. Influence of Noise Instrumental Noise: Fig. 7 shows the results of distributed dipoles simulations for a range of values of the noise level (1.25 ISIN,,,,, = (A,)-’ I00). For dipole position and dipole orientation the results for S/N,,,,, = 0.1 are shown to illustrate the effect of dominating noise. As expected all errors tend to increase with noise level. The average error on dipole orientation increases from 18 O to 58” (for A, = l), which is still a relatively good result, since it allows in the worst case to determine the dipole direction. The dipole amplitude is generally underestimated for low noise levels and becomes overestimated as the level of noise increases. Estimated dipole positions tend to be more central in a noisy environment; this corresponds to the fact that the noise field must be more or less uniform in amplitude and therefore has an effective center of energy near the center of the volume. DX and E D 0 tend to saturate, since their possible range of values is limited. Physiological Noise: Simulation results for distributed dipoles simulations with varying numbers and magnitudes

(b) Fig. 7. Average errors on dipole parameters as a function of the relative noise level. (a) Error on position, DX, and error on orientation, EDO; (b) difference in distance from center, DDV, and error on dipole magnitude, EDM . N = 32, V = 167 cm3, p‘ = (1, 1 , 1). The vertical bars represent the standard deviations of the errors for 1 0 0 distributed dipoles.

of noisy dipoles are shown in Table IV. The noisy dipoles are identical for all simulations with the same N,, (therefore the “noisy field” is identi~al)~. The results demonstrated that the proposed method gives valuable results for a [S/NIphys ratio of 0.5 (in the case of 20 noisy dipoles). This helps clarify the notion of a “dominant” source; more precisely, to which degree the field in the volume of exploration is dominated by the activity of the source of interest in a noisy environment. It is also of interest to note that the dependence of the errors on the total noisy dipole amplitude, N,,p,, is identical to the dependence on the relative noise level, A,, for the case of “instrumental” noise. We therefore conclude that the two types of simulated noises are equivalent with regard to the proposed method.

E. Simulation of Increase in Dipole Magnitude The simulation results illustrated in Fig. 8 are for series simulations: the calculation of the source parameters is performed for the same source but for different randomly generated noise values of the same amplitude. Volume I1 ’By taking the average error for two opposite source orientations, we ensure that the error on source orientation, DA, is 90” for a completely noisy field.

LEMIEUX AND LEDUC: SOURCE ESTIMATION BASED ON ELECTRIC FIELD CALCULATION

TABLE IV

EFFECTOF PHYSIOLOGICAL NOISE N”

DX

NnP. (A cm)

(#)

[SIN],

(cm)

ED0

DDV (cm)

(deg.)

EDM

0

-

OD

1.2

-0.46

14

0.32

5

0.1 1 .o 5.0

10 1 0.2

1.2 1.3 1.9

-0.46 -0.73 -1.4

14 22 70

0.31 0.28 -0.54

0.1

10 1 0.5 0.5

1.2 1.3 1.6 2.4

-0.46 -0.83 -1.2 -1.5

14 18 29 82

0.31 0.30 0.23 -2.8

10 1 0.25 0.025

1.2 1.4 2.2 2.9

-0.46 -0.65 -0.37 0.0

14 19 70 88

0.32 0.30 -0.23 -11.0

;:;

20

20.0 0.1

i:;

40

40.0

Irregular volume of exploration (11). N = 32. Average for ten distributed dipoles, two source orientations: p’ = (1, 1 , l ) , ( - 1 , - 1 , - 1). First row: noise.

2.5 C

z t--

r

J

I

2

1.5

0z

1

a

5

0.5

0

W

O

__I

9

LL

n

-0.5

815

tuations of the parameters when the magnitude of the dipole source increases above the noise level by a factor 10 (at “instant” #51). As seen in Fig. 8(b), (c) the relative amplitude of the “noise” fluctuations decreases in the case of dipole orientation and source position coordinates when the magnitude of the source increases. For dipole moment components [Fig. 8(a)], and therefore total dipole moment, the increase in the magnitude of the source is reflected by the sharp increase of these parameters, while the absolute amplitude of the fluctuations remains constant (as does the amplitude of the noise voltages). For the first half of the series the average E D 0 is 33” (SD = 22”), EDM = -0.86, DX = 1.5 cm; for the second half, ED0 = 12” (SD = 2”), EDM = 0.1, DX = 0.9 cm. Therefore, in this case even for a signal to noise ratio equal to 1, the method can give meaningful results. Here the main objective is not to attempt to reproduce actual behavior of the source but to show the sensitivity of the method for the various parameters. For example, we know that the frequency band width of the depth EEG is mainly in the 5-50 Hz region. We can therefore assume that the actual time behavior of the estimated parameters from real EEG should also be in that range and that the data reproduced in Fig. 7 varies much too fast considering a 200 s-’ sampling rate.

-1 -1.5

0

1 0

20

30

40

50

60

70

80

90 1 0 0

V . DISCUSSION AND CONCLUSION We have presented the theoretical basis as well as a 280 simulation study of a method for the calculation of source 240 parameters based on the calculation of the electric field U 200 from depth EEG data. We have demonstrated that valu1 6 0 able estimates can be obtained through the proposed U + z 1 2 0 method. Furthermore, the tetrahedral method for the calw 8 EO culation of the electric field is a simple device which is 40 well adapted for the present purpose. 0 0 1 0 2 0 30 40 50 60 70 80 90 1 0 0 The method we propose for the estimation of source parameters is independent of any assumption on the phys4 E I ical nature of the source beside it being of limited extent. v - 3 This is in contrast with other currently proposed methods CrJ s 2 which rely on the assumption of a point dipole. In effect, n z we propose to estimate directly the source parameters 8 ’ based on the assumption that a unique source is located S z O inside the volume generated by the implanted electrodes. e--1 0 m 0 A potential advantage of the proposed approach is the fact n -21 ’ ’ ’ ‘ ’ ’ ’ ’ ’ I that the calculation of the source parameters is nonitera0 1 0 20 J?,N450TA5NoT.60 loo tive and avoids problems of unicity of the solution. Al(C) though the implementation presented here is based on a Fig. 8. Sudden increase in dipole magnitude. (a) Dipole moment compoparticular geometry of stereotaxically implanted multiplenents. (b) Dipole orientation. (c) Dipole position coordinates. Volume 11: N = 32, V = 167 cm3. Simulated instrumental noise: SIN,,,,, is 1 for the electrode strands, that is, parallel strands on each side of 50 first calculations (“instants”), and 10 for the last 50. The true values the head, the principle applies to any configuration of of the dipole parameters are indi2ated by horizontal lines (suffix r): cr = electrodes within certain limits. ( - L O , I), e, = 4 5 , 4 , = MO, x, = ( - 1 . 0 , 2 . 5 , 2 . 5 ) . In summary, the simulation results have demonstrated that for sources located inside the volume of exploration was used with 32 electrodes. These serve to illustrate the the proposed approach can give valuable estimates of the influence of the noise on the estimated dipole parameters source parameters, even for a relatively small number of and more specifically the relative amplitude of the fluc- electrodes, electrode locations far from ideal, and rela(a)

320

h

v

e

h

U

,