In-vivo verified anatomically aware deep learning for real-time electric field simulation

Objective. Transcranial magnetic stimulation (TMS) has emerged as a prominent non-invasive technique for modulating brain function and treating mental disorders. By generating a high-precision magnetically evoked electric field (E-field) using a TMS coil, it enables targeted stimulation of specific brain regions. However, current computational methods employed for E-field simulations necessitate extensive preprocessing and simulation time, limiting their fast applications in the determining the optimal coil placement. Approach. We present an attentional deep learning network to simulate E-fields. This network takes individual magnetic resonance images and coil configurations as inputs, firstly transforming the images into explicit brain tissues and subsequently generating the local E-field distribution near the target brain region. Main results. Relative to the previous deep-learning simulation method, the presented method reduced the mean relative error in simulated E-field strength of gray matter by 21.1%, and increased the correlation between regional E-field strengths and corresponding electrophysiological responses by 35.0% when applied into another dataset. In-vivo TMS experiments further revealed that the optimal coil placements derived from presented method exhibit comparable stimulation performance on motor evoked potentials to those obtained using computational methods. The simplified preprocessing and increased simulation efficiency result in a significant reduction in the overall time cost of traditional TMS coil placement optimization, from several hours to mere minutes. Significance. The precision and efficiency of presented simulation method hold promise for its application in determining individualized coil placements in clinical practice, paving the way for personalized TMS treatments.


