Next Article in Journal
Task Cortical Connectivity Reveals Different Network Reorganizations between Mild Stroke Patients with Cortical and Subcortical Lesions
Previous Article in Journal
Combination Treatment of Locoregionally Aggressive Granulomatosis with Polyangiitis and Cranial Base Infiltration
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Review

A Review of Formulations, Boundary Value Problems and Solutions for Numerical Computation of Transcranial Magnetic Stimulation Fields

by
J. A. Pérez-Benítez
1,*,
P. Martínez-Ortiz
1 and
J. Aguila-Muñoz
2,*
1
Laboratorio de Bio-Electromagnetismo, ESIME-SEPI, Edif. Z-4, Instituto Politécnico Nacional, Mexico City 07738, CDMX, Mexico
2
CONAHCYT—Centro de Nanociencias y Nanotecnología, Universidad Nacional Autónoma de México, km 107 Carretera Tijuana-Ensenada, Apartado Postal 14, Ensenada 22800, BC, Mexico
*
Authors to whom correspondence should be addressed.
Brain Sci. 2023, 13(8), 1142; https://doi.org/10.3390/brainsci13081142
Submission received: 29 May 2023 / Revised: 22 July 2023 / Accepted: 24 July 2023 / Published: 29 July 2023
(This article belongs to the Section Computational Neuroscience and Neuroinformatics)

Abstract

:
Since the inception of the transcranial magnetic stimulation (TMS) technique, it has become imperative to numerically compute the distribution of the electric field induced in the brain. Various models of the coil-brain system have been proposed for this purpose. These models yield a set of formulations and boundary conditions that can be employed to calculate the induced electric field. However, the literature on TMS simulation presents several of these formulations, leading to potential confusion regarding the interpretation and contribution of each source of electric field. The present study undertakes an extensive compilation of widely utilized formulations, boundary value problems and numerical solutions employed in TMS fields simulations, analyzing the advantages and disadvantages associated with each used formulation and numerical method. Additionally, it explores the implementation strategies employed for their numerical computation. Furthermore, this work provides numerical expressions that can be utilized for the numerical computation of TMS fields using the finite difference and finite element methods. Notably, some of these expressions are deduced within the present study. Finally, an overview of some of the most significant results obtained from numerical computation of TMS fields is presented. The aim of this work is to serve as a guide for future research endeavors concerning the numerical simulation of TMS.

1. Introduction

Transcranial magnetic stimulation (TMS) is a method that modulates brain activity using an electric field induced in the brain by potent magnetic field pulses [1,2,3,4]. This method is based on the application of magnetic field pulses to the brain cortex using coils, inducing electric fields oriented along certain paths to stimulate or inhibit brain activity in these areas. Several works [5,6,7,8,9] have shown that after TMS treatments, patients show improvement in their conditions, such as depression [5], mania, obsessive–compulsive disorder, post-traumatic stress disorder and schizophrenia [6,7,8,9]. TMS in combination with EEG methods [10,11] can also be used as a tool too study brain functions. However, the efficacy of TMS for treatments is highly influenced by factors such as the efficiency [12,13], focality [14,15,16,17,18,19] and convenient distribution [15] of the induced electric fields. Furthermore, these factors are determined by several parameters, such as the quantity, [16,20], arrangement [16,17,20] and geometric parameters of excitation coils [16,18,21,22]; coil distance and orientation with respect to the brain [20,23,24]; excitation current waveform [25,26,27,28,29]; and electric properties of the brain tissues [30,31,32,33,34], among others.
Based on the studies of the influence of the aforementioned parameters on TMS, previous works [35,36,37,38,39,40,41] have analyzed the general mechanisms associated with TMS. From these works, it can be deduced that to optimize the excitation parameter for each TMS application, avoiding ethical conflicts, and for a practical and cost-effective estimation of the effect of the excitation electric field on neuronal activity, it is necessary to use simulation methods. In particular, Roth and Basser [42] proposed a method for computing the electric field and its gradient along a nerve induced by a single coil. They modeled the nerve as a passive cable whose electric potential is influenced by an excitatory electric field gradient. The excitatory electric field corresponds to the primary electric field [42,43,44], which depends only on the magnetic potential vector rate [45,46], without considering the effect of the scalar potential. This approximation was used in the some of initial TMS simulation works [42,43,44] and os still widely used but for coil parameter optimization [21,22]. However, it was later demonstrated [47,48,49,50,51] that a more precise estimation of the actual electric field induced in the brain tissues requires consideration of the secondary electric field equal to the negative gradient of the scalar potential. The scalar potential results from the electric charges formed at the tissue interfaces [47,48,49,50,51]. From this point, most of the works dedicated to the computation of the induced electric field [38,52,53,54,55,56,57] have proposed increasingly more realistic models for the coils, brain geometry and tissue properties and different formulations based on quasistatic approximation of the Maxwell equations [45,46]. In particular, quasistatic approximations have varied from magnetic quasi-static Laplace to magnetic quasistatic Poisson formulations, and with improvement of finite-element software, it was possible to include most of the source and potential contributions. However, the simultaneous use of complex coils, realistic brain geometry and complex formulations does not come for free, and even with the use of advanced software, the computational time increases substantially. This fact has brought forward questions such as: how much detail of the coil geometry is necessary to estimate the actual electric field induced by the real coil [58]? The development of new numerical methods such as neural networks [59] has also been motivated in order to increase the computation speed of TMS fields. Nevertheless, the most common approach has been the use of specific formulations of quasistatic approximation of Maxwell equations depending on the TMS applications, which is likely to remain the dominant approach for future TMS numerical computation, which, due to the microscopic size of neurons and nerves, requires increasing resolution and precision. Therefore, it is necessary to know the limitations and advantages of different formulations, as well as the different available numerical methods for their solutions.

2. Scope and Contributions

The TMS formulations analyzed in the present work correspond to the Maxwell equation potential representations and their quasistatic approximations, which are solved using numerical methods [1,2,3,4,5,6,7,8,9,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98]. Other methods such as the impedance method [99] and neural network methods [59] are not considered in this review. The methods that start from Maxwell equations but derive non-differential equation representations such as dipole-based methods [100] are also not considered in the present analysis.
The objective of this study is to serve as a comprehensive guide for the numerical computation of TMS fields by reviewing previous works in this field. To achieve this, a series of deductions leading to various quasistatic approximations for the computation of TMS fields are presented, followed by a review of works that have utilized these approximations. Furthermore, the numerical solutions using finite difference and finite element methods are discussed in detail, drawing upon existing methods proposed in the literature. Additionally, this study deduces methods not found in the consulted literature, such as finite difference for quasistatic magnetic A- ϕ and Darwin models and finite element for quasistatic magnetic A- ϕ and Darwin models using the Galerkin method, as well as their implementation of boundary conditions applied to TMS fields.
The papers included in this review pertaining to field simulations in TMS were systematically chosen in a chronological manner, with emphasis placed on commencing with highly impactful publications based on their citation count. Additionally, works that cited these influential papers while introducing novel variants of formulations and numerical methods were also considered. It is essential to acknowledge that not all the works encompassing TMS field simulations propose new formulations or numerical methods. In fact, a significant portion of the TMS field simulation literature references the original works that introduced these methods, with their contributions primarily concentrated in the domain of applications. Every effort was made to incorporate as many of these application-oriented works as feasible.

3. Field Theory of TMS

3.1. Maxwell Equations and Their Representation Using Vector and Scalar Potentials

To reduce the computational time and facilitate the numerical computation of TMS fields, some quasistatic approximations of the Maxwell equations [45,46] or, more specifically, their vector and scalar representation have been proposed [42,43,44,47,48,49,50,51]. Each of the formulations described below, some of which were used in previous TMS works, is derived from the potential representation of Maxwell equations. The classical forms of Maxwell equations are the Maxwell–Faraday, Maxwell–Ampere, Maxwell–Thompson and Maxwell–Gauss equations, expressed, in that order, as follows:
× E = B t
× H = J s + J i n d + D t
· B = 0
· D = ρ
The constitutive equations should also be considered:
D = ϵ E = ϵ o ϵ r E
B = μ H = μ o μ r H ,
as well as Ohm’s law
J = σ E
which establish the relation between the current density and the electric field. The Maxwell equations could be expressed using vector and scalar potentials [45,46]. Since the divergence of the curl of a vector is always zero, · × A = 0 . Then, B, the divergence of which is zero according to the Maxwell–Thompson equation ( · B = 0 ), could be replaced by the curl of an arbitrary vector such that:
B = × A
where A is the vector magnetic potential. Substituting this expression into the Maxwell–Faraday’s Equation (1a) and considering the distributive and commutative properties of the curl and the partial derivative, the following expression is obtained:
× E + A t = 0
According to Jackson [45], × ϕ = 0 for any scalar function as it is, for example, for the electric scalar potential ( ϕ ). Then, the term inside the curl can be expressed as E + A t = ϕ and resolving with respect to E :
E = A t ϕ
Substituting Equations (3) and (6) and the constitutive Equation (2) into the Maxwell–Ampere Equation (1b) results in:
× 1 μ × A + σ A t + ϕ + t ϵ A t + ϕ = J s
Applying the divergence operator to this equation and taking into account the aforementioned identity, i.e., that the divergence of the curl of a vector is always zero, gives:
· σ A t + ϕ + · t ϵ A t + ϕ = · J s
where J s is the external current density.
Using Equations (7) and (8), the full Maxwell equations are represented using vector and scalar potentials. However, the computational cost required to solve this set of differential equations could be significantly high depending on the geometry of the domains and their electric and magnetic properties. These equations could be simplified by making some considerations valid for TMS excitation characteristics and electric and magnetic properties of brain tissue. The formulations corresponding to this approximations are referred to in the literature as quasistatic approximations. The most well-known quasistatic approximations are [45,60,61] electro-quasistatic (EQS) approximation, which enables modeling of the potentials under the influence of conduction and displacement currents, including capacitive-resistive effects, but neglects the inductive phenomena; magnetic quasistatic (MQS) approximation, which neglects displacement current and capacitive effects and considers inductive phenomena; and the Darwin model, which comprises short-off electromagnetic quasistatic approximation and allows for the inclusion of inductive phenomena, as well as capacitive effects.

3.2. Poisson MQS ϕ -Formulation

The principle of TMS stimulation is that an electric field is induced in the brain by a time-varying magnetic field. The magnetic field is produced by a coil (or coils) fed by current peaks. Therefore, it is impossible to neglect inductive phenomena. Consequently, the behavior of potential can be computed using MQS. In MQS, the most restrictive formulation is the MQS- ϕ formulation, which is based on three approximations:
  • The wavelength of the excitation field is significantly higher than the size of the head. For TMS pulses with a duration of τ 0.1 ms corresponding to a frequency of f 10 kHz, the corresponding wavelength is λ 300 m, which is much higher than the head dimensions;
  • Diminishment of the capacitive effects in the brain tissue: Resulting from continuity conditions of electric currents in the interface between materials with different electric conductivities [45], the electric charges accumulate in this area, which could provoke capacitive effects. However, under this quasistatic approximation, induced charges are considered to move freely inside the brain, not allowing for static accumulation of charges that could generate capacitive effects. Furthermore, polarization effects are not considered;
  • Neglecting the skin-effect: A time-varying magnetic field induces electric currents that oppose the magnetic field. The amplitude of the induced current is proportional to the electric conductivity. Consequently, when a magnetic field enters a medium with an electric conductivity other than zero, it decays as it penetrates in the medium. However, the electric conductivity of brain tissue is low ( σ w h i t e m a t t e r = 0.15 S/m), and together with the paramagnetic magnetic properties of these tissues, δ = 1 / π f μ σ 5 m [62], which confirms this approximation.
