## 2 Materials and Methods
### Implementation of the DBIM-FISTA with a GPU-Based FDTD Forward Solver  
We have used CUDA toolkit to accelerate our FDTD forward solver in GPU similar to [49], which focused however on modeling RF wave interactions in high-field MRI. The FDTD algorithm consists of electric and magnetic field updates, and each can be viewed as a kernel on GPU. The kernel, executed as a grid in GPU, can be divided into multiple blocks of threads, whose number can be adjusted based on the GPU platform. We can therefore increase the computational efficiency significantly by computing the blocks of threads simultaneously. Our algorithm is designed to use the high performance of GPU without changing the MATLAB based DBIM code, to take advantage of MATLAB's capability for matrix computing. To this end, we have implemented the 3-D FDTD algorithm in C++ with CUDA and then incorporated it with MEX functions in MATLAB, which is also used to implement the FISTA inverse solver. The environment of our implementation is MATLAB 2020 (Mathworks Inc., Natick, MA, USA), CUDA 10.2 (Nvidia Co., Santa Clara, CA, USA) and Visual Studio 2015 (Microsoft Co., Redmond, WA, USA) with Intel(r) Xeon(r) CPU E5-2640 v3 @ 2.60 GHz and GPU Tesla K20c with 5 GB memory. The maximum number of threads for each block is 1024 and we use a 2-D block with size \(32\times 32=1024\) to make full use of each block. The number of threads for each block should be set to a multiple of 32 as the working unit in the GPU is a warp, which is a set of 32 threads [50]. Similarly, we use a 2-D grid GPU with size \(\frac{lmax}{32}\times\frac{lmax}{32}\) for each kernel, \(E_{i}\) (\(i=x,y\) or \(z\)) for example, and the loop inside each kernel for updating the fields is in the following structure shown in Algorithm 1, where Imax, Jmax and Kmax represent number of voxels for the \(x,y,\) and \(z\) dimension, respectively. ```
\(ti\gets blockIdx.x*blockDim.x+threadIdx.x\) \(tj\gets blockIdx.y*blockDim.y+threadIdx.x\) if\(ti<lmax\) and \(tj<lmax\)then fork\(tk\gets 1\) to \(Kmax\)do \(E_{i}[ti+tj*Imax+tk*(Imax*Jmax)]\gets UpdateE_{i}\) endfor endif
```  
**Algorithm 1** Kernel of updating electric field  
The GPU Tesla K20c has 13 streaming processors (SMs) and for each SM there are at most 2048 threads. Thus 13 \(\times\) 2048 threads can be run simultaneously. The flowchart of the 3-D FDTD algorithm is shown in Figure 1. To reduce computational times, most of the FDTD functions are run on GPU (kernels), including the source update. At each DBIM iteration, the FDTD variables are first allocated space on GPU after initialization by MEX functions. Thus all FDTD-related variables are on GPU memory. Then the FDTD part on GPU starts to simulate the wave propagation for each transmitter antenna, and the calculated data to be used in DBIM is transferred back to MATLAB via MEX functions, which can be used directly by the MATLAB code. Finally the ill-conditioned linear system is built by (3) and solved by FISTA to obtain the reconstructed values. The discretization size for FDTD is defined as  
\[\Delta h=\frac{\lambda_{f}}{n_{d}}, \tag{33}\]  
where \(\lambda_{f}\) is the wavelength and \(n_{d}\) is the discretization step. To ensure accuracy and reduce numerical dispersion, \(n_{d}>10\) is required. While \(\Delta h\) is 30 mm in air, its value inside the brain is defined by a much smaller wavelength (\(\lambda_{f}=\frac{300}{\sqrt{e}}\)), thus requiring a greater \(n_{d}\) and a smaller \(\Delta h\). We have therefore selected a cubic grid voxel size of 2 mm for each dimension in all our 2-D and 3-D FDTD simulations. This value can balance accuracy and computational burden. For example, using a 2 mm vs. a 1 mm cubic voxel side results in field values with a mean relative difference of \(10^{-3}\) but requires an 8 times bigger grid. The time step \(\Delta t\) defined as  
\[\Delta t=\frac{\Delta h}{2c}, \tag{34}\]  
where \(c\) is the speed of EM wave, will also increase for a higher resolution model and thus lead to more iterations.