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

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.

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.

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].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.

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: The constitutive equations should also be considered: as well as Ohm's law 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: 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: 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: Substituting Equations ( 3) and ( 6) and the constitutive Equation ( 2) into the Maxwell-Ampere Equation (1b) results in: 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: 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.

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: 1.
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; 2.
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; 3.
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 (σ whitematter = 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: where J ind = σ E = −σ ∂A ∂t + ∇φ are the induced conductive currents or eddy currents, and 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: 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 ind and J dis from Equation ( 9).This approximation does not mean that the induced current is zero but that J ind << J s .Furthermore, all the materials involved in the phenomenon are paramagnetic (µ r = 1).Therefore, Equation (9) becomes: 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: On the other hand, in the domain of the brain, neglecting the displacement current, J dis = 0, and since there is no external current, J s = 0, and Equation (10) becomes: which could also be expressed as: 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: is called the primary electric field and 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.

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.
According Ohm's Law, J p i = σ i E p , and considering that electric conductivities of the interfacing areas are different, (σ 1 = σ 2 ), at the interfaces: Therefore, 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]: 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 i s = σ i E i s ).In this case, Equation ( 21) can be expressed as: or 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, σ air = σ 2 = 0.Then, Equation (23) reduces to: or 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 (26) 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).

MQS A-φ Formulation
If only capacitive effects are neglected ( J dis = 0) using the general transformations and considerations described in Section 3.2, Equations ( 7) and ( 8) can be expressed as: or using the Coulomb gauge: 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): 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.

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: 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 irr + E sol , where ∇ × E irr = 0 and ∇ • E sol = 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: 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.

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.Table 1.Domains, their boundaries, and electric and magnetic properties.

