Magnetic Resonance-Electrical Properties Tomography by Directly Solving Maxwell’s Curl Equations

: Magnetic Resonance-Electrical Properties Tomography (MR-EPT) is a method to reconstruct the electrical properties (EPs) of bio-tissues from the measured radiofrequency (RF) ﬁeld in Magnetic Resonance Imaging (MRI). Current MR-EPT approaches reconstruct the EP proﬁle by solving a second-order partial di ﬀ erential wave equation problem, which is sensitive to noise and can induce large reconstruction artefacts near tissue boundaries and areas with inaccurate ﬁeld measurements. In this paper, a novel MR-EPT approach is proposed, which is based on a direct solution to Maxwell’s curl equations. The distribution of EPs is calculated by iteratively ﬁtting the RF ﬁeld calculated by the ﬁnite-di ﬀ erence-time-domain (FDTD) technique to the measured values. To solve the time-consuming problem of the iterative ﬁtting, a graphics processing unit (GPU) is used to accelerate the FDTD technique to process the ﬁeld calculation kernel. The new EPT method was evaluated by investigating a numerical head phantom, and it was found that EP values can be accurately calculated and were relatively insensitive to simulated RF ﬁeld errors. The feasibility of the proposed method was further validated using phantom-based experimental data collected from a 9.4 Tesla (T) Magnetic Resonance Imaging (MRI) system. The results also indicated that more accurate EP values could be reconstructed from the new method compared with conventional methods. Moreover, even in the absence of phase information of the RF ﬁeld, the proposed approach is still capable of o ﬀ ering robust EPT solutions, thus having improved practicality for potential clinical applications. electric ﬁeld components in both space and time domains. The explicit ﬁeld updating procedure was easily implemented in multiple, synchronized threads within the parallel-computing framework of a GPU. In this work, all computations have been executed on an Intel Xeon CPU E5-2680 V2 @ 2.80 GHz (128 GB RAM) equipped with an NVIDIA TITAN V GPU (12 GB RAM, 5120 CPUs, and 1455 MHz clock frequency) and more details about the GPU-FDTD computing can be found in Reference [29].


