## Highlights

- •FDM and 1st order FEM with 1.5 mm average mesh edge length have numerical errors >7%.
- •BEM or 2nd (and higher) order FEM with 1.5 mm average mesh edge length achieve numerical errors < 2%.
- •Coil wire cross-section must be accounted to achieve E-field errors below <2%.
- •Coil eddy currents can account for >2% of E-field when very brief pulses are used.

## Abstract

### Background

Computational simulations of the E-field induced by transcranial magnetic stimulation (TMS) are increasingly used to understand its mechanisms and to inform its administration. However, characterization of the accuracy of the simulation methods and the factors that affect it is lacking.

### Objective

To ensure the accuracy of TMS E-field simulations, we systematically quantify their numerical error and provide guidelines for their setup.

### Method

We benchmark the accuracy of computational approaches that are commonly used for TMS E-field simulations, including the finite element method (FEM) with and without superconvergent patch recovery (SPR), boundary element method (BEM), finite difference method (FDM), and coil modeling methods.

### Results

To achieve cortical E-field error levels below 2%, the commonly used FDM and 1st order FEM require meshes with an average edge length below 0.4 mm, 1st order SPR-FEM requires edge lengths below 0.8 mm, and BEM and 2nd (or higher) order FEM require edge lengths below 2.9 mm. Coil models employing magnetic and current dipoles require at least 200 and 3000 dipoles, respectively. For thick solid-conductor coils and frequencies above 3 kHz, winding eddy currents may have to be modeled.

### Conclusion

BEM, FDM, and FEM all converge to the same solution. Compared to the common FDM and 1st order FEM approaches, BEM and 2nd (or higher) order FEM require significantly lower mesh densities to achieve the same error level. In some cases, coil winding eddy-currents must be modeled. Both electric current dipole and magnetic dipole models of the coil current can be accurate with sufficiently fine discretization.

## Keywords

## Introduction

During transcranial magnetic stimulation (TMS) application, a coil placed on the scalp generates a magnetic field that induces an electric field (E-field) in the underlying head and brain [

[1]

]. The delivered E-field dose mediates the physiological effects of TMS; however, there are no established methods for directly measuring it. Therefore, computational modeling has become the dominant approach for TMS E-field dosimetry [[2]

,[3]

]. This work aims to enable the safe and effective use of E-field computer simulations for TMS research by providing a systematic study of their numerical accuracy.Computational methods for TMS simulations require a model of the TMS setup, which consists of the subject’s head, the TMS coil, and their placement relative to each other. Modeling errors in the predicted E-field can result from limited accuracy of the head model, which is typically derived from MRI images that are segmented into individual tissue compartments and using sparse conductivity data [

4

, 5

, 6

, 7

, 8

, 9

, 10

, 11

]. Furthermore, modeling errors can result from limited accuracy in coil geometry and placement representation [- Indahlastari A
- Chauhan M
- Sadleir R

Benchmarking transcranial electrical stimulation finite element simulations: a comparison study.

*J Neural Eng.*2019; 16 (2)https://doi.org/10.1088/1741-2552/aafbbd

[10]

]. Additional to modeling errors, the computational algorithms that calculate the E-field introduce numerical error. Modeling and numerical sources of errors are virtually independent, and the total error in computational simulations can be well approximated as the root sum square of their contributions [[12]

]. Good modeling practice requires that numerical methods do not contribute significantly to the total error [[13]

].The above is particularly vital for computational E-field dosimetry studies that aim to establish causal relationships between TMS induced E-fields and TMS setup parameters (e.g. coil position [

[14]

], coil orientation [[14]

,[15]

], brain size [[14]

], and conductivity [[14]

,[16]

]). These kinds of studies typically fix most other TMS setup parameters and vary only the ones that are being studied. If the numerical error is not negligibly small compared to the effects of the TMS setup parameters being studied, it can be a significant confounder. Correspondingly, the required numerical error level is study dependent. We recommend that simulation studies aim for a numerical error less than 2% to ensure that it is insignificant relative to current levels of modeling error [14

, 15

, 16

].We present a systematic study of the numerical error of the TMS-induced E-field estimated by computational methods, and we provide guidelines for sufficiently reducing it. Specifically, we study the convergence of predicted cortical E-fields with respect to the following numerical approximation variables: coil current approximation method, type of computational method, head mesh density, and computational method approximation order.

