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.

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–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–11]. To address this issue in part, many efforts have been pursued to pinpoint optimal TMS coil placements tailored to individual brain targets [12–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–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].

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.

Download figure:

Standard image High-resolution image

Recent research indicates that deep learning methods have the potential to simulate the E-field statistically [26–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 E-fields 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-learning-based method using explicit, supervised brain tissue representation as the intermediate variable (figure 1(a)), termed the anatomically aware, real-time 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.

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 E-field inference using an attentional U-net architecture (figure 2(c)).

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.

Download figure:

Standard image High-resolution image 2.1. Local sampling based on coil placements

In a TMS navigation system, the scalp position 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 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 , 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 . A transformation based on coil placements 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 .

2.2. 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 E-field 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 max-pooling 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 low-level 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].

2.3. 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 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 of coil placements is set in the range of coil angle () as well as in the range of coil position shown in figure 3(a). After that, the network could learn a mapping from the input T1w images 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 .

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.

Download figure:

Standard image High-resolution image

The entire deep-learning network can be described as a function to optimize the parameters in equation (1). The first term in the objective for the network was to generate E-fields ) that are closed to the E-fields 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 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)

The anatomical tissues and E-fields objectives were optimized in a two-step way. First, the anatomical objective was achieved by optimizing the cross-entropy loss () between the predicted anatomical structures () and anatomical distributions from realistic head model . Then the network achieved the regression objective through the optimization on the mean square error () between the network simulations ) and the high-precision simulations ). 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 . 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.

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 E-fields 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.

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

3.2. Brain targets and coil placement configurations3.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–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-to-cortex 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.

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