A fast direct solver for surface-based whole-head modeling of transcranial magnetic stimulation

Background: When modeling transcranial magnetic stimulation (TMS) in the brain, a fast and accurate electric field solver can support interactive neuronavigation tasks as well as comprehensive biophysical modeling. Objective: We formulate, test, and disseminate a direct (i.e., non-iterative) TMS solver that can accurately determine global TMS fields for any coil type everywhere in a high-resolution MRI-based surface model with ~200,000 or more arbitrarily selected observation points within approximately 5 sec, with the solution time itself of 3 sec. Method: The solver is based on the boundary element fast multipole method (BEM-FMM), which incorporates the latest mathematical advancement in the theory of fast multipole methods – an FMM-based LU decomposition. This decomposition is specific to the head model and needs to be computed only once per subject. Moreover, the solver offers unlimited spatial numerical resolution. Results: Despite the fast execution times, the present direct solution is numerically accurate for the default model resolution. In contrast, the widely used brain modeling software SimNIBS employs a first-order finite element method that necessitates additional mesh refinement, resulting in increased computational cost. However, excellent agreement between the two methods is observed for various practical test cases following mesh refinement, including a biophysical modeling task. Conclusion: The method can be readily applied to a wide range of TMS analyses involving multiple coil positions and orientations, including image-guided neuronavigation. It can even accommodate continuous variations in coil geometry, such as flexible H-type TMS coils. The FMM-LU direct solver is freely available to academic users.


Introduction 72
Transcranial magnetic stimulation (TMS) is a non-invasive technique used to stimulate specific 73 regions of the brain. For TMS to be a more effective and personalized therapeutic tool, it is crucial 74 to accurately simulate the induced electric field (E-field) distribution in the brain. For interactive E-75 field based TMS neuronavigation, the E-field should be updated and displayed fast enough to 76 provide the operator with a sense of real-time interaction. 77 There are several time scales to consider. First, most real-time visualization applications 78 have a minimum update rate of around 15-30 frames per second (fps) to have a smooth and 79 continuous experience; this allows for approximately 30-60 ms for coil position update, E-field 80 solution update, and image rendering. The second time scale of relevance relates to the pace at 81 which a human or a robot [1],[2] could move a TMS coil from one position to another. Given an 82 estimate for sequential robotic targeting speed of ~1cm/s [1], and if we consider coil position 83 change of a few millimeters to produce a meaningful change in the E-field, then the solution can 84 be updated on the time scale on the order of 100 ms. The final time scale is that of the TMS 85 interstimulus interval, if the E-field solution were to update synchronous with TMS pulse delivery. 86 For conventional motor mapping with single pulses, an approximately ~5 s interstimulus interval 87 is typical to avoid inhibitory effects (cf. [

3],[4]). 88
Over the past two years, significant progress has been made with regard to excellent real-89 time high-resolution solvers for TMS field modeling [5],[6], [7], which execute in approximately 20-90 50 ms. These simulation speeds are close to the standard screen refresh rate and should be very 91 useful in practice for continuous smooth field visualization. In general, these solvers do not involve 92 accelerating the main numerical algorithm itself but rather apply highly efficient and "smart" 93 interpolation methods to thousands of precomputed solution sets. 94 This study aims to introduce a direct accurate TMS solver that can output TMS fields from 95 any coil and anywhere in a high-resolution head model within ~3-5 sec. This solver is based on 96 the boundary element fast multipole method or BEM-FMM [8], [9], [10], [11] and the most recent 97 mathematical advance in the theory of fast multipole methodsan FMM-based LU factorization 98 or FMM-LU [12]. It does not use any interpolation routines or precomputed solutions, but does 99 require initial LU decomposition of a head-model matrix, which takes approximately 40 min. 100 One motivation for this study is a modern precise motor mapping protocol [13], which 101 would benefit from for fast accurate field calculations while the participant is still in the hospital or 102 facility so that a follow up TMS session could be avoided. In addition, the mapping process could 103 be further optimized through adaptive algorithms that will utilize a sequential approach for 104 estimating the "hot spot" [14]. On the other hand, the fast solver can be readily utilized to pre-105 compute a large number of E-field solutions covering "all possible coil positions/orientations" of 106 interest that can be subsequently used to stimulate a desired cortical target derived from 107 anatomical of functional connectivity data (see, e.g., [15]). The precalculated solutions could also 108 be used to guide the operator to the right target in an interactive neuronavigation setting. 109 In contrast to the finite element method, BEM formulations are restricted to the surface of 110 compartment domains, thereby reducing the dimensionality of the problem for brain modeling in 111 general. There has been a resurgence of interest in these surface-based methods due to the 112 availability of fast algorithms such as fast multipole methods (FMMs) [ [21]. 121 There is an obvious task where iterative solvers are not satisfactory. This includes near 122 real-time E-field predictions for TMS brain mapping or stimulation planning with a large number 123 of possible coil positions, orientations, or even coil geometry changes. Each such setup gives rise 124 to a unique right-hand side. Here, direct solvers would be preferred to solve the same system 125 matrix with multiple right-hand sides. Along with this, fast direct solvers could be useful when 126 exploring low-rank perturbations of the head geometry (and/or tissue conductivities) and, hence, 127 the system matrix. Updating the solution in such cases requires only a few applications of −1 or 128 a fast update of the inverse itself [22], [23]. 129 In the last few years, several algorithmic ideas have emerged which permit the 130 construction of a compressed approximation of −1 at a cost of the order ( ) or ( ),

131
for modest . To construct the first fast direct TMS solver, we apply one such scheme [12], which 132 is referred to as the FMM-LU method [12]. It uses FMM-type hierarchical compression strategies 133 to rapidly compute an LU-factorization of the large system matrix.

Charge-based BEM in the succinct form 137
Induced electric charges with a surface charge density ( ) will reside on every tissue conductivity 138 interface once an external electromagnetic stimulus or a primary E-field ( ) of a TMS coil is 139 applied. These induced surface charges will alter the primary field to fulfill the law of current 140 conservation across the boundaries. The E-field generated by all surface charges anywhere in 141 space except the charged interfaces themselves is governed by Coulomb's law. The total E-field 142 ( ) becomes the sum of the primary field and the secondary charge field i.e., 143 where 0 is dielectric permittivity of vacuum (a normalization constant). The E-field in Eq. (1) is 144 discontinuous at the interfaces. When approaching a charged interface with a certain normal 145 vector and assigning index into the medium from which is pointing and index out of the 146 medium toward which is pointing, the E-field close to the boundary is given by two limiting 147 Here, , are the conductivities just inside and outside with respect to the direction of the 159 normal vector. After combining similar terms, a Fredholm equation of the second kind is obtained 160 where is the electric conductivity contrast = − + . 161

Iterative (BEM-FMM) versus direct (FMM LU) solution of Eq. (4) 162
Assuming that the charge density ( ) has a constant value for every triangular facet with 163 area , Eq. (4) is discretized in standard matrix form ( is the primary field at the m-th facet) 164 where is Kronecker delta. Matrix is well conditioned. Therefore, an iterative solution of Eq. 165 iterations. FMM is used to compute the matrix-vector product in a "matrix free" fashion so that 167 matrix is never computed/stored [30], [26]. This solution was used previously [8], [9], [10], [11]. 168 Another way to solve Eq. (5) is to apply FMM to a direct LU-decomposition of matrix . In 169 other words, a compressed approximation of −1 is found directly so that a direct solution = 170 −1 is attempted, without using an iterative method and for multiple right -hand sides, . While 171 matrix is determined by the head model, the right-hand sides describe coil fields. This is the 172 compression algorithm of Ref.
[12], which is applied in this study. Indeed, the compressed inverse 173 differs from the iterative solution. We will show that this difference is vanishingly small. 174 the skin effect included) and the coil fields are also computed via the fast multipole method [34]. 184

