Real-time computation of brain E-field for enhanced transcranial magnetic stimulation neuronavigation and optimization

Transcranial magnetic stimulation (TMS) (Barker et al., 1985; Paulus et al., 2013) is a non-invasive brain stimulation approach that is widely used in neuroscience research to study brain function and has been approved by the Food and Drug Administration in the treatment of depression, obsessive-compulsive disorder, migraine, and smoking (Carmi et al., 2019; Cohen et al., 2022; Couturier, 2005; De Ridder et al., 2007; Hoffman et al., 2005; Kaptsan et al., 2003; Lan et al., 2017; Lipton et al., 2010; Nahas et al., 2003; O’Reardon et al., 2007; Pascual-Leone et al., 1996). TMS uses coils driven by low-frequency current pulses to stimulate targeted brain regions. Computational electric field (E-field) dosimetry can be used to 1) quantify the E-field strength and distribution to determine brain regions affected by TMS, 2) identify coil placements and orientation that would maximize the E-field induced at a prespecified target, and 3) design coils with a desired E-field profile (Beynel et al., 2020; Dannhauer et al., 2022; Gomez et al., 2018, 2021; Laakso & Hirata, 2012; D. Wang et al., 2023; Weise et al., 2023). All of these applications require repeated execution of E-field solvers to determine the E-field. Correspondingly, there is ongoing interest in fast and accurate E-field solvers for TMS (Daneshzand et al., 2021; Gomez et al., 2021; Stenroos & Koponen, 2019). In particular, there is an interest in incorporating accurate and precise E-field information in neuronavigation systems that use subject Magnetic Resonance Imaging (MRI) along with cameras to provide practitioners with precise information about the coil location relative to the head. Appending E-field information to neuronavigation tools would allow for on-the-fly E-field informed coil placement reconfiguration and intensity dosing to target multiple or even changing cortical targets during a single TMS session (Daneshzand et al., 2023; Nieminen et al., 2022). The integration of E-field information into neuronavigation requires a solver that can accurately determine the E-field in real-time.

Several approaches have been proposed for determining the E-field in real-time. For example, the most widely used approach locally fits a sphere to the head to estimate the E-field under the coil rapidly (Nummenmaa et al., 2013; Sollmann et al., 2016). However, this solver does not account for gyral folds of the cerebrospinal fluid and gray matter boundary surface, consequently, yielding less accurate local E-field estimates (Nummenmaa et al., 2013). Furthermore, deep learning approaches have been introduced, which have been promising. However, at this point, they still result in a high relative residual error of about 6% even in a small target region (Yokota et al., 2019) and 18% in the whole brain in (Xu et al., 2021). More recently, boundary element method (BEM)-based near real-time solvers that can account for local brain anatomy have been developed (Daneshzand et al., 2021; Stenroos & Koponen, 2019).

The elegant method of (Stenroos & Koponen, 2019) uses efficient quadratures for the coil and tissue boundary sources enabling real-time E-field solutions using GPUs. For example, the authors developed coil dipole models that achieved the same accuracy as those with thousands of dipoles but with only 42 dipoles per layer. Additionally, they created BEM meshes that properly accounted for gyral geometry while only requiring 21,052 nodes. These efficient quadratures coupled with a highly optimized hardware implementation of BEM enabled the approximation of the E-field in real time. Specifically, they use pre-computed boundary potentials on a mesh to compute the TMS-induced E-field in a cortical region of interest (ROI) in real-time via reciprocity. Using their method, the computation scales as the product of ‘the number of cortical E-field evaluation points’ times ‘the number of vertices of the surface meshes’ plus ‘the number of coil quadrature nodes’. As a result, a trade-off between the number of mesh nodes and the number of cortical evaluation points must be made to keep the computation low. These computations are done rapidly by leveraging massively parallel architectures like GPUs, enabling the mitigation of asymptotic computational costs for many practical scenarios.

To maintain fast computation, there must be a trade-off between mesh resolution and the number of ROI points. For example, to achieve 36 ms on a GPU they used a mesh consisting of 21,052 nodes and 20,324 cortical ROI locations. Recently, the authors introduced an improved GPU-accelerated implementation of the method (Stenroos et al., 2023), where it took 20.41 ms for a 42 dipole coil model and 22.73 ms with a 15,000 dipole coil model to compute the E-field for a single coil placement over the same head model. The Ernie mesh, a typical head mesh generated by SimNIBS software v3.2 (Thielscher et al., 2015), has more than 10 times the number of cortical surface nodes (216,130 nodes). Such a computation cannot be achieved using current hardware in real-time via the above method. Furthermore, newer more detailed head models generated by using SimNIBS v4.0 contain an even higher number of nodes. As a result, we aim to extend real-time solvers to analyze E-fields while using these denser and more detailed meshes that are typically used in the non-invasive brain stimulation community for modeling.

