Particle Swarm Optimization for Positioning the Coil of Transcranial Magnetic Stimulation

The distribution of the induced electric field (E-field) during transcranial magnetic stimulation (TMS) depends on the individual anatomical structure of the brain as well as coil positioning. Inappropriate stimulation may degrade the efficacy of TMS or even induce adverse effects. Therefore, optimizing the E-field according to individual anatomy and clinical need has become a research focus. In this paper, particle swarm optimization (PSO) was applied for the first time to the positioning of TMS coils with anatomical head models. We discuss the parameters of the PSO algorithm, which were optimized to achieve a reasonable convergence time suitable for in-time treatment planning. The optimizer improved the distribution of the induced E-field strength at the dedicated cortical region, with a mean value of 48.31% compared with that from the conventional treatment position. The optimization terminated after 4–11 iterations for 13 head models. The applicability and performance of the optimizer for a large population are discussed in terms of cortical complexity. This study could benefit not only clinics but also research on brain modulation.


Introduction
Transcranial magnetic stimulation (TMS) [1] is a widely used, noninvasive tool for modulating brain activity. Short pulses of magnetic fields are delivered to the cortex through current-carrying coils. e rapidly changing magnetic field induces current within the cortex and stimulates activity over a predefined cortical region [2]. Image-based navigation systems could provide accurate and repeatable coil positioning according to anatomical or contour-based landmarks. However, the distribution of the induced electric field (E-field) in the brain depends on individual anatomy [2,3], such as gyral orientation [4] and the dielectric properties of cortical tissues [2].
e maximum E-field strength is unnecessarily located beneath the focus of the coil, usually by several centimeters [5]. erefore, the coil's placement and orientation should be optimized to ensure appropriate stimulation of the target region [6][7][8][9].
e induced E-field in the target region can be maximized using extensive numerical simulations by traversing all possible configurations. Nevertheless, that process requires a huge time cost, making it unsuitable for in-time clinical planning.
Heuristic algorithms (e.g., genetic algorithm [10,11] and particle swarm optimization (PSO)) [12] have been reported to facilitate TMS coil design. Inspired by its efficiency at managing convergence rates and search diversity, we applied PSO [12] to the optimization of coil positioning and orientation. Tested using various anatomical head models, the proposed method could enhance local E-field strength by an average value of 48.31% (ranging between 10.98% and 116.57%) compared with the results of stimulation at conventional stimulation position. e number of simulations could be reduced to 3% of that used by the brute force method (traversing all possible configurations). Parametric optimization for PSO has been discussed. e relationship between optimization performance and cortical complexity has been investigated using the local gyrification index (lGI [13]) and local factual dimensionality (lFD [14]). e results have demonstrated the applicability of that method to a large population.
e proposed method facilitated individual treatment design with a reasonable time cost. Physicians can use this method to optimize the E-field distribution before the treatment almost automatically (the time cost can be around 3 to 10 hours as shown in our study by serial computation), without human intervention. Integrated with the navigation system, the proposed method will be beneficial to not only clinics but also brain studies. In addition, it aims to exploit the performance of an existing coil instead of designing a new one, which provides a low-cost option for researchers.

Positioning of the Head/Coil Models in the Simulations.
e purpose of optimization was to maximize the E-field strength in the region of interest (ROI) by selecting the position (x, y, z) and rotational angle (φ) of the coil center. e bottom-center of the model was positioned at the origin (0, 0, 0) of the Cartesian coordinate system, while the long axis of the head was aligned to the Z-axis. e coil moved around the scalp with a constant separation of 10 mm (to mimic the space occupied by the protective shell), and its plane was tangential to the surface. e rotational axis passed through the intersection point of the FOE coil and was perpendicular to the coil's surface. A rotational angle of 0°was defined as the coil's long axis being parallel to the X coordinate. Clinically, the stimulation was conventionally performed with the ROI beneath the intersection point of the FOE coil. Figures 2(a)-2(c) show the configurations.

