Abstract

The aim of this paper is to propose an approach for an accurate and fast (real-time) computation of the electric field induced inside the whole brain volume during a transcranial magnetic stimulation (TMS) procedure. The numerical solution implements the admittance method for a discretized realistic brain model derived from Magnetic Resonance Imaging (MRI). Results are in a good agreement with those obtained using commercial codes and require much less computational time. An integration of the developed code with neuronavigation tools will permit real-time evaluation of the stimulated brain regions during the TMS delivery, thus improving the efficacy of clinical applications.

1. Introduction

Transcranial magnetic stimulation (TMS) is a noninvasive and painless technique that delivers brain stimulations via externally applied magnetic fields generated by a coil positioned above the patient’s scalp surface [1]. The coil, fed by a current pulse, generates a time varying magnetic field that penetrates into the head’s tissues placed in the near field zone of the coil. According to Faraday’s induction law, rapid variations of magnetic field induce an electric () field in the brain interacting with the neuronal cells of the cortex. From a macroscopic point of view, brain activation occurs in the regions where the induced reaches a certain threshold.

In the last years, TMS has become an important tool not only in research and diagnostic fields, but also in clinical applications, being a promising alternative treatment for a broad range of neurological and psychiatric disorders including strokes, Parkinson’s disease, tinnitus, epilepsy, and depression [2].

Despite the efficacy of many TMS clinical applications, as for all stimulation techniques of the nervous system, more research is needed on dosimetry [3, 4]. In fact, knowledge of the induced field distribution is fundamental for evaluating position, level, and extent of the stimulation. In addition dosimetry, when coupled with biophysical neuronal models [57], is useful for interpreting experimental results, designing more efficient applicators [810], and clarifying possible mechanisms of action of electromagnetic fields on the brain [11].

The need for an accurate calculation has led to the development of increasingly realistic brain models, starting from simple homogenous spheres [1214] up to recent patient specific models obtained from Magnetic Resonance images. The patient specific models include a realistic number of tissues [1517] and even an anisotropic behavior of white matter [1820].

Various computational methods have been used to solve the problem, based on either differential or integral equations. They are the finite difference time domain (FDTD) method [21, 22], the finite element method (FEM) [19, 2327], the impedance method (IM) [15, 18, 28, 29], the boundary element method (BEM) [16, 30], and the surface integral equation (SIE) [31].

Due to the high wavelength of the stimulating field with respect to the head size, the majority of these methods applied to TMS use quasistatic approximation to determine the distribution of the field induced in the human brain [12, 1517, 27, 3240]. Conversely, some authors [31] recently proposed a more realistic physical description of the problem including the electromagnetic field propagation.

In recent years, the calculation of the field distribution inside the brain has assumed a further fundamental role. In fact, the increasing use of computer-aided neuronavigation tools requires a correct coil positioning with respect to specific brain region identified by MRI patient specific images [30, 4143], thus enhancing the necessity of exact field knowledge.

Such an application requires a very fast calculation, in order to furnish real-time results during the coil displacements over the scalp.

To fulfill the real-time requirement, most of the existing online TMS navigation tools estimate the induced field using spherical conductor models [41, 42]. Recently, Nummenmaa et al. [30] presented an efficient BEM approach for realistic brain model. If 1-layer BEM is considered, the computational time was comparable to that consumed when a locally fitting spherical model of the brain is considered [30]. However, the field was calculated only on the brain surface, neglecting its penetration in deeper regions. In fact, methods such as BEM and SIE deal with fully populated matrices resulting in high storage requirements and long computational time for fine spatial resolution [40].

In this paper we propose a tool for the calculation of the field within the whole brain based on an improved revision of the admittance method [22, 44, 45], aimed at solving the real-time requirement introduced by the neuronavigation computer-assisted technologies. The admittance method is fast and appropriate for this kind of near field, quasistatic approximation and does not require meshing codes, needed when using FEM methods, since the regular 3D grid for the brain model is directly extractable from MRI images [40].

