Automatized set-up procedure for transcranial magnetic stimulation protocols

NeuroImage 153 (2017) 307–318

Contents lists available at ScienceDirect

NeuroImage journal homepage: www.elsevier.com/locate/neuroimage

Automatized set-up procedure for transcranial magnetic stimulation protocols

MARK

⁎

S. Harquela,b,c,d, , J. Diarda,b, E. Raffina,c, B. Passeraa,b, G. Dall'Ignae,c, C. Marendaza,b, O. Davida,c, A. Chauvina,b a

Univ. Grenoble Alpes, F-38000 Grenoble, France CNRS, UMR 5105, Laboratoire de Psychologie et de Neurocognition, LPNC, F-38000 Grenoble, France c INSERM, U1216, Grenoble Institut des Neurosciences, GIN, F-38000 Grenoble, France d CNRS, INSERM, UMS 3552, IRMaGe, F-38000 Grenoble, France e Pôle de Psychiatrie et Neurologie, Centre Hospitalier Universitaire Grenoble Alpes, F-38000 Grenoble, France b

A R T I C L E I N F O

A BS T RAC T

Keywords: Robotized TMS Hotspot hunting Cortical excitability MEP Motor mapping

Transcranial Magnetic Stimulation (TMS) established itself as a powerful technique for probing and treating the human brain. Major technological evolutions, such as neuronavigation and robotized systems, have continuously increased the spatial reliability and reproducibility of TMS, by minimizing the influence of human and experimental factors. However, there is still a lack of efficient set-up procedure, which prevents the automation of TMS protocols. For example, the set-up procedure for defining the stimulation intensity specific to each subject is classically done manually by experienced practitioners, by assessing the motor cortical excitability level over the motor hotspot (HS) of a targeted muscle. This is time-consuming and introduces experimental variability. Therefore, we developed a probabilistic Bayesian model (AutoHS) that automatically identifies the HS position. Using virtual and real experiments, we compared the efficacy of the manual and automated procedures. AutoHS appeared to be more reproducible, faster, and at least as reliable as classical manual procedures. By combining AutoHS with robotized TMS and automated motor threshold estimation methods, our approach constitutes the first fully automated set-up procedure for TMS protocols. The use of this procedure decreases inter-experimenter variability while facilitating the handling of TMS protocols used for research and clinical routine.

1. Introduction Transcranial Magnetic Stimulation (TMS) is a non-invasive brain stimulation technique (Barker et al., 1985; Hallett, 2000). Applied alone or coupled with other neuroimaging techniques (Bestmann and Feredoes, 2013), its application are now numerous in both fundamental (Rogasch and Fitzgerald, 2013; Bortoletto et al., 2015) and clinical research (Ragazzoni et al., 2013; Lefaucheur et al., 2014; Lefaucheur and Picht, 2016). In order to standardize procedures and consequently reduce intersubject variability in response to TMS, the field has recently embraced major technological evolutions. Neuronavigation systems dedicated to TMS (Herwig et al., 2001) significantly improved its spatial precision and reproducibility (Julkunen et al., 2009; Weiss et al., 2013) and TMS-robotized systems enabled the automation of coil positioning (Finke et al., 2008; Kantelhardt et al., 2009; Ginhoux et al., 2013). In addition to improving spatial precision and reproducibility compared

⁎

to manual positioning (Ginhoux et al., 2013), robotized TMS paves the way for new acquisition protocols, such as automated cortical mapping procedures (Harquel et al., 2016a). It is thus likely that the future of TMS resides in the full automation of protocols, partly enabled by robotics. Every TMS protocol begins by a mandatory set-up procedure, which aims at defining the stimulation intensity to be employed on the cortical target (Rossi et al., 2009; Herbsman et al., 2009; Wassermann and Epstein, 2012). This intensity has to be defined specifically for each subject because it depends on individual neuroanatomy. The procedure consists in assessing the resting (or active) motor threshold (rMT, or aMT) by measuring the muscular activity evoked by the stimulation of the motor hotspot (HS) over the primary motor cortex (M1). Stimulation intensities are then most often expressed as a percentage of this threshold, in order to conform to safety guidances and to standardize stimulation power between subjects (Herbsman et al., 2009).

Correspondence to: Laboratoire de Psychologie et de Neurocognition, CNRS UMR 5105, Université Grenoble Alpes, BSHM, BP47, 38040 Grenoble Cedex 9, France. E-mail addresses: [email protected] (S. Harquel).

http://dx.doi.org/10.1016/j.neuroimage.2017.04.001 Received 16 November 2016; Accepted 1 April 2017 Available online 05 April 2017 1053-8119/ © 2017 Elsevier Inc. All rights reserved.

NeuroImage 153 (2017) 307–318

S. Harquel et al.

Fig. 1. Main components and hypotheses of the AutoHS model. a: stimulation grid used for hunting and for expressing point coordinates. b: Probabilistic prior concerning HS position, centered on the hand knob of the primary motor cortex M1 (precentral gyrus). c: Probability distribution of MEP amplitudes, estimated on N EMG recordings. d: MEP mean amplitude (μi) modulation as a function of the distances dX and dY between Si and HS along the X and Y-axis respectively.

the HS hunting procedure is the first interosseous muscle (FDI).

Depending on its definitions (Meincke et al., 2016), the HS corresponds to the cortical target over M1 where TMS evokes the lowest MT (Rossini et al., 1994), or the largest motor evoked potentials (MEPs) on the targeted muscle (van de Ruit et al., 2015). While efficient automated MT estimation methods have already been developed and are used since then (Awiszus, 2003; Awiszus and Borckardt, 2011), it is not the case for HS hunting. In practice, the HS position is manually set by experienced practitioners, even though Meincke et al. (2016) recently developed the first automated HS hunting procedure using the mapping of MTs. This method appeared to be effective in automatically assessing the HS position, and producing insightful data for motor mapping protocols. However, its substantial duration (over one hour) prevents its practical use in clinical settings and in most TMS experiments. Although quicker (about ten minutes), manual set-up procedures also have limitations: i. they represent an additional source of inter-experimenter variability (Gugino et al., 2001; Cincotta et al., 2010; Sollmann et al., 2013), ii. they require well-trained TMS practitioners, and iii. they rely on the observation of MEP mean amplitudes which are highly variable (Wassermann, 2002; Jung et al., 2010). In order to overcome these limitations, we propose here an automated HS hunting procedure (AutoHS) based on a Bayesian model. AutoHS aims at localizing the HS in a faster, more reliable and more reproducible way compared to a manual HS hunting performed by TMS experimenters. The present paper describes how we implemented HS hunting in a Bayesian model, and how we tested our method on virtual data and validated it against manual HS hunting performed by four TMS experts on 19 healthy volunteers. We finally discuss our method and its relevance for progressing towards fully automated set-up procedures for TMS protocols.

2.1. Bayesian model of AutoHS 2.1.1. Overview AutoHS is a probabilistic procedure, as is classical in the domain of multisensor data fusion in robotics (Bessière and Lebeltel, 2008). Its objective is to estimate the HS position, using the history of stimulated sites and recorded MEP amplitudes. AutoHS is built in two steps, applying the Bayesian programming methodology (Bessière et al., 2013). The first step consists in defining the joint probability distribution over five variables, including the HS position. From this joint probability distribution, the second step consists in applying Bayesian inference, so as to compute the probability distribution over HS positions, conditioned on previous observations. This procedure is complemented by a “smart” prospective method, that considers the most promising next cortical position to be stimulated, in terms of information gain (Baek et al., 2016), in order to find the HS as fast as possible. AutoHS automatically stops and settles on the HS position once enough information has been gathered. The method was implemented using Matlab (The Mathworks Inc., USA) and was run concurrently with the neuronavigation and the EMG recording systems. The default values of the model variables used in this work are reported throughout this section. These values have been estimated during pre-tests conducted on real motor mapping datasets not shown in this report. Their robustness are discussed later (see Discussion). 2.1.2. Model description The structure of the Bayesian model is defined by specifying the joint probability distribution over the five following variables.

