Artifact correction and source analysis of early electroencephalographic responses evoked by transcranial magnetic stimulation over primary motor cortex

www.elsevier.com/locate/ynimg NeuroImage 37 (2007) 56 – 70

Artifact correction and source analysis of early electroencephalographic responses evoked by transcranial magnetic stimulation over primary motor cortex Vladimir Litvak, a,b,⁎ Soile Komssi, c Michael Scherg, d Karsten Hoechstetter, d Joseph Classen, e Menashe Zaaroor, b Hillel Pratt, a,b and Seppo Kahkonen f,g a

Faculty of Biomedical Engineering, Technion-Israel Institute of Technology, Haifa, Israel Evoked Potentials Laboratory, Technion-Israel Institute of Technology, Haifa, Israel c Helsinki Medical Imaging Center, University of Helsinki, Helsinki, Finland d MEGIS Software GmbH, Grafelfing/Munich, Germany e Human Cortical Physiology and Motor Control Laboratory, Department of Neurology, University of Wuerzburg, Germany f BioMag Laboratory, Helsinki University Central Hospital, Finland g Cognitive Brain Research Unit, University of Helsinki, Helsinki, Finland b

Received 22 January 2007; revised 2 May 2007; accepted 7 May 2007 Available online 18 May 2007 Analyzing the brain responses to transcranial magnetic stimulation (TMS) using electroencephalography (EEG) is a promising method for the assessment of functional cortical connectivity and excitability of areas accessible to this stimulation. However, until now it has been difficult to analyze the EEG responses during the several tens of milliseconds immediately following the stimulus due to TMS-induced artifacts. In the present study we show that by combining a specially adapted recording system with software artifact correction it is possible to remove a major part of the artifact and analyze the cortical responses as early as 10 ms after TMS. We used this methodology to examine responses of left and right primary motor cortex (M1) to TMS at different intensities. Based on the artifact-corrected data we propose a model for the cortical activation following M1 stimulation. The model revealed the same basic response sequence for both hemispheres. A large part of the response could be accounted for by two sources: a source close to the stimulation site (peaking ∼ 15 ms after the stimulus) and a midline frontal source ipsilateral to the stimulus (peaking ∼ 25 ms). In addition the model suggests responses in ipsilateral temporo-parietal junction areas (∼ 35 ms) and ipsilateral (∼ 30 ms) and middle (∼ 50 ms) cerebellum. Statistical analysis revealed significant dependence on stimulation intensity for the ipsilateral midline frontal source. The methodology developed in the present study paves the way for the detailed study of early responses to TMS in a wide variety of brain areas. © 2007 Elsevier Inc. All rights reserved. Keywords: TMS; Human; EEG; Artifact correction; Source analysis

⁎ Corresponding author. Sobell Department of Motor Physiology and Movement Disorders, Institute of Neurology, 8–11 Queen Square, London WC1N 3BG, UK. Fax: +44 20 7278 9836. E-mail address: [email protected] (V. Litvak). Available online on ScienceDirect (www.sciencedirect.com). 1053-8119/$ - see front matter © 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2007.05.015

Transcranial magnetic stimulation (TMS) is a technique for noninvasive stimulation of the brain, based on the principle of electromagnetic induction (Barker et al., 1985; Hallett, 2000). In this method, short current pulses are driven through an induction coil, which is placed above the subject's head. The consequential transient magnetic field induces an electrical current in the tissues of the head. This current is strongest at the cortical surface. With the commonly used stimulation intensities and coils, the cortex is activated to a depth of 2–3 cm, and the activated surface area is several square centimeters (Ruohonen and Ilmoniemi, 1999; Thielscher and Kammer, 2002). The most prominent effect of a single TMS pulse is its ability to induce activity in contralateral muscles when applied over the primary motor cortex (M1). This has been the basis for the use of TMS in the diagnosis of motor nerve conduction velocity. When applied over the occipital visual areas TMS can induce visual sensations called phosphenes. However, these sensations are unreliable and have been used for research, but not for clinical diagnostics (Antal et al., 2004; Amedi et al., 2005). By applying series of stimuli, it is possible to induce longer-term effects. At the motor cortices, these can be detected as suppression or facilitation of corticospinal excitability (Pascual-Leone et al., 1994; Chen et al., 1997; Wolters et al., 2003; Huang et al., 2005). Repetitive TMS (rTMS) of association cortices can disturb performance of cognitive tasks, as the induced activation interferes with ongoing neuronal processes. rTMS of the dorsolateral prefrontal cortex has antidepressant effect and has, therefore, been used in the therapy of depression (Pascual-Leone et al., 1996; Klein et al., 1999). Combining TMS with neuroimaging methods such as positron emission tomography (PET), functional magnetic resonance imaging (fMRI) and electroencephalography (EEG) expands the

V. Litvak et al. / NeuroImage 37 (2007) 56–70

applicability of TMS to the study of cortical reactivity and connectivity. Studies with PET (Fox et al., 1997; Paus et al., 2001a; Fox et al., 2005) and fMRI (Bohning et al., 1999, 2000; Baudewig et al., 2001; Nahas et al., 2001; Bestmann et al., 2004, 2005; Denslow et al., 2005) revealed TMS-induced changes in blood flow and oxygenation both at the stimulation site and at remote sites. However, the low time resolution of these methods did not allow establishing the time course of activation of these areas. Even though multi-channel EEG could have provided some of the lacking information, TMS-evoked potentials recorded with standard EEG amplifiers are too distorted by TMS-induced artifacts. These slowly decaying transients or saturations in the recording do not allow examining the time period immediately following the stimulus which is interesting for studying functional connectivity. Extensive research to develop a TMS-compatible EEG system resulted in the construction of an amplifier with sample-and-hold circuits capable of suppressing a large part of the TMS-induced artifact (Virtanen et al., 1999) and allowing EEG recording as early as 7 ms after the TMS pulse (or even 2 ms in the first study (Ilmoniemi et al., 1997)). Studies using this system were able to follow the propagation of the TMS induced activation to other ipsilateral and contralateral cortical areas (Ilmoniemi et al., 1997, 1999; Komssi et al., 2002; Kahkonen et al., 2004, 2005a) and to examine the effect of stimulation intensity (Komssi et al., 2004; Kahkonen et al., 2005a) and movement (Nikulin et al., 2003) on TMS-evoked potentials. Recently, using a commercial version of the same EEG system it was possible to show differences in functional connectivity between wakefulness and deep stages of sleep (Massimini et al., 2005a), to study the effects of TMS-induced long-term potentiation (LTP) in the motor cortex (Esser et al., 2006) and to characterize abnormalities in EEG responses to TMS in schizophrenic patients (Massimini et al., 2005b). These studies have shown the potential of combining TMS and EEG for the study of cortical function in humans. The applicability of the TMS-EEG method is presently limited since the artifact cannot be completely suppressed using specific hardware. In the data of Komssi et al. (2004) and Kahkonen et al. (2005a,b), residual artifacts were present especially for high stimulation intensities. The authors dealt with this problem by excluding channels strongly affected by the artifact. This is problematic, however, since these channels are usually the ones closest to the stimulation site and, therefore, the most informative about the early phase of the response. In addition, exclusion of channels from a particular part of the scalp is likely to introduce biases in source estimation (Michel et al., 2004). An alternative approach to the use of an EEG amplifier with a sample-and-hold circuit is offline correction of the stimulus artifact. This requires that the EEG responses are recorded with an amplifier with wide operating range that does not become saturated by the TMS pulse. Thut et al. described a method where the artifact was removed by subtracting an event related potential (ERP) induced by TMS alone from an ERP induced by TMS during a task (Thut et al., 2003). This approach is suitable for studying the effect of TMS on potentials related to processing of sensory stimuli and cognitive tasks. However, it does not allow studying the responses evoked by TMS itself since these responses are eliminated by the subtraction procedure (Thut et al., 2005). In some cases in which only late components of the ERP are of interest, the parts of the recording affected by the artifact can be excised from the recording in order to prevent their interference with digital filtering and the remaining

57