Real-time calculation is achieved by decoupling the coil domain, where the magnetic potential is calculated only once to build up a library file for each considered coil, from the brain domain, where the field is calculated for each new coil position. This procedure allows us to recalculate in a few seconds the field due to the real-time positioning of the coil, using a rototranslation operation that dynamically maps one domain into the other. The field calculation, compliant with the strict time requirement, is obtained with acceptable spatial resolution and good accuracy, evaluated by comparison with commercial codes.

Such a strategy allows (i) considering different kinds of TMS coils, that is, an easily expandable “Coil Library”; (ii) obtaining a fast and reliable solution, in agreement with real-time requirements of clinical applications; (iii) calculating the field distribution within a whole realistic brain volume using a spatial resolution of a few mm. Therefore, the developed tool can be efficiently integrated with available neuronavigation systems.

2. Models and Methods

2.1. The Solution Approach: Decoupling and Calculation Domains

The coil is modeled as a discrete sequence of segments carrying a sinusoidal current (Section 2.2) in a discretized analysis volume (coil domain). Under the quasistatic hypothesis we obtain the magnetic vector potential by solving the discrete Poisson vector equation in the coil domain (Section 2.3).

As for the brain, we consider both homogeneous and inhomogeneous models. The homogeneous model is a realistically shaped volume completely filled with grey matter (GM) and immersed in the cerebrospinal fluid (CSF); the inhomogeneous one accounts for both GM and white matter (WM) (Section 2.4).

The brain domain is discretized using a 3 mm spatial step. This choice allows us to significantly reduce calculation time and memory occupation while maintaining an acceptable spatial resolution.

The electric field is calculated in the brain domain by the superimposition of the primary and the secondary fields, using the expression where and are the magnetic vector potential and the electric scalar potential, respectively. The electric scalar potential is obtained using the admittance method (Section 2.5) in the brain domain.

The overall strategy used to fulfill versatility and fast calculation, while maintaining an acceptable spatial resolution, is summarized in the flowchart of Figure 1. For each coil considered, the distribution is calculated in an analysis volume surrounding the coil, discretized using a Cartesian cubic grid 3 mm on a side. Since the analysis domain must be large enough (at least 1.75 times the brain size) to include the whole brain volume for any realistic position of the coil, this calculation is generally time-consuming, but it is carried out offline only once. The fields are then stored in “ad hoc” library files. This library, together with the file containing the brain model, becomes the input of the revisited admittance method whose engine is devoted to the real-time field calculation (Figure 1). To change the coil position, the brain domain is dynamically mapped (rototranslated) into the coil domain (Section 2.5) on the basis of the coordinates of three fiducial points of the coil given in the brain Cartesian reference system (Figure 3).

To further reduce the computational time, V and thus the field are evaluated in an analysis domain, whose external boundaries are placed at a “cut-off” distance from the brain, outside which the influence of the brain-CSF interface on becomes negligible. Here, this cut-off distance has been fixed to 9 mm, that is, 3 mesh cells (Figure 3).

Since the admittance method works in the frequency domain, to avoid time-consuming spectral decomposition, the biphasic TMS stimulus is replaced with the best approximating monochromatic signal (Section 2.2) [18].

The possibility of carrying out a real-time calculation relies on the adopted strategy that provides:(i)reading of the distribution in the coil domain from a library file;(ii) and calculation in the domain within the cut-off distance from the brain surface, suitably rototranslated with respect to the coil domain.In addition, a 3 mm spatial resolution and the monochromatic assumption have been considered to speed up calculation while maintaining an acceptable accuracy, as shown in Results. The procedure described in Figure 1 is implemented in C++ environment.

2.2. Source Model

In this study we consider the figure-of-eight commercial coil MAG-9925-00 (Magstim). However, thanks to the modular code implementation (Figure 1), as many coils as one wants can be added to the library without changing the fast engine calculation.

Under the assumption of dimensionless transversal section of the wire, the coil is discretized in a sequence of 2721 segments , describing the current path (Figure 2), using a MATLAB tool suitably developed. Each segment is identified by the coordinates of its central point in the coil reference system (Figure 3).