To overcome these challenges, Daneshzand et al. proposed the Magnetic Stimulation Profile (MSP) approach that approximates the TMS coil-induced E-fields in a cortical ROI as a linear combination of dipole-induced E-fields (Daneshzand et al., 2021). Dipole-induced E-fields are precomputed in a process that can take 5 to 18 hours depending on desired accuracy and mesh resolution. Then, these precomputed E-fields are used in real-time to determine the TMS coil-induced E-fields. To find E-field expansion coefficients, the method employs a least squares to match the primary E-field of the coil and the primary E-fields induced by a linear combination of dipoles. This method requires $∼$3000 dipole E-fields, requiring $∼32 GBs$ of memory and about 0.37 s using a Xeon E5-2360 CPU for 120,000 cortical triangles to match the full BEM solution with an average relative residual error of about 10% and a 5% error in predicting the peak E-field (Daneshzand et al., 2021). The memory and CPU time of this method scales as the number of dipole E-fields times the number of E-field evaluation points. As such, to lower the memory requirements and achieve real-time performance, trade-offs between accuracy and cortical E-field samples must be made.

In this paper, similar to the MSP approach, we approximate the E-field as a linear combination of precomputed basis E-fields. However, we use a novel method for determining an E-field basis instead of the dipole E-field basis. This results in a more accurate E-field solution for a fixed number of precomputed basis E-fields. To find an efficient basis of E-field modes on the cortex due to any source outside of the head, we apply an approach similar to the probabilistic matrix decomposition (PMD) based approach (Hasan et al., 2023). By using modes instead of dipole E-fields, we reduce the required number of modes more than 10-fold relative to the MSP approach. For example, to reach a 10% error, we require 110 modes relative to the 3000 modes required using the MSP approach. Furthermore, our approach requires under 400 basis modes to estimate the TMS-induced E-field with an error lower than 3%. Additionally, we introduce a novel method to determine the expansion coefficients from the primary E-fields and magnetic field (H-field) on a fictitious surface engulfing the head that minimizes the energy of the error of the predicted field. Using this method, we can estimate the E-field to 2% error within 4 ms from the primary E-field and H-field for 216,000 cortical surface targets. Finally, unlike the MSP, the expansion coefficients are found analytically using reciprocity and Huygens’s principle formalism to determine the coefficients without the need for numerical techniques requiring regularization to prevent over-fitting.

The notations used in this article are summarized in Appendix Table 1. This section describes the procedure for real-time determination of an approximate expansion for the E-field induced by the coil in the brain during TMS

Here, $r$ denotes a Cartesian location, $Ω$ denotes the brain region, $M(i)(r;t)=I′(t)M(i)(r)$ and $a(i)$ are one of the $Nm$ mode functions and expansion coefficients, respectively, $I′(t)$ is time-derivative of the driving current pulse waveform. Here, we assume that $I′(t)$ is normalized to have a maximum time derivative of one. Correspondingly, this results in $ETMS(Nm)(r)$ being the peak E-field observed. Note that above and in the rest of this paper, we assume the TMS-induced E-fields are quasi-static and this is standard for TMS E-field modeling (Daneshzand et al., 2021; Gomez et al., 2020b, 2021; Plonsey, 1972; Thielscher et al., 2015; D. Wang et al., 2023; B. Wang et al., 2024). As a result, we assume that the TMS coil driving current is separable $JTMS(r;t)=I(t)JTMS(r)$ and the H-fields and E-fields are proportional to the coil driving current and its temporal derivative, respectively (more details can be found in section 6.1 and 6.2 of the Supplementary File).

In Section 2.2, we describe a procedure for finding the spatial variation of the orthonormal mode functions $M(i)(r)$ (i.e., $〈M(i),M(j)〉=δi,j$⁠, where $〈f,g〉=∫Ωf(r)⋅g(r)dr$ and $δi,j$ is the Kronecker delta function) that can efficiently represent the E-fields induced in the brain by any TMS coil. Once these functions are determined, the coefficients $a(i)$ are chosen to minimize the $L2$ error of the expansion (i.e., $argmina∈ℝNm‖ETMS−ETMS(Nm)‖$⁠, where $ETMS$ is the E-field modeled using the in-house built field solver, $‖f‖=‖f,f‖$⁠, and $a=(a(1),a(2),...,a(Nm)))$⁠. Since the mode functions are orthonormal, the coefficients are $a(i)=〈M(i),ETMS〉$⁠. The above expression for determining each $a(i)$ requires a priori knowledge of the TMS-induced E-field in the brain. As such, it is not amenable to real-time solvers. Alternatively, in Section 2.3, we use reciprocity and Huygens’s equivalence principles (Balanis, 2012) to show that

