Fine multi-coil electronic control of transcranial magnetic stimulation: effects of stimulus orientation and intensity

Despite the evident benefit of transcranial magnetic stimulation (TMS) to probe human brain function and treat neurological diseases, current technology allows only a slow, mechanical adjustment of the electric filed orientation. Automated and fast control of the TMS orientation is critical to enable synchronizing the stimulation with the ongoing brain activity. We overcome these limitations with a two-coil electronically controlled TMS transducer to define precisely the pulse orientation (∼1° steps) without mechanical movement. We validated the technology by determining the dependency of motor evoked responses on the stimulus orientation and intensity with high angular resolution. The motor response was found to follow a logistic function of the stimulus orientation, which helps to disentangle the TMS neuronal effects. The electronic control of the TMS electric field is a decisive step towards automated brain stimulation protocols with enhanced accuracy of stimulus targeting and timing for better diagnostics and improved clinical efficacy.


Introduction
Transcranial magnetic stimulation (TMS) has been used extensively to study human brain function and to treat a multitude of neurological diseases 1 . In TMS, the intensity and orientation of the induced electric field (E-field) are critical factors to excite specific neuronal populations [2][3][4] and to activate distinct neurotransmitters 5,6 . Conventionally, the E-field orientation for optimal response is adjusted manuallya slow process required in every TMS application. Manual TMS control does not allow the manipulation of the E-field orientation at the millisecond-level time scale of cortical signalling. Overcoming this limitation with electronic control of the E-field orientation will enable finely defined modulation of neuronal populations, a methodology that may prove effective in enhancing clinical efficacy.
The induced E-field parallel to the longitudinal axis of pyramidal neurons, i.e., along the cortical columns, leads to maximal excitation of the targeted neurons, following qualitatively the so-called cortical column cosine model 7 . On a macroscopic scale, the cosine model implies that the stimulation effect is proportional to the cosine of the angle between the E-field and the normal of the cortical surface. Such a model has been largely employed in simulations to estimate the effect of TMS in the brain 8,9 . On the other hand, recent computational and experimental findings suggest that TMS excites neurons mostly at axon terminals and that all components of the E-field, i.e., not only its normal component, have a substantial contribution to the neuronal depolarization [10][11][12] . Although previous results have advanced our understanding of TMS-triggered cortical activation, empirical models that fully describe the relationship between the evoked response and the stimulus orientation and intensity are still missing.
We developed a multi-coil TMS (mTMS) transducer that allows an experimenter to rotate the peak Efield electronically, i.e., without manual coil movement. The mTMS concept was introduced by Ruohonen et al. 13,14 and recently elaborated by  to electronically shift the peak Efield along a line 16 . Our new transducer allows delivering consecutive pulses in different orientations at millisecond intervals, enabling an unprecedented non-invasive approach to probe distinct neuronal mechanisms. Furthermore, we determined quantitatively how the evoked motor responses depend on the stimulus orientation and intensity, emphasizing the importance of fine E-field adjustments to unveil the neuronal effects of TMS.