part of the ERP can be analyzed by standard methods. Obviously, this method is limited by its inability to examine the early part of the response (Fuggetta et al., 2006). In the present study, we used a methodology originally developed for correcting eye-movement artifacts based on their scalp topography (Berg and Scherg, 1994; Ille et al., 2002) in order to correct the residual TMS artifact. This method made it possible to substantially reduce the contamination of the early TMS-evoked ERPs by the stimulus artifact without distorting the physiological components of the response. Our goal was to apply this multiplesource artifact-correction method to previously published data (Komssi et al., 2004) in order to analyze the sources of the early EEG responses to magnetic stimulation of the left and the right motor cortex. Methods Experimental procedures A total of 7 healthy volunteers (2 female, 2 left handed, age 28 ± 3) participated in the study after giving a written informed consent. The study protocol was approved by the Ethics Committee of Helsinki Central University Hospital. Biphasic magnetic pulses (duration 385 μs) were delivered with a custom-made figure-of-eight coil consisting of two coplanar circular loops (40 mm Ø) so that the induced current was in the posterior–anterior direction. The electromyogram was recorded bilaterally from the abductor digiti minimi muscles (ADM). The optimal position for stimulation of the ADM was found separately in each hemisphere. The resting motor threshold (RMT) was determined separately for the left and right M1 by adjusting the TMS intensity until 50 μV motor evoked potentials were evoked in about 50% of the trials. During stimulation, the subject sat in a reclining chair with the neck and back supported by pillows, arms relaxed and eyes closed, wearing ear plugs. TMS was delivered to the left and right M1 in sequences of 50 pulses with an interstimulus interval of 1.5–2.5 s and with constant intensity. Four sequences, with intensities of 60, 80, 100, and 120% of RMT, were targeted in random order to the two stimulation sites. This was repeated after a 5–10 min break. The EEG was recorded from 60 scalp electrodes (Ag/AgCl C-shaped) referred to the forehead. The saturation of the EEG amplifiers during the TMS pulse was avoided by using a sampleand-hold circuit (Virtanen et al., 1999) that kept the outputs locked from 50 μs pre- to 7 ms post-stimulus to avoid saturation. The data from 100 ms pre- to 300 ms post-stimulus were bandpass-filtered between 0.1 and 500 Hz and sampled at the rate of 1450 Hz. Detailed description of the experimental procedures can be found in the previous publication (Komssi et al., 2004). Analysis of the EEG data The analysis presented here improved the previously presented analysis methods (Komssi et al., 2004; Kahkonen et al., 2005b) in three main respects: 1. We applied the artifact correction method previously suggested for common EEG artifacts such as eye blinks (Berg and Scherg, 1994) to the residual TMS artifact and were able to show that, in combination with the hardware artifact-blocking circuit (Virtanen et al., 1999) and averaging across subjects, this methodology

58

V. Litvak et al. / NeuroImage 37 (2007) 56–70

makes it possible to remove the artifact to a very large extent at all the stimulation intensities tested. 2. Based on the corrected data, we were able to derive a source model for the TMS-evoked potentials. This model allowed us to separate different components of the response and estimate propagation delays of TMS-evoked activation for different brain areas. 3. We applied a novel framework for the statistical analysis of neuroimaging data (Osipova et al., 2006; Takashima et al., 2006; Maris and Oostenveld, in press) in order to assess the significance of the effect of TMS intensity on the amplitudes of the different source activities. Thus, our analysis procedure consisted of 3 stages: correction of the TMS-induced artifact, computation of the cortical activity from corrected EEG data using source modeling and statistical analysis of the brain source activity estimated for each subject and condition. In the present study we report the analysis of the early part of the TMS-evoked responses between 10 and 60 ms after the stimulus. The later part of the response is in most cases not contaminated with the residual artifact and can be analyzed by standard methods (Bender et al., 2005; Fuggetta et al., 2006). Correction of the artifact The method used for artifact correction in the present study was developed by Berg and Scherg (Berg and Scherg, 1994). It is implemented in the BESA software (Version 5.1, MEGIS Software GmbH, Gräfelfing, Germany) that was used for the pre-processing and source analysis of the data. The artifact correction we applied is based on a multiple source approach: A source model is constructed that consists of a set of artifact topographies (obtained in the present study from a principal component analysis (PCA) decomposition of the averaged artifact) and a set of brain topographies (multiple dipoles that model brain activity). From this source model, a linear inverse operator is computed that decomposes the data into a linear combination of brain and artifact activities. The estimated artifact activities are then subtracted from the data. The modeled brain activity is not altered. This multiple source approach thus aims at maintaining the brain activity undistorted after the correction process. The major advantage over PCA- and independent component analysis (ICA) based artifact correction methods is that no assumption is made regarding spatial or temporal orthogonality or independence of brain and artifact activities. Given appropriate models for brain and artifact topographies, the discrete multiple source inverse allows for separation of brain and artifact activities despite spatial or temporal correlation. A mathematical description is given in the following. Under the assumption of linear summation the relationship between the scalp voltage and the activity of the underlying sources can be written as follows: V n
t ¼ Ln
m d S m
t

ð1Þ

Vn×t is the matrix of observed voltages (n—number of channels, t—number of time frames). Each row of Vn×t is a waveform recorded by one of the electrodes and each column is a voltage map on the scalp at a specific time point. Sm×t is the matrix of source activities, where each row is a source waveform. m is the number of sources.

Ln×m is called a ‘lead field matrix’ and each of its columns is a ‘scalp topography’ or ‘scalp signature’ of the corresponding source. This topography would be observed on the scalp if this source was the only active source contributing to the EEG. A specific case which is of interest in the present study is the use of discrete multiple source modeling for correcting artifacts. There are artifacts present in practically any EEG recording, such as eye movements, blinks and muscle artifacts and this study deals specifically with artifacts generated by TMS. Under the assumption of linear summation, artifacts can be presented as part of the sources in matrix S. 0 a1 1 N Sta1 S1 0 1 B v dd v C b d B ar ar C l1a1 N l1ar l1b1 N l1p B N S S t C B C 1 C v ddd v Ad B Vn
t ¼ @ v d d d v b b 1 B S1 N St 1 C bp B C a1 ar b1 ln N ln ln N ln @ v ddd v A b b S1 p N St p a a b b ð2Þ ¼ Ln
r d S r
t þ Ln
p d S p
t In this case, the topography of the artifacts is approximated by r spatial components and the brain activity is represented by a model with p sources. The equation can be solved by inverting the combined leadfield matrix L with the appropriate mathematical method (in BESA Moore–Penrose pseudoinverse is computed with Tikhonov regularization). The main difference between modeling artifacts and brain sources is that artifacts can be measured using their averaged topography on the head while brain sources require a forward model based on equivalent current dipoles (ECD) and a particular volume conductor model of the head to estimate the surface topography. The leadfield matrix is composed of both the artifact and dipole source topographies. Inverting this matrix intrinsically leads to the separation of the related artifact and brain source activities of the modeled regions (Berg and Scherg, 1994). The artifact topographies can be estimated by performing PCA of an artifact dominated time segment (Ille et al., 2002). PCA is a linear transformation that chooses an orthogonal coordinate system for the data topographies and waveforms such that the greatest variance is accounted for by the first component, the next largest variance by the second etc. Since the dimensionality of the blink and TMS artifact topographies is small, only one or a few PCA components are required to represent these artifacts, thus leaving other dimensions in the sensor space to represent brain activity, e.g. by model sources. Note that the rows of the inverse matrix related to the artifact topographies are dependent on the brain source topographies and vice versa. Therefore, the quality of separation between brain activity and the artifacts depends on both the model of the artifact and the model of the brain activity. This poses a problem, since an adequate model of the brain activity cannot be built in the presence of the artifact, whereas the artifact cannot be corrected well without such a model. In studies where the generators of evoked potentials in question are well known, standard models, based on previous studies, can approximate the brain activity rather well (Scherg et al., 2002). Since the present study is, to our knowledge, the first attempt to suggest a detailed model for the potentials evoked by TMS, we had to base our artifact correction on an iterative approach consisting of two stages. In the first stage the brain activity was represented by a surrogate model (Berg and Scherg, 1994) that covered all major cortical regions by 15 regional sources (Scherg

V. Litvak et al. / NeuroImage 37 (2007) 56–70

59

and Von Cramon, 1986; Scherg and Berg, 1996). This model is capable of reasonably approximating any activity generated in the cerebral cortex (Scherg et al., 2002). After artifact correction based on this model, we used the grand average of the corrected data to construct a source model specifically accounting for the TMSevoked activity. In the second stage, the correction of the original data was redone using the grand average source model in conjunction with the previously defined individual artifact topographies. Thus, artifact free averages could be obtained for each subject without distortion of the underlying model source topographies (Berg and Scherg, 1994). In summary, the analysis steps for Stage I of the analysis were as follows (flowchart of the analysis is shown in Fig. 1): 1. The EEG channels that were obviously disconnected during the recording or severely contaminated by noise were identified by eye and excluded from further analysis (see also (5) below). 2. For most subjects, the EEG epochs did not contain slow deflections associated with eye blinks. In cases where eye blinks were present, they were identified using the automatic search and average procedure of the BESA program. A topography of the eye blink artifact was then estimated based on these averages and was used in combination with TMS artifact topographies (see below) for artifact correction. 3. TMS-evoked ERPs were computed by averaging baseline corrected (baseline − 100 to 0 ms) data time-locked to the TMS triggers. The segments following the stimulus were examined for the presence of TMS artifacts, defined as sharp peaks that clearly exceeded normal ERP amplitudes (i.e. N 30 μV). The topographies of the artifact were estimated by performing PCA of a segment around these peaks, usually the first 15 ms after the TMS trigger. The data was then corrected using the defined artifact topographies with the method described above. At this stage, the surrogate model with 15 regional sources (Scherg et al., 2002) was used for artifact correction. If necessary, additional artifact topographies were added in a sequential procedure until there were no more obvious artifacts in the average. For experiments in which more than one location was stimulated and the coil was moved between the two repetitions of stimulating the same location, we found that the artifact topography differed between the two repetitions. We therefore applied our correction procedure separately for each repetition. 4. The corrected epoched data were then examined for other artifacts using the BESA artifact rejection interface (Scherg et al., 2004). Trials were considered artifactual and excluded from averaging if they satisfied one of the following two conditions: 1) The amplitude difference between the largest and the smallest sample within the trial in any channel exceeded the maximal amplitude threshold. 2) The voltage difference between two adjacent samples in any channel exceeded the maximal gradient threshold. The thresholds were adjusted individually using a graphical user interface to balance between the need to retain trials for analysis and the need to reject artifact-contaminated trials. The maximal amplitude threshold was 208.7 ± 47.5 μV. The maximal gradient threshold was 83.6 ± 28.9 μV. In some cases, additional channels were excluded from analysis at this stage in order to retain more trials (see also (5) below).

