A stereotactic method for image-guided transcranial magnetic stimulation validated with fMRI and motor-evoked potentials

www.elsevier.com/locate/ynimg NeuroImage 21 (2004) 1805 – 1817

A stereotactic method for image-guided transcranial magnetic stimulation validated with fMRI and motor-evoked potentials S.F.W. Neggers, a,* T.R. Langerak, b D.J.L.G. Schutter, a R.C.W. Mandl, c N.F. Ramsey, c P.J.J. Lemmens, a and A. Postma a a

Department of Psychonomics, Helmholtz Institute, Utrecht University, Utrecht, The Netherlands Department of Computer Science, Utrecht University, The Netherlands c Department of Psychiatry, Rudolf Magnus Institute, University Medical Center Utrecht, The Netherlands b

Received 8 October 2003; revised 3 December 2003; accepted 3 December 2003

Transcranial Magnetic Stimulation (TMS) delivers short magnetic pulses that penetrate the skull unattenuated, disrupting neural processing in a noninvasive, reversible way. To disrupt specific neural processes, coil placement over the proper site is critical. Therefore, a neural navigator (NeNa) was developed. NeNa is a frameless stereotactic device using structural and functional magnetic resonance imaging (fMRI) data to guide TMS coil placement. To coregister the participant’s head to his MRI, 3D cursors are moved to anatomical landmarks on a skin rendering of the participants MRI on a screen, and measured at the head with a position measurement device. A method is proposed to calculate a rigid body transformation that can coregister both sets of coordinates under realistic noise conditions. After coregistration, NeNa visualizes in real time where the device is located with respect to the head, brain structures, and activated areas, enabling precise placement of the TMS coil over a predefined target region. NeNa was validated by stimulating 5
5 positions around the ‘motor hotspot’ (thumb movement area), which was marked on the scalp guided by individual fMRI data, while recording motor-evoked potentials (MEPs) from the abductor pollicis brevis (APB). The distance between the center of gravity (CoG) of MEP responses and the location marked on the scalp overlying maximum fMRI activation was on average less then 5 mm. The present results demonstrate that NeNa is a reliable method for image-guided TMS coil placement. D 2004 Elsevier Inc. All rights reserved. Keywords: fMRI; Transcranial magnetic stimulation; Neural navigator

Introduction Transcranial magnetic stimulation (TMS) is a safe, noninvasive, and reversible method for intervening with neural processes in the human brain by a short magnetic pulse. In cognitive neuroscience, the latter method is used to study the specific involvement of brain

* Corresponding author. Department of Psychonomics, Helmholtz Institute, University of Utrecht, Heidelberglaan 2, 3584 CS Utrecht, The Netherlands. Fax: +31-30-253-4511. E-mail address: [email protected] (S.F.W. Neggers). Available online on ScienceDirect (www.sciencedirect.com.) 1053-8119/$ - see front matter D 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2003.12.006

areas in certain types of behavior or actions, by observing changes in the affected behavior shortly after stimulation. A TMS device consists of a coil through which a condensator can discharge a strong current, generating a short-lived magnetic field extending below the coil. The magnetic pulses up to 2.5 T in strength last less than a millisecond and can penetrate the skull and induce electrical currents in the underlying neural tissue that disrupts neural processing for a very brief period (Cowey and Walsh, 1998). At present, two methods for TMS administration are in use. The repetitive TMS (rTMS) method offers a train of pulses for an extended period of time (several minutes), a method which has proven to induce effects in relatively large areas outlasting the duration of actual stimulation up to 15 min (Chen et al., 1997). Testing the participant can reveal functional deficits similar to neuropsychological studies. In single pulse studies, a pulse is administered at a certain time during a behavioral trial, which will transiently induce changes in the brains neurophysiology for several tens of milliseconds (Paus et al., 2001). Single pulse TMS therefore allows one to investigate the time course in which a brain area is contributing to cognitive processes. The spatial resolution is approximately 1 – 2 cm (Bohning et al., 2001), probably considerably smaller than in rTMS studies where the induced electric fields might spread out over a larger area (Lorenzano et al., 2002). For single pulse TMS, the placement of the coil above the proper site of the cortex is critical. In the past, researchers often placed the TMS coil relative to the positions of the international 10 – 20 EEG system, a system based on average brain anatomy. In addition, it is common practice to first locate the motor hotspot on the scalp over the primary motor cortex where TMS can evoke small thumb movements (Paus, 1999). Evoking thumb movements with TMS requires the coil to be placed roughly within 15 mm of the motor hot spot, the size of the region slightly varies between individuals (Boroojerdi et al., 1999; Herwig et al., 2002; Lotze et al., 2003). Relative to this position, the coil is then shifted to a region of interest, whose position relative to the motor hot spot is again obtained from a brain atlas, based on average brain coordinates (Gershon et al., 2003)). Several serious problems arise when using both aforementioned methods, in that one does not know exactly which brain area is affected by the

1806

S.F.W. Neggers et al. / NeuroImage 21 (2004) 1805–1817

