key: cord-0319576-6pxcwfii
authors: Gomez, Luis J.; Dannhauer, Moritz; Peterchev, Angel V.
title: Fast computational optimization of TMS coil placement for individualized electric field targeting
date: 2020-05-30
journal: bioRxiv
DOI: 10.1101/2020.05.27.120022
sha: b9280715a4a1a9d828c6569dd4492273a7e2b638
doc_id: 319576
cord_uid: 6pxcwfii

Background During transcranial magnetic stimulation (TMS) a coil placed on the scalp is used to non-invasively modulate activity of targeted brain networks via a magnetically induced electric field (E-field). Ideally, the E-field induced during TMS is concentrated on a targeted cortical region of interest (ROI). Objective To improve the accuracy of TMS we have developed a fast computational auxiliary dipole method (ADM) for determining the optimum coil position and orientation. The optimum coil placement maximizes the E-field along a predetermined direction or the overall E-field magnitude in the targeted ROI. Furthermore, ADM can assess E-field uncertainty resulting from precision limitations of TMS coil placement protocols, enabling minimization and statistical analysis of E-field dose variability. Method ADM leverages the reciprocity principle to rapidly compute the TMS induced E-field in the ROI by using the E-field generated by a virtual constant current source residing in the ROI. The framework starts by solving for the conduction currents resulting from this ROI current source. Then, it rapidly determines the average E-field induced in the ROI for each coil position by using the conduction currents and a fast-multipole method. To further speed-up the computations, the coil is approximated using auxiliary dipoles enabling it to represent all coil orientations for a given coil position with less than 600 dipoles. Results Using ADM, the E-fields generated in an MRI-derived head model when the coil is placed at 5,900 different scalp positions and 360 coil orientations per position can be determined in under 15 minutes on a standard laptop computer. This enables rapid extraction of the optimum coil position and orientation as well as the E-field uncertainty resulting from coil positioning uncertainty. Conclusion ADM enables the rapid determination of coil placement that maximizes E-field delivery to a specific brain target. This method can find the optimum coil placement in under 15 minutes enabling its routine use for TMS. Furthermore, it enables the fast quantification of uncertainty in the induced E-field due to limited precision of TMS coil placement protocols. Highlights Auxiliary dipole method (ADM) optimizes TMS coil placement in under 8 minutes Optimum orientations are near normal to the sulcal wall TMS induced E-field is less sensitive to orientation than position errors

Transcranial magnetic stimulation (TMS) is a noninvasive brain stimulation technique [1, 2] , where a TMS coil placed on the scalp is used to generate a magnetic field that induces an electric field (E-field) in the head. This, in turn, directly modulates the activity of brain regions and network nodes exposed to a high intensity E-field [3, 4] . As such, computational E-field dosimetry has been identified by the National Institute of Mental Health as critical for identifying brain regions stimulated by TMS and for developing rigorous and reproducible TMS paradigms [5] . It is important to position and orient the coil to induce a maximal E-field on the targeted cortical region of interest (ROI) [6, 7] . Furthermore, since coil placement protocols have limited precision, it is important to quantify variability in the TMS induced ROI E-field due to potential errors in coil placement. This work proposes the novel auxiliary dipole method (ADM) for fast E-field-informed optimal placement of the TMS coil and for quantifying uncertainty in the TMS induced E-field due to possible coil placement errors.

Optimal TMS coil placement is often determined by using scalp landmarks that correlate with targeted cortical ROIs. For example, to stimulate the dorsolateral prefrontal cortex (DLPFC), the coil is often positioned 5 cm anterior to a position that elicits motor evoked potentials in the contralateral hand muscle [8, 9] . Alternatively, the coil is centered at a 10-20 coordinate location commonly used for EEG electrode positioning [10, 11] . Scalp landmark-based strategies can result in significant misalignments between the coil placement and the targeted cortical ROI [12] .