To apply these approximations to Equations (7) and (8), it is important to explain the cause–effect relations in these equations. According to Equation (6), Ohm’s Law and the constitutive equation for the electric flux density, Equation (7) can be expressed as:
× 1 μ × A = J s + J i n d + J d i s p
where J i n d = σ E = σ A t + ϕ are the induced conductive currents or eddy currents, and J d i s p = D t = t ϵ E = ϵ A t + ϕ are the displacement currents. This expression shows that the magnetic field is induced by an external current whose value decreases due to the induced currents. The induced currents can be computed according to Equation (8), which can be rewritten as:
· J i n d + · J d i s p = · J s
The simplification of the MQS ϕ -formulation implies neglect of the effects of eddy currents on the magnetic field (no skin effect) and of the displacement currents (no capacitive effects). Consequently, this approximation implies removing J i n d and J d i s from Equation (9). This approximation does not mean that the induced current is zero but that J i n d < < J s . Furthermore, all the materials involved in the phenomenon are paramagnetic ( μ r = 1 ). Therefore, Equation (9) becomes:
× × A = μ 0 J s
Since the vector magnetic potential is not unique, a restriction such as the Coulomb gauge ( · A = 0 ) could be used [45], which, together with the identity × × A = · A 2 A , modifies Equation (9) as follows:
2 A = μ 0 J s
On the other hand, in the domain of the brain, neglecting the displacement current, J d i s = 0 , and since there is no external current, J s = 0 , and Equation (10) becomes:
· J i n d = 0
or
· σ A t + ϕ = 0
which could also be expressed as:
· σ ϕ = · σ A t
Expressions (12) and (15) could be used to obtain A ( r , t ) and ϕ ( r , t ) , which allow for computation of the TMS-induced electric field ( E ) using Equation (6). For this purpose, the electric field is divided in two components: E = E p + E s , where:
E p = A t
is called the primary electric field and
E s = ϕ
is the secondary electric field.
This formulation can be called the Poisson MQS ϕ -formulation, since both expressions are Poisson equations. This is also a ϕ formulation because A and ϕ are only partially coupled. The scalar electric potential ( ϕ ) is produced by electric charges induced by the time variation of the vector magnetic potential ( A ), but the influence of the electric field on the magnetic potential is neglected. This approximation has the advantage that it could be easily solved using numerical methods with low computational resources. However, it is necessary to consider the aforementioned limitations of the model.

3.3. Laplace MQS ϕ Formulation

A further simplification of the Poisson MQS formulation consists of considering that the electric charges responsible for the secondary electric field are formed only at the brain-to-air interface. The principle of TMS stimulation is that an electric field is induced in the brain by a time-varying magnetic field. The magnetic field is produced by a coil (or coils) fed by current peaks. The accumulation of charges at the brain-to-air interface could be explained by considering the following principles: The primary electric field ( E p ) is computed as if the medium were isotropic and uniform. Therefore, E p is continuous across all interfaces.
E p 1 · n = E p 2 · n = E p · n
According Ohm’s Law, J i p = σ i E p , and considering that electric conductivities of the interfacing areas are different, ( σ 1 σ 2 ), at the interfaces:
σ 1 E p · n σ 2 E p · n
Therefore,
J p 1 · n J p 2 · n
In the quasistatic limit, the normal component of total current crossing an interface between tissues is continuous, and it is given by the continuity law [45]:
J 1 · n = J 2 · n
which seems to contradict (20). However, this law is satisfied by the formation of electric charges, which generate a secondary electric field ( E s ) and the corresponding current ( J s i = σ i E s i ). In this case, Equation (21) can be expressed as:
J p 1 + J s 1 · n = J p 2 + J s 2 · n
or
σ 1 E p + E s 1 · n = σ 2 E p + E s 2 · n
The fact that the secondary fields in the media of the interface are different compensates for the total electric field in such a way that the continuity condition holds. In the case of the brain-to-air interface, σ a i r = σ 2 = 0 . Then, Equation (23) reduces to:
σ 1 E p + E s 1 · n = 0
or
σ 1 ϕ · n + σ 1 A t · n = 0
Since the charges and the corresponding electric currents are restricted to the brain boundary, which corresponds to the case when the brain model is made of a uniform tissue, then σ is considered to be uniform and isotropic inside the brain, and Equation (15) for the brain domain becomes:
2 ϕ = 0
This Laplace equation can be used to compute the electric potential inside the brain domain. To obtain a unique and non-trivial solution from this equation, the Neumann boundary condition expressed by Equation (25) should be used. After obtaining ϕ from Equation (26) with the boundary condition (25) and A from Equation (12), the electric field is computed using Equation (6).

3.4. MQS A- ϕ Formulation

If only capacitive effects are neglected ( J d i s = 0 ) using the general transformations and considerations described in Section 3.2, Equations (7) and (8) can be expressed as:
σ A t + ϕ + 1 μ 0 × × A = J s
· σ ϕ + · σ A t = 0
or using the Coulomb gauge:
σ A t + ϕ 1 μ 0 2 A = J s
· σ ϕ + · σ A t = 0
which is the transient representation of the MQS A- ϕ formulation. In this case, A and ϕ are coupled, and the effects of the conductively induced currents on A are considered. However, due to the complexity of numerical computation of the time-dependent functions ( A ( r , t ) and ϕ ( r , t ) ), it is more usual to use the harmonic expressions of Equations (28a) and (28b):
j ω σ A + σ ϕ + 1 μ 0 × × A = J s
j ω · σ A + · σ ϕ = 0
The harmonic formulation can be used considering that the stimulation current presents a constant frequency ( ω ). This approximation should be used carefully, since the stimulation peak patterns are far from having a harmonic waveform. It seems unnecessary and a waste of computational resources to consider the induction currents for TMS. However, this approach could be useful for optimizing the effect of electric conductivity and geometry of the coil on the TMS. It could also be used for analyzing TMS in subjects with non-normative brains, such as those with brain implants.

3.5. TMS Full Maxwell Equation Formulation and the Darwin Model

The numerical computation of the electric field can be performed without any quasistatic approximation using Equations (7) and (8), but again, the complexity of the numerical computation of these equations is substantial. Therefore, the harmonic formulation of these equations is more commonly used:
j ω σ ω 2 ϵ o ϵ r A + σ + j ω ϵ o ϵ r ϕ + 1 μ 0 × × A = J s
· j ω σ ω 2 ϵ o ϵ r A + · σ + j ω ϵ o ϵ r ϕ = · J s
Nevertheless, apart from the harmonic stimulus limitation, this formulation still presents high computation complexity. An alternative that could be used for transient process computation is the Darwin model [61]. The Darwin model is not a full Maxwell equation formulation but allows for consideration of almost all the coupled effects, except for wave propagation. Since the domain scale in TMS is significantly smaller than the wavelength of the excitation magnetic field, simulation of wave propagation is completely unnecessary. The Darwin model is based on Helmholtz decomposition [45], in which the electric field is decomposed into an irrotational component and a solenoid component: E = E i r r + E s o l , where × E i r r = 0 and · E s o l = 0 . After replacing the electric field in the Maxwell equations with these approximations, neglecting the radiation effects and replacing the expression of the potentials yields:
× 1 μ × A + σ A t + ϕ + t ϵ ϕ = J s
· σ A t + ϕ + · t ϵ ϕ = · J s
This formulation takes into account all the effects relevant for TMS included in the full Maxwell equation formulation, and it is less computationally complex than integrating the full Maxwell equations. Moreover, in this case, it seems excessive to include capacitive effects for TMS simulations, but it could be useful to optimize the stimulation devices and for analysis of TMS in subjects with non-normative brains, such as those with brain implants.

3.6. Boundary Value Problems of TMS

Figure 1 shows the domain and boundary representation of a generic TMS setup. The meaning of the domains, their boundaries, and electric and magnetic properties are shown in Table 1.
The coil is a special domain, as the source of external current, the electric conductivity of which depends on the type of quasistatic formulation. In this case, it is not the boundary that is relevant but the integration path used to compute the magnetic potential vector, as shown latter.
Based on these assumptions related to the domains, the boundary value problems(BVP) can be expressed as follows:
The Laplace MQS-ϕ BVP:
2 A = μ 0 J s on Ω
2 ϕ = 0 on Ω b
× A · n = 0 on Γ
ϕ · n = A t · n on Γ b
The Poisson MQS-ϕ BVP:
2 A = μ 0 J s on Ω
· σ ϕ = · σ A t on Ω b
× A · n = 0 on Γ
ϕ · n = 0 on Γ
The MQS A-ϕ BVP:
σ A t + ϕ 1 μ 0 2 A = J s on Ω
· σ ϕ = · σ A t on Ω b
× A · n = 0 on Γ
ϕ · n = 0 on Γ
The Darwin model BVP:
× 1 μ × A + σ A t + ϕ + t ϵ ϕ = J s on Ω
· σ A t + ϕ + · t ϵ ϕ = · J s on Ω
× A · n = 0 on Γ
ϕ · n = 0 on Γ
Several works have used the MQS ϕ formulation. Table 2 lists some of these works and their models:

4. Solutions of TMS Field BVPs

4.1. Solution of MQS- ϕ BVP

Various methods have been proposed to address the MQS- ϕ formulation. In the case of the Laplace MQS- ϕ formulation, one of the most common solution is, first, to compute the primary field, E p = A t , where A is obtained by integrating the Poisson Equation (32a) over the volume of the coil ( Ω c o i l ), which gives:
A ( r o , t ) = μ o 4 π Ω c o i l J s ( r , t ) r r o d V
where J s ( r , t ) is the current density. This equation is easy to solve numerically. For example, in the case of a circular coil with uniform current density ( J s ( t ) ), a constant cross-section area (S) and a current path in a plane perpendicular to the z axis in the domain ( Ω ), each component of the magnetic potential vector at the position of r = r o , A ( r o , t ) could be computed as as:
A ( r o , t ) = μ o S J s ( t ) 4 π Σ c o i l 1 r r o d l
where Σ c o i l is the path. Equation (37) can be numerically computed as follows:
A x ( r o , t ) = μ o S J s ( t ) 4 π i = 0 N 1 Δ x i ( x i x o ) 2 + ( y i y o ) 2 + ( z i z o ) 2
A y ( r o , t ) = μ o S J s ( t ) 4 π i = 0 N 1 Δ y i ( x i x o ) 2 + ( y i y o ) 2 + ( z i z o ) 2
A z ( r o , t ) = 0 , r o ϵ Ω
r = ( x i , y i , z i ) ϵ Ω , r r o
where
x i = X m + R cos ( ψ i )
y i = Y m + R sin ( ψ i )
and
Δ x i = x i + 1 x i , Δ y i = y i + 1 y i
where S, X m and Y m are the coil cross-sectional area, the coil center position on the x axis and the coil center position on the y axis, respectively. It is possible to set up different types of excitation coils by adding coils with different parameters, positions and sizes.
Once the primary field is obtained, the secondary field can be computed using the following equation: E s = ϕ , where ϕ is obtained after solving the Laplace Equation (32b) and the boundary condition (32d).
In the case of the Poisson MQS- ϕ formulation, A is computed using the same procedure as in the Laplace formulation, but ϕ is computed by solving the Poisson Equation (33b) with boundary conditions (33c) and (33d). This equation can be solved using numerical, analytical or semianalytical methods depending on the complexity of the brain model. Two of the most commonly used numerical methods are the finite difference method (FDM) [63,64,65,66,67] and the finite element method (FEM) [64,68,69,70].
Equation (33b) for the scalar potential can be expressed as:
· σ ϕ = σ · ϕ + σ · 2 ϕ = · J i n d
where J i n d = σ E = σ A t is the conductively induced current.The FDM to solve this Poisson equation can be implemented using the interactive method:
ϕ i , j , k = 1 α [ σ i 1 2 , j , k ϕ i 1 , j , k + σ i + 1 2 , j , k ϕ i + 1 , j , k Δ x 2 + σ i , j 1 2 , k ϕ i , j 1 , k + σ i , j + 1 2 , k ϕ i , j + 1 , k Δ y 2 + σ i , j , k 1 2 ϕ i , j , k 1 + σ i , j , k + 1 2 ϕ i , j , k + 1 ] Δ z 2 + 1 2 α J i + 1 , j , k x J i 1 , j , k x Δ x + J i , j + 1 , k y J i , j 1 , k y Δ y + J i , j , k + 1 z J i , j , k 1 z Δ z
where α = σ i 1 2 , j , k + σ i + 1 2 , j , k Δ x 2 + σ i , j 1 2 , k + σ i , j + 1 2 , k Δ y 2 + σ i , j , k 1 2 + σ i , j , k + 1 2 Δ z 2 and J i , j , k x = σ i , j , k E i , j , k x , J i , j , k y = σ i , j , k E i , j , k y , J i , j , k z = σ i , j , k E i , j , k z are the currents densities due to the primary electric field ( E p = ( E p x , E p y , E p z ) ), and Δ x , Δ y and Δ z are the grid intervals. This method is valid for σ 0 .
It is clear that in the case of the Laplace MQS-V formulation, the same approach as before could be used, except that J i n d = 0 ,
ϕ i , j , k = 1 6 ϕ i 1 , j , k + ϕ i + 1 , j , k + ϕ i , j 1 , k + ϕ i , j + 1 , k + ϕ i , j , k 1 + ϕ i , j , k + 1
This expression should be solved using the boundary conditions (32d), which can be expressed as:
ϕ x n x + ϕ y n y + ϕ z n z = E x p n x + E y p n y + E x p n z
which is solved using the finite difference as: ϕ i , j , k = 1 n x Δ x + n y Δ y + n z Δ z E x p + ϕ i 1 , j , k Δ x n x + E y p + ϕ i , j 1 , k Δ y n y + E z p + ϕ i , j , k 1 Δ z n z
It is also important to note that the above numerical solution corresponds to an Euclidean coordinate system, and consequently, the brain should be circumscribed inside an air cube, and the boundary conditions in this case correspond to an interior boundary. From a computational coding perspective, the FDM is straightforward. However, the precision and the computation time required by this method depend on the complexity of the geometry.