2. Materials and methods

HS We describe first the Bayesian model of AutoHS in detail and second the experimental procedure used to test and compare AutoHS to manual methodology. Throughout this work, the targeted muscle for 308

represents the spatial position of the HS. The model makes the assumption that the 3D stimulation space can be projected to a 2D stimulation grid placed on the scalp surface over the motor cortex (Fig. 1a). HS is thus a two dimensional

NeuroImage 153 (2017) 307–318

S. Harquel et al.

P (HSμHS σsp S1: T μ1: T ),

variable expressed in a Cartesian coordinate system. Its x and y values are continuous over the stimulation grid area, so that (xHS , yHS ) ∈ ([1; Nx ], [1; Ny]), with Nx and Ny being the number of discrete points in the stimulation grid. This work uses a square grid, so that Nx = Ny = 7. The spacing between points is set to 7 mm, leading to a 4.2 cm sided square. This spacing affords a good spatial resolution regarding cortical mapping, while still making sure to be twice above navigation precision (around 3 mm, including camera measurement uncertainty and MRI registration error).AutoHS makes the initial assumption that the FDI HS is located on the hand knob of the precentral gyrus (primary motor cortex, M1). The prior probability distribution of HS over the stimulation space is thus a 2D-Gaussian distribution centered on this location:

⎛ ⎛ ( y − yHK )2 ⎞ ⎞⎟ (x − xHK )2 ⎟⎟ , P (HS ) ∝ exp ⎜⎜ −⎜⎜ + ⎟ 2 2 2σprior ⎠⎠ ⎝ ⎝ 2σprior

Si

μi

P (μi ) ∝

μHS

T

P (HSμHS σsp S1: T μ1: T ) = P (HS ) P (μHS ) P (σsp )

T

∏ P (S1: T ) ∏ P (μi |A), i =1

i =1

(5)

where A, for notational simplicity, is the conjunction of variables HS , μHS , σsp and Si. The priors on μHS and σsp are kept uninformative because of high inter-subject variability regarding these values. P (μHS ) and P (σsp ) thus are uniform distributions over their respective domains. The prior on the stimulation position Si is a Dirac distribution centered on the position selected during the previous iteration (see 2.1.5). S1 is initialized on the center of the stimulation grid. Finally, P (μi |A) corresponds to the probability of measuring the MEP mean amplitude μi given the coordinates of the stimulated site Si and HS , the estimated maximal amplitude μHS and spatial spread σsp of the muscle representation. We can write our knowledge μi given A in P (μi |A), by merging Eqs. (2) and (3), so that:

(1)

P (μi |A) ∝

⎛ (μ − μ ) 2 ⎞ 1 i exp ⎜ − ⎟ σi2 σi 2π ⎝ ⎠

(6)

⎛ ⎛ ⎛ ⎛ 2 ⎞⎞⎞ ⎞ 2 ⎜ ⎜μ − μ exp ⎜ −⎜ (xi − xHS ) + ( yi − yHS ) ⎟ ⎟ ⎟ ⎟ HS ⎟⎟⎟ ⎟ ⎜ ⎜ ⎜ ⎜ 2σsp2 2σsp2 ⎠⎠⎠ ⎟ ⎝ ⎝ ⎝ 1 exp ⎜⎜ − ∝ ⎟. σi2 σi 2π ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ 2.1.3. Question and inference The objective of AutoHS is to estimate the position of the HS, given all the MEP mean amplitudes μ1: T observed after stimulating cortical positions S1: T . However, this is not possible without jointly inferring μHS and σsp , since these variables are also conditioning the probability P (μi |A). Therefore, the joint probability P (HSμHS σsp |μ1: T S1: T ) is computed, which is a four-dimensional probability distribution. Bayesian inference yields (see appendix for the complete derivation):

(2)

where σi is the standard deviation measured over the N recorded MEP amplitudes (Fig. 1c). This work uses 5 MEP recordings per stimulation point, so that N=5. If no MEP was induced after 3 consecutive stimulations, μi was set to 0 and the model moved to the next iteration. represents the expected mean amplitude measured on the HS. Here, in accordance with the HS definition (see 1), AutoHS makes the assumption that μHS represents the maximal mean MEP amplitude that can be measured over the whole simulation space. μHS is a discrete variable over the domain of possible MEP max amplitudes Dmax (in μV). Here, we set Dmax = {250, 500, …, 4000}μV .Moreover, the model assumes that MEP amplitudes μi decrease as a function of the distance between the stimulated position Si and the HS position HS (Fig. 1d):

⎛ ⎛ ( y − y )2 ⎞ ⎞ (x − x ) 2 μi = μHS exp ⎜⎜ −⎜⎜ i 2HS + i 2HS ⎟⎟ ⎟⎟ . 2σsp 2σsp ⎠⎠ ⎝ ⎝ σsp

where S1: T and μ1: T denotes the series of variables Si and μi from the first to the last stimulated position. The joint probability distribution is simplified, by the appropriate conditional hypotheses assumptions, in order to obtain a product of simpler probability distributions:

where xHK and yHK are the coordinates of the hand knob, and σprior is the spatial spread of the distribution that quantifies the prior precision (Fig. 1b). We use here σprior = 2 (in grid units), leading to a spatial spread of 1.4 cm. is the position of the coil at the ith stimulation, i going from 1 to the total number of stimulated positions T. Si is also a two dimensional variable expressed in a Cartesian coordinate system. However, its x and y values are discrete and correspond to the nodes of the stimulation grid, so that (xi , yi) ∈ ({1, 2, …, Nx}, {1, 2, …, Ny}). is the mean amplitude observed from N MEPs on Si. μi is a discrete variable over the domain of possible MEP ampliDMEP tudes (in μV). Here, we set DMEP = {10, 20, …, 4000}μV .Here, μi is meant to follow a Gaussian distribution such that

⎛ (μ − μ ) 2 ⎞ 1 i exp ⎜ − ⎟, σi2 σi 2π ⎝ ⎠

(4)

T

P (HSμHS σsp |S1: T μ1: T ) ∝ P (HS )

∏ P (μi |A).

(7)

i =1

Once this distribution is evaluated, we compute its center of gravity, and the estimates for the HS x and y dimensions are its first two coordinates. Because P (HSμHS σsp |μ1: T S1: T ) is four dimensional, it is not easily displayed. For visualization, we thus compute and show the marginal distributions over each variable, that can be derived from Eq. (7): T

P (HS|μHS σsp S1: T μ1: T ) ∝ P (HS )

∑∑∏ μHS

σsp

P (μi |A)

i =1 T

P (μHS |HSσsp S1: T μ1: T ) ∝ P (HS )

∑∑∏ HS

(3)

σsp

P (μi |A)

i =1

(8)

T

is the standard deviation parameter of the 2D Gaussian pattern for MEP amplitude modulation (Eq. (3)). It represents the spatial spread of the muscle representation on the stimulation grid (Fig. 1d). σsp is a discrete variable over the domain of possible muscle representation sizes Dsp (in grid units). Here, we set Dsp = {0.2, 0.4, …, 2}.

P (σsp |μHS HSS1: T μ1: T ) ∝ P (HS )

∑ ∑ ∏ P (μi |A) μHS HS i =1

2.1.4. Stop criteria If the inference on HS can be done after each stimulation i, the total number of stimulations T is a priori unknown, and depends on the information gained during the whole process. The level of information at the ith iteration is assessed through the entropy hi of the joint distribution inferred in Eq. (7):

The joint probability distribution over the five above variables is then defined as follows: 309

NeuroImage 153 (2017) 307–318

S. Harquel et al.

hi = − ∑ ∑ ∑ P (HSμHS σsp |S1: i μ1: i ) ln (P (HSμHS σsp |S1: i μ1: i )).

