Accounting for Stimulations That Do Not Elicit Motor-Evoked Potentials When Mapping Cortical Representations of Multiple Muscles

The representation of muscles in the cortex can be mapped using navigated transcranial magnetic stimulation. The commonly employed measure to quantify the mapping are the center of gravity or the centroid of the region of excitability as well as its size. Determining these measures typically relies only on stimulation points that yield motor-evoked potentials (MEPs); stimulations that do not elicit an MEP, i.e., non-MEP points, are ignored entirely. In this study, we show how incorporating non-MEP points may affect the estimates of the size and centroid of the excitable area in eight hand and forearm muscles after mono-phasic single-pulse TMS. We performed test-retest assessments in twenty participants and estimated the reliability of centroids and sizes of the corresponding areas using inter-class correlation coefficients. For most muscles, the reliability turned out good. As expected, removing the non-MEP points significantly decreased area sizes and area weights, suggesting that conventional approaches that do not account for non-MEP points are likely to overestimate the regions of excitability.


INTRODUCTION
Single-pulse transcranial magnetic stimulation (TMS) is a non-invasive and painless technique that monitors neurophysiological alterations of the human motor cortex (Barker et al., 1985;Schambra et al., 2003). A TMS coil discharge at suitable intensity will induce transient currents and cause depolarization of axons of nerve cells (Rossini et al., 2015). When applied over the motor cortex, this can elicit a motor-evoked potential (MEP) that can be recorded in contralateral target muscles using conventional electromyography (EMG). Amplitudes and latencies of the MEPs reveal the excitability and conduction times of the cortical-spinal tract. Both have been conceived as valid outcomes of TMS motor mapping (Rossini et al., 1994). Neuroscientists and physicians alike utilized TMS motor mapping to evaluate muscle synergies and motor cortical plasticity (Siebner and Rothwell, 2003), to plan brain tumor surgery (Krieg et al., 2012), or to follow recovery after stroke (Mark et al., 2006;Sondergaard et al., 2021). There is ample evidence that the location at which TMS elicits the maximum MEP is particularly close to the location found using direct cortical stimulation, which is considered the gold standard in motor mapping. The localization clearly outperforms other modalities like magnetoencephalography (Tarapore et al., 2012) or functional magnetic resonance imaging (Forster et al., 2011).
TMS combined with neuro-navigation increases mapping accuracy (Krieg et al., 2017). A popular approach for coil positioning is a pseudo-random walk. Delivering stimulations at random locations roughly evenly spaced over the motor cortex is more efficient and potentially more accurate than timeconsuming, course-grained grid-based positioning. Likely, this random placement will also elicit MEPs in muscles other than the target muscle. Often considered a confounder, this isin fact-particularly useful when multiple muscles are being evaluated, presuming that the muscles have similar resting motor thresholds (Krieg et al., 2017) and close-by cortical representations (Schieber, 2001). Very recently, Tardelli et al. (2021) used the pseudo-random walk method to assess the cortical representation of abductor digiti minimi, flexor carpi radialis, and flexor pollicis brevis. As a rule of thumb: the more muscles are measured simultaneously, the more efficient assessments via pseudo-random coil positioning can be.
Irrespective of the experimental protocol, navigated TMS derived cortical map outcomes should have good reliability (Novikov et al., 2018). Nonetheless, "even most commonly used outcomes such as areas, volumes, the location of centers of gravity (CoGs), and hotspots have (hardly) been validated for being reliable measures in test-retest studies (Kraus and Gharabaghi, 2016)." We slightly modified this quote from Novikov et al. (2018), because they and other likewise recent reports did indeed test for the reliability of navigated TMS outcomes considered in the respective studies. For instance, Nazarova et al. (2021) evaluated the reliability of the CoG and the size of the area (volume) of excitability, next to the position of the MEP hotspot. In a grid-based approach, all measures displayed high relative but low absolute reliability, with the latter arguably reflecting between-subject variability.
The area of excitability can be defined as the cortical region within which TMS elicits an MEP. It is usually determined by projecting the focal point of the coil's magnetic field on a (re-)constructed spherical surface or volume and determining the resulting convex hull. State-of-the-art fine-tuning of this approach is to a priori concentrate on the cortical patch of interest. However, a mapping that agrees with the "real" anatomical structure generally provides better area estimates. This is particularly true when realizing that the gray matter border may have large curvatures along gyral ridges (Van Essen, 2004), where spherical approximations will be poor. We followed these lines and extracted subject-specific cortical surfaces at high resolution, projected the stimulation points to that surface and estimated the area spanned by the pair-wise shortest paths connecting the stimulation points. More importantly, we also projected the stimulation points where TMS did not elicit an MEP and removed these points from the estimated area. As will be shown, together these steps circumvent potential over-estimation of the area of excitability.