Fig. 1. Flowchart for artifact correction and analysis of TMS-evoked ERPs. See Methods for details.

5. Averaged responses were computed for each data file twice — with and without artifact correction. The two averages were compared in order to verify that only the artifact was eliminated and the activity at later times and at other scalp areas was not affected. In some cases, channels were excluded at this stage, if the activity after correction was clearly distorted and not con-

60

V. Litvak et al. / NeuroImage 37 (2007) 56–70

1988) of the model sources and the related brain areas are given in Table 1. For the grand average data used for constructing the model the goodness-of-fit of the model was 89.7% for the response to left M1 stimulation and 89.2% for the right M1 stimulation. Values of the goodness-of-fit of the model for individual subjects and all stimulation intensities are provided in supplementary Table S1. The final grand-average model with the 9 fixed dipoles was used to construct a fixed source montage, i.e. an inverse spatial filter (Scherg et al., 2002; Hoechstetter et al., 2004), that allowed to decompose the data of each subject and each condition into source waveforms for statistical analysis.

sistent with activity at the neighboring channels. The total number of channels excluded from analysis at stages 1, 4 and 5 (averaged across all analyzed files) was 4 ± 3, maximum 13. 6. Excluded channels were interpolated using spherical spline interpolation (Scherg et al., 2002), and for each subject, intensity and location an average was computed over the two repetitions of stimulating with the same parameters. 7. From these data grand averages across subjects were computed for each location and each intensity. For the creation of a grand average source model, one common grand average was computed over the intensities of 80%, 100% and 120%, separately for stimulation of the right and left M1.

Computing the individual activities of the model sources In Stage II of the analysis (see Fig. 1) the uncorrected individual averages were corrected using the grand average model combined with the individual artifact topographies to generate the individual source waveforms. This was done for each averaged data file separately. The corrected source waveforms were filtered with a forward 5-Hz high-pass filter. Averages were subsequently computed over the two files representing the same condition and subject. Grand averages across subjects were also computed for each condition.

