|Figure V.1: Position of the markers indicated in the MRI scans
In this chapter, a method will be described, providing a tool for combining structural information obtained from MRI-images with functional data from MEG and EEG. To be able to do so, one has to convert the coordinates of the location of the source area represented by a current dipole found by means of an inverse procedure into the coordinate system of the MRI-scan. Different models for the head that can be used in the source localisation are extracted from the MRI. It will be shown that the source location found, does depend on the model used. Moreover, integration of MRI with MEG/EEG provides us with a better insight into the location of an active region, but also allows the evaluation of the results. It will be shown that by combining MEG and EEG, one may attempt to estimate the size of the active area in the cortex. For further verification of the locations found, PET scans will also be matched with the MRI as well. The software to perform these tasks has been developed at the Biomagnetic Centre Twente, and comprises a single program which is called the Biomagnetic Imaging System (BIS).
In matching MEG and EEG with the MRI, differing coordinate systems have to be taken into account. The most straightforward of which is the MRI coordinate System. In an MRI-dataset each voxel has its own Cartesian coordinates, of which the Z-coordinate reflects the number of the slice. In the scans used, all coordinates can run from 0 to 255, in units of one mm. The origin is located at the lower left side of the head, at the back. The scanned area can be viewed as a cube which contains the head.
A second coordinate system is related to the head itself, the head coordinate system. It is defined by X- and Y-axis running through anatomical landmarks. The Z-axis is taken perpendicular to the X- and Y-axis, and always points to the top of the head. The origin is always located in the head, somewhere between the ears. Different landmarks are in use. De Munck (1989) used the inion and the external auditory meatus, with the origin halfway between the line connecting the left and right auditory meatus, with the Y-axis running towards the right. The X-axis runs from the origin to the inion. For convenience often the preauricular points are used instead of the meatus. We will call this the modified head system. Also, it is possible to use the nasion instead of the inion. The X-axis runs in this case towards the nasion, and the Y-axis points toward the left-ear. This will be called the reversed head-system. In the reversed head system it is assumed that also the preauricular points are used. The program called BIS can deal with these three variants of the head coordinate system.
To match the head coordinate system with the MRI coordinates, the MRI scan is performed, with markers attached at the nasion and preauricular points of the subject. These markers consist of small plexiglass holders, filled with the fluid Gadolinium, which shows up very well in the MRI-scans. If they have been found in the scans, they have to be identified for the computer, so it can set up the coordinate system (see figure V.1).
|
Figure V.1: Position of the markers indicated in the MRI scans
The origin, and the axis of the head coordinate system are subsequently expressed as vectors in MRI-coordinates. Any coordinate in the head coordinate system is given in MRI-coordinates by:
where Hi, i=x,y,z, is the location vector in Cartesian head coordinates, the vectors Xi, Yi and Zi represent the three axis of the head coordinate system, expressed in MRI-coordinates, and Mi is the resulting location vector in MRI-coordinates.
Since it is very uncomfortable for the subject to be scanned with a marker on the inion, we decided not to use such a marker. Therefore, if the non-reversed head coordinate system needs to be used, the inion has to be located afterwards in the scan. Although a rough estimate can be made of the inion on visual inspection, it is preferable to measure the distances of the inion to the preauricular points and the nasion. With these distances known, it is possible to calculate the location of the inion in the MRI-scans. This is the same method De Munck used to calculate the position of electrodes with respect to the head-coordinate system (De Munck et al., 1991). He derived the following formula:
where a is the distance from the origin of the head coordinate system, to the inion or nasion, whichever is used, b is the distance from the origin to one of the earpoints used, where it is assumed that the distance between the earpoints is 2b. dl, d2 and d3 are the measured distances from the point of interest, to respectively the inion, the right earpoint and the left earpoint in the case of the non-reversed coordinate system (see figure V.2).
For the reversed coordinate system the distances are measured to the nasion, left earpoint and the right earpoint respectively. One can use this method to find the inion. The distances of the inion to the landmarks of the reversed-coordinate system are used for this purpose. Once the position of the inion has been established, the non-reversed coordinate system can be used. Formula V.2 has two solutions. The user will have to decide which of both solutions is the proper one.
V.2.1 Sphere fitting
The inverse procedure based on EEG, described by De Munck(1989) uses concentric spheres as a model of the head. The procedure has been adapted to handle MEG as well. In this case one sphere suffices. Usually, De Munck's procedure uses a sphere fitted to the electrode positions as the outer one. The electrode positions are calculated from the distances measured from each electrode to the anatomical landmarks used. The resulting dipole(s) are expressed with respect to the centre of the spheres. To reproduce the locations of these dipoles in the MRI scans, this centre has to be expressed in MRI-coordinates. The position of the centre of the sphere is given in the head coordinate system, and can therefore be translated to the MRI-coordinates as described above. For dipole positions to be translated to MRI-coordinates, the centre of the sphere has to be taken as the origin of the coordinate system. We will call this the sphere coordinate system. Obviously, this coordinate system is only displaced with respect to the head coordinate system. The directions of the axis remain the same.
Instead of fitting the sphere to the electrodes, a sphere can be fitted locally in the MRI-scan. Various regions of interest can be selected, like the back of the head, the side of the head or the top of the head. Points are taken from this region using the same method as described for generating a realistically shaped model of the head in chapter IV, except that points are only taken from a limited region of the head. Once points are selected, a sphere can be fitted to these points, using an iterative procedure which finds the best fit of the origin and the radius of the sphere using the method of steepest descend. Once the sphere has been found, an inverse coordinate transformation can be used, to express the location of its centre in the head coordinate system. These coordinates can be used in the inverse procedure. An example of a sphere fitted in the MRI scan can be seen in figure V.3.
Figure V.3: Example of a sphere fitted in the MRI scan at the back of
the head
V.2.2 The measurement positions
For any good source localisations, it is important to know the measurement
positions with respect to the head. For EEG-measurements, the position of
the electrodes of the head can be determined by measuring the distance of an
electrode to each of the three points which are used to define the head
coordinate system. Equation V.2 then gives the coordinates of the
electrode for the head coordinate system.
For an array of magnetometers the situation is quite different. Although a
mechanical solution is easy to implement, its accuracy leaves much to be
desired. Therefore many groups are using coil sets which are attached to
the head (Ernč et al., 1987; Ahlfors and Ilmoniemi, 1989; Hulsteyn et al, 1989,
Higuchi et al., 1989; Incardona et al., 1992) . A current is sent through each
coil in turn, which generates a magnetic field which is measured by the
magnetometer. From these measurements, the position of each of the
coilsets can be determined with respect to the pick-up coils of the
magnetometer. Since the position of the coils with respect to the
head-coordinate system ran be determined in the same way as the electrodes, the
magnetic measurement position is known with respect to the head.
Such a system is in the final stages of development at the Biomagnetic Centre Twente.
V.2.3 Displaying the dipole in the MRI
Once an inverse solution has been calculated, usually with respect to a sphere, the coordinates of the equivalent current dipole have to be entered into the program. The BIS program can accept dipole locations in spherical or Cartesian coordinates, with the origin being either the centre of the sphere, or the origin of the head coordinate system. It is also possible to enter the dipole directly in MRI-coordinates. This is useful in case the inverse procedure is carried out using the boundary element method, because the elements are derived from the MRI, and their coordinates correspond to the MRI-coordinate system. The location of the resulting dipole would therefore also be expressed in MRI-coordinates. The program transforms the dipole position to MRI-coordinates and retrieves the proper MRI-slice, indicated by the Z-coordinate of the dipole. The dipole is indicated with a dot in this slice. A sagital slice of this MRI-scan is also retrieved to indicate the vertical position of the slice, and the dipole location in this slice. Figure V.4 shows the position of an equivalent dipole indicated in the MRI.
Figure V.4: Example of an equivalent dipole shown in the MRI
As indicated earlier, one can use various models for the description of the volume conductor, i.e. the head. The most often used model, the sphere, has the advantage that there is an analytical solution of the forward problem. Something which significantly reduces the computational effort in solving the inverse problem. Furthermore, the magnetic field is not influenced by volume currents. For the magnetic case, this model is equivalent to the concentric spheres model. The more realistically shaped descriptions of the head need numerical calculations to solve the forward problem. To illustrate that the volume conductor influences the location of an estimated equivalent dipole, the estimation based on a sphere model will be compared to one based on a more realistically shaped model of the head. The equivalent current dipole location will be calculated from an actual magnetic measurement.
The realistically shaped volume conductor model is a bounded homogeneous conductor with the shape of the inner surface of the skull. This should be equivalent to a compartment with the shape of the brain, since the layer of fluid between the brain and the skull is very thin. This model is within an error of 5% equivalent to a model consisting of three realistically shaped homogeneous, isotropic compartments, describing the brain, the skull and the scalp, as argued by Hamalainen and Sarvas (1987). The inverse problem is solved by means of the boundary element method (BEM). The BEM is based on the following integral equations derived by Barnard et al. (1967) and Geselowitz (1970) as discussed in Chapter III.
The first two equations describe the electric potential 
,
while the second two are valid for the magnetic field B. Ns denotes the number
of surfaces. 
is the electric
conductivity outside the j th surface, while 
is the electric conductivity on the inside of that surface. a. is the
conductivity in the source region. The source is located at position r0
.
These equations show that the contribution of the volume currents in a homogeneous, isotropic volume conductor can be considered equivalent to the influence of secondary sources which lie at the boundary. The orientation of these secondary sources is normal to the boundary and their value is linearly proportional to the electric potential at that location. The magnetic field due to these secondary sources has to be added to the magnetic field due to the current dipole. To compute magnetic fields and electric potentials numerically, the integral equations mentioned above must be approximated by summations. The surface involved is presented by many automatically generated triangles and the potential is assumed to be constant over each triangle.
In our algorithm the fields are computed with the so-called centre of mass approach in which the potential of the centre of mass of a triangle is ascribed to the potential of that triangle. The set of linear equations resulting from the numerical expression is singular. This singularity is first removed by deflation of the associated matrix (Lynn and Timlake, 1968). Usually, the matrix equation obtained is solved by iterative procedures. However, especially in the case that the volume conductor consists of a single compartment, it is more convenient to solve the system of linear equations by computing the inverse of the matrix (Oostendorp and Van Oosterom, 1991). For the magnetic case the matrix for the discretised integral equation has to be created and multiplied by the inverted potential matrix. This System matrix remains the same for any source within the volume conductor chosen. Consequently, once this inverse is obtained, the forward problem can be solved very fast for any choice of the location of the dipole within that volume conductor. The inverse solution is based on the algorithm of Marquardt (1963). That dipole is estimated for which the model predicted distribution of the field, Bcomp is closest to the measured distribution of the field, Bmeas The function which is minimized reads:
The
forward solution by means of the BEM is tested by computing the magnetic field
distribution due to a current dipole p within a homogeneous sphere
and by comparing the results obtained with those obtained by the analytical
expression. The surface of the sphere is discretized with 320
triangles. Magnetic fields components are computed in 19 grid points,
which are chosen such that they coincide with those measured with our
magnetometer system (see figure V.5). The distances of the observation points
are 2.5 cm from the surface of the sphere with radius r = 10 cm .
The numerical error in the computed magnetic field distribution (in least-squared error sense) turned out to be smaller than 5% for tangential dipoles located in a range of 10 to 90 mm from the Centre of the sphere.
The same sphere is used to test the inverse solution. Because radial dipoles do not contribute to the outer magnetic field, a penalty function is added to the sum of quadratic differences function F. This penalty function P reads
where Pr is the instantaneous radial component of the current dipole.
The function to be minimized is F', where
F'=F+P
The fitted dipoles are found within 0.5 mm from the known dipole locations and the strength of fitted dipoles differ about 1% from the original ones.
As an example the inverse procedure is used to localise a source from visually evoked magnetic fields. These measurements are part of a study on the effects of visual spatial attention on early components in Visual Evoked Potentials and Fields.
Eight subjects with a mean age of 23 years, with normal or corrected to normal vision, were investigated. The stimuli consisted of sets of two letters (consonants), presented to the left or the right of a central fixation dot that was continuously present. The distance between the subject and the monitor was kept constant at 1.5 meters. The visual angle between the dot and the letters was 2.5 degrees, the letters were 0.9 degrees in height. The subject was instructed to pay attention to the left or the right visual field while fixating on the central dot. They had to react with a button press, as fast and accurately as possible when a letter from a memory set was presented in the attended field. Magnetic measurements were performed in 36 grid points. Memory load did not affect the early components like P1 and N1.
The
source analysis described in this paper is carried out on the first large
component in the VEF the P1m, the magnetic counterpart of the electrical
P1-component, which shows up for all stimulus categories, either attended or
unattended. The comparison between head-models as described here was done
for one of the subjects of this experiment.
The BEM method described above is used to estimate the source from the P1m component of the VEF. The model of the head is shown in figure V.6. The head is also described by a sphere which is locally fitted to the scalp. The results obtained are displayed in figure V.7. Although it is the common opinion that both models (i.e., the sphere locally fitted to the scalp and the inner surface of the skull) are adequate as models for the head to localize activity from MEG measurements, we observe a difference in dipole location of about 1 cm.
Figure V.7: Equivalent dipole locations calculated using a
realistically shaped and a spherical model
In order to see if this distance can be explained by noise in the measurements or by the assumption that there is only one dipole present, we did the following simulation: The dipole found was used in a forward solution using the realistically shaped model. The resulting field distribution was used in an inverse solution with the sphere model. The dipole obtained in this way was located at the same position as the one which was obtained with the sphere model from the measured fields. This location is at a distance of 1 cm from the original dipole which was used in the forward solution.
The different dipole positions found by using either the sphere or the realistically shaped model of the head illustrate the strong influence of the volume conductor on the measured magnetic fields. Furthermore, it is impossible to determine from the MRI which solution is better. Both dipoles are found in reasonable locations. The application of the BEM, combined with MRI based models, enables us to improve the accuracy of dipole estimations with acceptable computational effort by enabling us to use more detailed descriptions of the head.
Although the generation of the realistically shaped homogeneous model is fully automated, the localization procedure based on such a model is still time consuming and costly. Consequently, it will be advantageous if not the actual geometry of the inner surface of the skull of the individual subject is needed to solve the inverse problem for MEG. For this reason the adequacy of a standard realistically shaped model has to be tested. Such a standard model might be obtained by scaling the inner surface of the skull of an arbitrarily chosen subject. From a computational point of view it is to be preferred that the same scaling factor is used for three orthogonal directions, because in such an operation the Maxwell equations are conserved and so is the matrix describing the linear system of equations given by equation III.7. As a consequence the inversion of the (deflated) matrix used in the inverse solution based on the boundary element method has to be carried out only once. The variation in the lengths and widths of skulls is rather limited. Van Veenendaal (1982), a coworker at our university, measured the skulls of 29 women and 64 men who died in the age between 20 and 60 years. The average distance between inion and nasion was 182 ± 8 mm, the average distance between the ear canals was 144 ± 7 mm
In order to test the adequacy of a transformation of a standard model by
scaling, simulations were carried out. Two models were derived from the
MRI of two different subjects. A dipole was placed within one of the
models. The magnetic field was calculated in 57 points, based on one model
and the inverse solution was carried out using the scaled model. The
distance in localization is called the error. The forward and inverse
procedures were based on the boundary element method. The number of
triangles used was about 650. The errors due to the discretization are of
the order of 1 mm. From the calculated magnetic field distribution the
dipole was localized using the scaled model. The relative position of the
pick-up coils with respect to the head was the same for the two models
used. When the scaling factor was taken to be the same for all three
orthogonal directions, the errors (dependent on the position and orientation of
the dipole) could be quite large, i.e. up to 2 cm. Taking a different scaling
factor for each of the three orthogonal directions led to better results.
The three scaling factors used in our example were less than 10%. The two
triangulated models are shown in figure V.8. 
The
results for a tangentially oriented dipole at different depths in the occipital
region near the midsagital plans are shown in figure V.9.
Also the results are shown for the case that instead of the scaled model a
sphere was used. For deeper sources the scaled model gives better results
than the spherical model. It turns out that for superficial sources the sphere
gives rise to errors which are of the same order or even smaller than the
scaled model dependent on the orientation of the dipole. This is only true for
the occipital region. From the results mentioned above we
have to conclude that an individual realistically shaped model is superior to a
standard realistically shaped model or a spherical model.
V.4.1 Determination of the strength of a dipole
To be able to estimate the area of cortex in the brain, activated at a
certain instance, the strength of the equivalent dipole has to be known.
In a spherical model, the inverse solution of an MEG measurement only yields the
tangential component of the equivalent dipole. However, using only EEG to
estimate the active patch of cortex, does not give reliable results, as ran been
seen from an investigation by Stok (1986). In order to study the influence
of the head model parameters, he calculated the potential distribution due to a
current dipole, and used the data obtained to solve the inverse problem.
The inverse solution was based on a reference model. The conductivities
used for the scalp, skull, fluid and brain compartments were respectively 0.33,
0.0042, 1.0 and 0.33 S/m as given by Cohen et al. (1990). Values of the
radii of the four concentric shells were 75, 69, 66 and 64 mm. Stok used
current dipoles which were either radially or tangentially oriented. He
found that both the sphere radii and the conductivities especially influence the
strength of the EEG dipole, and much less its location. To study the
influence of the changes in conductivities on dipoles with an arbitrary
orientation, we performed similar simulations. In these simulations the
conductivities were allowed to change by a total of 30% in steps of 5% and the
radii by 2.5% in steps of 0.5%. From all possible combinations the worst-case
maximum and worst-case minimum potential were selected. The influence of
the changes in the conductivity and radii upon the ratio between the tangential
and radial dipole were computed for the worst-case combinations of the
conductivities, for sources at various depths. The results are given in
table V.I. R is the distance from the centre of the sphere to the dipole.
As shown in this table the influence is in all cases smaller than 5%. In
other words, the ratio between the radial and tangential components of the
EEG-based equivalent dipoles is hardly influenced by the model parameters,
although both the strength and the position of the equivalent dipoles are
influenced. The influence of uncertainties in the radii turned out to be
such that for both the worst-case maximum and minimum the orientation of the
equivalent dipole did not change more than 0.5% for all depths and all angles
between the radial and tangential components. The tangential component of
the equivalent dipole can be derived from MEG data. As shown by the
simulations, the angle 
can be reliably derived from EEG data. The fact
that the orientation of the dipole is preserved was also found by Cuffin et al.
(1991), who studied the influence of the real head by means of activated
implanted electrodes. As a consequence MEG measurements can be used to
localize the current dipoles and to estimate the strength of the tangential
component. EEG measurements can be used to compute the orientation of the
current dipole. The strength is equal to the tangential component times 1/cos
). In other words, using both EEG and MEG, all six parameters describing
the dipole can be found more reliably than with EEG or MEG alone. Once the
strength and the orientation of the dipole are known, the patch of cortex which
is represented by this equivalent dipole can be depicted in the MRI, using data
from electrophysiology.
| relative errors in
orientation (1/cos ) in
% |
||||||
| R = 30 mm | R = 50 mm | R = 62 mm | ||||
| Worst case | ||||||
| |
Min . | Max. | Min. | Max. | Min. | Max. |
| 0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
| 15 | 0.1 | 0.2 | -0.5 | 0.1 | 2.8 | 0.1 |
| 30 | -0.2 | 0.1 | -0.8 | 0.3 | -4.6 | -0.1 |
| 45 | -0.3 | 0.3 | -0.9 | 0.3 | -4.1 | -1.1 |
| 60 | -0.6 | 0.4 | -1.2 | 0.4 | -3.6 | -2.0 |
| 75 | -1.5 | 0.5 | -2.6 | 0.8 | -4.2 | -2-2 |
| 90 | - | - | - | - | - | - |
Table V.1:Influence conductivities on EEG-equivalent
dipoles. R-distance centre sphere to dipole. Errors in 1/cos
in %
Williamson et al. (1991) determined the strength of the dipole from MEG en
MRI, by determining the position of the dipole from the MEG, and determining the
orientation and thus the strength of the dipole from the fact that the direction
of the dipole should be perpendicular to the cortical surface. The
correctness of the strength depends in this case entirely on the accuracy of the
dipole localisation. As shown by the comparison of the head models above,
different localisations may yield different but both quite possible
solutions. Therefore this method does probably not yield reliable results.
V.4.2 Determination of the active area of cortex
In the cortex, the pyramidal neurons form the main source of electric activity. Due to the axial alignment of their apical dendrites and to the lamellar organization of the synaptic inputs, the population of pyramidal neurons generates dipole fields (Lopes de Silva and Van Rotterdam, 1987). It is interesting to establish a quantitative relationship between the strength of the estimated dipole and the dimensions of the cortical dipole layer generated by the local pyramidal cells. In order to realize this objective, some assumptions have to be made about the quantitative characteristics of the cortical dipole layer. While recording electrocorticograms, Freeman found the largest transcortical potential differences to be I mV over a distance of 1.5 mm. Assuming the specific resistance of the cortical tissue to be, on average, 25 ohm-mm, one can estimate the transcortical dipole current density to be about 270 nAmm-2 (Freeman, 1975, p.454). In the case of visually evoked fields, Stok (1986) found for the largest equivalent dipoles measured, having a strength of 100 nAm, that the area of the cortical layer activated should be about 400 MM2, assuming that the poles were separated by a distance of 1 mm, and that the dipole current density was 250 nAm M-2 . Thus taking into account these basic quantities, we can estimate the active cortical area corresponding to a given equivalent dipole. As mentioned before, the location of the estimated dipole is derived from MEG data, as well as the strength of its tangential component. The orientation of this dipole is derived from EEG data. Combining both yields location, strength and orientation of the equivalent dipole. The direction of the dipole corresponds to the direction of the intracellular current. However, the direction of this current can only be known if the pattern of synaptic activation underlying a given field is also known from basic electrophysiological studies. For example, if we assume that an evoked field component is primarily caused by excitatory synapses in the middle or deep layers of the cortex, the main intracellular current is oriented towards the cortical surface. Knowledge of the direction of the dipole can be used to reposition it within the error estimate in order to obtain the correct solution.
It follows from simulation studies that the electric potential and the magnetic induction, for any homogeneous dipole layer with a circular rim (e.g. a disc, a hemisphere or a homogeneous dipole layer having the shape of a half prolate spheroid) is equivalent to a single current dipole at the centre of the circle of the ground plans, at distances comparable to the dimensions of the layer. Only a dipole layer, where all the conditions mentioned above are met within the error estimate, can be a solution.
The segmentation of the MRI-images ran often discriminate between white matter (cortex) and grey matter within the brain. This allows the computer to indicate a patch of cortex corresponding to the dipole strength, under the assumption that that patch of cortex is rotational symmetric around an axis formed by the direction of the dipole. The computer generates an enlarged view of the brain area around the location of the equivalent dipole and tests to see if the dipole location is within the grey matter (cortex) or fluid. In the latter case the dipole could be located in a gyrus, which can be a valid solution. The computer algorithm searches in the direction of the dipole orientation for grey matter first. When found, a region growing operation is started. The size of the region is dependent on the strength of the dipole. The choice to search in the direction of the dipole orientation is arbitrary, since it is not known what direction the current inside the neurons have, towards or away from the cortical surface. An example is shown in figure V.10.
Figure V.10: Example of an estimated area of cortex from the strength
of the equivalent dipole
It has been shown that MRI and MEG/EEG can be of mutual benefit to each other. The structural information contained in the MRI has been used to generate models, to display and evaluate the location of equivalent current dipoles and to estimate the active region of the cortex. The strength of the dipole has to be determined from both MEG and EEG. Determining the strength from MEG and MRI, as proposed by Williamson et al. (1991) does probably not lead to reliable results. The MRI can not be used to discriminate between two possible solutions, but it can be used to modify a solution to satisfy electrophysiological restraints. Although one has to take care not to suggest more knowledge or precision than available, the combination of these three techniques may prove very fruitful in future. Another technique may be added to these to add further information, the PET scan. This will be discussed in the next chapter.
(c) MEG, EEG and the integration with Magnetic Resonance Images, H.J. Wieringa, 1993
[<<Chapter IV][Contents][Home][Chapter VI: Matching MRI with PET data>>]