mTMS transducer to rotate the E-field electronically
Our mTMS transducer has two overlapping perpendicular figure-of-eight coils to control the induced Efield orientation electronically. In a single TMS pulse, we manipulate the net E-field orientation by driving a monophasic current waveform simultaneously through each coil with pre-calibrated capacitor voltages. The E-field intensity decreases rapidly with the distance from the cortex 17 ; thus, we aimed at a thin and compact transducer case (150-mm outer diameter) to improve the coil-cortex coupling and reached up to 129 V/m measured with our probe at 70 mm from the centre of an 85-mm-radius spherical head model.
The two sets of coil windings were computed using a minimum-energy optimization procedure developed by   (Online Methods). This procedure computes the minimum-energy surface current density distributions that reproduce a reference E-field distribution generated by a widely used commercial figure-of-eight coil (Magstim 70mm Double Coil; The Magstim Co Ltd, UK; Fig. 1a-b) 18 . The discretized surface current densities provide the winding paths for manufacturing the coils (Fig. 1c). We wound two layers of 12 turns of litz wire (1.6 mm thick and 2.4 mm wide; Rudolf Pack GmbH & Co. KG, Germany) to two 3D-printed formers (Fig. 1d). The formers were printed by selective laser sintering of glass-filled polyamide (Maker 3D, Finland) which has a tensile strength of 38 MPa and dielectric strength of 15 kV/mm, and are resistant to the pressure of about 10 MPa from the Lorentz forces during the TMS pulse (derived from Ilmoniemi et al. (1999)). The coils were potted with epoxy for further mechanical support and electric insulation 20 .
The full width at half maximum (FWHM) of the E-field profiles parallel and perpendicular ( Fig. 1e-f) to the peak E-field below the coil centre were computed to verify the stimulation focality. The FWHM in the direction perpendicular to the peak E-field of the bottom and top coils were 25.6 and 26.5 mm, respectively. The FWHM in the direction parallel to the peak E-field was 44.8 and 46.4 mm for the bottom and top coils, respectively. The capacitor voltage required by each coil follows the cosine and sine function of the desired E-field orientation and norm, as expected, given that the coils are orthogonal. The bottom and top coils required 868 V and 1128 V of capacitor voltage, respectively, to induce a 100 V/m E-field measured by our probe 21 (Online methods). The average deviations from the theoretical values in the E-field norm for stimulation at 25 V/m and orientation from 0° to 180°, were 0.1 V/m and 1.3°, respectively. The measured E-field distributions with the stimulus orientation electronically rotated to 0°, 45°, and 90° are illustrated in Fig. 1g-i. Figure 1: mTMS transducer and E-field distributions. a, Example of a target E-field induced on a spherical surface by a model of a commercial figure-of-eight coil. b, Optimized minimum-energy surface current density distributions for the bottom and top coils. c, Surface current distributions of (b) discretized in 12 turns. The induced E-field has a similar focality to that of the figure-of-eight coil in (a), but the induced E-field can be rotated by adjusting the currents in the top (orange) and bottom (black) coils. d, mTMS transducer with litz wire wound in the bottom and top 3D-printed formers. e-f, The profile along the direction parallel and perpendicular to the peak induced E-field, respectively. The solid and dashed lines refer to the bottom and top coils, respectively. g-i, Measured E-field distribution induced on a spherical surface with a 70-mm radius for a TMS pulse at 0°, 45°, and 90°. The shaded grey outer boundary represents a spherical scalp with a radius of 85 mm on top of which the coil assembly was placed.