Introduction
Magnetic Resonance-Electrical Properties Tomography (MR-EPT) is an imaging technique which maps the electrical properties (permittivity and conductivity) of a subject by analyzing the radiofrequency (RF) magnetic field measured during an MR scan. This technique has shown promise for clinical diagnosis and interests in ultra-high field Magnetic Resonance Imaging (MRI), because of its potential in the evaluation of patient-specific absorption rate (SAR) [1]. MR-EPT has also 2. Theory and Method Figure 1 shows the frameworks of different EPT methods. Conventional differential-based methods (Figure 1a,b) rely on inverse processes involving second-order differential kernels in the EP reconstruction. In the proposed approach (Figure 1d), a forward process is used, which investigates the relationship between the EPs and the amplitude of B + 1 data, and only involves first-order derivative; thus, it is Laplacian-operator free. Details of differential equations for the description of radiofrequency problems in MRI can be found in the literature [18]. When the other conditions are pre-known (such as the structures of the coil and subject), the amplitude of the B + 1 field is mainly determined by the EPs of the imaging subject. Therefore, it is possible to optimize the EPs during the forward process so that the calculated amplitude of the B + 1 field can match the measured one. The efficiency and accuracy of this approach rely on the forward solver and the optimizer. Because the EPs are iteratively optimized during each forward process, an efficient forward EM solver is required to achieve a practical solution. As opposed to conventional methods, the proposed method calculates the EP values using the full set of magnetic field components available from the EM solver, without neglecting the z-component of the magnetic field (B z ) [18]. The details of the forward solver and optimization process are given below.
Appl. Sci. 2020, 10, x 3 of 13 radiofrequency problems in MRI can be found in the literature [18]. When the other conditions are pre-known (such as the structures of the coil and subject), the amplitude of the B + 1 field is mainly determined by the EPs of the imaging subject. Therefore, it is possible to optimize the EPs during the forward process so that the calculated amplitude of the B + 1 field can match the measured one. The efficiency and accuracy of this approach rely on the forward solver and the optimizer. Because the EPs are iteratively optimized during each forward process, an efficient forward EM solver is required to achieve a practical solution. As opposed to conventional methods, the proposed method calculates the EP values using the full set of magnetic field components available from the EM solver, without neglecting the z-component of the magnetic field ( z B ) [18]. The details of the forward solver and optimization process are given below.  [7,8], (b) the second-order differential wave equation-based solution [16][17][18][19][20], (c) integral equation-based solution [11,[23][24][25][26], and (d) the proposed EPT approach. It is noted that the proposed method finds the EPT solution without accessing the phase information of the radiofrequency (RF) field.

GPU-Accelerated FDTD Solver
The in-house built FDTD solver [29] was constructed based on the differential form of the Maxwell equations (ME). As the first-order central differential scheme was used to implement the partial derivatives in the ME numerically, the noise effects were only amplified in the first-order rate,  [7,8], (b) the second-order differential wave equation-based solution [16][17][18][19][20], (c) integral equation-based solution [11,[23][24][25][26], and (d) the proposed EPT approach. It is noted that the proposed method finds the EPT solution without accessing the phase information of the radiofrequency (RF) field.

GPU-Accelerated FDTD Solver
The in-house built FDTD solver [29] was constructed based on the differential form of the Maxwell equations (ME). As the first-order central differential scheme was used to implement the partial derivatives in the ME numerically, the noise effects were only amplified in the first-order rate, which was more robust than conventional EPT methods based on second-order wave equation. In the FDTD computational domain, a uniform mesh was used to discretize the subject-coil space, and specific cell sizes were problem-dependent. The calculation volume was truncated with perfectly matched layers (PML). The PML was composed of 8 to 10 layers, and the conductivities of each layer were calculated based on the method presented in Reference [30].
In the FDTD solver, the vector forms of the electric and magnetic fields were calculated for the evaluation of B + 1 = B x + iB y /2 [31]. Here, as Maxwell's curl equations state, the calculated B + 1 fields were a function of EPs of the imaged subject. For more details about the in-house developed FDTD solver, please refer to our previous work [29].
In the constructed optimization problem, the forward process was implemented in each iteration step; thus the FDTD solver was required to run efficiently for the evaluation of the B + 1 field. To this end, a GPU-accelerated FDTD solver [29] was applied to the proposed EPT solution procedure. In this study, a 12-rung birdcage coil was modelled as an EM source and loaded with a subject (phantom or head), as shown in Figure 2. A perfect electric conductor (PEC) was placed around the birdcage coil to simulate the RF shield, and PML boundary conditions were placed at the top and bottom of the computation domain. In the FDTD simulation, Maxwell's curl equations are approximated with Yee's 'leap-frog' finite-difference algorithm, and the entire computational space is translated from standard Cartesian space to Yee-cell space, where the FDTD equations exhibit inherent strong parallelism. At each Yee-cell, this parallelism permits the simultaneous calculation of magnetic and electric field components in both space and time domains. The explicit field updating procedure was easily implemented in multiple, synchronized threads within the parallel-computing framework of a GPU. In this work, all computations have been executed on an Intel Xeon CPU E5-2680 V2 @ 2.80 GHz (128 GB RAM) equipped with an NVIDIA TITAN V GPU (12 GB RAM, 5120 CPUs, and 1455 MHz clock frequency) and more details about the GPU-FDTD computing can be found in Reference [29].
Appl. Sci. 2020, 10, x 4 of 13 which was more robust than conventional EPT methods based on second-order wave equation. In the FDTD computational domain, a uniform mesh was used to discretize the subject-coil space, and specific cell sizes were problem-dependent. The calculation volume was truncated with perfectly matched layers (PML). The PML was composed of 8 to 10 layers, and the conductivities of each layer were calculated based on the method presented in Reference [30].
In the FDTD solver, the vector forms of the electric and magnetic fields were calculated for the evaluation of + 1 = + 2 x y B B iB [31]. Here, as Maxwell's curl equations state, the calculated B + 1 fields were a function of EPs of the imaged subject. For more details about the in-house developed FDTD solver, please refer to our previous work [29].
In the constructed optimization problem, the forward process was implemented in each iteration step; thus the FDTD solver was required to run efficiently for the evaluation of the B + 1 field. To this end, a GPU-accelerated FDTD solver [29] was applied to the proposed EPT solution procedure. In this study, a 12-rung birdcage coil was modelled as an EM source and loaded with a subject (phantom or head), as shown in Figure 2. A perfect electric conductor (PEC) was placed around the birdcage coil to simulate the RF shield, and PML boundary conditions were placed at the top and bottom of the computation domain. In the FDTD simulation, Maxwell's curl equations are approximated with Yee's 'leap-frog' finite-difference algorithm, and the entire computational space is translated from standard Cartesian space to Yee-cell space, where the FDTD equations exhibit inherent strong parallelism. At each Yee-cell, this parallelism permits the simultaneous calculation of magnetic and electric field components in both space and time domains. The explicit field updating procedure was easily implemented in multiple, synchronized threads within the parallel-computing framework of a GPU. In this work, all computations have been executed on an Intel Xeon CPU E5-2680 V2 @ 2.80 GHz (128 GB RAM) equipped with an NVIDIA TITAN V GPU (12 GB RAM, 5120 CPUs, and 1455 MHz clock frequency) and more details about the GPU-FDTD computing can be found in Reference [29].

Optimization
EPs were estimated through an iterative optimization procedure.

Optimization
EPs were estimated through an iterative optimization procedure. The optimization variables (EPs) are x EP = [ε tissue σ tissue ], that is, the permittivity and conductivity values of each voxel of the imaging subject. Using the developed FDTD solver described in Section 2.1, the tissue dielectric properties variable x EP were mapped to the B + 1 field as: where F FDTD denotes the FDTD process.

of 13
Then, a cost function which describes the residual between the mapped B + 1 fields and the measured B + 1 field can be built as: Therefore, an unconstrained optimization problem which aims to minimize the cost function can be formed as: min The problem in Equation (3) can be efficiently solved by using optimization methods such as quasi-Newton or trust-region techniques. The EPs of the imaging subject can be obtained after the optimum x EP is found from Equation (3).
In this work, the genetic algorithm (GA, standard MATLAB function) was used to find the optimal solution of x EP . In the optimization procedure, an initial EP distribution denoted as x O was set based on the LHA. Inside the uniform region, the EP distribution was directly derived from the LHA-based result. On the tissue boundaries, average EPs of the neighboring tissues were used. This operation avoided the boundary artifacts in the LHA-based results. The variable range in the GA was set as Alternatively, the initial EPs of each tissue can be set based on the data from Reference [32], and the variable range was set as [33]. The GA-based optimization was implemented to satisfy one of the following stop criteria: (1) the fitness, C, was less than the expected fitness (least-square residual error), C desired , and (2) the number of iterations exceeds certain pre-set limits, N desired . Once condition (1) or (2) was satisfied, the optimization process could be terminated. The parameter setup of C desired and N desired is problem-dependent. In this study, it was observed that the variation of residual error was negligible after 50 iterations, and when it reached a value of around 10 −3 . Therefore, we chose C desired and N desired as 10 −3 and 100, respectively.

Simulation Study
The proposed EPT method was first evaluated using a simple phantom. The phantom was composed of an inner compartment with a diameter of 15 mm and an outer compartment with a diameter of 40 mm. The height of the inner and outer compartments was 50 mm. The relative permittivity for the inner and outer compartments was 78.5, and the conductivities were 0.43 S/m and 0.78 S/m, respectively. After the phantom evaluation, the head of the realistic human model Duke, from the virtual family [32], was used to study the proposed EPT method. The initial model resolution of 1 mm 3 was re-meshed as 2 mm 3 , as a trade-off between accuracy and computational efficiency of the forward solver. The electrical properties of the brain tissues at 298 MHz (the operating frequency of 7T MRI) were determined based on the reports from Reference [33]. The head model was placed at the center of the coil ( Figure 2). To assess the noise robustness of the proposed method, white Gaussian noise (GWN) was added to the simulated B + 1 field data to achieve signal-to-noise ratios (SNR) between 15 and 55 dB [34].

Experimental Setup
To evaluate the practicality of the proposed EPT method, a phantom-based experiment at 9.4T preclinical MRI system (BioSpec 94/30 USR, Bruker, Germany) located at the Centre for Advanced Imaging, the University of Queensland, was conducted. A Bruker birdcage coil (Øouter/Øinner = 112/86 mm) was used to acquire images and the B + 1 field. The phantom described in Section 2.3 was fabricated, with conductivities of 0.43 S/m and 0.78 S/m for the inner and outer compartments respectively, and permittivity of 78 for both compartments. The two compartments were separated by a wall (thickness = 2 mm) made of acrylic material. The EPs of the liquid solution Appl. Sci. 2020, 10, 3318 6 of 13 were measured using a dielectric probe (Kit DAK-12, SPEAG Ltd., Zurich, Switzerland). The probe uncertainty and temperature-associated variation in dielectric properties were 3.3% and 2%, respectively. After careful calibration, the total measurement error can be considered less than 6.6%. The aqueous solutions consist of distilled water, NaCl, and CuSO4·5H2O with mass ratios of 100:0.18:0.025 and 100:0.5:0.025, respectively. Figure 3 shows the phantom and its cross-section, acquired with a gradient echo sequence. Two MR images with flip angles equal to 45 • and 90 • were obtained by using the fast low-angle shot (FLASH) sequence (echo time (TE) equals 5 ms and repetition time (TR) equals 5000 ms). The field-of-view in the experiment was set to 48 × 48 × 48 mm, and a 0.5 mm isotropic resolution was used in the image acquisition. The double-angle-method (DAM) [35,36] was applied to calculate the amplitude of the B + 1 field, while the phase was calculated by using the transceive phase assumption (TPA) [6,37]. To implement the TPA, two images with different TEs (TE1 = 2.5 ms and TE2 = 3.5 ms) were used to estimate the phase difference due to the inhomogeneous B o field. Another two images, generated using identical TE/TR but reversed gradient directions, were used to cancel the gradient-induced phase (e.g., eddy current) [6]. It should be noted that the B + 1 phase was not used in the EP reconstruction but only for the purpose of evaluating the accuracy of the forward solver, as illustrated in the following section.
Appl. Sci. 2020, 10, x 6 of 13 Imaging, the University of Queensland, was conducted. A Bruker birdcage coil (Øouter/Øinner = 112/86 mm) was used to acquire images and the B + 1 field. The phantom described in Section 2.3 was fabricated, with conductivities of 0.43 S/m and 0.78 S/m for the inner and outer compartments respectively, and permittivity of 78 for both compartments. The two compartments were separated by a wall (thickness = 2 mm) made of acrylic material. The EPs of the liquid solution were measured using a dielectric probe (Kit DAK-12, SPEAG Ltd., Zurich, Switzerland). The probe uncertainty and temperature-associated variation in dielectric properties were 3.3% and 2%, respectively. After careful calibration, the total measurement error can be considered less than 6.6%. The aqueous solutions consist of distilled water, NaCl, and CuSO4•5H2O with mass ratios of 100:0.18:0.025 and 100:0.5:0.025, respectively. Figure 3 shows the phantom and its cross-section, acquired with a gradient echo sequence. Two MR images with flip angles equal to 45° and 90° were obtained by using the fast low-angle shot (FLASH) sequence (echo time (TE) equals 5 ms and repetition time (TR) equals 5000 ms). The field-of-view in the experiment was set to 48 × 48 × 48 mm, and a 0.5 mm isotropic resolution was used in the image acquisition. The double-angle-method (DAM) [35,36] was applied to calculate the amplitude of the + 1 B field, while the phase was calculated by using the transceive phase assumption (TPA) [6,37]. To implement the TPA, two images with different TEs (TE1 = 2.5 ms and TE2 = 3.5 ms) were used to estimate the phase difference due to the inhomogeneous o B field. Another two images, generated using identical TE/TR but reversed gradient directions, were used to cancel the gradient-induced phase (e.g., eddy current) [6]. It should be noted that the + 1 B phase was not used in the EP reconstruction but only for the purpose of evaluating the accuracy of the forward solver, as illustrated in the following section.

Simulation Results
The simple phantom described in Section 2.3 was first used to evaluate the performance of the proposed EPT method. To implement the FDTD calculation, the Yee-cell size was 1 mm, and the entire computational domain was divided into Dx × Dy × Dz = 112 × 112 × 112 ≈ 1,404,928 cells. Each GPU-FDTD iteration time was approximately 18 seconds, and the total time for the optimization was around 30 min. In the evaluation, six different noise levels were added to the simulated + 1 B field.
As shown in Table 1, the EPs can be retrieved with high accuracy using the proposed approach, even under an extreme scenario in which the SNR level declined to 15 dB (the lowest SNR level used in current MR-EPT numerical evaluations was 35 dB). The maximum errors were 3% and 0.63% for the permittivity and conductivity reconstruction, respectively. These results show the strong robustness of the method to noise.

Simulation Results
The simple phantom described in Section 2.3 was first used to evaluate the performance of the proposed EPT method. To implement the FDTD calculation, the Yee-cell size was 1 mm, and the entire computational domain was divided into D x × D y × D z = 112 × 112 × 112 ≈ 1,404,928 cells. Each GPU-FDTD iteration time was approximately 18 seconds, and the total time for the optimization was around 30 min. In the evaluation, six different noise levels were added to the simulated B + 1 field. As shown in Table 1, the EPs can be retrieved with high accuracy using the proposed approach, even under an extreme scenario in which the SNR level declined to 15 dB (the lowest SNR level used in current MR-EPT numerical evaluations was 35 dB). The maximum errors were 3% and 0.63% for the permittivity and conductivity reconstruction, respectively. These results show the strong robustness of the method to noise. Following the simple phantom evaluation, the numerical head model, as described in Section 2.3, was used to assess the performance of the proposed EPT method. The target and reconstructed EP profiles are shown in Figure 4. It can be seen that accurate EPs can be retrieved by using the proposed EPT method.
Appl. Sci. 2020, 10, x 7 of 13 Following the simple phantom evaluation, the numerical head model, as described in Section 2.3, was used to assess the performance of the proposed EPT method. The target and reconstructed EP profiles are shown in Figure 4. It can be seen that accurate EPs can be retrieved by using the proposed EPT method. To further illustrate the effectiveness of the new method, the traditional algorithm based on LHA and a recently developed EPT solution using a modified Finite-Difference (mFD) scheme [18] were both used for comparison. As a recently developed second-order differential equation-based approach, the mFD approach is a good candidate for comparison between the conventional EPT method with the proposed one. The root-mean-square-error (RMSE) [16] was used to evaluate the performance quantitatively, and the corresponding results are shown in Figure 5. It can be seen from this figure that the new EPT method achieved significantly better results than those from compared methods across SNR levels from 15 to 55 dB. In the extreme case, in which the SNR level was 15 dB, both comparison methods obtained RMSE values larger than 1.5 for the permittivity reconstruction and 1.2 for the conductivity reconstruction. However, the proposed method was capable of retrieving more accurate EP's values with an RMSE less than 0.2 for both the permittivity and conductivity for all SNR levels. To further illustrate the effectiveness of the new method, the traditional algorithm based on LHA and a recently developed EPT solution using a modified Finite-Difference (mFD) scheme [18] were both used for comparison. As a recently developed second-order differential equation-based approach, the mFD approach is a good candidate for comparison between the conventional EPT method with the proposed one. The root-mean-square-error (RMSE) [16] was used to evaluate the performance quantitatively, and the corresponding results are shown in Figure 5. It can be seen from this figure that the new EPT method achieved significantly better results than those from compared methods across SNR levels from 15 to 55 dB. In the extreme case, in which the SNR level was 15 dB, both comparison methods obtained RMSE values larger than 1.5 for the permittivity reconstruction and 1.2 for the conductivity reconstruction. However, the proposed method was capable of retrieving more accurate EP's values with an RMSE less than 0.2 for both the permittivity and conductivity for all SNR levels. Appl. Sci. 2020, 10, x 8 of 13

Experimental Results
The new EPT method was further verified by using experimental data obtained from the 9.4T MRI system. The measured + 1 B field data were compared with that from the GPU-FDTD solver, as shown in Figure 6, indicating that an accurate + 1 B field can be calculated for both the magnitude and phase. These results confirmed the accuracy of the forward solver in the new EPT approach, thus showing that stable and fast convergence can be expected. The measured + 1 B data (including the magnitude and phase information) was firstly used to generate the EP profile by using the conventional method with the LHA. As shown in Figure 7b,f, large reconstruction artefacts can be clearly observed in the two solutions. In the inner and outer regions of the acrylic container, more distortions are observed

Experimental Results
The new EPT method was further verified by using experimental data obtained from the 9.4T MRI system. The measured B + 1 field data were compared with that from the GPU-FDTD solver, as shown in Figure 6, indicating that an accurate B + 1 field can be calculated for both the magnitude and phase. These results confirmed the accuracy of the forward solver in the new EPT approach, thus showing that stable and fast convergence can be expected. The measured B + 1 data (including the magnitude and phase information) was firstly used to generate the EP profile by using the conventional method with the LHA.

Experimental Results
The new EPT method was further verified by using experimental data obtained from the 9.4T MRI system. The measured + 1 B field data were compared with that from the GPU-FDTD solver, as shown in Figure 6, indicating that an accurate + 1 B field can be calculated for both the magnitude and phase. These results confirmed the accuracy of the forward solver in the new EPT approach, thus showing that stable and fast convergence can be expected. The measured + 1 B data (including the magnitude and phase information) was firstly used to generate the EP profile by using the conventional method with the LHA. As shown in Figure 7b,f, large reconstruction artefacts can be clearly observed in the two solutions. In the inner and outer regions of the acrylic container, more distortions are observed As shown in Figure 7b,f, large reconstruction artefacts can be clearly observed in the two solutions. In the inner and outer regions of the acrylic container, more distortions are observed compared with the Appl. Sci. 2020, 10, 3318 9 of 13 results from mFD and the proposed method. This is because the LHA-based method is susceptible to measurement noise; hence EP distortions still appear in those homogenous regions. Figure 7c,g showed the reconstruction results by using the mFD and Figure 7d,h are the solutions of the new EPT method. It can be observed that although the artifacts were reduced when using the mFD method, and the specific EPs were still deviated from the measured values. With the proposed method, an accurate estimation of the EPs in the two solutions can be achieved. The proposed method attempted to seek the optimum EP values of the inner and outer regions by fitting the measured and calculated B + 1 field. This process mainly considers the global B + 1 field distribution, thereby mitigating the local measurement noise effect and significantly improving the robustness of the algorithm. Owing to this feature, accurate EP values of the inner and outer regions can be calculated. Table 2 shows that permittivity and conductivity values can be calculated with less than 0.6% errors, which indicated the superiority of the proposed EPT in a practical utilization. The excellent results shown in Figure 7 are obtained with the available exact acrylic boundaries derived from the acquired MR image, and the EPs in segmented regions are taken as unknown constants to reconstruct through optimization. The precise position of the acrylic pixels is also known from segmentation and excluded in the EP reconstruction.
Appl. Sci. 2020, 10, x 9 of 13 compared with the results from mFD and the proposed method. This is because the LHA-based method is susceptible to measurement noise; hence EP distortions still appear in those homogenous regions. Figure 7c,g showed the reconstruction results by using the mFD and Figure 7d,h are the solutions of the new EPT method. It can be observed that although the artifacts were reduced when using the mFD method, and the specific EPs were still deviated from the measured values. With the proposed method, an accurate estimation of the EPs in the two solutions can be achieved. The proposed method attempted to seek the optimum EP values of the inner and outer regions by fitting the measured and calculated + 1 B field. This process mainly considers the global + 1 B field distribution, thereby mitigating the local measurement noise effect and significantly improving the robustness of the algorithm. Owing to this feature, accurate EP values of the inner and outer regions can be calculated. Table 2 shows that permittivity and conductivity values can be calculated with less than 0.6% errors, which indicated the superiority of the proposed EPT in a practical utilization. The excellent results shown in Figure 7 are obtained with the available exact acrylic boundaries derived from the acquired MR image, and the EPs in segmented regions are taken as unknown constants to reconstruct through optimization. The precise position of the acrylic pixels is also known from segmentation and excluded in the EP reconstruction.  In addition to EP reconstruction, the new EPT can also provide information about the z B field, which is usually neglected in current differential and integral-based MR-EPT methods. Figure 8 shows the exported z B components from the proposed solution procedure. As expected, in the central plane of the birdcage coil, the z B was quite small. Close to the coil ends, the z B value was elevated.
The availability of these magnetic field components should be helpful to improve the EPT solution and its applications in MRI.  In addition to EP reconstruction, the new EPT can also provide information about the B z field, which is usually neglected in current differential and integral-based MR-EPT methods. Figure 8 shows the exported B z components from the proposed solution procedure. As expected, in the central plane of the birdcage coil, the |B z | was quite small. Close to the coil ends, the |B z | value was elevated. The availability of these magnetic field components should be helpful to improve the EPT solution and its applications in MRI.

Advantages of the Proposed Method
It has been well recognized that most EPT methods developed so far are based on the secondorder differential wave equation. These methods are theoretically sound but practically un-robust in real imaging applications. In particular, these methods all rely on the accurate phase information of the + 1 B data, which is particularly challenging to obtain in ultra-high field MRI. To find EPT solutions, assumptions and simplifications have to be made and cause system errors. By directly solving Maxwell's equations related to EPT (in which first-order differential operation only is involved) and probing dielectric properties via analysis of the B + 1 data, the proposed method inversely finds solutions [38,39] with improved robustness and practicability for ultra-high field applications. This has been demonstrated in both simulation studies and phantom-based experimental testing at 9.4T. Results indicated that the proposed EPT offers a significant enhancement in the reconstruction quality compared with existing EPT approaches. Specifically, lower tissue boundary artefacts and less sensitive to measurement noise are observed, and acceptable EP reconstruction can be achieved with less than a 0.2 RMSE value when the SNR level reached 15 dB.
As shown in the experimental results, the proposed method can provide EPT solutions, as well as export full components of the magnetic and electrical fields. With the available z B profiles, one can see that near the birdcage center, the EPT reconstruction should be quite accurate without using the z B information. However, in the regions close to the birdcage end rings, neglecting the z B provide much better solutions when surface coils are used. Furthermore, the proposed EPT method can also export the electric fields, which cannot be measured inside patients. The availability of these field components will be particularly useful to the estimation of patient-specific SAR and also new development of imaging methods in ultra-high field MRI.

Limitations of the Proposed Method
Two main limitations can be identified regarding the proposed method. Firstly, the imaged subject needs to be modelled in the EM field solver. As it is challenging to find solutions to the formed optimization problem with millions of unknowns, the following assumption was made: one type of tissue has uniform EP values. In electromagnetic computing, a segmented MR image is required as prior information to help assign initial EP values of corresponding FDTD cells representing the

Advantages of the Proposed Method
It has been well recognized that most EPT methods developed so far are based on the second-order differential wave equation. These methods are theoretically sound but practically un-robust in real imaging applications. In particular, these methods all rely on the accurate phase information of the B + 1 data, which is particularly challenging to obtain in ultra-high field MRI. To find EPT solutions, assumptions and simplifications have to be made and cause system errors. By directly solving Maxwell's equations related to EPT (in which first-order differential operation only is involved) and probing dielectric properties via analysis of the B + 1 data, the proposed method inversely finds solutions [38,39] with improved robustness and practicability for ultra-high field applications. This has been demonstrated in both simulation studies and phantom-based experimental testing at 9.4T. Results indicated that the proposed EPT offers a significant enhancement in the reconstruction quality compared with existing EPT approaches. Specifically, lower tissue boundary artefacts and less sensitive to measurement noise are observed, and acceptable EP reconstruction can be achieved with less than a 0.2 RMSE value when the SNR level reached 15 dB.
As shown in the experimental results, the proposed method can provide EPT solutions, as well as export full components of the magnetic and electrical fields. With the available B z profiles, one can see that near the birdcage center, the EPT reconstruction should be quite accurate without using the B z information. However, in the regions close to the birdcage end rings, neglecting the B z component would introduce system errors as the intensity of B z component becomes comparable with that of B x and B y . This B z issue will be more prominent in array coil systems because of the larger value of |B z |/ B x,y in such cases, and it is therefore expected that the new EPT method can provide much better solutions when surface coils are used. Furthermore, the proposed EPT method can also export the electric fields, which cannot be measured inside patients. The availability of these field components will be particularly useful to the estimation of patient-specific SAR and also new development of imaging methods in ultra-high field MRI.

Limitations of the Proposed Method
Two main limitations can be identified regarding the proposed method. Firstly, the imaged subject needs to be modelled in the EM field solver. As it is challenging to find solutions to the formed optimization problem with millions of unknowns, the following assumption was made: one type of tissue has uniform EP values. In electromagnetic computing, a segmented MR image is required as prior information to help assign initial EP values of corresponding FDTD cells representing the imaged subject. The structure of the imaged subject can be constructed by using image registration techniques or implementing segmentation of the obtained MR images [40,41]. During optimization, several or tens of unknown variables are only required to update. With a significant dimension reduction, superior convergence can thus be achieved. Owing to an ideal experimental setup involving birdcage coil and sample phantom, the FDTD model can match the experiment shown in this work. Compared with conventional EPT approaches with z-field independency and LHA, the full-wave EM calculation is capable of regularizing the EPT solution closer to measurements. Even with the B + 1 data only, the strong coupling mechanism of E-and H-fields in Maxwell's curl equation can still faithfully find the EP values. Lastly, the proposed EP solution is sensitive to the initial setup of the EP values. We expect that the accuracy of the proposed method for human tissue mapping could be much lower than the given simple phantom and simulation results, and further validation is thus required. On the other hand, the coil structure can be obtained using a reference-based approach [42,43]. That is, using a phantom with pre-known dielectric properties, one can determine the coil information. The coil model can then be used for EPT studies. Secondly, Maxwell's curl equations are solved by the FDTD technique, which is usually time-consuming, and may take up to several hours to converge. In this study, to balance the spatial resolution and computational cost, an mm-level grid size for the FDTD solver was found to be effective to resolve the MRI-EPT problem. In this work, we have successfully employed GPU acceleration to enable an efficient implementation of the FDTD computing (in only tens of seconds). Through appropriately applying these techniques, it is then feasible to use a standard optimization technique to retrieve accurate EP distributions in practice.

Conclusions
In conclusion, this work proposed a novel EPT method that reconstructs tissue dielectric properties by directly solving Maxwell's curl equations. This method has several advantages over existing methods: (a) Without the use of second-order differential operations, the EPT solutions are more robust against B + 1 measurement errors and noise, (b) the GPU-enabled FDTD solver allows for the efficient implementation of an iterative optimization procedure for finding the EPT solutions, and the proposed method could therefore be used in clinical environments, and (c) since B + 1 phase information is not required in the solution procedure, the proposed method is more straightforward than many of the existing methods, particularly for the cases involving coil arrays in high-field MRI. These advantages have been validated through simulation and phantom-based experimental studies. The feasibility and effectiveness of the proposed EPT method will be further tested in human imaging at 7T, using both a birdcage coil and phased array coils.