The current flowing in the MAG-9925-00 coil energized with the biphasic stimulator Magstim Rapid2 was measured using a Tektronix DPO 2024 (200 MHz; 1 GS/s) digital oscilloscope equipped with amperometric probes Tektronix TCT 0030 A (120 MHz; 30 A).

In agreement with theoretical predictions, the measured current time course is well fitted by a damped sinusoid with period μs, amplitude 100 A (stimulator set at the 1% of the maximum deliverable power), and attenuation constant  s−1.

Since the field induced in the brain reaches its maximum value at the time instants where the time derivative of the current is maximum (i.e., at the beginning of the stimulus), generally the current signal is assimilated to a pure sinusoid that has the same slope of the damped sinusoid at the origin [18]. This slope, calculated in the first 20 μs, is perfectly approximated (error < 1%) by a pure sine wave with frequency  kHz and amplitude  A.

Therefore, the code for the field calculation works at the approximated frequency of 3 kHz; the current amplitude is an input variable that can be changed depending on the power level set on the stimulator.

2.3. and Fields Calculation

Under the quasistatic approximation, assuming a uniform current density over the coil cross section [31] and considering the meshed coil, the magnetic vector potential at any space point can be calculated as the superimposition of all contributions from linear segments used to discretize the coil model:where is the free space permeability, I and are the intensity and direction of the current flowing through the coil, and is the distance from the observation point to the source point on the coil.

The uniform current density is a plausible assumption since the wire thickness does not exceed two times the skin depth at the considered frequencies.

The distributions stored in the library files are calculated for  A. Due to the linear relationship between and , the actual distributions are simply multiplied by a factor equal to the current intensity I.

values, stored in the library files, are used to calculate the distribution.

Using the same assumption, magnetic flux density has been also calculated for comparison with results from commercial codes, in order to validate the coil modeling.

at any space point is calculated using the discretized version of the Biot-Savart law:where is the unit vector of .

2.4. Human Brain Model

The developed code is able to operate with any patient specific voxeled brain model coming from MRI.

In this paper the brain model is obtained from the Virtual Population member Duke (v.1.0, Zurich Med Tech AG, Zurich, Switzerland, [46]), that is, a 1 mm resolution voxeled male model in which a number from 1 to 77 encodes for 77 different tissues. We consider only a box of  mm3 containing all the brain structures; this box is processed in order to have a simplified model consisting of only three materials, that is, GM (that takes into account the cortex and all the nuclei like thalamus, hypothalamus, etc.), WM (that includes the structures that can be assimilated to nerve fibers), and CSF (that incorporates all the other tissues).

In this paper we consider both the three-tissue model, referred to as inhomogeneous, and a further simplified model where all tissues in the brain are treated as GM (homogeneous).

To speed up the solution process while keeping an acceptable spatial resolution in a real-time clinic framework, the brain model is undersampled and discretized in cubic cells 3 mm on a side.

Due to the fact that biological tissues do not possess magnetic properties, the permeability values are taken to be , that is, the free space permeability.

All materials are considered as lossy ones with complex conductivity equal to where σ and are the electric conductivity and the relative permittivity, respectively, is the permittivity of the free space, and is the operating frequency of 3 kHz.

The used values are = 6.7 × 104,  S/m for the GM and ,  S/m for the CSF [47]. WM is considered isotropic with = 3.01 × 104, = 0.065 S/m [47].

2.5. Calculation of the Induced Field
2.5.1. The Admittance Method

Among different numerical approaches for the field calculation inside the brain we chose to adopt the admittance method since it allowed us to explore the whole brain volume with reduced computational time [22, 44, 45].

The admittance method is based on a finite difference (FD) approach to the solution of Maxwell’s equations in the frequency domain in quasistatic condition [48], which is considered valid at the typical frequency of a biphasic TMS pulse (≤3 kHz).

