Minimum-norm Estimation of TMS-activated Motor Cortical Sites in Realistic Head and Brain Geometry

—Navigated transcranial magnetic stimulation (nTMS) is a widely used tool for motor cortex mapping. However, the full details of the activated cortical area during the mapping remain unknown due to the spread of the stimulating electric ﬁeld (E-ﬁeld). Computational tools, which combine the E-ﬁeld with physiological responses, have potential for revealing the activated source area. We applied the minimum-norm estimate (MNE) method in a realistic head geometry to estimate the activated cortical area in nTMS motor mappings of the leg and hand muscles. We cal-culatedthe MNE alsoin a sphericalheadgeometryto assess the effect of the head model on the MNE maps. Finally, Academy Foundation for the Vilho, the Academy of Science J. Ilmoniemi effect on the distance between the MNE CoG and the motor map (p < 0.05). The optimized coil locations were < 1 cm from the initial hotspot in 7/10 subjects. Further research is required to determine the level of anatomical detail and the optimal mapping parameters required for robust and accurate localization.


Minimum-Norm Estimation of TMS-Activated
Motor Cortical Sites in Realistic Head and Brain Geometry we determined optimized coil placements based on the MNE map maxima and compared these placements with the initial hotspot placement. The MNE maps generally agreed well with the original motor maps: in the realistic head geometry, the distance from the MNE map maximum to the motor map center of gravity (CoG) was 8.8 ± 4.6 mm in the leg motor area and 6.6 ± 2.5 mm in the hand motor area. The head model did not have a significant effect on these distances; however, it had a significant effect on the distance between the MNE CoG and the motor map (p < 0.05). The optimized coil locations were <1 cm from the initial hotspot in 7/10 subjects. Further research is required to determine the level of anatomical detail and the optimal mapping parameters required for robust and accurate localization.

I. INTRODUCTION
T RANSCRANIAL magnetic stimulation (TMS) is a non-invasive method for stimulating the cerebral cortex [1]. It is based on electromagnetic induction: a rapid current pulse in a stimulation coil produces a time-varying magnetic field outside the coil, inducing an electric field (E-field) in the cortex [2]. A large proportion of conventional TMS studies stimulated the motor cortex, since the resulting response is straightforward to detect as a motorevoked potential (MEP) in an electromyogram (EMG). These studies revealed the basic physiological principles of TMS [3], which paved the way for a variety of diagnostic [4] and therapeutic [5], [6] applications. Navigated TMS (nTMS) enables stimulation with high localization accuracy [7], [8], which is crucial for modern neurosurgery in pre-surgical mappings that aim to delineate the eloquent motor areas [9]- [13]. Conventionally, these mappings stimulate the motor cortex with a slightly supra-threshold stimulation intensity (SI), e.g., 105% or 110% of the resting motor threshold (rMT), to outline the cortical sites associated with MEPs [13]- [15]. The pre-surgical nTMS maps have been shown to correspond well with the gold standard of direct electrical stimulation (DES) based on comparisons between hotspots, defined as the cortical sites eliciting the highest This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ MEP amplitudes, or centers of gravity (CoGs) [16]- [18]. However, the supra-threshold SI in nTMS mapping increases the spatial extent of the cortical E-field sufficient to excite the cortical neurons, increasing also the extent of the resulting motor map [19]. This tends to cause overestimation of the true motor representation area [20]- [22]. Furthermore, selecting the SI without the prior consideration of individual input-output (IO) characteristics may hamper inter-individual comparisons [13], [23].
The TMS-induced E-field strength is considered a fundamental factor in defining the activated cortex and must be included in the assessment of the motor maps to gain more accurate localization and delineation of the motor representation area. The cortical E-field exhibits a non-uniform distribution, which depends on individual anatomy [24]. Several studies performed E-field modeling using a realistic head geometry and combined the individualized E-field distributions with MEP data for more accurate localization [25]- [30]. Recent studies also implemented multi-scale computational modeling, combining individualized E-field calculations with morphologically realistic cortical neurons, to take into account activation mechanisms at the neuronal level [31]- [33].
These studies focused on the hand motor area, typically concluding that the most likely activated area resides around the crown of the precentral gyrus. In the leg motor area, the homunculus extends to the interhemispheric fissure, requiring higher SI to induce a sufficiently strong E-field in the cortex. However, some leg muscles, such as the tibialis anterior (TA), have more superficial representation areas, which can be effectively stimulated with a figure-of-eight coil [34], [35]. Few studies have combined nTMS mapping with E-field modeling to localize motor function in the leg motor area [36].
Pitkänen et al. applied the minimum-norm estimate (MNE) [37] in a spherical head geometry to take into account the effect of the E-field distribution on the nTMS-mapped motor areas [38]. The MNE method, when applied to nTMS mapping, aims to solve the inverse problem by finding the motor representation area based on the nTMS-MEP trials. In the current study, we applied the MNE method in a realistic head geometry to assess the activated cortical area in nTMS motor maps of the leg and hand motor areas. We also included calculations in the spherical head model to assess the effect of head model geometry on the outcome of MNE. Furthermore, to compare the MNE map maxima with the hotspots determined during the nTMS procedure, we calculated optimized coil placements by applying the computational approach presented by Gomez et al. [39].