4.1.1. Solutions Using FEM

The FEM usually requires much less computational time than the FDM but can be more tricky to implement. The Poisson MQS- ϕ BVP can be solved with FEM using different approaches: the direct minimization of energy functional and the weighted residuals or the Galerkin method. For example, Weiping Wang and Solomon R. Eisenber [95] considered the following energy functional:
W ( ϕ ) = J · E d V ,
which represents the energy dissipated in the conduction media during the induction process. Using the expressions of the induced current and electric fields, the resultant equation is:
W ( ϕ ) = σ ϕ + A t · ϕ + A t d V
The objective is to find the potentials ( ϕ and A ) that minimize the energy functional, i.e., W ( ϕ ) . The potentials at any coordinate ( x , y , z ) inside the domain can be obtained from the potential at the nodes of the finite elements (usually tetrahedral elements), in which the domain is divided using interpolation functions:
Φ e ( x , y , z ) = i = 1 n n ϕ i e N i ( x , y , z )
A e ( x , y , z ) = i = 1 n n A e i N i ( x , y , z )
where n n is the number of nodes of the finite elements, n n = 4 represents tetrahedral elements, N i ( x , y , z ) are the interpolation functions and ϕ i e , A e i = ( A x ) i e i + ( A y ) i e j + ( A z ) i e k are the electric scalar potential and the spacial part of the magnetic vector potential, respectively, corresponding to node i in element e. In this case, variable separation for the magnetic potential vector ( A ( x , y , z , t ) = A ( x , y , z ) · f ( t ) ) is assumed. The interpolation functions are usually linear functions of the form N i ( x 1 , x 2 , , x n n ) = j = 1 n n a i j x j + d . For example, in the case of tetrahedral elements ( N i ( x , y , z ) = a i x + b i y + c i z + d ), the coefficients are obtained by solving a linear equations system, considering that N i ( x , y , z ) = 1 for x = x i , y = y i , z = z i and N i ( x , y , z ) = 0 for x = x j , y = y j , z = z j , where x i , y i , z i are the coordinates of node i, and x j , y j , z j are the coordinates of the other nodes of the element.
The integral corresponding to total dissipated power for each element given by Equation (46) is divided in two parts: a voltage-dependent power and a constant term:
W e ( Φ e ) = W N e ( Φ e ) + W c o n s t e
Considering that A t = f ˙ A ( x , y , z ) , f ˙ = f ( t ) t .
W N e ( Φ e ) = σ ( Φ e ) 2 + 2 σ f ˙ A · Φ e d V
W c o n s t e ( Φ e ) = ( f ˙ ) 2 σ ( A ) 2 d V
The constant term does not influence the minimum value of the functional and is therefore removed. The potential dependent term after replacing Equations (47a) and (47b) gives:
W N e ( Φ e ) = ( Φ e ) T P e Φ e + ( Φ e ) T Q e
where Φ e = ( ϕ 1 e , ϕ 2 e , ϕ 3 e , ϕ 4 e , , ϕ n n e ) T is the vector of scalar potentials corresponding to the nodes of the element (e), P e is a square matrix and Q e is a column vector given by:
P i j e = σ N i · N j d V
and
Q i e = 2 f ˙ j = 1 n n σ N i · A j N j d V
These integrals are very easy to compute, since the interpolation functions are linear functions. For example, in the case of tetrahedral elements, the integrated expressions are:
P i j e = V e Λ i x a j + Λ i y b j + Λ i z c j
and
Q i e = 2 f ˙ V e 2 j = 1 4 Λ i x A j x + Λ i y A j y + Λ i z A j z
where V e is the volume of the element, and Λ i k = σ k x a i + σ k y b i + σ k z c i , k = x , y , z . Equation (51) corresponds to the power dissipated in one element. The total power is obtained as the sum of the power of all elements:
W N ( Φ ) = e W N e ( Φ e )
Nevertheless, it is worth noting that to obtain this sum, the dimensions of the vectors ( Φ e and Q e ) and the matrix ( P e ) should be increased such that d i m ( Φ e ) = n n × N e , d i m ( Q e ) = n n × N e and d i m ( P e ) = ( n n × N e ) 2 , filling the rest of the components with zero and shifting the position of the vector or matrix according to the element index (e) and the number of nodes n n . For example, the new vector ( Φ e ) should look like:
Φ e = 0 , , 0 ( e 1 ) × n n , ϕ 1 e , ϕ 2 e , ϕ 3 e , ϕ 4 e , , ϕ n n e , 0 , , 0 T
Similarly to the vector ( Q e ), the matrix ( P e ) should look like:
P e = 0 0 ( e 1 ) × n n + 1 0 0 0 P 11 e P 14 e 0 0 P 41 e P 44 e 0 0 0         ( e 1 ) × n n + 1
In the new space, which is called disjoint space, Equation (56) can be written as:
W d i s ( Φ d i s ) = Φ d i s T P d i s Φ d i s + Φ d i s T Q d i s
The potentials in the disjointed space are represented by variables refereed to each element, which gives a large number of variables. The number of variables can be reduced significantly by considering that neighbor elements share nodes, and therefore, it is possible to use the same variable for the potential of the nodes shared by different elements. In other words, the variables can be renamed using a common global notation for all the nodes inside the domain. This process is called assembly and can be performed using a connectivity matrix ( C ) such that:
Φ d i s = C Φ a s e m
Substituting this equation into Equation (69) gives:
W a s e m ( Φ a s e m ) = Φ a s e m T P Φ a s e m + Φ a s e m T Q
where P = C P d i s C T and Q = C Q d i s .
The linear equation system used to obtain the potentials at each node ( Φ i ) is formed by finding the minimum of the energy with respect to each node potential:
W ( Φ ) Φ i = 0
which results in the following linear equation system:
2 P Φ + Q = 0
Furthermore, to obtain a unique solution, it is necessary to impose an additional condition, which is usually achieved nu setting the potential of the Nth node to zero ( Φ N = 0 ). Thus, the equation system is updated by removing the Nth row and column of matrix P and the last element of Q, resulting in a matrix ( P ( d i m ( P ) = ( N 1 ) × ( N 1 ) )) and a vector ( Q ( d i m ( Q ) = 1 × ( N 1 ) )). The solution for the equations system is:
Φ = 1 2 P 1 Q
After computing the scalar potential, it is possible to obtain the electric field using Equation (6).

4.1.2. Weighted Residual Galerkin Method

The Weighted residual Galerkin method [68,69] is one of the most widely used FEMs. In this case, the starting point is Equation (32), which considers the domain of a single element ( Ω e ), and the weak formulation [68,69] and the residual for the element (e) can be written as:
r e = · σ ϕ + · σ A t
If the numerical solution were exact, the residual would be zero. However, since the numerical solutions are not ideal, it is necessary to find a solution that minimizes the residual or, in the case of this method, the weighted residual. The procedure involves multiplying the residual by a weight (w), integrating this product over the volume of the element and setting the integral to zero.
Ω e w r e d V = Ω e w · σ ϕ + · σ A t d V = 0
using the identity · w A = w · A + A · ( w ) , which implies that:
w · A = · w A A · ( w )
Using this identity, Equation (67) can be modified as:
Ω e · w σ ϕ d V Ω e σ ϕ · w d V + Ω e w · σ A t d V = 0
Applying the divergence theorem [45] ( Ω e · F d V = Γ e F · n d S ) to the first term of Equation (68) and rearranging the equation yields:
Ω e σ ϕ · w d V + Ω e w · σ A t d V + Γ e w σ ϕ · n d S = 0
The three components of this equation are the potential-dependent term, the source term and the boundary condition term, in that order. The source term can be modified by considering the induced current density ( J = σ · A t ):
Ω e σ ϕ · w d V Ω e w · J d V + Γ e w σ ϕ · n d S = 0
According to the Garlerkin formulation [68], the weight function (w) should be set equal to the interpolation functions:
w = N i e ( x , y , z ) = a i x + b i y + c i z + d
Moreover, using Equation (47a) for the interpolation of the scalar potential ( ϕ e ) and considering the interpolation for the induced current density values:
J e ( x , y , z ) = i = 1 n n J i e N i ( x , y , z )
Equation (70) becomes a system of n n equations, where n n is the number of element nodes:
P e Φ e + Q e J e + B e = 0
where
P i j e = Ω e σ x x e N i x N j x + σ y y e N i y N j y + σ z z e N i z N j d z d V
Q i j e = Ω e N i N j x + N j y + N j z d V
B i e = Γ e N i σ x x e Φ x n x + σ y y e Φ y n y + σ z z e Φ z n z d V
d i m ( P e ) = n n × n n , d i m ( Q e ) = n n × n n , d i m ( J e ) = 1 × n n , d i m ( B e ) = 1 × n n ,
Due to the fact that the normal vector of the faces shared by two neighbor elements has opposite directions [68,69], the integral B e cancels for interior element boundaries, and B e is different from zero only for the outer boundaries. Furthermore, if the simulation is performed considering an air box circumscribing the brain, in the air, σ = 0 ; therefore, B e ( x , y , z ) = 0 , ( x , y , z ) ϵ Ω s y s . Accordingly, this term can be removed from the equation, and the equation system becomes:
P e Φ e + Q e J e = 0
On the other hand, if the brain is considered the outer boundary and the approximation of having only one tissue interface, for example, the interface between the brain and the air, the second term can also be modified using the divergence theorem, which gives:
Ω e σ ϕ · w d V + Γ e w σ A t n d S + Γ e w σ ϕ · n d S = 0
or
Ω e σ ϕ · w d V + Γ e w σ A t + σ ϕ · n d S = 0
Considering the boundary condition corresponding to Equation (32d), it is evident that the second term is zero. Therefore, for this particular case, Equation (77) becomes:
Ω e σ ϕ · w d V = 0
which gives the following linear equation system:
P Φ = 0
The trivial solution ( Φ = 0 ) is avoided, considering the boundary condition expressed by Equation (32d) for boundary nodes. This can be implemented in several forms, for example, by considering the following equation system that results from condition Equation (32d):
Φ x = A x t = E x p
Φ y = E y p
Φ z = E z p
The value of Φ for the boundary nodes can be obtained from these equations for each element separately, forming a 3 × 3 linear equation system:
i = 1 3 Φ i e N i s x = ( E x p ) e
i = 1 3 Φ i e N i s y = ( E y p ) e
i = 1 3 Φ i e N i s z = ( E z p ) e
where N i s represents the interpolation functions for the potential of the nodes forming the triangular face, that is, facing the outer boundary, and ( E x p ) e , ( E y p ) e , ( E z p ) e are the components of the primary electric field in the boundary space.
Equation (75) corresponds to one element. To solve the global equation system, it is necessary to assembly the matrices for the entire domain ( Ω ). This can be done again by applying a procedure similar to that explained in the previous section. The dimensions of the matrices ( P e and Q e ) and the vectors ( J e and Φ e ) are extended as shown in Equations (57) and (58), generating a system of linear equations in the disjointed space such that:
P dis = e ϵ Ω P e
Q dis = e ϵ Ω Q e
and
J dis = e ϵ Ω J e .
Therefore,
P dis Φ d i s + Q dis J d i s = 0
There is also a redundancy of variables of the disjointed space due to the presence of shared nodes between neighbor elements. To remove this redundancy, it is necessary to assemble the matrices using an assembly process as in the previous section, which can be achieved using the connectivity matrix:
Φ d i s = C Φ a s e m
J d i s = C J a s e m
and substituting in (85):
P asem Φ a s e m + K a s e m = 0
where P a s e m = P d i s C and K a s e m = Q d i s C J a s e m , and the potential can be obtained from:
Φ a s e m = P asem 1 K a s e m