Effect of E-field orientation and intensity on motor response
The 2-coil mTMS transducer allowed us to model the effect of the E-field orientation and intensity on the motor evoked potential (MEP) amplitude and latency with high angular resolution. We collected MEPs on the abductor pollicis brevis (APB) muscle on 16 healthy subjects at every 3° of stimulus orientation (coordinate system illustrated in Fig. 2a). The neuronal excitation properties critically depend on the stimulation strength 22 , which potentially affects the orientation dependency curve. To investigate this effect, we tested four stimulation intensities: 110, 120, and 140% of the resting motor threshold at 0° (MT0°), and 120% of the resting motor threshold at 90° (MT90°).
Changing the stimulus orientation for the excitation of a neuronal population, for instance pyramidal neurons with anisotropic dendritic arborization, is perceived fundamentally as a variation in the effective stimulation intensity 23 . The logistic (sigmoid) equation describes well the relationship between the stimulus strength (input) and the neuronal output 22 , making it a suitable candidate to parametrize the orientation dependency of the cortico-motor response. We visually assessed how other conventional models, such as the cosine 7 and the Gaussian equation 24 , fit the measured MEP amplitude versus the stimulus orientation. These models failed to fully capture the essential characteristics of the entire peak, as shown in Fig. 3a. Note that the cosine does not fit well to the exponential increase in the MEP amplitude between 45° and 90°; the Gaussian model does not match the inflection around ±45° and the slightly broad plateau around the peak. Therefore, we modelled the log-transformed MEP amplitude versus the stimulus orientation with the generalized logistic equation: where y = log 10 MEP , the MEP is the MEP amplitude, b is the baseline MEP amplitude, max is the maximum MEP amplitude, θ is the stimulus orientation, and θ 0 and σ are the center and slope of the logistic equation, respectively. These four parameters can provide information on the structuralfunctional relationship of the primary motor cortex. Given a set of pre-defined physical constrains, such as the pulse shape, coil-to-cortex distance, and intensity, max represents the neuronal capacity as the maximal response that can be obtained from the stimulated population and b the corresponding minimum neuronal response. The σ represents the orientation sensitivity of the neuronal population, high values indicate a neuronal population preferentially aligned in a narrow range of orientations resulting in a sharp change between the baseline and maximum amplitudes. Lastly, θ 0 characterizes the intrinsic properties of the neuronal population combining max , b , and σ.
The orientation dependency of the MEP amplitude was greatly affected by the increase in the stimulation intensity. We observed the same effect in all five subjects, illustrated with data from subject 16 (S16) in Fig. 2b-d. The traces in Fig. 2b show two clusters of motor responses around 0° and 180° with increasing response amplitude for the stimulus intensities 110%, 120%, and 140% of the MT0°. At  Fig. 3b). The maximum MEP amplitude reached the saturation value at a lower intensity (110% of the MT0°) than the baseline amplitude which only increased with the tested stimulation intensities (Fig. 3b-c). This was expected due to the 28% higher MT at 90° compared with 0° orientation (one-tailed paired t-test, p < 0.001). Thus, the difference between the 0° and the baseline amplitude, at ±90°, greatly reduced with the increase in the stimulation intensity ( Fig. 3b-c). The highest stimulation intensities also led to similar MEP amplitudes across all orientations, evidenced by the change from an eight-shaped curve at 110% MT0° to a circle with constant amplitude at 120% of MT90° (Fig. 2c). The centre orientation increased linearly with the intensity (Fig. 3d), which corresponds to the broader peaks in the orientation dependency curve, as illustrated for S16 (Figs. 2c). In turn, the logistic equation slope (sigma) did not show a clear variation with the stimulation intensity apart from a small difference between the 0° and 180° peaks at 120% of MT90° (df = 34.0; tratio = −3.27; p < 0.039; Fig. 3e). Figure 2: MEP amplitude and latency as functions of the stimulus orientation and intensity. a, A schematic representation of the stimulus orientations relative to the subject's head. The 0° orientation refers to the first phase of the induced E-field pointing to the anteromedial orientation, approximately perpendicular to the central sulcus. b, The MEP epochs (median across three trials), from 20 to 50 ms after the TMS pulse, at each stimulus orientation and intensity. The colormap encodes the stimulus orientation depicted in a. c,d, MEP amplitude (radial axis in logarithmic scale) and latency as functions of the stimulus orientation. The dots represent the median MEP amplitude and latency across the repeated trials for each stimulus orientation and intensity. The solid lines represent the logistic and the trigonometric polynomial fits of the (c) MEP amplitude and (d) latency, respectively. The colormap indicates the stimulation intensity in % of the resting motor threshold at 0° (MT0°) and 90° (MT90°). The data in (b-d) was recorded from subject 16 (S16).
The MEP latency is commonly used as a proxy to the mechanisms involved in the action potential generation 25 . A shorter latency may indicate direct depolarization of the neuron's axon, whereas a longer latency is associated with transsynaptic or indirect neuronal excitations. We fit a 2 nd -order trigonometric polynomial to the MEP latency for each subject, which was also affected by the TMS orientation and intensity (Figs. 2d and 3f). For all stimulation intensities, MEP latencies at 0° and 180° were about 2 ms shorter compared to those at 90° and 270°. Moreover, the lowest stimulation intensity The grey markers are the median MEP amplitudes across three trials recorded from one subject (S15) with stimulation at 110% MT0°. b-e, Linear mixed model analysis of the logistic equation parameters maximum amplitude ( max ), baseline amplitude ( b ), centre orientation (θ 0 ), and slope sigma (σ) across subjects. The error bars represent the 95% confidence interval of the estimated marginal means. The difference in the error bar sizes is due to the distinct number of subjects in each intensity (16 subjects in 110% MT0°, 2 subjects in 120% MT0°, and 5 subjects in 140% MT0° and 120% MT90°). f, Linear mixed model analysis of the MEP latency extracted from the four peak orientations (−90°, 0°, 90°, and 180°) in the trigonometric polynomial fit. We removed the MEP latency from 120% MT0° because only subjects showing latencies across the entire range of orientations for each intensity were included in this analysis (7 subjects in 110% MT0°, and 5 subjects in 140% MT0° and 120% MT90°).