Here, $S$ denotes the Huygens’s surface separating the head and the coil, $JS(i)(r;t)=I(t)JS(i)(r)$ and $KS(i)(r;t)=I′(t)KS(i)(r)$ are fictitious equivalent electric and magnetic current densities associated with the $ith$ mode function, respectively, and $ETMSP(r;t)=I′(t)ETMSP(r)$ and $HTMSP(r;t)=I(t)HTMSP(r)$ are the primary E-field and H-field, respectively. Equation (2) only requires primary fields on Huygens’s surface and here it is used to determine expansion coefficients in real-time. In section 2.5, we provide a method for rapidly determining $ETMSP(r)$ and $HTMSP(r)$ on the Huygens’s surface. Section 2.6 summarizes the pre-processing and real-time stages for an easy implementation. Finally, sections 2.7 and 2.8 describe the head and coil models for validating the algorithm and error quantization metrics, respectively.

The mode functions $M(i)(r;t)$⁠, where $i=1,2,...,Nm$⁠, are chosen as an orthonormal basis to E-fields induced in the brain by $Nm$ magnetic surface current density distributions residing on a fictitious surface 1 mm directly outside of the scalp (Fig. 1). The reason we chose magnetic surface current distributions is that they do not need to be divergence-free, in fact, they are oriented randomly on the fictitious surface at each point in space with a normal distribution. The $Nm$ current density distributions are chosen as distinct realizations of Gaussian white noise $W(i)(r;t)=I′(t)W(i)(r)$⁠. The $Nm$ Gaussian white noise magnetic current density realizations are continuous analogous to the Gaussian white noise vectors that we previously successfully used to compress matrices of TMS-induced brain E-field samples for many coil placements and brain locations (Hasan et al., 2023).

To generate each Gaussian white noise distributed current density, we first create a triangle mesh by extruding the nodes of the scalp surface mesh 1 mm normally outward. Then, we approximate the current as a piece-wise constant within each of the $Nd$ triangles of the triangle mesh. The value of the current density is determined using a Gaussian distributed random number generator. To determine the primary (or free space) E-field $EWP(i)(r;t)=I′(t)EWP(i)(r)$ for the $ith$ surface current realization, we apply a single point quadrature to obtain

where $Aj$ and $W(i)(r′j)$ are the area and the $ith$ realization of magnetic surface current at the center of the $jth$ triangle, respectively. The value of $W(i)(r′j)$ is determined using a normally distributed random number generator. ${W(i)(r;t);r∈S}$ can be simply considered as $ith$ realization of random magnetic dipoles on the surface $S$⁠.

The presence of the head will distort the primary E-field. As a result, the total E-field induced in the head by current $W(i)(r;t)$ is $EW(i)(r;t)=I′(t){EWP(i)(r)−∇ϕ(i)(r)}$⁠, where $ϕ(i)(r)$ is the scalar potential. Applying the current continuity to the conduction current results in the Poisson’s Equation (4) that can be solved for scalar potential-

where $σ(r)$ is the conductivity at Cartesian location $r$ in the head and the temporal variations are canceled from both sides. The region surrounding the head is assumed to be insulating (i.e., $σ(r)=0$⁠). The mesh is truncated on the outermost boundary of the head (i.e., the scalp boundary) where the normal component of the current density is continuous across the boundary. In other words, $n^⋅∇ϕ(i)(r)=n^⋅EWP(i)(r)$ on the scalp boundary. Here, $n^$ is a normal vector on the scalp surface pointing outward. To solve Equation (4), we approximate the head by a tetrahedral mesh where each tetrahedron is assigned a homogeneous conductivity.

Then, Equation (4) is solved using either our $1st$- or 2^nd-order in-house finite element solvers that have been cross-verified (Gomez et al., 2020b) and available online (Gomez et al., 2020a).

The spatial variation of the total E-field in the brain is approximated as a piece-wise constant within each tetrahedron. The resulting expansion is $EW(i)(r)=$$∑k=1NeLk(r)(e3(k−1)+1(i)x^+e3(k−1)+2(i)y^+e3(k−1)+3(i)z^)$⁠, where $Ne$ is the number of tetrahedrons in the brain, $Lk(r)=1$ in the $kth$ tetrahedron and zero outside it and $e$ is a vector of expansion coefficients and $e3(k−1)+1$⁠, $e3(k−1)+2$⁠, and $e3(k−1)+3$ are the average $Ex$⁠, $Ey$⁠, and $Ez$ in the $kth$ tetrahedron.

To find the spatial variation of $M(i)(r)$⁠, where $i=1,2,...,$$Nm$⁠, we must orthogonalize the set ${EW(1)(r),EW(2)(r),...,$$EW(Nm)(r)}$ with respect to their inner-product

The above establishes an equivalence between the inner-product of two E-fields and the dot-product of two column vectors of a matrix $Z3Ne×​Nm$ with entries $Z3(k−1)+α,i=Vke3(k−1)+α(i)$⁠, where $k=1,2,...,Ne$⁠, $α=1,2,3$⁠, $i=1,2,...,Nm$⁠, and $Vk$ is the volume of the $kth$ tetrahedron. To find orthonormal mode functions from $Z$ we first compute its singular value decomposition (SVD) $Z=USVT$⁠. The modes are defined as