A. Subjects
We recruited five healthy, right-handed volunteers (two female, age: 24-40 years) for the measurements of the leg motor area. Furthermore, we analyzed our previous data [20] performed on five healthy, right-handed volunteers (two female, age: 25-48 years) in the hand motor area. The study was approved by the Research Ethics Committee of the Northern Savo Hospital District (permissions 72/2016 and 140/2017); the subjects signed a written informed consent.

B. Measurements
T1-weighted MRIs were acquired with two different 3-T scanners (Philips Achieva 3.0 T X, Philips, Eindhoven, The Netherlands or Siemens MAGNETOM Vida, Siemens Healthcare, Erlangen, Germany) to enable nTMS. In addition, T2-weighted MRIs were acquired to be used in generating the head models. TMS was applied with an nTMS system (NBS 4.3, Nexstim Plc, Helsinki, Finland) that produced biphasic pulses through a figure-of-eight coil (35-mm inner diameter, 75-mm outer diameter, 10 turns). EMG was measured with a system-integrated device and disposable Ag-AgCl electrodes from the contralateral muscles of interest (TA in the leg motor area and first dorsal interosseous (FDI) in the hand motor area) with a sampling frequency of 3 kHz. In the leg motor area, the measurements were conducted on both hemispheres, while only the left hemisphere was stimulated in the hand motor area.
To study the motor representation of the leg muscle, the nTMS mapping protocol started with an initial search for the motor hotspot by applying supra-threshold stimuli on the precentral gyrus near the longitudinal fissure. The SI was adjusted to induce MEPs with amplitudes of approximately 1 mV in the TA muscle, and the coil was kept tangentially to the scalp. The coil was oriented in the lateral direction, approximately perpendicularly to the longitudinal fissure. The initial hotspot was identified as the cortical site (estimated as the E-field maximum location by the nTMS system) that repeatedly produced the highest-amplitude MEPs, and the optimal orientation was further refined at this hotspot with the targeting tool of the nTMS system [40]. The rMT was determined (as a percentage of the maximum stimulator output, %-MSO) in this location and orientation with a systemintegrated threshold hunting algorithm [41], [42].
By applying stimuli with SI of 105% of the rMT, motor mapping was conducted with the help of a rectangular grid (5 mm × 5 mm cells) overlaid around the hotspot on the cortical view of the nTMS software. One stimulus was applied in each cell with a minimum inter-stimulus interval of 5 s, starting from the grid center and continuing to the surrounding grid cells until no MEPs were induced. The mapping protocol in the hand motor area was identical, except that two stimuli were applied in each grid cell [20].
Finally, IO curves [43] were measured by stimulating the hotspot with varying SIs at 10-%-MSO steps and by applying ten stimuli at each SI. The order of the ten-stimulus batches was randomized, and the SI range was 90-150% rMT for the TA and 90-130% rMT for the FDI muscle.