4.2. Solution of MQS Φ -V BVP

4.2.1. Solution Using FDM

The FDM can be applied to solve the BVP given by Equation (34) using several approaches. One of the most straightforward methods is the explicit interactive method [63], which can also be applied using different considerations. For example, in the case of isotropic conductivity, Equation (34a) can be divided in three equations corresponding to each of three coordinates in the Euclidean space. If Cartesian coordinates are used:
σ A x t + ϕ x x 1 μ 0 2 A x = J x s
σ A y t + ϕ y y 1 μ 0 2 A y = J y s
σ A z t + ϕ z z 1 μ 0 2 A z = J z s
The procedure to obtain the electric field is described as follows:
  • For t = 0 , J = 0 ; therefore, Φ , A = 0 ;
  • For t = t 1 , A is obtained from Equation (90) using the explicit interactive method in an isotropic grid ( Δ x = Δ y = Δ z = Δ h ) as follows:
    σ i , j , k ( A x t ) i , j , k ( A x t Δ t ) i , j , k Δ t 1 μ 0 ( A x t ) i 1 , j , k + ( A x t ) i + 1 , j , k + ( A x t ) i , j 1 , k ( Δ h ) 2 + 6 μ 0 ( A x t ) i , j , k ( Δ h ) 2 1 μ 0 ( A x t ) i , j + 1 , k + ( A x t ) i , j , k 1 + ( A x t ) i , j , k + 1 ( Δ h ) 2 + σ i , j , k ( Φ x t Δ t ) i + 1 , j , k ( Φ x t Δ t ) i 1 , j , k 2 Δ h = ( J x s ) i , j , k
    which, resolved with respect to ( A x t ) i , j , k , gives:
    ( A x t ) i , j , k = 1 α [ β ( A x t ) i 1 , j , k + ( A x t ) i + 1 , j , k + ( A x t ) i , j 1 , k + ( A x t ) i , j + 1 , k + ( A x t ) i , j , k 1 + ( A x t ) i , j , k + 1 ] + γ α ( Φ x t Δ t ) i + 1 , j , k ( Φ x t Δ t ) i 1 , j , k + θ α ( A x t Δ t ) i , j , k + 1 α ( J x s ) i , j , k
    where α = σ i , j , k Δ t + 6 μ o ( Δ h ) 2 , β = 1 μ o ( Δ h ) 2 , γ = σ i , j , k 2 Δ h and θ = σ i , j , k Δ t . This equation is applied recursively to all nodes, except for the boundary nodes (where the solution is already given) and given the initial condition ( A = 0 ) until the solution converges. The obtained value of A is replaced in Equation (41), which can be resolved for Φ using Equation (42);
  • The same procedure is applied to compute A y and A z ;
  • The values of A and Φ are used to compute the electric field using Equation (6). For the next time instant, the value of Φ is replaced in step 2, and the process is repeated.

4.2.2. Solution Using FEM (Galerkin Method)

The FEM solution can also be obtained using different approaches. Considering that electric conductivity is isotropic, Equation (90) can be represented as a weighted residual as follows:
Ω e w σ A x t + ϕ x 1 μ 0 2 A x J x s d V = 0
or
Ω e w σ A x t + w σ ϕ x 1 μ 0 w 2 A x w J x s d V = 0
which, applying the identity a 2 b = · a b a · b to the third term, becomes:
Ω e w σ A x t + w σ ϕ x 1 μ 0 · w A x + 1 μ 0 A x · w w J x s d V = 0
Using the divergence theorem as in the previous sections applied to the third term and rearranging the equation yields:
1 μ 0 Ω e A x · w d V + Ω e w σ A x t + w σ ϕ x d V Ω e w J x s d V 1 μ 0 Γ e w A x · n d S = 0
The potential at any position ( ( x , y , z ) ϵ Ω s y s ) is computed using the potential at the nodes and the interpolation functions:
A x e ( x , y , z ) = i = 1 n n ( A x e ) i N i ( x , y , z )
Φ e ( x , y , z ) = i = 1 n n ( Φ e ) i N i ( x , y , z )
Replacing the weight function, (w) with the interpolation functions in Equation (96) and A x , ϕ with expression (97) gives:
P e ( A x t ) e + T e ( A x t ) e + R e Φ + T e ( A x t Δ t ) e + F e + B e = 0
where P i j e = 1 μ o Ω e N i x N j x + N i y N j y + N i z N j d z d V , T i j e = σ e Δ t Ω e N i N j d V , R i j e = σ e Ω e N i x N j d V , F x e = ( J x s ) e Ω e N i d V and B x e = 1 μ o Ω e N i A x x n x + A x y n y + A x z n z d V is the boundary condition as in the previous section.
Equation (98) corresponds to one element and can be expanded for the disjoint space using the same procedure explained in the previous two sections, after which the equation system is assembled using the connectivity matrix. The resulting equation system considering all components of A and the potential ( Φ ) is:
P + T A x t + R Φ = T A x t Δ t S x
P + T A y t + R Φ = T A y t Δ t S y
P + T A z t + R Φ = T A z t Δ t S z
P Φ = Q J
or
P + T R 0 P A t Φ = T A t Δ t S Q J
where S v = F v + B v for v = { x , y , z } , J is a row of vectors formed by the values of induced currents for each of the nodes of each element ( J i = σ i A x i t + A y i t + A y i t and P i j = μ o σ i P i j ).

4.2.3. Solution of the Darwin Model BVP

It can be seen from Equations (35a) and (35b) or their versions considering the Coulomb gauge that they can be integrated using the procedure explained in the two previous sections by adding a new source term to Equations (34a) and (34b) and considering the change of potential due to the presence of displacement currents in the second equation. Considering the Coulomb gauge and that the electric permittivity ( ϵ ) is time-independent, Equations (35a) and (35b) can be expressed as:
1 μ o 2 A + σ A t + ϕ + ϵ ϕ t = J s
· σ A t + ϕ + · ϵ ϕ t = · J s
where the new source term for the first equation is ϵ ϕ t , and the new scalar potential term in the second equation is · ϵ ϕ t . Equations (101a) and (101b) can be further modified by approximating the scalar potential using a finite difference scheme:
1 μ o 2 A + σ A t + ϕ t + ϵ Δ t ϕ t ϵ Δ t ϕ t Δ t = J s
· σ A t + ϕ + · ϵ Δ t ϕ t · ϵ Δ t ϕ t Δ t = · J s
The term J d i s = ϵ Δ t ϕ t Δ t is a displacement current density that comes from the potential induced in the previous time instant. The term σ A t is the conductive current density ( J i n d ), which is magnetically induced. Therefore, these equations can be rearranged as:
σ A t 1 μ o 2 A + σ + ϵ Δ t ϕ t = J s + J d i s
· σ + ϵ Δ t ϕ t = · J s + J d i s + J i n d
which essentially have the same form as the A- Φ formulation, with an additional current source. Therefore, this equation set can be solved using the finite difference method as described in Section 3.1.

4.2.4. Solution of the Darwin Model BVP Using FEM

Houssein Taha [60] proposed the solution of the Darwin model directly using Equations (35a) and (35b) applied to electric machines. However, working with the numerical curl of vectors usually requires a high computational cost if applied to a complex geometry mesh such as the brain cortex. Based on the modification of the equations described in the previous section, the Darwin model for TMS applications can be reduced to a type of Φ -V formulation. Therefore, comparing Equations (34a) and (34b) and their FEM numerical solution (Equation (100)) with Equation (103), the FEM formulation for the Darwin model can be expressed as:
P + T R 0 P A t Φ = T A t Δ t S Q J
where P , T and Q are the same as in the Φ -V formulation solution, and R i j = β σ ϵ i Ω N i x N j d V , where β σ ϵ i = σ i + ϵ i Δ t , is constant for the nodes of each element. P i j = μ o β σ ϵ i P i j , S v = F v + B v for v = { x , y , z } , with F = ( J s + J d i s ) Ω N i d V and J = J s + J d i s .

4.3. Analytical and Semianalytical Solutions