To improve coil placement, MRI imaging data is also sometimes used for 'neuronavigated' coil positioning. Standard neuronavigated protocols identify the location on the scalp directly over the targeted cortical site's center of mass (CM) as the optimal coil center position [6, 7, 12] . This is because commonly used TMS coils have a figure-8 winding configuration that generates a primary E-field (E-field in the absence of the subject's head) that is strongly concentrated underneath its center [13] . However, the optimum coil placement site on the scalp can be shifted up to 12 mm (5.5 mm on average) away from the scalp location directly above the CM [14] by a secondary E-field generated inside the subject's head due to charge build-up on tissue interfaces [15] [16] [17] .

The orientation of the TMS coil is typically chosen so that the direction of the induced E-field is perpendicular to the ROI's sulcal wall. This orientation is known to maximize the magnitude of the E-field in the ROI [14, 18, 19] and therefore corresponds to the lowest threshold for cortical activation, which is reinforced by the perpendicular orientation of the main axon of pyramidal neurons relative to the sulcal wall [4] . Indeed, E-field directed into the ROI sulcal wall requires, on average, lowest TMS coil currents to evoke a motor potential [20] and is close to the optimal orientation for targeting each region of the hand motor cortical area [14, [21] [22] [23] . Therefore, in the absence of an explicit model of neural activation, choosing a coil position and orientation that maximizes the E-field strength is a suitable objective.

A further limitation of TMS procedures is that they have a limited, and often unquantified, precision and accuracy of determining and maintaining the coil placement. Even with gold-standard neuronavigation and robotic coil placement, coil position error can exceed 5 mm [24] [25] [26] . As such, there is an uncertainty in the TMS induced E-field resulting from uncertainty in the coil placement.

This may result in variability in the outcomes of TMS interventions and needs to be quantified.

Therefore, in addition to linking the external coil placement and current to the E-field induced in the brain, computational models should ideally account for coil placement uncertainty.