C. Motor-Evoked Potential Analysis
MEPs were detected automatically by the nTMS software and further visually inspected offline. Responses with amplitudes of at least 50 μV were accepted as MEPs, and responses with preceding muscle activity were rejected from further analysis. The CoG was calculated for each motor map by utilizing the stimulus coordinates and the corresponding MEP amplitudes [20]: where x i , y i , and z i are the coil coordinates in the motor mapping, and M i is the corresponding MEP amplitude. The cortical sites associated with the MEPs were also visualized by projecting the corresponding coil coordinates along the coil normal direction to the gray matter surface. This was done with the understanding that the accuracy of this procedure relies on the coil being very nearly tangential with respect to the local head surface [8].

D. Electric Field Modeling
The individual head models with realistic geometries were generated from the T1-and T2-weighted MRIs by using the headreco pipeline of the SimNIBS software (version 3.2) [44]. This pipeline created a tetrahedral volume mesh divided into five tissues (skin, skull, cerebrospinal fluid, gray matter, and white matter), which were associated with standard, isotropic conductivity values [24], [45].
The TMS coil dimensions were estimated from X-ray images ( Fig. 1) to form a coil model based on magnetic dipoles [46]. The coil coordinates of the motor mapping were transformed from the MRI coordinates to the SimNIBS coordinates, and SI (reported originally in %-MSO) was converted into a time derivative of the coil current. The coil model, coordinates, and current enabled the calculation of its magnetic vector potential (A), which was used as input to the finite element method (FEM) for the calculation of the electric scalar potential (V ) [24]. These potentials were then used for calculating the final E-field distribution ( Fig. 2A): The E-fields were also estimated in a spherically symmetric conductor (i.e., head model with a spherical geometry) by computing the mutual inductance between the nTMS coil and the triangle construction [47] in a grid of points on the cortical surface (Fig. 2B) [38]. This approach, including the simplified coil model, was similar to the one presented by Pitkänen et al. [38]; however, instead of a regular grid of points based on the nTMS coordinates, we formed the grid by utilizing the nodes of the realistic head mesh on the gray matter surface. For each TMS pulse, the E-field was calculated separately at the grid points. Therefore, by adjusting the calculation points based on anatomical information, our approach considered the gyrification pattern of the cortex.

E. Minimum-Norm Estimation
The cortical E-field distributions were utilized for calculating the MNE of the motor representation area (Fig. 7); the computational approach was equivalent to the method presented by Pitkänen et al. [38]. In summary, a sigmoidal function was fitted to the IO curve to transform the E-field distribution E(r) into an activating function E (r), and the motor representation area excitabilities x were solved from where m is a vector containing the MEP amplitudes obtained from different stimulation locations, the i th row of E is the activating function used to evoke the i th MEP in m, and x is a vector containing the (unknown) excitabilities in the locations where the E-fields were estimated. These excitabilities can be estimated by the MNE of x: where E + is the pseudoinverse of E , calculated by applying the singular value decomposition. The noise is taken into account by regularization; similarly to the original approach of Pitkänen et al., the present study applied the Wiener regularization, where E and m are orthonormalized [38]. The resulting MNE map consists of the values ofx in the cortical locations where the E-field was estimated.

F. Coil Placement Optimization
We estimated the optimized coil placement (i.e., location and orientation) for both head models (i.e., with realistic and spherical geometries) and compared it with the coil placement corresponding to the initial hotspot determined in the mapping procedure. For comparison, we selected the maximum of the MNE map in the precentral gyrus separately in both head models as an estimate of the activated cortical area. We then determined the optimal coil location and orientation, which maximized the E-field strength at this location, by applying the optimization routine developed by Gomez et al. [39]. This routine compared 13896 different coil placements in an area of 20-mm radius around the initial target with the following parameters: spatial resolution = 2.5 mm, angle resolution = 5.0 • , search angle = 360 • .