It is clear that it is very difficult—if not impossible—to obtain an analytic solution of the TMS field formulations on realistic brain cortex geometries. However, it is possible to obtain analytic solutions considering simplified geometric models of the excitation coil and the brain. The advantage of these solutions is that they serve as a calibration reference for numerical simulations. Esselle et al. [49] computed the electric field produced by an arbitrarily shaped coil above half-space tissue with a planar interface with the air. The formulation corresponds to the Laplace MQS- ϕ formulation:
2 A = μ o J s , / A Ω , × A Γ = 0
2 ϕ = 0 , / ϕ Ω b , n 1 · J 1 Γ = n 2 · J 2 Γ
The solution of the magnetic vector potential generated by a coil formed by N thin wires carrying a current (I) is:
A ( r o ) = μ o N I 4 π Σ c o i l 1 r r o d l
The primary electric field can be computed from the magnetic potential vector; therefore:
d E p = μ o N ( d I / d t ) 4 π r r o d l
The secondary electric field ( E s ) of a coil element ( d l ) is computed from the scalar potential as:
d E s = ϕ = ϕ x i ϕ y j ϕ z k
Consequently, the total electric field of the coil element ( d l ) is given as:
d E = d E p + d E s = d E p x ϕ x i + d E p y ϕ y j + d E p z ϕ z k
Considering the condition of the continuity of the current in the interfaces, since the electric conductivity of the air is zero, then the current in the air–tissue interface produced by the d l element of the coil perpendicular to the tissue interface is in the z direction; therefore:
J b = σ b d E = σ b d E z = 0
or
d E z = 0
Therefore, it is only necessary to compute d E x and d E y . Furthermore, the boundary condition Equation (111) can be used to compute these components. From the boundary condition ( d E z = 0 ), it is deduced that at the interface is z = 0 :
μ o N ( d I / d t ) 4 π r r o d z ϕ z = 0
A general solution of the Laplace Equation (105b) can also be expressed as the Bessel integral:
ϕ = 0 U ( λ ) e λ z J o ( λ ρ ) d λ
where ρ = x 2 + y 2 , and J o is the Bessel function of the first kind. Replacing Equation (113) in Equation (112) and also using the expansion of 1 r r o as a function of the Bessel function ( J o ), the value of ϕ is obtained as:
ϕ = μ o N ( d I / d t ) 4 π d l z 0 1 λ e λ ( z o z ) J o ( λ ρ ) d λ
This expression is replaced in Equation (109) to obtain components d E x and d E y , which are given by:
d E x = d E p x ϕ x = μ o N ( d I / d t ) 4 π d l x R + x d l z ρ 2 1 z z o R
and
d E y = d E p y ϕ y = μ o N ( d I / d t ) 4 π d l y R + y d l z ρ 2 1 z z o R
The differential of the total electric field is obtained as d E = d E x i + d E y j .
Although the work by Esselle et al. [49] starts from the approximation that the tissue interface is plane, they also conclude that the result was independent of the tissue inhomogeneities inside the tissue half-plane if the conductivity of the inhomogeneous tissue changes perpendicular to the location of the interface.
A more precise geometrical model was used in the H. Eaton model [48], which proposed an analytical expression to compute the electric field and current density induced by a coil inside a homogeneous spherical coil. The author used the Laplace MSQ-V formulation as described in Section 3.3. Therefore, the solution of the magnetic potential vector is given by Equation (36). To analytically integrate this expression, the inverse of the distance term is expanded in terms of spherical harmonics:
1 r r o = 4 π l = 0 m = l l r l 2 l + 1 r o l + 1 Y l m * ( θ o , ϕ o ) Y l m ( θ , ϕ ) , r o > r
where Y l m ( θ , ϕ ) represents the spherical harmonic functions, and * indicates a complex conjugate. After replacing this expression in Equation (36), the value of the magnetic potential vector is:
A = μ o I l = 0 m = l l r l 2 l + 1 C l m Y l m ( θ , ϕ )
where
C l m = Ω c o i l Y l m * ( θ o , ϕ o ) r o l + 1 d l o
The values of the coefficient ( C l m ) are obtained by dividing the domain ( Ω b ) into rectangular components such that:
C l m = C l m x i + C l m y j + C l m z k
Moreover, ϕ is obtained after integrating the Laplace equation ( 2 ϕ = 0 ), which is the other equation of the the Laplace MQS- ϕ formulation, with boundary conditions given by the current continuity Equation (32d) and the electric potential outside the brain (sphere of radius (R)), vanishing ( ϕ = 0 ) for r :
ϕ = l = 0 m = l l F l m + G l m Y l m ( θ , ϕ )
where
G l m = F l m R 2 l + 1
F l m = j ω μ o I σ + j ω ( ϵ ϵ o ) l σ + j ω ϵ + j ω ϵ o l + 1 α l m D l 1 , m 1 + β l m E l 1 , m + 1 + γ l m C l m z
with α l m = l + m 1 l + m 2 l + 1 2 l 1 , β l m = l m 1 l m 2 l + 1 2 l 1 , γ l m = l m l + m 2 l + 1 2 l 1 , D l m = C l m x j C l m y 2 and E l m = C l m x + j C l m y 2 for l > 0 and F o o = 0 . Until this point, it is possible to obtain a semianalytical solution of the formulation by numerically integrating Equation (119) and computing A and ϕ using Equations (118) and (121), respectively, and the electric field using Equation (6). However, the authors took a step forward by replacing the analytic expressions with their corresponding spherical harmonic functions and obtaining analytical expressions for the electric field component in spherical coordinates, as described in [48].

5. Implementation of Numerical Solutions

Review of Solutions and Implementations Presented in the Literature

Table 3 provides a comprehensive overview of the predominant BVPs employed in the simulation of TMS fields. Alongside these BVPs, we present the corresponding Simulation Tools and their respective authors who have utilized them. The term "Direct implementation" pertains to cases where the authors have custom-built their simulation codes using programming languages, tailoring them to the specific requirements of their studies. On the other hand, "Commercial general-purpose software" refers to computer simulation programs, frequently employing the finite element method, capable of simulating various multiphysics problems represented by partial differential equations and boundary conditions.

6. Overview of Some of the Main TMS Simulation Results

The effectiveness and applications of TMS depend on several parameters of the coil-brain system, as mentioned in the Introduction. Some analyses using simulation methods have revealed how these parameters influence the electric field acting on the nerves. The use of different formulations, brain models and levels of numerical method precision can also influence the accuracy of estimating the field. A brief overview of some of these results is presented below.
Although it is difficult to establish a level of priority of elements that influence the electric field induced in the brain by TMS and the computation of the actual field distribution in the brain, there seems to be a consensus among works [16,18,20,21,22,23,24] related to the simulation of the electric field induced by TMS, which analyzed the influence of the coil parameters (geometry, positioning, current waveform, etc.). Their main conclusion is that the adequate selection of these parameters is crucial for inducing the value of the electric field. In other words, the excitation parameters are the most fundamental elements for achieving the necessary value and distribution of the electric field in the brain, and the smallest variation of these parameters could have a significant impact on the excitation field.
Research conducted by Mills [20] demonstrated that the utilization of single coils for the purpose of stimulating long nerves is not advantageous. This is attributed to the presence of “hot spots” within the induced electric field generated by such coils, which possess both positive and negative fields. These hot spots can elicit contradictory effects on the nerve. Additionally, through the simulation of various coil geometries and configurations, Deng [17] established that figure-8-type coils exhibit greater focal properties compared than circular-type coils. None of the simulated coil configurations surpassed the figure-8-type coils in terms of achieving a tradeoff between depth of penetration and focal precision. Furthermore, Deng’s research [17] indicated that an increase in the number of coil turns (N) amplifies coil inductance and peak voltage while reducing the necessary peak current. Theoretically, augmenting N enhances the accuracy of field replication, although this is constrained by the coil’s dimensional requirements.
The impact of TMS is heavily influenced by the parameters employed for excitation. Specifically, the electric field distribution, which crucial for effective stimulation, may not always achieve the desired level of focus due to the configuration of the excitation coil in TMS. Consequently, achieving highly localized stimulation through the appropriate selection and arrangement of coils can prove challenging for specific applications. In this regard, a novel emerging technique called transcranial ultrasound stimulation (TUS) [101,102] offers several advantages, including the ability to provide highly focused excitation using an ultrasonic transducer. Moreover, the utilization of ultrasonic waves allows for deeper penetration of the stimulus into the brain, enabling access to central brain regions. Nevertheless, despite being considered a safe technique, further research is necessary to establish the long-term safety and potential risks associated with TUS. Additionally, the use of coupling gel is required for excitation in TUS. Furthermore, the indirect mechanism by which ultrasound stimulates neural activity through mechanical oscillations is more intricate compared to the direct electrical stimulation employed by TMS.
In reference to the precision of the estimation of the induced electric field, the accuracy of the coil and current parameters used in the model is very important, as well as an adequate formulation and BVP [20,23,24]. In order of significance, the correct assignment of electric conductivity to each point in the brain model [30,31,32,33,34] and the selection of an appropriate numerical method for the simulations [73,90] are the subsequent key factors to be considered. Lastly, there are the details of brain model geometry, provided the use of a basic brain model including most of the fundamental parts (gyrus and sulcus) of the brain is used, although there is no clear consensus about the importance of the later. In this review, the effect of the shape and length of neurons and nerves, [103,104,105] which is also of fundamental importance, is not considered.
Petrov’s study [58] reported that in the specific region where neurostimulation is typically performed, spiral wire coil models exhibited higher precision (RE < 5%) in predicting the actual induced field compared to a single circular wire coil model, which deviated from field measurements by up to 10% RE. These findings emphasize the necessity of employing a spiral winding-turns model to accurately approximate the induced field of a typical TMS coil [58]. Moreover, the utilization of realistic magnetic resonance imaging (MRI)-derived head models demonstrated that differences in tissue conductivities had a negligible impact on the magnitudes of the electric field (E-field). However, the position and orientation of the coil, as well as the size of the brain, exerted significant effects on the magnitude of the E-field [24]. Notably, substantial deviations in E-field magnitude were observed with respect to coil placement when considering the evaluation of effects over a 2 cm range in each direction of a two-dimensional plane of the TMS coil [24]. This estimation may be regarded as conservative when assessing the disparity between the scalp location and the intended cortical target. The results highlight the extensive cortical area that may be affected when accounting for this level of positioning uncertainty, encompassing a significant portion of the dorsolateral prefrontal cortex [24]. Certainly, it is feasible to reduce this uncertainty within an individual patient by developing reliable coil placement procedures combined with neuronavigation techniques [24,58].
On the other hand, in a recent study [83], it was discovered that the standard deviations in maximum electric field values resulting from uncertainties in tissue conductivities accounted for approximately 5% of the mean value in TMS. Tissue conductivity emerged as the primary factor determining the magnitudes of current density when the displacement current was insignificant [83]. Conversely, with the presence of displacement currents, the permittivity became the principal determinant of current density magnitude. Furthermore, the study revealed that the existence of displacement currents could potentially increase the maximum cortical current density by an order of magnitude, assuming the extreme permittivity values reported by certain researchers prove to be accurate [83]. Notably, all solutions exhibited currents that were perpendicular to the cortical interface, challenging models that propose preferential stimulation of tangentially oriented neurons based on the assumption that fields normal to the interface are negligible. Consequently, these models warrant re-examination [83].
The impact of brain models on simulation accuracy remains a contentious topic. Salvador et al. [96] demonstrated that the most consistent correspondence between the estimated electric field and motor activation was observed in realistic brain models, specifically in the crown region of the precentral gyrus. Conversely, the use of simplified head models, such as spherical head models, yielded spatially nonspecific findings [96]. This study also proposed that the maximum electric field strength serves as the optimal parameter for exploring the applicability of field calculations in quantitative dosing [62,96]. However, Janssen et al. [89] showed that in a highly realistic head model, variations in sulcus width (up to 3 mm) did not result in significant differences in the calculated electric field values for most brain areas. Thus, for a global estimation of the electric field, an accurate representation of sulci is of limited importance [89]. Nevertheless, considerable overestimation of sulcus width led to an overestimation of local field strength in locations distant from the cortical hot spot. In contrast, Saturnino et al. [85] demonstrated that accurate anatomical representation of tissue boundaries, including sulci, significantly improved the numerical accuracy of TMS field estimation in full-head models and sulcus models. Specifically, the sulcus model highlighted the importance of higher mesh density around highly curved regions of the gray matter/cerebrospinal fluid boundary, where electric potential undergoes rapid changes. This approach helped prevent local simulation errors and emphasized that the under-representation of sulci could lead to substantial errors in the electric field [85]. Moreover, modifications to cortical geometry were shown to disrupt stimulating fields, potentially impairing cortex targeting in non-normal populations [83].
Regarding numerical methods employed in simulations, it was found that the finite difference method (FDM) requires high spatial resolution in regular grids or hexahedral meshes to achieve detail comparable to that of tetrahedral meshes. However, this results in a large number of elements, necessitating lengthy computation times and significant computational resources. Therefore, the finite element method (FEM) based on tetrahedral meshes was suggested as better-suited for these calculations [90].

7. Conclusions