2.2.2. Protocol design Subjects went through two stimulation sessions, separated by a delay of one week. Each session consisted of the two tested TMS set-up procedures: manual and automatized. The order of the procedures was counterbalanced between subjects and sessions. Each procedure consisted in one hotspot-hunting step followed by a resting motor threshold (rMT) assessment. Ten baseline MEPs were also recorded using an intensity of 120% rMT, in order to asses the influence of the procedures on its variability. In one of the two sessions, an extensive motor mapping was done either before, between or after the two procedures.

(9)

HS μHS σsp

By definition, the entropy hi starts at its theoretical maximum hmax and decreases throughout the HS hunting process as the model gains information. An increase of entropy during the process would incidentally reveal the violation of at least one of the model hypothesis. A first stop criterion can then be set up by defining a threshold level to be reached by hi. This threshold is expressed as a proportion of hmax. The process thus stops after the ith iteration whenever:

hi < αh max ,

(10) 2.2.3. Anatomical MRI processing Prior to the stimulation sessions, a cerebral anatomical T1-weighted MRI was acquired at 3T (Achieva 3.0T TX, Philips, Netherlands) for each subject. The MRI data were entered in the TMS neuronavigation system (Localite GmbH, Germany) for post-processing. First, the anatomical hotspot (HSanat ) for the right FDI muscle was defined according to its classical location in the hand knob of the left precentral gyrus (primary motor cortex) (Yousry et al., 1997). Second, a grid of stimulation positions was defined as a square of seven by seven points centered on HSanat and oriented with respect to the central sulcus (Fig. 1a). This grid was automatically generated on the 3D segmented cortex using the dedicated function of the neuronavigation software. The distance between points on the cortical surface was set to 7 mm, resulting in a 17.64 cm2 stimulation area. The X-axis was in the lateral-median direction, its origin being at the most lateral point, while the Y-axis ranged from the most anterior point to the most posterior one. All spatial coordinates reported throughout this manuscript were projected to and expressed in this 2D-Cartesian coordinate system.

with α ∈ (0; 1). However, the situation where hi never reaches this threshold also has to be considered. The occurrence of this bad-case scenario depends on the match between the recorded data and the model hypothesis, as well as on the variability of the data (σi). The second stop criterion is thus based on the first derivative dhi / di . The process stops whenever the model does not get any further information from the data, i.e.:

dhi < β, di

(11)

where β is a real number close to 0. To prevent false-alarms, i.e. temporary decreases of information gain, the criterion has to be satisfied on two consecutive iterations. The default values used in this work are α = 0.3 and β = 0.2 . 2.1.5. Smart planning of Si+1 AutoHS integrates a smart selection of the next stimulation position Si+1. Although any procedure of position selection within the stimulation grid would be possible, a prospective choice of Si+1 allows the model to converge more rapidly to the HS. To do so, the model estimates future entropy resulting from the virtual stimulation of each position:

2.2.4. Manual set-up procedure The manual set-up procedure (mSUP) was performed by an experienced TMS practitioner. Four local experimenters (three neurophysiologists and one psychiatrist) participated in this study. They were randomly distributed and counterbalanced between subjects and sessions, so that they never processed the same subject twice. The experimenters were blind to the other set-up procedures. In order to reduce inter-experimenter variability, normalized guidelines were followed throughout the experiment. First, the HS hunting session had to start on HSanat using an a priori stimulation intensity of 55% of the maximal stimulator output. A minimum of three MEPs had to be recorded on each tested stimulation point. The experimenters were free to increase or decrease the stimulation power during HS hunting. Finally, they were asked to continuously control the coil angles so that the coil was tangential to the scalp surface and the stimulation angle was perpendicular to the central sulcus. Once experimenters considered they found the HS position, they were asked to maintain the coil in its position during rMT and baseline measurements using the neuronavigation marker.

∀ Si +1 |(xi +1, yi +1) ∈ ({1, …, Nx}, {1, …, Ny}), hi +1 = − ∑ ∑ ∑ P (HSμHS σsp |S1: i +1 μ1: i +1) ln (P (HSμHS σsp |S1: i +1 μ1: i +1)), HS μHS σsp

(12) where μi+1 is directly estimated using Eq. (3) and the current knowledge at iteration i about HS , μHS and σsp , and σi+1 is equal to σi. AutoHS finally chooses the cortical position Si+1 inducing the maximal entropy decrease between iterations i and i + 1:

Si +1 | max(hi − hi +1). Si+1

(13)

In order to avoid intensive stimulation over the same cortical positions, each point of the grid was limited to two non-consecutive iterations during the whole hunting procedure. The model selected the point leading to the next best entropy decrease in case of conflict with the previous rule.

2.2.5. Automatized set-up procedure The automatized set-up procedure (aSUP) was completed by the robotized TMS system together with AutoHS (Fig. 2). The procedure began with an initial threshold rMTanat estimated on HSanat . HS hunting then started on HSanat , using an exploration intensity I of 110% of rMTanat . Since the measured rMT increases while moving away from the HS, this intensity always ensured a supraliminal stimulation level while preventing from stimulating too much, in the case where HSanat was too far from the actual HS position. On each stimulated position, five MEP amplitudes (see Section 2.1) were recorded and entered into AutoHS. Upon computation completion, AutoHS provided the next cortical position to be stimulated. This position was then manually selected in the neuronavigation system (see Section 2.2.7), in order to put the robotized coil on it. Once over the HS position defined by AutoHS, the robot maintained the coil on it during rMT and baseline measurements.

2.2. Experimental procedure 2.2.1. Subjects Nineteen right-handed healthy volunteers (14 males, aged 29.6 ± 10.1 years old) participated in the protocol. They all gave their informed and written consent prior to the experiment and received payment for their participation. None of them had either history of psychiatric or neurological disorders, or history of alcohol or substance abuse. They were free of any medicinal treatment likely to modulate their motor cortical excitability levels. MRI and TMS acquisitions were performed at IRMaGe MRI and Neurophysiology facilities (Grenoble, France). The ethical committee of Grenoble University Hospital (ID RCB: 2013-A01734-41) approved this study, which has been registered on ClinicalTrials.gov (number NCT02168413). 310

NeuroImage 153 (2017) 307–318

S. Harquel et al.

• • • • •

the 2D-coordinates of the HS (in x and y) for each session, reported using a doubled spatial resolution grid (3.5 mm between points instead of 7 mm), the HS shift between the first and second session (in mm), defined as the Euclidean distance between the two HS: ∑x, y (HS1 − HS2 )2 , the rMT (in % of maximal stimulator output), the procedure duration, the mean and standard deviation of the measured baseline.

2.3.2. EMG data EMG data were processed using CortExTool, a Matlab toolbox developed in our lab and freely available to the community (Harquel et al., 2016b). Data were first band-pass filtered (50–600 Hz), and trials containing any muscular activity in the pre-stimulus period were rejected. MEPs were then semi-automatically detected, and their peak to peak amplitudes were extracted. Motor maps showing the MEP mean amplitudes on each stimulation point were generated using a spline 2D interpolation between points (seven iterations). 2.3.3. Virtual data Virtual data were generated using the real MEP distributions recorded during motor mapping. For each stimulation position Si, N MEPs were drawn from a Gaussian distribution N (μSi , σS2i ), where μSi and σSi corresponded respectively to the mean and standard deviation of the 10 MEPs recorded on Si during motor mapping. In order to assess the influence of the parametrical form of the real MEP distribution on the HS hunting performance, amplitudes were multiplied by random values drawn from an uniform distribution U (1 − ϵ, 1 + ϵ) where ϵ represented a noise factor. Three noise factors were tested: 0 (no noise), 0.5, and 0.9 (strong noise). High values modified the mean and standard deviation, while significantly increasing the skewness and kurtosis of the distribution. Each MEP amplitude yi was thus calculated as follows:

Fig. 2. Automatized set-up procedure. The supraliminal stimulation intensity I used for HS hunting is first derived from an initial rMT estimated on HSanat . HS hunting is then iteratively performed by AutoHS, in closed loop with the robotized TMS and the EMG recording systems. Once one of the stop criteria has been satisfied, the model settles the TMS coil over the estimated HS where the actual rMT is assessed.