Human models and coil models used to test the algorithm
Finally, we test a healthy subject scanned at Max Planck Inst. for Human Cong. & Brain 185 Sciences Leipzig, Germany with the MagVenture CB65 coil and the same left 1 HAND target. 186 However, this model is manually refined in the region of interest (ROI) so that its overall size is 187 1,8 M facets. In all cases, the default SimNIBS conductivity parameters [24] have been used. 188

Testing method 189
Both accuracy and speed are tested. Two error types are considered: the relative error in the 190 vector E-field, , and the error in the magnitude of the E-field, 191 over an arbitrary domain of interest where ‖•‖ is a 2-norm for a vector or scalar field. Index 192 denotes the total E-field in the head which is a sum of the primary field and the secondary field.

Method accuracy versus SimNIBS FEM accuracy 198
For the default Ernie model example [24] of the FEM SimNIBS software [25], Table 1 gives the 199 field error values for a mid-surface between the gray and white matter with ~200,000 observation FMM solution, we also provide results for the uniformly refined models. 203 The iterative BEM-FMM solution in Table 1 Table 2 gives the field error values for the mid-surface between gray and white matter. Two

Method speed 232
The method speed was tested using four nearly identical workstations with Intel(R) Xeon(R) Gold 233 6348 CPU @ 2.60GHz, 512 GB RAM, 56 cores, OS Windows Server 2022 Standard and the 234 common MATLAB 2021b platform. Fig. 1 outlines the corresponding execution times including 235 precomputations. In Fig. 1, initial