The analysis domain is divided into homogeneous cells () centered at the point of coordinates () in a Cartesian grid that can be represented as a network of lumped electrical elements (see Figure of [45] for a simplified 2D representation) [45, 49, 50].

Applying the Kirchhoff law at each node () of the network leads to a linear system whose equations are given bywhere is the unknown scalar electric potential at each network node, and () are the admittance and the vector potential components at each surface of the discretizing cell, and ω is the operating angular frequency.

The admittance values are calculated from the complex conductivity of (4) as The linear system of (5) is solved using the overrelaxation iterative technique [22, 44, 45]. This choice is based on the best compromise among easy implementation, computer memory occupation, and speed of convergence [51, 52]: with convergence factor α between 1 and 2.

The iterative procedure ends when the error at the th step (left side of (8)) falls down below a defined tolerance level e:In this work is fixed to 1 × 10−7.

2.5.2. Rototranslation

Rototranslation is a fundamental operation for the field calculation: it permits a dynamical mapping of the brain domain into the coil domain.

If () are the coordinates of the origin of the coil Cartesian reference system (local system) with respect to the brain Cartesian reference system (global system), and (), (), and () are the coordinates of the unitary vectors identifying the local system with respect to the global one, the one-to-one correspondence between the coordinates of any point in the local system () and the correspondent point in the global system () is given by the linear relationship: where the coefficients , , and are calculated from the coordinates of three fiducial points () on the coil (Figure 3 shows a simplified sketch of the operation).

In the two discretized domains the correspondence between cells is not necessary one to one, since the rototranslation operation is not conformal, but the maximum committed error is equal to the cell size (3 mm); thus, it is inside the accepted tolerance.

3. Numerical Results

3.1. Accuracy Evaluation

In this section we want to evaluate how much the assumption of homogeneous brain, the spatial resolution of 3 mm, and the thickness of the shell around the brain (3 mesh cells), adopted in our code to speed up calculation, can affect the accuracy of results and how much it can be considered well acceptable within the approximation of a real-time clinical use. Moreover, to validate the used method, we compare our results with those obtained using a commercial software (SIM4LIFE, ZMT Zurich MedTech AG), based on a different solution method, on an inhomogeneous brain model with 1 mm spatial resolution. In this last comparison, to reduce possible truncation errors at the domain boundary, we solved the admittance method on a volume including a CSF shell of 10 mesh cells around the brain.

SIM4LIFE is a simulation platform developed for modeling of interactions between physical stimuli and the human body. In particular, this software can deal with medical image data obtained from MRI. The simulation platform includes many physics solvers; we chose to work with the Magneto Quasistatic module included in the Quasistatic EM Solvers (P-EM-QS) that solves problems in static or quasistatic EM regimes applying the FEM method on graded voxel meshes. To compare results on the same mesh grid we decided to work with SIM4LIFE by setting a regular mesh constituted of cubic voxels 1 mm on a side.

The stimulation is the MAG-9925-00 coil placed horizontally with center coordinates in cm (9.1, 10.9, and 18.4) and handle orthogonal to the coronal plane (Figure 3), modeled in SIM4LIFE as a bidimensional sketch, using the same MATLAB file developed for the C++ code.

A first comparison on field distributions generated by the coil within the brain volume reveals that differences between our method and SIM4LIFE solution are always below 0.1%.

Moving to the calculated field, Figure 4 reports the distributions of the induced field on the top surface of the brain (a, b, c, and d) and on a coronal section at  cm (e, f, g, and h).

Figure 4 shows that in all cases the same stimulated regions are predicted under the coil center, with distributions dependent on the shape of the cortical surface. There is no significant difference between homogeneous and inhomogeneous brain with 3 mm of resolution, whereas induced is higher under the coil and goes deeper inside the cortex for 1 mm resolution.

