If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Department of Psychiatry and Behavioral Sciences, Duke University, Durham, NC, 27710, USADepartment of Electrical and Computer Engineering, Duke University, Durham, NC, 27708, USADepartment of Neurosurgery, Duke University, Durham, NC, 27710, USADepartment of Biomedical Engineering, Duke University, Durham, NC, 27708, USA
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.
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 [
]. 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 [
Electric field calculations in brain stimulation based on finite elements: an optimized processing pipeline for the generation and usage of accurate individual head models.
]. 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 [
]. 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 [
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 [
]). 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 [
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 [
]. 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 [
]. 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 [
]. 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 [
] 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’ [
]. 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 [
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
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 [
Comparative performance of the finite element method and the boundary element fast multipole method for problems mimicking transcranial magnetic stimulation (TMS).
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) [
]. 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.
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
]. 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,
(1)
where , , and are the conductivity, scalar potential, and primary E-field at position , 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 [
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
] 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 [
]. 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 [
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 [
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.
] 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 [
]. 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 [
], 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 [
Fig. 1TMS E-field simulation setups. (A) The coil was placed directly above the apex of the spherical head model. (B) The coil was placed above the primary motor cortex hand-knob representation of the MRI-derived head model.
]. (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 [
]. 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 [
]. 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.
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 [
]. 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 [
] 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 [
]. 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. [
]. First, both circular loops of the figure- 8 coil were divided into rings of uniform width. Then, to obtain most uniform set of dipoles, each ring was discretized with equispaced dipoles (where is the inner radius of the ring and — its outer radius). Lastly, 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 L2 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 L2-norm error,
(2)
where denotes magnitude, and are the reference and approximate solution, respectively, and is the domain of integration. For head simulations the integration was done over the GM/WM matter region (i.e., ) and for coil current error—in the coil winding wire (i.e., ). The L2 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 L2 norm error was estimated by sampling the E-field on a grid with points 2 mm apart. For the MRI-derived head, the brain L2 norm error was estimated by applying a 2nd order accurate integration rule [
The second metric is the pointwise magnitude error
(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” L2 errors are shown in Fig. 2. For the homogenous sphere, the maximum pointwise and L2 errors are similar. FDM, 0th 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 L2 error by a factor of 6–13 compared to 1st order FEM. Furthermore, increasing the order of FEM or BEM decreases the L2 error by 1–2 orders of magnitude.
Fig. 2Maximum pointwise error of ‘cortical’ E-field in (A) homogenous and (B) inhomogeneous spherical head model. Global L2 norm error of ‘cortical’ E-field in (C) homogenous and (D) inhomogeneous spherical head model.
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 0th 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 0th to 1st, respectively, disappear. SPR-FEM again performs between 1st and 2nd order FEM, decreasing the L2 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 0th 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 0th for BEM, it is likely necessary to use higher-order mesh element shape functions [
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 0th 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 () is given below each plot. The reference has 100 V/m. For FDM, 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 [
]. The higher-order FEM methods (2nd and 3rd order) and 0th order BEM method predict an 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.
Fig. 3E-field on cortical surface. Results are arranged in columns by method. Rows are arranged with increasing number of mesh refinements top to bottom. Maximum simulated E-field strength is given in V/m in the right bottom corner of each image. The reference solution has a maximum E-field strength of 100 V/m.
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.
Fig. 4Pointwise error of cortical E-field. Results are arranged in columns by method. Rows are arranged with increasing number of mesh refinements top to bottom.
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 [
]. SPR decreases the L2 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 [
]. For mesh 0, SPR-FEM has L2 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 L2 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 0st 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, 0th order BEM, and 3rd order FEM are marginally different in error.
Fig. 5Global L2 norm error of cortical E-field in MRI-derived head model.
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 0th 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. [
]. More detailed discussion of CPU runtime, memory usage and possible speed-ups is given in the Supplement.
Fig. 6CPU runtimes of simulations with MRI-derived head model. (Note codes have only been optimized for accuracy. Speeds are reported here only for relative comparison and can be significantly improved upon in absolute terms.)
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 [
]). 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 L2 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 L2 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 [
] 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.
Fig. 7Primary E-field (A) pointwise error and (B) L2 norm error using current discretization and equivalent magnetic dipole approach.
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.
Fig. 8Top view of each coil model along with current distribution at dc, 3 KHz, and 6 KHz in the wire-cross section denoted by red dashed line. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)
The L2 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 L2 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 [
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 [
], 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 [
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% L2 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 [
]. 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 [
]. 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 0th 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 [
]. Because of the complicated cortical geometry, a dense mesh with average edge length ≤1.5 mm is required to sufficiently resolve the geometric features [
]. The use of a 2nd order FEM or 0th 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 [
Where does TMS stimulate the motor cortex? Combining electrophysiological measurements and realistic field estimates to reveal the affected cortex position.
], 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 0th 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:
Electric field calculations in brain stimulation based on finite elements: an optimized processing pipeline for the generation and usage of accurate individual head models.
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
Comparative performance of the finite element method and the boundary element fast multipole method for problems mimicking transcranial magnetic stimulation (TMS).
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.
Where does TMS stimulate the motor cortex? Combining electrophysiological measurements and realistic field estimates to reveal the affected cortex position.
30378-X&id=gr1.jpg){kind=link}
30378-X&id=gr2.jpg){kind=link}
30378-X&id=gr3.jpg){kind=link}
30378-X&id=gr4.jpg){kind=link}
30378-X&id=gr5.jpg){kind=link}
30378-X&id=gr6.jpg){kind=link}
30378-X&id=gr7.jpg){kind=link}
30378-X&id=gr8.jpg){kind=link}