where $U3(k−1)+α,i$ is the $ith$ column and $(3(k−1)+α)th$ row entry of $U$⁠.

In this section, we provide details of the numerical implementation used to compute the coefficients $a(i)$⁠.

The expansion coefficients can be computed from $ETMS(r)$ in the brain and $M(i)(r)$ by evaluating the following inner product integral

However, this requires apriori knowledge of TMS-induced E-field in the head. An alternate form of Equation. (7) can be found by invoking the reciprocity principle. First, we define a new scenario with a mode (⁠$M(i)(r)$⁠) as an electric current source (Fig. 2B). As shown in section 6.3 of the Supplementary File and Gomez et al. (2021), the reciprocity principle dictates that

where, $JTMS(r;t)=I(t)JTMS(r)$ is the electric current density in the coil outside the head and $EM(i)(r;t)=I′(t)EM(i)(r)$ is the E-field generated by the impressed mode current sources $I(t)M(i)(r)$⁠. Equation (8) implies that the integral of the TMS-induced E-field times the mode function is equal to the integral of the E-fields induced by the modes as currents times the TMS coil electric currents. For a detailed mathematical proof and discussion about the reciprocity principle in quasi-static scenarios, the readers are referred to section 6.3 in the Supplementary File and Plonsey (1972); Hasan et al. (2023); Gomez et al. (2021) and the Supplementary File of Gomez et al. (2021).

Evaluating $a(i)$ by using the reciprocity principle in Equation (8) still requires the calculation of $EM(i)(r)$ due to all mode current sources that are inside the head. To reduce the computation, we utilize Huygens’s principle (Balanis, 2012) to replace the mode current distributions $M(i)(r;t)$ by equivalent electric current distributions $JS(i)(r;t)=I(t)JS(i)(r)$ and magnetic current distributions $KS(i)(r;t)=I′(t)KS(i)(r)$ on the Huygens’s surface (⁠$S$⁠). Both of these current distributions generate E-field $EJS,KS(i)(r;t)$ everywhere. Outside of the head, $EM(i)(r;t)$ is equal to $EJS,KS(i)(r;t)$ (Fig. 2B, 2C). Therefore, Equation (8) becomes

As a next step, we invoke the reciprocity principle again (Fig. 2C, 2D), and this results in the following

The last equality involves the TMS coil primary fields because the equivalent scenario (Fig. 2C) has the equivalent surface (Huygens’s surface) current distributions generating fields in empty space (details in section 6.5 of the Supplementary File). For detailed proof of Equation (10), the readers are referred to section 6.6 of the Supplementary File. The Huygens’s surface (⁠$S$⁠) is chosen to be the same as the fictitious surface in section 2.2. A single-point Gaussian quadrature is used to determine $a(i)$⁠. This results in the following

where $Aj$ and $rj$ are the area and the center position of the $jth$ triangle. To evaluate $a(i)$ numerically, we need to first determine the equivalent electric and magnetic currents on the Huygens’s surface for each of the impressed currents $M(i)(r)$⁠, which are described in the next section.

According to the surface equivalence principle (section 7.8 in Balanis (2012)), the fields generated by $M(i)(r;t)$ outside of the surface $S$ can be represented as being generated by $JS(i)(r;t)$ and $KS(i)(r;t)$⁠. These current densities will generate electric fields

and magnetic fields

Note that, the fields inside the surface $S$ in Equations (12) and (13) are chosen to be zero inside the head to enable the replacement of the conductive head with free-space Furthermore, the fields have to satisfy the electromagnetic boundary conditions. Plugging in the electromagnetic boundary conditions on $S$ for the equivalent scenario (Fig. 2C) results in (Equation (7-43a) and (7-43b) in Balanis (2012)) 

The spatial variation of the currents $JS(i)(r)$ and $KS(i)(r)$ can be determined from $HM(i)(r)$ and $EM(i)(r)$⁠, respectively. Balanis (2012) and Peterson et al. (1998) provide more detailed explanations of the surface equivalence principle and Huygens’s principle for interested readers (section 7.8 in Balanis (2012)). A detailed proof of Equation (14) that follows Tai (1998)’s approach is provided in sections 6.4-6.5 of the Supplementary File for interested readers.

Each mode current generates an E-field inside and outside the head. The E-field inside the head is due to a scalar potential generated by the mode current. The E-field outside the head has a negligible scalar potential component and is equal to the time-derivative of the magnetic vector potential generated by the conduction currents and mode currents. Therefore, the spatial variation of the E-field on the surface $S$ is (Gomez et al., 2020a, 2021; Hasan et al., 2023; D. Wang et al., 2023)

where $Head$ is the whole head region. The H-field is just the curl of the magnetic vector potential (i.e. Biot-Savart applied to the conduction and mode currents) and its spatial variation is