Moving to a quantitative analysis, we will consider observables suitable to compare intensity, spread, and localization of the estimated stimulation. The stimulation intensity is measured by the maximum field () induced in the brain. To have a robust localization of the stimulation position, useful in clinical practice, we define the stimulation center (SC) as a sort of center of mass of the volume extension of the stimulation, where geometric coordinates are weighted by the intensity at each node [30]. Such a volume has been considered as the space where overcomes the 80% of its maximum value (VOL80). Moreover, another comparison between the results obtained with the admittance method and SIM4LIFE was performed computing the Symmetric Mean Absolute Percentage Error (SMAPE) [53, 54] between the values inside the same volume VOL80.

The SMAPE is calculated as follows: where denotes the th cell, is the number of cells inside VOL80, and and are the field magnitudes inside VOL80 calculated with the admittance method and SIM4LIFE, respectively.

Table 1 summarizes , SC coordinates, and SMAPE, obtained for the distributions calculated using SIM4LIFE and the admittance code in the inhomogeneous brain model, 1 mm of spatial resolution.

Table 1 shows that results on obtained with the same brain model (inhomogeneous 1 mm) but using the two different approaches (admittance and SIM4LIFE) differ only by 1.2%, with a distance below 1 mm between the estimated SCs. The similarity of the results obtained with the two methods is also confirmed by the SMAPE that, for the volume chosen, is less than 7%.

The observed variability of results is well below that obtained using the same SIM4LIFE software with different mesh settings; thus, such results confirm the goodness of our solution method.

When considering the solution obtained with the admittance method in the brain domain used to speed up calculation (homogeneous, 3 mm of resolution, and CSF layer of 3 mesh cells), the calculated is underestimated by 26% with respect to the solutions obtained with SIM4LIFE, gridded with the same 3 mm resolution. Still the results, at this spatial resolution, remain consistent and the SMAPE at 80% is calculated as 14.58%.

It seems evident that the implemented methodology can estimate the stimulated regions in good agreement with those calculated by the commercial code SIM4LIFE using a realistic inhomogeneous brain model with resolution of 1 mm. Therefore when more accurate offline computations are desired (e.g., for planning TMS therapy sessions) our methodology used with the detailed brain model is confirmed to be affordable.

Moreover, when real-time results are required for navigation during experiments or clinic therapy, we have to calculate the field in a brain model using 3 mm of resolution and a surrounding CSF layer of only 3 cells. In this case, results are obtained within 20 seconds on a PC 2.8 GHz Intel Core i7, both for the homogeneous and the inhomogeneous models. Nevertheless this coarse resolution does not affect significantly the brain region targeting, introducing an underestimation less than 30% of the maximum level that can still be acceptable in real-time applications [30].

The use of homogeneous or inhomogeneous brain models with 3 mm resolution affects neither computational times nor accuracy of the stimulation estimate. Thus one can choose to use a simple homogeneous model to represent different patients brains (general purpose) or a more realistic inhomogeneous one for patient specific applications [18].

3.2. Application of the Solution Strategy

In this section we want to show an example on how the developed code can operate online during a TMS section, following the flowchart of Figure 1. To fulfill the real-time requirement we use the homogeneous brain model discretized using 3 mm spatial step and we apply the admittance method for the field calculation inside a conductor volume including the brain and a CSF surrounding layer of 3 mesh cells thick (Section 2.1).

At the beginning of the TMS session the operator chooses the coil and the current intensity. Then, after loading the library files corresponding to the brain model and the vector potential generated by the coil, the code asks for the coordinates of three fiducial points identifying the coil positioning. After giving results on distribution inside the brain volume and on the stimulation center, the code asks for a new coil positioning until the end of the section.

In this example the coil is the MAG-9925-00 with a current of 3 kA flowing through it. The operator moves the coil in four positions identified by the following coordinates (in cm) of the three fiducial points (, , and in Figure 3) with respect to the brain:(1)C (), H (), and L (), corresponding to the MAG-9925-00 coil placed horizontally at the top of the brain (Figures 5(a) and 5(e));(2)C (), H (), and L (), corresponding to the coil translation of  cm,  cm, and  cm with respect to position () and a rotation of 30° around the axis through the center of the coil parallel to the -axis (Figures 5(b) and 5(f));(3)C (), H (), and L (), corresponding to the coil translation of  cm,  cm, and  cm with respect to position () and a rotation of 60° around the axis through the center of the coil parallel to the -axis (Figures 5(c) and 5(g));(4)C (), H (), and L (), corresponding to the coil translation of  cm,  cm, and  cm with respect to position () and a rotation of 90° around the axis through the center of the coil parallel to the -axis (Figures 5(d) and 5(h)).During the offline step, loading of files containing the distribution in the coil domain takes 6 s. This time depends only on the size of the coil domain.