Particle Swarm Optimization.
In PSO, a population of candidate particles is moved along the search surface, and measurements are made according to a given measure of quality (mathematical formula) that regulates the particle's solution (representing the coil's position and orientational angle in our study) and velocity [27]. Each particle's movement is influenced by its known position and that of the population in the search space. As such, the swarm is expected to move toward the best solutions.
PSO was initialized with a number of particles to search for optima. e solution of the particle at the i th iteration was where x i , y i , and z i represent the coordinates of the position of the coil center at the i th iteration, where φ i is the coil's orientational angle at the i th iteration. e particle updates its velocity and position according to where rand() is a random number between 0 and 1. c 1 and c 2 are learning factors (exploration and exploitation abilities), and usually, c 1 � c 2 � 2 to balance cognitive and social influences. ω is the inertial weight factor, which is between 0.9 and 1.2 [27]. e self-adaptive method is ω � ω max − t/ t max (ω max − ω min ) [28], with ω min and ω max being the minimum and maximum weights, respectively. t and t max are the current iterative number and the maximum iterative number, respectively. pos pi is the best solution that a particle has ever achieved, and pos pg is the best solution of the particle swarm. e coil was rotated along the axis (φ) with a step of 15°. e particles were initiated on a curved surface of 4 × 4 cm 2 spread over the cortex (with a spacing of 10 mm), and the ROI was located beneath its center. e optimization program was coded in Python.

Numerical Solver and
Hardware. e in-house ELF scalar potential finite difference (SPFD) solver was used to calculate the induced E-field distribution in the brains [29]. e solver used the incomplete lower-and uppermatrix preconditioner to speed up solution of the derived septa-diagonal matrix, where block Forward-Elimination and Backward-Substitution algorithms were developed to facilitate GPU-based multithread parallelization. is solver has been validated with commercial software and has been demonstrated to have high computational efficiency when processing ELF MF problems [29]. e solver's code is free to download at https://github.com/ licongsheng/OpenSPFD. e simulation volume was discretized into 1 × 1 × 1 mm 3 voxels. e hardware configuration was as follows: CPU: 2 × Xeon E5-2630, 2.2 GHz; memory on board: 256 GB; GPU: 2 × NVIDIA Tesla K40c with 24 GB total memory. e numerical simulation of each head under TMS with the given configuration took about 6 min.

Numerical Experiment.
In this study, stimulation of the motor cortex was used to evaluate the effects of optimization.
e motor cortex is involved in the planning, control, and execution of voluntary movements, and it is one of the most frequently stimulated regions in diagnostic and therapeutic applications [30]. e motor cortex can be divided into the primary motor cortex, premotor cortex, and supplementary motor area, with specific sites corresponding to various body movements [31,32]. Previous experimental and modeling-based reports have suggested that the size of a given muscle's representation is rather narrow [33]. In this numerical trial, a surface of 2 × 2 mm 2 in the middle of the precentral gyrus (Brodmann area 4, [34]) was selected as the ROI (Figure 3). e area belongs to the primary motor cortex. e initial position (IP) was set to C4 electrode point of international 10-20 electroencephalography system so that the ROI was directly beneath the intersection of the FOE coil when the initial rotational angle was 0°. is positioning method is frequently used in clinics [9].

Local Cortical Quantification.
As the induced E-field distribution depends on local anatomy, we investigated the effects of cortical geometry and complexity on the convergence of the optimization. Two metrics were used in the local cortical measurements: lGI and lFD.
lGI quantifies the amount of cortex buried within the sulcal folds as compared with the amount on the outer visible cortex. A cortex with extensive folding has a large gyrification index, whereas one with limited folding has a small gyrification index. lGI is the ratio of the total pial cortical surface over the perimeter of the brain delineated on 2-D coronal sections. Some neuroimaging processing tools, such as Freesurfer [35], provide pipelines for calculating lGI. However, Freesurfer's automatic pipeline uses the raw T1-weighted MRI as the input. e head models used in this study were segmented, and the raw data were unavailable. erefore, we developed tools to use the segmented models. First, the cortical surface was discretized on an adaptive triangular mesh using Amira ( ermo Fisher Scientific, Waltham, MA). Second, the local gyral vertex was connected on a smoothed surface and subsequently discretized on an adaptive triangular mesh. Accordingly, lGI can be obtained by the division of the two surface areas.
Cortical FD covaries with gyrification [36], and analyses of this relationship are useful for quantifying the convolutional properties of the cortex across multiple scales [37].
ere exist several algorithms to calculate the dimensionality measure [38]. We applied a dilation algorithm to measure the 3D structure. Using this method, each voxel of the 3D structure was replaced with a cube of given volume. e cube sizes can be dilated, usually by a multiple of 2, while the number of cubes is changed to fill the 3D volume. After taking the logarithm of both cube size and the count of cubes filled, the FD value was derived as in the following equation: e 3D cortical surface under the search surface of the particles was extracted for FD measurement. Calculation of 3D fractal dimensionality was done by a MATLAB program provided by Madan and Kensinger [39].
Both lGI and FD were calculated based on a search surface of 4 × 4 cm 2 .