2.2.6. Motor mapping An extensive motor mapping procedure was done on each subject during one of the two stimulation sessions. Motor mapping was randomly done before, between or after the two procedures. It consisted in the recording of ten MEPs on each point of the 7*7 stimulation grid using the robotized system. If the stimulation failed to evoke a consistent MEP (>50 μV ) in three consecutive trials, the coil was moved on to the next site. The stimulation intensity was set to 120% of the lowest rMT found (by either mSUP or aSUP) for the corresponding subject.

y′i ∼ N (μSi , σS2i ) ni ∼ U (1 − ϵ, 1 + ϵ) yi = y′i × ni .

2.2.7. TMS parameters TMS was delivered by a butterfly coil (2*75 mm) MagPro Cool B65RO (MagVenture A/S, Denmark) plugged into a MagPro X100 stimulator system (MagVenture A/S, Denmark). Biphasic pulses were delivered using a randomized inter-stimuli interval (ISI) between three and five seconds during mSUP, aSUP and motor mapping. A tensecond ISI was used during the baseline assessment of mSUP and aSUP. The coil position and orientation with respect to the individual cortical anatomy were controlled using a neuronavigation system (Localite GmbH, Germany). Every target and entry points were defined to be normal to the scalp surface. A robotized system (Axilum Robotics, France) was used during aSUP and motor mapping. All rMTs were assessed automatically using the PEST procedure (Awiszus and Borckardt, 2011) on 30 trials.

2.3.4. Statistics The results were analyzed through Bayesian statistics. Data were processed using JASP software (JASP Team (2016)), as well as modified scripts from Kruschke (2014) using R language (R Core Team, 2016) together with the rjags library. Three basic and standard tests were used: the Bayesian tests equivalent to paired t-tests, repeated measures ANOVAs, and linear regressions. Prior and likelihood functions were modeled using Gaussian or Cauchy distributions (the normality of processed data was systematically checked using Shapiro-Wilk tests). Prior distributions for paired comparison between mSUP and aSUP were kept wide (i.e. uninformative) and centered on 0, as no previous comparison between aSUP and mSUP was ever done before. Posterior distributions were computed using a Markov Chain Monte Carlo procedure (10,000 iterations). The convergence and length of the chains were systematically inspected (following the method detailed in Kruschke (2014)). The posteriors were reported in the text using 95% highest density intervals (HDI95). HDI95 indicates the most probable values for a tested parameter. Evidence of observed effects was provided using Bayesian factors (see Mulder and Wagenmakers, 2016 for an extensive review on this topic): BF01, BF10 (which is the multiplicative inverse of BF01), and BFincl denoted the level of evidence of the null hypothesis, the alternate hypothesis (non-signed difference), and the inclusion of a specific parameter in a model (ANOVA) respectively. Bayes factors were interpreted following the cut-offs proposed by Jeffreys (1998) (values between 1 and 3 denoted absence of evidence to anecdotal evidence, values between 3 and 10 denoted moderate evidence, etc.).

2.2.8. EMG parameters The electromyographic activity of the first dorsal interosseous muscle of the dominant hand (right) was recorded using EMG electrodes placed in a belly-tendon montage. EMG data were amplified, band-pass filtered (50– 6000 Hz), and finally sampled at 12 KHz using a Dantec Keypoint portable EMG recording system (Natus Medical Inc., USA). The EMG activity was recorded in a −200 ms to +600 ms window surrounding stimulation onset. 2.3. Data processing 2.3.1. Data of interest The following data were investigated for each procedure: 311

NeuroImage 153 (2017) 307–318

S. Harquel et al.

the dispersion of HS positions found by global maximum detection procedure and AutoHS, together with the COG, for all subjects using 10 noiseless MEPs per point (N=10 and ϵ = 0 ).

3. Results All subjects went through the entire experiment without any major problem. The average MRI co-registration error, as calculated by the neuronavigation software, was 1.97 ± 0.30 mm . Two subjects were excluded from the experimental comparison, as experimenters did not follow the stimulation angle constraint. Results regarding the reliability and reproducibility of AutoHS on virtual data are first described, followed by the experimental comparison between automatized and manual procedures.

3.2. Experimental comparison between automatized and manual procedures HS shift Fig. 5a shows all the HS positions inferred by both methods on all subjects and sessions, superimposed on corresponding motor maps. The mean HS shifts were 10.1 ± 6.7 mm for mSUP and 8.2 ± 5.4 mm for aSUP (Fig. 5b). The paired comparison of these two distributions failed to draw any conclusion regarding their difference (BF10 = 0.4 , BF01 = 2.4 ). The average distributions of MEP amplitude surrounding HS were computed for each procedure by averaging the normalized motor maps (i.e. by dividing each point of the maps by their maximum) centered on the HS inferred within the same session. The average maps were consequently scaled between 0, no MEP activity found in any subject, and 1, maximum MEP activity found in all subjects (Fig. 5a). aSUP tended to define HS in the center of a 2D Gaussian shape peaking at 0.71, whereas the HS map inferred by mSUP was more blurry and spread out towards the anterior direction, leading to a smaller peak value of 0.55.