Domain Boundary
Applicable Properties Description System domain that includes all other domains: Its boundary is limited by the air exterior boundary: The brain, composed of several brain tissues Brain tissues.Some of these domains are considered to, partially share boundaries Based on these assumptions related to the domains, the boundary value problems(BVP) can be expressed as follows: The Laplace MQS-φ BVP: The Poisson MQS-φ BVP: The MQS A-φ BVP: The Darwin model BVP: Several works have used the MQS φ formulation.Table 2 lists some of these works and their models: A neurocortical neuron model and a coil [77] E p + skin effect A 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]  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

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 (Ω coil ), which gives: 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: where Σ coil is the path.Equation ( 37) can be numerically computed as follows: where and 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: where The FDM to solve this Poisson equation can be implemented using the interactive method: where α = ,k E z i,j,k are the currents densities due to the primary electric field (E p = (E x p , E y p , E z p )), 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 ind = 0, This expression should be solved using the boundary conditions (32d), which can be expressed as: which is solved using the finite difference as: φ i,j,k = 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.

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: 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: 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: A e (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 φ e i , A e i = (A x ) e i i + (A y ) e i j + (A z ) e i 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 ) = ∑ n n j=1 a j i 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: 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 e N (Φ e ) = (Φ e ) T P e Φ e + (Φ e ) T Q e (51 where Φ e = (φ e 1 , φ e 2 , φ e 3 , φ e 4 , . . ., φ e n n ) 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: and 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: and where V e is the volume of the element, and Equation (51) corresponds to the power dissipated in one element.The total power is obtained as the sum of the power of all elements: 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 dim(Φ e ) = n n × N e , dim(Q e ) = n n × N e and dim(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: Similarly to the vector (Q e ), the matrix (P e ) should look like: In the new space, which is called disjoint space, Equation (56) can be written as: 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: Substituting this equation into Equation (69) gives: where P = CP dis C T and Q = CQ dis .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: which results in the following linear equation system: 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 (dim(P ) = (N − 1) × (N − 1))) and a vector (Q (dim(Q ) = 1 × (N − 1))).The solution for the equations system is: After computing the scalar potential, it is possible to obtain the electric field using Equation (6).

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 : 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.
using the identity ∇ • w A = w ∇ • A + A • ( ∇w), which implies that: Using this identity, Equation ( 67) can be modified as: Applying the divergence theorem [45] ( Ω e ∇ • F dV = Γ e F • n dS) to the first term of Equation ( 68) and rearranging the equation yields: 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 ( − According to the Garlerkin formulation [68], the weight function (w) should be set equal to the interpolation functions: Moreover, using Equation (47a) for the interpolation of the scalar potential (φ e ) and considering the interpolation for the induced current density values: Equation (70) becomes a system of n n equations, where n n is the number of element nodes: where 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) Ω sys .Accordingly, this term can be removed from the equation, and the equation system becomes: 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: 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: which gives the following linear equation system: 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): The value of Φ for the boundary nodes can be obtained from these equations for each element separately, forming a 3 × 3 linear equation system: where N s i represents the interpolation functions for the potential of the nodes forming the triangular face, that is, facing the outer boundary, and (E 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: and Therefore, 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: J dis = CJ asem (87) and substituting in (85): where P asem = P dis C and K asem = Q dis CJ asem , and the potential can be obtained from: 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: The procedure to obtain the electric field is described as follows: 1.
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: which, resolved with respect to (A t x ) i,j,k , gives: where α = ∆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); 3.
The same procedure is applied to compute A y and A z ; 4.
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.

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: or which, applying the identity a∇ 2 b = ∇ • (a∇b) − ∇a • ∇b to the third term, becomes: Using the divergence theorem as in the previous sections applied to the third term and rearranging the equation yields: The potential at any position ((x, y, z) Ω sys ) is computed using the potential at the nodes and the interpolation functions: Replacing the weight function, (w) with the interpolation functions in Equation ( 96) and A x , φ with expression (97)  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: where and

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 timeindependent, Equations (35a) and (35b) can be expressed as: 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: The term J dis = ∆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 ind ), which is magnetically induced.Therefore, these equations can be rearranged as: 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.

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: where P, T and Q are the same as in the Φ-V formulation solution, and , is constant for the nodes of each element.

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: The solution of the magnetic vector potential generated by a coil formed by N thin wires carrying a current (I) is: The primary electric field can be computed from the magnetic potential vector; therefore: The secondary electric field ( E s ) of a coil element (d l) is computed from the scalar potential as: Consequently, the total electric field of the coil element (d l) is given as: 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: Therefore, it is only necessary to compute dE x and dE 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: A general solution of the Laplace Equation (105b) can also be expressed as the Bessel integral: 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: This expression is replaced in Equation ( 109) to obtain components dE x and dE y , which are given by: and The differential of the total electric field is obtained as dE = dE x i + dE 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: where Y lm (θ, φ) 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: where The values of the coefficient ( C lm ) are obtained by dividing the domain (Ω b ) into rectangular components such that: 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 → ∞: where for l > 0 and F oo = 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].

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.

Overview of Some of the Main TMS Simulation Results
The effectiveness and applications of TMS depend on several parameters of the coilbrain 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].

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.

Figure 1 .
Figure 1.Domains and domain boundaries for TMS boundary value problems.
are the components of the primary electric field in the boundary space .
)with α lm = (l+m−1)(l+m) (2l+1)(2l−1) , β lm = (l−m−1)(l−m) (2l+1)(2l−1) , γ lm = (l−m)(l+m) (2l+1)(2l−1) , D lm = Electric scalar potential vector of element e Φ Electric scalar potential vector for all elements in the domain Γ J k i,j,k Nodal current density with the component k = {x, y, z} J e Current density vector of the element e J Current density vector of all the elements in the domain Γ A k i,j,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 Γ

Table 2 .
Studies of TMS numerical simulation-based MQS-φ boundary value problems.
Ω e N i dV and B e x = 1 µ o Ω e Ni ∂Ax ∂x n x + ∂Ax ∂y n y + ∂Ax ∂z n z dV is the boundary condition as in the previous section.
gives: P e (A t x ) e + T e (A t x ) e + R e Φ + T e (A t−∆t x ) e + F e + B e = 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 i

Table 3 .
TMS numerical simulation studies based on different BVPs.