The total E-field induced in cortical ROIs can be determined by using MRI-derived subject-specific volume conductor models and finite element method (FEM) or boundary element method (BEM) [23, [27] [28] [29] [30] [31] [32] [33] . Evaluating the TMS induced E-field for a single coil position using a standard resolution head model currently requires 35 seconds using FEM [34] or 104 seconds using a fastmultipole method accelerated BEM [35] . Furthermore, E-field-informed optimal coil placement would require iterative execution of such simulations until an optimal coil position and orientation are found [36, 37] . As such, computational requirements limit the routine use of E-field-informed optimization of coil placement. For example, for this reason we restricted the individual modelbased dosing in a recent study to selecting only the coil current setting [37] . This paper introduces a computational approach that enables fast E-field-informed optimization of coil placement using high resolution individual head models. Our framework is exceptionally computationally efficient as it can evaluate the E-fields generated in a targeted cortical ROI using a standard MRI-derived head model when the coil is placed at 5,900 scalp different positions and 360 coil orientations per position in under 15 minutes using a laptop computer.

To determine the optimum coil placement on the scalp it is unnecessary to evaluate the E-field induced outside of the cortical ROI. This enables the use of the reciprocity principle to compute only the E-field induced in the ROI. Reciprocity has been used previously in the context of boundary element simulations of TMS induced E-fields [38, 39] . These previous uses of reciprocity have two bottlenecks: First, they require the determination of either the magnetic field or E-field due to many electromagnetic point sources, which limits their use to low-resolution head discretizations. Second, they are limited to isotropic head models, and thus they cannot account for white matter tissue anisotropy. Here we avoid the low-resolution bottleneck by using a fastmultipole method (FMM) [40] , which enables the rapid calculation of fields due to many electromagnetic point sources. Furthermore, to enable the use of anisotropic conductivity head models, we apply reciprocity by using directly conduction currents in the brain [41] . These modifications to previous uses of reciprocity in TMS modeling enable the use of head models from commonly adopted transcranial brain stimulation simulation pipelines [17, 31] in our E-fieldinformed coil placement framework. Finally, we introduce a method for approximating the coil currents by auxiliary dipoles enabling it to represent all coil orientations for a given coil position with less than 600 dipoles. This enables the rapid generation of maps that quantify the dependence of the E-field on both coil position and orientation. Furthermore, we post-process these maps to quantify the E-field uncertainty due to the limited accuracy and precision of coil placement methods.

The rest of this paper is structured as follows. First, we describe the proposed approachs to optimally locate and orient a TMS coil on the scalp to maximize E-field components in a brain ROI. Second, the approach is benchmarked in terms of accuracy relative to analytical solutions of a spherical head model and in terms of run-time and memory requirement relative to results obtained by directly using FEM to determine TMS induced E-fields. Third, we use a number of detailed realistic head models and target cortical ROIs from TMS experiments to compare the proposed method to placing the coil in a position that minimizes the scalp to cortical ROI distance.

Finally, we use this approach to estimate quickly the E-field variation related to coil placement errors that inevitably occur with any TMS procedures.

Software implementations of the methods proposed in this paper are available for download [42] .

for use in MATLAB and Python environments enabling easy integration with existing transcranial brain stimulation software packages (e.g., SimNIBS [43] and SCIRun [30] ) as add-on toolkits.

During TMS the coil placement is specified by the position of the coil center and its orientation, which is indicated by a unit vector ̂ that specifies the direction of the E-field directly under its center (Fig.1A) . The average E-field along ̂ (i.e., 〈 ( ) ⋅〉, where 〈⋅〉 computes the average over the ROI) and the average E-field magnitude (i.e., 〈‖ ( )‖〉) in a specified brain region of interest (ROI) are both functions of coil position and orientation. We formulate the goal of optimal TMS coil placement to be the determination of a coil position and orientation ̂ that maximizes either 〈 ( ) ⋅〉 or 〈‖ ( )‖〉 . For example, 〈 ( ) ⋅〉 can be maximized if a preferential field alignment for the targeted neural population is known, whereas 〈‖ ( )‖〉 is maximized whenever this information is missing. Reciprocity method for determining 〈 ( ) ⋅〉 Reciprocity is an equivalence relationship between two scenarios (Fig. 2) . In one scenario, the 

Since for the TMS pulse waveform frequencies the E-field spatial and temporal components are separable (quasi-stationary) [44] , the temporal variation of all currents and E-fields-( ) and

′( ), respectively-is omitted in this and subsequent notation.

Computing 〈 ( ) ⋅〉 using Eq. (2) requires the evaluation of ( ) outside of the conductive head. This is done by computing the primary E-field due to the cortical current ( ) as well as the secondary E-field due to the conduction currents ̿( ) ( ), where ̿( ) is the conductivity tensor inside the head. Thus,

Here 0 is the permeability of free space, and the integration is done over the head domain. In the following section, we describe the FEM approach we used to compute ( ) inside the head, and numerical integration rules used to estimate Eqs. (2)-(3).

The E-field inside the head due to ( ) satisfies Laplace's equation,

where ( ) = ∇ ( ). To solve for and ∇ we used an in-house 1 nd order FEM method [33] .

First, the head model is discretized into a tetrahedral mesh having nodes and tetrahedrons, and each tetrahedron is assigned a constant tissue conductivity tensor. Second, the scalar potential is approximated by piece-wise linear nodal elements ( ) (where = 1, 2, . . . , ) [45] . Third, weak forms of Laplace's equation are sampled also using piecewise-linear nodal elements as testing functions in a standard Galerkin procedure. This results in = ,

Entries ( ) , are computed analytically using expressions provided in [45] and ROI denotes the support of ( ). The system of equations Error! Reference source not found. is solved using a transpose-free quasi-minimal residual [46] iterative solver to a relative residual of 10 −10 .

Samples of ( ) outside of the head are obtained via Eq. (3) and using the FEM solution. The volume cortical and conduction currents are approximated by current dipoles on the centroid of each tetrahedron. The current dipoles are computed as ( ) = (̿ ( ( ) ) ( ( ) ) + ( ( ) )) , where ( ) is the centroid, and is the volume of the j th tetrahedron. When the above conduction current approximation is applied to Eq. (3), it yields the following

Typically coil models consist of dipoles ( ) at locations ( ) , where = 1,2, … , . As a result, using a coil dipole model and Eq. (5) to evaluate Eq. (2) results in

Evaluating Eq. (5) Conversely, if we would like to evaluate the E-field for many coil placements, the resolution of the head model has to be lowered. This is a common limitation of reciprocity based E-field BEM solvers for TMS [38, 39] . To lower the computational cost, we use the fast-multipole method (FMM), which was designed to lower the computational complexity of the calculation of Eq. (6).

Specifically, it reduces the total computation from scaling as × × to scaling as the maximum between and × . Using FMM enables us to compute Eq. (6) with both dense head meshes (i.e., large ) and for many coil placements (i.e., large × ). Furthermore, in the supplemental material we describe an auxiliary dipole model generation used here that leverages the same set of dipole locations for multiple coil orientations further speeding up our method. Note the method as presented here does not work with magnetic dipole TMS coil models. Its extension to magnetic sources requires minor modifications given in the supplemental material. Software implementations for both electric and magnetic dipole sources are available in [42] . Oftentimes, it is desired to maximize the average E-field magnitude in the ROI

For a given coil position and ̂ the best possible approximation of 〈‖ ( )‖〉 that can be obtained while replacing unit vector instead of computing the magnitude and then taking the average over the ROI, the average is first computed and then the magnitude is taken. This approximation is valid whenever the E-field in the ROI is nearly unidirectional. As such, the approximation will be less accurate for larger ROIs relative to smaller ones. We compared results obtained by directly computing 〈‖ ( )‖〉 with ones for ‖〈 ( )〉‖ for ROIs of varying sizes to assess the accuracy of this approximation. Furthermore, we ran SimNIBS [47] coil placement optimization to maximize 〈‖ ( )‖〉 and compared with our results obtained by maximizing ‖〈 ( )〉‖. 

where ‖ − ‖ is the geodesic distance along the scalp between and , and parameters Δ and Δ are radii of the positional and angular support of the distribution. The values of Δ and Δ are specific to the TMS methods used and should be based on available data [24] [25] [26] .

Furthermore, other distributions (e.g., Gaussians) can also be used to discretize the coil position and orientation uncertainty whenever appropriate. The geodesic distance on the scalp is computed using the heat method described in [49] .

For the spherical head model, we compare the average E-field in the ROI predicted by our inhouse direct FEM, reciprocity, and ADM with the analytically computed value. The average Efield analytically computed for each coil placement is shown in Fig. 3 . As expected, the maximum occurs when the coil is placed directly above the ROI and oriented to induce a maximum primary E-field along . The errors normalized to the maximum E-field across all simulations is shown in For the Ernie head model, we compare the average E-field in the ROI predicted by the in-house direct FEM, reciprocity, and ADM. The value of ̂ is chosen (as described above) as (−0.57, −0.72, −0.39) by evaluating the E-field for coil orientations varying from 0° to 359° in 5° intervals at each of the coil positions. The average E-field computed by the direct FEM for each coil position is shown in Fig. 4 , and a maximum E-field is obtained when the coil is placed 4.5 mm away from the point on the grid closest to the ROI center of mass. The normalized absolute error of the different approaches, when using direct FEM as reference, is shown in Fig. 4 . All approaches agree with the direct FEM and have relative absolute errors below 0.13%. The maximum normalized absolute errors across all simulations of the ADM relative to the reciprocity approach is shown in Fig. 5 . The ADM converges to the reciprocity results exponentially with increasing number of dipoles, and an error below 0.5% is observed using the configuration of 578 dipoles. Additional comparisons with SimNIBS results given in the supplemental material have a maximum relative error below 0.11 %. For the four models (M1-M4) we compare metric 〈‖ ( )‖〉 to ‖〈 ( )〉‖ for ROIs of varying sizes. Fig. 6 shows the maximum magnitude error obtained by comparing the ADM to compute ‖〈 ( )〉‖ relative to using direct FEM to compute 〈‖ ( )‖〉 for ROIs of varying sizes. The maximum error increases with increasing ROI size. From what is known of the spatial resolution of TMS experimentally, the diameter of a typical ROI is below 1 cm diameter, and for these sizes the errors are below 0.5 %. For larger ROIs with diameters between 1 cm and 2 cm the error is still below 2.0 %. These results indicate that 〈‖ ( )‖〉 can be approximated by ‖〈 ( )〉‖ in many practical TMS settings. Figure 6 . Maximum E-field magnitude estimate error (Eq. (9)) using ADM to the reciprocity results as a function of ROI diameter.

Coil placement optimizations were run each time using a different number of candidate designs.

The optimizations were run using a high performance computation (HPC) system that has an hours, whereas ADM is 165 times faster, taking a little over 7 minutes. The memory required to run the SimNIBS optimization ranged 11-11.5 GB, whereas with ADM it was 6.3 -6.5 GB. Note the memory requirements of SimNIBS can be lowered to 3.5 GB by using their iterative solver for the optimization [43] , which, however, takes significantly longer. We decided not to use the cm diameter, the discrepancy between ADM and simNIBS optimum TMS coil position is likely due to inaccuracies in the approximation of the E-field magnitude using ADM. For regions with diameter below 4 cm, the orientation differences are likely because ADM was run with an angular step size of 1°, whereas the SimNIBS simulations were run with an angular step size of 5° (due to long computation times). As such, the observed differences in coil orientation between ADM and SimNIBS are below the resolution of our optimization.

In Fig. 8 , ̂ is shown for each ROI. ̂ is always perpendicular to a sulcal wall of the ROI. This

is consistent with what has been found experimentally [14, 18, 19] . Whenever there are multiple bends in the ROI it is difficult to identify the precise normal to the sulcal wall. Our framework helps disambiguate a precise coil orientation that should be used.

In a conventional neuronavigated coil placement protocol the coil would have been placed directly above the center of mass ( ) and oriented to generate a primary E-field along ̂ (i.e. ̂= . Distance between coil placements (D) and , (E) and , and (F) and . Additional coil optimization comparisons given in Tables S1 and S2 in the supplemental material.

As an example, we computed the uncertainty in |〈 ( ) ⋅〉| for the 10 mm diameter ROI in the head model M2 resulting from coil position and orientation uncertainty. Fig. 10 shows results for various levels of coil position uncertainty. For fixed the minimum of the standard deviation is seen where the expected value of the E-field is highest, and the standard deviation at each position increases with . For fixed the expected value decreases monotonically away from its optimum, and this decrease accelerates with increasing . In Fig. 11 , we plotted the 90% confidence interval and expected value for |〈 ( ) ⋅〉| for various levels of position and orientation uncertainty. The expected value of |〈 ( ) ⋅〉| decreases with increased . The range of the 90% confidence interval of |〈 ( ) ⋅〉| is nearly doubled if the coil is placed 5 mm away from its optimal position ( Fig. 11A-E) . Furthermore, coil orientation uncertainty does not significantly contribute to total uncertainty (Fig. 11F ). All of these results indicate that errors in identifying optimum coil position both decreases efficacy of TMS and also increases uncertainty in the TMS induced Efield due possible coil positioning errors. Furthermore, it is more critical to identify the optimum coil positioning whenever there are large coil positioning uncertainties associated with the TMS coil placement protocol. 

ADM, which uses auxiliary dipoles along with electromagnetic reciprocity, enables rapid extraction of optimal coil position and orientation to target the TMS E-field to specific brain regions.

Furthermore, ADM enables the rapid quantification of uncertainty of TMS induced ROI E-fields resulting from uncertainty in the coil position and orientation. The validation results indicate that the average E-field along a given direction in an ROI can be computed with a numerical error below 1% using ADM with 578 dipoles. For ROIs with diameters below 2 cm, the average magnitude of the E-field can be estimated to an error below 2% by executing ADM only three times, corresponding to three orthogonal spatial directions of the E-field. As such, ADM can accurately determine the optimal coil placement based on maximization of either the E-field along a given direction or the E-field magnitude for typical cortical ROIs.

The optimal coil positions and orientations determined via ADM and the direct FEM approach for maximizing the average E-field are virtually identical for ROIs with diameter below 2 cm. To optimize the coil placement for one subject, computations take over two days using FEM, whereas ADM runs in under 15 minutes on a laptop. Thus, ADM is a practical rapid method for E-fieldinformed coil placement.

The difference between the location on the scalp closest to the targeted ROI center of mass and the E-field-informed optimum coil placement is between 1 to 14 mm consistent with prior results using direct simulations that observed an average distance of 5.5 mm and as high as 1.2 cm for ROIs on the temporal brain region [14] . When the objective is to maximize the E-field magnitude in the ROI, the optimal coil orientation is ambiguous for ROIs containing highly curved sulcal walls.

The reciprocity-based approach and ADM provide a fast-computational method that can address all of these issues by enabling E-field-informed coil placement to be done on a standard laptop. Table S2 : Additional comparisons between maximizing E-field magnitude (i.e. maximize 〈‖⋅‖〉) and using ADM to approximately maximize E-field magnitude (i.e. maximize ‖〈⋅〉‖).

Unlike in the main text and Table S1 we compare E-field dose predictions using their respective methods to emphasize that the E-field magnitude approximation ‖〈⋅〉‖ gives a good estimate of actual E-field magnitude.

Transcranial electric and magnetic stimulation: technique and paradigms

Non-invasive magnetic stimulation of human motor cortex

Where does TMS stimulate the motor cortex? Combining electrophysiological measurements and realistic field estimates to reveal the affected cortex position

Simulation of transcranial magnetic stimulation in head model with morphologically-realistic cortical neurons

17-600: early stage testing of pharmacologic or device-based interventions for the treatment of mental disorders (R61/R33)

Transcranial magnetic stimulation in therapy studies: examination of the reliability of "standard" coil positioning by neuronavigation

Optimizing functional accuracy of TMS in cognitive studies: a comparison of methods

Daily repetitive transcranial magnetic stimulation (rTMS) improves mood in depression. Neuroreport: An International Journal for the Rapid Communication of Research in Neuroscience

Rapid-rate transcranial magnetic stimulation of left dorsolateral prefrontal cortex in drug-resistant depression

Using the International 10-20 EEG System for Positioning of Transcranial Magnetic Stimulation

Stimulation over the human supplementary motor area interferes with the organization of future elements in complex motor sequences

Optimal transcranial magnetic stimulation coil placement for targeting the dorsolateral prefrontal cortex using novel magnetic resonance image-guided neuronavigation

Electric field depth-focality tradeoff in transcranial magnetic stimulation: simulation comparison of 50 coil designs

Atlas of optimal coil orientation and position for TMS: A computational study

3D modeling of the total electric field induced by transcranial magnetic stimulation using the boundary element method

Functional localization in the human brain: Gradient-echo, spin-echo, and arterial spin-labeling fMRI compared with neuronavigated TMS

Impact of the gyral geometry on the electric field induced by transcranial magnetic stimulation

The effect of local anatomy on the electric field induced by TMS: evaluation at 14 different target sites

The coil orientation dependency of the electric field induced by TMS for M1 and other brain areas

Optimal Coil Orientation for Transcranial Magnetic Stimulation

Inter-individual variability in optimal current direction for transcranial magnetic stimulation of the motor cortex

Optimal focal transcranial magnetic activation of the human motor cortex: effects of coil orientation, shape of the induced current pulse, and stimulus intensity

Bringing transcranial mapping into shape: sulcus-aligned mapping captures motor somatotopy in human primary motor hand area

Accuracy of robotic coil positioning during transcranial magnetic stimulation

Transcranial magnetic stimulation and the challenge of coil placement: A comparison of conventional and stereotaxic neuronavigational strategies

Navigated transcranial magnetic stimulation

The development and modelling of devices and paradigms for transcranial magnetic stimulation

Fast multigrid-based computation of the induced electric field for transcranial magnetic stimulation

Software Toolkit for Fast High-Resolution TMS Modeling

A pipeline for the simulation of transcranial direct current stimulation for realistic human head models using

2012 Annual International Conference of the IEEE

Realistic vOlumetric-Approach to Simulate Transcranial Electric Stimulation--ROAST--a fully automated open-source pipeline

The ICVSIE: A General Purpose Integral Equation Method for Bio-Electromagnetic Analysis

Conditions for numerically accurate TMS electric field simulation

Field modeling for transcranial magnetic stimulation: A useful tool to understand the physiological effects of TMS?

A software toolkit for TMS electric-field modeling with boundary element fast multipole method: An efficient MATLAB implementation

Online repetitive transcranial magnetic stimulation during working memory in younger and older adults: A randomized within-subject comparison

Site-specific effects of online rTMS during a working memory task in healthy older adults

Real-time computation of the TMS-induced electric field in a realistic head model

Comparison of spherical and realistically shaped boundary element head models for transcranial magnetic stimulation navigation

Computational software: Simple fmm libraries for electrostatics, slow viscous flow, and frequency-domain wave propagation

Efficient computation of lead field bases and influence matrix for the FEM-based EEG and MEG inverse problem

Electric field simulations for transcranial brain stimulation using FEM: an efficient implementation and error analysis

Considerations of quasi-stationarity in electrophysiological systems

The finite element method in electromagnetics

A transpose-free quasi-minimal residual algorithm for non-hermitian linear systems

A novel approach to localize cortical TMS effects

Electric field properties of two commercial figure-8 coils in TMS: calculation of focality and efficiency

Geodesics in heat: A new approach to computing distance based on heat flow

The authors declare that there is no conflict of interest regarding the publication of this article. 

Reciprocity method applied to magnetic dipole coil models Here the reciprocity section of the main text is repeated with appropriate field related quantities replaced for magnetic dipole coil model application. Reciprocity is an equivalence relationship between two scenarios (Fig. S1 ). In one scenario, the TMS coil, modeled as impressed magnetic currents ( ; ) = ′( ) ( ), generates and E-field ( ; ) = ′( ) ( ) inside the head (Fig. S1A ). Computing 〈 ( ) ⋅〉 using Eq. (S2) requires the evaluation of ( ) outside of the conductive head. This is done by computing the primary H-field due to the cortical current ( ), and secondary H-field due to conduction currents ̿( ) ( ), where ̿( ) is the conductivity tensor inside the head. This results in the followingFinally, to estimate Eqs. (S2)-(S3) we use the same approach as described in the manuscript to estimate Eqs. (2)-(3). Specifically, in Eqs. (5)-(6) becomeand, assuming the coil model consists of magnetic dipoles ( ) at locations ( ) , where = 1,2, … , ,(S5) Here we describe the procedure used to approximate all coil orientations using a few dipoles located at auxiliary nodes (Fig. S2) . The dipoles are chosen as Gauss-Legendre tensor product nodes of a rectangular box enclosing the coil (Fig. S2B) . Assume that the x-y axis is chosen as the coil plane, and it consists of dipoles ( ) at locations ( ) , where i = 1,2, … (Fig. S2A ).

From (i) and (i) , we generate coil dipole models ( ( ) ) at locations ( ( ) ) and = 1,2, … , by rotating ( ) and ( ) (Fig. S2B) . Assume that all the coil models can be bound by a Since Ω box is outside of the head, (( ( ) ) ) is interpolated using a polynomial defined from samples in Ω box with an error that decreases exponentially with number of points. We 

Here additional validation results are given indicating that our implementation achieves comparable error levels to SimNIBS, a commonly used non-invasive brain stimulation E-field dosimetry toolbox. In Fig. S3 we show Spherical head model results. The errors observed for the spherical head are slightly higher than our direct FEM. These small differences in accuracy are likely due to differences in the way that the primary E-field is computed. In Fig. S4 we compare the results obtained for the Ernie head model via our direct FEM with SimNIBS. The maximum relative difference across all simulations was 0.11 %. In Fig. S5 we compare results obtained for the four additional subjects obtained via our direct FEM and SimNIBS. The agreement between the two sets of solutions increases with ROI size. This is likely because we compute the average the E-field on the ROI, which is less sensitive to numerical errors with increasing ROI size. In all cases SimNIBS was in agreement with our direct FEM to a fraction of a percentage indicating that both codes result in virtually the same solution.