Example #1. Cortical and subcortical fields computed in 4.7 sec (MRiB91 coil) 249
As the first computational example, Fig. 2 shows the simulation output for Connectome Young 250 Adult subject 120111 including: a) the total E-field just inside the gray matter interface; b) the total 251 E-field just outside the white matter interface; c) the total E-field at the cortical midsurface and; d) 252 the total E-field in a transverse plane beneath the coil. The MRiB91 coil of MagVenture is driven 253 with / = 9.4 7 A/s. The computational sequence for any type of output from Fig. 2 takes  254 approximately 4.7 seconds to run, including graphical rendering in MATLAB. 255

Example #2. Cortical and subcortical fields computed in 4.8 sec (authentic H1 coil) 256
As a second example, Fig. 3 shows a simulation output for the same Connectome subject 120111 257 but when an H1 coil of BrainsWay is used. In this case, the coil is made flexible, and it changes 258 its shape when continuously aligned with the patient's head. Therefore, step I in Fig. 1 cannot be  259 replaced by an interpolation of the precomputed coil field. The computational sequence for any 260 type of output from Fig. 3 Table 2 indicates that the accuracy of FMM LU is non-distinguishable from the accuracy of BEM-280 FMM given the FMM precision level of 1e-4 or better used for constructing the compressed 281 inverse. In this case, the factorization data stored in MATLAB workspace have the size of 282 approximately 50 Gigabytes; their creation requires approximately 70 min. The direct solution -283 step II in Fig. 1executes in 3.0 s. The FMM precision level used for the field computations -284 steps I and III in Fig. 1 may be as low as 1e-1 or 1e-2 without deteriorating the solution accuracy. 285

Iterative and direct FMM solutions produce nearly identical results 279
When the FMM precision level is reduced to 1e-3, the factorization data will have the size 286 of 30 Gigabytes, their generation requires approximately 40 min, and the direct solution (step II in 287 Fig. 1) executes in 2.4 s. This option could likely be preferred since its accuracy reported in the 288 last four rows of Table 2 is approximately 0.4% on average for the critical ROI fields. 289 As Table 2

Solid agreement with SimNIBS is achieved for a refined head model 295
The lower 1 st -order FEM accuracy from Table 1 can be improved using selective mesh refinement 296 in the ROI domain as illustrated in Fig. 4a,b. Here, the default headreco segmentation [24] of 297 the healthy subject was manually refined for cerebrospinal fluid, gray matter, and white matter 298 (both surface and volume meshes were refined) as shown in Fig. 4b. One thousand different coil 299 positions (using a MagVenture CB65 coil) have been tested while scanning the ROI. 300 Solutions using the refined meshes from SimNIBS and BEM-FMM were compared with 301 each other as shown in Fig. 4c,d, respectively. The average vector field and magnitude field 302 differences from Eq. (6) within the ROI in Fig. 4c,d now attain the sub-percent values, which 303 confirms the convergence of both methods to the same result. Due to the differences in the 2D 304 vs. 3D numerical formulations of the problems for the BEM-FMM and SIMNIBS, respectively, 305 some discrepancies remain close to the conductivity boundaries where the normal electric field is 306 discontinuous. The practical implications of these differences require further investigation. 307

Solid agreement with SimNIBS is achieved for activating thresholds of intracortical 308 neural cells for a refined head model 309
Three intracortical neural cells (#4, 6, and 9) from a multi-scale toolbox Neuron Modeling for TMS 310 (NeMo-TMS) [37] were placed 1 mm below the grey matter surface of the refined Ernie model 311 (0.4 nodes per mm 2 ) on the coil axis as shown in Fig. 5a and their activating thresholds were 312 computed as a function of the E-field intensity, / , of a MagStim 70 mm coil (the SimNIBS coil 313 model is used) for the default NeMo-TMS pulse form. Fig. 5b)-e) below compares the activating 314 thresholds computed with BEM-FMM and SimNIBS, respectively, using the corresponding quasi-315 potentials. It also illustrates the action potentials at the neuron activation for cells 6 and 9, 316 13 respectively, at 0.7 ms. The activating thresholds differ by 3.2% (cell 4), 2.6% (cell 6), and 4.2% 317 (cell 9) with BEM-FMM always predicting slightly lower threshold values while the action potentials 318 (membrane voltages) at 0.7 ms are hardly distinguishable. For the non-refined model, the 319 agreement is worse (~20% on average). The results indicate that 'inside' the gray matter the BEM-320 FMM and SIMNIBS based modeling of the neuronal activation thresholds show robust agreement. 321 On the other hand, axonal activation mechanisms may be influenced by the tissue heterogeneity 322 (that causes charge accumulation on the conductivity boundaries) that may require more detailed

.4. Major limitations of the present direct solution 338
The present method has been tested on four multicore (32 to 56 cores) 2.6 GHz workstations.