Discussion
The developed electronic control of E-field orientation opens the possibility for exploring new cortical stimulation paradigms and probing neuronal mechanisms. The traditional manual coil placement is timeconsuming and prone to localization errors on the order of several millimetres and degrees in orientation [26][27][28] , being a major source of variability across sessions 27 . Our millisecond-control of stimulus orientation enables the development of automated algorithms 29 for closed-loop paradigms triggered by neurophysiological recordings, increasing the efficacy for TMS in both research and clinical applications 30,31 . It also enables novel studies on intracortical inhibition and facilitation mechanisms 5,6 ; for instance, in paired-pulse protocols changing the E-field orientation within a millisecond interval without the need for the mechanical movement of the transducer. Moreover, we demonstrated with high angular resolution that the generalized logistic equation can be a simple model to parametrize the motor response dependency on the stimulus orientation and intensity, benefiting future realistic computational models to explain the TMS neuronal effects.

mTMS transducer
The developed energy-efficiency-optimized mTMS transducer allows one to change the induced E-field orientation by adjusting the ratio of the currents applied to the two overlapping coils. The transducer diameter is similar to conventional TMS coils, and only half the size (150 mm) of previously manufactured minimum-energy coils 15,32 . The transducer was built with a compact design for more comfortable positioning over the scalp and better handling. The reduced size requires a higher voltage (bottom coil, 868 V; top coil, 1128 V) to generate a desired pulse than a previously built large minimumenergy figure-of-eight coil (570 V) 32 . The required energy depends on the stimulation orientation, with the bottom coil requiring 45 J, which is 12% higher than the energy required by a previous model (40 J) 32 and 41% lower than the energy required by the top coil (76 J). Furthermore, the transducer's bottom coil required 23% lower capacitor voltage than the top coil to generate a 100 V/m E-field at the depth of 15 mm from the surface of an 85-mm radius sphere. This difference is due to the top coil's 5-mm extra distance from the cortical surface, which leads to a weaker coil-cortex coupling 32 . The larger distance also leads to 3% (1-2 mm) poorer focality in perpendicular and parallel directions for the top compared to the bottom coil, a negligible difference for standard TMS applications. If desired, this difference can be eliminated in future models by designing the bottom coil so that its focality matches the other coil's pattern. The possibility to use 3D-printed formers might ease the manufacturing of transducers with more complicated winding patterns, such as those reported by Deng et al. 33 and future mTMS transducers with more coils 15 .
The manufactured mTMS transducer requires lower energy and current to generate a desired E-field than a conventional TMS coil. To induce a 60 µs 100 V/m pulse, the developed bottom coil requires 2.6 kA and 45 J whereas a traditional TMS coil (Magstim 70mm Double Coil) requires 3.8 kA and 120 J 34 . The energy and maximum current were estimated using the surface current density distribution, as reported by Koponen et al. 32 . It is worth noting that the E-field values reported in this study were verified by measurements with a triangular probe 21 that accurately mimics E-field recordings at 70-mm radius in the spherical head geometry. The induced E-field distribution in a real brain is influenced also by cortical folding, conductivity profiles, and other deviations of the head from spherical symmetry 3,35 .