Efficiency of PSO.
Using eight particles, we repeated the optimization of each head model five times. e difference between the five simulations was the initial particle positions, which were randomly generated. e E-field distribution with each optimization is shown in Figure 4. e detailed results are shown in Table 3 and summarized in Table 4.
e optimization converged at 4-11 iterations, theoretically corresponding to 192-528 min by serial computing (the actual time cost using the above-mentioned hardware was approximately 220-560 min). In comparison, the induced E-field strength can be enhanced by up to 116% (Table 3, Korean Adult, #2 optimization), with an overall improvement of about 43% for the 13 head models. e spatial deviation from the IP was up to 18 mm (Table 3). e maximum shift by optimization was 18 mm: the researchers alternately conducted the brute force simulations (simulations by traversing all possible configurations) on a surface of about 18 × 18 mm 2 to search for the optimal value, resulting in 3,888 candidate configurations (1 mm resolution and 12 rotational angles for each point). We conducted a validation experiment using a Chinese male head model as an example. e histogram of the calculated maximum E-field strength is shown in Figure 5. e calculated maximum E-field strength was 2.46 V/m, compared with 2.35 V/m by PSO with less than 70 simulations (8.4 iterations × 8 particles � 64). More than 98% of the simulations were saved, with a difference of only 4%.
Admittedly, the focality of the FOE coil was approximately a few cm 2 [40]. However, using our proposed method, we could manage the E-field distribution in a much finer region. ese notions are not contradictory because the present study aimed to direct the peak values to the predefined brain regions. In this approach, the operators could reduce the power delivered to the coil so that the above-threshold stimulation was achieved only in the ROI. As such, the coil's performance was exploited to achieve fit with a very narrow stimulation. e results listed in Tables 3 and 4 were averaged over ROIs of 2 × 2 mm 2 (i.e., four voxels). Hence, the statistical results for the 99 th percentile of E-field strength, which is prescribed by ICNIRP to reduce stair-case errors, are unavailable. e absolute percentage increase may be subject to change when other statistical metrics are applied. As E-field enhancement was found for all of the optimizations, the proposed optimization effectively improved the localized E-field strength. It should be note that the ROI size could potentially affect the performance of the algorithm method. e different ROIs may be used according to the clinical application so that the realistic performance improved need be investigated for various ROI definitions. e threshold for brain stimulation was above 100 V/m. In the study, we obtained the induced E-field strength at the level of several V/m. ere were mainly two reasons for it. Firstly, the coil used in the simulation was about half the diameter of a conventional FOE coil. e advantage was that it had better focality, but with reduced induced E-field strength. To validate our results, we can refer to existing literature using similar coils for rodent stimulation [41]. Secondly, the coil in the study had only one turn. In contrast, the clinical coils usually had 10 to 15 turns. As such, the induced E-field strength was further lowered. However, the purpose of the study was to present an optimization algorithm for the induced E-field strength. e results from this simplified FOE coil were still representative.

Optimization of PSO Parameters.
Previous studies concluded that the best approach to optimize PSO parameters is the rule of thumb, i.e., fixing the inertial weight while carefully selecting c 1 and c 2 [42] or vice versa. In general, parameter selection was empirical. We conducted numerical simulations to investigate parameter selection.
In this study, parametric optimization was initiated in terms of population size. Larger population size can accelerate convergence, but at the cost of an increased number of simulations per iteration. In addition, hardware parallelization should be taken into consideration when deciding the number of particles. Some studies have reported that PSO was not sensitive to population size [27] and that a population size of 20-30 is a conventional choice. We conducted trial simulations using 4, 6, 8, 12, 16, and 20 particles. With 8 particles, a 100% success rate was achieved for all head models. is particle count was selected because it was appropriate for the core size of the current CPU, which facilitated parallelization. e present study adopted the c 1 � c 2 � 2.05. It was proposed by the previous empirical studies on PSO [27]. Some studies have also proposed to fix the sum of c 1 and c 2 to 4.1 while adjusting the ratio of c 1 /c 2 from 2.8/1.3 to 1.3/2.8 [43]. We applied these coefficients to our simulations with a step of 0.2, using the Chinese adult female and male head models.
e results indicate that there is no significant difference. Prior knowledge could help to further reduce the number of simulations. For example, we may design the initial rotational angles of the FOE so that they are close to perpendicular to the local gyral orientation. It has been reported that this layout could lead to higher E-field strength [9].
Besides PSO, other methods based on brain Atlas [7] and deep neural networks [44] have also shown promise in facilitating accurate brain stimulation.

Relation between Local Anatomical Complexity and
Convergence Rate. Cortical quantification of the local cortical regions from various head models is shown in Table 5. As the induced E-field distribution was influenced by local anatomy, the convergence rate was assumed to change with the cortical complexity of the local cortex beneath the search surface. We conducted correlational analysis of the results from the 13 head models using Spearman's rank correlation coefficient and estimated the resulting statistical significance (α � 0.05). e results were as follows: r � 0.19 and p � 0.53 for lFD vs. mean E-field strength enhancement, r � 0.08 and p � 0.81 for lGI vs. mean E-field strength enhancement, r � − 0.07 and p � 0.82 for lGI vs. mean iteration, and r � 0.41, p � 0.16 for lGI vs. mean iteration. e results indicated that there was no significant correlation between local cortical complexity and either the convergence rate or the enhancement to the expected E-field strength. e findings need further investigation using more anatomical head models, which will be performed in our future work. e measured lGL and FD values fell in the average range for the population, according to anatomical/radiological reports [13,45]. In addition, the anatomical head models were generated with various segmentation tools and protocols. Accordingly, optimal performance could be expected from a large population.
In addition to anatomical factors, the most effective coil orientation depends on the shape of the induced current pulse. Further, when the first and second phases of the pulse are of similar size, it depends on the intensity of stimulation. Optimal mapping of the human motor cortex with magnetic stimulation requires knowledge of the influences on all these factors [2,46].  Figure 2: Positioning of the coil and the head model. During the stimulation, the distance between the coil and scalp was 10 mm, and its plane was tangential to the surface (a). e coil can move around the surface with a constant separation of 10 mm to the scalp (b). e rotational axis passed through the intersection point of the 8-shape coil and was perpendicular to the coil's surface (c). 0°of the rotational angle was defined when the long axis of the coil was parallel to X coordinate. e rotation angle φ was defined hereafter.

Conclusion
is study proposed numerical methods to optimize TMS coil positioning according to clinical needs. e application of PSO-based optimization to TMS coil positioning was first studied in this work. e versatility and efficiency of the optimizer have been numerically demonstrated. is study confirmed that the proposed algorithm is valid and efficient for providing optimal plans (in terms of induced E-field strength in the ROI) within a clinically acceptable period. A FOE coil was used in this study, and the proposed method further exploited its performance in brain stimulation. In conclusion, PSO can enhance the E-field strength in an ROI by a mean value of about 48%. e derived parameters can benefit robotic neuroimaging navigation systems by facilitating stable and desired cortical activation.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors report no conflicts of interest in this work. None of the authors have any personal or financial involvement with organizations that have a financial interest in its content.

Authors' Contributions
Congsheng Li and Chang Liu contributed equally to this work.