To evaluate Equations (15 and 16), we first need to determine $EM(i)(r)$ in the head (procedure shown in Fig. 3). This is done by applying an FEM procedure to $∇⋅σ(r)EM(i)(r)=∇⋅M(i)(r)$ (available online at Gomez et al. (2020a)). The integrals in Equations (15 and 16) are approximated by applying a single point Gaussian quadrature rule for each tetrahedron in the head mesh and are rapidly evaluated using the Fast Multipole Method library (Greengard & Gimbutas, 2012).

Evaluation of $a(i)$ and the TMS-induced E-field requires computation of the primary fields generated by the TMS coil on Huygens’s surface. As such, these must be computed in real-time each time the coil moves. Here, we adopt the approach proposed in Daneshzand et al. (2021); Thielscher et al. (2015) that leverages the fact that the primary fields are functions of position relative to the coil (i.e., translational invariance). As a result, the primary fields are rapidly computed by interpolating samples on a 3D Cartesian grid using a process described next.

Here, we assume a reference coil placement flat on and centered about the x-y plane. The primary fields are sampled on a 3D grid with 4 mm grid spacing with 822,000 grid points (Fig. 4). This grid spacing empirically was found to result in an interpolation error of the order of $10−3%$⁠. Results of interpolation errors for various grid spacing are given in the Supplementary File.

The grid is large enough to contain the whole Huygens’s surface for any relative placement with respect to the coil. Coil placements are defined as coordinate transformation (i.e., a rotation and translation) to the coil relative to the reference coil placement as

where $Rc$ is a $3×3$ rotation matrix and $T0$ is a $3×1$ translation vector. Applying the transformation $T$ to the coil is equivalent to applying $T−1$ to the Huygens’s surface. We apply the transformation $T−1$ to the subjects’ Huygens’s surface instead of the transformation $T$ to the coil. This avoids the need to translate the grid that has $822,000$ grid points and only requires translating $120,000$ Huygens’s surface nodes. Furthermore, this enables us to use the computationally efficient standard multi-linear interpolation method (Kreyszig, 1972). In the standard multi-linear interpolation scheme, the E-field is interpolated at the center of Huygens’s surface triangular facets using the nearby grid point values (where the E-fields have already been pre-computed), and the interpolation is performed along $x^,y^,$ and $z^$⁠, separately, using the ‘interpn()’ functionality in MATLAB (MATLAB, 2022).

This section summarizes the offine mode and surface equivalent current calculation stage and the real-time E-field calculation stage. Algorithm 1 summarizes the critical steps for computing the mode functions. Algorithm 2 describes the four fundamental steps to calculate the TMS-induced E-field in the $ROI$ in real-time while the modes and primary fields are already pre-computed. The algorithms were implemented in MATLAB 2022a (MATLAB, 2022) with built-in GPU functionalities from the ‘Parallel Computing Toolbox’. The current implementation of the real-time stage is on NVIDIA GPUs or any GPU that MATLAB supports. However, the real-time stage requires only two dense matrix-vector multiplications and a multilinear interpolation. As such, it is easily portable to any GPU package with those capabilities.

Note that the induced E-fields $(EW(i)(r))$ in the brain due to random magnetic sources $(W(i)(r))$ span the TMS-induced E-field space. To find an orthonormal basis set of $EW(i)(r)$⁠, we must perform an SVD on $Z$⁠. Instead of performing SVD on the matrix $Z$ itself, we implement the SVD by doing an economic QR decomposition of $Z$ followed by an SVD on $R$ because it is more efficient than doing an SVD on the original matrix. Both methods provide the same result, but the latter is more computationally efficient as we perform SVD on a much smaller matrix $R$ than $Z$⁠.

The algorithm is tested on 14 MRI-derived heads and three distinct TMS coil models. The MRI-derived head models are generated from eight distinct subject MRIs. One is the ‘Ernie’ subject included in SimNIBS-3.2 (Thielscher et al., 2015) and the seven others were collected from Lee et al. (2016). The 16 head models are generated using the ‘mri2mesh’ and ‘headreco’ tools in SimNIBS (Nielsen et al., 2018; Windhoff et al., 2013). The ‘mri2mesh’ models comprise 668,000–742,000 nodes and 3.73–4.16 million tetrahedrons. On the other hand, the ‘headreco’ models consist of 528,000–886,000 nodes and 2.87–4.92 million tetrahedrons. The average edge lengths on the cortex, scalp, and Huygens’s surface are 1.4 mm, 1.65 mm, and 1.66 mm, respectively.

We consider only the five homogeneous concentric compartments such as (from inner to outer) white matter, grey matter, cerebrospinal fluid (CSF), skull, and scalp with corresponding conductivity of 0.126, 0.275, 1.654, 0.01, and 0.465 S/m, respectively (Wagner et al., 2004) in all head models. The time required for the generation of each head model was between 20–24 hours using the ‘mri2mesh’ tool and 1.5–2 hours using the ‘headreco’ tool.