administration of TMS pulses. The anatomy of the central nervous system differs substantially between individuals (Luders et al., 2003), and therefore the methods described above only give a rough estimation of the spot being searched for. Therefore, one might not stimulate the actual area of interest. The individual differences in brain anatomy are also a problem in statistically comparing (f)MRI data between subjects or groups, a problem that spurred the development of algorithms that can take these differences into account, map (functional) images to a standard space (Luders et al., 2003; Mazziotta et al., 1991; Woods et al., 1999). Finally, an accurate search for the desired stimulus location might, besides knowledge of individual brain anatomy, also require knowledge of the exact functional representations of the test subject’s brain. In recent decades, new and noninvasive radiological techniques have been developed that are able to measure and visualize the human brain in a volumetric way (an intensity value for each 3D location), such as computer tomography (CT) and magnetic resonance imaging (MRI), providing anatomical information on brain structures. More recently, several functional imaging methods became available that can measure the location and intensity of neurophysiological processes reflecting brain activity, such as functional MRI (fMRI), positron emission tomography (PET), and magnetoencephalography (MEG). Even changes over time can be observed, albeit with varying degrees of temporal resolution (around 1 s for fMRI and PET, tens of milliseconds for MEG). In the present study, we use fMRI data, where the blood oxygenation level dependent (BOLD) signal can be used to indicate regions associated with a certain task or event. In order to use such neuroanatomical and functional information from neuroimaging techniques for coil placement, we constructed a neural navigator (NeNa), a frameless stereotactic device that allows us to use a structural MRI scan of a person’s head to guide a TMS coil to the proper region on the skull. On top of that, it is possible to superimpose fMRI activation maps on the structural brain scan, allowing researchers to target individual functional areas. By doing so, the probability of affecting the intended neuroanatomical structure, and even the intended functions, becomes much higher. The main problem in stereotaxy is to align the 3D space (as measured with a pointer containing a 3D position measurement device, called physical space coordinates throughout this study) of the participant’s head to the MRI coordinates (from now on called image space coordinates) obtained earlier. In other words, one wants to know where the TMS coil is located with respect to the brain structure of the participant, visible on the computer screen. In order to be able to guide the placement of the coil by the MRI image, one has to measure the 3D position of a pointer (used to mark the stimulation target on the participant’s head) in physical space, and transform the 3D coordinates to image space. A few stereotactic methods that are already available use structural and functional data (for examples, see http://www.magstim.com/Brainsight.html), but the specific requirements and circumstances for image-guided TMS experiments using stereotaxy remain to be investigated to implement an optimal solution in a new method. Such a method to calculate the aforementioned transformation is presented in this study. When the transformation is calculated, one can then visualize a 3D representation of the coil together with a 3D representation of the brain surface in image space on the computer screen in real time, that is, one can see a 3D drawing of the pointer ‘move’ with

respect to the brain on the screen so that a location overlying a brain region of interest can be marked on the scalp. In order to be able to calculate the transformation of TMS-coil position to image space, we adopted a method that uses anatomical landmarks that are visible in both the MRI scan and on the surface of the participant’s head, such as the bridge of the nose, the ears, etc. There is no need to place markers on the participant’s head during the scans. The method to calculate the transformation from physical space to image space uses Jacobian coordinate transforms for a first estimate, and a gradient descent procedure for further fine tuning. It is designed to deal with noise properties of the marker positioning and measurements realistic to MRI-guided TMS experiments. The algorithms that were used are described in detail in this paper. The accuracy of the method is estimated based on the position measurement data of six participants. In addition, in a simulation, 25 marker positions were created, with head and image space positions drawn from a similar spatial distribution as the marker positions from a real subject. Noise patterns of increasing magnitude were added to these marker positions, and the robustness of the transformation from head to image space was tested on interpolated positions, as a function of the magnitude of the noise. The transformation was calculated using 15 out of 25 marker positions, and the quality of the transformation was tested on the remaining 10 markers (interpolated positions). The latter reflects the situation in an actual TMS experiment, where usually locations are stimulated that are not included as stereotactic markers. Finally, a validation study of the method was performed, in which the region with increased activation during thumb movements, as measured with fMRI, was targeted with TMS using NeNa. During TMS stimulation, electromyograms (EMG) were recorded on the abductor pollicis brevis (APB) hand muscle, and from the EMG, the motor-evoked potentials (MEPs) were calculated. Motorevoked potentials are short EMG responses starting after approximately 21 ms, and are directly caused by currents induced in the motor cortex that are transmitted to the hand muscles. MEP intensity as a function of coil position, or a cortical MEP map, was constructed to assess the accuracy of the location indicated by NeNa as directly overlying the motor hotspot.

Methods Mapping algorithms The main problem in stereotaxy is to calculate the spatial transformation that is able to express the position and orientation coordinates of a 3D position measurement device, which is held at a certain position with respect to the head of a participant, to the 3D voxel coordinates in an MRI data set of that same person. In order to calculate an estimate of such a transformation, a method is adopted that uses the 3D positions of anatomical landmarks that are visible on a skin rendering of the participants MRI scan as well as measurable with a 3D position measurement device directly at the head of the participant. The problem can be defined by M

XI ! XH þ E

ð1Þ

In Eq. (1), XI and XH denote 3
N matrices containing N rowvectors containing the landmark positions pi (i a [1,N]) in image space I (as selected on the skin-surface in the MRI scan) and qi in

S.F.W. Neggers et al. / NeuroImage 21 (2004) 1805–1817

1807

It is reasonable to assume that with a rigid body transformation (translating + rotating) and scaling (with the same factor for all 3 cardinal directions), the error E can be sufficiently small, since both sets of markers should be located on an isomorphic object. See Fig. 1 for a 2D example. By doing so, the sheer component in the transformation A is assumed to be negligible. Furthermore, it is assumed that the scale factors for the three cardinal directions are identical. Here, we choose A to be a rotation matrix R (expressed in Euler angles of rotation: /, h, and w) in addition to a scale operation (uniform for all directions):

Fig. 1. A 2D example of a rigid body transformation and scale transformation. With this transformation, a shape can be transformed to any other isomorphic shape. T: translation. R: rotation. S: scale. In NeNa, the same transformations are used in 3D.

0

0

ð2Þ

T and Tc are 3
N matrices with the same translation vector t and tc in each row. tc = (tx,ty,tz) is the centroid (average position vector, tc = p¯) of image space marker positions, a fixed parameter. It is included to bring the coordinates XI to the origin, around which the rotations in A operate, whereas t is a free parameter to be estimated by an error minimization procedure (see below). Tc is included only to find a reasonable starting position for the numerical error minimization procedure described below, thereby making it less likely that the latter procedure gets trapped in local minima, and could in principle be omitted.

cos/

sin/

0

1 0

1 0

cosh 0 sinh

C B C B B sinw C C  B0 A @ sinh cosw

s

0 0

1

0

1 C C C C A

0 cosh

1

B C B C B C B C C B0 s 0C B sin/ cos/ 0 B C B C @ A @ A 0 0 1 0 0 s

ð3Þ

There are several definitions for expressing a rotation by three sequential rotations using Euler angles, since a different order of the same rotation yields different results (rotation matrices are not commutative). In Eq. (3), we adopt the NASA standard Aeroplane set of rotations, around the Z, Y, X axes, respectively. The scale factor s is determined by the ratio of the distance between two marker positions in physical space divided by the distance between the same points in image space, averaged over all possible combinations of marker positions:

s¼ XI  T c ¼ AðXH  TÞ þ E

0

B B A ¼ Rð/; h; wÞ  S ¼ B B 0 cosw @ 0 sinw

physical space H (positions of landmarks measured directly at participants head). E denotes the remaining errors that should be minimized by choosing the proper transformation M, mapping the coordinates qi to pi. The latter is not straightforward, since the transformation algorithms should be robust with respect to small deviations in the measured coordinates of anatomical landmarks. Rigid body and scale transformation for coregistering real cranial marker positions with MRI space It is assumed that an affine (linear) transformation M is capable of mapping the landmark positions in physical space to the corresponding positions in image space. The latter implies that nonlinear deformations in both physical and image space can be assumed to be negligible and do not need to be taken into account in the mapping transformation. This is probably safe to assume for MR images that are corrected for geometrical distortions. The affine transformation M in expression (1) can be decomposed into a 3
3 transformation matrix A describing rotation, scaling and shearing, and a separate translation component:

1 0

N 1 X N X Api  pj A 2 N ðN  1Þ i¼1 j¼iþ1 Aqi  qj A

ð4Þ

The latter method of calculating the scaling factor was only included to avoid a calibration procedure of the measurement units used by 3D position measurement device (that can be altered by the experimentator) with the measurements units used to represent space in the MRI (that can vary between images). In principle, the scale factor could of course be calculated and fixed when the scale of the units is known beforehand, to exclude biases in the scale factor caused by marker placement. For the sake of flexibility and ease of use, the scale factor was calculated as in Eq. (4).

Fig. 2. In three steps, the point set qi (i = 1,2,3) in physical space is used to construct an orthogonal set of coordinate axis, as in Eq. (7). The two sets of three points in physical space and image space are both orthogonalized in the same way for all possible combinations of markers.

1808

S.F.W. Neggers et al. / NeuroImage 21 (2004) 1805–1817

Gradient descent minimization: the Powell method The final definition of E is obtained by the substitution of Eq. (3) into Eq. (2): Eðtx ; ty ; tz ; /; h; wÞ ¼ X I  Rð/; h; wÞSðX H  Tðtx ; ty ; tz ÞÞ  T c SSE ¼

N X 3 X

E2ji

ð5Þ

i¼1 j¼1

The remaining six parameters, tx, ty, tz, /, h, w, are free parameters and will be determined by a numerical minimization of the sum of squared errors (SSE). To find a minimum of Eq. (5), a numerical procedure was chosen instead of trying to solve dE/dt = 0, dE/dA = 0 analytically. This was mainly to improve ease of implementation and flexibility, since in the future more complicated mapping transformations might be adopted with more free parameters, possibly including deformation transformations such as shear (that were now left out since it was assumed no deformations of space are present in the acquired marker positions) or even nonlinear deformations that are possibly not solvable in an analytical way. The SSE as a function of six variables will be minimized using a ‘steepest descent’ algorithm that iteratively searches for the direction in which the 6D vector of free parameters has to be altered to maximally decrease SSE (Powell, 1964). Iterative numerical minimization procedures typically could end up with parameters in a local minimum, where the function value does not decrease anymore in any direction of the parameter vector, although a smaller value of the function to be minimized exists. Typically, one wants to find the global minimum of a function, that is, the set of free parameters yielding the smallest possible function value. One minimizes the risk of finding a local minimum, by starting the search at a set of six parameters that lies as close as possible to the global minimum. Therefore, one of the main problems that has to be encountered when using numerical minimization procedures is finding an informed guess of a starting position that is believed to be not too far from the global minimum.

Fig. 3. The NeNa rendering and navigation window, with a surface of the skin, brain, and activated regions. Also, the six markers used to coregister physical space to image space are plotted as small spheres.

Fig. 4. With the tip of this probe (T), the landmark positions on the participants skin can be measured. Inside the probe, a MiniBIRD marker is mounted, of which the position and orientation can be read every 10 ms. Since the dimensions of the probe are known, the position of the tip can be calculated.

Rough estimate of rotation based on Jacobian coordinate transformation of paired positions For the minimization problem at hand, a method is developed that roughly coregisters the set of marker positions in physical space to those in image space, that is, finds a starting vector (tx, ty, tz, /, h, w) with values for the translation and rotation components that already yield relatively small values for Eq. (4). To do so, first the fixed translation matrix Tc (a matrix with the vector tc in all rows) and the scaling factor s are calculated as in Eqs. (5) and (4), respectively. A good estimate for the translation T would also bring the marker positions in physical space around the origin, as Tc does for the markers in image space. Therefore, the centroid (average position) of all markers in physical space is therefore taken as a starting position for the translation parameters in Eq. (5): t = (tx,ty,tx) = q¯. Finding a good approximation of the Euler angles /, h, and w, determining rotation matrix R that minimizes Eq. (5) is not as straightforward and discussed in the following. Here we are first deriving an expression for the rotation matrix R, from which the Euler angles can be calculated. Consider XHV , the translated (to the origin using an estimate of t as discussed above) and scaled (as in Eq. (4)) set of physical space markers XH from the right-hand side of Eq. (5) and the set of image

Fig. 5. The positions of the 5
5 grid (green lines, intersections were stimulation positions) that were stimulated during the validation experiment, visualized on a skin rendering of one of the subjects (DB). The T-map representing neural activity is also rendered (in red).

S.F.W. Neggers et al. / NeuroImage 21 (2004) 1805–1817

markers X IV from the left hand side of Eq. (5) also translated to the origin: XHV ¼ SðXH  Tðtx ; ty ; tz ÞÞ XIV ¼ X I  T c

ð6Þ

XHV only needs to get rotated around its centroid (now lying at the origin) to roughly fit the image space markers X IV , since

1809

both sets have the same scale (distance between marker positions) and have their centroids at the same position (origin). An exact match is not possible, since both sets have some noise (position errors due to measurements). To determine the appropriate estimate of the rotation angles /, h, and w, the following approach is adopted. First, two local coordinate frames F and G spanned by the orthogonal vectors fx, fy, fz and gx, gy, gz are defined. F is constructed from two markers in translated image

Fig. 8. The results for each participant are plotted in rows. The columns, from left to right, represent the following results. (A) T-maps of bR – bL contrast (see Methods), or right-thumb movement vs. left-thumb movement, overlaid on a structural T1 image (coronal view). (B) NeNa visualization of the same data, with renderings of the skin, cortical surface, and T-maps. The cross (
) denotes the position T, the location with the highest T value targeted by NeNa and used to center the 5
5 stimulation grid. (C) Pseudocolor map of MEP amplitude as a function of stimulated position, relative to the point T (0,0) that was indicated by NeNa. Red indicates high and blue indicates low amplitudes. The cross (
) indicates the position of the center of gravity (CoG) as defined in Eq. (10). The color bars give the exact EMG amplitudes in millivolt corresponding to a color in the map.

1810

S.F.W. Neggers et al. / NeuroImage 21 (2004) 1805–1817

F and G, the coordinate transformation JF!G takes the form of the following equation: 0

gx  f x

B B JF!G ¼ ðgx ; gy ; gz ÞT ðfx ; fy ; fz Þ ¼ B B gy  f x @ gz  f x ¼ RF!G ðGÞ

gx  f y gy  f y gz  f x

gx  f z

1

C C gy  f z C C A gz  f z ð8Þ

J is also the required rotation matrix R rotating the coordinate axis G to F; note that this rotation matrix in Eq. (8) is expressed relative to subspace G. Finally, the rotation matrix R( G) in Eq. (8) has to be expressed in global space S in which the coordinates p and q are stored. R(S) can be derived from the Jacobian transformations between global space S and G, JS!G, and JF!G: T RðSÞ ¼ JS!G RðGÞJS!G

Fig. 6. The average error per marker, CRE (=M(SSE)/6) from the first estimate of coregistration until after the minimization procedure is plotted, as a function of the number of iterations, for six participants. Note that the minimization slightly reduces the average errors with respect to the first estimates. Except for one participant, all resulting errors are less then 5 mm.

space XIV, say piV and pjV, and the centroid in image space p¯ (now at the origin), and G from two marker positions representing the same landmarks in XHV , say qiV and q Vj , and the centroid position, now at the origin (q¯ = 0). The indices i and j could be any combination of markers. Space G is constructed as follows: gx ¼

qiV AqiVA

gy ¼

qiV
qjV AqiV
qjVA

gz ¼ gy
gx

For all possible combinations of markers (indices i and j choosing piV, pjV and qi, qj from which the local coordinate systems described above were constructed), the matrices R1(S) . . . RN(S) were calculated, and then the component-wise average matrix Mmn P ðSÞ ¼ N1 Ni¼1 Rimn was constructed by averaging the matrix components Rnm. Since M is not a pure rotation matrix anymore (for pure rotation matrices det(R) = 1 and RRT = I), we use singular value decomposition (SVD) to find the pure rotation matrix closest to M that can be expected to be close to the average rotation matrix ¯ . SVD can be used to decompose an affine matrix M into two R unitary matrices U and V, and a ‘stretching’ matrix S with only numbers on the diagonal, such that M = USVT. The ‘average’ ¯ (with det(R ¯ ) = 0 and R ¯R ¯ T = I) close to M can rotation matrix R ¯ = UVT (eliminating the stretching component than be derived as R of M). One cannot, as in coordinates with an Euclidian metric,

ð7Þ

see Fig. 2 for a graphical representation of the construction of the coordinate system. In the same way, the local coordinate system F is constructed from pi and pj. Mathematically, this comes down to Gram – Schmidt orthogonalization for each of the normalized set of vectors (piV, pjV, piV
pjV) and (qiV, q jV , q iV
qjV) to obtain an orthogonal basis. Note that G and F share the same origin, and hence F can be rotated to exactly match G. Roughly, the rotation matrix bringing G to F is close to the rotation matrix bringing the set of markers XHV to XIV , since the local coordinate systems G and F are constructed from the same set of markers in XHV to XIV and hence F has the same relative orientation with respect to space I as G has with respect to H. It is probably not the rotation matrix that optimally matches XHV to XIV since the markers pi, pj, qi, and qj that were used to construct the subspaces F and G contain some errors as well. The matrix transforming coordinates of a position in space expressed in the axes of F to the coordinates of that same position expressed in G is the Jacobian transformation matrix JF!G. For two systems of orthogonal normalized coordinate axes

Fig. 7. The performance (errors) of the coregistration algorithm on the 10 markers in the test set (TRE) as a function of noise, of which the image space positions are known but not used to determine the coregistration transformations. The performance on these interpolated (i.e., not used for coregistration) positions under noisy conditions is representative for experimental situations, in which usually positions are targeted with TMS (at the scalp) that are not used for coregistration.

S.F.W. Neggers et al. / NeuroImage 21 (2004) 1805–1817

directly average rotation matrices, due to their Riemannian metric (for more details, see Moahker, 2002). Finally, the Euler angles (according to the NASA standard Airplane definition) /, h, w as in Eq. (3) corresponding to the ¯ can be directly derived from R ¯: ‘average’ rotation matrix R ¯ 11 Þ ¯ 12 =R / ¼ arctanðR ¯ 13 Þ h ¼ arcsinðR ¯ 33 Þ ¯ 23 =R w ¼ arctanðR

ð9Þ

The ‘average’ rotation matrix R¯ can be expected to be close to the rotation bringing XHV to XIV, since the rotations between the ¯ correctly rotate local coordinate system used above to create R subsets of markers from space H to I, and can therefore serve as a starting point for /, h, w for the minimization procedure from Eq. (5). Visualization techniques The NeNa software in its present form uses the Visualization Toolkit C++ library VTK (http://public.kitware.com/VTK/) to create the surface renderings of skin, brain, and brain activity. Preprocessing of MRI volumes In the present study, representations (surface rendering) of the skin, brain, and functional maps were needed (see Fig. 3). To construct these surfaces, a T1 weighted, segmented gray matter image and a statistical T-contrast (T-map) derived from T2* weighted time series images were used. In order to generate the skin rendering from a T1 structural image or a statistical T-map, the data does not need to be preprocessed; only a threshold value has to be entered (see next paragraph). The rendering of the brain surface requires a segmented gray matter image that can be generated manually or with many MRI data analysis packages, such as statistical parametrical mapping (SPM). Triangulation of skin, brain, and T-map surface based on intensity threshold: marching cubes The surfaces of the skin, brain, and T-maps (see Fig. 3) are generated using the marching cubes algorithm (Lorensen and Cline, 1987). The marching cubes algorithm produces a polygonal model from the voxel data by traversing all the voxels and using the values of their corner points to replace them by polygons. The placement of these polygons depends on the value of an intensity threshold. The value of the intensity threshold depends on the image modality that was used and can be changed interactively. In the present study, the skin rendering was created from a T1 weighted image using a threshold between the value for voxels outside the head and skin voxel values. For T-map rendering, a statistical threshold was chosen, corresponding to T = 4. For segmented brain image renderings, a threshold slightly below that of gray matter voxel values was used. After subjecting the voxel data to the marching cubes algorithm, it needs to be post processed in two ways. First, since the marching cubes algorithm creates a very inefficient polygonal model, it needs to be triangulated to enable a faster rendering.

1811

Second, the marching cubes algorithm may produce unconnected parts in the polygonal model, which can be deleted. Marker positioning and navigation In the main rendering window of the NeNa software package, the marker positions used for the coregistration procedure can be chosen and moved around along the X, Y, and Z coordinates to the proper locations on the skin rendering. In Fig. 3, the six markers that were used in the present study can be seen. Many more markers could be entered when necessary, as long as they are visible on both a skin rendering derived from a structural image, as well as directly at the participant’s head in the experimental setup. From experience, it appears that the best results for marker placement (and hence coregistration) can usually be obtained when the same person selects the marker positions on the skin rendering as well as at the head of the participant in the setup. In order to facilitate navigation, the image in the rendering window can be viewed from all directions and distances, which can be done in real time with adequate graphics hardware. Furthermore, it is also possible to view 2D slices from the structural image volume, in three separate sagittal, coronal, and transversal projections. This is particularly useful for markers that lie in notches in the surface, which are difficult to locate precisely in a three-dimensional surface rendering. Finally, the position and orientation of the MiniBIRD probe used to measure marker positions and target areas of interest are visualized in the main rendering window after all positions are measured and mapping procedure has been calculated. It can be seen hovering over the head with a thin line extending into the head with the same orientation as the real probe has with respect to the head. This line can be used to target a certain region of interest. Even more useful is the feature to align the camera-view (from which perspective the images are rendered) with that line, allowing a close-up top view of the target region and/or activity maps as if one were looking through the MiniBIRD probe. The center of the screen (which is indicated by a circle) can then be used to target a region; when a certain area is visible through the circle, it is directly pointed at by the MiniBIRD probe. The latter procedure was used in this study. Cranial marker position measurement The anatomical marker positions used to perform the coregistration of a participants head with his corresponding structural MRI data have to be clearly visible in a comparable way both in the skin rendering obtained from the MRI data and directly at the participant in the experimental setup. In the current experiments, six anatomical landmarks are used: the crest of helix of the left and right ear, the bridge and the tip of the nose, and the medial commissure of the left and right eye. The latter locations are often used in medicine for image-guided neurosurgery relying on comparable coregistration with imaging data and have proven to yield reliable coregistration. To improve coregistration accuracy, one could place much more than six fiducial markers on the head during the MRI scan, and use them as visible markers in a TMS experiment directly afterwards. Acquisition of probe tip position In order to measure the exact location of a position on the participants skin, a small probe has been constructed with a 2-mmwide spherical tip (see Fig. 4). Inside the probe, a MiniBIRD marker is located. The MiniBIRD system of Acension Technolo-

1812

S.F.W. Neggers et al. / NeuroImage 21 (2004) 1805–1817

gies is capable of measuring 3D position (x, y, z) coordinates with about 1-mm precision, as well as the 3D orientation angles, of a small marker in real time (110 Hz). The MiniBIRD position tracker consists of an electronics unit, a transmitter, and a receiver (the marker inside the probe). The transmitter generates a pulsed DC magnetic field. In a pulsed DC magnetic field, stainless steel and other non-ferromagnetic metal objects do not interfere much with the field strength and do not disturb the position measurements. The changes in magnetic field strength induce a current in the MiniBIRD marker, which is depending on the orientation and distance from the transmitter. This current is interpreted by the electronics unit, and the position and orientation measurements are sent directly to the NeNa software. Since the dimensions of the probe have been calibrated, it is possible to calculate the position of the 2-mm wide tip in real time from the position and orientation of the MiniBIRD marker. Validation experiment: data acquisition, design, and analysis Procedure For a typical NeNa-guided TMS study, a separate fMRI study needs to be conducted first. For the motor mapping study presented here, the entire procedure was as follows (see below for more details on each step): Four participants made thumb opposition movements in an fMRI scanner with their left or right thumb in alternating blocks, during which the blood oxygenation in their brain was measured. A T1 weighted structural MRI scan was obtained as well. After the fMRI experiment, the statistical parametrical maps (SPMs) reflecting thumb movement activity in the motor cortex were constructed (see below). The same participants were measured in a TMS experiment that took place several days after the fMRI experiment. The structural scan as well as the statistical parametrical map (SPM) reflecting the T-contrast (T-map) between right and left thumb movements (visible as a ‘blob’ of activity on the motor cortex, see visualization) were loaded into NeNa before the experiment, and the cranial markers were set on the skin rendering. The participant would then be seated, the head fixed in a head support, and a rubber head cap placed over his head. The positions of the same cranial landmarks were measured with the MiniBIRD directly at the head of the participant, and the mapping described in Mapping algorithms section was calculated. With the NeNa probe, the region of increased activity obtained from the SPM was targeted by ‘looking’ through the probe to the renderings of the skin and activity blob on the screen (see Marker positioning and navigation section), and the position directly overlying the maximum value in the activity blob (the T-threshold was set so high that only a maximum activity map of about 4 mm remained) was marked on the head-cap. Around this position, a rectangular grid of 5
5 positions, 5 mm apart, was drawn (see Fig. 5). TMS was applied over all 25 positions, during which electromyograms (EMGs) were recorded from the thumb muscle abduction pollicis brevis (ABP) contralateral to the stimulated site. When stimulating the motor cortex with TMS, so-called motorevoked potentials (MEP) are known to be visible in the EMG signal measured at the contralateral thumb, that is, peaks of short duration (about 10 ms) at a fixed latency after pulse administration (about 22 ms), directly caused by the neural signal that was electrically induced by the TMS pulse in the motor cortex and conducted to the thumb muscle (Lorenzano et al., 2002; Ray et al., 2002). Finally, a map of the magnitude of the MEPs was constructed and the center of gravity (MEP maximum) was calculated.

When the coregistration algorithm described here worked perfectly, and a perfect functional match would exist between the MEP maximum and the fMRI region of maximum intensity, the MEP maximum should be in the fMRI region of maximum intensity (see Fig. 8 where the MEP maximum is located at coordinate 0). fMRI thumb movement experiment: acquisition and analysis Functional imaging during the thumb movement task was performed with a Philips ACSNT 1.5-T clinical scanner, using the blood oxygen level dependent (BOLD) sensitive, navigated 3D PRESTO pulse sequence (Ramsey et al., 1998; van Gelderen et al., 1995), with the following parameter settings: TE/TR, 35/24 ms; flip angle, 9j; FOV, 256
192
96 mm; voxel size, 4 mm isotropic; scan time per fMRI volume, 1.49 s; 345 scans per session. After the functional sequence, an anatomical scan was acquired (TE/TR, 4.6/30 ms; flip angle, 30j, FOV, 256
180
208 mm; voxel size, 1
1
1.2 mm). Importantly, in the 3DPRESTO technique as well as the other 3D MR pulse sequences, other than in an EPI sequence, the TR does not reflect scan duration (Liu et al., 1993; van Gelderen et al., 1995). The PRESTO pulse sequence uses crusher gradients to enable fast acquisition, and a navigator echo for movement correction (Ramsey et al., 1998). The crusher gradients destroy much of the signal in larger draining veins, thereby reducing the emergence of ‘activity’ in veins. The three-dimensional acquisition with a low flip angle ensures that there are also no ‘activations’ caused by blood flowing into the scanned volume (Duyn et al., 1996). All these properties of PRESTO improve the localization of activity in neuronal tissue, which has been validated with positron emission tomography (Ramsey et al., 1996) and with intraoperative electrical stimulation (Rutten et al., 1999). Another advantage of PRESTO is the fact that spatial distortion is negligible, further improving colocalization of activity with the anatomical image (Ramsey, 1999). Four healthy participants without a history of neurological treatment were included in the experiment. All were informed beforehand on the experimental procedure and gave their written informed consent. The Medical Ethical Committee (METC) of the Utrecht University Medical center approved the experiments, which were in accordance with the TMS safety guidelines. During functional imaging, the participants looked at a screen visible through an overhead mirror while lying in the MRI scanner. A beamer, controlled by a laptop PC running the Presentation visual stimulation software package from Neurobehavioral Systems, projected the word LEFT or RIGHT on the screen, and after 23 s the word STOP was presented. After 20 s of rest, the next word LEFT or RIGHT was projected on the screen, etc. Each word was presented four times in a random sequence. Also, three extra rest blocks of 23 s were included. Participants were instructed to move their left or right thumb, respectively, at a pace of about 2 Hz, while their arms were lying stretched and relaxed besides their bodies. Data was analyzed with SPM99, a software package that is developed and freely distributed by the Wellcome Trust in London, UK, using MATLAB 6.5. During preprocessing, the functional time series were coregistered to the high-resolution anatomical scan, and corrected for movement in the scanner by realigning them to the first scan using a least square minimization. Also, the functional data was smoothed with a kernel of 8 mm to increase the signal to noise ratio. The loss in localization accuracy is supposedly small with respect to the expected size of the area activated by thumb opposition movements (in the order of 1 – 2 cm). From here on,

S.F.W. Neggers et al. / NeuroImage 21 (2004) 1805–1817

two different sets of functional time series data were generated during spatial preprocessing that were then both analyzed the same way. In the first set, the data was not normalized to standard space (MNI) since the resulting activity maps will be used for TMS stimulation using stereotaxy, requiring the data to be in native space (not deformed). In the second set, the functional time series data were normalized to MNI space based on the T1 structural images of each subject (using nonlinear warping techniques available in SPM99) to classify the activated areas (corresponding to the activated areas from the unnormalized set that were targeted with TMS) in a standard stereotactic space. SPM uses a General Linear Model (GLM) approach, assessing brain activation in a statistical parametrical map (SPM) to a series of regressors using a least square linear fit method (Friston et al., 1995). Per subject, a GLM was fitted to each voxel using three regressors: Y¯ ¼ X¯¯ b¯ þ e¯ where Y is the data for a voxel, X = {XL, XR, XI, XC} is the 3
345 design matrix with three regressors for left- and right-thumb movements: the moment of written instruction on the screen and a constant term, b = {bL, bR, bI, bC}. e is the remaining error variance that could not be explained by the model. A set of cosine regressors was added for high-pass filtering the data with a cutoff period of 150 s (wavelength). A box-car regressor was used for left- and rightthumb movements separately, where 1 = moving and 0 = rest. One regressor was added (effect of no interest) for the instruction (event related): delta function at moment of word presentation convolved with canonical hemodynamic response function. The latter regressor serves to regress out any activity related to the presentation and interpretation of the instruction written on the screen. The b’s associated with each regressor were tested with 2 one-tailed t tests for the following hypotheses: bL – bR > 0 and bR – bL > 0, reflecting voxels that are more active for left than for right thumb movements and vice versa. Also, by calculating the difference between activity associated with right and left movements instead of directly analyzing the activity associated with one effector (bR), the resulting areas can expected to be more specific for the movement of that specific effector, and not for additional processes related to making a movement in general. The reported regions are corrected for multiple comparisons using the random field correction method (Worsley et al., 1992, 1996). See Fig. 8 in the Results for the T-maps for each participant. Measuring motor-evoked potentials in thumb abductor during TMS The T-maps constructed from the fMRI data described above were loaded into the NeNa software and thresholded at T = 4.27 ( P < 0.05 corrected) for all four subjects. The rendered surface of this activation is then targeted with the MiniBIRD probe by aligning the viewing direction with the probe axis and navigating with the probe (during which it is held perpendicular to the head surface) until the target blob is directly in the center of the field of view, on the NeNa crosshair. The position T is defined as the position where the central axis of the MiniBIRD probe intersected the surface of the head (see Fig. 5 for an example), when directed at the center of the rendered surface representing maximum activation. Subjects were wearing rubber caps on which T was marked. A matrix with a grid of 5
5 positions with 5 mm between the positions was drawn on the cap as well, aligned with the midsagittal plane of the head (see Fig. 5). During the TMS experiment, first the motor threshold was deter-

1813

mined for that particular subject as the lowest intensity that induced visible finger movements in 50% of the trials when TMS was applied to T (Pridmore et al., 1998). The intensity of TMS stimulation was then set at 110% of the aforementioned individual motor threshold throughout the experiment. A Neopulse TMS device (Neotonus Inc, Atlanta) with an iron core coil was used (Epstein and Davey, 2002), in contrast to many other TMS studies that used a figure of 8 coil without a ferromagnetic coil. The iron core is embedded at the center of a rectangular figure of 8 coil, parallel to the wiring at the center. Ferromagnetic cores cause the generated magnetic field to be stronger and to penetrate deeper into the brain (Epstein and Davey, 2002). In addition, iron core coils generate more focal magnetic fields than coils without such cores. The induced currents in the brain form a figure 8, whose transition point for current direction is constrained by the outer edges of the iron core. This makes for a somewhat smaller boundary to the induced field than a figure 8 coil without an iron core. Summarized, the core acts to project the field deeper into the brain, with the net effect that the field is more focused than for typical figure 8 coils (personal communications with Dr. Epstein). During the actual experiment, 10 TMS pulses were administered at approximately 0.25 Hz, at each of these 25 positions, yielding 250 pulses per subject. A random delay of 500 ms was introduced before each pulse to prevent preparation effects or movements produced by the participant. The TMS device was triggered by a PC using the Presentation stimulus software package (Neurobehavioral Systems), by a 5 V TTL pulse over the parallel (LPT) port. During stimulation, the EMG potentials were measured at the APB using the PsyLab EMG acquisition system, sampled at 1000 Hz, 3 kHz lowpass, and 0.3 Hz highpass filtering and stored on disk for offline analyses (MEP map construction). Motor-evoked potentials and center of gravity The EMG signals stored on disk are used to make maps of motor-evoked potentials MEPi at i = 5
5 = 25 positions pi for each subject. Per stimulated location, the 10 raw data traces (a vector with 1000 Hz samples) were prepared as follows. First, all traces were aligned at t = 0, the time of pulse administration. Then for each location, an average signal was calculated from the 10 traces that were measured at that location. The difference between the maximum and the minimum value, or peak-to-peak distance, between 20 and 35 ms after TMS pulse administration was taken as the amplitude of the MEP, since according to literature the main MEP response can be expected in this time interval (Lorenzano et al., 2002; Ray et al., 2002). The averaging and peak detection was done in home made Matlab 6.5 scripts. The center of gravity (CoG) of a map of MEPs, or the average of stimulated position coordinates weighted by MEP intensity, is known to yield a more accurate estimate of the location on the scalp directly overlying the region of maximum excitability of thumb potentials than simply taking the location from the 5
5 matrix where the MEP is maximal (Boroojerdi et al., 1999). For each participant, the CoG vector was determined from the MEP matrix according to Eq. (10): N X

CoG ¼

pi MEPi

i¼1 N X i¼1

ð10Þ MEPi

1814

S.F.W. Neggers et al. / NeuroImage 21 (2004) 1805–1817

The latter definition of CoG had been used in other studies as well (Boroojerdi et al., 1999; Lotze et al., 2003), and is therefore suitable to compare the present method to other studies.

Results Coregistration The coregistration procedure described above was performed for six participants, of which four also participated in the fMRI study. First, the six cranial landmarks pi, described in Cranial marker position measurement section, were set in the MRI data using NeNa. The same six landmark positions qi were measured directly at the participant. First, the matrices S, R, and translation T were estimated as described in the Methods. The scaling factor s was calculated as in Eq. (4) from the relative distances between corresponding markers for the sake of flexibility, and not fitted using a gradient descent method as was done for the other parameters. The values of s were fairly constant over the six subjects (0.874 F 0.032 mm/voxel), indicating that no substantial bias was introduced at this level (see Methods). Finally, during the gradient descent process, the free parameters (tx, ty, tz, /, h, w) startingpwith ffiffiffiffiffiffiffiffithe informed guess values described above, were monitored. SSE =N , the estimated average error per marker or correspondence registration error (CRE) as in Eq. (4) was plotted against time (see Fig. 6) for all six participants, and for all except one participant, the SSE nicely converged to a value less than 5 mm. The starting estimate of the free parameters is already below 8 mm, demonstrating that the informed guess of the starting parameters is already close to the optimal value, and hence the risk of ending up in a local minimum is small. The coregistration process as described in the Methods was also tested in a simulation using 25 simulated marker positions, resembling the marker position coordinates in physical space and image space from a real subject. The 25 simulated marker positions were formed by 25 sets of image and physical space coordinates pi and qi, with a known spatial relationship (translation, scale, and rotation) between them, and constructed as follows. For one of the participants above, the free parameters (tx, ty, tz, /, h, w) and the fixed parameter s were calculated from the six marker positions pi and qi that were measured in image and in physical space as described in the Methods. In physical space, the average R¯ and standard deviation r of the distances Ri between the centroid position q¯ and all qi’s were calculated for these six markers. From a normal distribution with the same mean and standard deviation, 25 samples of radii Ri were drawn. Also, 25 corresponding sets of angles /i and hi were randomly drawn from a uniform distribution (p < / < p and p/2 < h < p/2). These 25 sets of points qi with coordinates /, h, and Ri in a spherical coordinate system are thus centered around q¯. They therefore occupied roughly the same space as the original six measured positions in physical space. Finally, the coordinates were expressed in Cartesian coordinates x, y, and z. With the rotation and translation matrices corresponding to the fitted parameters tx, ty, tz, /, h, w of the subject on which the distribution was based from which the 25 coordinates were selected, all the 25 points qi were transformed to the simulated ‘image space’ coordinates pi. After this, some random noise was added to the marker positions in image space,

with different magnitudes. To the x, y, and z coordinates, a number L was added independently, where L was drawn from a uniform random distribution with mean k R¯, k equals (0, 0.01, 0.02, . . .. 0.14). R¯ now equals the average distance of the image space markers to the centroid, here 61.1 mm. To estimate the spatial dimension of the noise for each k, this procedure was repeated 10 000 times for each noise level k, and the amplitude of the noise added to each marker was averaged (see Fig. 7). The optimal spatial transformation between both point sets is known for these 25 sets of markers, and it will be investigated how well NeNa is capable of finding it, under different noise conditions. First, the set (pi, qi) of 25 markers is divided into some so-called training set of 15 markers, and test sets of 10 markers. All possible combinations of individual markers over the two sets are being used. The mapping algorithm as described in the Methods section will be performed using the 15 markers in a training set, under one of the noise conditions. Then, the physical space coordinates qi of the 10 markers in the corresponding test set that were not used to determine the coregistration transformation were transformed to image space using the coregistration matrices mentioned above (based on the 15 training set markers), and the average distance in millimeter to their image space counterparts, the so-called target registration error (TRE), will be calculated for that noise condition. This will be repeated for all the test and training sets (with all combinations of the 25 markers) obtained, and the TRE will be averaged. The TRE reflects the error when targeting a region of interest in a real TMS experiment, since usually regions will be targeted (in the brain) that were not used as markers in the coregistration. The average TRE in millimeter as a function of the noise parameter is plotted in Fig. 7. Clearly, the coregistration found by the algorithm is perfect (TRE is zero) when no noise is added. The TRE for noisy conditions roughly increases linearly with noise, with a slope smaller than 1, that is, the TRE is somewhat smaller than the spatial extent of the noise that was added to the complete set. Fig. 7 is reflecting a robust mapping under realistic noise conditions. fMRI thumb tapping experiment The T-maps for contrast bR – bL are plotted on a structural T1 scan (unnormalized data set) in Fig. 8, left column, for all four participants in the fMRI study. In the middle column, the skin, cortex, and T-map rendering from the NeNa navigation window are plotted, the left motor strip activation is clearly visible. Also, the position with the highest T-value is indicated with a cross (
). This position was used to target the focus of the activated area with TMS (with the viewing direction aligned with the pointer) and mark the overlying point T on the scalp. The corresponding MNI coordinates of the maximum T value for the contrast bR – bL activation for all four subjects (normalized data set) are given in Table 1, thresholded at P < 0.05 corrected for multiple comparisons using a random field correction (Worsley et al., 1992, 1996), together with the associated Brodmann areas. As expected, clear activation can be seen around the left motor cortex (Brodmann area 4). MEPs For each of the 25 positions in the grid around the position T marked by the NeNa pointer, 10 traces of EMG were acquired. The signal was then preprocessed as described in the Methods to

S.F.W. Neggers et al. / NeuroImage 21 (2004) 1805–1817

1815

Table 1 For each participant, this table lists (from left to right) an anonymous subject ID, the MNI coordinates in millimeter of the maximum fMRI activation (T value) for right-thumb movements, the corresponding T and Z values, the Brodmann area in which the latter maximum was located, the location of the center of gravity (CoG) of the MEP map in millimeter, and the gridpoint location in millimeter yielding the largest MEP response after TMS stimulation Participant

MNI coord

DB EH MN BN

32, 48, 40, 40,

28, 20, 20, 32,

68 60 64 72

T (DF)

Z

Area

CoG (x,y) (mm)

Max. loc. (mm)

6.38 6.29 5.44 6.63

6.19 6.10 5.32 6.42

4 3/4 4 4

(1.50, 1.25) (6.38, 5.67) (3.17, 2.20) (0.17, 2.18)

(10,0) (10,10) (0,0) (5,5)

(330) (330) (330) (330)

calculate an average per position. In Fig. 9, the voltage of 10 EMG recordings is plotted against time, aligned at the time of TMS pulse administration (t¯ = 0), for a stimulation position that yielded significant MEPs. In addition, the estimated MEP amplitude for that stimulation position is indicated in Fig. 9. Other stimulation positions resulting in significant MEP responses from the 5
5 grid yield EMG signals that are similar to those shown in Fig. 9, for all four subjects (although the absolute numbers of measured voltages differ). A clear first peak can be seen after 23 ms, representing the MEP directly caused by TMS. The MEP amplitude of the first peak was calculated as described in the methods for each stimulated position. The magnitude of the MEP response as a function of the stimulated position is displayed as color-coded maps in the right panels of Fig. 8 (next to the fMRI results), for each subject that participated in the fMRI study. A cross (
) denotes the CoG as defined in Eq. (10). The CoGs of each subject are also given in Table 1. The average CoG was 1.13 and 2.20 mm, and the average distance of the CoG to the point T was 4.14 mm. The location of the grid-point yielding a maximum MEP is also given in Table 1. On average, the maximum MEP tends to be evoked when stimulating a region slightly posterior of the position T indicated to overly the fMRI maximum.

Discussion The present study proposed a neuro-navigational method using point-based coregistration. Based on an fMRI study that was able

to activate the thumb-movement area in Brodmann area 4, four participants were stimulated with TMS using the method proposed here to place the TMS coil. The validation study using imageguided TMS and simultaneous EMG recordings demonstrated that the center of gravity of the MEPs was always within 4 mm of the position T marked by the neuro-navigator to be right above the maximum fMRI activation, except for subject EH where the CoG was located 8.5 mm from the position T. In the central panel of Fig. 8, one can see that this might indicate that the region in the statistical map of fMRI activation with the highest T value is not always the best indicator for maximum EMG responses, since relative to the target used for subject EH, the active area extends posterior-medially, just as the location of his CoG. For EH, the location yielding maximum MEP tends to be located up to 10-mm posterior to the fMRI maximum (see Table 1). Interestingly, studies by Herwig et al. (2002) and Lotze et al. (2003) that used a similar technique observed CoGs 18.8 and 10 mm anterior to fMRI activation for thumb movements, respectively. It was speculated (Lotze et al., 2003) that this might have been due to the figure of 8 TMS coil that was used in both aforementioned studies. In figure of 8 coils, the location of neural tissue affected by TMS pulses has been shown to depend on coil orientation (around a line perpendicular to the surface of the figure 8 coil) and design (Kammer et al., 2001). In our study, an iron core TMS coil was used, where the field directly underneath the center of the coil is stronger, and penetrates deeper into the brain (Epstein and Davey, 2002), which may explain the differences (depending on how the coil was held in the mentioned studies). Another reason for the different results might be that the

Fig. 9. Single EMG traces (thin gray lines) and their average (thick black line), as a function of time since TMS administration, at a stimulated location that yielded significant MEPs. The vertical dashed lines denote the time of maximum and minimum EMG values, as determined by the algorithm described in the Methods. The traces are representative for all recordings during which MEPs were observed. The difference between the maximum and minimum value in the average signal is the MEP amplitude in the maps in Fig. 8 (column C).

1816

S.F.W. Neggers et al. / NeuroImage 21 (2004) 1805–1817

two studies mentioned before, which also investigated spatial congruency of MEP maximum and the maximum BOLD signal as measured with fMRI, both used the Echo Planar Imaging (EPI) technique to measure the BOLD signal, whereas the present study uses the PRESTO technique, which is less sensitive to signal change caused by inflowing as well as draining veins, which can shift the BOLD contrast away from the site of real neural activation (Ramsey et al., 1996, 1998). On another note, several studies have reported that fMRI-based localization of brain function does not map perfectly onto localization based on direct cortical stimulation during surgical procedures in epilepsy patients (Rutten et al., 2002). This mismatch (in the order of 5 mm) may be because fMRI detects activity in neurons that are involved but not critical for motor function (as opposed to TMS), or to complexity of the source of fMRI signal change (Attwell and Iadecola, 2002). Furthermore, simulations with coregistrations based on 15 marker positions revealed that under realistic noise conditions, an estimated error of 3 – 4 mm could be reached in theory, depending on the noise in marker placement. The magnitude of these simulated errors was indeed close to the observed CoG accuracy of 4.14 mm. Realistic noise is most likely to be generated by slight errors in skin surface renderings (on which the 3D markers were placed on the screen), measurement inaccuracies of the magnetic position measurement device, and errors generated by the investigator when placing the markers. When one has an estimate of the noise in measuring marker positions, one can get an idea of the expected accuracy of TMS coil placement using NeNa by taking the corresponding value in Fig. 7 obtained from the noise simulations. For TMS experiments, the accuracy of 4.14 mm obtained from the CoGs in the present validation experiment is a promising result. In single pulse TMS, the neural tissue affected by stimulation is most likely 1 – 2 cm in diameter (Bohning et al., 2001), and the intended targeted area should therefore be affected when the present method is used to determine the locus of the TMS stimulation. In summary, the device for fMRI-guided positioning of a TMS coil yields a promising degree of spatial accuracy. The device is inexpensive and relatively easy to apply. In combination with fMRI-derived brain activity maps, this device enables investigators to manipulate brain regions accurately for individual subjects.

Acknowledgments The authors thank Erno Hermans and and Matthijs Noordzij for assistance in the EMG acquisition and TMS pulse administration, and Dr. Epstein for helpful comments on TMS coil design. Furthermore, we would like to thank two anonymous reviewers for very helpful comments on an earlier version of this manuscript. This work was supported by a ‘Pionier’ grant from the Netherlands Foundation for Scientific Research (NWO, nr: 440-20-00).

References Attwell, D., Iadecola, C., 2002. The neural basis of functional brain imaging signals. Trends Neurosci. 25 (12), 621 – 625. Bohning, D.E., He, L., George, M.S., Epstein, C.M., 2001. Deconvolution of transcranial magnetic stimulation (TMS) maps. J. Neural Transm. 108 (1), 35 – 52.

Boroojerdi, B., Foltys, H., Krings, T., Spetzger, U., Thron, A., Topper, R., 1999. Localization of the motor hand area using transcranial magnetic stimulation and functional magnetic resonance imaging. Clin. Neurophysiol. 110 (4), 699 – 704. Chen, R., Classen, J., Gerloff, C., Celnik, P., Wassermann, E.M., Hallett, M., Cohen, L.G., 1997. Depression of motor cortex excitability by low-frequency transcranial magnetic stimulation. Neurology 48 (5), 1398 – 1403. Cowey, A., Walsh, V., 1998. Magnetic stimulation studies of visual cognition. Trends Cogn. Sci. 2 (3), 103 – 110. Duyn, J.H., Yang, Y., Frank, J.A., Mattay, V.S., Hou, L., 1996. Functional magnetic resonance neuroimaging data acquisition techniques. NeuroImage 4 (3 Pt 3), S76 – S83. Epstein, C.M., Davey, K.R., 2002. Iron-core coils for transcranial magnetic stimulation. J. Clin. Neurophysiol. 19 (4), 376 – 381. Friston, K.J., Holmes, A.P., Worsley, K.J., Poline, J.B., Frith, C.D., Frackowiak, R.S., 1995. Statistical parametric maps in functional imaging: a general linear approach. Hum. Brain Mapp. 2, 189 – 210. Gershon, A.A., Dannon, P.N., Grunhaus, L., 2003. Transcranial magnetic stimulation in the treatment of depression. Am. J. Psychiatry 160 (5), 835 – 845. Herwig, U., Kolbel, K., Wunderlich, A.P., Thielscher, A., von Tiesenhausen, C., Spitzer, M., Schonfeldt-Lecuona, C., 2002. Spatial congruence of neuronavigated transcranial magnetic stimulation and functional neuroimaging. Clin. Neurophysiol. 113 (4), 462 – 468. Kammer, T., Beck, S., Thielscher, A., Laubis-Herrmann, U., Topka, H., 2001. Motor thresholds in humans: a transcranial magnetic stimulation study comparing different pulse waveforms, current directions and stimulator types. Clin. Neurophysiol. 112 (2), 250 – 258. Liu, G., Sobering, G., Olson, A.W., van Gelderen, P., Moonen, C.T., 1993. Fast echo-shifted gradient-recalled MRI: combining a short repetition time with variable T2* weighting. Magn. Reson. Med. 30 (1), 68 – 75. Lorensen, W.E., Cline, H.E., 1987. Marching cubes: a high resolution 3D surface reconstruction algorithm. Comput. Graph. 21 (4), 163 – 169. Lorenzano, C., Gilio, F., Inghilleri, M., Conte, A., Fofi, L., Manfredi, M., Berardelli, A., 2002. Spread of electrical activity at cortical level after repetitive magnetic stimulation in normal subjects. Exp. Brain Res. 147 (2), 186 – 192. Lotze, M., Kaethner, R.J., Erb, M., Cohen, L.G., Grodd, W., Topka, H., 2003. Comparison of representational maps using functional magnetic resonance imaging and transcranial magnetic stimulation. Clin. Neurophysiol. 114 (2), 306 – 312. Luders, E., Rex, D.E., Narr, K.L., Woods, R.P., Jancke, L., Thompson, P.M., Mazziotta, J.C., Toga, A.W., 2003. Relationships between sulcal asymmetries and corpus callosum size: gender and handedness effects. Cereb. Cortex 13 (10), 1084 – 1093. Mazziotta, J.C., Valentino, D., Grafton, S., Bookstein, F., Pelizzari, C., Chen, G., Toga, A.W., 1991. Relating structure to function in vivo with tomographicimaging. CibaFound.Symp.163, 93 – 101 (discussion 101 – 112). Moahker, M., 2002. Means and averaging in the group of rotations. SIAM J. Matrix Anal. Appl. 24 (1), 1 – 16. Paus, T., 1999. Imaging the brain before, during, and after transcranial magnetic stimulation. Neuropsychologia 37 (2), 219 – 224. Paus, T., Sipila, P.K., Strafella, A.P., 2001. Synchronization of neuronal activity in the human primary motor cortex by transcranial magnetic stimulation: an EEG study. J. Neurophysiol. 86 (4), 1983 – 1990. Powell, M., 1964. An efficient method for finding the minimum of a function of several variables without calculating derivatives. Comput. J. 7, 55 – 162. Pridmore, S., Fernandes Filho, J.A., Nahas, Z., Liberatos, C., George, M.S., 1998. Motor threshold in transcranial magnetic stimulation: a comparison of a neurophysiological method and a visualization of movement method. J. ECT 14, 25 – 27. Ramsey, N.F., 1999. Direct comparison of functional MRI and PET. In: Moonen, C.T.W., Bandettini, P.A. (Eds.), Functional MRI. Springer, Berlin, pp. 421 – 432. Ramsey, N.F., Kirkby, B.S., Van Gelderen, P., Berman, K.F., Duyn, J.H., Frank, J.A., Mattay, V.S., Van Horn, J.D., Esposito, G., Moonen, C.T., et

S.F.W. Neggers et al. / NeuroImage 21 (2004) 1805–1817 al., 1996. Functional mapping of human sensorimotor cortex with 3D BOLD fMRI correlates highly with H2(15)O PET rCBF. J. Cereb. Blood Flow Metab. 16 (5), 755 – 764. Ramsey, N.F., van den Brink, J.S., van Muiswinkel, A.M., Folkers, P.J., Moonen, C.T., Jansma, J.M., Kahn, R.S., 1998. Phase navigator correction in 3D fMRI improves detection of brain activation: quantitative assessment with a graded motor activation procedure. NeuroImage 8 (3), 240 – 248. Ray, J.L., McNamara, B., Priest, A., Boniface, S.J., 2002. Measuring TMS stimulus/response characteristics from both hemispheres simultaneously for proximal and distal upper limb muscles. J. Clin. Neurophysiol. 19 (4), 371 – 375. Rutten, G.J., van Rijen, P.C., van Veelen, C.W., Ramsey, N.F., 1999. Language area localization with three-dimensional functional magnetic resonance imaging matches intrasulcal electrostimulation in Broca’s area. Ann. Neurol. 46 (3), 405 – 408. Rutten, G.J., Ramsey, N.F., van Rijen, P.C., Noordmans, H.J., van Veelen,

1817

C.W., 2002. Development of a functional magnetic resonance imaging protocol for intraoperative localization of critical temporoparietal language areas. Ann. Neurol. 51 (3), 350 – 360. van Gelderen, P., Ramsey, N.F., Liu, G., Duyn, J.H., Frank, J.A., Weinberger, D.R., Moonen, C.T., 1995. Three-dimensional functional magnetic resonance imaging of human brain on a clinical 1.5-T scanner. Proc. Natl. Acad. Sci. U. S. A. 92 (15), 6906 – 6910. Woods, R.P., Dapretto, M., Sicotte, N.L., Toga, A.W., Mazziotta, J.C., 1999. Creation and use of a Talairach-compatible atlas for accurate, automated, nonlinear intersubject registration, and analysis of functional imaging data. Hum. Brain Mapp. 8 (2 – 3), 73 – 79. Worsley, K.J., Evans, A.C., Marrett, S., Neelin, P., 1992. A three dimensional statistical analysis for CBF activation studies in the human brain. J. Cereb. Blood Flow Metab. 12, 900 – 918. Worsley, K.J., Marrett, S., Neelin, P., Vandal, A.C., Friston, K.J., Evans, A.C., 1996. A unified statistical approach for determining significant voxels in images of cerebral activation. Hum. Brain Mapp. 4, 58 – 73.