Herein, we present a survey of the most used formulations, boundary value problems and implementations for the numerical computation of TMS fields. These aspects were analyzed, showing theirs limitations and advantages. The deduction of the formulations from the Maxwell equations and the numerical solution of the corresponding BVP were described in detail. We also induced some numerical solutions not found in the consulted literature, such as finite difference for quasistatic magnetic A- ϕ and Darwin models and finite element for quasistatic magnetic A- ϕ and Darwin models using the Galerkin method, as well as their implementation of boundary conditions applied to TMS fields. The aim of the present work was to serve as a guide for reproduction and future implementation of TMS simulations.
Several of the main results deduced from the TMS fields simulations emphasize the importance of the stimulation parameters, in particular, the coil parameters and arrangement, followed by the precision of the specification of the electrical properties of the tissue and details of the brain geometry model. Our literature analysis highlights that the most used approach is the Poisson MQS- ϕ formulation applied to a model consisting a figure-8-type coil and a brain model. The advantage of this method is that it does not require large computational resources and enables consideration of the effects of the formation of electric charges in the inner interfaces of the brain. In this case, computational resources can also be redirected to increase the detail of the brain model geometry. Laplace MQS approximation requires significantly fewer computational resources. However, apart from its limitation in considering the charges formed at the inner interfaces, it is practically limited in terms of setting the boundary condition in the complex brain-to-air interface, which usually requires an advance graphical interface. The use of commercial software is common practice in the reviewed works. Nevertheless, the use of software has its drawbacks, as it limits the flexibility of connection in the simulated electric field output with the input for neuron equations and models. The use of a direct implementation could solve these limitations, but there is still a considerable range of code implementations and source code styles, which limits the widening of their use. The recent development of TMS custom platforms such as SimNIBS represents a step forward in solving these drawbacks, although the tendency of the development of graphical interfaces that should be mastered could affect the effectiveness of such platforms. In that sense, it may be more useful to work towards the development of flexible toolboxes and libraries for TMS simulations.

Author Contributions

Conceptualization, J.A.P.-B. and P.M.-O.; methodology, J.A.P.-B.; writing and original draft preparation, J.A.P.-B., P.M.-O. and J.A.-M.; writing, review and editing, P.M.-O. and J.A.P.-B.; funding acquisition, J.A.-M. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The comments provided by the reviewers are deeply appreciated.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

for Maxwell Equations
E Electric field
D Electric flux density
H Magnetic field
B Magnetic flux density
A Magnetic potential vector
ϕ Electric scalar potential
ρ Electric charge density
J s Vector of external (excitation) current density
J i n d Vector of induced current density
J d i s p Vector of displacement current density
σ Electric conductivity
ε 0 Electric permittivity of free space
ε Material electric permittivity
μ 0 Magnetic permeability of free space
μ Material magnetic permeability
ω Frequency of the excitation current
j = 1 Complex number imaginary unit
for Numerical Methods
ϕ i , j , k Electric scalar potential
Φ e Electric scalar potential vector of element e
Φ Electric scalar potential vector for all elements in the domain Γ
J i , j , k k Nodal current density with the component k = { x , y , z }
J e Current density vector of the element e
JCurrent density vector of all the elements in the domain Γ
A i , j , k k         Nodal magnetic potential of the component k = { x , y , z }
A e Magnetic potential vector for the element e
A k Magnetic potential k-component for all elements in the domain Γ