For computational modeling of TMS, the induced E-field can be partitioned into a primary E-field due to the coil currents and a secondary E-field due to charge accumulation on regions where there is a change in conductivity [

[17]

]. Numerical error of the primary E-field is dependent on the method used for modeling the TMS coil windings and its currents. Numerical error of the secondary E-field is dependent on the numerical method used—finite element method (FEM) or boundary element method (BEM)—and the specific type of FEM or BEM discretization.Typical TMS coils consist of circular windings of wire with rectangular cross-section. For example, a commonly used Magstim 70 mm Double Coil (P/N 9925-00, Magstim Ltd., UK) consists of a pair of 70 mm average diameter, nine-turn circular spiral windings of wire with cross-section of 1.9 mm by 6 mm. Because of their small cross-section most TMS modeling frameworks assume that the currents flow uniformly through the wire [

18

, 19

, 20

]. However, the eddy currents cause the currents to redistribute and concentrate toward the wire surface [- Wang B
- Shen MR
- Deng Z-D
- Smith JE
- Tharayil JJ
- Gurrey CJ
- et al.

Redesigning existing transcranial magnetic stimulation coils to reduce energy: application to low field magnetic stimulation.

*J Neural Eng.*2018; 15https://doi.org/10.1088/1741-2552/aaa505

[21]

]. Neglecting this eddy-current redistribution can potentially result in errors of the predicted primary E-field. Furthermore, some modeling frameworks use Faraday’s law to model the coil as magnetic dipoles, thereby, introducing another possible source of numerical error [[18]

,[22]

,[23]

]. We evaluate different approximations of modeling the coil primary E-field, including the impact of eddy currents and explicit representation of the coil current flow in contrast to approximation through magnetic dipoles.In FEM, the head is discretized into a mesh consisting of polyhedral volumetric elements [

[24]

] (e.g., tetrahedrons or hexahedrons). Then, a weak form of Poisson’s equation [[18]

] is solved numerically to determine the secondary E-field. A commonly used subset of FEM called finite difference method (FDM) approximates the head by a mesh consisting of equal sized cuboid elements, thereby approximating tissue boundaries by ‘staircases’ [[15]

,[25]

]. These staircase approximations can be problematic since accurately resolving tissue boundaries requires many elements [[15]

], and artifactual outliers are typically observed in the FDM solution [[15]

,25

, 26

, 27

, 28

]. In contrast, tetrahedral meshes can be used to better conform to the smooth tissue boundaries, and many existing software packages use them for TMS modeling [[18]

,29

, 30

, 31

].For BEM, the subject’s head tissue-boundary surfaces are discretized by two-dimensional polyhedral elements (triangles or quadrilaterals). Weak forms of surface integral equations are solved to determine the TMS induced E-field. FEM and BEM solve a system of equations for a piecewise polynomial approximation of a quantity related to the secondary E-field (e.g., charge distributions or scalar potential) on the mesh. Because BEM only discretizes tissue surfaces, it typically requires less elements than FEM. Nevertheless, in BEM computer memory requirements increase quadratically with the number of mesh elements, limiting its use [

[32]

]. Recently, an implementation of the adjoint double layer BEM [33

, 34

, 35

, 36

] using fast-multipole-method (FMM) accelerators, which reduces memory requirements to linearly increasing with the number of mesh elements [[37]

,[38]

], has been shown to be computationally tractable for analyzing realistic TMS scenarios [- Htet AT
- Saturnino GB
- Burnham EH
- Noetscher G
- Nummenmaa A
- Makarov SN

Comparative performance of the finite element method and the boundary element fast multipole method for problems mimicking transcranial magnetic stimulation (TMS).

*J Neural Eng.*2019; 16https://doi.org/10.1088/1741-2552/aafbb9

[39]

]. A remaining limitation of BEM is that, unlike FEM or FDM, it is not effective at modeling highly heterogeneous media such as anisotropic tissue [[40]

]. For highly heterogeneous media, extensions of BEM termed volume integral equation methods can be used [[31]

,[41]

,[42]

]. However, like FEM and FDM, these methods need volumetric meshes, which makes them more computationally expensive than BEM [[40]

].FDM, FEM, and BEM are known to converge to the actual solution by either refining the mesh (H-refinement) and/or increasing the polynomial approximation order (P-refinement) [

[43]

]. Because of the relatively complex geometry of the human brain, it is difficult to estimate the P- and H-refinement necessary to compute accurate E-field representations in the brain with FDM, FEM, and BEM. Therefore, we also evaluate via simulations the impact of P- and H-refinements on the TMS E-field modeling accuracy with FDM, FEM, and BEM, and compare them to a reference solution.In summary, in this paper we benchmark various TMS E-field simulation methods and determine the conditions for achieving specific error levels. We compare the absolute accuracy and computation time of FDM, FEM, and BEM implementations with respect to analytical solutions of a spherical head model and a high-fidelity simulation of a realistic head model. We also explore the impact of the TMS coil representation on the solution accuracy.

## Methods and materials

### Approximation of electromagnetic equations