Modeling the TMS-evoked responses The grand averages across all subjects and intensities from 80% to 120% were analyzed using the sequential multiple source analysis strategy (Scherg and Berg, 1996). In order to remove lowfrequency activities, the data were filtered with a forward 5-Hz high-pass filter. First, a regional source was fitted to the earliest epoch that was free of residual artifact (11–20 ms). The regional source was oriented, and the predominant dipole retained as the first source of the multiple source model. Then, symmetric pairs of regional sources were added to the model in a sequential procedure to account for the residual activities in the next baseline-to-peak interval. Again, after each pair of regional sources was added and fitted, the sources were rotated at their earliest peak and converted into single dipoles retaining the predominant orientation. The right and left TMS data were fitted in alternation. This showed a fair symmetry of the sources with some of the ipsilateral activities being matched by contralateral activities. Since the contralateral activities were weak and showed a similar orientation as compared to the ipsilateral activity when stimulating the other side, the symmetry constraint was released and the ipsilateral dipoles were refitted in their respective intervals of the ipsilateral stimulation conditions where they had a larger signal. Thus, a grand average model could be created that contained comparable sources ipsiand contralateral to each stimulation while the fitting optimization was only performed on the ipsilateral, i.e. stronger source activities. Thereby, the uncertainty of fitting a small contralateral component with unstable location was overcome. The grand average model comprised three nearly symmetric pairs of dipoles in anterior motor (AM), midline-frontal-motor (MFM) and temporoparietal junction areas (TPJ) as well as three dipoles in the cerebellum (Cb). The Talairach coordinates (Talairach and Tournoux,

Statistical analysis The source waveforms were analyzed statistically using the open source Fieldtrip toolbox for Matlab (http://www.ru.nl/fcdonders/ fieldtrip/, F.C. Donders Centre for Cognitive Neuroimaging, Nijmegen, The Netherlands). The nonparametric statistical test implemented in this toolbox is described in detail in Maris and Oostenveld (in press) and is similar to the well established nonparametric approach used for the analysis of functional imaging data (Nichols and Holmes, 2002). A similar approach has previously been used in several published studies (Osipova et al., 2006; Takashima et al., 2006). The null hypothesis of this nonparametric statistical test is that the EEG data in the different experimental conditions come from the same probability distribution. To maximize the sensitivity of the statistical test, a test statistic is used that takes advantage of the fact that the effects produced by physiological sources usually span several consecutive time frames. Conversely, differences between experimental conditions that are due to noise have a much smaller temporal correlation. Thus it is possible to increase the sensitivity for discovering physiological effects by using a test statistic that is based on clusters of neighboring time frames that exhibit the same effect.

Table 1 The model fitted to TMS evoked ERPs (Talairach coordinates) Source name

Side

Anterior motor (AM)

Left Right Left Right Left Middle Right Left Right

Midline frontal motor (MFM) Cerebellum (Cb)

Temporo-parietal junction (TPJ)

Location

Orientation

Brain areas within 1 cm radius

X

Y

Z

X

Y

Z

−24.5 19.7 −12.4 8.5 −36.9 −3.7 30.1 −42.2 31.0

7.8 20.6 14.5 23.8 − 57.7 − 69.2 − 63.6 − 22.1 − 29.1

60.3 54.8 40.2 34.1 −27.3 −21.7 −28.0 28.8 38.0

0.4 − 0.4 − 0.5 0.1 − 0.4 0.1 0.4 − 0.2 0.2

0.9 0.9 0.1 0.2 0.5 0.0 0.0 0.5 0.4

−0.2 −0.2 0.9 1.0 −0.8 −1.0 −0.9 0.8 0.9

BA 6 BA 6, 8 BA 32 BA 32, 6, 8, 9, 24 Left cerebellum Left cerebellum, right cerebellum Right cerebellum BA 2 BA 2, 3, 40

V. Litvak et al. / NeuroImage 37 (2007) 56–70

The nonparametric cluster-based statistical test is performed in two steps. In the first step, the effect of interest is evaluated for all time frames in the analysis window. This is performed by means of a time-specific statistic. Time frames whose time-specific statistic exceeds some prior threshold are selected for clustering in the second step. As will become clear later on, the type of statistic and the value of the prior threshold does not affect the Type-I error rate of the statistical test. In the second step, the consecutive time-frames where the timespecific statistic crossed the threshold are grouped into clusters. After clusters are identified, each cluster is assigned a cluster-level statistic, whose value equals the sum of the point-specific statistics for all the point in the cluster. Thus, the cluster-level statistic depends on both the extent of the cluster and the magnitudes of the point-specific statistics at the points that belong to this cluster (Takashima et al., 2006). The test statistic for which a statistical reference distribution is constructed, is the maximum over all cluster-level statistics. The statistical reference distribution is the so-called permutation distribution. This distribution is obtained by randomly permuting the condition labels over the condition-specific data. For each of these permutations, the maximum cluster-level statistic is recorded and these are used to construct a histogram that approximates the permutation distribution. In the present study, this histogram was created from 1000 random draws. To evaluate the significance of the observed maximum cluster-level statistic (i.e. the maximum cluster-level statistic that was obtained with the original labeling), this statistic was compared against the approximate permutation distribution. This comparison results in a p-value: the proportion of permutations in which the maximum cluster-level test statistic exceeds the observed maximum cluster-level test statistic. This proportion is called a Monte Carlo p-value in the statistics literature. A thousand random draws provide an accurate estimate of the true p-value. These are the p-values reported for the nonparametric tests in the present study. The time-specific statistic used in the present study was derived from regression analysis for finding temporal patterns that are significantly dependent on TMS intensity. Since our data consists of repeated measurements for the same subjects, classic linear regression would not be appropriate. Our analysis is thus based on dependent-samples regression statistic suggested by Lorch and Myers (1990). In the first stage, the idea of this statistic is to compute the regression coefficient β between the independent variable X and dependent variable Y (in our case, the TMS intensity and the neural response, respectively) separately for each subject: K X

b ̂s ¼

ðYi  Y¯ ÞðXi  X¯ Þ

i¼1 K X

ð3Þ ðXi  X¯ Þ2

i¼1

where K is the number of samples for each subject (in our case K = 4 because four TMS intensities were tested in all experiments). In the second stage, the dependent-regression t statistic is computed from the individual regression coefficients by dividing the mean regression coefficient by its standard error across subjects: Tdsr ¼

b¯ SEðbÞ

ð4Þ

61

In order to find the aspects of the response significantly dependent on TMS intensity, this statistic can be tested against a t-distribution with N − 1 degrees of freedom (N is the number of subjects). Results Correction of the residual artifact The data was recorded with an artifact blocking sample-and-hold circuit. Fig. 2A illustrates the dependence of the magnitude of the residual artifacts not blocked by this circuit on stimulus intensity. There was a general trend for increase in artifact magnitude with TMS intensity. However, at the same time the artifact amplitude varied greatly between subjects and in one case high artifact amplitudes were observed even for relatively low TMS intensities. For the specific case of average between responses evoked by 80%, 100% and 120% RMT stimulation, used in subsequent analysis, artifact amplitudes are shown for all individual cases (Fig. 2B). In total, 112 raw data files were corrected separately using the procedure described in Methods. In 2 cases, no correction was necessary, since no clear artifacts were found. In 94 cases, single artifact topography was sufficient to correct the artifact. In 11 cases, two topographies were used, in 3 cases, three topographies, and five and six topographies were necessary in a single case each. These numbers include both TMS artifact and eye-blink artifact topographies. In order to evaluate the usefulness of artifact correction we calculated the percentage of trials and channels excluded from analysis due to artifacts for uncorrected and corrected data using the same rejection criteria (see Methods). The results are presented in supplementary Fig. S1. An example of the effect of software correction is shown for a case with intermediate artifact amplitude (Fig. 2C). In the uncorrected data (Fig. 2C, Uncorrected), sharp peaks can be seen around the stimulation site (right M1), exceeding by far the magnitude typical of evoked potentials, as can be seen by their comparison with later peaks. The response returns to baseline only about 30 ms after the stimulus. The residual artifact, therefore, completely masks the early ERP response. After correction, the sharp peaks are completely removed (Fig. 2C, Corrected). Prominent scalp distribution patterns in the corrected data To reveal the most robust features of the response and to reduce noise, the responses to TMS at 80%, 100%, and 120% RMT were corrected using the surrogate model and averaged. The grand average of these individual averages was examined and subjected to multiple source analysis. Fig. 3 shows superimposed multi-channel plots (“butterfly plots”) of the grand-average responses to stimulation of left and right M1 after Stage I correction. Scalp potential distribution at the prominent response peaks is shown. Note the remarkable symmetry between the scalp maps at closely corresponding times between the left and right M1 stimulation. The earliest pattern observed after the period immediately following the amplifier reconnection (and still contaminated by artifacts) is negativity ipsilateral to the stimulation around 10–12 ms after stimulus. This pattern is followed by a dipolar pattern centered over the side of stimulation with a midfrontal positivity and an ipsilateral parietal negativity (15–25 ms). The next pattern is widespread anterior-central positivity lateralized to the side of stimulation (25–

62

V. Litvak et al. / NeuroImage 37 (2007) 56–70

Fig. 2. Distribution of artifact peak amplitude and the effect of correction. (A) TMS intensity-dependence of the distribution of residual artifact amplitude. The residual artifact amplitude was defined as the maximal absolute difference between the software corrected and uncorrected TMS-evoked ERPs across all latencies between − 20 and 60 ms poststimulus and across all channels. Data was pooled across stimulation sites. Box and whisker plot (McGill et al., 1978) was used to show the distribution. The boxes have lines at the lower quartile, median, and upper quartile values. The whiskers are lines extending from each end of the box to the extreme data values not defined as outliers (see below). Points lying beyond 1.5 times the interquartile range from the median are shown as outliers (empty circles). In this case all outliers belong to a single subject (Subject 1). (B) Residual artifact amplitude (defined as in A) for all subjects, data averaged between 80%, 100% and 120% RMT and stimulation of right M1 (black bars) and left M1 (grey bars). (C) Example of the effect of computational correction for a single subject with intermediate residual artifact amplitude (subject 5, TMS applied to right M1). The data was averaged between 80%, 100% and 120% RMT stimulation intensities. Here, as in most other cases, the maximal artifact amplitude was attained at channels close to the stimulation site.

35 ms). This is followed at 35–50 ms by a dipolar pattern similar to that of 15–25 ms with reversed polarity. In addition, a dipolar pattern appears around 30 ms over the occipital region with an occipito-parietal negativity and an ipsilateral inferior positivity, suggesting cerebellar origin (Fig. 3).

AMi (first peaking at 15 for left and at 12 ms for right M1 stimulation), MFMi (peaking at 28 ms—left, 26 ms—right), MFMc (peaking at 28 ms—left, 30 ms—right), Cbi (31 ms—left, 30 ms— right), TPJi (36 ms—left, 34 ms—right) and Cbm (48 ms left, 50 ms right).

Source model reveals similar response sequence for left and right M1 stimulation

Dependence of the source waveforms on stimulation intensity

Fig. 4 shows the source waveforms and scalp projections of equivalent current dipoles of the model fitted to the grand average of the responses to 80%, 100% and 120% RMT, combined. The names, locations and orientations of model sources are shown in Table 1. To facilitate the comparison between left and right M1 stimulation, in the following discussion, we use suffixes i and c to denote brain structures ipsilateral and contralateral to stimulation, and m to denote structures close to the midline. Clear symmetry can be observed between the responses to the left and right M1 stimulation. In addition, the model reveals a likely sequence of signal propagation following magnetic stimulation of M1. The sequence is

To analyze the intensity dependence of the response of each of the model sources on TMS intensity, non-parametric statistical analysis of the source waveforms was performed using dependentsamples regression statistic described in Methods. The results of the analysis are shown in Fig. 5. The sources included in the analysis were the six sources shown in Fig. 4. Only the MFMi sources showed significant intensity dependence at similar latencies for both left and right M1 stimulation (p b 0.05 for stimulation of left M1 (26–32 ms) and p b 0.01 for stimulation of right M1 (27– 35 ms)). For right M1 stimulation intensity-dependent temporal clusters were also observed for AMi (28–54 ms p b 0.01), Cbi (48– 54 ms p b 0.05), and Cbm (35–45 ms p b 0.01). For left M1

V. Litvak et al. / NeuroImage 37 (2007) 56–70

63

Fig. 3. Prominent scalp potential patterns in the artifact corrected data. The data were averaged across 80%, 100% and 120% RMT intensities and grand averaged across subjects. The arrow marks the end of the time window where the amplifiers were disconnected. Note the symmetry between the responses to left and right M1 stimulation.

stimulation a significant cluster was found for Cbm (11–14 ms, p b 0.05) but the short latencies indicate that this cluster was probably a case of a false positive. Discussion In the present study we showed that a previously published artifact correction method (Berg and Scherg, 1994) makes it possible to almost completely remove the TMS-induced artifact and reveal physiological EEG responses to TMS, starting as early as 10 ms after the magnetic stimulus. We applied the artifact correction to responses evoked by magnetic stimulation of the left and right M1. Based on the corrected data we suggest a model for the generation of these responses and their dependence on stimulation intensity. The physiological significance of the findings Temporal sequence of the response The exploration of the brain's functional connectivity has always been an important motivation behind the attempts to combine magnetic stimulation with EEG recording (Komssi and Kahkonen, 2006). The idea is to stimulate a particular cortical area and to

determine other areas to which the activity propagated at the early stages of the response. The first study that utilized the TMScompatible EEG amplifier developed by Virtanen et al. (1999) showed, using a minimum norm inverse solution, that after stimulation of M1, the activity propagates in about 22 ms to the homologous contralateral area (Ilmoniemi et al., 1997). Since then, further attempts to present the spatiotemporal sequence of TMSevoked activation with distributed currents have been scarce (Komssi et al., 2002; Massimini et al., 2005a). In the present study, we use current dipoles to model the source activity evoked by TMS of the left and right M1. Our source distribution partially disagrees with and partially extends the previous findings. As could be expected, the earliest response was observed close to the stimulation site. The early negativity seen at 11 ms for left M1 and at 10 ms for right M1 stimulation (Fig. 3) is probably the immediate response of the tissue directly stimulated by TMS. This response was brief and largely obstructed by the disconnection of the recording system and the residual artifact and therefore was not modeled in the present study. In the next stage, activation could be well represented by a tangentially oriented AMi source peaking around 15 ms for stimulation of left M1 and 12 ms for stimulation of right M1. The coordinates of both left and right AM sources are

64

V. Litvak et al. / NeuroImage 37 (2007) 56–70

Fig. 4. The multiple source model fitted to the grand averaged data and the corresponding source waveforms. The side panels show the model dipoles from different view angles (AM — red, MFM — blue, TPJ — green, Cb — magenta (left, right), cyan (mid-line)). In the middle — source waveforms separately for left and right M1 stimulation and scalp projections for the sources for which clear peaks were observed. The black bars mark peaks mentioned in the text. The magnitude of the scalp potentials shown in the scalp maps corresponds to the marked times. Note the similarity of the scalp projections to raw data scalp patterns of Fig. 3.

consistent with a location close to the middle frontal gyrus, at Brodmann area 6. The tangential orientation of the ECD indicates that the cortical generator is located in the wall of a sulcus. Based on various studies of sensory cortical responses using intracortical current-source-density estimates, Mitzdorf (1985) suggested that the initial response of a sensory area to thalamic input recorded above the cortical surface consists of a positivity produced by current flow from the cell bodies to apical dendrites of pyramidal cells in layer 4, followed by a negativity produced by current flow in the opposite direction as a result of synaptic input at the apical dendrites. If this basic sequence also applies to the observed tangential dipole pattern of the AMi source, then its generator should be located at the posterior wall of a sulcus. Thus, location and orientation of this source point to the posterior wall of the precentral sulcus. The following response from the midfontal region (MFMi) exhibited a very prominent positive peak around 28 ms in the scalp

maps. The ECD fitted to this response seems to be relatively deep below the cortical surface and its location and orientation is consistent with a source in the superior wall of the ipsilateral cingulate gyrus or, possibly, supplementary motor area (SMA). However, these interpretations should be taken with caution, because it is possible that these regions are activated via corticocortical connections with a different sink-source distribution over the cortical layers. Therefore, the polarity of the dipole fields should not be over-interpreted. The sequence of scalp patterns explained in our model by the AM and MFM sources has been described in the original publication of Komssi et al. (2004) and recently independently confirmed by Bonato et al. (2006) using an fMRI-compatible EEG system where TMS artifact is relatively brief. The responses in the cerebellum modeled by the sources Cbi and Cbm have not been previously described in EEG studies. On one hand, existence of these responses should not be surprising since the

V. Litvak et al. / NeuroImage 37 (2007) 56–70

65

Fig. 5. EEG responses' dependence on stimulus intensity of M1 stimulation. Grand average source waveforms for the four intensities are shown (120% — solid black, 100% — dotted black, 80% — solid grey, 60% — dotted grey). The black bars indicate time periods where significant intensity dependence was found for a particular source. The p-values are not corrected for multiple comparisons between sources.

motor cortex projects to neo-cerebellum via the pontine nuclei (Ramnani, 2006). On the other hand, the responses observed in the present study were mainly accounted for by the ipsilateral and the median cerebellar sources, whereas M1 projects to the contralateral cerebellum. Furthermore, in previous studies, the contralateral cerebellum was also found to participate in the motor planning network based on the imaging of cortical sources showing coherence with muscle activity (Schnitzler et al., 2006). However, the observed occipital-inferior dipole field around 30 ms was very prominent in the 3D topographic maps over the ipsilateral side and clearly separated in timing from the other responses. Thus, the

model clearly required sources in the inferior occipital region. Further studies are necessary to confirm the existence of these cerebellar sources and to characterize them in more detail. The responses in the TPJ sources might be either a direct consequence of the magnetic stimulation coming perhaps from the somatosensory system or physiological artifacts such as middle latency auditory responses to the coil click or somatosensory responses to scalp stimulation. Location and orientation of the two relatively weak TPJ sources suggest an origin in the auditory cortex. Also, the peak latencies of 33 ms are compatible with the middle latency auditory P30 components of responses to the coil click. A

66

V. Litvak et al. / NeuroImage 37 (2007) 56–70

physiological somatosensory evoked potential due to scalp stimulation may have also contributed to these components, but their latency appears too long for such activity. The model revealed remarkable symmetry between the responses to stimulation at either side. Although there are both structural and functional differences between the two hemispheres (Haaland and Harrington, 1996), the basic sequence of events following M1 stimulation seems to be similar in both cases. The latencies of the responses in the right hemisphere seem to be a few milliseconds shorter than the latencies of the analogous responses in the left hemisphere. This difference could possibly be due to some remaining artifact at high intensities, e.g. for the MFMi source that had a larger difference and stronger residual artifact as compared to the other source activities. The response latencies revealed by our model can be compared with propagation delays obtained in paired-pulse stimulation studies in monkeys (Cerri et al., 2003; Shimazu et al., 2004) and in humans (Civardi et al., 2001; Daskalakis et al., 2004; Baumer et al., 2006; Koch et al., 2006). In the study of Civardi et al. (2001) conditioning TMS stimuli were applied 3–5 cm anterior to the hand motor area (this location might correspond to our AMi source) or 6 cm anterior to the vertex on the nasion–inion line (this location might correspond to our MFMi and MFMc sources). In both cases the effect on M1 was inhibitory and the optimal delay between the conditioning stimulus and M1 stimulus was about 6 ms. Slightly longer delays were found between left and right M1 (Baumer et al., 2006). In the study of Daskalakis et al. (2004) conditioning at the cerebellum affected M1 with a delay of about 5–7 ms. These values are consistent with the fact that corticocortical connections conduct at a velocity of about 10–20 m/s (Kandel et al., 2000). Therefore, theoretically, excitation can propagate from M1 to any other cortical area in the first 10 ms and we would expect that the earliest EEG responses to TMS are completely missed by our recording method due to the blocking of recording immediately after the stimulus. There is no way to refute this suggestion in the framework of the present study. However, in practice, it is not likely to be the case since delays between ERP peaks are comprised not only of the time of propagation in axons, but also of the processing delays at the cellular and network levels. Thus, it may take some time from the arrival of the first action potentials from another brain area until a potential gradient generated by the pyramidal cells is detectable on the scalp. Therefore, it would be surprising to find delays between activations of ECD model sources which are shorter than expected conduction times between the corresponding brain regions, but it would not be surprising to find delays which are much longer. In this sense our results are consistent with paired pulse stimulation data. Some constraints can also be drawn on the possible propagation pathways from the known conduction delays. For instance, since in the case of left M1 stimulation MFMi and Cbi were activated with a very short delay of about 4 ms, it is unlikely that the input to cerebellum comes from the SMA or CMA and more likely that it comes from M1. Intensity dependence of the EEG responses to TMS Significant intensity dependence was found in several aspects of the response. Although the waveforms of the corresponding sources seemed qualitatively similar for left and right M1 stimulation, the only case for which significant intensity dependence was found for both stimulation sites was the positive phase of the MFMi response. In addition, with right M1 stimulation significant dependence was

found at the late phase of AMi and the cerebellar responses. For the other response peaks one should differentiate between the cases where there was no intensity dependence at all and cases where intensity dependence appeared to exist in the grand average (Fig. 5), but could not be shown statistically. The latter is true for the early peaks of the AMi and MFMi sources. This might have to do with residual artifact at this time period that attenuated in the grand average but still affected the statistical analysis. There seem to be several modes of intensity dependence. In some cases the response increased linearly with the stimulation strength (positive peak of the MFMi, positive peak of AMi in the right M1 stimulation), whereas in other cases there was a threshold phenomenon—a response only to high intensity stimulation. In some cases the threshold was between 100% and 120% RMT (AMi left, CBi left, Cbm left) and in other cases between 80% and 100% RMT (TPJi right, Cbm right). Although the limitations of the present study do not allow detailed statistical analysis of these phenomena, the ECD modeling of the response combined with statistical analysis of the source waveforms are a promising method for performing such detailed analyses in the future. The use of significance threshold uncorrected for multiple comparisons between sources was expected to result in some false positives and the cluster in Cbm between 11–14 ms seen for left M1 stimulation is likely to be such a case, since its latency was too short compared to other cerebellar responses. Relation with the results of other TMS–EEG studies Several studies published in recent years examined the multichannel-EEG responses to TMS of M1 (Ilmoniemi et al., 1997; Paus et al., 2001b; Komssi et al., 2002, 2004; Bonato et al., 2006). The studies differ with respect to the earliest time that was examined but, roughly, all report a dipolar pattern ipsilateral to stimulation with an anterior negativity and a posterior positivity between 10 and 20 ms after the stimulus (sometimes referred to as P15), the broad vertex-positive pattern around 30 ms after the stimulus (P30) and a dipolar pattern ipsilateral to stimulation with opposite polarity to that of P15 around 45 ms after the stimulus (N45). The results of these studies have been remarkably consistent and the same sequence of scalp patterns was also observed in the present study (Fig. 3). The novel contribution of the present study is in suggesting a detailed source model that improves our understanding of the mechanisms generating the scalp patterns described previously. In our model P15 and N45 are accounted for by the AMi source and P30 by the MFMi source with a small contribution of MFMc. In the original study based on the same series of experiments as the present study (Komssi et al., 2004) source reconstruction was performed using minimum-norm method (Hamalainen and Ilmoniemi, 1994). This method does not allow relating scalp patterns occurring at different times to the same source and analyzing the time course and stimulus intensity dependence of each source separately as was done in the present study. Paus et al. (2001b) discussed the possible origins of P30 and N45 and their conclusion was that only N45 originates in the motor cortex—consistent with our model. In a recent study from the same group (Van Der Werf and Paus, 2006) N45 amplitude was shown to be affected by 1 Hz rTMS delivered at M1 but not at dorsal premotor cortex, which again points at M1 as a likely generator. Another finding of Paus et al. (2001b) was that P30 amplitude does not correlate with TMS intensity. This contradicts our finding of a significant correlation of the MFMi response with stimulus

V. Litvak et al. / NeuroImage 37 (2007) 56–70

strength (Fig. 5). Perhaps the discrepancy originates in a difference in the set of stimulation intensities tested. Similarly to Komssi et al. (2002), contralateral responses were weak compared to the ipsilateral ones and the dipoles of these responses were mainly located close to the midline. Response was observed in the MFMc source primarily for the left M1 stimulation. However, the fact that it coincided in time with the MFMi response and both sources have overlapping scalp projections does not allow ruling out the possibility of crosstalk. We did not observe a contralateral dipolar pattern at about 24 ms post-stimulus as reported in Ilmoniemi et al. (1997). Such a pattern was also absent in several other studies (Paus et al., 2001b; Komssi et al., 2004; Bonato et al., 2006). Moreover, there is experimental evidence based on TMS indicating that the transcallosal connections are mainly inhibitory (Gerloff et al., 1998). Thus, the question whether responses in contralateral M1 can indeed be recorded at least in some subjects requires further clarification. Relation with the results of studies using other imaging methods Several previous studies examined the effect of TMS of M1 using imaging methods other than EEG. Overall there is great disparity between the results which might be related to differences in stimulation protocols and the imaging methods used. The common elements seen in all the studies are activation at the stimulation site and auditory activation due to the coil click. In the present study this would correspond to the AMi source and possibly the TPJ sources respectively. A question of special relevance to our results is whether there is activation that could correspond to our MFM sources i.e. relatively deep activation close to the midline. This is indeed found in the fMRI study of Li et al. (2004) where activation at the cingulate gyrus (Brodmann area 32) was found in response to single TMS pulses delivered at the left M1, in the study of Bestmann et al. (2004) where activation was found at both SMA and CMA in response to 3.125 Hz repetitive TMS and in the study of Baudewig et al. (2001) where SMA was activated by 10 Hz rTMS. Interestingly, the response in both areas was found to increase with TMS intensity (Bestmann et al., 2004) which would be consistent with the properties of the MFMi source in our study. Finally, activation in the cerebellum was also found in the study of Bestmann et al. (2004) and Baudewig et al. (2001). In the latter study cerebellar activation ipsilateral to stimulation was described, consistent with our results. For many of the areas found to be activated in fMRI studies there are no clear correlates in the EEG. This is not surprising given the difference in the experimental protocols and the types of brain activity that can be detected by the different imaging methods (Nunez and Silberstein, 2000). Nevertheless, some utility can be drawn from the comparison between the two kinds of results. Limitations of the present study The residual artifact The existence of a residual artifact in the data recorded with the sample-and-hold circuit is surprising since the main purpose of this circuit is TMS artifact suppression. However, this artifact has also been observed in previous studies using the same recording system (Paus et al., 2001b). Although in some cases the problem was not very severe because of low stimulation intensity and carefully chosen stimulation sites, the artifact is still a severe enough problem to limit the applicability of TMS–EEG combination.

67

There are several suggestions as to reasons for this artifact, although to our knowledge none has been experimentally tested. One of them is that the origin of the artifact is electrode movement caused by the TMS pulse (Virtanen et al., 1999). While this factor is a possible contributor for electrodes located in the immediate proximity to the coil, the electrodes distant from the stimulation site are unlikely to be affected. A very plausible theory put forward by Thut et al. (2005) distinguishes between the part of the artifact induced directly by the brief magnetic pulse and the long-lasting part attributed to the decay of electrical charges induced in the hardware components by the brief high-voltage peak. Since in our case the hardware components were isolated at the time of the pulse, it is more likely that charges of this sort accumulated at the interface between the electrodes and the skin. The time required for decay of this charge to physiological range was longer than the time of amplifier disconnection (7 ms) and therefore after the reconnection the still ongoing discharge was registered by the recording system as the artifact. This suggestion is supported by our observations that the artifact got larger as the electrode paste became dry during the measurement and that in the majority of cases the artifact time course (i.e. the time course corresponding to the artifact sources of the model in Stage II of the correction) could be well fitted by exponential decay (data not shown). The exact scalp topography of the artifact may vary between subjects and also as a result of repositioning the coil or head movements within one experiment. The primary reason for that is that while the position of the coil is optimized by examining motor evoked responses, the angle of the coil with respect to the electrodes is influenced by the head size and shape and local skull curvature under the coil. Even small differences in coil orientation can have profound effects on the effective magnetic field strength near a particular electrode and perhaps the exact angle of the electrode with respect to the spatial gradient of the field also makes a difference in the amount of charge that can be accumulated at the skin/electrode junction as a result of TMS pulse. Furthermore, the initial amplitude and time constant of the artifact are influenced by the impedance of the electrode/skin junction i.e. the local conductive properties of the subject's skin, the quality of impedance reduction during electrode application, the amount of conductive gel in the electrode and the degree to which the gel has dried out during the experiment. In the present study for each combination of stimulation site and intensity two data files were separately acquired with the coil being repositioned in between. Each file contained 50 stimuli and with random ISI's between 1.5 and 2.5 s (Komssi et al., 2004) was acquired in about 2 min. The position of the coil relative to the subject's head was fixed by a coil holder. Therefore, neither changes in angle nor drying out of gel could be an issue in such a short time period. Indeed it is our impression that the artifact topography stayed constant within a file. In contrast, we observed that the topography changed between the two files and representation of the artifact created based on one file could not be used for correcting the other file. This was probably due to minor differences in the coil angle. Furthermore, as mentioned above, the artifact amplitude increased over the whole experiment as the electrode paste became dry. All the sources of variability discussed above were accounted for in our analysis due to the fact that we created artifact representations separately for each data file while keeping the brain activity representation constant. This made it possible to accom-

68

V. Litvak et al. / NeuroImage 37 (2007) 56–70

modate variable artifact topographies and still get comparable brain source activities for different data files within and between subjects. The artifact correction method of Berg and Scherg (1994) works well when the artifact can be represented by a linear combination of a relatively small number of basis topographies. While this is usually true for the late decaying part of the artifact, the early part mostly blocked by the sample-and-hold circuit might contain additional rapidly decaying spatial components. For this reason small residual artifacts are also seen in the very early part of the corrected data (Fig. 3) between 7 and 10 ms. Thus, this part of the response was not analyzed. Furthermore, for the same reason it might not be possible to correct the early part of the artifact with the software method alone in the case when no sample-and-hold circuit is used. The correction method can handle artifact topographies that change over time within a trial as long as the artifact can be represented by a (time-varying) linear combination of the selected basis topographies. Theoretically, this cannot be guaranteed since the different magnetic field strength at electrodes with different distances to the coil might lead to different charges in the capacitors which are present in the biological and/or recording system. Oscillating circuits might present different eigenfrequencies, which might affect the initial artifact peaks in a different way than later time intervals. Our study could not completely rule out these possibilities. However, our iterative correction procedure guaranteed that there were no artifact residuals exceeding the physiological signal amplitudes in the epochs that were analyzed. In our view the fact that the artifact could be successfully corrected and the amplitude of the signals could be brought into physiological range using a small number of topographical components for representing the artifact is an empirical result indicating that the abovementioned property holds in our case. It is possible that due to variations in the artifact topography during the decay there were artifact residuals in individual data files. However, since (a) these residuals were small, in the physiological amplitude range and (b) their topographies can be assumed to be different and independent in each subject (since these effects largely depend on individual distribution of electrode impedances), it can be assumed that they averaged out in the grand average computed across all subjects and three stimulation intensities, on which our source model was based. We hypothesize that the artifact residuals left in individual subjects could be the reason for the fact that no significant intensity dependence could be shown for early AMi and MFMi peaks. The practical conclusion from the above discussion is that averaging of corrected evoked responses over multiple subjects or multiple experimental sessions in each subject is an essential component of the correction procedure. Reduction of the artifact cannot be achieved effectively by averaging alone even with a large number of subjects because the topographies of uncorrected residual artifacts depend primarily on the stimulation site and are therefore expected to be highly correlated across subjects, whereas topographies of residuals left after correction depend on distribution of electrode impedances and are therefore expected to be uncorrelated. Possible bias in source locations Since electrode positions for each subject were not measured in the original experiments, the registration of the electrode cap to the brain anatomy might be imprecise in the present study. Thus caution should be taken when interpreting model source coordinates and

possible bias, especially in the anterior–posterior direction should be taken into account. Modeling approach In the present study we used an equivalent current dipole model using a sequential fitting strategy. This approach is only one of a variety of methods for solving the so-called ‘inverse problem’ i.e. computing the current source density in the cortex from EEG data recorded on the scalp (Dale and Sereno, 1993; Hamalainen and Ilmoniemi, 1994; Pascual Marqui et al., 1994; Scherg and Berg, 1996; Gross et al., 2001). Since the problem is underdetermined, the solution is not well defined by the data and additional assumptions are necessary. Some of the source localization methods compute the current source density from the scalp potentials according to a pre-defined formula based on method-specific assumptions. Such methods have the advantage of objectivity because they do not require that the experimenter manually constructs a source model. However, in the case of EEG responses to TMS they might produce erroneous results since they do not differentiate between the physiological cortical response and the TMS artifact. The process of manually constructing a multiple source model makes it possible to keep track of the relation between the scalp patterns in the raw data and the model sources accounting for these patterns and to avoid modeling patterns that are clearly artifacts. In the particular case of the present study, the symmetry of responses to the left and right M1 stimulation that clearly emerged from the grand average data (Fig. 3) could be used as a constraint in the model building. Taking this symmetry into account as a constraint for a solution such as minimum norm (Dale and Sereno, 1993; Hamalainen and Ilmoniemi, 1994) or LORETA (Pascual Marqui et al., 1994) is a highly non-trivial task and would not have general applicability. An additional subjective decision that should be made while building a multiple ECD model is about the number of sources. While this parameter is not unambiguously determined by the data, the decision can be based on examining the residual activity not accounted for by the initial model with only one source or two symmetric sources. Additional sources are added as long as the residual contains clear peaks. The advantage of this process is the ability to roughly estimate the number of distinct generators involved in the response. Such information cannot be readily obtained from current source density estimation as discussed above. Conclusions The methodology developed in the present study made it possible to substantially reduce the TMS stimulation artifact even for high stimulation intensities. This paves the way to applying TMS–EEG imaging to a variety of stimulation sites and conditions many of which were not accessible until now. Based on artifactcorrected data we propose for the first time a multiple source model for the early phase of the response evoked by magnetic stimulation of left and right M1. A remarkable symmetry was found between left and right stimulation. The main sources contributing to the response according to our model are the stimulated motor cortex and a deep ipsilateral source close to the midline, possibly cingulate gyrus or SMA. We also discovered responses in the cerebellum which have not been described previously. The results of the study were largely consistent with previous EEG and fMRI studies.

V. Litvak et al. / NeuroImage 37 (2007) 56–70

Acknowledgments VL was funded by Technion stipend for graduate students, by the German Academic Exchange Service (DAAD) and by Boehringer Ingelheim Fonds. We thank Robert Oostenveld and Eric Maris for their assistance with nonparametric statistical analysis and Naomi Bleich for graphical work. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.neuroimage.2007.05.015. References Amedi, A., Malach, R., Pascual-Leone, A., 2005. Negative BOLD differentiates visual imagery and perception. Neuron 48, 859–872. Antal, A., Nitsche, M.A., Kincses, T.Z., Lampe, C., Paulus, W., 2004. No correlation between moving phosphene and motor thresholds: a transcranial magnetic stimulation study. NeuroReport 15, 297–302. Barker, A.T., Jalinous, R., Freeston, I.L., 1985. Non-invasive magnetic stimulation of human motor cortex. Lancet 1, 1106–1107. Baudewig, J., Siebner, H.R., Bestmann, S., Tergau, F., Tings, T., Paulus, W., Frahm, J., 2001. Functional MRI of cortical activations induced by transcranial magnetic stimulation (TMS). NeuroReport 12, 3543–3548. Baumer, T., Bock, F., Koch, G., Lange, R., Rothwell, J.C., Siebner, H.R., Munchau, A., 2006. Magnetic stimulation of human premotor or motor cortex produces interhemispheric facilitation through distinct pathways. J. Physiol. 572, 857–868. Bender, S., Basseler, K., Sebastian, I., Resch, F., Kammer, T., Oelkers-Ax, R., Weisbrod, M., 2005. Electroencephalographic response to transcranial magnetic stimulation in children: Evidence for giant inhibitory potentials. Ann. Neurol. 58, 58–67. Berg, P., Scherg, M., 1994. A multiple source approach to the correction of eye artifacts. Electroencephalogr. Clin. Neurophysiol. 90, 229–241. Bestmann, S., Baudewig, J., Siebner, H.R., Rothwell, J.C., Frahm, J., 2004. Functional MRI of the immediate impact of transcranial magnetic stimulation on cortical and subcortical motor circuits. Eur. J. Neurosci. 19, 1950–1962. Bestmann, S., Baudewig, J., Siebner, H.R., Rothwell, J.C., Frahm, J., 2005. BOLD MRI responses to repetitive TMS over human dorsal premotor cortex. NeuroImage 28, 22–29. Bohning, D.E., Shastri, A., McConnell, K.A., Nahas, Z., Lorberbaum, J.P., Roberts, D.R., Teneback, C., Vincent, D.J., George, M.S., 1999. A combined TMS/fMRI study of intensity-dependent TMS over motor cortex. Biol. Psychiatry 45, 385–394. Bohning, D.E., Shastri, A., Wassermann, E.M., Ziemann, U., Lorberbaum, J.P., Nahas, Z., Lomarev, M.P., George, M.S., 2000. BOLD-f MRI response to single-pulse transcranial magnetic stimulation (TMS). J. Magn. Reson. Imaging 11, 569–574. Bonato, C., Miniussi, C., Rossini, P.M., 2006. Transcranial magnetic stimulation and cortical evoked potentials: a TMS/EEG co-registration study. Clin. Neurophysiol. 117 (8), 1699–1707. Cerri, G., Shimazu, H., Maier, M.A., Lemon, R.N., 2003. Facilitation from ventral premotor cortex of primary motor cortex outputs to macaque hand muscles. J. Neurophysiol. 90, 832–842. Chen, R., Classen, J., Gerloff, C., Celnik, P., Wassermann, E.M., Hallett, M., Cohen, L.G., 1997. Depression of motor cortex excitability by lowfrequency transcranial magnetic stimulation. Neurology 48, 1398–1403. Civardi, C., Cantello, R., Asselman, P., Rothwell, J.C., 2001. Transcranial magnetic stimulation can be used to test connections to primary motor areas from frontal and medial cortex in humans. NeuroImage 14, 1444–1453. Dale, A.M., Sereno, M.I., 1993. Improved localization of cortical activity by

69

combining EEG and MEG with MRI cortical surface reconstruction: a linear approach. J. Cogn. Neurosci. 5, 162–176. Daskalakis, Z.J., Paradiso, G.O., Christensen, B.K., Fitzgerald, P.B., Gunraj, C., Chen, R., 2004. Exploring the connectivity between the cerebellum and motor cortex in humans. J. Physiol. 557, 689–700. Denslow, S., Lomarev, M., George, M.S., Bohning, D.E., 2005. Cortical and subcortical brain effects of transcranial magnetic stimulation (TMS)induced movement: an interleaved TMS/functional magnetic resonance imaging study. Biol. Psychiatry 57, 752–760. Esser, S.K., Huber, R., Massimini, M., Peterson, M.J., Ferrarelli, F., Tononi, G., 2006. A direct demonstration of cortical LTP in humans: a combined TMS/EEG study. Brain Res. Bull. 69, 86–94. Fox, P., Ingham, R., George, M.S., Mayberg, H., Ingham, J., Roby, J., Martin, C., Jerabek, P., 1997. Imaging human intra-cerebral connectivity by PET during TMS. NeuroReport 8, 2787–2791. Fox, P.T., Narayana, S., Tandon, N., Fox, S.P., Sandoval, H., Kochunov, P., Capaday, C., Lancaster, J.L., 2005. Intensity modulation of TMSinduced cortical excitation: Primary motor cortex. Hum. Brain Mapp. 27 (6), 478–487. Fuggetta, G., Pavone, E.F., Walsh, V., Kiss, M., Eimer, M., 2006. Corticocortical interactions in spatial attention: a combined ERP/TMS study. J. Neurophysiol. 95 (5), 3277–3280. Gerloff, C., Cohen, L.G., Floeter, M.K., Chen, R., Corwell, B., Hallett, M., 1998. Inhibitory influence of the ipsilateral motor cortex on responses to stimulation of the human cortex and pyramidal tract. J. Physiol. 510 (Pt. 1), 249–259. Gross, J., Kujala, J., Hamalainen, M., Timmermann, L., Schnitzler, A., Salmelin, R., 2001. Dynamic imaging of coherent sources: studying neural interactions in the human brain. Proc. Natl. Acad. Sci. U. S. A. 98, 694–699. Haaland, K.Y., Harrington, D.L., 1996. Hemispheric asymmetry of movement. Curr. Opin. Neurobiol. 6, 796–800. Hallett, M., 2000. Transcranial magnetic stimulation and the human brain. Nature 406, 147–150. Hamalainen, M.S., Ilmoniemi, R.J., 1994. Interpreting magnetic fields of the brain: minimum norm estimates. Med. Biol. Eng. Comput. 32, 35–42. Hoechstetter, K., Bornfleth, H., Weckesser, D., Ille, N., Berg, P., Scherg, M., 2004. BESA source coherence: a new method to study cortical oscillatory coupling. Brain Topogr. 16, 233–238. Huang, Y.Z., Edwards, M.J., Rounis, E., Bhatia, K.P., Rothwell, J.C., 2005. Theta burst stimulation of the human motor cortex. Neuron 45, 201–206. Ille, N., Berg, P., Scherg, M., 2002. Artifact correction of the ongoing EEG using spatial filters based on artifact and brain signal topographies. J. Clin. Neurophysiol. 19, 113–124. Ilmoniemi, R.J., Virtanen, J., Ruohonen, J., Karhu, J., Aronen, H.J., Naatanen, R., Katila, T., 1997. Neuronal responses to magnetic stimulation reveal cortical reactivity and connectivity. NeuroReport 8, 3537–3540. Ilmoniemi, R.J., Ruohonen, J., Virtanen, J., Aronen, H.J., Karhu, J., 1999. EEG responses evoked by transcranial magnetic stimulation. Electroencephalogr. Clin. Neurophysiol. Suppl. 51, 22–29. Kahkonen, S., Wilenius, J., Komssi, S., Ilmoniemi, R.J., 2004. Distinct differences in cortical reactivity of motor and prefrontal cortices to magnetic stimulation. Clin. Neurophysiol. 115, 583–588. Kahkonen, S., Komssi, S., Wilenius, J., Ilmoniemi, R.J., 2005a. Prefrontal TMS produces smaller EEG responses than motor-cortex TMS: implications for rTMS treatment in depression. Psychopharmacology (Berl.) 181, 16–20. Kahkonen, S., Komssi, S., Wilenius, J., Ilmoniemi, R.J., 2005b. Prefrontal transcranial magnetic stimulation produces intensity-dependent EEG responses in humans. NeuroImage 24, 955–960. Kandel, E.R., Schwartz, J.H., Jessell, T.M., 2000. Principles of Neural Science, 4th Edition. McGraw-Hill Health Professions Division, New York. Klein, E., Kreinin, I., Chistyakov, A., Koren, D., Mecz, L., Marmur, S., BenShachar, D., Feinsod, M., 1999. Therapeutic efficacy of right prefrontal slow repetitive transcranial magnetic stimulation in major depression: a double-blind controlled study. Arch. Gen. Psychiatry 56, 315–320.

70

V. Litvak et al. / NeuroImage 37 (2007) 56–70

Koch, G., Franca, M., Del Olmo, M.F., Cheeran, B., Milton, R., Alvarez Sauco, M., Rothwell, J.C., 2006. Time course of functional connectivity between dorsal premotor and contralateral motor cortex during movement selection. J. Neurosci. 26, 7452–7459. Komssi, S., Kahkonen, S., 2006. The novelty value of the combined use of electroencephalography and transcranial magnetic stimulation for neuroscience research. Brain Res. Brain Res. Rev. 52, 183–192. Komssi, S., Aronen, H.J., Huttunen, J., Kesaniemi, M., Soinne, L., Nikouline, V.V., Ollikainen, M., Roine, R.O., Karhu, J., Savolainen, S., Ilmoniemi, R.J., 2002. Ipsi- and contralateral EEG reactions to transcranial magnetic stimulation. Clin. Neurophysiol. 113, 175–184. Komssi, S., Kahkonen, S., Ilmoniemi, R.J., 2004. The effect of stimulus intensity on brain responses evoked by transcranial magnetic stimulation. Hum. Brain Mapp. 21, 154–164. Li, X., Teneback, C.C., Nahas, Z., Kozel, F.A., Large, C., Cohn, J., Bohning, D.E., George, M.S., 2004. Interleaved transcranial magnetic stimulation/ functional MRI confirms that lamotrigine inhibits cortical excitability in healthy young men. Neuropsychopharmacology 29, 1395–1407. Lorch Jr., R.F., Myers, J.L., 1990. Regression analyses of repeated measures data in cognitive research. J. Exp. Psychol., Learn. Mem. Cogn. 16, 149–157. Maris, E., Oostenveld, R. (in press). Nonparametric Statistical Testing of EEG- and MEG-data. J. Neurosci. Methods. Massimini, M., Ferrarelli, F., Huber, R., Esser, S.K., Singh, H., Tononi, G., 2005a. Breakdown of cortical effective connectivity during sleep. Science 309, 2228–2232. Massimini, M., Ferrarelli, F., Peterson, M.J., Esser, S.K., Huber, R., Lazar, M., Alexander, A., Tononi, G, 2005b. Cortical EEG response to TMS in schizophrenia patients and normal controls. In: Society for Neuroscience, pp Program No. 443.445. Washington, DC. McGill, R., Tukey, J.W., Larsen, W.A., 1978. Variations of Boxplots. Am. Stat. 32, 12–16. Michel, C.M., Murray, M.M., Lantz, G., Gonzalez, S., Spinelli, L., Grave de Peralta, R., 2004. EEG source imaging. Clin. Neurophysiol. 115, 2195–2222. Mitzdorf, U., 1985. Current source-density method and application in cat cerebral cortex: investigation of evoked potentials and EEG phenomena. Physiol. Rev. 65, 37–100. Nahas, Z., Lomarev, M., Roberts, D.R., Shastri, A., Lorberbaum, J.P., Teneback, C., McConnell, K., Vincent, D.J., Li, X., George, M.S., Bohning, D.E., 2001. Unilateral left prefrontal transcranial magnetic stimulation (TMS) produces intensity-dependent bilateral effects as measured by interleaved BOLD fMRI. Biol. Psychiatry 50, 712–720. Nichols, T.E., Holmes, A.P., 2002. Nonparametric permutation tests for functional neuroimaging: a primer with examples. Hum. Brain Mapp. 15, 1–25. Nikulin, V.V., Kicic, D., Kahkonen, S., Ilmoniemi, R.J., 2003. Modulation of electroencephalographic responses to transcranial magnetic stimulation: evidence for changes in cortical excitability related to movement. Eur. J. Neurosci. 18, 1206–1212. Nunez, P.L., Silberstein, R.B., 2000. On the relationship of synaptic activity to macroscopic measurements: does co-registration of EEG with fMRI make sense? Brain Topogr. 13, 79–96. Osipova, D., Takashima, A., Oostenveld, R., Fernandez, G., Maris, E., Jensen, O., 2006. Theta and gamma oscillations predict encoding and retrieval of declarative memory. J. Neurosci. 26, 7523–7531. Pascual-Leone, A., Valls-Sole, J., Wassermann, E.M., Hallett, M., 1994. Responses to rapid-rate transcranial magnetic stimulation of the human motor cortex. Brain 117 (Pt. 4), 847–858. Pascual-Leone, A., Rubio, B., Pallardo, F., Catala, M.D., 1996. Rapid-rate transcranial magnetic stimulation of left dorsolateral prefrontal cortex in drug-resistant depression. Lancet 348, 233–237.

Pascual Marqui, R.D., Michel, C.M., Lehmann, D., 1994. Low resolution electromagnetic tomography: a new method for localizing electrical activity in the brain. Int. J. Psychophysiol. 18, 49–65. Paus, T., Castro-Alamancos, M.A., Petrides, M., 2001a. Cortico-cortical connectivity of the human mid-dorsolateral frontal cortex and its modulation by repetitive transcranial magnetic stimulation. Eur. J. Neurosci. 14, 1405–1411. Paus, T., Sipila, P.K., Strafella, A.P., 2001b. Synchronization of neuronal activity in the human primary motor cortex by transcranial magnetic stimulation: an EEG study. J. Neurophysiol. 86, 1983–1990. Ramnani, N., 2006. The primate cortico-cerebellar system: anatomy and function. Nat. Rev., Neurosci. 7, 511–522. Ruohonen, J., Ilmoniemi, R.J., 1999. Modeling of the stimulating field generation in TMS. Electroencephalogr. Clin. Neurophysiol. Suppl. 51, 30–40. Scherg, M., Berg, P., 1996. New concepts of brain source imaging and localization. Electroencephalogr. Clin. Neurophysiol. Suppl. 46, 127–137. Scherg, M., Von Cramon, D., 1986. Evoked dipole source potentials of the human auditory cortex. Electroencephalogr. Clin. Neurophysiol. 65, 344–360. Scherg, M., Ille, N., Bornfleth, H., Berg, P., 2002. Advanced tools for digital EEG review: virtual source montages, whole-head mapping, correlation, and phase analysis. J. Clin. Neurophysiol. 19, 91–112. Scherg, M., Hoechstetter, K., 2004. BESA—Brain Electrical Source Analysis. Software manual, Graefelfing/Munich. Schnitzler, A., Timmermann, L., Gross, J., 2006. Physiological and pathological oscillatory networks in the human motor system. J. Physiol. (Paris) 99, 3–7. Shimazu, H., Maier, M.A., Cerri, G., Kirkwood, P.A., Lemon, R.N., 2004. Macaque ventral premotor cortex exerts powerful facilitation of motor cortex outputs to upper limb motoneurons. J. Neurosci. 24, 1200–1211. Takashima, A., Jensen, O., Oostenveld, R., Maris, E., van de Coevering, M., Fernandez, G., 2006. Successful declarative memory formation is associated with ongoing activity during encoding in a distributed neocortical network related to working memory: a magnetoencephalography study. Neuroscience 139, 291–297. Talairach, J., Tournoux, P., 1988. Co-planar Stereotaxic Atlas of the Human Brain: an Approach to Medical Cerebral Imaging. G. Thieme; Thieme Medical Publishers, Stuttgart. Thielscher, A., Kammer, T., 2002. Linking physics with physiology in TMS: a sphere field model to determine the cortical stimulation site in TMS. NeuroImage 17, 1117–1130. Thut, G., Northoff, G., Ives, J.R., Kamitani, Y., Pfennig, A., Kampmann, F., Schomer, D.L., Pascual-Leone, A., 2003. Effects of single-pulse transcranial magnetic stimulation (TMS) on functional brain activity: a combined event-related TMS and evoked potential study. Clin. Neurophysiol. 114, 2071–2080. Thut, G., Ives, J.R., Kampmann, F., Pastor, M.A., Pascual-Leone, A., 2005. A new device and protocol for combining TMS and online recordings of EEG and evoked potentials. J. Neurosci. Methods 141, 207–217. Van Der Werf, Y.D., Paus, T., 2006. The neural response to transcranial magnetic stimulation of the human motor cortex. I. Intracortical and cortico-cortical contributions. Exp. Brain Res. 175 (2), 231–245. Virtanen, J., Ruohonen, J., Naatanen, R., Ilmoniemi, R.J., 1999. Instrumentation for the measurement of electric brain responses to transcranial magnetic stimulation. Med. Biol. Eng. Comput. 37, 322–326. Wolters, A., Sandbrink, F., Schlottmann, A., Kunesch, E., Stefan, K., Cohen, L.G., Benecke, R., Classen, J., 2003. A temporally asymmetric Hebbian rule governing plasticity in the human motor cortex. J. Neurophysiol. 89, 2339–2345.