References

  1. Merton, P.A.; Morton, H.B. Stimulation of the cerebral cortex in the intact human subject. Nature 1980, 285, 227. [Google Scholar] [CrossRef]
  2. Barker, A.T.; Jalinous, R.; Freeston, I.L. Non-invasive magnetic stimulation of human motor cortex. Lancet 1985, 1, 1106–1107. [Google Scholar] [CrossRef] [PubMed]
  3. Kobayashi, M.; Pascual-Leone, A. Transcranial magnetic stimulation in neurology. Lancet Neurol. 2003, 2, 145–156. [Google Scholar] [CrossRef] [PubMed]
  4. Paulus, W.; Peterchev, A.V.; Ridding, M. Transcranial electric and magnetic stimulation: Technique and paradigms. Handb. Clin. Neurol. 2013, 116, 329–342. [Google Scholar] [PubMed]
  5. Gershon, A.A.; Dannon, P.N.; Grunhaus, L. Transcranial magnetic stimulation in the treatment of depression. Am. J. Psychiatry 2003, 160, 835–845. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Simons, W.; Dierick, M. Transcranial magnetic stimulation as a therapeutic tool in psychiatry. World J. Biol. Psychiatry 2005, 6, 6–25. [Google Scholar] [CrossRef]
  7. Hallett, M. Transcranial magnetic stimulation and the human brain. Nature 2000, 406, 147–150. [Google Scholar] [CrossRef]
  8. Burt, T.; Lisanby, S.H.; Sackeim, H.A. Neuropsychiatric applications of transcranial magnetic stimulation: A meta analysis. Int. J. Neuropsychopharmacol. 2002, 5, 73–103. [Google Scholar] [CrossRef] [Green Version]
  9. Kahkonen, S.; Ilmoniemi, R.J. Transcranial magnetic stimulation: Applications for neuropsychopharmacology. J. Psychopharmacol. 2004, 18, 257–261. [Google Scholar] [CrossRef]
  10. Maymandi, H.; Perez Benitez, J.L.; Gallegos-Funes, F.; Perez Benitez, J.A. A novel monitor for practical brain-computer interface applications based on visual evoked potential. Brain-Comput. Interfaces 2021, 8, 1–13. [Google Scholar] [CrossRef]
  11. Kawala-Sterniuk, A.; Browarska, N.; Al-Bakri, A.; Pelc, M.; Zygarlicki, J.; Sidikova, M.; Martinek, R.; Gorzelanczyk, E.J. Summary of over Fifty Years with Brain-Computer Interfaces—A Review. Brain Sci. 2021, 11, 43. [Google Scholar] [CrossRef]
  12. Erbay, M.F.; Zayman, E.P.; Erbay, L.; Ünal, S. Evaluation of Transcranial Magnetic Stimulation Efficiency in Major Depressive Disorder Patients: A Magnetic Resonance Spectroscopy Study. Psychiatry Investig. 2019, 16, 745–750. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Wang, B.; Shen, M.R.; Deng, Z.-D.; Smith, J.E.; Tharayil, J.J.; Gurrey, C.J.; Gomez, L.J.; Peterchev, A.V. Redesigning existing transcranial magnetic stimulation coils to reduce energy: Application to low field magnetic stimulation. J. Neural Eng. 2018, 15, 036022. [Google Scholar] [CrossRef] [PubMed]
  14. Yunokuchi, K.; Koyoshi, R.; Wang, G.; Yoshino, T.; Tamari, Y.; Hosaka, H.; Saito, M. Estimation of focus of electric field in an inhomogenous medium exposed by pulsed magnetic field. In Proceedings of the IEEE 1st Joint BMES/EMBS Conference Serving Humanity and Advancing Technology, Atlanta, GA, USA, 13–16 October 1999. [Google Scholar]
  15. Wilson, S.A.; Thickbroom, G.W.; Mastaglia, F.L. Transcranial magnetic stimulation mapping of the motor cortex in normal subjects. The representation of two intrinsic hand muscles. J. Neurol. Sci. 1993, 118, 134–144. [Google Scholar] [CrossRef] [PubMed]
  16. Thielscher, A.; Kammer, T. Electric field properties of two commercial Figure-8 coils in TMS: Calculation of focality and efficiency. Clin. Neurophysiol. 2004, 115, 1697–1708. [Google Scholar] [CrossRef] [PubMed]
  17. Deng, Z.-D.; Lisanby, S.H.; Peterchev, A.V. Electric field depth-focality tradeoff in transcranial magnetic stimulation: Simulation comparison of 50 coil designs. Brain Stimul. 2013, 6, 1–13. [Google Scholar] [CrossRef] [Green Version]
  18. Cohen, L.G.; Roth, B.J.; Nilsson, J.; Dang, N.; Panizza, M.; Bandinelli, S.; Friauf, W.; Hallett, M. Effects of coil design on delivery of focal magnetic stimulation. Technical considerations. Electroencephalogr. Clin. Neurophysiol. 1990, 75, 350–362. [Google Scholar] [CrossRef]
  19. Toschi, N.; Welt, T.; Guerrisi, M.; Keck, M.E. Transcranial magnetic stimulation in heterogeneous brain tissue: Clinical impact on focality, reproducibility and true sham stimulation. J. Psychiatr. Res. 2009, 43, 255–264. [Google Scholar] [CrossRef]
  20. Mills, K.R.; Boniface, S.J.; Schubert, M. Magnetic brain stimulation with a double coil: The importance of coil orientation. Electroencephalogr. Clin. Neurophysiol. 1992, 85, 17–21. [Google Scholar] [CrossRef]
  21. Salinas, F.S.; Lancaster, J.L.; Fox, P.T. Detailed 3D models of the induced electric field of transcranial magnetic stimulation coils. Phys. Med. Biol. 2007, 52, 2879–2892. [Google Scholar] [CrossRef]
  22. Thielscher, A.; Opitz, A.; Windhoff, M. Impact of the gyral geometry on the electric field induced by transcranial magnetic stimulation. Neuroimage 2010, 54, 234–243. [Google Scholar] [CrossRef] [PubMed]
  23. Knecht, S.; Sommer, J.; Deppe, M.; Steinstrater, O. Scalp position and efficacy of transcranial magnetic stimulation. Clin. Neurophysiol. 2005, 116, 1988–1993. [Google Scholar] [CrossRef] [PubMed]
  24. Sparing, R.; Buelte, D.; Meister, I.G.; Paus, T.; Fink, G.R. Transcranial magnetic stimulation and the challenge of coil placement: A comparison of conventional and stereotaxic neuronavigational strategies. Hum. Brain Mapp. 2008, 29, 82–96. [Google Scholar] [CrossRef]
  25. Maccabee, P.J.; Nagarajan, S.S.; Amassian, V.E.; Dur, D.M.; Szabo, A.Z.; Ahad, A.B.; Cracco, R.Q.; Lai, K.S.; Eberle, L.P. Influence of pulse sequence, polarity and amplitude on magnetic stimulation of human and porcine peripheral nerve. J. Physiol. 1998, 513, 571–585. [Google Scholar] [CrossRef] [PubMed]
  26. Peterchev, A.V.; Jalinous, R.; Lisanby, S.H. A transcranial magnetic stimulator inducing near-rectangular pulses with controllable pulse width (cTMS). IEEE Trans Biomed. Eng. 2008, 55, 257–266. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Casula, E.P.; Rocchi, L.; Hannah, R.; Rothwell, J.C. Effects of pulse width, waveform and current direction in the cortex: A combined cTMS-EEG study. Brain Stimul. 2018, 11, 1063–1070. [Google Scholar] [CrossRef] [Green Version]
  28. Geddes, L.A. Optimal stimulus duration for extracranial cortical stimulation. Neurosurgery 1987, 20, 94–99. [Google Scholar] [CrossRef]
  29. Taylor, J.L.; Loo, C.K. Stimulus waveform influences the efficacy of repetitive transcranial magnetic stimulation. J. Affect. Disord. 2007, 97, 271–276. [Google Scholar] [CrossRef]
  30. Foster, K.R.; Schwan, H.P. Dielectric properties of tissues. In Biological Effects of Electromagnetic Fields; CRC Press: Boca Raton, FL, USA, 1996; pp. 25–102. [Google Scholar]
  31. Geddes, L.A.; Baker, L.E. The specific resistance of biological material: A compendium of data for the biomedical engineer and physiologist. Med. Biol. Eng. 1967, 5, 271–293. [Google Scholar] [CrossRef]
  32. Pethig, R.; Kell, D.B. The passive electrical properties of biological systems: Their significance in physiology, biophysics, and biotechnology. Phys. Med. Biol. 1987, 32, 933–970. [Google Scholar] [CrossRef] [Green Version]
  33. Martinsen, O.G.; Grimmes, S.; Schwan, H.P. Interface phenomena and dielectric properties of biological tissue. In Encyclopedia of Surface and Colloid Science; Marcel Dekker: New York, NY, USA, 2002. [Google Scholar]
  34. Gabriel, S.; Lau, R.; Gabriel, C. The dielectric properties of biological tissues: II. Measurements in the frequency range 10 Hz to 20 GHz. Phys. Med. Biol. 1996, 41, 2251–2269. [Google Scholar] [CrossRef] [Green Version]
  35. Roth, B.J. Mechanisms for electrical stimulation of excitable tissue. Crit. Rev. Biomed. Eng. 1994, 22, 253–305. [Google Scholar] [PubMed]
  36. Rossi, S.; Hallett, M.; Rossini, P.M.; Pascual-Leone, A. The Safety of TMS Consensus Group, Safety, ethical considerations, and application guidelines for the use of transcranial magnetic stimulation in clinical practice and research. Clin. Neurophysiol. 2009, 120, 2008–2039. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Peterchev, A.V.; Wagner, T.A.; Miranda, P.C.; Nitsche, M.A.; Paulus, W.; Lisanby, S.H.; Pascual-Leone, A.; Bikson, M. Fundamentals of transcranial electric and magnetic stimulation dose: Definition, selection, and reporting practices. Brain Stimul. 2012, 5, 435–453. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  38. Goetz, S.M.; Deng, Z.D. The development and modelling of devices and paradigms for transcranial magnetic stimulation. Int. Rev. Psychiatry 2017, 29, 115–145. [Google Scholar] [CrossRef] [Green Version]
  39. Silva, S.; Basser, P.; Miranda, P. Elucidating the mechanisms and loci of neuronal excitation by transcranial magnetic stimulation using a finite element model of a cortical sulcus. Clin. Neurophysiol. 2008, 119, 2405–2413. [Google Scholar] [CrossRef] [Green Version]
  40. Ruohonen, J.; Ilmoniemi, R.J. Physical principles for transcranial magnetic stimulation. In Handbook of Transcranial Magnetic Stimulation; Arnold: London, UK, 2002; pp. 18–30. [Google Scholar]
  41. Heller, L.; van Hulsteyn, D.B. Brain stimulation using electromagnetic sources: Theoretical aspects. Biophys. J. 1992, 63, 129–138. [Google Scholar] [CrossRef] [Green Version]
  42. Roth, B.J.; Basser, P.J. A Model of the Stimulation of a nerve fiber by Electromagnetic Induction. IEEE Trans. Biomed. Eng. 1990, 37, 588–597. [Google Scholar] [CrossRef]
  43. Roth, B.J.; Saypol, J.M.; Hallett, M.; Cohen, L.G. A theoretical calculation of the electric field induced in the cortex during magnetic stimulation. Muscle Nerve 1991, 81, 47–56. [Google Scholar] [CrossRef]
  44. Abdeen, M.A.; Stuchly, M.A. Modeling of magnetic field stimulation of bent neurons. IEEE Trans. Biomed. Eng. 1994, 41, 1092–1095. [Google Scholar] [CrossRef]
  45. Jackson, J.D. Classical Electrodynamics; John Wiley & Sons: Hoboken, NJ, USA, 2012. [Google Scholar]
  46. Plonus, M.A. Applied Electromagnetics; John Wiley & Sons: Hoboken, NJ, USA, 2012. [Google Scholar]
  47. Davey, K.; Cheng, C.H.; Epstein, C.M. Prediction of magnetically induced electric fields in biological tissue. IEEE Trans. Biomed. Eng. 1991, 38, 418422. [Google Scholar] [CrossRef] [PubMed]
  48. Eaton, H. Electric field induced in a spherical volume conductor from arbitrary coils: Application to magnetic stimulation and MEG. Med. Biol. Eng. Comput. 1992, 30, 433–440. [Google Scholar] [CrossRef] [PubMed]
  49. Esselle, K.P.; Stuchly, M.A. Neural Stimulation with Magnetic Fields: Analysis of Induced Electric Fields. IEEE Trans. Biomed. Eng. 1992, 39, 693–700. [Google Scholar] [CrossRef] [PubMed]
  50. Grandori, F.; Ravazzani, P. Magnetic-stimulation of the motor cortex-Theoretical considerations. IEEE Trans. Biomed. Eng. 1991, 38, 180–191. [Google Scholar] [CrossRef]
  51. Durand, D.; Ferguson, A.S.; Dalbasti, T. Effects of surface boundary on neuronal magnetic stimulation. IEEE Trans. Biomed. Eng. 1992, 37, 588–597. [Google Scholar] [CrossRef]
  52. Kim, D.H.; Sykulski, J.K.; Loukaides, N.; Georghiou, G.E. Numerical investigation of the electric field distribution induced in the brain by transcranial magnetic stimulation. IEE Proc. Sci. Meas. Technol. 2004, 151, 479–483. [Google Scholar] [CrossRef] [Green Version]
  53. Kowalski, T.; Silny, J.; Buchner, H. Current density threshold for the stimulation of neurons in the motor cortex area. Bioelectromagnetics 2002, 23, 421–428. [Google Scholar] [CrossRef]
  54. Krasteva, V.T.; Papazov, S.P.; Daskalov, I.K. Magnetic stimulation for non-homogeneous biological structures. Biomed. Eng. 2002, 1, 3. [Google Scholar] [CrossRef] [Green Version]
  55. Miranda, P.C.; Hallett, M.; Basser, P.J. The electric field induced in the brain by magnetic stimulation: A 3-D finite-element analysis of the effect of tissue heterogeneity and anisotropy. IEEE Trans. Biomed. Eng. 2003, 50, 1074–1085. [Google Scholar] [CrossRef]
  56. Bungert, A.; Antunes, A.; Espenhahn, S.; Thielscher, A. Where does TMS stimulate the motor cortex? Combining electrophysiological measurements and realistic field estimates to reveal the affected cortex position. Cereb. Cortex 2017, 27, 5083–5094. [Google Scholar] [CrossRef] [Green Version]
  57. Gomez-Tames, J.; Laakso, I.; Hirata, A. Review on biophysical modelling and simulation studies for transcranial magnetic stimulation. Phys. Med. Biol. 2020, 65, 24TR03. [Google Scholar] [CrossRef] [PubMed]
  58. Petrov, P.I.; Mandija, S.; Sommer, I.E.; Van Den Berg, C.A.; Neggers, S.F. How much detail is needed in modeling a transcranial magnetic stimulation figure-8 coil: Measurements and brain simulations. PLoS ONE 2017, 12, e0178952. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  59. Li, H.; Deng, Z.-D.; Oathes, D.; Fan, Y. Computation of transcranial magnetic stimulation electric fields using self-supervised deep learning. NeuroImage 2022, 264, 119705. [Google Scholar] [CrossRef] [PubMed]
  60. Taha, H. Implementation of Darwin’s Model by the Finite Element Method in Order to Model Electrical Machines at Intermediate Frequencies. Ph.D. Thesis, Université de Lille, Lille, France, 2021. [Google Scholar]
  61. Larsson, J. Electromagnetics from a quasistatic perspective. Am. J. Phys. 2006, 75, 230–239. [Google Scholar] [CrossRef] [Green Version]
  62. Salvador, R.N.B.F. Numerical Modelling in Transcranial Magnetic Stimulation. Ph.D. Thesis, Universidade de Lisboa, Lisboa, Portugal, 2009. [Google Scholar]
  63. LeVeque, R.J. Finite Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2007. [Google Scholar]
  64. Gholami, A.; Malhotra, D.; Sundar, H.; Biros, G. FFT, FMM, or multigrid? A comparative study of state-of-the-art Poisson solvers for uniform and nonuniform grids in the unit cube. SIAM J. Sci. Comput. 2016, 38, C280–C306. [Google Scholar] [CrossRef] [Green Version]
  65. Toschi, N.; Welt, T.; Guerrisi, M.; Keck, M.E. A reconstruction of the conductive phenomena elicited by transcranial magnetic stimulation in heterogeneous brain tissue. Phys. Medica 2008, 24, 80–86. [Google Scholar] [CrossRef]
  66. Cerri, G.; Leo, R.D.; Moglie, F.; Schiavoni, A. An accurate 3-D model for magnetic stimulation of the brain cortex. J. Med. Eng. Technol. 1995, 19, 7–16. [Google Scholar] [CrossRef]
  67. Woolfson, M.M.; Pert, G.J. An introduction to Computer Simulation; Oxford University Press: Oxford, UK, 1999. [Google Scholar]
  68. Polycarpou, A.C. Introduction to the Finite Element Method in Electromagnetics; Synthesis Lectures on Computational Electromagnetics Series; Morgan & Claypool Publishers: San Rafael, CA, USA, 2022. [Google Scholar]
  69. Jin, J.-M. The Finite Element Method in Electromagnetics, 3rd ed.; John Wiley & Sons: Hoboken, NJ, USA, 2014. [Google Scholar]
  70. Webb, J.P. Application of the finite-element method to electromagnetic and electrical topics. Rep. Prog. Phys. 1995, 58, 1673–1712. [Google Scholar] [CrossRef]
  71. Ueno, S.; Tashiro, T.; Harada, K. Localized stimulation of neural tissues in the brain by means of a paired configuration of time-varying magnetic fields. J. Appl. Phys. 1988, 64, 5862–5864. [Google Scholar] [CrossRef] [Green Version]
  72. Branston, N.M.; Tofts, P.S. Analysis of the distribution of currents induced by a changing magnetic field in a volume conductor. Phys. Med. Biol. 1991, 36, 161–168. [Google Scholar] [CrossRef]
  73. Gomez, L.J.; Yucel, A.C.; Hernandez-Garcia, L.; Taylor, S.F.; Michielssen, E. Uncertainty quantification in transcranial magnetic stimulation via high dimensional model representation. IEEE Trans Biomed. Eng. 2015, 62, 361–372. [Google Scholar] [CrossRef] [PubMed]
  74. Güllmar, D.; Haueisen, J.; Reichenbach, J.R. Influence of anisotropic electrical conductivity in white matter tissue on the EEG/MEG forward and inverse solution.A high-resolution whole head simulation study. Neuroimage 2010, 51, 145–163. [Google Scholar] [CrossRef] [PubMed]
  75. Opitz, A.; Windhoff, M.; Heidemann, R.M.; Turner, R.; Thielscher, A. How the brain tissue shapes the electric field induced by transcranial magnetic stimulation. NeuroImage 2011, 58, 849–859. [Google Scholar] [CrossRef] [PubMed]
  76. Wu, Y.X.; Yu, H.Y.; Liu, Z.W. Numerical Investigation of the Magnetic and Electric Field Distributions Produced by Biconical Transcranial Magnetic Stimulation Coil for Optimal Design. IEEE Trans. Magn. 2018, 54, 11. [Google Scholar] [CrossRef]
  77. Kamitani, Y.; Bhalodia, V.M.; Kubota, Y.; Shimojo, S. A model of magnetic stimulation of neocortical neurons. Neurocomputing 2001, 38, 697–703. [Google Scholar] [CrossRef] [Green Version]
  78. Seo, H.; Jun, S.C. Relation between the electric field and activation of cortical neurons in transcranial electrical stimulation. Brain Stimul. 2019, 12, 275–289. [Google Scholar] [CrossRef]
  79. Romero, M.C.; Davare, M.; Armendariz, M.; Janssen, P. Neural effects of transcranial magnetic stimulation at the single-cell level. Nat. Commun. 2019, 10, 2642. [Google Scholar] [CrossRef] [Green Version]
  80. Laakso, I.; Murakami, T.; Hirata, A.; Ugawa, Y. Where and what TMS activates: Experiments and modeling. Brain Stimul. 2018, 11, 166–174. [Google Scholar] [CrossRef]
  81. Aberra, A.S.; Wang, B.; Grill, W.M.; Peterchev, A.V. Simulation of transcranial magnetic stimulation in head model with morphologically-realistic cortical neurons. Brain Stimul. 2020, 13, 175–189. [Google Scholar] [CrossRef] [Green Version]
  82. Htet, A.T.; Saturnino, G.B.; Burnham, E.H.; Noetscher, G.; Nummenmaa, A.; Makarov, S.N. Comparative performance of the finite element method and the boundary element fast multipole method for problems mimicking transcranial magnetic stimulation (TMS). J. Neural. Eng. 2019, 16, 024001. [Google Scholar] [CrossRef] [Green Version]
  83. Wagner, T.A.; Zahn, M.; Grodzinsky, A.J.; Pascual-Leone, A.P. Three-dimensional head model simulation of transcranial magnetic stimulation. IEEE Trans. Biomed. Eng. 2004, 51, 1586–1598. [Google Scholar] [CrossRef] [PubMed]
  84. Tofts, P.S. The distribution of induced currents in magnetic stimulation of the nervous system. Phys. Med. Biol. 1990, 35, 1119–1128. [Google Scholar] [CrossRef]
  85. Saturnino, G.B.; Madsen, K.H.; Thielscher, A. Electric field simulations for transcranial brain stimulation using FEM: An efficient implementation and error analysis. J. Neural. Eng. 2019, 16, 066032. [Google Scholar] [CrossRef] [PubMed]
  86. Chen, M.; Mogul, D.J. A structurally detailed finite element human head model for simulation of transcranial magnetic stimulation. J. Neurosci. Meth. 2009, 179, 111–120. [Google Scholar] [CrossRef] [PubMed]
  87. Granula, B.; Porzig, K.; Toepfer, H.; Gacanovic, M. A Comparison between an A-V and V Formulation in Transcranial Magnetic Stimulation. In Proceedings of the 2013 COMSOL Conference, Rotterdam, The Netherlands, 23–25 October 2013. [Google Scholar]
  88. Roth, B.J.; Cohen, L.G.; Hallett, M.; Friauf, W.; Basser, A.P.J. A Theoretical Calculation of the Electric Field Induced by Magnetic Stimulation of a Peripheral Nerve. Muscle Nerve 1990, 13, 734–741. [Google Scholar] [CrossRef] [PubMed]
  89. Janssen, A.M.; Rampersad, S.M.; Lucka, F.; Lanfer, B.; Lew, S.; Aydin; Wolters, C.H.; Stegeman, D.F.; Oostendorp, T.F. The influence of sulcus width on simulated electric fields induced by transcranial magnetic stimulation. Phys. Med. Biol. 2013, 58, 4881–4896. [Google Scholar] [CrossRef] [Green Version]
  90. Windhoff, M.; Opitz, A.; Thielscher, A. Electric Field Calculations in Brain Stimulation Based on Finite Elements: An Optimized Processing Pipeline for the Generation and Usage of Accurate Individual Head Models. Hum. Brain Mapp. 2013, 34, 923–935. [Google Scholar] [CrossRef]
  91. Laakso, I.; Hirata, A. Fast multigrid-based computation of the induced electric field for transcranial magnetic stimulation. Phys. Med. Biol. 2012, 57, 7753–7765. [Google Scholar] [CrossRef]
  92. Yamamoto, K.; Takiyama, Y.; Saitoh, Y.; Sekino, M. Numerical Analyses of Transcranial Magnetic Stimulation Based on Individual Brain Models by Using a Scalar-Potential Finite-Difference Method. IEEE Trans. Magn. 2016, 52, 7. [Google Scholar] [CrossRef]
  93. Gomez, L.J.; Dannhauer, M.; Koponen, L.M.; Peterchev, A.V. Conditions for numerically accurate TMS electric field simulation. Brain Stimul. 2020, 13, 157–166. [Google Scholar] [CrossRef] [Green Version]
  94. Thielscher, A.; Antunes, A.; Saturnino, G.B. Field modeling for transcranial magnetic stimulation: A useful tool to understand the physiological effects of TMS? In Proceedings of the 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), Milan, Italy, 25–29 August 2015; pp. 222–225. [Google Scholar]
  95. Wang, W.; Eisenberg, S.R. A Three-Dimensional Finite Element Method for Computing Magnetically Induced Currents in Tissues. IEEE Trans. Magn. 1994, 30, 6. [Google Scholar]
  96. Salvador, R.; Mekonnen, A.; Ruffini, G.; Miranda, P.C. Modeling the electric field induced in a high resolution realistic head model during transcranial current stimulation. In Proceedings of the 2010 Annual International Conference of the IEEE Engineering in Medicine and Biology, Buenos Aires, Argentina, 31 August–4 September 2010; pp. 2073–2076. [Google Scholar]
  97. Miranda, P.C.; Correia, L.; Salvador, R.; Basser, P.J. Tissue heterogeneity as a mechanism for localized neural stimulation by applied electric fields. Phys. Med. Biol. 2007, 52, 5603–5617. [Google Scholar] [CrossRef]
  98. Saypol, J.M.; Roth, B.J.; Cohen, L.G.; Hallett, M. A theoretical comparison of electric and magnetic stimulation of the brain. Ann. Biomed. Eng. 1991, 19, 317–328. [Google Scholar] [CrossRef]
  99. Nadeem, M.; Thorlin, T.; Gandhi, O.; Persson, M. Computation of electric and magnetic stimulation in human head using the 3-D impedance method. IEEE Trans. Biomed. Eng. 2003, 50, 900–907. [Google Scholar] [CrossRef]
  100. Daneshzanda, M.; Makarova, S.N.; de Lara, L.I.N.; Guerin, B.; McNab, J.; Rosen, B.R.; Hamalainen, M.S.; Raij, T.; Nummenmaa, A. Rapid computation of TMS-induced E-fields using a dipole-based magnetic stimulation profile approach. NeuroImage 2021, 237, 1180897. [Google Scholar] [CrossRef]
  101. Ai, L.; Bansal, P.; Mueller, J.K.; Legon, W. Efects of transcranial focused ultrasound on human primary motor cortex using 7T fMRI: A pilot study. BMC Neurosci. 2018, 19, 56. [Google Scholar] [CrossRef] [PubMed]
  102. Zhang, D.; Li, H.; Sun, J.; Hu, W.; Jin, W.; Li, S.; Tong, S. Antidepressant-Like Effect of Low-Intensity Transcranial Ultrasound Stimulation. IEEE Trans. Biomed. Eng. 2019, 66, 411–420. [Google Scholar] [CrossRef] [PubMed]
  103. Nagarajan, S.S.; Dur, D.M. A generalized cable equation for magnetic stimulation of axons. IEEE Trans. Biomed. Eng. 1996, 43, 304–312. [Google Scholar] [CrossRef] [PubMed]
  104. Wang, B.; Grill, W.M.; Peterchev, A.V. Coupling magnetically induced electric fields to neurons: Longitudinal and transverse activation. Biophys. J. 2018, 115, 95–107. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  105. Wu, T.; Fan, J.; Lee, K.S.; Li, X. Cortical neuron activation induced by electromagnetic stimulation: A quantitative analysis via modeling and simulation. J. Comput. Neurosci. 2016, 40, 51–64. [Google Scholar] [CrossRef]