During the preprocessing stage (offline stage or, mode calculation stage), we compute the E-fields in the brain (grey matter and white matter) consisting of 1.31–1.84 million tetrahedrons for ‘headreco’ models and 1.55–1.65 million tetrahedrons for ‘mri2mesh’ models. During the real-time stage, we compute the E-field at the barycenter of each triangular facet on the middle grey matter surface (a surface approximately midway into the grey matter (Stenroos & Koponen, 2019)) consisting of 122,000–289,000 triangular elements for ‘headreco’ models and 241,000–284,000 triangular elements for ‘mri2mesh’ models.

Our method is suitable for modeling any TMS coil; to illustrate, we included different coil types and sizes in this study. The Figure-8 coil consists of two $9$ turn concentric circular loops with inner and outer loop diameters of 53 mm and 88 mm, respectively, which matches the 70-mm Figure-8 #31 in Deng et al. (2013) and is approximated with the coil model described in Gomez et al. (2020b). The circular and double cone coil models were obtained from Makarov (2020) and are models of the MagVenture Cool-40 Rat coil and D-B80 coil, respectively. The coils are modeled with electric dipoles. The coil is centered 5 mm directly above and oriented tangent to a scalp landmark. The number of dipoles in Figure-8, MagVenture Cool-40, and D-B80 coils are 193,536, 57,024, and 22,400, respectively. The Huygens’s surface is just 1 mm outward from the scalp. During the grid interpolation, the primary fields are pre-computed on the grid and the grid is right below the coil but large enough to hold the head and the Huygens’s surface.

We performed a benchmark comparison between our real-time algorithm with the conventional $1st$-order FEM solver. We consider global vector error (GVE) and global magnitude error (GME) as means of comparison between the FEM-calculated E-field $(ETMS(r))$ and the real-time computed E-field $(ETMS(Nm)(r))$ in the $ROI$⁠, defined as follows:

where $‖·‖$ denotes the $L2$ norm and $|·|$ is the magnitude for the E-field. The real-time solver aims to recreate the FEM predicted field, as a result, the reference E-field (⁠$ETMS(r)$⁠) is computed by the in-house built $1st$-order FEM solver (Gomez et al., 2020a). Unlike the inverse relative transformation and multilinear interpolation in the real-time solver, the primary E-field of the reference solver was computed by direct transformation of the coil with respect to the head. Additionally, for a visual comparison of the E-fields, we consider the point-wise relative errors: local vector error (LVE) and local magnitude error (LME), normalized by the largest E-field magnitude in the $ROI$ from the FEM solver defined as follows:

For all of our analyses, we perform simulations for modes 100 to 500, with a step size of 50.

In this section, we compare observed errors for all 16 head models and the 70-mm Figure-8 coil model. We randomly select 1000 coil placements for each head model and calculate the errors by comparing them with the reference $1st$-order FEM solution. The convergence of GME and GVE is shown in Figure 5 as a function of the number of modes. For ‘mri2mesh’ models, the mean GME and mean GVE are below 2% at modes (equals matrix rank) of $325$ and $450$⁠, respectively. For ‘headreco’ models, the required modes are $350$ and $475$ for mean GME and GVE, respectively. There are some outlier errors, primarily corresponding to specific coil placements, which nonetheless remain below 3% for GVE and under 2% for GME at rank $500$⁠. Here, we used the default ‘headreco’ and ‘mri2mesh’ mesh models whose FEM solution is known to have a GVE near 5% (Nielsen et al., 2018). With 400 modes we observed a maximum GVE and GME error of 4% and 3%, respectively, across all simulations. To ensure that the real-time results are just as accurate as the $1st$-order FEM, we estimated the error of the $1st$-order FEM and real-time solutions by using a $2nd$-order FEM as reference (results are given in the Supplementary File). For 400 modes, the resulting difference in GVE and GME between real-time and $1st$-order FEM was, on average, 0.17% and 0.14%, respectively (Fig. S4). Furthermore, the GVE and GME across the 16,000 simulations showed that the real-time solution with 400 modes is up to 1.3% (1.7%) and 0.7% (1.1%) more (less) in agreement, respectively, than the $1st$-order FEM to the $2nd$-order FEM solution (Fig. S5 and S6), indicating that the real-time results are as accurate as the $1st$-order FEM ones.

Here, consider the relative accuracy performance of the method for three distinct coil models 70-mm Figure-8, MagVenture D-B80 coil, and Cool-40 Rat coil. We consider the sample mean error over 1000 random coil placements over the scalp on each of the 16 head models and use the $1st$-order FEM solution as a reference. The mean GME and the mean GVE are shown in Figure 6. The mean GME is below 2% for ‘mri2mesh’ (as well as ‘headreco’) head models at ranks above $325$⁠, $425$⁠, and $375$ (⁠$350$⁠, $375$⁠, and $375$⁠) for the Figure-8, D-B80, and Cool-40 coils, respectively. Due to the unique bending shape, the D-B80 coil has an E-field that has more fine features relative to the others, thereby, requiring comparatively more modes for its expansion. The mean GVE is below 2% for ‘mri2mesh’ (‘headreco’) head models at ranks above $450$⁠, $550$⁠, and $475$ (⁠$475$⁠, $525$⁠, and $500$⁠) for the Figure-8, D-B80, and Cool-40 coils, respectively. All coils exhibit similar errors. However, compared to the D-B80 and Cool-40 coils, the Figure-8 coil model converges to the 2% error limit with fewer modes.