For each position, the fast engine calculation takes 14 s to calculate the distribution inside the brain and the coordinates of the stimulation center (Figure 5). This time does not depend on the specific position but only on the number of cells discretizing the calculation volume.

During the whole TMS session, the developed code requires only 62 s to calculate the maps and the stimulation centers for all the coil positions, around 15 seconds per position (Figure 5 and Table 2). As expected, the stimulated region is always under the coil center, with distributions dependent on the geometry of the brain surface.

Table 2 summarizes and SC calculated in each position.

Results of Table 2 confirm that the stimulation centers are located in the brain cortex under the axis of the coil passing through its center. The higher values of calculated in positions 30° and 60° are due to coil closer to the brain surface.

4. Discussion

Here we present a fast and robust computational approach for field calculation inside realistic brain models. The computational approach is based on an improved revision of the admittance method. In addition, the strategy adopted to fulfill the real-time requirement decouples the domains where and fields are calculated. This allows us to calculate the field inside the brain for each new coil position by dynamically mapping the brain model into the coil domain avoiding to be recalculated each time. Both domains are discretized with  mm3 cubic cells resolution to obtain results in less than 20 s.

Results were compared with those obtained using a more detailed brain model (i.e., 1 mm resolution) and a different solution method implemented by the commercial software SIM4LIFE. Comparisons show good agreement between results obtained with the two methods using the detailed brain model (i.e., 1 mm of resolution). The less refined resolution (i.e., 3 mm) determines an acceptable underestimation of the induced maximum field while maintaining a good accuracy in the stimulation positioning.

The presented methodology could benefit clinical TMS application in several ways. It may be used with the inhomogeneous  mm3 resolution brain to improve the individual targeting and dosing of TMS during the offline planning of therapeutic TMS sessions. However, the actual advantage is the possibility of performing a real-time and accurate prediction (i.e., a few seconds) of the field, using a 3 mm resolution, inside realistic brain regions with different coil models and any coil position and orientation.

Some limitations of the proposed approach rely on the head model simplification. The brain model inside the CSF accounts for the tissue interface closest to the stimulated regions but does not represent a realistic scenario. The real brain is surrounded by various tissues such as adipose, skull, and scalp. In practice, the presence of these tissues may affect the overall distribution of the fields in a way that could be better investigated and quantified in a future research. Indeed, the analysis in [55] showed that the inclusion of the skull would not affect the distribution of currents in the adjacent cortex of the brain, while other authors claim the contrary [38]. Concluding, a direct consequence of this work would lead to including the main of the surrounding tissues such as skull and scalp.

5. Conclusions

This paper presents a fast, versatile, and robust approach to calculate the field induced by different TMS coils within a realistic brain model obtained from MRI images.

The presented approach is practical in the sense that it can be readily integrated with neuronavigation tools to accommodate online TMS, leading to possible improvement in the efficacy of TSM clinic applications.

Future work will be related to the development of a realistic model of the human head, taking into account different biological tissues (bone, muscle, and skin), to evaluate the effect of the adopted head simplification.

Conflict of Interests

Gianluigi Rubino is a full-time employee and Paolo Tampieri is the CEO of EMS s.r.l. (Bologna, Italy). This does not alter professional judgment concerning the validity of the research conducted for the publication of this paper. All the other authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors want to thank Rosanna Pinto and Simone Chicarella for assistance and availability in performing measurements on the coil with the amperometric probes and Stefano Pisa for having introduced them to the smart potentialities of the admittance method and for fruitful discussions.