We used quasistatic FEM as in Refs. [

where $\text{\sigma}\left(\mathbf{r}\right)$, $\phi \left(\mathbf{r}\right)$, and ${\mathbf{E}}_{\mathbf{p}}\left(\mathbf{r}\right)$ are the conductivity, scalar potential, and primary E-field at position $\mathbf{r}$, respectively. Both FEM and FDM solve Eq. (1) by approximating the head by a mesh consisting of elements (tetrahedra for FEM and identical cuboidal elements for FDM), assuming a piecewise polynomial approximation to the scalar potential within each element, and determining the piecewise polynomial expansion by applying weak forms of Eq. (1). The convergence of the FEM and FDM approximation with respect to the true solution depends on the size of the individual elements used to approximate the head and the polynomial order of the solution [

[18]

,[29]

,[30]

] and a 27 point stencil FDM, which is equivalent to the FEM with multilinear basis functions and identical cuboidal elements, as in Ref. [[25]

]. In this formulation, displacement currents are neglected, as they are much smaller than conduction currents. Magnetic fields scattered by the human head are also neglected, as they are small compared to the magnetic field produced by the TMS coil. Within these assumptions, the primary E-field can be determined from the coil currents via the Biot–Savart law [], and the secondary E-field arises purely from a scalar potential field. The scalar potential can be determined by solving Poisson’s differential equation from the primary E-field,$$\nabla \cdot \sigma \left(\mathbf{r}\right)\nabla \phi \left(\mathbf{r}\right)=\nabla \cdot \sigma \left(\mathbf{r}\right){\mathbf{E}}_{\mathbf{p}}\left(\mathbf{r}\right)\text{,}$$

(1)

where $\text{\sigma}\left(\mathbf{r}\right)$, $\phi \left(\mathbf{r}\right)$, and ${\mathbf{E}}_{\mathbf{p}}\left(\mathbf{r}\right)$ are the conductivity, scalar potential, and primary E-field at position $\mathbf{r}$, respectively. Both FEM and FDM solve Eq. (1) by approximating the head by a mesh consisting of elements (tetrahedra for FEM and identical cuboidal elements for FDM), assuming a piecewise polynomial approximation to the scalar potential within each element, and determining the piecewise polynomial expansion by applying weak forms of Eq. (1). The convergence of the FEM and FDM approximation with respect to the true solution depends on the size of the individual elements used to approximate the head and the polynomial order of the solution [

[43]

]. TMS modeling packages [[18]

,[29]

,[30]

] typically employ piecewise linear (1st order) approximations for the scalar potential, resulting in piecewise constant approximations of the E-field within the head. Relative to higher order methods, these piecewise constant approximations converge slowly to the true solution [[43]

]. We additionally evaluated 2nd and 3rd order FEM formulations that result respectively in piecewise linear and piecewise quadratic approximations of the E-field within the head. We also used superconvergent patch recovery (SPR) to increase the accuracy of 1st order FEM [[13]

,[45]

].- Saturnino GB
- Madsen KH
- Thielscher A

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

*J Neural Eng.*2019; https://doi.org/10.1088/1741-2552/ab41ba

For BEM, we used an FMM-accelerated BEM formulation for analysis of TMS, which solves integral equations in terms of surface charges derived by imposing weak forms of current continuity equations on tissue boundaries [

[39]

]. Then, the surface-charge induced E-field is evaluated using Coulomb’s law. We implemented the formulation with a piecewise constant (0^{th}order) [[39]

] as well as a piecewise linear (1st order) approximation of the charges.All solvers where developed in-house, with implementation details provided in the Supplement. All of the solvers are freely available online at [

[46]

]. The FMM acceleration of the BEM is achieved using FMM libraries [[47]

].### Head modeling

Spherical head models: we employed concentric spherical shell models [

[48]

,[49]

] composed of a spherical core and three spherical shell compartments with inner to outer layers representing brain, cerebrospinal fluid (CSF), skull, and scalp, with each with respective outer boundary radius of 78, 80, 86, and 92 mm. We considered an inhomogeneous sphere model where the brain, CSF, skull, and scalp layer conductivity is set to 0.33, 1.79, 0.01, and 0.43 S/m, respectively [- Wolters C.H.
- Anwander A.
- Tricoche X.
- Weinstein D.
- Koch M.A.
- MacLeod R.S.

Influence of tissue conductivity anisotropy on EEG/MEG field and return current computation in a realistic head model: a simulation and visualization study using high-resolution finite element modeling.

*Neuroimage.*2006; 30: 813-826

[8]

]. Additionally, we considered a homogenous sphere model comprised of a single-shell sphere with radius of 92 mm and conductivity of 1 S/m. For both versions of sphere models we used the same FEM meshes, and utilize the same meshing software, GMSH [[50]

], for discretization. The resulting mesh parameters are given in Table 1, with additional details in supplemental Table S.3. Plots of each mesh are provided in Supplemental Figs. S13–S19. Additionally, we generated FDM meshes consisting of cubic voxels, and centered about the spherical head by assigning each voxel the head model conductivity corresponding to the voxel center. For each mesh, we ran a simulation scenario where the sphere head is centered about the origin, and the bottom of the coil windings is centered 5 mm above its apex at location (97, 0, 0) mm, as shown in Fig. 1A. Analytical expressions for the E-field generated inside the spherical head model [[51]

] were used to determine accuracy of all solvers for the sphere head at each mesh density level.Table 1Spherical head model mesh parameters.

Mesh refinement level | Number of tetrahedrons | Edge length (mm) | Number of homogeneous sphere BEM triangles | Number of inhomogenous sphere BEM triangles | |
---|---|---|---|---|---|

Average | Maximum | ||||

0 | 11,496 | 14.6 | 32.4 | 1202 | 3972 |

1 | 54,258 | 8.5 | 17.5 | 4668 | 15,642 |

2 | 319,304 | 4.8 | 9.0 | 20,758 | 68,220 |

3 | 1,613,233 | 2.8 | 5.1 | 62,706 | 210,714 |

4 | 11,286,205 | 1.5 | 2.7 | 233,524 | 785,144 |

5 | 35,108,346 | 1.0 | 1.8 | 501,914 | 1,682,698 |

MRI-derived head model: As a head mesh, we used the one provided in the popular SimNIBS v1.0 package [

[18]

]. (In SimNIBS v2.0 this mesh was updated with one having nearly doubled resolution; however, we employed the lower-resolution version since it enables more barycentric refinements while maintaining computational tractability.) The mesh consisted of five homogenous compartments including white/gray matter (WM/GM), CSF, skull, and skin with conductivity of 0.126, 0.276, 1.65, 0.01 and 0.465 S/m, respectively [[52]

]. Starting from the SimNIBS v1.0 example mesh, we generated refined meshes, each time by subdividing each tetrahedron of the earlier mesh into eight equal sized tetrahedrons using the GMSH ‘refine by splitting’ option [- Gomez L
- Goetz S
- Peterchev AV

Design of transcranial magnetic stimulation coils with optimal trade-off between depth, focality, and energy.

*J Neural Eng.*2018; 15https://doi.org/10.1088/1741-2552/aac967

[50]

]. Mesh parameters for each refinement level are given in Table 2. The coil was placed 5 mm above the primary motor cortex hand-knob region of the brain as shown in Fig. 1B. We generated FDM meshes 0, 1, 2, 3, and 4 consisting of 2.8 mm, 1.4 mm, 0.8 mm, 0.4 mm, and 0.2 mm cubic voxels, respectively. A conductivity was assigned to each voxel from the MRI-derived head mesh by assigning each voxel the conductivity of its voxel center.Table 2MRI-derived head model mesh parameters.

Mesh refinement level | Number of tetrahedrons | Edge length (mm) | Number of BEM triangles | ||
---|---|---|---|---|---|

Minimum | Average | Maximum | |||

0 | 877,106 | 0.2 | 2.9 | 15.7 | 251,888 |

1 | 7,016,848 | 0.1 | 1.5 | 14.1 | 1,007,552 |

2 | 56,134,784 | 0.06 | 0.8 | 7.2 | 4,030,208 |

3 | 449,078,272 | 0.03 | 0.4 | 5.3 | 16,120,832 |

### Coil modeling

We considered two models of 70 mm loop diameter figure- 8 coils. Both models assume that the coil is made of rectangular cross-section copper wire and consist of two circular windings side by side, where each winding has 9 concentric turns and inner/outer diameter of 52/88 mm, respectively. One model has a ‘thick’ cross-section of 1.9 mm by 6 mm, as measured from a Magstim coil winding. The other model has a ‘thin’ cross-section of 1 mm by 7 mm [

[23]

]. Models of the coils were constructed and meshed in GMSH. Same as for the head models, we constructed a hierarchy of meshes each time by using the ‘refine by splitting’ option in GMSH.The currents flowing through the coil windings were determined by solving the current continuity equation via FEM. Additionally, eddy-current effects were modeled using a time-harmonic partial-element equivalent circuit method [

[21]

] developed in-house (details provided in the Supplement). We computed the primary E-field generated by the coil inside the spherical head region (i.e. on a spherical region with apex 5 mm below the center of the coil surface) using a first order accurate quadrature for each tetrahedron of the coil mesh. Primary E-field samples were taken on a grid with points 2 mm apart. Furthermore, we compared the total E-field predicted by the 2nd order FEM with the MRI-derived head mesh with one refinement level.Additionally, we modeled the thick cross-section coil with equivalent magnetic dipoles [

[43]

]. To generate sets of dipoles with varying density, where each dipole represents a small subregion of the coil, we used dipole-placement rules generalized from the layout used in Ref. [[19]

]. First, both circular loops of the figure- 8 coil were divided into ${N}_{ring}$ rings of uniform width. Then, to obtain most uniform set of dipoles, each ring was discretized with $\frac{{r}_{o}+{r}_{i}}{{r}_{o}-{r}_{i}}\pi +\frac{1}{2}$ equispaced dipoles (where ${r}_{i}$ is the inner radius of the ring and ${r}_{o}$— its outer radius). Lastly, ${N}_{layer}$ identical copies of the dipole rings were stacked uniformly between the bottom and top of the coil windings.For the results comparing head E-field modeling methods, we did not include eddy-effects and computed the primary E-field of the figure- 8 coil using an integration rule with 193,536 sample points, which resulted in L

^{2}norm error below 0.1% for the primary E-field generated in the head. We used the same coil model for these simulations to remove coil modeling as a contributor to numerical error when comparing E-field solvers.### Accuracy benchmarks

To quantify the accuracy of our solution we considered two metrics. The first metric is relative L

where $\cdot $ denotes magnitude, ${\mathbf{E}}_{\mathbf{r}\mathbf{e}\mathbf{f}}\left(\mathbf{r}\right)$ and ${\mathbf{E}}_{\mathbf{a}\mathbf{p}\mathbf{p}\mathbf{r}\mathbf{o}\mathbf{x}}\left(\mathbf{r}\right)$ are the reference and approximate solution, respectively, and $\text{\Omega}$ is the domain of integration. For head simulations the integration was done over the GM/WM matter region (i.e., $\text{\Omega}=Brain$) and for coil current error—in the coil winding wire (i.e., $\text{\Omega}=Wire$). The L

^{2}-norm error,$$er{r}_{{L}^{2}\left(\Omega \right)}=\sqrt{\frac{{\iiint}_{\Omega}{\mathbf{E}}_{\mathbf{approx}}\left(\mathbf{r}\right)-{\mathbf{E}}_{\mathbf{ref}}{\left(\mathbf{r}\right)}^{2}d\mathbf{r}}{{\iiint}_{\Omega}{\mathbf{E}}_{\mathbf{ref}}{\left(\mathbf{r}\right)}^{2}d\mathbf{r}}}\phantom{\rule{0.25em}{0ex}}\text{,}$$

(2)

where $\cdot $ denotes magnitude, ${\mathbf{E}}_{\mathbf{r}\mathbf{e}\mathbf{f}}\left(\mathbf{r}\right)$ and ${\mathbf{E}}_{\mathbf{a}\mathbf{p}\mathbf{p}\mathbf{r}\mathbf{o}\mathbf{x}}\left(\mathbf{r}\right)$ are the reference and approximate solution, respectively, and $\text{\Omega}$ is the domain of integration. For head simulations the integration was done over the GM/WM matter region (i.e., $\text{\Omega}=Brain$) and for coil current error—in the coil winding wire (i.e., $\text{\Omega}=Wire$). The L

^{2}norm error is a measure of the relative energy of the error of the E-field (i.e., the residue), thereby quantifying the global convergence of the numerical solution. For the spherical head model, the brain L^{2}norm error was estimated by sampling the E-field on a grid with points 2 mm apart. For the MRI-derived head, the brain L^{2}norm error was estimated by applying a 2nd order accurate integration rule [[53]

] in each tetrahedron of the original mesh.The second metric is the pointwise magnitude error

$$err\left(\mathbf{r}\right)=\frac{{\mathbf{E}}_{\mathbf{approx}}\left(\mathbf{r}\right)-{\mathbf{E}}_{\mathbf{ref}}\left(\mathbf{r}\right)}{\underset{\mathbf{r}\in \Omega}{\mathrm{max}}{\mathbf{E}}_{\mathbf{ref}}\left(\mathbf{r}\right)}\phantom{\rule{0.25em}{0ex}}\text{.}$$

(3)

This error quantifies the magnitude of the E-field difference at each point relative to the E-field maximum, thereby, characterizing the local convergence of the numerical solution.

## Results

### Spherical head model

For the inhomogeneous and homogenous sphere simulations, the maximum pointwise E-field error on a spherical shell 0.5 mm into the cortex (i.e., with radius 77.5 mm) and ‘‘brain” L

^{2}errors are shown in Fig. 2. For the homogenous sphere, the maximum pointwise and L^{2}errors are similar. FDM, 0^{th}order BEM, and 1st order BEM achieve similar accuracy to 1st order FEM, 2nd order FEM, and 3rd order FEM, respectively. SPR applied to 1st order FEM (SPR-FEM) results in error performance between 1st and 2nd order FEM, decreasing L^{2}error by a factor of 6–13 compared to 1st order FEM. Furthermore, increasing the order of FEM or BEM decreases the L^{2}error by 1–2 orders of magnitude.Compared to the homogenous sphere (Fig. 2A–C), for the inhomogeneous sphere (Fig. 2B–D) FDM errors increase by at least a factor of 3; 1st order FEM remains virtually unchanged; 2nd order FEM and 0

^{th}order BEM errors increase by at most a factor of 2; and accuracy gains from increasing FEM and BEM order from 2nd to 3rd and 0^{th}to 1st, respectively, disappear. SPR-FEM again performs between 1st and 2nd order FEM, decreasing the L^{2}error of 1st order FEM by a factor of 6–13. For FDM a systematic error due to the staircase approximation of each compartment boundary causes the decreased accuracy in the inhomogeneous case and a relatively constant maximum pointwise error with increasing mesh resolution. The low error differences in the inhomogeneous case indicate that for orders beyond 2nd for FEM and 0^{th}for BEM, the numerical error near tissue interfaces is dominated by spatial discrepancies between the mesh boundary and the ideal sphere geometry. To achieve significant accuracy gains by increasing order beyond 2nd for FEM or 0^{th}for BEM, it is likely necessary to use higher-order mesh element shape functions [[54]

] to approximate the head model geometry.The sphere is a smooth geometry with minor geometrical features, resulting in a smooth E-field that rapidly converges to the true solution. As a result, we observe errors below 2% even for mesh edge lengths as high as 1.4 cm. The next section analyzes solution convergence in a more realistic head model.

### MRI-derived head model

We consider the E-field simulations for the MRI-derived head models using 0

^{th}order BEM and 3rd refinement mesh as a reference solution, since it was the most accurate test case we could run within our computational resources. The E-field magnitude 0.5 mm below the cortical surface is shown in Fig. 3. Other than the lowest resolution FDM result, which predicts far larger stimulated area extending to adjacent gyri in frontal and central directions, all simulations exhibit a visually identical E-field pattern enabling their use for qualitative analysis of TMS induced E-fields. In Fig. 3, the E-field maximum (${\left|\mathbf{E}\right|}_{max}$) is given below each plot. The reference has ${\left|\mathbf{E}\right|}_{max}=$100 V/m. For FDM, ${\left|\mathbf{E}\right|}_{max}$ for lowest to highest resolution is 168, 118, 132, 109, and 101 V/m, respectively. Thus, all FDM results with resolution above 0.2 mm have a significant (>9%) error. This also highlights the necessity of removing outlier E-field values (e.g. above 99% percentile) in FDM [[25]

]. The higher-order FEM methods (2nd and 3rd order) and 0th order BEM method predict an ${\left|\mathbf{E}\right|}_{max}$ between 98 and 102 V/m (i.e. with less than 2% error) already at the initial mesh resolution. With 1st order BEM, SPR–FEM, 1st order FEM, and FDM similar accuracy is only obtained for much denser meshes with average edge lengths reduced from 2.9 mm to 1.5 mm, 0.8 mm, 0.4 mm, or below the densest tested mesh at 0.4 mm, respectively.The pointwise E-field error 0.5 mm into the cortical surface is shown in Fig. 4. For the first three refinements, the FDM pointwise error exhibits large regions with errors above 2%; as such, this solution is likely not reliable on this surface. The fourth refinement of FDM improved relative to the other three indicating that, like FEM and BEM, it converges to the correct solution with a sufficiently dense mesh.

The 1st order FEM solution slowly improves as the mesh is refined. After two refinements, significant regions near the sulci have a pointwise error above 2%, consistent with prior observations [

[55]

]. SPR decreases the L^{2}norm of the pointwise errors of 1st order FEM on the cortical surface by a factor of 1.6 and 2.1 for mesh 1 and 2, respectively, without introducing significant computational overhead [[13]

]. For mesh 0, SPR-FEM has L^{2}norm of the pointwise errors that is only 10% higher than that of 1st order FEM. SPR-FEM did not exhibit significant regions with error above 2% after two mesh refinements, whereas three refinements are needed for this accuracy for 1st order FEM. The two and three mesh refinements required 56 and 400 million elements, respectively, in contrast to the typical 3–5 million elements used. All other methods exhibit errors below 2% after one refinement, and after two refinements most of the cortical shell have errors below 0.1%. Additional results of the E-field error on a cut slice into the brain (given in the Supplement) indicate that the above findings about the relative performance of the various methods are also valid for deep brain regions.The L

^{2}norm errors are shown in Fig. 5. From least to most accurate: FDM, 1st order FEM, SPR-FEM, 2nd order FEM, 1st order BEM, 3rd order FEM, and 0^{st}order BEM. The accuracy of FDM lags the accuracy of 1st order FEM by one refinement. Both FDM and 1st order FEM lag all other methods in accuracy by at least two refinements. SPR exhibits super-convergence; i.e. with decreasing average edge length, SPR converges more rapidly than 1st order FEM to the reference solution. 2nd order FEM has an error that is nearly twice that of 1st order BEM. Like the inhomogenous sphere case, 1st order BEM, 0^{th}order BEM, and 3rd order FEM are marginally different in error.Finally, to compare the computational demands within simulation method, each simulation was run using a single core of an Intel(R) Xeon(R) Gold 6154 CPU computer node with 3.0 GHz 36 core processor and 768 GB of memory. Fig. 6 shows the CPU runtime for each of the methods. We divide the CPU runtime into two phases: the setup phase, where a linear system of equations is assembled, and the solution phase, where the system of equations is solved iteratively. For FEM methods, increasing the order or refining the mesh result in similar CPU runtime increases. For BEM methods, the setup phase of the 1st order BEM requires roughly three times the CPU runtime relative to 0

^{th}order BEM, and the solution phase requires nearly the same amount of CPU runtime for both methods. The CPU runtimes of all methods could be significantly reduced by using lower level languages and parallelizing the computational routines. Since FEM/FDM versus BEM implementations are significantly different, the time comparisons here are intended for illustration of their dependence on the mesh density. A valid comparison of timing across solvers would require optimizing the implementations for CPU resources, which has not been done here. For example, we observed a tenfold reduction of the CPU runtime of the solution phase of BEM using shared-memory parallelism with ten cores. Further, FDM runtimes and memory requirements could be reduced using techniques described in Ref. [[56]

]. More detailed discussion of CPU runtime, memory usage and possible speed-ups is given in the Supplement.With each mesh refinement the CPU runtime for FEM/FDM increases by a factor of 16 and for FMM compressed BEM it increases by a factor of 4. In FEM/FDM, this factor can be reduced to its optimal value of 8 by using preconditioning techniques (e.g. multigrid methods [

[56]

]). Nevertheless, BEM implementations that leverage FMM are still asymptotically superior to FEM/FDM in terms of CPU runtime.### Coil modeling comparison

We consider the error of the primary E-field resulting from numerical errors in the current modeling for the coil with thick wire cross-section. As reference we use a high-resolution simulation result from the coil discretized with 12,124,160 tetrahedron elements having an average edge length of 0.5 mm. The L

^{2}norm and maximum pointwise error of the primary E-field induced in the spherical head model for each coil winding discretization is shown in Fig. 7. The magnetic dipole approach is more efficient than the current discretization approach, as it requires approximately an order of magnitude fewer dipoles for comparable error levels. Nevertheless, the magnetic dipole approach used here does not account for gaps between coil turns, and the L^{2}error stagnates near 0.2%. Modeling fine geometrical details and more complicated coils requires partitioning the flux integral domain, which is not always straightforward. For the current discretization approach, the primary E-field can be determined directly using a computer aided design model of the coil. Thus, the current discretization approach provides flexibility at the cost of slightly increased computation. This cost can be virtually removed by precomputing the primary E-field and storing it [[18]

] or using FMM methods. To achieve errors below 2%, it is recommended to use >3000 uniformly distributed current dipoles or >200 uniformly distributed magnetic dipoles.Next, we consider numerical errors in the primary E-field stemming from ignoring eddy current effects in the coil current modeling. We consider the figure- 8 coils with both thin (1 mm by 7 mm) and thick (1.9 mm by 6 mm) cross-section wires at a frequency of 3 kHz and 6 kHz. The currents along a coil cross-section are shown in Fig. 8. There is significant current redistribution along the vertical direction for both coils. In the horizontal direction, there is significantly more variation in the distribution of the current for the thick than thin wire coil.

The L

^{2}and pointwise differences of the primary E-field sampled in the spherical model were both below 0.4% for 3 kHz and 1.0% for 6 kHz for the thin wire coil, but increased for the thick model respectively to 1.7% and 2.6% for 3 kHz, and 3.2% and 5.0% for 6 kHz. Comparing the total E-field generated in the MRI-derived head model brain by the thick wire coil with and without eddy current redistribution revealed L^{2}and maximum pointwise differences of 1.6% and 1.9% for 3 kHz, and 3.1% and 3.8% for 6 kHz, respectively. Therefore, when the TMS pulse has frequencies beyond about 3 kHz, the introduction of eddy currents is likely important to achieve a numerical error below 2%. The primary E-fields induced by the thin and thick wire coils (without including eddy effects) differed by more than 3%, confirming that accurate representation of the wire geometry is important [[57]

].Finally, in the Supplement we compare a Magstim figure- 8 spiral coil model with feed and without one, and the primary E-field differs by 2.6%. This indicates that the feed may have to be modeled to achieve errors below 2% in the coil primary E-field. Including the coil feed in the analysis is relatively straightforward using the approaches given in the Supplement.

## Discussion

FEM, BEM, and FDM are commonly used for TMS E-field dosimetry and to study E-field sensitivity to TMS setup parameters. The accuracy of these methods is limited by modeling and numerical errors. Unlike modeling error, numerical error can be easily made arbitrarily small, and ideally it should not significantly affect solution accuracy. Based on simulation studies [

14

, 15

, 16

], we consider numerical error below 2% to be insignificant relative to other possible contributors to error. We systematically studied error in E-field simulations and determined solver parameters to sufficiently suppress numerical error. As TMS coil positioning systems, brain imaging techniques, and in-vivo electrical conductivity data are improved [[58]

], lower numerical errors might become necessary.In TMS modeling there are three main sources of numerical error: coil modeling fidelity, mesh resolution, and accuracy of the numerical method. Most methods typically employed for modeling coil currents sufficiently suppress numerical errors. TMS coil wire usually has a small cross-section or is of a litz wire, and therefore eddy current effects are typically negligible for TMS simulation. However, to attain TMS coil E-field errors below 2%, modeling eddy currents might be necessary in extreme cases when the wire is made of a solid conductor with a relatively large cross-section, and/or very brief pulses are used.

The mesh resolution mediates both how well the mesh approximates the tissue geometry and how accurate the solution is. FDM meshes typically do not conform well to tissue boundaries, resulting in additional errors relative to FEM and BEM methods. Consequently, to achieve accurate representation of the cortical features, FDM requires denser meshes relative to FEM and BEM. For example, in our results FDM required a voxel edge length of 0.2 mm compared to an average element edge length of 1 mm for second-order FEM to reach the 2% L

^{2}accuracy target. Nevertheless, FDM methods can operate directly with segmented MRI images, which are easily constructed. Furthermore, there are many FDM algorithms that efficiently leverage massively parallel systems and speed up iterative convergence by using multi-grid, Fourier transforms, or FMM [[59]

]. Thus, solving FDM equations on fine grids can be computationally tractable.In this study we did not consider E-fields in the scalp or other body parts. As such, our numerical error estimates are not directly applicable to scalp or peripheral stimulation, although one would expect similar relative performance of the various numerical methods there. Furthermore, additional results (given in the Supplement) indicate that low pointwise errors for E-fields in the cortical surface translate to low pointwise errors for E-fields in deep brain regions. However, these results also indicate that to obtain low absolute E-field errors deep into the brain will require an even lower overall numerical error level than the one recommended here.

Head models with spherical geometry cannot be used to estimate the numerical errors in a realistic TMS setup; for a given mesh resolution, the numerical errors even in the inhomogeneous spherical geometry are at least a factor of ten smaller than in a realistic head geometry. Unsurprisingly, 1st order FEM methods converge slowly with respect to increasing mesh resolution. This is because they employ a piecewise constant approximation of the E-field within each mesh element. This slow convergence can be partially improved by interpolating the 1st FEM solution using the super-convergent patch recovery method [

[60]

]. Our 1st order SPR-FEM errors are similar to those reported in another study [[61]

]. SPR-FEM does not improve accuracy whenever tissue compartments are one tetrahedron thick, and the relative accuracy gains of SPR-FEM increase with the number of tetrahedrons contained within the thickness of a compartment. For example, the MRI-derived head model mesh 0 has a GM compartment that is only one tetrahedron between its CSF and WM boundaries, and we did not observe accuracy improvements from using SPR-FEM in this case.The generation of the BEM system of equations requires the approximation of integral by quadrature, and the use of low-order quadratures can result in errors similar to those of 1st order FEM methods (see Supplement). Increasing the order or refining the mesh of FEM results in similar CPU cost increases. Increasing the order in FEM reduces the error more rapidly than increasing the mesh resolution. 2nd and 3rd order FEM, and 0

^{th}and 1st order BEM exhibit similar levels of accuracy near tissue interfaces. The lack of improvement through P-refinement of the BEM methods could be due to excessive regularity imposed through the use of piece-wise linear basis functions, as discretizations of the adjoint BEM operator should admit piece-wise discontinuous functions [[33]

]. Because of the complicated cortical geometry, a dense mesh with average edge length ≤1.5 mm is required to sufficiently resolve the geometric features [[61]

]. The use of a 2nd order FEM or 0^{th}order BEM method may diminish the numerical error to acceptable levels without the need of additional mesh refinements in practical applications. Furthermore, SPR-FEM increases the accuracy of 1st order FEM without adding significant computational overhead making it a possible alternative to higher order methods. Although achieving 2% error will require additional computational resources, this increase is smaller than the computational resources already required to generate a tetrahedral mesh from an MRI image. For TMS E-field analysis applications requiring ensembles of simulations like uncertainty quantification [[9]

,[16]

], functional mapping [[62]

], and neuro-navigation [[63]

], improving the computational cost to achieve a low numerical error remains an open research area.## Conclusions

We benchmarked the accuracy of various methods for simulating the TMS induced E-field as a function of the density of spatial discretization. Whereas at present 1st order FEM is most commonly used, 2nd order FEM or 0

^{th}order BEM appear to be judicious strategies for achieving negligible numerical error relative to modeling error, while maintaining tractable levels of computation, and SPR-FEM offers some advantages as well. FDM requires significantly higher mesh densities than FEM and BEM to achieve a given accuracy level, likely due to the staircase approximation of tissue boundaries. In addition to the numerical errors arising from the solver, if the TMS coil windings use a large cross-section solid wire, the eddy currents within the coil windings may have to be modeled to reach error levels below 2%.## Funding

Research reported in this publication was supported by the National Institute of Mental Health and the National Institute of Neurological Disorders and Stroke of the National Institutes of Health under Award Numbers K99MH120046 , RF1MH114268 , R01NS088674-S1 , RF1MH114253 , and U01AG050618 . The content of current research is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

## Declaration of competing interest

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

## Acknowledgement

Research reported in this publication was supported by the National Institute of Mental Health and the National Institute of Neurological Disorders and Stroke of the National Institutes of Health under Award Numbers K99MH120046 , RF1MH114268 , R01NS088674-S1 , RF1MH114253 , and U01AG050618 . The content of current research is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

## Appendix A. Supplementary data

The following are the Supplementary data to this article:

- Multimedia component 1

## References

- Transcranial electric and magnetic stimulation: technique and paradigms.in: Handbook of clinical neurology. vol. 116. Elsevier, 2013: 329-342
- Fundamentals of transcranial electric and magnetic stimulation dose: definition, selection, and reporting practices.
*Brain Stimul: Basic Transl Clin Res Neuromodulation.*2012; 5: 435-453 - Electric field calculations in brain stimulation based on finite elements: an optimized processing pipeline for the generation and usage of accurate individual head models.
*Hum Brain Mapp.*2013; 34: 923-935 - Accuracy and reproducibility study of automatic MRI brain tissue segmentation methods.
*Neuroimage.*2010; 51: 1047-1056 - Accuracy and reliability of automated gray matter segmentation pathways on real and simulated structural magnetic resonance images of the human brain.
*PLoS One.*2012; 7e45081 - Evaluation of automated brain MR image segmentation and volumetry methods.
*Hum Brain Mapp.*2009; 30: 1310-1327 - Influences of skull segmentation inaccuracies on EEG source analysis.
*Neuroimage.*2012; 62: 418-431 - Modeling of the human skull in EEG source analysis.
*Hum Brain Mapp.*2011; 32: 1383-1399 - Uncertainty quantification in transcranial magnetic stimulation via high dimensional model representation.
*IEEE Trans Biomed Eng.*2015; 62: 361-372 - Development and characterization of the InVesalius Navigator software for navigated transcranial magnetic stimulation.
*J Neurosci Methods.*2018; 309: 109-120 - Benchmarking transcranial electrical stimulation finite element simulations: a comparison study.
*J Neural Eng.*2019; 16 (2)https://doi.org/10.1088/1741-2552/aafbbd - Navigated transcranial magnetic stimulation.
*Neurophysiol Clin/Clin Neurophysiol.*2010; 40: 7-17 - A posteriori error estimation in finite element analysis. vol. 37. John Wiley & Sons, 2011
- Uncertainty quantification in transcranial magnetic stimulation via high-dimensional model representation.
*IEEE Trans Biomed Eng.*2015; 62: 361-372 - Effects of coil orientation on the electric field induced by TMS over the hand motor area.
*Phys Med Biol.*2013; 59: 203https://doi.org/10.1088/0031-9155/59/1/203 - A principled approach to conductivity uncertainty analysis in electric field calculations.
*Neuroimage.*2018; - The development and modelling of devices and paradigms for transcranial magnetic stimulation.
*Int Rev Psychiatry.*2017; 29: 115-145 - Field modeling for transcranial magnetic stimulation: a useful tool to understand the physiological effects of TMS?.in: 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). 2015: 222-225 (2015)
- Electric field depth–focality tradeoff in transcranial magnetic stimulation: simulation comparison of 50 coil designs.
*Brain Stimul.*2013; 6: 1-13 - Redesigning existing transcranial magnetic stimulation coils to reduce energy: application to low field magnetic stimulation.
*J Neural Eng.*2018; 15https://doi.org/10.1088/1741-2552/aaa505 - VoxHenry: FFT-accelerated inductance extraction for voxelized geometries.
*IEEE Trans Microw Theory Tech.*2018; 66: 1723-1735 - Bioelectromagnetic forward problem: isolated source approach revis (it).
*Phys Med Biol.*2012; 57: 3517 - Linking physics with physiology in TMS: a sphere field model to determine the cortical stimulation site in TMS.
*Neuroimage.*2002; 17: 1117-1130 - The finite element method in electromagnetics.third ed. ed. John Wiley & Sons, 2014
- Fast multigrid-based computation of the induced electric field for transcranial magnetic stimulation.
*Phys Med Biol.*2012; 57: 7753 - Guidelines for limiting exposure to time-varying electric and magnetic fields (1 Hz to 100 kHz).
*Health Phys.*2010; 99: 818-836 - Where does transcranial magnetic stimulation (TMS) stimulate? Modelling of induced field maps for some common cortical and cerebellar targets.
*Med Biol Eng Comput.*2012; 50: 671-681 - Reducing the staircasing error in computational dosimetry of low-frequency electromagnetic fields.
*Phys Med Biol.*2012; 57: N25 - A pipeline for the simulation of transcranial direct current stimulation for realistic human head models using SCIRun/BioMesh3D. Engineering in Medicine and Biology Society (EMBC).in: 2012 annual international conference of the IEEE. 2012: 5486-5489
- Realistic vOlumetric-Approach to Simulate Transcranial Electric Stimulation-ROAST-a fully automated open-source pipeline.
*bioRxiv.*2017; : 217331 - The ICVSIE: a general purpose integral equation method for bio-electromagnetic analysis.
*IEEE (Inst Electr Electron Eng) Trans Biomed Eng.*2018; 65: 565-574 - Comparison of spherical and realistically shaped boundary element head models for transcranial magnetic stimulation navigation.
*Clin Neurophysiol.*2013; 124: 1995-2007 - Conforming discretizations of boundary element solutions to the electroencephalography forward problem.
*Compt Rendus Phys.*2018; 19: 7-25 - A common formalism for the integral formulations of the forward EEG problem.
*IEEE Trans Med Imaging.*2005; 24: 12-28 - The application of electromagnetic theory to electrocardiology: I. Derivation of the integral equations.
*Biophys J.*1967; 7: 443-462 - Software toolkit for fast high-resolution TMS modeling.
*bioRxiv.*2019; : 643346 - A fast algorithm for particle simulations.
*J Comput Phys.*1987; 73: 325-348 - Comparative performance of the finite element method and the boundary element fast multipole method for problems mimicking transcranial magnetic stimulation (TMS).
*J Neural Eng.*2019; 16https://doi.org/10.1088/1741-2552/aafbb9 - A quasi-static boundary element approach with fast multipole acceleration for high-resolution bioelectromagnetic models.
*IEEE Trans Biomed Eng.*2018; 65: 2675-2683https://doi.org/10.1109/TBME.2018.2813261 - On an indirect boundary element method for the anisotropic EEG forward problem.in: Antennas and propagation (EuCAP), 2015 9th European conference on. 2015: 1-3
- Two volume integral equations for the inhomogeneous and anisotropic forward problem in electroencephalography.
*J Comput Phys.*2017; 348: 732-743 - Improved testing of the magnetic-field integral equation.
*IEEE Microw Wirel Compon Lett.*2005; 15: 615-617 - The p-version of the finite element method.
*SIAM J Numer Anal.*1981; 18: 515-545 - Classical electrodynamics.John Wiley & Sons, 2012
- Electric field simulations for transcranial brain stimulation using FEM: an efficient implementation and error analysis.
*J Neural Eng.*2019; https://doi.org/10.1088/1741-2552/ab41ba - Computational software: simple fmm libraries for electrostatics, slow viscous flow, and frequency-domain wave propagation.
*Commun Comput Phys.*2015; 18: 516-528 - The discontinuous galerkin finite element method for solving the MEG and the combined MEG/EEG forward problem.
*Front Neurosci.*2018; 12 - Influence of tissue conductivity anisotropy on EEG/MEG field and return current computation in a realistic head model: a simulation and visualization study using high-resolution finite element modeling.
*Neuroimage.*2006; 30: 813-826 - Gmsh: a 3-D finite element mesh generator with built-in pre-and post-processing facilities.
*Int J Numer Methods Eng.*2009; 79: 1309-1331 - Minimum-energy coils for transcranial magnetic stimulation: application to focal stimulation.
*Brain Stimul.*2015; 8: 124-134 - Design of transcranial magnetic stimulation coils with optimal trade-off between depth, focality, and energy.
*J Neural Eng.*2018; 15https://doi.org/10.1088/1741-2552/aac967 - A numerical algorithm for the construction of efficient quadrature rules in two and higher dimensions.
*Comput Math Appl.*2010; 59: 663-676 - Electrical Io, engineers E. Computational methods for electromagnetics. vol. 2. IEEE press, New York1998
- The influence of sulcus width on simulated electric fields induced by transcranial magnetic stimulation.
*Phys Med Biol.*2013; 58: 4881 - Multigrid finite element methods for electromagnetic field modeling. vol. 28. John Wiley & Sons, 2006
- How much detail is needed in modeling a transcranial magnetic stimulation figure- 8 coil: measurements and brain simulations.
*PLoS One.*2017; 12e0178952 - Low-frequency conductivity tensor imaging of the human Head in vivo using DT-MREIT: first study.
*IEEE Trans Med Imaging.*2018; 37: 966-976 - FFT, FMM, or multigrid? A comparative study of state-of-the-art Poisson solvers for uniform and nonuniform grids in the unit cube.
*SIAM J Sci Comput.*2016; 38: C280-C306 - The superconvergent patch recovery and a posteriori error estimates. Part 1: the recovery technique.
*Int J Numer Methods Eng.*1992; 33: 1331-1364 - Efficient electric field simulations for transcranial brain stimulation.
*bioRxiv.*2019; : 541409 - Where does TMS stimulate the motor cortex? Combining electrophysiological measurements and realistic field estimates to reveal the affected cortex position.
*Cerebr Cortex.*2016; 27: 5083-5094 - Real-time computation of the TMS-induced electric field in a realistic head model.
*bioRxiv.*2019; : 547315

## Article Info

### Publication History

Published online: October 03, 2019

Accepted:
September 29,
2019

Received in revised form:
September 25,
2019

Received:
January 15,
2019

### Identification

### Copyright

© 2019 Elsevier Inc.

### User License

Creative Commons Attribution – NonCommercial – NoDerivs (CC BY-NC-ND 4.0) | How you can reuse

Elsevier's open access license policy

Creative Commons Attribution – NonCommercial – NoDerivs (CC BY-NC-ND 4.0)

## Permitted

### For non-commercial purposes:

- Read, print & download
- Redistribute or republish the final article
- Text & data mine
- Translate the article (private use only, not for distribution)
- Reuse portions or extracts from the article in other works

## Not Permitted

- Sell or re-use for commercial purposes
- Distribute translations or adaptations of the article

Elsevier's open access license policy