Orientation selectivity of the TMS effect
Our measurements revealed that the MEP amplitude abruptly decreases when moving away from the optimal orientation, with a stronger dependency than just the first-degree cosine of the angle. For instance, on average for all the subjects, at 110% of MT0°, a 46.2° deviation from the 0° orientation led to a reduction by half relative to the maximum amplitude. In addition, the peak amplitude at 0° was almost twice the amplitude at 180°. This strong asymmetry has been recorded at the single-cell level and is possibly related to the excitation of distinct neuronal populations 36 . Moreover, we observed a single, smooth amplitude peak at the optimal orientation for all the 16 participants. We did not observe any bimodal response curves at the 0° orientation as reported for most participants by Kallioniemi et al. (2015). The bimodal response might have been caused by inaccuracies inherent to the manual coil adjustments which are eliminated by the electronic rotation of the E-field with the 2-coil mTMS.
The increase in stimulation intensity severely reduced the orientation selectivity of the TMS effect. One likely explanation is that distinct neuronal populations are recruited depending on the stimulus orientation 25,36 . For instance, if we assume that on the hand knob area of the primary motor cortex a neuronal population A is optimally aligned with the 0° stimulus while a population B, slightly further away from the hotspot, is aligned at 90°, illustrated in Fig. 4a; a low-intensity stimulation would selectively excite the population A. An increase in the intensity would rapidly lead to the maximal MEP amplitude of population A (saturation level in the input-output curve) while only slightly above the excitation threshold of population B. Further increasing the stimulation intensity would excite more the population B without any further effect on A, which is already at maximum capacity. Our findings support this activation model from multiple perspectives. First, the MEP amplitude induced with stimuli at ±90° rapidly increased to a similar level than those obtained at 0°, as evidenced from the steeper increase in the baseline amplitude compared to the maximum amplitude (Fig. 3b-c). Second, the centre orientation increased linearly with the intensity, indicating a reduction in the orientation selectivity (Fig.  3d). Lastly, the MEP latency was about 1.5-2.6 ms shorter at the 0° compared with the ±90° orientations (Fig. 3f). Such delay would justify the excitation of a neuronal population B in the vicinity of the cortical hotspot, projecting to the hotspot through one synaptic connection. In fact, a recent study we showed that neuronal populations in the vicinity of the cortical target site seem to have an effect on the generated neural drive 16 . The induced E-field along the somatodendritic axis gives rise to an extracellular current on an ensemble of neurons firing synchronously. When seen from a large distance, the phenomenon can be approximated with an equivalent dipole. The illustration objects are not scaled to each other.
The logistic function provides a simple mathematical formulation based on neurophysiological principles 37,38 to describe the neuronal orientation selectivity to the TMS-induced E-field, encompassing all cortical mechanisms engaged in the MEP generation. A neuronal population firing action potentials synchronously through parallel dendrites can be approximated by an equivalent electric current dipole. The obtained sensitivity profile would then indicate that the TMS pulse induces extracellular currents along with the somatodentritic axis of neuronal populations that are tangential to the scalp surface in the precentral sulcus wall, illustrated in Fig. 4b, reciprocal to the principles in magnetoencephalography 39 . Also, the probability of neurons firing increases exponentially as the stimulation reaches the membrane depolarization threshold. Second, as the stimulus approaches the optimal orientation, the neuronal firing rate reaches a maximum level defined by the membrane hyperpolarization. These transitions are smooth due to the specific oscillations in the resting membrane potential of each neuron in the ensemble. The MEP results from combined cortical and subcortical (cerebellar and spinal) modulation 4,40 . It may be an oversimplification to assume that the neuronal or muscle activation is proportional only to the cosine between the induced E-field and the pyramidal neurons' somatodendritic axis, as previously reported 7,8 . This view is supported by recent multi-scale realistic simulations, in which all E-field components seem to contribute significantly to the neuronal excitation due to the widespread axonal ramification in different orientations, with the normal component leading to a 30% higher contribution 11 . Furthermore, the excitatory and inhibitory interneurons contributing to the cortical excitation and its control exhibit specific alignments relative to the cortical columns 41,42 . For instance, excitatory neurons mainly project from layers 2/3 to pyramidal neurons in layer 5 while inhibitory interneurons exhibit a stellate arborization with mainly horizontal projections within layer 2/3 43 . In this case, inhibitory mechanisms have a lower excitation threshold and are less affected by the E-field orientation 5 , possibly explaining the stronger decay in the muscle response on suboptimal directions. Thus, our models provide evidence on the structure-function relationship of neuronal populations from the primary motor cortex and the TMS mechanisms of neuronal excitation.
We found that MEP latencies are shortest for the 0° stimulus, i.e., when the E-field is perpendicular to the cortex, and somewhat longer for the ±90° stimulus for all stimulation intensities. Changes in latency presumably reflect differences in the generation of cortical direct and indirect epidural descending volleys. In contrast to our observations, previous studies showed that monophasic TMS pulses delivered perpendicular to the central sulcus elicited stronger indirect waves (I-waves), i.e., long-latency MEPs, whereas pulses along the central sulcus have been shown to elicit stronger direct (D-) waves resulting in short-latency MEPs 25,44 . One possible explanation is that the optimal orientation (around 0°) evoked earlier descending volleys I1-I3 (Day at al., 1989), while 90° pulses preferentially recruited later waves, e.g., I3, or neuronal populations in the vicinity of the hotspot. Thus, our findings suggest that axonal activation occurred preferentially with E-field aligned parallel to neuronal bundles at bending terminals, supported by early simulation studies 23 . Besides, we observed slightly longer MEP latencies for 180° than the optimal 0° only at the lowest stimulation intensity, as previously detected 45,46 . The selective activation of distinct circuitry and synaptic inputs in each specific E-field direction likely explains the difference in latency between opposing directions 47,48 . This is in line with the observation that an increase in stimulation intensity seems to reduce the neuronal selectivity of the TMS pulse.
We should note that the relation between MEP and stimulus orientation and intensity might differ depending on several other factors, such as the pulse waveform 22,47 and studied muscle 2,49 . The neurophysiological implications from our results should be carefully interpreted as the estimated E-field orientation is based on the spherical head geometry used to design and calibrate the mTMS transducer, not on a realistic morphology. However, due to the symmetric and smooth change in MEP amplitude relative to the stimulus orientation (Fig. 2c), it is unlikely that spurious changes in the E-field maximum 8,10 in the cortex would affect the recorded MEPs.