G. Statistical Methods
The distances from the motor hotspot and CoG to the MNE map maximum and CoG were calculated and compared between the head models. Furthermore, the distances from the hotspot coil location to the optimized coil locations were compared between the head models. Finally, the coil orientations were compared between the hotspot coil orientation in the initial mapping and the optimized orientations. All these comparisons applied the Wilcoxon signed-rank test, and the threshold of statistical significance was p = 0.05. MATLAB (version: 2015b, MathWorks Inc., Natick, MA) was used for statistical analyses.

A. Motor Maps, Hotspots, and Centers of Gravity
The motor mapping protocol was successfully completed for all the participants. Fig. 3 shows the motor maps of the TA muscle together with the cortical projections of hotspot and CoG. In the leg motor area, the total number of stimuli in the motor map was 24 ± 14 (range = 7-59). The average distance between the TA hotspot and CoG was 6.1 mm (95% confidence interval (CI) = 2.8-9.3 mm). Fig. 4 shows the motor maps of the FDI muscle and the corresponding hotspots and CoGs. In the hand motor area, the total number of stimuli was 55 ± 41 (range = 20-101). The average  distance between the FDI hotspot and CoG was 7.2 mm (95% CI = 0.9-13.5 mm). Fig. 8 illustrates the E-field distributions corresponding to hotspot stimulation. These maps were thresholded at 50% of their maximum for better visualization of the MNE maxima. Fig. 9 illustrates the effect of thresholding on the MNE maps. The MNE map maxima and CoGs were close to the motor map CoGs and hotspots (Table I). When comparing the distances in different head models (Table I), the head model had a statistically significant effect only on the distance from the motor map to the MNE CoG in the leg motor area ( p = 0.020). Table II Fig. 5. Minimum-norm estimations of the leg (tibialis anterior muscle, S1-S5) and hand (first dorsal interosseous muscle, S6-S10) motor representations in realistic (red) and spherical (blue) head geometries. The maps are thresholded at 50% of their maxima. The maxima are marked with crosses.

C. Optimized Coil Placements
The optimized coil placements and hotspot coil locations are shown in Fig 6. In general, the optimized coil locations were in  good agreement with the hotspot coil locations (Table III): the mean distance from the hotspot coil location to the optimized coil location was between 7.1 and 15.6 mm. The difference between the distances in the different head models was not statistically significant. In three subjects (S1, S5, and S9), the optimized coil location was >1 cm away from the initial hotspot coil location. Table IV lists the optimized coil orientations. The range of the mean orientations was 19-76 degrees in the leg motor area (with respect to the longitudinal fissure) and 51-70 degrees in the hand motor area (with respect to the direction perpendicular to the longitudinal fissure). In the leg motor area, the coil orientations optimized for the spherical head geometry were systematically different from the orientations optimized for the realistic geometry ( p = 0.001) and the hotspot coil orientation ( p = 0.002).

IV. DISCUSSION
We applied MNE in head models with realistic and spherical head geometries to assess nTMS motor maps. The MNE maps of the TA and FDI muscles agreed with the original motor maps based on their distance to the motor map CoGs and hotspots. The head model had a significant effect on this distance only for the MNE CoG in the leg motor area. Furthermore, we determined optimized coil placements for the MNE maxima in the different head models; the optimized coil locations were generally close to the initial hotspot coil location.