Figure 7 shows the convergence of error metrics as a function of modes for the $ROI$ regions where $ETMS(r)≥0.7maxr∈ROI(ETMS(r))$⁠.)). In this scenario, it requires 350 (500) modes for mean GME (GVE) to converge to the 2% error limit (considering all coil models and all head model types). This is in contrast to 425 (550) (Fig. 6) for E-fields across the whole cortex.

In this section, several exemplary simulation results of coil placements over the scalp of the ‘Ernie’ head model are shown. Figure 8 shows the specific coil placement over the scalp, the associated E-field induced on the middle grey matter surface computed in Real-time and FEM solvers, and the corresponding LME and LVE distributions. The E-field distributions predicted in real-time and FEM are visually indistinguishable. Furthermore, the peak E-field is the same up to 0.65 V/m in all cases shown. The maximum LME for each scenario (top to bottom) are 3.7%, 3.6%, 2.7%, 3.1%, 2.9%, and 3.2%, whereas the corresponding maximum LVE is 4%, 3.8%, 3.9%, 3.2%, 3.4%, and 4.5%, respectively. These results indicate that the real-time predicted E-field distributions are equally valid to the FEM $1st$-order solutions, which have an LME of about 5%.

Here we consider single-precision arithmetic during real-time computation, as the results remained unchanged up to seven digits relative to double-precision. Figure 9A and 9B show the mean mode calculation time (pre-processing time) as a function of the number of modes across eight ‘mri2mesh’ models and eight ‘headreco’ models, respectively. The pre-processing stage is computed using an AMD Rome 2.0 GHz CPU. On average the pre-processing time required to generate $400$ modes is $38$ hours and $34$ hours for ‘mri2mesh’ and ‘headreco’ models, respectively. The total computational runtime is 2 FEM runs per mode. Our single-threaded implementation requires, on average, 3 minutes of computation time per simulation in an AMD Rome 2.0 GHz processor. The required pre-processing computations could be significantly sped up by using multiple threaded FEM solvers.

The core FEM and reciprocity integrals are implemented in C complied for MATLAB and the real-time stage was implemented in MATLAB with GPU functionalities. Figure 10 shows the computational reconstruction time as a function of the number of modes for a GPU and CPU, respectively. The above result was obtained by running 48,000 simulations across the 16 head models and 3 coil models. The mean reconstruction time to predict the E-field over the intermediate grey matter-white matter surface for a fixed coil placement is 2.2 ms for 400 modes in an NVIDIA RTX 3080 GPU with a maximum time of 3.8 ms. The first step of the real-time TMS (coordinate transformation of the Huygens’s surface) takes only 0.03 ms, whereas the CPU takes 0.9 ms (a $30$ times speed-up). The second step (multi-linear interpolation of primary fields) takes 1.7 ms in a GPU and 37.40 ms in a CPU. In other words, the GPU requires $22$ times less run-time than the CPU. The mode coefficient calculation is completed in 0.4 ms in a GPU and 1100 ms in a CPU (i.e., a $2750$ times faster). Finally, the fourth step (TMS-induced E-field computation) is rapidly computed in 0.03 ms in a GPU versus 105 ms in the CPU (⁠$3500$ times faster). Overall, the mean total time for estimating the TMS-induced E-field using $400$ modes is 2.2 ms in a GPU and 1200 ms in the CPU (⁠$550$ times faster). This ratio could be improved by accelerating Huygens’s surface coordinate transformation.

Figure 11 shows the required memory in the reconstruction stage for both CPU (AMD Rome CPU, 2.0 GHz) and GPU (NVIDIA RTX 3080-10 GB) across 14 head models. The total required memory in the reconstruction stage is the same for both the GPU and the CPU. The differences in GPU and CPU memory requirements stem from the fact that the Matlab environment requires overhead that is not accounted for in the GPU memory. In other words, the GPU only has all required data structures (e.g., modes, surface currents, and interpolatory primary fields). When the real-time computation is performed in the GPU, the required mean CPU and the GPU memory for 400 modes are 1.3 GigaBytes (GB) and 3 GB, respectively. Additionally, the required mean CPU memory during real-time computation in the same CPU is 4.3 GB. Additionally, section 6.7 in the Supplementary File provides an estimation of the floating point operations (FLOPS) required for the real-time stage.