3.1. Model validation AutoHS was tested on virtual data using either 3, 5, 7 or 10 recorded MEPs per stimulation position, and 3 noise factors (ϵ = 0, 0.2, 0.9). 100 virtual motor maps were generated for each subject and parameter set, leading to a total of 20,400 maps. Overall, it appeared that the model found the HS position in a fast and reliable way, while being more reproducible and less prone to noise than global maximum detection procedures. Model reliability The model converged to the HS position for all virtual maps and parameter sets after 13.5 ± 3.3 iterations, i.e. stimulation positions (10.7 ± 1.4 for N=5 and ϵ = 0 , using the model's default value on noiseless data). Fig. 3a and 3b shows in detail the proceedings of the model on one representative virtual motor map. As expected, the entropy h rapidly decreased during HS hunting, while the marginal distributions peaked toward their final estimations of HS position, μHS and σsp values. The motor mapping data used in this example is plotted in Fig. 3c. The estimated MEP maximum (μHS = 1, 750 μV ) was close to the recorded maximum (1, 689 μV ), and the spatial spread of the muscle representation was correctly estimated (σsp = 8.4 mm ) by the model.

rMT Fig. 6a presents the results regarding rMTs assessed from mSUP and aSUP, for all subjects and sessions. The Bayesian repeated measure ANOVA indicated that rMTs were moderately similar between procedures (BF01 = 3.31). The mean observed rMT was 46.4 ± 7.5% for mSUP and 45.7 ± 7.1% for aSUP. rMTs were strongly correlated between the two sessions for both methods (Fig. 6c). However, the most probable values of the correlation and its evidence were higher for aSUP (HDI95 (ρ) = [0.73; 0.97], BF10 = 6e 4 ) than for mSUP (HDI95 (ρ) = [0.43; 0.91], BF10 = 185). A Bayesian ANOVA on rMTs assessed during mSUP using experimenters as fixed factor, and subjects and sessions as noise factors, provided moderate evidence that the experimenters did not influence rMT (BF01 = 3.45).

Spatial reproducibility A comparison between the spatial reproducibility obtained with AutoHS and with the classical HS hunting procedure was performed by computing the spatial spread of HS positioning on 100 virtual motor maps for each subject (Fig. 4). The classical HS hunting procedure is usually based on the observed MEP amplitude, and defines the HS as the cortical position where maximal muscle contraction is induced. To simulate such a procedure, virtual data were analyzed using a global maximum detection procedure, which settled the HS on the point showing the absolute maximum MEP mean amplitude over the complete motor maps (49 points). The spatial spread was defined as the mean standard deviation of HS positions in both axes. The center of gravity (COG) of the motor maps was also computed using x = ∑ (xμ)/ ∑ μ , and y = ∑ ( yμ)/ ∑ μ as coordinates, with μ being the mean MEP amplitude for each point of the grid. Overall, the HS spatial spread was smaller using AutoHS (2.8 ± 0.8 mm ) than the global maximum detection procedure (3.9 ± 1.6 mm ), as stated by the main effect of the method (BFincl =+ inf , HDI95 = [1.0; 1.4] mm). As expected, the number of recorded MEPs by stimulation points as well as the noise factor modulated the reproducibility of both methods. HS hunting procedures performed better using more data (main effect of N, BFincl =+ inf ) and noise-free data (main effect of ϵ, BFincl = 4e+11). It appeared that the global maximum search procedure was more sensitive to the noise parameter than AutoHS (interaction effect, BFincl = 1e+7), whereas both methods seemed to be equally influenced by the number of MEPs (BF01 = 5.56 ). Fig. 4a shows two representative examples of this finding. On the top row, the size of the dispersion ellipses was modulated by the number of recorded MEPs for both procedures. On the bottom row, the HSs defined as global MEP maxima tended to spread away from their original position (i.e. without noise) as the noise factor increased, whereas those inferred by AutoHS remained the same. Fig. 4c shows

Estimated duration Fig. 6a shows the comparison between the duration of mSUP and aSUP, for all subjects and sessions. The theoretical duration of aSUP was simulated by considering AutoHS as fully implemented in the TMS acquisition loop (Fig. 3), using an ISI of 5 seconds and an interval of 3 seconds between stimulation positions, compensating for the robotized arm movement. Both durations included the rMT assessment time and thus represented the time of the whole set-up procedures, from the first TMS pulse to the final rMT outcome. It appeared that aSUP lasted shorter than mSUP (HDI95 = [1.6; 4.6] min, BF10 = 273). On average, mSUP lasted 11.5 ± 3.7 min while aSUP lasted 8.3 ± 1.7 min. A Bayesian ANOVA was run on the mSUP duration using experimenters as fixed factor, and subjects and sessions as random factors. No conclusion could be drawn regarding the influence of the experimenters on the procedure time (BFincl = 1.16, BF01 = 0.87). Baseline measurements The baseline means were 999 ± 848μV for mSUP and 852 ± 609μV for aSUP. There was a moderate evidence that they did not differ between methods (BF01 = 3.6 for the method effect). The average baseline standard deviations were 526 ± 363μV for mSUP and 462 ± 369 μV for aSUP, and did not differ between methods (BF01 = 4.6 ). 4. Discussion This work provides the first proof of concept for a practical and fully automatized set-up procedure dedicated to TMS protocols. The core element of this procedure is a Bayesian model, AutoHS, which guides 312

NeuroImage 153 (2017) 307–318

S. Harquel et al.

Fig. 3. Inferred probability distributions throughout AutoHS, for one virtual motor map of a representative subject (S16). a: (top) Evolution of the marginal distribution Pm (HS ) (as defined by Eq. (8)) over time. Iteration 0 represents the prior distribution P (HS ). (bottom) Entropy as a function of iterations i. AutoHS settled on the HS position after the 7th stimulation position, inducing an entropy under the defined threshold αhmax (see text). b: From left to right, evolution over time of the marginal distributions Pm (σsp ) and Pm (μHS ) respectively. For this example, AutoHS estimated μHS = 1, 750 μV and σsp = 8.4 mm . c: Final motor mapping. The black area corresponds to the spatial distribution (1 σ ellipse) of the HS positions found during the 100 simulations (N=5, ϵ = 0 ).

2007; Desmurget et al., 2014; Giszter, 2015). The conception of M1 as a simple body plan is shifting to a more complex picture that includes action representations (Graziano, 2016). Our motor mapping data (see S1, S4, S5, S9 and S10), as well as many examples in the literature, suggest that the FDI representation is far from being uni-modal in some subjects (see the motor or MT maps shown in Gugino et al. (2001), Neggers et al. (2004), Sparing et al. (2008), Meesen et al. (2011), van de Ruit et al. (2015), Meincke et al. (2016)). Additionally, significant muscular contractions can also be induced when stimulating pre-motor areas (Ahdab et al., 2016), which are sometimes included in the most anterior part of the stimulation space, depending on the gyri configuration (see S3, S6, S11 and S15). Despite the fact that the spatial spread of the electrical field induced by TMS might be limiting for revealing such complex multiple representations in most of subjects, choosing multi-modal distributions for muscle representations may optimize the method for the rest of them. Testing our procedure with more focal coils could also allow to better observe and take into account this phenomenon. The assumption of Gaussianity regarding MEP distributions has been extensively tested on virtual data. The model appeared to be very robust to its violation, i.e. when MEP distributions were not Gaussian. As shown in Fig. 4a and b, the reproducibility of the method stayed unaffected by increasing noise factors, contrary to the global maximum search procedure. As a trade-off, the convergence speed was however reduced from 11 (perfect Gaussian distributions) to 17 iterations (noisy distributions) on average. The prior regarding the HS position had little influence on model performance. By choosing a large value for σprior (14 mm), the prior was kept uninformative so that its influence was

the HS hunting. Results of both simulations and experimental comparison against classical methodologies are promising. The automatized procedure appears to be at least as reliable as the manual one, while being more reproducible and significantly shorter than the necessary processing time of this step. Moreover, this work challenges the classical hotspot definition, and supports the need for future developments, optimizations, and extensions of the model to other experimental fields in TMS. 4.1. Robustness of the model The hypothesis regarding the relationship between recorded MEP amplitudes and distance to the HS may be one of the most influential factor on model performance. It is based on the assumption that muscle representations are centered symmetrically around the HS. Even though this hypothesis appeared to be violated with varying degrees in most of the subjects, it did not prevent the model to systematically converge towards reliable and reproducible solutions. However, assuming that targeted muscles are restricted to uni-modal representations may be erroneous for any kind of HS hunting procedure. The functional organization of the primary motor cortex is more complex than a successive arrangement of muscle representations along the central sulcus (Schieber, 2001). While S1 somatotopy has a discrete and segregated layout of body part representations, M1 in contrast displays a more integrated and overlapping functional organization. Recent invasive stimulation studies on both animals and humans clearly identified representations of several synergies and multijoint movements over the motor cortex (Graziano and Aflalo, 313

NeuroImage 153 (2017) 307–318

S. Harquel et al.

Fig. 4. Reproducibility of HS positioning on virtual motor maps, using AutoHS or maximal mean MEP amplitude (Max). a: Results obtained on two representative subjects (S10 and S13). The black and white areas represent the spatial distribution (1 σ ellipse) of the HS found by AutoHS and Max procedures for 100 simulations respectively. The distributions are superimposed on motor maps showing the mean MEP amplitude on each point. The top and bottom rows show the influence of N (number of recorded MEPs per stimulation site) and ϵ (noise factor) parameters on the generated motor maps respectively. b: Group result showing the mean HS spatial spread (size of the spatial distribution, see text) observed using AutoHS (black bars) and Max (white bars) procedures for 100 simulations. Data are grouped in X-axis by N and ϵ parameters. c: Spatial distribution of the HS found by AutoHS (black) and Max (white) for all subjects, using 10 noise-free MEPs per point. The center of gravity (COG) is represented by a white cross.

4.2. Redefining the HS

minimized compared to the iterative exploration process. Recorded data supported this choice, as only one subject (S14) showed maximal contractions exactly on the anatomical prior for HS position. For these particular cases, choosing a lower value for σprior could then increase the convergence speed (about 3 iterations for σprior = 4 mm ). On average, choosing a wider or a narrower prior distribution would slightly decrease or increase the convergence speed respectively, depending on the distance between the real and anatomical HS, i.e. depending on the goodness of the prior. Variable domains should always cover the most probable physiological values regarding MEP amplitudes, or spatial spreads of muscle representations. Choosing domains larger than the proposed defaults would enable the model to cover more particular cases. However, recording maximum MEP amplitudes higher than 4 μV , or spatial muscle representations larger than 14 mm, are highly improbable using the proposed stimulation intensity. Choosing more precise domains, i.e. using smaller steps between values, would theoretically make the probabilistic inference more precise, while consequently increasing computing time significantly. Finally, the values defining the stop criteria are always representative of the trade-off between procedure precision and duration. Higher values for the entropy thresholds (α and β) would allow decreasing procedure time, as they would be reached faster, while lowering the reproducibility of the calculation (see the example depicted in Fig. 3a).

By developing AutoHS, we sightly revisited the interpretation of the most commonly used HS definition. Theoretically, the HS is the site where TMS induces maximal contractions over a targeted muscle. Most TMS experimenters generally assess HS by finding the position of the global maximum within the stimulation space, by computing the mean of several trials for each site. In contrast, AutoHS looks for the probabilistic position of this maximum. By fully taking inter-trial variability into account, probabilistic methods are more reliable than global maximum search procedures and are less prone to outliers. This is especially important given the high trial-to-trial variability of MEP amplitude and our results on virtual data strongly support the beneficial contribution of the probabilistic approach to HS determination. We found that the average HS shift was significantly reduced between virtual sessions with AutoHS, especially when data were noisy. Interestingly, our procedure led to a good solution between the global maximum of the motor maps and its center of gravity. Like the COG, it tended to position the HS in the center of large muscle representations (see for example Fig. 4c, subjects S3, S4, S9, S11), while it rather converged to global maximum position for narrower motor maps (see S16 and S17). The optimization in HS positioning had a direct effect regarding the reproducibility of outcomes during the experimental comparison with manual methods. The automatized set-up procedure showed improved inter-session reproducibility regarding the measured rMT. Experimental

314

NeuroImage 153 (2017) 307–318

S. Harquel et al.

Fig. 5. Experimental comparison of mSUP and aSUP regarding the reproducibility of HS positioning. a: HS positions assessed by mSUP (blue) and aSUP (orange) for all subjects and sessions. The positions are superimposed on motor maps showing the mean MEP amplitude on each point. The black stars indicate positions inferred within the same session as motor mapping (MM session). The two supplementary maps (bottom right) represent the average distributions of MEP amplitude surrounding HS inferred by aSUP (left) and mSUP (right). b: HS shift (Euclidean distance between HS found in the two sessions) as a function of procedure. mSUP and aSUP distributions are colored in blue and orange respectively. Pair-wise comparisons are highlighted using segments.

Fig. 6. Experimental comparison of mSUP and aSUP performances based on rMTs and procedure durations. mSUP and aSUP distributions are colored in blue and orange respectively. Paired-wise comparisons are highlighted using segments. a: distributions of inferred rMT (top) and procedure durations (bottom) using mSUP and aSUP, for all subjects and sessions. b: distributions of mSUP values for each of the four experimenters compared to the aSUP values (top: rMT, bottom: procedure duration). c: inter session correlation between rMTs for aSUP (top) and mSUP (bottom), where a sample of credible regression lines are plotted.

315

NeuroImage 153 (2017) 307–318

S. Harquel et al.

have a strong influence over the inferred rMT on some subjects (Richter et al., 2013). Using a single stimulation angle (perpendicular to the precentral gyrus) might be limiting here, even if it is theoretically supposed to be the most adapted. Two factors might at least influence the optimal stimulation angle within an actual experiment: errors in the MRI registration and complex modulations of the induced electrical field geometry by neuro-anatomical properties (Thielscher et al., 2011; Janssen et al., 2015). In this study, two subjects were excluded of the experimental comparison as experimenters did not follow the angle constraint (disparities of 30 and 50°), inducing a rMT difference of ± 10% between the two procedures (opposite in direction). It appeared then that the pre-defined stimulation angle was not optimal in at least one case. In future developments of the method, the stimulation angle should become a variable of the model and be estimated during HS hunting. Finally, the core principle of the model can be applied beyond the specific HS hunting problem, i.e. in any experimental procedure aiming at finding a cortical target optimized using TMS readouts. More specifically, it can be adapted to find the probabilistic origin of any a priori known distribution of TMS-induced readouts over a stimulation space. In principle, the TMS readouts can be either behavioral or neurophysiological. For example, the model could be modified for developing a phosphene tresholding procedure over the occipital cortex (Dugué et al., 2011; Herring et al., 2015). It could also be turned into a functional area hunting procedure, by estimating the stimulated cortical target best able to modulate the performance of a specific cognitive task. Finally, the model could also be applied in concurrence with other brain imaging technique, like electroencephalography (EEG). Under this condition, it would allow to find stimulation areas optimizing the emergence of any EEG components, such as evoked potentials or induced oscillations, in order to help the characterization of specific outcomes or biomarkers in various clinical or fundamental issues (see Farzan et al., 2016 for a recent review). Such procedure would also be especially useful in future closed-loop approaches combining TMS and electrophysiology (Bergmann et al., 2016), that aim at adapting the stimulation parameters based on readouts in real time.

results also confirmed that it was more reliable than the manual procedure in finding a position inducing maximal contractions. The average motor maps presented in the bottom right part of Fig. 5a are especially meaningful. Overall, the automatized method sets the HS over the maximum contraction point at a higher rate than the manual one. However, our work failed to bring strong evidence for an improvement of the measured outcomes reliability, i.e. smaller rMT and baseline variability, due to a high inter-subject variability between procedure performance. Three non-exclusive hypotheses can be posited to explain the lack of result regarding reliability. First, as suggested by the small BF values, further data have to be acquired in the future in order to accumulate more evidence in favor of one or another hypothesis. Second, the current neuronavigation system precision, about 3 mm when adding the registration error and the camera precision (1 mm), and TMS spatial resolution might not have the ability to detect fine differences in HS shift induced by AutoHS compared to classical approaches. Using smaller coils and more precise camera might be two solutions for solving these issues respectively. Finally, our model is still not fully optimized regarding its hypotheses and stop criteria. Developing better and more adaptable priors about MEP spatial distribution, as well as lowering the thresholds for stop criteria could improve the reliability of HS positioning and measured outcomes, while slightly increasing the number of needed stimulation targets, as a trade-off (see Section 4.1). 4.3. Future developments This initial study paves the way for future optimizations and developments of AutoHS, by reconsidering some of the hypotheses of the model, including angle estimation in its process, and extending its principle to other experimental fields. First, the subject anatomy and gyral geometry should be better taken into account in the stimulation space. The model should first include the scalp-to-cortex distance and adapt the stimulation intensity to it (Stokes et al., 2005), as it may vary within the stimulation space. Second, even if the square grid used in this work covered most of the precentral gyrus, the exploration space would be optimized using a non-cartesian stimulation grid exactly following the gyrus shape towards the X-axis. A recent study showed that TMS could better discriminate between muscles when following the gyral geometry compared to following its main direction (Raffin and Pellegrino, 2015). Second, the parametrical form of the distribution of MEP amplitudes knowing the positions of HS and the stimulation point (P (μi |A)) is simplistic and non-optimal for several subjects, especially for those presenting a bi-modal distribution (see above). Less constrained parametrical forms could be tested. For instance, introducing covariance terms in the 2D Gaussian would allow capturing decrease of MEP amplitudes not geometrically aligned with the grid; using Gaussian mixtures would allow capturing multi-modality; etc. Second, the model should estimate the optimized stimulation angle within its process. As shown in previous studies, the stimulation angle may

5. Authors contribution S.H., J.D., E.R., and A.C. designed research. S.H., E.R., B.P., G.I., C.M. and A.C. performed research. S.H. and B.P. analysed data. S.H., J.D., E.R., O.D. and A.C. wrote the paper. Acknowledgments This work was funded by the Pôle Grenoble Cognition grant “PGC2014” and the Agence Nationale pour la Recherche grant “ANR-15CE37-0015-02”. Data were acquired on a platform of France Life Imaging Network partly funded by the grant “ANR-11-INBS-0006”. The authors declare no conflict of interest.

Appendix The derivation leading to Eq. (7) used during Bayesian inference is:

316

NeuroImage 153 (2017) 307–318

S. Harquel et al.

P (HSμHS σsp |S1: T μ1: T ) P (HSμHS σsp S1: T μ1: T ) = P (S1: T μ1: T ) 1 = P (HSμHS σsp S1: T μ1: T ) Z1 = = =

1 P (HS ) P (μHS ) P (σsp ) Z1 T

[Bayes theorem] [P (S1: T μ1: T ) = Z1 (constant)] T

T

∏ P (Si ) ∏ P (μi |HSμHS σsp Si ) i =1 T

1 P (HS ) Z1 Z 2

∏ P (Si ) ∏ P (μi |HSμHS σsp Si )

1 P (HS ) Z1 Z 2

∏ P (μi |HSμHS σsp Si )

i =1 T

[by definition of the joint]

i =1

[P (μHS ) P (σsp ) = Z 2 (uniforms)]

i =1

[P (Si ) = 1(Dirac)]

i =1

T

∝ P (HS )

∏ P (μi |HSμHS σsp Si )

□

i =1 T

∝ P (HS )

∏ P (μi |A)

(14)

i =1

recalling that A is the conjunction of variables HS , μHS , σsp and Si. Gugino, L.D., Rafael Romero, J., Aglio, L., Titone, D., Ramirez, M., Pascual-Leone, A., Grimson, E., Weisenfeld, N., Kikinis, R., Shenton, M.-E., 2001. Transcranial magnetic stimulation coregistered with MRI a comparison of a guided versus blind stimulation technique and its effect on evoked compound muscle action potentials. Clin. Neurophysiol. 112 (October (10)), 1781–1792. Hallett, M., 2000. Transcranial magnetic stimulation and the human brain. Nature 406 (6792), 147–150. Harquel, S., Bacle, T., Beynel, L., Marendaz, C., Chauvin, A., David, O., 2016a. Mapping dynamical properties of cortical microcircuits using robotized TMS and EEG Towards functional cytoarchitectonics. NeuroImage 135 (May), 115–124. Harquel, S., Beynel, L., Guyader, N., Marendaz, C., David, O., Chauvin, A., October 2016b. CortExTool: a toolbox for processing motor cortical excitability measurements by transcranial magnetic stimulation. URL 〈https://hal.archivesouvertes.fr/hal-01390016/document〉. Herbsman, T., Forster, L., Molnar, C., Dougherty, R., Christie, D., Koola, J., Ramsey, D., Morgan, P.S., Bohning, D.E., George, M.S., Nahas, Z., 2009. Motor Threshold in Transcranial Magnetic Stimulation The Impact of White Matter Fiber Orientation and Skull-to-Cortex Distance. Human. brain Mapp. 30 (July (7)), 2044–2055. Herring, J.D., Thut, G., Jensen, O., Bergmann, T.O., 2015. Attention Modulates TMSLocked Alpha Oscillations in the Visual Cortex. J. Neurosci. 35 (October (43)), 14435–14447. Herwig, U., Schönfeldt-Lecuona, C., Wunderlich, A.P., von Tiesenhausen, C., Thielscher, A., Walter, H., Spitzer, M., 2001. The navigation of transcranial magnetic stimulation. Psychiatry Res.: Neuroimaging 108 (November (2)), 123–131. Janssen, A.M., Oostendorp, T.F., Stegeman, D.F., 2015. The coil orientation dependency of the electric field induced by TMS for M1 and other brain areas. J. Neuroeng. Rehabil. (May), 12. Jeffreys, H., August 1998. The Theory of Probability. OUP Oxford. Julkunen, P., Säisänen, L., Danner, N., Niskanen, E., Hukkanen, T., Mervaala, E., Könönen, M., 2009. Comparison of navigated and non-navigated transcranial magnetic stimulation for motor cortex mapping, motor threshold and motor evoked potentials. NeuroImage 44 (February (3)), 790–795. Jung, N.H., Delvendahl, I., Kuhnke, N.G., Hauschke, D., Stolle, S., Mall, V., 2010. Navigated transcranial magnetic stimulation does not decrease the variability of motor-evoked potentials. Brain Stimul. 3 (April (2)), 87–94. Kantelhardt, S.R., Fadini, T., Finke, M., Kallenberg, K., Siemerkus, J., Bockermann, V., Matthaeus, L., Paulus, W., Schweikard, A., Rohde, V., Giese, A., 2009. Robot-assisted image-guided transcranial magnetic stimulation for somatotopic mapping of the motor cortex a clinical pilot study. Acta Neurochir. 152 (November (2)), 333–343. Kruschke, J., November 2014. Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan. Academic Press. Lefaucheur, J.-P., André-Obadia, N., Antal, A., Ayache, S.S., Baeken, C., Benninger, D.H., Cantello, R.M., Cincotta, M., de Carvalho, M., De Ridder, D., Devanne, H., Di Lazzaro, V., Filipović, S.R., Hummel, F.C., Jääskeläinen, S.K., Kimiskidis, V.K., Koch, G., Langguth, B., Nyffeler, T., Oliviero, A., Padberg, F., Poulet, E., Rossi, S., Rossini, P.M., Rothwell, J.C., Schönfeldt-Lecuona, C., Siebner, H.R., Slotema, C.W., Stagg, C.J., Valls-Sole, J., Ziemann, U., Paulus, W., Garcia-Larrea, L., 2014. Evidence-based guidelines on the therapeutic use of repetitive transcranial magnetic stimulation (rTMS). Clin. Neurophysiol. 125 (November (11)), 2150–2206. Lefaucheur, J.-P., Picht, T., 2016. The value of preoperative functional cortical mapping using navigated TMS. Neurophysiol. Clin./Clin. Neurophysiol. 46 (April (2)), 125–133. Meesen, R.L., Cuypers, K., Rothwell, J.C., Swinnen, S.P., Levin, O., 2011. The effect of long-term TENS on persistent neuroplastic changes in the human cerebral cortex. Human. Brain Mapp. 32 (June (6)), 872–882. Meincke, J., Hewitt, M., Batsikadze, G., Liebetanz, D., 2016. Automated TMS hotspothunting using a closed loop threshold-based algorithm. NeuroImage 124 (January), 509–517, (Part A). Mulder, J., Wagenmakers, E.-J., 2016. Editors introduction to the special issue Bayes

References Ahdab, R., Ayache, S.S., Brugières, P., Farhat, W.H., Lefaucheur, J.-P., 2016. The Hand Motor Hotspot is not Always Located in the Hand Knob A Neuronavigated Transcranial Magnetic Stimulation Study. Brain Topogr. (March), 1–8. Awiszus, F., 2003. TMS and threshold hunting. Suppl. Clin. Neurophysiol. 56, 13–23. Awiszus, F., Borckardt, J., 2011. TMS motor threshold assessment tool (MTAT 2.0). Baek, J., Lesmes, L.A., Lu, Z.-L., 2016. qPR An adaptive partial-report procedure based on Bayesian inference. J. Vision. 16 (August (10)), 25. Barker, A.T., Jalinous, R., Freeston, I.L., 1985. Non-invasive magnetic stimulation of human motor cortex. Lancet 325 (May (8437)), 1106–1107. Bergmann, T.O., Karabanov, A., Hartwigsen, G., Thielscher, A., Siebner, H.R., 2016. Combining non-invasive transcranial brain stimulation with neuroimaging and electrophysiology Current approaches and future perspectives. NeuroImage, (URL 〈http://www.sciencedirect.com/science/article/pii/S1053811916001191〉). Bessière, P., Mazer, E., Ahuactzin, J.M., Mekhnacha, K., 2013. Bayesian Programming. CRC Press, Boca Raton, Florida. Bessière, P., Lebeltel, O., 2008. Basic Concepts of Bayesian Programming. In: Bessière, P., Laugier, C., Siegwart, R. (Eds.), Probabilistic Reasoning and Decision Making in Sensory-Motor Systems. No. 46 in Springer Tracts in Advanced Robotics. Springer Berlin Heidelberg, pp. 19–48, http://dx.doi.org/10.1007/978-3-540-79007-5_2. Bestmann, S., Feredoes, E., 2013. Combined neurostimulation and neuroimaging in cognitive neuroscience: past, present, and future Combined neurostimulation and neuroimaging. Ann. New Y. Acad. Sci. 1296 (August (1)), 11–30. Bortoletto, M., Veniero, D., Thut, G., Miniussi, C., 2015. The contribution of TMS EEG coregistration in the exploration of the human cortical connectome. Neurosci. Biobehav. Rev. 49 (February), 114–124. Cincotta, M., Giovannelli, F., Borgheresi, A., Balestrieri, F., Toscani, L., Zaccara, G., Carducci, F., Viggiano, M.P., Rossi, S., 2010. Optically tracked neuronavigation increases the stability of hand-held focal coil positioning evidence from transcranial magnetic stimulation-induced electrical field measurements. Brain Stimul. 3 (April (2)), 119–123. Desmurget, M., Richard, N., Harquel, S., Baraduc, P., Szathmari, A., Mottolese, C., Sirigu, A., 2014. Neural representations of ethologically relevant hand/mouth synergies in the human precentral gyrus. Proc. Natl. Acad. Sci. 111 (April (15)), 5718–5722. Dugué, L., Marque, P., VanRullen, R., 2011. The phase of ongoing oscillations mediates the causal relation between brain excitation and visual perception. J. Neurosci. 31 (August (33)), 11889–11893. Farzan, F., Vernet, M., Shafi, M.M.D., Rotenberg, A., Daskalakis, Z.J., Pascual-Leone, A., September 2016. Characterizing and Modulating Brain Circuitry through Transcranial Magnetic Stimulation Combined with Electroencephalography. Frontiers in Neural Circuits 10.URl 〈http://journal.frontiersin.org/Article/10.3389/ fncir.2016.00073/abstract〉. Finke, M., Fadini, T., Kantelhardt, S., Giese, A., Matthaus, L., Schweikard, A., 2008. Brain-mapping using robotized TMS. In: Engineering in Medicine and Biology Society, 2008. EMBS 2008. Proceedings of the 30th Annual International Conference of the IEEE, pp. 3929–3932. Ginhoux, R., Renaud, P., Zorn, L., Goffin, L., Bayle, B., Foucher, J., Lamy, J., Armspach, J.P., de Mathelin, M., 2013. A custom robot for transcranial magnetic stimulation: first assessment on healthy subjects. In: 2013 Proceedings of the 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). IEEE, pp. 5352–5355. Giszter, S.F., 2015. Motor primitives new data and future questions. Curr. Opin. Neurobiol. 33 (August), 156–165. Graziano, M.S.A., 2016. Ethological action maps a paradigm shift for the motor cortex. Trends Cogn. Sci. 20 (February (2)), 121–132. Graziano, M.S.A., Aflalo, T.N., 2007. Mapping Behavioral Repertoire onto the Cortex. Neuron 56 (October (2)), 239–251.

317

NeuroImage 153 (2017) 307–318

S. Harquel et al.

Schieber, M.H., 2001. Constraints on Somatotopic Organization in the Primary Motor Cortex. J. Neurophysiol. 86 (November (5)), 2125–2143. Sollmann, N., Hauck, T., Obermüller, T., Hapfelmeier, A., Meyer, B., Ringel, F., Krieg, S.M., 2013. Inter- and intraobserver variability in motor mapping of the hotspot for the abductor policis brevis muscle. BMC Neurosci. 14 (1), 94. Sparing, R., Buelte, D., Meister, I.G., Pauš, T., Fink, G.R., 2008. Transcranial magnetic stimulation and the challenge of coil placement A comparison of conventional and stereotaxic neuronavigational strategies. Human. Brain Mapp. 29 (January (1)), 82–96. Stokes, M.G., Chambers, C.D., Gould, I.C., Henderson, T.R., Janko, N.E., Allen, N.B., Mattingley, J.B., 2005. Simple Metric For Scaling Motor Threshold Based on ScalpCortex Distance Application to Studies Using Transcranial Magnetic Stimulation. J. Neurophysiol. 94 (December (6)), 4520–4527. Thielscher, A., Opitz, A., Windhoff, M., 2011. Impact of the gyral geometry on the electric field induced by transcranial magnetic stimulation. NeuroImage 54 (January (1)), 234–243. van de Ruit, M., Perenboom, M.J.L., Grey, M.J., 2015. TMS Brain Mapping in Less Than Two Minutes. Brain Stimul. 8, 231–239. Wassermann, E., Epstein, C. (Eds.), November 2012. The Oxford handbook of transcranial stimulation. Oxford University Press, Great Clarendon Street, Oxford, OX2 6DP. Wassermann, E.M., 2002. Variation in the response to transcranial magnetic brain stimulation in the general population. Clin. Neurophysiol. 113 (July (7)), 1165–1171. Weiss, C., Nettekoven, C., Rehme, A.K., Neuschmelting, V., Eisenbeis, A., Goldbrunner, R., Grefkes, C., 2013. Mapping the hand, foot and face representations in the primary motor cortex Retest reliability of neuronavigated TMS versus functional MRI. NeuroImage 66 (February), 531–542. Yousry, T.A., Schmid, U.D., Alkadhi, H., Schmidt, D., Peraud, A., Buettner, A., Winkler, P., 1997. Localization of the motor hand area to a knob on the precentral gyrus. A new landmark. Brain 120 (1), 141–157.

factors for testing hypotheses in psychological research Practical relevance and new developments. J. Math. Psychol. 72 (June), 1–5. Neggers, S.F.W., Langerak, T.R., Schutter, D.J.L.G., Mandl, R.C.W., Ramsey, N.F., Lemmens, P.J.J., Postma, A., 2004. A stereotactic method for image-guided transcranial magnetic stimulation validated with fMRI and motor-evoked potentials. NeuroImage 21 (April (4)), 1805–1817. R Core Team, 2016. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. URL 〈https://www.R-project. org〉. Raffin, E., Pellegrino, G., Di Lazzaro, V., Thielscher, A., Siebner, H.R., 2015. Bringing transcranial mapping into shape Sulcus-aligned mapping captures motor somatotopy in human primary motor hand area. NeuroImage 120 (October), 164–175. Ragazzoni, A., Pirulli, C., Veniero, D., Feurra, M., Cincotta, M., Giovannelli, F., Chiaramonti, R., Lino, M., Rossi, S., Miniussi, C., 2013. Vegetative versus minimally conscious states a study using TMS-EEG, sensory and event-related potentials. PLoS One 8 (2), e57069. Richter, L., Neumann, G., Oung, S., Schweikard, A., Trillenberg, P., 2013. Optimal Coil Orientation for Transcranial Magnetic Stimulation. PLoS One 8 (April (4)), e60358. Rogasch, N.C., Fitzgerald, P.B., 2013. Assessing cortical network properties using TMSEEG. Human. Brain Mapp. 34 (July (7)), 1652–1669. Rossi, S., Hallett, M., Rossini, P.M., Pascual-Leone, A., 2009. Safety, ethical considerations, and application guidelines for the use of transcranial magnetic stimulation in clinical practice and research. Clin. Neurophysiol. 120 (December (12)), 2008–2039. Rossini, P.M., Barker, A.T., Berardelli, A., Caramia, M.D., Caruso, G., Cracco, R.Q., Dimitrijević, M.R., Hallett, M., Katayama, Y., Lücking, C.H., August 1994. Noninvasive electrical and magnetic stimulation of the brain, spinal cord and roots: basic principles and procedures for routine clinical application. Report of an IFCN committee. Electroencephalography and Clinical Neurophysiology, 91 (2), pp. 79–92.

318