Preprocessing General Head Models for BEM-FMM Modeling Pertinent to Brain Stimulation

Introduction: Transcranial magnetic stimulation (TMS) is a major noninvasive neurostimulation method in which a coil placed near the head employs electromagnetic induction to produce electric fields and currents within the brain. To predict the actual site of stimulation, numerical simulation of the electric fields within the head using high-resolution subject-specific head models is required. A TMS modeling software toolkit has been developed based on the boundary element fast multipole method (BEM-FMM), which has several advantages over conventional finite element method (FEM) solvers.Objective: To extend the applicability of the BEM-FMM TMS simulation toolkit to head models whose meshing scheme produces a single mesh for every unique tissue instead of producing a single mesh for every unique tissue/tissue boundary.Method: The MIDA model of the IT’IS Foundation, Switzerland, comprises 115 high-resolution tissue models in the form that the BEM-FMM toolkit is modified to accept. The updated BEM-FMM toolkit is tested using this head model.Results: The BEM-FMM toolkit has been successfully modified to accept head models consisting of one unique mesh per unique tissue while still supporting its initial model format of one unique mesh per boundary between two specific tissues. Performance impacts occur in the preprocessing phase only, meaning that the charge computation method performs equally well regardless of model format.


Introduction
Transcranial magnetic stimulation (TMS) is a noninvasive neurostimulation method wherein a coil placed near the subject's head induces electric currents within the brain [7,10,13]. However, intermediate tissues between the coil and the cortex strongly affect the induced electric field (and thus the induced current). Numerical simulation of the interaction between the primary electric field and tissues of the head is necessary to predict the behavior of the total induced electric field and find the ultimate activation site(s). Further, wide intersubject variations cause the actual fields to deviate strongly from expected fields calculated using a generic head model. To minimize deviation between the simulated and actual fields, the simulated fields must be calculated using an accurate, high-resolution, subject-specific head model.
The TMS toolkit (complete computational code and supporting documentation) available for academic use at the Dropbox repository [2] is one such TMS simulator, which utilizes the boundary element fast multipole method (BEM-FMM) described in [4,9]. The toolkit is written for MATLAB R2019a and has dependencies on the Image Processing Toolbox, Partial Differential Equations Toolbox/Antenna Toolbox, and Statistics and Machine Learning Toolbox. Its core FMM method is that of [3], included with permission in the redistributable software package.
This toolkit enables users to simulate TMS behavior using predefined or custom coil CAD models and subject-specific head models. These head models consist of a set of nested 3D triangular meshes, where each mesh marks the boundary between two tissues with different electrical properties (e.g., one mesh follows the skin/skull boundary, and another mesh follows the gray matter (GM)/white matter (WM) boundary) [5]. Because the BEM-FMM algorithm operates directly in terms of induced charges on these interfaces [8], it is robust against several common mesh defects that would hinder conventional volumetric finite element method (FEM) simulations, including intersecting meshes. The BEM-FMM algorithm further supports computation of the net electric field at locations arbitrarily close to tissue interfaces, where FEM routines cannot provide field resolution that exceeds the resolution of the underlying volumetric mesh.
Despite the BEM-FMM algorithm's robustness against common mesh defects, the current implementation of the software toolkit is applicable only to one specific meshing scheme: one in which each mesh represents a boundary between exactly two tissues. This is the standard output format of the SimNIBS v2.1 pipeline [11,12,[14][15][16][17] in particular, and it produces meshes that are layered one inside the other. For example, the boundary between the skull and cerebrospinal fluid (CSF) completely surrounds and encloses the boundary between the CSF and gray matter (GM), which in turn completely surrounds and encloses the boundary between gray matter and white matter (WM). The goal of this exercise was to add support for a second meshing scheme, in which each mesh represents the entire outer boundary of a single tissue. One model that employs this meshing scheme is the MIDA model, produced by the IT'IS Foundation [6] The MIDA head model comprises 115 CAD tissue models with more than 11 M triangular facets total. The model was produced from scans of a healthy 29-year-old female volunteer. Data was compiled from several medical imaging methods, including MRI, MRA, and DTI. These diverse imaging methods ensured that high-contrast images of most tissues existed in at least one of the image sets and image resolution approached 500 μm. Special care was taken to obtain high-contrast images of nerve tissue and vasculature. The entire data set was segmented independently by three experts using both manual and automated segmentation techniques, and their individual segmentations were combined to produce a highly accurate final segmentation. A triangulation algorithm was then applied to the resulting voxel model to extract triangular mesh surfaces for every tissue [6].

Mesh Preprocessing
Because the MIDA model's tissue meshes are not nested in general (i.e., a given MIDA tissue mesh explicitly segments every boundary between that mesh and any other tissue), adjacent tissue meshes each contain their own copies of the facets that form the border between them. When two adjacent tissue meshes are loaded simultaneously, their shared border comprises two sets of coincident facets, one set contributed by each tissue mesh. These coincident facets necessarily share coincident centroids, which in turn create singularities that invalidate simulation results. Figure 11 depicts this case for three hypothetical meshes, Object 1, Object 2, and Object 3. Though Object 1 and Object 2 both segment their shared boundary, Object 3 does not explicitly segment its boundary with Objects 1 and 2 for  simplicity. Every mesh is assigned a default exterior conductivity derived from manual inspection of the surrounding tissues. The function knnsearch of MATLAB's Statistics and Machine Learning Toolbox is first used to pair facets that have coincident centroids. One facet of each pair is designated as the facet to be kept, and the other is designated as the facet to be deleted. The outer conductivity of the facet to be kept is set equal to the inner conductivity of the facet to be deleted, and associated contrast information is updated for the facet to be kept. After this process has been completed for all coincident facet pairs, all information related to the facets to be deleted (e.g., centroid, area, connectivity) is removed, and any now-unreferenced vertices are cleared from the list of vertices (and face connectivity information is updated as appropriate).

Results
The final test setup modeled a TMS configuration intended to target the motor hand area of the precentral gyrus (the hand knob area, [18]) of the MIDA model. The coil model used was a generic figure-eight coil with circular cross-sectional wire, as shown in Fig. 12. The coil was approximated by 16,000 elementary current segments driven by time-varying current dI dt ¼ 9:4e7 Amperes= sec . Preprocessing of the MIDA model for simulation using the BEM-FMM algorithm took approximately 525 s in total. Of those 525 s, 138 s were required to resolve coincident facets. Of the original 11 M facets, approximately 5.4 M were removed, and approximately 5.6 M remain. Table 1 lists the times associated with each preprocessing step.
The coil model was positioned above the head model according to four simple geometric rules: 1. The coil's centerline passes through a selected point on the hand knob area. 2. The coil's centerline is perpendicular to the skin surface. 3. The distance from the coil to the skin surface along the coil's centerline is 10 mm. 4. The dominant field direction (the y-axis of the coil coordinate system) is approximately perpendicular to the gyral crown and associated sulcal walls of the precentral gyrus pattern at the target point. Fig. 12 The coil model employed for this test Figure 13 shows the BEM-FMM convergence curve for this test setup after 100 GMRES iterations, and Fig. 14 shows the convergence curve for the first 15 iterations. The relative residual falls well below the threshold 10 À3 within 15 iterations, indicating that 15 iterations produce results within an acceptable error margin. The test was run on a 32-core Intel® Xeon® E5-2683 v4 CPU operating at 2.1 GHz with 256 GB RAM. On this machine, the total computational time required with 5.6 M facets for 15 GMRES iterations was 373 s.   Figure 15 shows simulation results for the electric field at the gray matter/CSF interface (a, c) and the white matter/gray matter interface (b, d). Figure 15 (a, b) shows a heat map of the magnitude of the electric field at the respective surfaces, scaled in V/m. Figure 15 (c, d) shows a focality estimate of the total electric field. In these figures, small blue balls are drawn at every facet for which the total field magnitude is within the range 80% to 100% of the maximum field magnitude observed for that particular surface. We see that the naïve geometric coil positioning rules barely stimulate the desired region at all and instead produce local maxima at distant sulci rather than the targeted motor hand area. This further reinforces the necessity of subject-specific head modeling for TMS applications. Figures 16, 18 and 20 depict cross sections of the tissue meshes coregistered with T1 MRI data for the MIDA subject. The planes of these cross sections pass through the point on the white matter surface where the maximum E-field magnitude occurs. Pink spheres are drawn at the center of every WM facet that experiences a field with magnitude within 80-100% of the maximum E-field magnitude on this surface. Figures 17 and 19 show contour plots of the electric field magnitude in the immediate vicinity (i.e., AE10 mm) of the maximum field locations of Figs. 16 and 18, respectively.

Conclusion
The BEM-FMM TMS modeling toolkit has been made compatible with a previously unsupported head mesh scheme, in which each mesh corresponds to one tissue's entire outer surface. The modifications were tested using the MIDA head model, which employs the newly supported mesh scheme. Simulation was executed successfully with the MIDA model, achieving convergence within 15 GMRES iterations. The performance penalty associated with the new mesh format occurs solely in the preprocessing stage, and there is little to no effect on field calculation performance. The toolkit is now applicable to a wider range of head models and is more robust against models whose meshes have coincident facets in general.