mTMS transducer
We extracted the winding paths for the construction of the two-coil mTMS transducer from the surface current density distributions of a minimum-energy optimization solution 15 . The coil windings were optimized to reproduce reference E-field distributions induced by a widely used commercial figure-ofeight coil (Magstim 70mm Double Coil; Magstim Co Ltd, UK; Fig. 1a) 18 rotated to multiple orientation over a spherical head geometry (70-mm cortical radius; 2562-vertex triangular mesh). We placed this coil model 15 mm above the spherical cortex and rotated it from 0 to 180° in steps of 10° around the normal of the cortex and the coil face and calculated the E-field distribution for each case. Next, for each reference E-field, we calculated the surface current density in an octagonal plane (15-cm outer diameter; 1089-vertex triangular mesh) that produces a minimum-energy magnetic field inducing an Efield distribution that has the same focality and peak intensity as the reference E-field. The resulting set of surface current densities was decomposed with singular value decomposition with each component representing a single coil. The optimization was performed separately with two coil planes placed 15 and 20 mm above the spherical cortical surface, respectively, as the manufactured coils could not occupy the same physical space. The first two components were selected to construct two physical coils because the peak E-field rotation around the normal direction of the cortex is a problem with one degree of freedom (Fig. 1b), the first and second components corresponded to the bottom and top coils, respectively. Finally, we computed the contour lines of the surface current distribution stream function of these two components to obtain the paths for the coil windings (Fig. 1c).
The coil formers were designed in SolidWorks 2016 (Dassault Systèmes SA, France). The formers were 5 mm thick (1-mm solid bottom and 4-mm-deep grooves) printed by selective laser sintering of 30% glass-filled polyamide (Maker 3D, Finland). Glass-filled polyamide was selected due to its high stiffness (tensile modulus of 6 GPa), good mechanical strength (tensile strength of 38 MPa), and suitable electrical insulation (dielectric strength of 15 kV/mm). Each coil comprised of two layers of copper litz wire (1.6 mm thick and 2.4 mm wide; Rudolf Pack GmbH & Co. KG, Germany) in series and glued to the former grooves with epoxy. The remaining space inside each coil former was filled with epoxy for increased strength. Each coil was soldered to a cable provided by Nexstim Plc (Finland).
The finalized transducer was calibrated to allow accurate electronic control of the peak E-field orientation by adjusting the voltages applied to each coil. The E-field distribution was measured with a probe having two perpendicular triangular coils with 5-mm tangential and 70-mm radial lengths 21 . The spatial distribution of the E-field was measured at 1000 points in a hemisphere for stimulus orientations of 0°, 45°, and 90°. The self-inductance and resistance of each coil were determined by connecting the coil in series with a half-bridge circuit with a 99.3 Ω resistor and measuring the phase difference between the voltage across the coil and the half-bridge. The circuit was powered with a sinewave voltage at frequencies of 1, 2, 3, 5, 10, 20, and 50 kHz (AFG1062, Tektronix, USA); voltages were measured with an oscilloscope (InfiniiVision MSOX3034T, Keysight, USA

Participants
Sixteen healthy subjects (mean age: 29 years, range 22-41; five women) participated in the study. All participants gave written informed consent before their participation. The study was performed in accordance with the Declaration of Helsinki and approved by the Coordinating Ethics Committee of the Hospital District of Helsinki and Uusimaa.

Experimental procedure
Subjects, sitting in a reclining chair, were instructed to stay fully relaxed during the TMS sessions. EMG was recorded from the right APB, abductor digiti minimi, and first dorsal interosseous muscles using surface electrodes in a belly-tendon montage (Neuroline 720, Model 72001-K/12; Ambu, Denmark). Data from the abductor digiti minimi and first dorsal interosseous were used to monitor that hand muscles were relaxed during TMS and were not considered for further analysis. The skin surface (epidermis) was scrapped with sandpaper and cleaned with alcohol to reduce the electrode-skin impedance. EMG was amplified and recorded with a Nexstim eXimia EMG device (500 Hz low-pass filtering, 3000 Hz sampling frequency; Nexstim Plc). The mTMS transducer was placed relative to the subject's cortical anatomy using a neuronavigation system (NBS 3.2, Nexstim Plc). All subjects underwent an anatomical T1 magnetic resonance imaging session before the mTMS experiments (a gradient recalled echo sequence with voxel dimensions less than or equal to 1 mm). First, for each participant, the APB hotspot was identified as the cortical site beneath the transducer centre resulting in MEPs with maximum peak-to-peak amplitude for single-pulse TMS. The hotspot was obtained with the peak E-field induced by the bottom coil being approximately perpendicular to the central sulcus and with its first phase inducing an E-field from posterolateral to anteromedial direction in the cortex. After defining the hotspot, the transducer was manually rotated to obtain the highest MEP amplitudes with a fixed suprathreshold intensity. Then, the resting motor threshold (MT) was estimated as the minimum stimulation intensity eliciting at least 10 out of 20 MEPs with at least 50 µV of peak-to-peak amplitude 1 . The MT was estimated for stimulation at 0° (MT0°; posterolateral to anteromedial) and 90° (MT90°; posterolateral to anteromedial) E-field orientations (see Fig. 2a for the orientations relative to the head).
The mTMS pulses were monophasic and delivered by a custom-made electronics 15,50 . The single-pulse stimulation intensity was adjusted by varying the capacitor voltage with the waveform phases lasting for 60.0 µs (positive E-field), 30.0 µs (near-zero E-field), and 43.2 µs (negative E-field). The 60 µs rise-time monophasic TMS pulses have been shown to elicit more consistent MEPs across different orientations 47,51 . The transducer temperature was monitored with a thermal infrared camera (i3, FLIR Systems, USA). If the surface temperature reached 41 °C, the transducer was cooled to about 20-30 °C in approximately 5 min.
To measure the E-field orientation effect on the MEP amplitude and onset latency, we applied to the APB hotspot five single pulses at each orientation in 3° steps, i.e., 120 orientations, on 11 subjects. In addition, 20 single pulses were delivered in the 0° (anteromedial) orientation to obtain a reliable reference for each subject. The pulses were divided into ten sequences of 62 pulses each, followed by short breaks of about 5 min. The stimulation intensity was 110% of the MT0°; the interstimulus interval was pseudo-randomized from 4 to 6 s. To evaluate the effect of stimulation intensity on the orientation dependency curves, we applied a similar experimental protocol on five other subjects. MEPs were collected from three single pulses at the APB hotspot at each orientation in 3° steps (i.e., 120 orientations) in each of the three stimulation intensities 110% and 140% of the MT0°, and 120% of the MT90°. For two out of the five subjects, we measured the MEPs for an extra intensity of 120% of the MT0°. The pulses were divided into nine or twelve (extra intensity) sequences of 120 pulses each, followed by short breaks of about 5 min. The interstimulus interval was pseudo-randomized from 2.4 to 2.7 s. The order of stimulation orientation and intensity of the pulses was pseudo-randomized.

Data processing
Data were pre-processed using MATLAB R2017a (MathWorks Inc, USA). MEPs were extracted from the continuous EMG recordings and trials showing muscle pre-activation or movement artifacts greater than ±15 µV within 1000 ms before the TMS pulse were removed from the analysis. For each trial, we computed the MEP peak-to-peak amplitude in a time window 15-60 ms after the TMS pulse and manually annotated the MEP latency. Given the inherent uncertainty in assigning the latency to small MEPs, the latencies from trials with a peak-to-peak amplitude below 50 µV were rejected from the analysis (0.9% of the trials were rejected and 51% of the MEP latencies were not annotated). As expected, the relatively high amount of non-annotated latencies was due to the small MEPs in the ±45-±135° quadrants.
To align the data across subjects, we smoothed the orientation-amplitude curve of each subject with a moving-average filter with a window size of 7 (i.e., considering data for pulses within 21°) and fit a twopeak composite Gaussian to extract the peak orientation. The angle corresponding to the maximum amplitude was subtracted from the original set of angles, ensuring that all subjects had the maximum amplitude at 0°. The average absolute correction across subjects was 8.5° (standard deviation: 6.5°). The correction obtained from the amplitude data was applied to the latency data, maintaining the correspondence between both measures. For the MEP latency, only subjects showing latencies across the entire range of orientations for each intensity were included in this analysis (7 subjects in 110% MT0°, and 5 subjects in 140% MT0° and 120% MT90°). Then, we fit a 2 nd -order real trigonometric polynomial equation on the latency data set of each stimulation intensity and subject, accounting for the periodicity and visible asymmetry of the data across the closed circle: where ( ) is either the measured latency as a function of the stimulation angle , = 2 is the degree of the trigonometric polynomial, and and (0 ≤ ≤ and ≠ 0 or ≠ 0) are the regression coefficients. The regression coefficients and model parameters were computed using the least-square solver implemented in the lmfit 1.0 package. The regression residuals were visualized in a Q-Q plot and inspected for deviations from normality. The data pre-processing and visualization was performed using custom-made scripts written in the Python 3.7 language (Python Software Foundation, USA).

Statistical analysis
We applied a linear mixed-effects model to assess the effect of stimulus intensity and peak orientation on the MEP amplitude model parameters and latency. The intensity and peak orientation were modelled as fixed effects and subjects were modelled as a random effect using restricted maximum likelihood estimation. The linear mixed-effects model accounts for the spurious variations in the measured data due to the inherent inter-subject variability and is less affected by the unequal sample sizes across the tested categories. The p-values of the fixed effects were derived with Satterthwaite approximations in a Type III Analysis of Variance table. Then, we performed post-hoc multiple comparisons using the estimated marginal means with p-value correction for the false discovery rate. The model residuals were inspected for critical deviations from normality with Q-Q plots and homoscedasticity was inspected with a standard versus fitted values plot. Statistical analysis was performed in R 4.0 (R Core Team, Austria) using the lme4 1.1, afex 0.28 packages, and emmeans 1.5 packages. The level of statistical significance was set at 0.05.

Supplementary Material
Supplementary Table 1: Type III Analysis of Variance Table with Satterthwaite's method for the mixedeffect model of the maximum motor evoked potential (MEP) amplitude. Fixed factors were the stimulation intensity (SI) and peak. Interaction between factors is represented as "×" and p-values in bold are smaller than the statistical significance level of 0.05.

Effect
Degrees  Table 12: Multiple comparisons on sigma between the peak orientations. The results are presented as follows, for a fixed stimulation intensity (SI), we compute the sigma difference between the peak orientations. Each comparison has a standard error (SE), degrees of freedom (DoF), t-ratio, and p-value. The SI is given as a percentage of the orientation-specific resting motor threshold at 0° (MT0°) and 90° (MT90°). p-values in bold are smaller than the statistical significance level of 0.05.  Table 15: Multiple comparisons on motor evoked potential (MEP) latency between the peak orientations. The results are presented as follows, for a fixed stimulation intensity (SI), we compute the sigma difference between the peak orientations. Each comparison has a standard error (SE), degrees of freedom (DoF), t-ratio, and p-value. The SI is given as a percentage of the orientation-specific resting motor threshold at 0° (MT0°) and 90° (MT90°). p-values in bold are smaller than the statistical significance level of 0.05.