Participants
Twenty healthy, right-handed volunteers (average age: 29.6 ± 7.5, eight females) participated in the study. Prior to the experiment, all participants were screened for contraindications of MRI and TMS through questionnaires (Rossi et al., 2011). All of them provided signed informed consent prior to joining the experimental sessions. The Edinburgh Handedness Inventory served to determine hand dominance (Oldfield, 1971). The study had been approved by the medical ethics committee of Amsterdam University Medical Center (VUmc, 018.213-NL65023.029.18).

Materials
Our set-up consisted of three devices: a TMS system, an EMG amplifier, and a neural navigation system. Single-pulse TMS was delivered by a Magstim 200 2 stimulator (Magstim Company Ltd., Whitland, Dyfed, United Kingdom) using a figure-ofeight coil with 70 mm windings. Eight bipolar EMG signals were recorded using a 16-channel EMG amplifier (Porti, TMSi, Oldenzaal, Netherlands) and continuously sampled at a rate of 2 kHz. The EMG recordings were triggered by the TMS to allow for online EMG assessments using a custom-made Labview-program with embedded Matlab functions (designed at our department using Labview 2016, National Instruments, Austin, TX, and Matlab 2018b, The MathWorks, Natick, MA). In brief, upon receiving a trigger, peak-to-peak amplitudes and latencies of MEPs were estimated from all EMG signals during the following 500 ms. These outcomes, as well as the original EMG signals (duration = 500 ms), were sent to the neural navigation system (Neural Navigator, Brain Science Tools, De Bilt, Netherlands) 1 for online monitoring and storage. The neural navigation software also stored the position and orientation of the coil with respect to the head.
Prior to running the TMS protocol, we acquired the participants' anatomical T1-weighted MRI (3 Tesla Philips Achieva System, Philips, Best, Netherlands; matrix size 256 × 256 × 211, voxel size 1.0 × 1.0 × 1.0 mm 3 , TR/TE 6.40/2.94 ms). For online neuro-navigation, gray matter was segmented using SPM 2 ; note that for offline analysis, we employed a more detailed segmentation via Freesurfer. 3 We considered the first dorsal interosseous (FDI), abductor digiti minimi (ADM), abductor pollicis brevis (APB), flexor pollicis brevis (FPB), extensor digitorum communis (EDC), flexor digitorum superficialis (FDS), extensor carpi radialis (ECR), and flexor carpi radialis (FCR) muscles, which were measured using bipolar electrodes (Blue Sensor N-00-S, Ambu, Ballerup, Denmark), placed after cleaning the skin with alcohol (cf. Figure 1). The ground electrode was attached to the ulnar styloid process. We monitored and kept the electrode impedance below 5 k . During the experiment, the orientation of the TMS coil was held 45 degrees to the sagittal plane, and tangential to the scalp. By this, we meant to induce currents in the cortex along the posterior-to-anterior direction. To control the TMS output, we used the Matlab-toolbox Rapid2. 4

Experimental Procedures
Participants were seated comfortably in an armchair, relaxing muscles of hands and arms. The experiment consisted of two identical sessions, Session 1 and Session 2. These sessions were separated by 1h and served to test for test-retest reliability of our outcomes. EMG electrodes were kept fixed to minimize placement errors. The interval of 1h was set to prevent drying of the conductive electrolyte gel.
In each session, we searched for the hotspot positions for FDI, EDC and FCR before testing the RMTs. First, the stimulation intensity was identified that yielded MEPs for all three muscles when stimulating in the omega-shaped area ("hand knob") of the precentral gyrus. We started at 45% of the maximum stimulator output and increased or decreased the intensity until a consistent MEP was present. Then, we performed thirty stimulations around the hand knob region along the precentral gyrus. From these stimulations, we determined the position with the largest peakto-peak amplitude for every muscle and labeled that position as the hotspot. Next, the RMT for every muscle was determined at the muscle-specific hotspot as the minimum stimulator output at which peak-to-peak amplitudes exceeded 50 µV in five out of ten stimuli. This was followed by the actual mapping procedure.
The TMS coil was pseudo-randomly positioned such that stimulations covered the entire left precentral gyrus. We applied 120 stimulations (Cavaleri et al., 2017) and repeated this at three intensities: 105% RMT of FDI, EDC, and FCR, respectively. In total, we performed 360 stimulations in every session. We chose 105% RMT because previous studies suggested it to be the lowest possible intensity for upper limb muscles mapping (Krieg et al., 2017), thus leading to the least stimulation cross-talk. Finally, we estimated the hotspots of the other five muscles (ADM, APB, FPB, FDS, and ECR) and determined the respective RMTs.

Motor-Evoked Potentials Definition
We discriminated between TMS with and without eliciting MEPs, i.e., MEP and non-MEP points. MEPs were considered proper if 4 https://github.com/armanabraham/Rapid2 their amplitude exceeded 20 times the EMG-baseline's standard deviation (defined over 100 ms prior to each stimulation). While on average these thresholds were [51, 51, 76, 63, 59, 55, 70, and 64] µV for FDI, ADM, APB, FPB, EDC, FDS, ECR, and FCR, respectively, the baseline's standard deviations differed substantially over the group rendering a subject-specific threshold definition appropriate-see Supplementary Figure 2 for the corresponding boxplots and median values. Amplitudes were also required to stay below 10 mV (to exclude movement and cable artifacts) and the peak's latency had to fall within the range of 5-50 ms after stimulation. All other stimulations were marked as non-MEP points; see below under Outcome measures for further details.

Area Estimate
The area estimates were based on triangulated cortical surface meshes that we extracted using Freesurfer (see text footnote 3; version 7). We imported the meshes into Brainstorm (version 3) 5 to ease converting between world and subject-specific MRI coordinates. Next to the original meshes with about 230,000-340,000 vertices dependent on the participant, we also generated low-resolution version by downsampling the mesh to either 15,000 or 100,000 vertices. This enabled us to test for effects of surface resolution. In all cases, we assigned the Mindboggle anatomical atlas (version 6; see also Klein et al., 2017) 6 to select left primary motor cortex. The area construction consisted of four steps: (i) Stimulation points were projected to the triangulated cortical surface mesh yielding a set of vertices as illustrated in Figure 2 and further detailed under Outcome measures. (ii) Vertices that did not fall in left primary motor cortex were excluded (label "precentral L" in the Mindboggle6 atlas). (iii) We connected the vertices along their shortest connecting paths. In brief, we converted the mesh into a sparse, weighted graph. The edges of triangularization served as adjacencies that we weighted by the Euclidean distance between the corresponding vertices. Then, we searched for the shortest paths between all MEP points (Dijkstra, 1959). We repeated this iteratively for all points of the connecting paths until no points were added. Figure 3A briefly summarizes this iteration. Further details can be found at github. 7 (iv) Finally, we excluded the vertices (and triangles) corresponding to non-MEP points from the resulting area as sketched in Figure 3B. Note that if a triangle contained both, one or more stimulations that elicited a MEP and one or more stimulations that did not, we kept the triangle. By this we limited the risk of underestimation that may stem only from falsely considering stimulation points as non-MEP points.
By construction, removing non-MEP points will reduce the size of the active areas. To appreciate the benefits of non-MEP point removal, consider the case in which the true area of excitability has a non-convex boundary, e.g., if the area is U-shaped. Conventional estimates, in particular ones based on estimating the convex hull of the cloud of stimulation points, will clearly provide an overestimate of the excitable area. Rather than opting for non-convex hull estimates, we used a more general approach that also allows for removing points that are scattered across the area spanned by MEP points. We briefly illustrate this in Figure 4 showing data of a single subject where isolated triangles that are being removed.
From here on, we refer to the reconstruction without accounting for non-MEP points as method M1, whilst the removal of non-MEP points will be method M2 (with examples in Figures 4A,B, respectively).

Outcome Measures
Most of the outcome measures were based on the MEP amplitudes and latencies. We quantified them for every TMS pulse and for every muscle using the original EMG signals, from which we removed the stimulation artifact via linear interpolation (-1 to +2 ms around stimulation) followed by high pass filtering at 10 Hz (2nd order, bi-directional Butterworth design). We defined the epoch -100 to -1 ms before the stimulus as a baseline and determined its mean value µ and 7 https://github.com/marlow17/surfaceanalysis the standard deviation σ. A peak in the interval 5-100 ms after the stimulus was considered an MEP if its value exceeded µ ± 20σ. Its latency was set as the first sample after the stimulus at which the signal exceeded µ + 2σ (µ − 2σ) if the first peak is a maximum (minimum). Finally, we removed MEPs with peak-to-peak amplitude larger than 10 mV as we considered them artifacts. Except for the latter artifact definition, we opted for relative-to-baseline changes when identifying MEPs to circumvent between-subject variability in skin conductance; the choices for ± 20σ and ± 2σ were based on visually inspecting the EMG traces.
The stimulation points were mapped onto cortical surfaces given as triangulated meshes. The triangle of the surface mesh closest to a stimulation point along the direction of coil orientation was determined following the approach by Möller and Trumbore (1997) (see Figure 2). We assigned the vertices v 0 , v 1 , and v 2 of the closest triangle the corresponding amplitude value, a. If two or more stimulations with an MEP shared a vertex, we averaged their amplitudes at the shared point. Since the total area of excitability possibly covered points that were not projected directly, we set all amplitude values via natural interpolation with C 1 continuity toâ; note that if interpolation was not needed, i.e., at the original vertices v 0 , v 1 , and v 2 , thenâ = a. For every triangle we defined the length between their vertices as ||· · · || denoting the Euclidean distance between vertices. That is, given Cartesian coordinates v i = x i , y i , z i , we used, e.g., The total area A of k = 1, ..., M triangles weighted by the MEP amplitudes was computed via a slight modification of Heron's formula, namely the triangular prism, that reads: Frontiers in Human Neuroscience | www.frontiersin.org  and a k = 1 3 2 i = 0âi,k being the mean value of (interpolated) MEP amplitudes at the three vertices of triangle k.
Given all i = 1, ..., N area vertices, we further defined the centroid of the total area in line with the conventional form of the center of gravity (Opitz et al., 2014) as: (2)

Statistics
We first estimated the reliability between Session 1 and Session 2 via the intraclass correlation coefficients (ICC) of the centroid (C), the weighted area size (A) for every muscle and intensity level. In more detail, we use a two-way mixed-effects model for single measurement type and estimated the absolute agreement, i.e., ICC (2,1) conform the Shrout and Fleiss convention (Koo and Li, 2016). We ran this analysis separately for the three representations of the cortex (i.e., three mesh resolutions) and for the three intensities. While the highest resolution for the intensity of 105% RMT of FCR will be reported below, all the other results can be found as Supplementary Material. There we also report the ICCs for centers of gravity (CoG) and for both the MEP amplitudes and the latencies. To further confirm the absence of significant differences between Sessions 1 and 2, we performed a two-way ANOVA with repeated measures including factors of intensity and session. This also allowed for assessing effects of stimulation intensity. We applied a Bonferroni correction for multiple comparisons. Again, we restrict ourselves to reporting "only" the findings of the estimates at maximum resolution in the body text and refer to Supplementary Material for all other cases.
Finally, to assess effects of cortex mesh resolution and of ignoring/removing non-MEP stimulation points, we used a twoway repeated ANOVA with factors method and resolution (again with Bonferroni correction for multiple comparisons).
Prior to conducting the ANOVAs, sphericity was verified via Mauchly's test. A Greenhouse-Geisser correction was performed if necessary. Throughout hypothesis testing, we used a significance threshold of α = 0.05. All statistical analyses were conducted using Matlab (The Mathworks Inc., Natwick MA, version 2020b).

RESULTS
All N = 20 participants completed the experimental procedure without adverse reactions. Of all mappings (subjects × muscle × intensity × session = 960, each containing 120 stimulations) 2% did not contain any valid MEP, and thus did not enter further analyses. In 11/320 (subjects × muscles × session) cases this was for 105% RMT of FDI, 5/320 for 105% RMT of EDC, and 1/320 for 105% RMT of FCR; see Supplementary Table 1). For five subjects, we could not detect any MEPs for ADM when using the second intensity in both sessions.
The ICCs appeared consistent between methods M1 (ignoring non-MEP points) and M2 (removing non-MEP points). Most of them were moderate to good. Good-excellent reliability was found for estimated centroids in the anterior/posterior and superior/inferior directions (x and z coordinates, respectively). While the area sizes' ICCs of FDI, FDS, and FCR were poor, the ANOVA did not reveal any significant differences between the area estimates between sessions. We illustrate this in Table 2 for the highest cortex resolution and refer to Supplementary Tables 3A,B for the ANOVA results for the other cortex     Supplementary Tables 4-6A-C, we also provide the results for the corresponding centroid positions. In a nutshell there were hardly any significant effects of session or intensity (let alone their interaction) on the centroids; when correcting for multiple comparisons all effects will turn out not significant.

resolutions. In
Here we would like to note that this dependency on stimulation intensity can be understood when looking at the effects of intensity on the mere MEP amplitudes, i.e., without projecting them onto the cortex. The corresponding results can be found as Supplementary Table 7. In a nutshell, the amplitudes of FDI, EDC, FCR, FDS, and FCR significantly increased with increasing stimulation intensity.
As expected, ignoring non-MEP points (M1) consistently resulted in larger area sizes when compared to the case when non-MEP points were removed (M2). Our second ANOVA confirmed this. We summarized this in Table 3 where we highlighted the main effects of method. Yet, we also would like to note the interaction effect with resolution, suggesting that the correction for non-MEP points is especially relevant when incorporating low-resolution cortical meshes (see also Figure 5, upper row).
The main effect of method (ignoring non-MEPs vs. removing them) is also illustrated in Figure 6 where we show the relative change in the estimated area sizes. Irrespective of resolution, not removing the non-MEP points yields an overestimation of the active areas, though this effect appears particularly pronounced at low resolution (top panel in Figure 6).
We finally illustrate the effect of removing non-MEP points in Figure 5, where it can be clearly seen that higher resolutions lead to less area being removed.

DISCUSSION
We assessed the reliability of the cortical representation of eight muscles mapped simultaneously using navigated TMS. We distinguished two methods to estimate the active area of a muscle. In the first, more conventional one (M1), we included all stimulation points that elicited an MEP. In the second method (M2), we included the same points but also excluded all stimulation points that did not elicit an MEP. We tested for the effects of the type of measure with the obvious expectation that the latter will yield smaller active areas. We also tested for effects of stimulation intensity and cortical mesh resolution in two consecutive sessions. By and large, we found that the reliabilities of the size and the centroids of the active areas for all the muscles were excellent, good, or moderate. Exceptions were the area size estimates in three muscles ( Table 1) that came with small areas sizes but strong outliers when looking at their representation at high-resolution cortical surface meshes (Figure 6, lower panel). The ICCs of amplitude and latency were excellent or good for all the muscles, again supporting the reliability of our experimental approach (cf. Supplementary Table 8).
One must realize that designing multiple muscle mapping experiments can-in general-be problematic as the RMT of a single muscle must be considered a reference when setting the stimulation intensity. In our case, the difference of RMTs values between different stimulation intensities was small (on average 3.1% of stimulator output; when looking at the individual subjects we found maximum differences of | RMT FDI -RMT EDC | = | RMT FCR -RMT EDC | = 9%, and | RMT FDI -RMT FCR | = 8%). Intensities of 105%, 110-120% (Akiyama et al., 2006) RMT have been widely used in motor mapping (Bohning et al., 2001;Akiyama et al., 2006;Tarapore et al., 2012), suggesting that the here-observed difference is acceptable if not negligible. Hence, forearm and hand muscles might be pooled in a group of muscles with "similar RMTs" and may be evaluated at the same intensity.
For all the muscles, the ICCs of the centroids' positions were moderate to excellent. In the Supplementary Tables 2A-C we show the likewise good results for the more conventional CoGs. Both the estimated centroid as well as the centers-ofgravity hence appeared very consistent and should be considered reliable outcomes in motor mapping, in particular also the anterior/posterior and superior/inferior directions, in line with previous studies (Weiss et al., 2013;Cavaleri et al., 2018). The CoG is commonly employed to quantify the cortical representation of muscles (Massé-Alarie et al., 2017;Nazarova et al., 2021). However, there are several issues with the notion of "CoG" itself. For instance, for many shapes (of cortical representations), the CoG will lie outside the actual stimulation area itself (consider a banana, whose CoG will not be inside the banana itself). The CoG may hence be a tricky measure to give an estimate of the cortical representation of a muscle, especially when the true cortical representation is non-trivially shaped and on a curved surface. Supplementing the CoG, or in our case the  centroid, by the area of excitability is clearly needed, especially when the area estimate is weighted by the MEP amplitude. Again, we advocate incorporating the non-MEP stimulation points in these estimates. When removing non-MEP points, the areas of excitability became significantly smaller than when non-MEP points were simply ignored. We argue that by ignoring the stimulation points that do not elicit MEPs one runs the risk of overestimating the area of excitability and thus to mis-represent muscles in the cortex. Our results show that accounting for non-MEP points does not jeopardize the reliability of assessments. As such we advocate for correcting any potential structural error and provide the tools to do so. Of course, one may counter the fear for structural errors by subsuming that the neuronal population that ought to be covered by our cortical map are likely to be homogeneously distributed. However, several invasive studies already speak against this (e.g., Schieber, 2001). By using intracortical micro-stimulation, Nudo et al. (1996) revealed that the cortical representation of distal forelimb muscles is quite complicated and clearly not uniform. Moreover, to date most area measures rely on estimating convex hulls that clear yield weak approximations if the excitable area has a non-convex boundary-when looking at precentral gyrus that might be the rule rather than the exception.

CONCLUSION
Estimating the active area can be improved when incorporation points at which TMS does not elicit an MEP. Navigated TMS and a pseudo-random coil placement allow for correcting area estimates post-hoc and hence reduce the risk of overestimating the cortical representation of active areas. As such, the very fact that at certain points, a stimulation does not yield a measurable response appears informative. And, even when assessing multiple muscles in unison, this approach comes with high reliability, albeit under the provison that stimulation intensity has been chosen properly.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Medisch Ethische Toetsingscommissie VUmc.
The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
FJ designed and conducted the experiment, analyzed the data, and wrote the draft. SB designed and conducted the experiment, analyzed the data, and modified the draft. AD designed the experiment, analyzed the data, and modified the draft. All authors contributed to the article and approved the submitted version.