The real-time TMS-induced E-field approximation method excels in rapidly calculating the E-field over the middle grey matter surface using GPU acceleration, achieving this computation within 4 ms for any coil placement over the scalp of the subject. Our detailed benchmarking analysis indicates that $450$ modes are enough to achieve a mean GME and a mean GVE below 2%. Additionally, with only $400$ modes, the maximum GME and maximum GVE are 2.8% and 4.2%, respectively. Therefore, with $400$ modes, the agreement between the real-time predicted E-fields and $1st$-order FEM is closer than the expected 5% error of the FEM solution. These results are consistent across multiple subjects, head model construction pipelines, coil types, and coil placements.

In this study, we conducted a benchmark comparison between our real-time algorithm and the $1st$-order FEM solver, which has a known relative error of about $5%$ (Thielscher et al., 2015). Correspondingly, we observed that the convergence of the mode expansion greatly decelerated below $5%$⁠. This is likely because, below the FEM solution error threshold, we are increasingly spanning the erroneous part of the solution, which is not expected to be spanned by a small number of modes. We additionally used $2nd$-order FEM, which results in a more accurate E-field prediction, for a single subject to generate the mode expansion. The results given in the Supplementary File indicate that using a more accurate solver (i.e., FEM $2nd$-order) could result in slightly improved convergence. At the current time, running the $2nd$-order FEM requires excessive preprocessing time, thereby limiting the applicability of this method. Although not pursued here, this tool could be implemented using efficient implementations of BEM (Makarov et al., 2020) that are more accurate than $1st$-order FEM (Gomez et al., 2020b) to achieve improved convergence. Furthermore, $1st$-order FEM is currently the most commonly used tool for TMS E-field dosimetry.

Our method relies on the existence of a small set of modes that span the possible range of E-fields induced by a TMS coil. From spherical solutions, we expect that the total required modes are larger for superficial regions relative to deeper ones. As such, this method could in principle be applied to determine E-fields on more superficial compartments like the CSF, skull, and skin. However, a larger number of modes would likely be necessary in such cases, which could potentially limit its applicability.

The preprocessing stage can be done using any available FEM or BEM code. Using our $1st$-order in-house FEM solver (Gomez et al., 2020a, 2020b), which is optimized for accuracy, requires less than $40$ hours for $400$ modes at 180 s per FEM run. However, a faster implementation of the underlying FEM can speed up the process. For example, the SimNIBS implementation of FEM requires a runtime of 29.8 s, which reduces the total time to 6.6 hours.

The real-time stage requires only under 4 ms (assuming 400 modes) in an NVIDIA RTX 3080 GPU and 1.2 s in an AMD Rome 2.0 GHz CPU to predict the TMS-induced E-field. The reconstruction time is almost independent of the number of modes in a GPU. Therefore, for better accuracy, a higher number of modes can be accommodated in a higher-end GPU.

We found empirically that the pre-processing stage should be computed in double-precision arithmetic whereas the real-time stage can be performed in single-precision arithmetic for accommodating more modes inside the GPU and a larger $ROI$ (i.e., the whole brain/head). A GPU with only 4 GB of dedicated memory is large enough for this process accommodating up to 400 modes, pre-computed primary fields, Huygens’s surface nodes, and an $ROI$ of size 220,000 samples on the middle grey matter surface.

Future work includes connecting this GPU-based E-field computation and GPU-based visualization within one screen refresh in real-time. The integration of real-time E-field modeling in combination with neuronavigation and TMS stimulator control will enable several avenues of research and clinical intervention. For example, this methodology will enable support for more reliable dosing throughout a TMS session by automatically changing TMS intensities accounting for coil placement variability. It will also allow the development of methods for faster motor threshold determination applying fewer TMS pulses by using E-field-based instead of grid search approaches. Furthermore, this approach will allow the development of rapid coil placement optimization for updated brain targets and E-field constraints during the TMS procedure. Updated brain targeting may be derived from behavioral (e.g., task performance, or treatment indicators) or physiological (e.g., heart rate variability, or electroencephalography) response even in closed-loop settings.

We have developed a method for rapidly calculating the TMS-induced E-field. This method is based on a functional generalization of the probabilistic matrix decomposition and generalized Huygens’s and reciprocity principles. The initial preprocessing stage takes approximately 40 hours to complete given the algorithm requires solving 400 FEM simulations and 400 ADM simulations. However, subsequent E-field computations can be done within 4 ms in a GPU and 1.2 s in a CPU over the middle grey matter surface (containing, on average, 216,000 barycentric points). The resultant E-field has accuracy comparable to standard FEM solutions. Notably, this computational performance can be achieved using a standard GPU with a dedicated memory of only 4 GB, making it practical for many users. This framework enables real-time E-field calculation in the cortex for arbitrary coil design and coil placement, and generalizes well for distinct head model pipelines, underscoring its adaptability and suitability for a wide range of applications.

The pre-print was submitted to BioRxiv. However, the article is not submitted elsewhere. No AI tools have been utilized in the preparation of the article. There is no experiment conducted for this research involving live animal or human subjects.