Figure 1. Domains and domain boundaries for TMS boundary value problems.
Figure 1. Domains and domain boundaries for TMS boundary value problems.
Brainsci 13 01142 g001
Table 1. Domains, their boundaries, and electric and magnetic properties.
Table 1. Domains, their boundaries, and electric and magnetic properties.
DomainBoundaryApplicable PropertiesDescription
Ω Γ σ 0 , ϵ r 1 , μ r = 1 System domain that includes all other domains: Ω = Ω a Ω c Ω b . Its boundary is limited by the air exterior boundary: Γ Γ a = Γ
Ω a Γ a σ = 0 , ϵ r = 1 , μ r = 1 Air
Ω c Σ c (integration path) σ = { 0 , σ c u p p e r } , ϵ r = 1 , μ r = 1 Coil
Ω b Γ b σ 0 , ϵ r 1 , μ r = 1 The brain, composed of several brain tissues ( Ω b = i = 0 n Ω b i ); its boundary is limited by the air ( Γ a Γ b = Γ b )
Ω b i Γ b i σ 0 , ϵ r 1 , μ r = 1 Brain tissues. Some of these domains are considered to, partially share boundaries ( Γ b i j = 0 n Γ b j O , i j )
Table 2. Studies of TMS numerical simulation-based MQS- ϕ boundary value problems.
Table 2. Studies of TMS numerical simulation-based MQS- ϕ boundary value problems.
BVPModelPublications
MQS only E p Several coils and coil configurations in an air box[21,50]
A coil over a passive cable[42]
Figure-8 coil over a tissue planar interface[44]
A quasispherical volume conductor and a paired coil[71]
A neurocortical neuron model and a coil[77]
E p + skin effectA coil over an non-homogeneous volume conductor[54]
Laplace MQS- ϕ A coil over a cylindrical volume conductor[43]
Three types coils over a spherical volume conductor[48]
An arbitrarily shaped coil over a half-plane conductor[49]
A circular coil over a spherical conductor[72,88,97]
A coil over a half-plane tissue[84]
Poisson MQS- ϕ Figure-8 coil over a realistic brain model[22,56,58,75,79,80,81,82,85,87,90,91,92,94]
A circular coil over a parallelepiped volume conductor[47,51,95]
A coil over an approximate brain model[52,93]
Figure-8 coil over several brain models[73]
Figure-8 coil over a high-resolution brain model[74,86]
Uniform and realistic E-fields and a realistic brain model[78]
MQS A- ϕ An 8-shaped coil over a cortical sulcus[39,98]
A circular coil over a realist head model[53]
A custom coil over three concentric spheres[76]
Figure-8 coil over an approximate head model[83]
MQS full Maxwell equationsFigure-8 coil over a brain approximate model[62]
Table 3. TMS numerical simulation studies based on different BVPs.
Table 3. TMS numerical simulation studies based on different BVPs.
Boundary Value ProblemSimulation ToolTypePublications
MQS only EpMatlabDirect implementation[21]
Equations solved with VAX 750 and FortranDirect implementation[42]
Laplace MQS- ϕ FEM using version 4.7 of SCIRun: A Scientific Computing Problem Solving Environment (SCI), Utah, USACustom software[58]
FEM model using Matlab and the GetFEM++ libraryDirect implementation[75,89]
Comsol multiphysicsCommercial general purpose software[96]
Poisson MQS- ϕ FEM custom-written Matlab and C++ routines together with Getfem++ functionsDirect implementation[22]
FEM softwareUnknown[52]
Ansoft finite-element packageCommercial general-purpose software[97]
SimNIBS pipelineOpen-source software for TMS simulation[56,94]
SimBio software environmentOpen-source software for TMS simulation[74,81,82]
SimNIBSOpen-source software for TMS simulation[78,79,85,93]
Matlab and C++Direct implementation[80,90,91]
Matlab FDMDirect implementation[92]
I-DEAS FEM packageCommercial general-purpose software[95]
MQS A- ϕ FEM implemented by Comsol MultiphysicsCommercial general-purpose software[39]
Eddy current solver from the commercial FEM package Maxwell3D from AnsoftCommercial software[53,83,86]
Full Maxwell equationsFEM implemented by Comsol MultiphysicsCommercial general-purpose software[62]
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Pérez-Benítez, J.A.; Martínez-Ortiz, P.; Aguila-Muñoz, J. A Review of Formulations, Boundary Value Problems and Solutions for Numerical Computation of Transcranial Magnetic Stimulation Fields. Brain Sci. 2023, 13, 1142. https://doi.org/10.3390/brainsci13081142

AMA Style

Pérez-Benítez JA, Martínez-Ortiz P, Aguila-Muñoz J. A Review of Formulations, Boundary Value Problems and Solutions for Numerical Computation of Transcranial Magnetic Stimulation Fields. Brain Sciences. 2023; 13(8):1142. https://doi.org/10.3390/brainsci13081142

Chicago/Turabian Style

Pérez-Benítez, J. A., P. Martínez-Ortiz, and J. Aguila-Muñoz. 2023. "A Review of Formulations, Boundary Value Problems and Solutions for Numerical Computation of Transcranial Magnetic Stimulation Fields" Brain Sciences 13, no. 8: 1142. https://doi.org/10.3390/brainsci13081142

APA Style

Pérez-Benítez, J. A., Martínez-Ortiz, P., & Aguila-Muñoz, J. (2023). A Review of Formulations, Boundary Value Problems and Solutions for Numerical Computation of Transcranial Magnetic Stimulation Fields. Brain Sciences, 13(8), 1142. https://doi.org/10.3390/brainsci13081142

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop