# Automatized set-up procedure for transcranial magnetic stimulation protocols

Sylvain Harquel, Julien Diard, Estelle Raffin, Passera Brice, Dall'Ignà Gaelle, Christian Marendaz, Olivier David, Alan Chauvin

► To cite this version:

Sylvain Harquel, Julien Diard, Estelle Raffin, Passera Brice, Dall'Ignà Gaelle, et al. Automatized set-up procedure for transcranial magnetic stimulation protocols. NeuroImage, 2017, Epub ahead of print. 10.1016/j.neuroimage.2017.04.001. inserm-01504574

HAL Id: inserm-01504574

https://inserm.hal.science/inserm-01504574

Submitted on 10 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Automatized set-up procedure for transcranial magnetic stimulation protocols

S. Harquel<sup>a,b,c,d,*</sup>, J. Diard<sup>a,b</sup>, E. Raffin<sup>a,c</sup>, B. Passera<sup>a,b</sup>, G. Dall'Ignà<sup>e,c</sup>, C. Marendaz<sup>a,b</sup>, O. David<sup>a,c</sup>, A. Chauvin<sup>a,b</sup>

<sup>a</sup>Univ. Grenoble Alpes, F-38000 Grenoble, France

<sup>b</sup>CNRS, UMR 5105, Laboratoire de Psychologie et de Neurocognition, LPNC, F-38000 Grenoble, France

<sup>c</sup>INSERM, U1216, Grenoble Institut des Neurosciences, GIN, F-38000 Grenoble, France

<sup>d</sup>CNRS, INSERM, UMS 3552, IRMaGe, F-38000 Grenoble, France

<sup>e</sup>Pôle de Psychiatrie et Neurologie, Centre Hospitalier Universitaire Grenoble Alpes, F-38000 Grenoble, France

Abstract

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.

**Keywords:** robotized TMS, hotspot hunting, cortical excitability, MEP, motor mapping

Highlights

* Automatized set-up procedures would facilitate TMS experiments and increase reproducibility

* We developed a Bayesian model aiming at automatically finding the motor hotspot

* Implementation of this model in a robotized TMS system allows the automatic search for the motor hotspot

* Definition of the motor hotspot is significantly improved in terms of speed and reproducibility

*Corresponding author. Laboratoire de Psychologie et de Neurocognition, CNRS UMR 5105. Université Grenoble Alpes, BSHM. BP47, 38040 Grenoble Cedex 9, France. Tel: +334 76 82 58 91. Mail: sylvain.harquel@univ-grenoble-alpes.fr. Webpage: http://lpnc.univ-grenoble-alpes.fr/Sylvain-Harquel?lang=fr
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 inter-subject 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).

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. Materials & Methods

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 the HS hunting procedure is the first interosseous muscle (FDI).

### 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.

HS 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 (Figure 1a). HS is thus a two dimensional variable expressed in a Cartesian coordinate system. Its $x$ and $y$ values are continuous over the stimulation grid area, so that $(x_{HS}, y_{HS}) \in ([1; N_x], [1; N_y])$, with $N_x$ and $N_y$ being the number of discrete points in the stimulation grid. This work uses a square grid, so that $N_x = N_y = 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:

$$P(\text{HS}) \propto \exp\left(-\left(\frac{(x - x_{\text{HK}})^2}{2\sigma_{\text{prior}}^2} + \frac{(y - y_{\text{HK}})^2}{2\sigma_{\text{prior}}^2}\right)\right), \quad (1)$$

where $x_{HK}$ and $y_{HK}$ are the coordinates of the hand knob, and $\sigma_{prior}$ is the spatial spread of the distribution that quantifies the prior precision (Figure 1b). We use here $\sigma_{prior} = 2$ (in grid units), leading to a spatial spread of 1.4 cm.

$S_i$ is the position of the coil at the $i^{th}$ stimulation, $i$ going from 1 to the total number of stimulated positions $T$. $S_i$ 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 $(x_i, y_i) \in (\{1, 2, \dots, N_x\}, \{1, 2, \dots, N_y\})$.

$\mu_i$ is the mean amplitude observed from $N$ MEPs on $S_i$. $\mu_i$ is a discrete variable over the domain of possible MEP amplitudes $D_{MEP}$ (in $\mu$V). Here, we set $D_{MEP} = \{10, 20, \dots, 4000\}$ $\mu$V.
Here, $\mu_i$ is meant to follow a Gaussian distribution such that

$$P(\mu_i) \propto \frac{1}{\sigma_i \sqrt{2\pi}} \exp\left(-\frac{(\mu - \mu_i)^2}{\sigma_i^2}\right), \qquad (2)$$

where $\sigma_i$ is the standard deviation measured over the N recorded MEP amplitudes (Figure 1c). This work uses 5 MEP recordings per stimulation point, so that $N = 5$. If no MEP was induced after 3 consecutive stimulations, $\mu_i$ was set to 0 and the model moved to the next iteration.

$\mu_{HS}$ represents the expected mean amplitude measured on the HS. Here, in accordance with the HS definition (see 1), AutoHS makes the assumption that $\mu_{HS}$ represents the maximal mean MEP amplitude that can be measured over the whole simulation space. $\mu_{HS}$ is a discrete variable over the domain of possible MEP max amplitudes $D_{max}$ (in $\mu V$). Here, we set $D_{max} = \{250, 500, \dots, 4000\} \mu V$.

Moreover, the model assumes that MEP amplitudes $\mu_i$ decrease as a function of the distance between the stimulated position $S_i$ and the HS position HS (Figure 1d):

$$\mu_i = \mu_{HS} \exp\left(-\left(\frac{(x_i - x_{HS})^2}{2\sigma_{sp}^2} + \frac{(y_i - y_{HS})^2}{2\sigma_{sp}^2}\right)\right). \qquad (3)$$

$\sigma_{sp}$ 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 (Figure 1d). $\sigma_{sp}$ is a discrete variable over the domain of possible muscle representation sizes $D_{sp}$ (in grid units). Here, we set $D_{sp} = \{0.2, 0.4, \dots, 2\}$.

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

$$P(HS \ \mu_{HS} \ \sigma_{sp} \ S_{1:T} \ \mu_{1:T}), \qquad (4)$$

where $S_{1:T}$ and $\mu_{1:T}$ denotes the series of variables $S_i$ and $\mu_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:

$$P(HS \ \mu_{HS} \ \sigma_{sp} \ S_{1:T} \ \mu_{1:T}) = P(HS) P(\mu_{HS}) P(\sigma_{sp}) \prod_{i=1}^{T} P(S_{1:T}) \prod_{i=1}^{T} P(\mu_i | A), \qquad (5)$$

where $A$, for notational simplicity, is the conjunction of variables HS, $\mu_{HS}$, $\sigma_{sp}$ and $S_i$. The priors on $\mu_{HS}$ and $\sigma_{sp}$ are kept uninformative because of high inter-subject variability regarding these values. $P(\mu_{HS})$
<figure><img src="image_1.png" /><figcaption>Figure 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 ($\mu_i$) modulation as a function of the distances $dX$ and $dY$ between $S_i$ and HS along the X and Y-axis respectively.</figcaption></figure>

and $P(\sigma_{sp})$ thus are uniform distributions over their respective domains. The prior on the stimulation position $S_i$ is a Dirac distribution centered on the position selected during the previous iteration (see 2.1.5). $S_1$ is initialized on the center of the stimulation grid. Finally, $P(\mu_i | A)$ corresponds to the probability of measuring the MEP mean amplitude $\mu_i$ given the coordinates of the stimulated site $S_i$ and HS, the estimated maximal amplitude $\mu_{HS}$ and spatial spread $\sigma_{sp}$ of the muscle representation. We can write our knowledge $\mu_i$ given A in $P(\mu_i | A)$, by merging Eq. (2) and (3), so that:

$$ 
\begin{align} 
P(\mu_i | A) &\propto \frac{1}{\sigma_i \sqrt{2\pi}} \exp\left(-\frac{(\mu - \mu_i)^2}{\sigma_i^2}\right) \tag{6} \\ 
&\propto \frac{1}{\sigma_i \sqrt{2\pi}} \exp\left(-\frac{\left(\mu - \mu_{HS} \exp\left(-\left(\frac{(x_i - x_{HS})^2}{2\sigma_{sp}^2} + \frac{(y_i - y_{HS})^2}{2\sigma_{sp}^2}\right)\right)\right)}{\sigma_i^2}\right). \nonumber 
\end{align} 
$$
### 2.1.3. Question and inference

The objective of AutoHS is to estimate the position of the HS, given all the MEP mean amplitudes $\mu_{1:T}$ observed after stimulating cortical positions $S_{1:T}$. However, this is not possible without jointly inferring $\mu_{HS}$ and $\sigma_{sp}$, since these variables are also conditioning the probability $P(\mu_i | A)$. Therefore, the joint probability $P(HS \mu_{HS} \sigma_{sp} | \mu_{1:T} S_{1:T})$ is computed, which is a four-dimensional probability distribution. Bayesian inference yields (see appendix for the complete derivation):

$$P(HS \mu_{HS} \sigma_{sp} | S_{1:T} \mu_{1:T}) \propto P(HS) \prod_{i=1}^{T} P(\mu_i | A). \quad (7)$$

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 \mu_{HS} \sigma_{sp} | \mu_{1:T} S_{1: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):

$$ 
\begin{align}
P(HS | \mu_{HS} \sigma_{sp} S_{1:T} \mu_{1:T}) &\propto P(HS) \sum_{\mu_{HS}} \sum_{\sigma_{sp}} \prod_{i=1}^{T} P(\mu_i | A) \nonumber \\
P(\mu_{HS} | HS \sigma_{sp} S_{1:T} \mu_{1:T}) &\propto P(HS) \sum_{HS} \sum_{\sigma_{sp}} \prod_{i=1}^{T} P(\mu_i | A) \tag{8} \\
P(\sigma_{sp} | \mu_{HS} HS S_{1:T} \mu_{1:T}) &\propto P(HS) \sum_{\mu_{HS}} \sum_{HS} \prod_{i=1}^{T} P(\mu_i | A) \nonumber
\end{align}
$$

### 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 $i^{th}$ iteration is assessed through the entropy $h_i$ of the joint distribution inferred in Eq. (7):

$$h_i = - \sum_{HS} \sum_{\mu_{HS}} \sum_{\sigma_{sp}} P(HS \mu_{HS} \sigma_{sp} | S_{1:i} \mu_{1:i}) \ln(P(HS \mu_{HS} \sigma_{sp} | S_{1:i} \mu_{1:i})). \quad (9)$$

By definition, the entropy $h_i$ starts at its theoretical maximum $h_{max}$ 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 $h_i$. This threshold
is expressed as a proportion of $h_{max}$. The process thus stops after the $i^{th}$ iteration whenever:

$$h_i < \alpha h_{max}, \qquad (10)$$

with $\alpha \in (0; 1)$.

However, the situation where $h_i$ 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 ($\sigma_i$). The second stop criterion is thus based on the first derivative $dh_i/di$.

The process stops whenever the model does not get any further information from the data, i.e.:

$$\frac{dh_i}{di} < \beta, \qquad (11)$$

where $\beta$ 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 $\alpha = 0.3$ and $\beta = 0.2$.

### 2.1.5. Smart planning of $S_{i+1}$

AutoHS integrates a smart selection of the next stimulation position $S_{i+1}$. Although any procedure of position selection within the stimulation grid would be possible, a prospective choice of $S_{i+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:

$$\forall S_{i+1} | (x_{i+1}, y_{i+1}) \in (\{1, \dots, N_x\}, \{1, \dots, N_y\}),$$

$$h_{i+1} = - \sum_{HS} \sum_{\mu_{HS}} \sum_{\sigma_{sp}} P(HS \ \mu_{HS} \ \sigma_{sp} | S_{1:i+1} \ \mu_{1:i+1}) \ln(P(HS \ \mu_{HS} \ \sigma_{sp} | S_{1:i+1} \ \mu_{1:i+1})), \quad (12)$$

where $\mu_{i+1}$ is directly estimated using Eq. (3) and the current knowledge at iteration $i$ about HS, $\mu_{HS}$ and $\sigma_{sp}$, and $\sigma_{i+1}$ is equal to $\sigma_i$. AutoHS finally chooses the cortical position $S_{i+1}$ inducing the maximal entropy decrease between iterations $i$ and $i+1$:

$$S_{i+1} | \max_{S_{i+1}}(h_i - h_{i+1}). \qquad (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. *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).

### 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 assess 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.

### 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 ($HS_{anat}$) 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 $HS_{anat}$ and oriented with respect to the central sulcus (Figure 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 cm$^2$ 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.

### 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 $HS_{anat}$ 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.

### 2.2.5. Automatized set-up procedure

The automatized set-up procedure (aSUP) was completed by the robotized TMS system together with AutoHS (Figure 2). The procedure began with an initial threshold $rMT_{anat}$ estimated on $HS_{anat}$. HS hunting then started on $HS_{anat}$, using an exploration intensity $I$ of 110% of $rMT_{anat}$. 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 $HS_{anat}$ was too far from the actual HS position. On each stimulated position, five MEP amplitudes (see 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 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.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
<figure><img src="image_2.png" /><figcaption>Figure 2: Automatized set-up procedure. The supraliminal stimulation intensity *I* used for HS hunting is first derived from an initial rMT estimated on *HS*<sub>anat</sub>. 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.</figcaption></figure>

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.

### 2.2.7. TMS parameters

TMS was delivered by a butterfly coil (2\*75 mm) MagPro Cool B65-RO (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 ten-second 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.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 -200ms to +600ms window surrounding stimulation onset.

## 2.3. Data processing

### 2.3.1. Data of interest

The following data were investigated for each procedure:

* 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: $\sqrt{\sum_{x,y} (HS_1 - HS_2)^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 $S_i$, *N* MEPs were drawn from a Gaussian distribution $N(\mu_{S_i}, \sigma_{S_i}^2)$, where $\mu_{S_i}$ and $\sigma_{S_i}$ corresponded respectively to the mean and standard deviation of the 10 MEPs recorded on $S_i$ 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 - \epsilon, 1 + \epsilon)$ where $\epsilon$ 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 $y_i$ was thus calculated as follows:

$$y'_{i} \sim N(\mu_{S_{i}}, \sigma_{S_{i}}^{2})$$

$$n_{i} \sim U(1 - \epsilon, 1 + \epsilon)$$

$$y_{i} = y'_{i} \times n_{i}.$$

### 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 ($HDI_{95}$). $HDI_{95}$ 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): $BF_{01}$, $BF_{10}$ (which is the multiplicative inverse of $BF_{01}$), and $BF_{incl}$ 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.).

## 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
274

excluded from the experimental comparison, as experimenters did not follow the stimulation angle constraint.

275

Results regarding the reliability and reproducibility of AutoHS on virtual data are first described, followed
by the experimental comparison between automatized and manual procedures.

277

3.1. *Model validation*

278

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.

280

282

283

Model reliability

284

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).
Figure 3a and 3b show 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, *μ*<sub>HS</sub> and *σ*<sub>sp</sub> values. The motor mapping data used in this example is
plotted in Figure 3c. The estimated MEP maximum (*μ*<sub>HS</sub> = 1,750 *μ*V) was close to the recorded maximum
(1,689 *μ*V), and the spatial spread of the muscle representation was correctly estimated (*σ*<sub>sp</sub> = 8.4 mm) by
the model.

285

286

287

288

289

290

291

292

Spatial reproducibility

293

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 (Figure 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 = \frac{\sum (x\mu)}{\sum \mu}$, and
$y = \frac{\sum (y\mu)}{\sum \mu}$ as coordinates, with $\mu$ being the mean MEP amplitude for each point of the grid.

294

295

296

297

298

299

300

301

302

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 (BF<sub>incl</sub> = +inf, HDI<sub>95</sub> =

303
<figure><img src="image_3.png" /><figcaption>Figure 3: Inferred probability distributions throughout AutoHS, for one virtual motor map of a representative subject (S16). a: (top) Evolution of the marginal distribution $P_m(\text{HS})$ (as defined by Eq.(8)) over time. Iteration 0 represents the prior distribution $P(\text{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 $\alpha h_{\text{max}}$ (see text). b: From left to right, evolution over time of the marginal distributions $P_m(\sigma_{\text{sp}})$ and $P_m(\mu_{\text{HS}})$ respectively. For this example, AutoHS estimated $\mu_{\text{HS}} = 1,750$ µV and $\sigma_{\text{sp}} = 8.4$ mm. c: Final motor mapping. The black area corresponds to the spatial distribution (1 $\sigma$ ellipse) of the HS positions found during the 100 simulations ($N = 5, \epsilon = 0$).</figcaption></figure>

[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$, $BF_{\text{incl}} = +\inf$) and noise-free data (main effect of $\epsilon$, $BF_{\text{incl}} = 4e^{+11}$). It appeared that the global maximum search procedure was more sensitive to the noise parameter than AutoHS (interaction effect, $BF_{\text{incl}} = 1e^{+7}$), whereas both methods seemed to be equally influenced by the number of MEPs ($BF_{01} = 5.56$). Figure 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. Figure 4c shows
<figure><img src="image_4.png" /><figcaption>Figure 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.</figcaption></figure>

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.2. *Experimental comparison between automatized and manual procedures*

### HS shift

Figure 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 (Figure 5b). The paired comparison of these two distributions failed to draw any conclusion regarding their difference ($BF_{10} = 0.4, BF_{01} = 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
a

<figure><img src="image_5.png" /><figcaption>Figure 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.</figcaption></figure>

b

<figure><img src="image_6.png" /></figure>

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 (Figure 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.

rMT

Figure 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 ($BF_{01}$ = 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 (Figure 6c). However, the most probable values of the correlation and its evidence were higher for aSUP ($HDI_{95}(\rho)$ = [0.73; 0.97], $BF_{10}$ = $6e^4$) than for mSUP ($HDI_{95}(\rho)$ = [0.43; 0.91], $BF_{10}$ = 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 ($BF_{01}$ = 3.45).
<figure><img src="image_7.png" /><figcaption>Figure 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.</figcaption></figure>

Estimated duration

Figure 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 (Figure 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 ($HDI_{95} = [1.6; 4.6]$ min, $BF_{10} = 273$). On average, mSUP lasted $11.5 \pm 3.7$ min while aSUP lasted $8.3 \pm 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 ($BF_{incl} = 1.16$, $BF_{01} = 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 (BF<sub>01</sub> = 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 (BF<sub>01</sub> = 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 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, 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 Figure 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 $\sigma_{prior}$ (14 mm), the prior was kept uninformative so that its influence was 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 $\sigma_{prior}$ could then increase the convergence speed (about 3 iterations for $\sigma_{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 ($\alpha$ and $\beta$) would allow decreasing procedure time, as they would be reached faster, while lowering the reproducibility of the calculation (see the example depicted in Figure 3a).
## 4.2. Redefining the HS

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 Figure 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 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 Figure 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 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 et al., 2015). Second, the parametrical form of the distribution of MEP amplitudes knowing the positions of HS and the stimulation point ($P(\mu_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 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.

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.

6. Acknowledgments

This work was funded by the Pôle Grenoble Cognition grant “PGC-2014” and the Agence Nationale
pour la Recherche grant “ANR-15-CE37-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.
## 7. Appendix

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

$$
\begin{align*}
& P(HS \ \mu_{HS} \ \sigma_{sp} | S_{1:T} \ \mu_{1:T}) \\
&= \frac{P(HS \ \mu_{HS} \ \sigma_{sp} \ S_{1:T} \ \mu_{1:T})}{P(S_{1:T} \ \mu_{1:T})} \tag*{[Bayes theorem]} \\
&= \frac{1}{Z_1} P(HS \ \mu_{HS} \ \sigma_{sp} \ S_{1:T} \ \mu_{1:T}) \tag*{[P(S_{1:T} \ \mu_{1:T}) = Z_1 \ (constant)]} \\
&= \frac{1}{Z_1} P(HS) \ P(\mu_{HS}) \ P(\sigma_{sp}) \prod_{i=1}^{T} P(S_i) \prod_{i=1}^{T} P(\mu_i | HS \ \mu_{HS} \ \sigma_{sp} \ S_i) \tag*{[by definition of the joint]} \\
&= \frac{1}{Z_1 Z_2} P(HS) \prod_{i=1}^{T} P(S_i) \prod_{i=1}^{T} P(\mu_i | HS \ \mu_{HS} \ \sigma_{sp} \ S_i) \tag*{[P(\mu_{HS})P(\sigma_{sp}) = Z_2 \ (uniforms)]} \\
&= \frac{1}{Z_1 Z_2} P(HS) \prod_{i=1}^{T} P(\mu_i | HS \ \mu_{HS} \ \sigma_{sp} \ S_i) \tag*{[P(S_i) = 1 \ (Dirac)]} \\
&\propto P(HS) \prod_{i=1}^{T} P(\mu_i | HS \ \mu_{HS} \ \sigma_{sp} \ S_i) \tag*{[\hspace{2.5cm} \square]} \\
&\propto P(HS) \prod_{i=1}^{T} P(\mu_i | A) \tag{14}
\end{align*}
$$

recalling that A is the conjunction of variables HS, $\mu_{HS}$, $\sigma_{sp}$ and $S_i$.

## 8. References

Ahdab, R., Ayache, S. S., Brugières, P., Farhat, W. H., Lefaucheur, J.-P., Mar. 2016. The Hand Motor Hotspot is not Always Located in the Hand Knob: A Neuronavigated Transcranial Magnetic Stimulation Study. Brain Topography, 1–8.

Awiszus, F., 2003. TMS and threshold hunting. Supplements to Clinical Neurophysiology 56, 13–23.

Awiszus, F., Borckardt, J., 2011. TMS motor threshold assessment tool (MTAT 2.0).

Baek, J., Lesmes, L. A., Lu, Z.-L., Aug. 2016. qPR: An adaptive partial-report procedure based on Bayesian inference. Journal of Vision 16 (10), 25.

Barker, A. T., Jalinous, R., Freeston, I. L., May 1985. Non-invasive magnetic stimulation of human motor cortex. The Lancet 325 (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, doi: 10.1007/978-3-540-79007-5_2.
Bestmann, S., Feredoes, E., Aug. 2013. Combined neurostimulation and neuroimaging in cognitive neuroscience: past, present, and future: Combined neurostimulation and neuroimaging. Annals of the New York Academy of Sciences 1296 (1), 11–30.

Bortoletto, M., Veniero, D., Thut, G., Miniussi, C., Feb. 2015. The contribution of TMS-EEG coregistration in the exploration of the human cortical connectome. Neuroscience & Biobehavioral Reviews 49, 114–124.

Cincotta, M., Giovannelli, F., Borgheresi, A., Balestrieri, F., Toscani, L., Zaccara, G., Carducci, F., Viggiano, M. P., Rossi, S., Apr. 2010. Optically tracked neuronavigation increases the stability of hand-held focal coil positioning: Evidence from “transcranial” magnetic stimulation-induced electrical field measurements. Brain Stimulation 3 (2), 119–123.

Desmurget, M., Richard, N., Harquel, S., Baraduc, P., Szathmari, A., Mottolese, C., Sirigu, A., Apr. 2014. Neural representations of ethologically relevant hand/mouth synergies in the human precentral gyrus. Proceedings of the National Academy of Sciences 111 (15), 5718–5722.

Dugué, L., Marque, P., VanRullen, R., Aug. 2011. The Phase of Ongoing Oscillations Mediates the Causal Relation between Brain Excitation and Visual Perception. The Journal of Neuroscience 31 (33), 11889–11893.

Farzan, F., Vernet, M., Shafi, M. M. D., Rotenberg, A., Daskalakis, Z. J., Pascual-Leone, A., Sep. 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. 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 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). IEEE, pp. 5352–5355.

Giszter, S. F., Aug. 2015. Motor primitives — new data and future questions. Current Opinion in Neurobiology 33, 156–165.

Graziano, M. S. A., Feb. 2016. Ethological Action Maps: A Paradigm Shift for the Motor Cortex. Trends in Cognitive Sciences 20 (2), 121–132.

Graziano, M. S. A., Aflalo, T. N., Oct. 2007. Mapping Behavioral Repertoire onto the Cortex. Neuron 56 (2), 239–251.

Gugino, L. D., Rafael Romero, J., Aglio, L., Titone, D., Ramirez, M., Pascual-Leone, A., Grimson, E., Weisenfeld, N., Kikinis, R., Shenton, M.-E., Oct. 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. Clinical Neurophysiology 112 (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., May 2016a. Mapping dynamical properties of cortical microcircuits using robotized TMS and EEG: Towards functional cytoarchitectonics. NeuroImage 135, 115–124.

Harquel, S., Beynel, L., Guyader, N., Marendaz, C., David, O., Chauvin, A., Oct. 2016b. CortExTool: a toolbox for processing motor cortical excitability measurements by transcranial magnetic stimulation.

URL https://hal.archives-ouvertes.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., Jul. 2009. Motor Threshold in Transcranial Magnetic Stimulation: The Impact of White Matter Fiber Orientation and Skull-to-Cortex Distance. Human brain mapping 30 (7), 2044–2055.
Herring, J. D., Thut, G., Jensen, O., Bergmann, T. O., Oct. 2015. Attention Modulates TMS-Locked Alpha Oscillations in the Visual Cortex. The Journal of Neuroscience 35 (43), 14435-14447.

Herwig, U., Schönfeldt-Lecuona, C., Wunderlich, A. P., von Tiesenhausen, C., Thielscher, A., Walter, H., Spitzer, M., Nov. 2001. The navigation of transcranial magnetic stimulation. Psychiatry Research: Neuroimaging 108 (2), 123-131.

Janssen, A. M., Oostendorp, T. F., Stegeman, D. F., May 2015. The coil orientation dependency of the electric field induced by TMS for M1 and other brain areas. Journal of NeuroEngineering and Rehabilitation 12.

Jeffreys, H., Aug. 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., Feb. 2009. Comparison of navigated and non-navigated transcranial magnetic stimulation for motor cortex mapping, motor threshold and motor evoked potentials. NeuroImage 44 (3), 790-795.

Jung, N. H., Delvendahl, I., Kuhnke, N. G., Hauschke, D., Stolle, S., Mall, V., Apr. 2010. Navigated transcranial magnetic stimulation does not decrease the variability of motor-evoked potentials. Brain Stimulation 3 (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., Nov. 2009. Robot-assisted image-guided transcranial magnetic stimulation for somatotopic mapping of the motor cortex: a clinical pilot study. Acta Neurochirurgica 152 (2), 333-343.

Kruschke, J., Nov. 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., Nov. 2014. Evidence-based guidelines on the therapeutic use of repetitive transcranial magnetic stimulation (rTMS). Clinical Neurophysiology 125 (11), 2150-2206.

Lefaucheur, J.-P., Picht, T., Apr. 2016. The value of preoperative functional cortical mapping using navigated TMS. Neurophysiologie Clinique/Clinical Neurophysiology 46 (2), 125-133.

Meesen, R. L., Cuypers, K., Rothwell, J. C., Swinnen, S. P., Levin, O., Jun. 2011. The effect of long-term TENS on persistent neuroplastic changes in the human cerebral cortex. Human Brain Mapping 32 (6), 872-882.

Meincke, J., Hewitt, M., Batsikadze, G., Liebetanz, D., Jan. 2016. Automated TMS hotspot-hunting using a closed loop threshold-based algorithm. NeuroImage 124, Part A, 509-517.

Mulder, J., Wagenmakers, E.-J., Jun. 2016. Editors' introduction to the special issue "Bayes factors for testing hypotheses in psychological research: Practical relevance and new developments". Journal of Mathematical Psychology 72, 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., Apr. 2004. A stereotactic method for image-guided transcranial magnetic stimulation validated with fMRI and motor-evoked potentials. NeuroImage 21 (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., Oct. 2015. Bringing transcranial mapping into shape: Sulcus-aligned mapping captures motor somatotopy in human primary motor hand area. NeuroImage 120, 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., Apr. 2013. Optimal Coil Orientation for Transcranial Magnetic Stimulation. PLOS ONE 8 (4), e60358.

Rogasch, N. C., Fitzgerald, P. B., Jul. 2013. Assessing cortical network properties using TMS-EEG. Human Brain Mapping 34 (7), 1652-1669.

Rossi, S., Hallett, M., Rossini, P. M., Pascual-Leone, A., Dec. 2009. Safety, ethical considerations, and application guidelines for the use of transcranial magnetic stimulation in clinical practice and research. Clinical Neurophysiology 120 (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., Aug. 1994. Non-invasive 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), 79-92.

Schieber, M. H., Nov. 2001. Constraints on Somatotopic Organization in the Primary Motor Cortex. Journal of Neurophysiology 86 (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 Neuroscience 14 (1), 94.

Sparing, R., Buelte, D., Meister, I. G., Pauš, T., Fink, G. R., Jan. 2008. Transcranial magnetic stimulation and the challenge of coil placement: A comparison of conventional and stereotaxic neuronavigational strategies. Human Brain Mapping 29 (1), 82-96.

Stokes, M. G., Chambers, C. D., Gould, I. C., Henderson, T. R., Janko, N. E., Allen, N. B., Mattingley, J. B., Dec. 2005. Simple Metric For Scaling Motor Threshold Based on Scalp-Cortex Distance: Application to Studies Using Transcranial Magnetic Stimulation. Journal of Neurophysiology 94 (6), 4520-4527.

Thielscher, A., Opitz, A., Windhoff, M., Jan. 2011. Impact of the gyral geometry on the electric field induced by transcranial magnetic stimulation. NeuroImage 54 (1), 234-243.

van de Ruit, M., Perenboom, M. J. L., Grey, M. J., 2015. TMS Brain Mapping in Less Than Two Minutes. Brain Stimulation 8, 231-239.

Wassermann, E., Epstein, C. (Eds.), Nov. 2012. The Oxford handbook of transcranial stimulation. Oxford University Press, Great Clarendon Street, Oxford, OX2 6DP.

Wassermann, E. M., Jul. 2002. Variation in the response to transcranial magnetic brain stimulation in the general population. Clinical Neurophysiology 113 (7), 1165-1171.

Weiss, C., Nettekoven, C., Rehme, A. K., Neuschmelting, V., Eisenbeis, A., Goldbrunner, R., Grefkes, C., Feb. 2013. Mapping the hand, foot and face representations in the primary motor cortex — Retest reliability of neuronavigated TMS versus functional MRI. NeuroImage 66, 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.