Introduction
Transcranial magnetic stimulation (TMS) is a noninvasive brain stimulation technique widely applied in the modulation of brain functions such as motor and working memory, as well as in the treatment of psychiatric and neurological disorders, including depression, anxiety, and mild cognitive impairments [1][2][3][4][5][6].By generating magnetically induced electric fields (E-fields), neural activities within specific brain regions can be effectively modulated with a TMS coil [7].Nevertheless, considerable variability exists in the therapeutic outcomes of TMS treatments [8][9][10][11].To address this issue in part, many efforts have been pursued to pinpoint optimal TMS coil placements tailored to individual brain targets [12][13][14].Such optimization requires quantifying the distributions of the evoked E-fields for a massive number of coil placement choices [12,15].The implementation of an efficient E-field simulation methods would facilitate the development of personalized TMS interventions [8,9,[16][17][18].
The conventional approach to E-field simulation is to use numerical computations to approximate the solution of the partial differential equations (PDEs) within a subject-specific volume conductor model.Specifically, a realistic individualized head model is initially segmented and constructed from individual T1-or T2-weighted images to represent the conductivities of various anatomical tissues (figure 1(a)) within a finely triangulated mesh [19,20].Subsequently, the E-field distribution under one coil placement is simulated using a variety of numerical calculation methods, such as the finite element model (FEM) and the boundary element method [12,21,22].The precision and computation duration of these solutions are contingent upon the resolution of head meshes and the order of polynomial approximation [23].The considerable computational time (tens of seconds per coil placement) of E-field simulation poses a challenge when simulating a vast number of coil placements and pursuing high-precision numerical modelling.Recently accelerated solutions have been developed using a dipole-based magnetic stimulation profile approach [24,25] with a pre-computed magnetic stimulation profile or the auxiliary dipole method (ADM) that computes the mean of the E-field within targets [14].Besides, long pre-processing time is also needed for the creation of an accurate head model for these methods, indicating several preprocessing hours after acquiring a subject's magnetic resonance imagings (MRIs) [15,19].
Recent research indicates that deep learning methods have the potential to simulate the E-field statistically [26][27][28][29].Using individual MRI as input, deep neural networks (DNNs) can produce corresponding E-field distributions in real time, leveraging their efficient extraction of implicit tissue features and sidestepping the dimensional curse encountered in earlier machine learning methods [30,31].However, the implicit tissue features (as shown in figure 1(a)) extracted by DNNs might result in region-dependent performance when these networks are trained only on specific brain regions (e.g. the motor or Broca's areas) [26,27].Furthermore, directly mapping Efields from single-site MRIs can make DNN inferences more vulnerable to the acquisition conditions of the input MRI images [31,32].Specifically, in addition to brain tissue features, site-specific MRI features might also be extracted and utilized in the E-field simulation.This could undermine the model's generalizability to datasets for which sufficient re-training data is lacking.Therefore, additional constraints on model features and large region coverage on the brain cortex should be considered to ensure reliability during deep-learning-based E-field simulation.Moreover, previous studies most investigated the model performance (e.g.E-field simulation errors) through in-silico experiments.However, in-vivo validation of deep learning models for ascertaining TMS coil placement is also of significant importance.This verification may provide essential support for selecting effective simulation methods and facilitate future practical application on personalized medicine treatment.
Inspired by the generalization of computational simulation models, we designed a deep-learningbased method using explicit, supervised brain tissue representation as the intermediate variable (figure 1(a)), termed the anatomically aware, realtime e-field simulation (AnaRES) method.Given a coil placement, the presented method could generate evoked E-fields using individual T1-weighted imaging.The AnaRES method employs an attentional deep learning architecture and is trained with an extensive dataset (over 100TB) encompassing a majority of brain scalp regions.Given the large dataset of different coil configurations, trained model could be directly conducted on new subject in real time.To test the model precision, we implemented AnaRES-based coil placement optimization for more than 200 cortical areas from the brainnetome atlas (BNA) [32], Shen's atlas [33] and human connectome project (HCP)-styled parcellations [34].We validated the model using both in-silico experiments on E-field simulation and in-vivo TMS stimulation experiments on motor evoked potentials (MEPs).The results have demonstrated our proposed model generates highly accurate E-field simulation outcomes, comparable to computational methods but with considerably greater efficiency.

Methods
The deep-learning architecture of AnaRES (figure 2(c)) learns a mapping from individual MRI images to the E-field, with the explicit anatomical architecture as an intermediate variable.By optimizing E-fields with explicit anatomical distribution, the model could improve the accuracy of E-fields, acting as a good estimator to solve the governing PDE of TMS fields.The main diagram can be subdivided into local sampling (figure 2(a)) and anatomical-based Efield inference using an attentional U-net architecture (figure 2(c)).

Local sampling based on coil placements
In a TMS navigation system, the scalp position P (x, y, z) and orientation (roll, pitch, yaw) of a coil (considering the coil as a rigid body with six degrees of freedom) will influence the distribution of the evoked E-field within an individual brain.However, in a normal situation, the coil position P (x, y, z) is located on the patient's scalp and the axes of roll and pitch are always set to be tangential to the brain surface on the scalp point P (x, y, z), leaving the yaw orientation θ to be determined.
To further simplify the complexity of E-field regression, the first thing the presented method did was to unify the estimated scope of the E-field simulations under different coil placements (P, θ).A transformation based on coil placements (P, θ) was applied on original MRI coordinates into a new XYZ coordinate.In the transformation, the coil position was set as the origin of the new coordinate.The tangential plane of the coil position on scalp surface was set as the XY plane, with the coil direction as the positive orientation of the X-axis.The Y-axis is perpendicular to X-axis in the plane.The direction of the Z-axis was along the negative normal direction of the XY plane (shown in figure 2(d)).To reduce the simulation time, the range of the estimated scope is limited within the local region but includes the effective stimulation range.For the double coil, the sampling sizes were 80, 80, 40 mm for the X, Y, and Z axes (shown in figure 2(d)).After local sampling, the estimated scopes of the E-field simulations under various coil placements were unified to the same condition p 0 .

Attentional network architecture for the E-field simulation
The presented deep learning used two attentional U-Net architectures.The first attentional U-Net used sampled images from individual T1w MRIs, to infer the local distributions of five individual brain tissues (skin, bone, cerebrospinal fluid (CSF), gray matter (GM), and white matter (WM)) and null space.After a softmax layer, the output of the first attentional U-Net is supervised by the realistic head models adopted in the FEM computational methods.The second attentional U-Net is used to learns a mapping from different brain tissues to the E-fields.The output of the second attentional U-Net is the distribution of Efield strength.
The attentional U-Net architecture is a four-level encoder-decoder architecture (figure 2(b)), in which both the input and output have four dimensions (width, height, depth, and channel).Each level of the encoder layer is the combination of convolutional layers (3 × 3 × 3 kernel size), a batch normalization layer, a rectified linear unit (ReLU), and a maxpooling layer (2 × 2 × 2 or 2 × 2 × 1 kernel size), while each level of the decoder layer is a combination of convolutional layers (kernel size 3 × 3 × 3), a batch normalization layer, a ReLU, and an up-sampling layer (2 × 2 × 2 or 2 × 2 × 1 kernel size), as shown in figure 2(b).The number of channels in the high-level convolutional layers is double that of the lower-level layer, with a value of 8 for the lowest level.The attention layers [35] in the skip connection consist of two separate inputs from the same-level encoder and lowlevel decoder layers, and the output is the same-layer decoder layer.It uses the 1 × 1 × 1 convolutional layers to generate spatial attention weights on original input features based on the back-propagation process of prediction errors in the network training.Studies have showed the model with attentional layers learns to suppress irrelevant regions in an input image while highlighting salient features useful for a specific task, leading to an improvement on medical-image-based predictions [35].

The optimization of the deep-learning model
The optimization of the deep learning network is described below.Before network training, the estimation situations are unified into the same situation p 0 through local sampling technique.The local sampling technique aims to unify the scope of different coil placements to reduce the complexity of model training.In the network optimization, the domain (D) of coil placements is set in the range of coil angle (0 • < θ < 180 • ) as well as in the range of coil position (P) shown in figure 3(a).After that, the network could learn a mapping from the input T1w images (T) to E-fields ( Ẽ), with the explicit anatomical architecture ( Ã) as the intermediate variable to gain the generalizability and precision.Then E-field simulation ( Ẽ) is computed with the input of the anatomical structures ( Ã) and T1w images (T).
The entire deep-learning network can be described as a function f (W) to optimize the parameters W in equation (1).The first term in the objective for the network was to generate E-fields ( Ẽ) that are closed to the E-fields (E) computed by the high-precision computational method using the simulation of non-invasive brain stimulation (SimNIBS) software.The second term in the objective (equation ( 2)) is to regulate intermediate variables used to reference the E-filed distribution, by obtaining explicit brain tissues ( Ã) that could provide a basis distribution of different current conductivities.Explicit brain tissues were supervised by the results in realistic head model (A) which are reconstructed through headreco pipeline using T1w and T2w images [19].The parameter γ is the relative importance for two objectives.The network objective can be described as equation ( 2) minimize The anatomical tissues and E-fields objectives were optimized in a two-step way.First, the anatomical objective was achieved by optimizing the crossentropy loss (Loss anat ) between the predicted anatomical structures ( Ã) and anatomical distributions from realistic head model (A).Then the network achieved the regression objective through the optimization on the mean square error (Loss reg ) between the network simulations ( Ẽ) and the high-precision simulations (E).The whole loss function is: To optimize this network, the ground truth of the E-field was simulated using FEM simulations 486 000 times on 15 subjects.In the optimization, 10 subjects were used to train and validate the network parameters (80% training, 20% validation), and the other 5 subjects were used to test the network performance.The deep-learning networks were trained using 50 epochs with the NVIDIA GeForce RTX 3090 24G and the best one with the minimum loss was chosen.The exponential-decay learning rate was used with an initial value of 0.005.The Adam optimizer for the network parameters had the values of β 1 = 0.5, β 2 = 0.99.To increase the number of training samples, we adopted augmentation techniques by flipping all images (the T1w MRI, anatomical structures, and induced E-field images) on the X-axis, Y-axis, or both axes with an augmentation probability of 0.5 during the training process.

Implementation
To verify the performance of presented method, both in-silico and in-vivo experiments were implemented.First in the in-silico experiment, E-field simulation performance on different TMS coil placements targeting various brain regions was investigated by applying the trained deep-learning model across validation subjects.Both the simulation errors of Efields and optimized deviations of coil placement were evaluated.Then, to further confirm whether there is a practical influence on the stimulation performance, the in-vivo stimulation was designed and tested for the TMS stimulation of the right-hand first dorsal interosseous (FDI) muscles.The MEPs recorded from electromyography (EMG) were used as alternative proof for simulation performance.The datasets used for different experiments are shown in the supplementary figure S1.

MRI acquisition and preprocessing
In this study, 15 subjects (10 females, right-handed, aged 24-33 years) was selected from the HCP S1200 dataset with the selection criterion for the MRI scan parameters [36].Each subject has T1w, T2w, and diffusion MRIs acquired using Siemens 3T Skyra scanners.The T1 MRIs were scanned with TR = 2400 ms, TE = 2.14 ms, FOV = 224 × 224 mm with 0.7 mm isotropic resolution.The diffusion MRI data consisted of three shells of b = 1000, 2000, and 3000 s mm −2 with 90 diffusion weighted directions plus another 6 b = 0 acquisitions with 1.25 mm slice thickness.All datasets were first analyzed using the HCP minimal preprocessing pipeline [36,37].The structural MRI has been preprocessed for gradient distortion correction, co-registration, brain extraction, field map distortion correction, and bias field correction.The diffusion MRI data was primarily preprocessed by motion correction, eddy current distortion correction, and echo-planar images (EPI) susceptibility-induced field distortion correction.Further preprocessing details of the structural and functional data can be found in the HCP preprocessing pipeline [38,39] (https://github.com/Washington-University/HCPpipelines).
Besides, we also recruited 27 healthy subjects (14 males and 13 females, right-handed, aged 23-29 years) to test the MEPs on the FDI using different coil placements.All subjects provided informed consent to accept MRI acquisitions and motor stimulation experiments.The TMS stimulation experiments were approved by the ethical committee of the Institute of Automation, Chinese Academy of Sciences.Each subject was scanned to obtain T1w, T2w and diffusion MRIs using a GE 3 T Signa UHP scanner.Both T1w MRIs (TE = 4.66 ms, TR = 4666 ms, FOV = 25.2 × 25.2 mm) and T2w MRIs (TE = 114 ms, TR = 4002 ms, FOV = 25.2 × 25.2 mm) were scanned at a 0.7 mm isotropic resolution.Diffusion MRIs were scanned with a single shell (b = 1000, 60 directions) and another 6 non-diffusion acquisitions (b = 0) with 1.5 mm slice thickness.Corresponding non-diffusion images with reverse phase-encoding were also scanned for eddy current distortion correction and EPI susceptibility-induced field distortion correction.The density distribution of the T1w MRIs was normalized ('fslmaths') into the density distribution of the HCP training subjects.The diffusion images were preprocessed for eddy current distortion correction and EPI susceptibility-induced field distortion correction.

Brain targets and coil placement configurations 3.2.1. In-silico configurations
To better evaluate the optimization performance using different E-field simulation methods, cortical targets were selected from three parcellations based on different MRI modalities: (1) BNA atlas (based on anatomical connectivity), (2) Shen's atlas (based on functional connectivity) and (3) HCP-styled parcellation (based on multi-modality information) as optimization targets [32][33][34].For the BNA atlas and Shen's atlas, these volume-based atlases were first transformed into fsaverage_LR32k surface templates using Workbench software to pursue a high delineation of individual targets.All area labels on the surface templates were remapped onto individual surfaces by surface-based registration [40] and assigned the area label of the individual cortical voxels based on the nearest surface labels, aiming to remove individual region not belonging to the GM.The individual region was not considered if it was located in a cortical area with a greater than 35 mm scalp-tocortex distance from the area center because it exceeds the effective depth of the coil.
The generation scheme for the coil placements is provided in the Supplementary Method.The yaw orientation (θ) of a TMS coil was analyzed from 0 to 180 degrees in the posterior-anterior (PA) directions using a 15-degree interval.For the Magstim 70 mm double coil, the induced E-field strength was assumed to be the same as that obtained from an additional 180-degree rotation of the TMS coil in the scalp tangent plane.For each subject, there were 2790 coil positions (P) covering the majority of the left hemisphere, excluding the regions near the eyes and ears.

In-vivo configurations
In the in-vivo experiment, single-pulse TMS stimulations were applied by a Magstim Rapid2 stimulator (the Magstim Co. Ltd) with a Magstim 70 mm double coil.The stimulation frequency was restricted to a maximum of 0.2 Hz.Different coil placements were controlled and navigated by an industrial robotic system (Robot: UR5e cobot, Universal Robots, Teradyne Inc.; Camera: MicronTracker H3-60, Claron Technology Inc).The navigation system ensured that the figure-8 coil was placed on the tangential plane of a scalp position.Through the face keypoint registration [41], the registration errors of navigation system were no more than 3 mm distance and no more than 2 degrees on the tangential plane of the brain scalps, otherwise such registration will be repeated manually.The EMG electrodes were placed on the FDI muscle.EMG recordings of the FDI muscles were detected with an integrated dry electrode (SX230, Biometrics Ltd) and recorded by a data acquisition system (DataLINK DLK900, Biometrics Ltd).The EMG signals were amplified (1000x), band pass filtered (20-460 Hz), and digitally sampled at 5000 Hz.The peak-to-peak amplitudes of the MEPs were recorded after 10-50 ms of the TMS pulses.We conducted three successive sessions on 27 participants from two sites (17 participants did not receive the second session due to the long stimulation time).The details of the tasks for the three sessions can be seen in the section entitled 'in-vivo TMS stimulation validations for the motor region' .
For the brain target corresponding to right-hand FDI muscle, previous studies have investigated the relationship between the electrophysiological signals and average E-field strength of the motor region [42,43].The variations in the FDI MEPs of the right hand have been reported to be highly consistent with the variations in the simulated E-field strengths around the left gyral crown using computational methods [42].Based on the conclusions of these studies, we chose the corresponding brain target in the standard image as the region around the gyral crown, which was centered on the MNI [−36.0,−17.3, 65.8] coordinates within a 7 mm radius covering the left pre and postcentral gyrus (shown in supplementary figure S2).The individual targets were obtained based on non-linear registration from the standard template.The regions that did not belong to the GM of a given subject were removed from the brain targets.

E-fields simulation from finite element method
To train and test the generalizability of the deep learning methods, we simulated the high-resolution Efield on 15 HCP subjects with scalp locations across the left hemisphere using the SimNIBS software V3 [20,44,45] (see supplementary figure S3).The head model of each subject was first constructed by the headreco pipeline from the T1 and T2 weighted images of each subject, including 5 tissue types, WM, GM, CSF, skull, and skin [19].To avoid the automated pipeline to highly influence the tissue segmentation of subject-specific head model, we have made the strict manual inspection for the brain segmentation results.The segmented results are double-checked by two brain researchers with background in brain imaging and brain anatomy.Then, the TMS-induced Efield distribution of a Magstim 70 mm double coil was simulated throughout the whole brain by 16 FEM models [46] under different coil poses, including 2790 scalp positions for the left hemisphere and 12 yaw angles for each position.To lessen the influence of hair, the coil scalp positions were moved up 6 mm along its normal direction.The conductivities of the brain tissues were set to σ skin = 0.465S/M, σ skull = 0.010S/M, σ CSF = 1.654S/M, σ GM = 0.275S/M and σ WM = 0.126S/M [44,47].The GM and WM were assigned anisotropic conductivity values derived from diffusion tensor imaging using the volume normalized mapping approach [48].The induced E-fields were calculated assuming a quasi-static regime.The change rate of the coil current was set to 1 A/µs.The computational simulations provided precise Efield results that can be regarded as the ground truth and thus used to train and validate the precision of the deep learning networks.We provide an example for the FEM E-field distribution in the same coil position with varying coil directions in supplementary figure S4.There were 2790 × 12 × 15 subjects = 502 200 computational E-field simulations using over 3000 × 4.5 h per subject with 2.5 GHz 2000 Hygon C86 cores (6 GB RAM per core) and over 100 TB of storage resources.

Quantitative indicators to evaluate simulated E-fields
We quantified three indicators to evaluate the overall performance of the predicted distribution of Efield strength ( Ẽ) for each scalp position (P mn , m = 50, . . ., 80; n = 11, . . ., 100) across the following test measures: mean absolute error (MAE), mean relative errors (MREs) for different anatomical structures, and anatomy-related absolute error (AAE).Specifically, the anatomy-related error was defined as the MAE that occurs on the edge of the brain anatomical structures (D edge ) across subjects, as shown in equation (6).This measure described the accuracy of the simulated E-field strength at the anatomical boundary.We calculated not only the general performance but also the spatial distribution and their differences on the three quantitative indicators Ẽ (P mn ) − E (P mn ) Ẽi (P mn ) − E i (P mn ) . (6)

Determination and evaluation of optimal coil placement 3.5.1. Traditional optimization of TMS coil placement
We use traditional coil placement optimization to locate the optimal coil placement by maximizing the Efield strength within targets.This optimization consists of three steps to optimize the coil placement for a brain target: (1) a nearby search for potential coil placements (D p , D θ ) of a target; (2) the E-field simulations under different coil placements; (3) the qualification of E-field strength E ROI within target and the selection of the optimum P opt , θ opt based on the E-field strength criterion.Specifically, the domain of potential coil positions D p was a 20 mm circle of scalp positions centered on the scalp position nearest to a target area.The potential coil orientations D θ was chosen from 12 orientations ranging from 0 to 180 degrees.The optimal coil position P opt and orientation θ opt were selected from the coil placement that provided the maximum E-field strength within a target (equation ( 7)).
In the second step of traditional coil optimization, three E-field stimulation methods have been implemented.The computational results with Efield simulation using direct FEM method were first investigated, named as FEM-CPO.The optimal results with two deep-learning E-field simulation methods are then investigated, including the Yokota et al [26] method and the presented method, named as RawRES-CPO and AnaRES-CPO, respectively.

Deviation quantifications from optimal coil placement
To investigate deviations of deep-learning coil optimization results from the computational results, we quantified such deviations in three aspects: the optimal coil positions, the optimal coil orientations, and the optimal E-field strengths across testing subjects.The deviations in the optimal coil scalp position and orientation were obtained by computing the differences from those in the FEM-CPO experiments.The deviation of the optimal E-field strength computes the deviations on the mean E-field strength in other methods from computational Efield strength under the coil placement determined by the other method.These deviation quantifications can be described as below equations: L Ma et al Moreover, even though the deep learning methods may not be able to duplicate the performance of the computational optimization perfectly, it may be possible to have a relatively large E-field strength in the optimal configuration.One way to qualify the deviation in the E-field strength is to evaluate whether the percentage of the mean E-field strength in different optimal choices would exceed 90% of the optimal E-field strength from the FEM-CPO results.Another way is to qualify the percent to include the same optimal coil placement as the FEM-CPO results in the top N choices from the deep learning methods.We added these two aspects in analyzing the optimal coil placement.

In-vivo TMS stimulation experiments for the motor region
To verify on the practical effectiveness of presented methods, in-vivo TMS stimulation on the righthand FDI muscle was conducted and to find out the electrophysiology responses of the motor target [42,43] under different optimal coil placements.27 recruited participants (17 participants did not receive the second session) received three successive TMS stimulation sessions.
The purpose of the first session was to determine the resting motor threshold (rMT) for a subject on individual hotspot location, the location with the largest MEP responses.The TMS coil was oriented in the PA45 • direction (45 • relative to the posterior-toanterior direction), which was approximately perpendicular to the central sulcus.The individual hotspot location was searched by a regular grid of 5 × 5 coil positions with 1 cm spacing.This grid was centered on the nearest scalp position above the target region [42].The average MEP amplitudes were calculated for each grid position.The hotspot location was defined as the position with the maximum MEP amplitude obtained from surface fitting the amplitude of MEPs on a 5 × 5 grid using the 'fit' command in Matlab software.Based on this position, the rMT was searched under the standard criterion that the minimum simulation intensity needed to evoke would be no less than 5 out of 10 MEPs with an amplitude exceeding 50 µV [49].
The second session aimed to test the associations between the spatial distribution of the FDI MEPs on a 5 × 5 grid with the spatial distribution of the Efield predicted by different methods.Only three orientations (PA0 • , PA45 • and PA90 • ) were considered in this stimulation experiment.Because they might not tolerate the process if the stimulation time was too lengthy for subjects.Based on the same coil positions in the 5 × 5 grid that were used to search the hotspot, each coil position received three single-pulse stimuli at an intensity of 120% rMT.Random stimuli were applied to the grid positions to avoid close positions from receiving successive stimulations [50].The average MEP amplitude was calculated for each coil position.Then, the averaged MEP amplitudes of a subject on each grid locations in all three directions were concatenated into a 75-dimension vector.These MEP vectors were used to tested the correspondence between target E-field strengths and the magnitude of the neural responses from the FDI brain target.
The third session was designed to test the optimization performance of the coil placements derived from four different groups: computational methods using (FEM-CPO) and (ADM-CPO) and deeplearning methods (RawRES-CPO and AnaRES-CPO).Top five coil placements with larger E-field strength within the target regions were selected for the deep learning methods (RawRES-CPO and AnaRES-CPO).For a subject, these coil placements were stimulated with 10 successive stimuli at an intensity of 120% rMT (110% rMT, based on the pain reactions of the participants).The stimulus sequence was arranged in random order with one-minute breaks to avoid as much as possible the influence of previous TMS stimulations [50].

High-precision real-time E-field simulations
After training on the HCP dataset, we first compared pointwise E-field strength differences between two deep learning (AnaRES and RawRES) and computation FEM methods under different coil placements across the scalp region (shown in figure 3).As shown in figure 3(b), the MAE (MAE = 0.059 ± 0.019 V m −1 ) of the AnaRES simulations under different coil poses was significantly (p < .001)lower than that (MAE = 0.081 ± 0.021 V m −1 ) of the RawRES simulations.The spatial distribution of the MAE was relatively stable, except for parts of the frontal scalp regions and the vicinity of the ear.The region with the maximum MAE reduction (figure 3(c)) for the AnaRES simulations occurred in the region with the higher absolute errors in the RawRES simulations.The E-field estimations in all five anatomical tissues had significantly (p < .001)lower MREs in the AnaRES method (WM: 11.2%, GM: 11.6%, CSF: 21.6%, bone: 8.5%, skin: 9.4%, average: 12.4%) than in the RawRES simulations (WM: 13.6%, GM: 14.7%, CSF: 27.9%, bone: 12.1%, skin: 12.4%, average: 16.1%), shown in figure 3(d).A spatial error visualization indicates that the reduction in MRE for GM occurs for almost all coil position on scalp, not for certain scalp regions in figure 3(e).Above all, AnaRES E-field simulations could provide a better precision in E-field simulation.
Another aspect is to evaluate the distribution of MAE in one simulation.Because the E-field simulation computed by computational method exists a sharp transition between adjacent brain tissues, it might cause anatomically related distribution of estimation errors in a deep-learning E-field estimations [26].These uneven errors may cause a large shift of optimal coil placements from results of the FEM method.Based on the results, brain structures learned by the AnaRES method match well with ground truth.The overlap between learned segmentation and ground truth is Dice = 0.974 on average within the local brain region (examples on HCP testing subjects are shown in supplementary figure S5).There is also a significant reduction in anatomically related errors (AAE = 0.107 ± 0.029 for the AnaRES and AAE = 0.137 ± 0.034 for the RawRES; paired t test, p < .001),as shown in figures 3(f) and (g).The reduction indicated that E-field simulation errors between the AnaRES and FEM results were more uniform and less related to brain anatomy.

Lower deviations in optimal coil placement from the FEM method
We show an example of coil placement optimization on the inferior frontal junction (IFJ) area (figure 4) in the BNA [32].We optimized the coil placement for the IFJ area on five testing subjects, independent from the training stage.The visualization shows that AnaRES-CPO had more consistent results with the optimal FEM positions and orientations across the subjects than the RawRES-CPO (figure 4).
Brain targets from three parcellations were used to test the reliability of different methods in the coil placement optimization.The spatial distribution of deviations (figure 5(a)) in both coil orientation and position were relatively evenly distributed for the cortical areas in the AnaRES-CPO results.Some of the areas with large cortex-to-scalp distances showed large deviations in orientation, position, and E-field strength.These deviations are shown for the RawRES-CPO method in supplementary figure S6.For the optimal coil orientation aspect (figure 5(b)), AnaRES-CPO had a significantly lower orientation deviation (12.53 • ± 16.64 • ) than RawRES (16.11 • ± 19.08 • ; paired t test, t1005 = 5.82, p < .001) on average.In more than 80% of the optimization cases, the orientation deviations in the AnaRES-CPO results were less than or equal to 15 • .For the optimal coil position aspect (figure 5(c)), AnaRES-CPO also showed a significantly lower position deviation (6.08 mm ± 6.69 mm) than RawRES-CPO (6.92 mm ± 6.74 mm; paired t test, t1005 = 3.37, p < .001) on average.We showed that the AnaRES-CPO results had fewer differences in orientation and position from the FEM-CPO results.
To assess the E-field strength within targets, we first took the FEM-CPO results as the gold standard, and compared the FEM E-field strength of the optimal coil placements from three different schemes.As shown in figure 5(d) the AnaRES-CPO had lower deviations in optimal E-field strength (MRE = 3.90%) than RawRES-CPO (MRE = 4.32%; paired t test, t1005 = 3.12, p < .001) on average.We also quantified the other two scores to show the similarity between the presented deep learning methods and the FEM method.In the top 5 optimal coil placements, there were 46% (RawRES-CPO) and 51.2% (AnaRES-CPO) chances, respectively, including the same optimal coil placement in the FEM-CPO results (shown in figure 5(e)).Moreover, even if the E-field strengths optimized by the AnaRES-CPO were not the largest compared with the FEM-CPO result, it still remains relatively higher E-field strengths for TMS stimulation.As shown in figure 5(f) there were more than 93.8% (95.5%, and 95.5%) chances to exceed 90% optimal E-field strength of the FEM-CPO in the top 1 (top 3 and top 5) choice in the AnaRES-CPO results.In other words, the AnaRES-CPO had more likelihood to reach a satisfactory E-field strength under its optimized coil placements.

Accurate E-field simulation contributes to a higher amplitude of MEPs
To further investigate the real-time estimation in practical application, in-vivo motor stimulation experiments were conducted on 27 subjects.By targeting the individual brain region corresponding to right-hand FDI muscles, optimal coil placements were generated by different deep-learning and computational schemes, and placed by a robotic arm in the TMS navigation system.The corresponding EMG signals were recorded during the stimulations.
We first investigated the associations between the distribution of the E-field strengths predicted by deep learning methods within the target and the distribution of corresponding electrophysiological responses recorded from EMG signals.We stimulated 120% rMT in three directions (PA0, PA45, PA90) on a 5 × 5 grid position in the experiment.As  In terms of working memory, the AnaRES-CPO method could be implemented with less than 2 GB memory on the laptop.Parallel computations using multiple threads can reduce the total time consumption but will have a higher working memory at cost.

Discussion
We developed an AnaRES method to simulate Efield distribution across the majority of scalp regions.The presented method directly uses individual T1w images as input to generate local E-fields near the coil position without time-consuming preprocessing.The method largely accelerates the traditional coil optimization process from several hours using FEM-based computational methods to several minutes with the presented deep learning method.By incorporating attentional layers in deep learning architecture and robust brain tissue representations, In in-silico experiments, the precision of E-field simulation is significant increased than previous deep-learning method based on the U-Net network.The in-vivo right-hand experiment demonstrated no significant difference in MEP magnitudes from AnaRES-optimized coil placements compared to those from computational methods, but higher than those from the compared deep-learning method.To the best of our knowledge, our study is the first to verify the effectiveness of deep-learning E-field simulation methods by examining individual neural responses to TMS stimulations on the optimized coil position.The study provides a robust foundation for the application of deep-learning methods in the future TMS modulation and clinical treatments.Conventional methods for choosing TMS coil placements (both coil position and orientation) are mostly based on empirical rules.There are empirical coil position rules including the 5 cm empirical rule approved by the U.S. food and drug administration [51,52], and an international 10-20 system-guided localization protocol recommended by the clinical TMS society [53][54][55][56][57] as well as empirical coil orientation rules that the TMS coil orientation is often set in the PA45 direction when evoking motor potentials [12,58,59].However, the topography of the cerebral cortex, especially for the dorsal prefrontal cortex is quite complicated and exhibits considerable inter-individual variability [60].Modulation on these regions necessitates more precise individual coil placements, given the large variation on individual targets, rather than employing fixed TMS coil settings [61][62][63].To achieve a personalized coil position for individual brain target, faster E-field simulation presented by AnatRES method is necessary and hopefully to mitigate the cross-subject variations on the stimulation performance.For instance, the magnitude of MEPs in the in-vivo experiment (shown in supplementary figure S8) using the 5 cm empirical rule (std = 405 µV) is higher variation than AnaRES-CPO results (std = 301 µV) and maintains high stimulation performance.Such E-field-based optimization is not limited by the motor target, and can be replicated to any arbitrary brain region which may lack region-specific stimulation experiences and rules.
The development of the E-field computational modeling allows to qualify individual E-field measurements under a particular coil placement prospectively [22,44].However, the FEM method to approximately solve the PDEs requires extensive image preprocessing and simulation time under different coil placements [15] which prevents the rapid application of TMS stimulation after MRI acquisition for a subject.Despite improvements in time efficiency of E-field simulations by the ADM method [14], its consideration only for E-field strength of brain targets, might have difficulties in the application on other E-field measurements, such as the on-target value [10] which also considers non-target areas to provide more stimulation specificity.Moreover, it is inconvenient for both subjects and medical staff to wait several hours after MRI acquisitions.In contrast, the presented method can directly optimize the coil placements of a brain target in less than several minutes (e.g.960 choices in less than 3 min using onethread parallel computing).Furthermore, it is capable of being generalized to other E-field measurements to fulfill optimization objectives because local brain region around a target is completely There is compatible performance between the magnitudes of evoked MEPs in optimal coil placements from the presented method and the computational method in our motor stimulation experiment.The uncertainty of the AnaRES-CPO E-field strength showed similar tendencies for the FEM-based computational method with increasing registration errors of coil position and orientation (supplementary figure S9).The presented method provides a time-efficient, precise and convenient way for multiple E-field measurements as well as multiple brain targets, especially those sharing a particular brain circuit [11,64,65].
Several deep-learning methods were built to supplant the time-consuming TMS pipeline in preprocessing steps or the FEM-based simulations [19,28].For example, Yokota et al [26] has significantly accelerated the simulation time through a U-Net architecture to directly map T1w images into E-field distribution.Skipping the reconstruction of individual brain anatomies with distinct conductivities allows the model to achieve greater time efficiency in E-field simulation.However, it poses a risk to the model precision under various image acquisition conditions as well as to the subsequent determination of TMS coil placement.In our study, the presented method adopts the explicit reconstruction (similar to computational method) and attentional architecture.The AnaRES method could learn a precise and stable E-field mapping to approximately solve to the PDEs under different coil settings.Comparison showed a significant reduction in E-field errors (21.1% reduction in GM MREs) as well as in anatomically related errors (21.8% reduction in anatomically related errors).Elevated precision in Efield simulation contributes to reduced deviations in optimal coil orientation, position, and E-field intensity.When directly applied deep-learning models into individual MRIs which acquisition is different from those in the training dataset, in-vivo experiment indicated the AnaRES method has shown an increasing correlation (r = 0.27) between target Efield strength and corresponding electrophysiological responses, 35% higher than the RawRES method (r = 0.20).It might be the reason to yield a significant higher value in evoked MEPs under the AnaRES optimal coil placements during motor stimulation.
The optimal E-filed strengths on the same target estimated from HCP dataset and the motor stimulation dataset (shown in supplementary figure 10) also indicates that there are fewer estimating biases resulting from the influence of image acquisitions in the present method.But it still leaves much space to improve the model generalization.
While the AnaRES method offers greater convenience for clinical practice and functions as on-device software for a navigation system, it also presents several limitations.Firstly, we only simulated the Efield at a 40 mm depth based on the effective distance of a Magstim double coil.The depth needed to changed and to be re-trained for other kinds of TMS coils, such as the H7 coil for deep TMS [18].Additionally, although the simulated E-field strength in the presented method demonstrates strong consistency with computational methods, the direction of the E-field is not considered, which could impede further exploration of direction-related E-field influence on the neural responses.Furthermore, the brain connectome and microstructure features (e.g.varying directions of fiber axons) from individual diffusion MRIs which may affect the local and remote stimulation performance, are not included in the presented deep-learning method.And the E-field simulation precision using deep learning method is still lower than new developed method with a realistic head model, such as ACA-ADM [66], leaving more space for the improvement of deep-learning-based method.
As we move forward, there still have a large space for image-based method to choose optimal coil placement.A by-product finding is that under the same stimulation intensity, the stimulated MEP magnitudes in experimental hotspot location are much higher than those in current image-based methods across ten subjects (shown in supplementary figure S11).Several key areas of development can be envisaged for both deep-learning-based or computational E-field methods.Firstly, improving the method to account for high-focality TMS coil types [67], highprecise head models [68], non-target region [10] will enable a more comprehensive understanding of the stimulation effects and their relationship with various coil configurations.Secondly, incorporating biophysical mechanisms, such as the effects of evoked E-field orientations and the diverse directions of fiber axons, the distribution of cell types might enhance the TMS stimulating performance and provide the explanatory power of the model [69,70].The neural responses for stimulations may also help to inspire efficient MRI image features for better coil placement in return.Thirdly, expanding the model to accommodate other modalities, such as functional and diffusion MRIs, and using the brain connectome to precisely delineate the individual brain target could provide a more holistic understanding of the individual brain's structural and functional organization.This, in turn, could pave the way for new advances in the treatment of various neurological and psychiatric disorders.

Figure 1 .
Figure 1.The workflow of various E-field simulation methods.The explicit and implicit tissue representations were extracted in the E-field simulation based on the computational and deep-learning methods.

Figure 2 .
Figure 2. Schematic diagram of the AnaRES method.(a): Local sampling based on coil placements (scalp positions and orientations).(b): The attentional U-Net architecture.(c): The construction and training of the attentional architectures.(d): The application of AnaRES E-field modeling for coil placement optimization.

Figure 3 .
Figure 3. E-field simulation errors for the RawRES and AnaRES methods.(a): Coil locations used for network training and validation.(b): General performance of the two methods in terms of E-field simulation results.(c): Method comparison of the spatial distribution of MAE scores.(d): Differences in MRE scores for brain tissues.(e): Spatial distribution of gray matter MRE scores for the two methods.(f): Average anatomically related error for the RawRES and AnaRES methods.(g): An example of anatomically related error.Border of the anatomical structure (left) and distribution of MAEs for the RawRES (middle) and AnaRES (right) methods.The * * * indicates the significance p < 0.001.

Figure 4 .
Figure 4. Five optimization examples using the inferior frontal junction (IFJ) regions from five individual subjects to illustrate the deviations in coil position (P) and orientation (θ) between the AnaRES-CPO and FEM-CPO methods.The first and the second rows show the coil placement optimized by RawRES-CPO and AnaRES-CPO.The third row shows the results from FEM-CPO.The fourth row presents a comparison of the optimal coil placements for all the methods on the cortical surface; the shaded color indicates the individual's IFJ region; the arrow represents the coil orientations.AnaRES-CPO: Red; RawRES-CPO: Green; FEM-CPO: Black.

Figure 5 .
Figure 5. Deviations in optimal coil position and orientations between the FEM-CPO and the other two deep-learning methods.(a): Spatial distribution of the optimal coil placements for AnaRES-CPO.(b): Orientation deviations for the two methods and their distributions.(c): Position deviations for the two methods over different targets.(d): Deviations in the optimal E-field strength.(e): Percent of the inclusion of the FEM-CPO optimal coil placement in the top N coil placements of the two deep learning methods.(f): Percent exceeding the 90% maximum E-field strength of the FEM-CPO method.* * * : p < .001,* * : p < .01,* : p < .05.

Figure 6 .
Figure 6.TMS stimulation experiments evoking the FDI MEPs based on different coil placements.(a).Illustration of the 5 × 5 grid positions and three coil orientations on a subject.(b).Pearson's correlation between the amplitude of the MEPs and the E-field strength within the target predicted by the AnaRES and RawRES simulations separately.(c).Illustration of four groups of coil placements on a subject.(d).Average amplitude of the FDI MEPs under different coil placements.The highest average amplitude is shown among the top N placements.ADM-CPO (red), FEM-CPO (orange), AnaRES-CPO (green), RawRES-CPO (blue).