TABLE IV HOTSPOT COIL ORIENTATION AND OPTIMIZED COIL ORIENTATIONS
In the leg motor area, the hotspots of the TA muscle were close to the interhemispheric fissure, which is in agreement with previous studies [35], [48]. In the hand motor area, the FDI hotspots were located in the "hand knob" [49]. This implies that the initial mapping was accurate in localizing the muscle-specific hotspot, which is crucial in the rMT determination and the subsequent motor mapping [13]. In most subjects, the motor maps were elongated in the direction of the E-field (Figs. 3 and 4), which is a characteristic finding in nTMS motor maps [50].
Previous studies have demonstrated that nTMS motor mapping is repeatable in terms of hotspots and CoGs between sessions [51], [52]. The present study indicated that the methods for estimating the source of the TMS-induced muscle activation may reveal independent characteristics, explaining the differences observed in hotspot and CoG locations, as well as in the optimized coil location and the initial hotspot. Therefore, the consistency and accuracy of these methods in nTMS motor mapping needs to be further studied to estimate their feasibility.
The MNE values were consistently highest at the crowns of the precentral gyrus. This is because the MNE solution, as its name indicates, provides a solution with the minimum L2 norm, which is best accomplished with superficial solutions. The MNE maps had a spatially limited extent, supporting the existence of a common activation site during motor mapping, which has been suggested by earlier modeling studies [20], [21], [27], [30]. In both the leg and hand motor areas, the MNE map maxima and CoGs were close to the motor map hotspots and CoGs (Table I): the range of the mean distances was 3-12 mm. Interestingly and importantly, these distances are similar to the reported distances between the nTMS hotspots and DES hotspots [16], [17], [53], [54].
The group average coordinates (Table II) of the TA hotspot matched with previously reported coordinates in the MNI space [35]; however, the group average MNE maximum in the realistic head geometry was more laterally located in the right hemisphere, where subjects S2 and S5 had more lateral MNE maxima when compared with the left hemisphere (Fig. 5). Importantly, the group average coordinates of the FDI MNE maximum agreed well with a recent modeling study, which regressed the E-fields against the elicited MEPs in each cortical element to determine the effectively stimulated cortical site [55].
Most previous studies that have combined the TMS-induced cortical E-fields with MEP measurements have focused on the hand motor area. These studies have mainly suggested the E-field strength at the gyral crowns as the primary predictor of activation [25], [26], [30], while some studies have highlighted the E-field normal at the sulcal walls [27], [56], [57]. The MNE method prefers the gyral crown as the origin of activation, since the MNE maps attenuate strongly when moving deeper in the sulcus.
The methodological approach of the present study differs significantly from those presented in previous modeling studies that aimed to outline the activated cortex. Opitz et al. [28], [29] mapped the motor cortex at 25 different locations with a fixed coil orientation and calculated the E-field CoG (weighted by MEP amplitude) as an estimate of the excitation area. This estimate exhibited relatively large extent around the "hand knob", requiring strong thresholding to match the DES-based area [29]. Aonuma et al. [25] multiplied the E-field strengths corresponding to the five stimulations associated with the highest MEP amplitudes. This estimate, after strong thresholding, was in good agreement with DES. However, these approaches, which were based on additive or multiplicative processing of the E-fields, may cause biases in the resulting estimate towards the areas with consistently higher E-field strengths during the mapping [30].
Another approach, which is probably not as susceptible to the above-mentioned bias, is based on the relative standard deviation (RSD) of the product between the MT and E-field strength for all the coil positions/orientations in the mapping. As this estimate requires the measurement of the MT, relatively few coil positions and orientations were measured in the studies applying this approach [26], [27]. While this may hamper the discriminative power at the individual level, the group-level activations were located focally at the "hand knob". The current study resembles the approach of Weise et al. [30], and especially the method applied by Numssen et al. [55], who also performed the mapping in a high number of different coil positions and considered the IO properties of the cortex in the calculation. However, Fig. 7. The minimum-norm estimation method. A: A sigmoidal function is fitted to the input-output curve, which is measured at the hotspot (marked with an asterisk in the original cortical mapping). M is the response maximum, erf is the error function, I is the stimulation intensity, and μ and σ define the midpoint and transition width of the sigmoidal function, respectively. B: The original motor mapping consists of i stimulations and motor-evoked potential amplitudes (the elements of vector m). These stimulations are modelled to obtain the corresponding E-field strength values in j cortical elements. The E-field values E(r) are transformed with the previously determined sigmoidal function to get an "activating function" E'(r). C: The excitability values x can be estimated by applying the pseudoinverse of E. their method was based on a different mathematical method: the goodness-of-fit of the IO curve at the different cortical elements. The MNE-based method presented in the current study provides a mathematically intuitive inverse approach for solving the source (excitability) distribution in TMS motor mapping. Direct validation is a challenge in modeling studies that combine E-field distributions with MEP measurements to localize the activated cortex. Valuable validation data can be obtained with DES [25], [29], [36], which is not an option for healthy volunteers. In the present study, we cannot provide direct validation of the MNE method, but provide data on Fig. 9. Minimum-norm estimations of the leg (tibialis anterior muscle, S1-S5) and hand (first dorsal interosseous muscle, S6-S10) motor representations in realistic (red) and spherical (blue) head geometries. The maps are thresholded at 50%, 70%, and 90% of their maxima. The maxima are marked with crosses. feasibility by comparing different methods. The feasibility of the MNE maps was evaluated in a manner that resembles that used by Numssen et al. [55]; they determined the optimal coil placement for stimulating the suggested activation site, and measured the rMTs for the optimal and adjacent placements. In the current study, we determined the optimized coil placement and compared it with the initial hotspot coil placement. The optimized coil location was close to the hotspot coil location in most subjects (Fig. 6), which suggests that our MNE approach located the muscle representation area effectively. However, subjects S1 (right hemisphere), S5, and S9 appeared as outliers. Their hotspot stimulation (Fig. 8) induced a high E-field strength in a neighboring perpendicularly oriented gyrus instead of the underlying gyrus, which was most strongly stimulated in the other subjects. Since the coil orientation remained relatively stable during the mapping, the neighboring gyrus probably accumulated the highest E-field values consistently throughout the measurement, resulting also in the highest MNE value. In these subjects, the larger distance from the optimized coil location to the hotspot coil location implies that MNE map may misrepresent the activated area, especially in the realistic head geometry. The coil orientations that were optimized based on the realistic head geometry matched relatively well with the initial hotspot orientations. In the leg motor area, the orientations deviated from the standard lateral orientation, which is consistent with the findings of Richter et al. [58].
Increasing the range of coil orientations during the motor mapping may improve localization accuracy. Weise et al. [30] noticed that higher variability between E-fields induced by the stimuli of the mapping protocol sharpened the localization results in their approach, which was based on the correlation between MEPs and the induced E-field. Moreover, sets of experimental conditions with selectively varied coil positions or orientations, especially the 45 • orientation in the hand M1, did not contain sufficient information for unique localization. Therefore, our approach, albeit based on a different computational method, could benefit from a mapping protocol utilizing more varied coil positions and orientations; reducing the correlation of the E-fields across these coil placements could reduce the number of stimuli required for accurate localization. The effect of mapping protocol on the resulting MNE map should be investigated in future studies.
The main limitations of the current study are the low number of subjects and the lack of proper validation data. The number of subjects included in the present study was justified by the study design, in which the intra-individual variations were of major interest, while the inter-individual variations provided an estimate of the general level of potential systematic differences in the evaluated parameters. In addition, the IO curves were measured in the hotspot to obtain an activating function, which was assumed to be the same for all cortical locations in the MNE calculation (Fig. 7). This assumption is valid based on the measurements of Thickbroom et al. [59], who showed that the shape of the IO curve is very similar at different cortical sites. While these limitations hinder generalized conclusions and induce relatively wide confidence intervals in the results, the MNE method appears a promising tool for locating the TMS-activated cortical sites as it considers the spread of the E-field. The MNE, computed with very fast linear matrix multiplications, may also prove useful in real-time adjustments of TMS targeting in closed-loop multi-locus TMS [60]. However, limitations still include the lack of proper quantification of representation area size based on MNE. From the perspective of pre-surgical mapping, the motor maps obtained with nTMS exhibit higher accuracy than other non-invasive methods, such as functional MRI, when compared with the gold standard of DES [18]. The MNE maps have potential in narrowing down the delineated cortical area,  Tables V-VII  list the individual data corresponding to Tables I, III, and IV,  respectively.