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

Abstract

Magnetic Resonance-Electrical Properties Tomography (MR-EPT) is a method to reconstruct the electrical properties (EPs) of bio-tissues from the measured radiofrequency (RF) field in Magnetic Resonance Imaging (MRI). Current MR-EPT approaches reconstruct the EP profile by solving a second-order partial differential wave equation problem, which is sensitive to noise and can induce large reconstruction artefacts near tissue boundaries and areas with inaccurate field 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 fitting the RF field calculated by the finite-difference-time-domain (FDTD) technique to the measured values. To solve the time-consuming problem of the iterative fitting, a graphics processing unit (GPU) is used to accelerate the FDTD technique to process the field 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 field 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 field, the proposed approach is still capable of offering robust EPT solutions, thus having improved practicality for potential clinical applications.

1. 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 shown promise for clinical diagnosis and guiding the treatment of malignant tissues [2,3] and improving the patient-specific electromagnetic (EM) modelling of hyperthermia [4]. The concept of MR-EPT was introduced by Haacke et al. in 1991 [5], and Wen et al. [6] verified this concept through experiments in 2003. The initial development of MR-EPT was based on the utilization of the Helmholtz equation, in which the permittivity and conductivity were assumed constant (also known as the “local homogeneity assumption” (LHA)) [7,8]. The methods based on this LHA have been found to be sensitive to noise, and severe distortions at the boundaries of different tissues may occur [9,10]. To reduce the noise amplification in the utilization of second-order derivatives, several methods have been proposed. Additional regularization terms were added to solve the EPT inverse problem [11,12], the Laplacian operator was replaced with a surface integral [8,13], the k-space data were weighted properly [14], and the region of interest (ROI) was expanded with more voxels [15]. Numerical and experimental results in References [11,12,13,14,15] demonstrated that these methods improved the noise robustness but failed to correct the edge artefacts caused by the LHA. To eliminate these edge artefacts, a new set of partial differential equations (PDE) was introduced [16,17,18,19,20,21]. In References [16,17,18], the special gradients of EPs were considered as unknown variables in the PDE, which were resolved by using multi-channel parallel transmit coils. In this scheme, the condition number of the constructed PDE was reduced, and the inverse solution was thus improved. An alternative method to resolve the unknown EP variables is to transform the PDE into a linearly regularized system and search for its global solution [19,20,21]. Although the methods in References [16,17,18,19,20,21] have been tested in experiments, it is noteworthy that the EP reconstructions were implemented in an inverse process, in which the solution robustness could be weakened by the ill-conditioned property of the linear system. For these reasons, the development of new EPT approaches without the involvement of second-order PDE (i.e., Laplacian-operator free) is of interest.

Recently, EP reconstruction in MRI has been alternatively investigated using the integral equations (IE). The IE-based EPT methods rely on a global convolution between the Green function and the contrast source [22,23,24,25,26] With the integral of global convolution, the noise robustness is theoretically improved in those IE-based methods. However, the development of IE-EPT is still at an initial stage and facing several practical problems. The technical problems include the requirement of the estimation of the incident electromagnetic fields, low spatial resolution, the accuracy issue associated with the use of the free-space Green function to construct the contrast source and the necessity of using global convolution.

In this paper, a novel MR-EPT method is proposed. In the presented method, the EP reconstruction was achieved by directly solving Maxwell’s equations [27,28]. That is, the first-order PDE was numerically solved by using an in-house finite-difference-time-domain (FDTD) method with the graphics processing unit (GPU)-based acceleration. The unknown EP profile was involved in the FDTD kernel solver, and an optimization procedure was performed to refine the EP profile to minimize the difference between the calculated and measured transmit magnetic fields ($B_1^{+}$). To demonstrate the feasibility of the proposed EPT method, human brain simulation and phantom-based experimental studies have been carried out, and EP reconstruction results were compared with conventional approaches.

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.

2.1. 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 $\left|B_1^{+}\right|=\left|B_x+i{B}_y\right|/2$ [31]. Here, as Maxwell’s curl equations state, the calculated $\left|B_1^{+}\right|$ 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 $\left|B_1^{+}\right|$ 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].

2.2. Optimization

EPs were estimated through an iterative optimization procedure. The optimization variables (EPs) are ${\mathit{x}}_{\mathrm{EP}}=\left[{\mathit{\epsilon }}_{\mathrm{tissue}}{\mathit{\sigma }}_{\mathrm{tissue}}\right]$, 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 ${\mathit{x}}_{\mathrm{EP}}$ were mapped to the $B_1^{+}$ field as:

$$B_1^{+}={ℱ}_{\mathrm{FDTD}}\left({\mathit{x}}_{\mathrm{EP}}\right),$$

(1)

where ${ℱ}_{\mathrm{FDTD}}$ denotes the FDTD process.

Then, a cost function which describes the residual between the mapped $B_1^{+}$ fields and the measured $B_1^{+}$ field can be built as:

$$C={\Vert \left|B_1^{+}\right|-\left|B_{1;meas}^{+}\right|\Vert }_2,$$

(2)

Therefore, an unconstrained optimization problem which aims to minimize the cost function can be formed as:

$$\underset{\mathit{x}}{\min }{\Vert \left|{ℱ}_{\mathrm{FDTD}}\left({\mathit{x}}_{\mathrm{EP}}\right)\right|-\left|B_{1;meas}^{+}\right|\Vert }_2,$$

(3)

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 ${\mathit{x}}_{\mathrm{EP}}$ is found from Equation (3).

In this work, the genetic algorithm (GA, standard MATLAB function) was used to find the optimal solution of ${\mathit{x}}_{\mathrm{EP}}$. In the optimization procedure, an initial EP distribution denoted as ${\mathit{x}}_{\mathbf{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 ${\mathit{x}}_{\mathbf{O}}\pm 0.3{\overline{\mathit{x}}}_{\mathbf{O}}$, where ${\overline{\mathit{x}}}_{\mathbf{O}}$ is the mean value of ${\mathit{x}}_{\mathbf{O}}$. Alternatively, the initial EPs of each tissue can be set based on the data from Reference [32], and the variable range was set as ${\mathit{x}}_{\mathbf{O}}\pm 0.3{\mathit{x}}_{\mathbf{O}}$, where ${\mathit{x}}_{\mathbf{O}}$ are the reported values in Reference [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_{\mathrm{desired}}$, and (2) the number of iterations exceeds certain pre-set limits, $N_{\mathrm{desired}}$. Once condition (1) or (2) was satisfied, the optimization process could be terminated. The parameter setup of $C_{\mathrm{desired}}$ and $N_{\mathrm{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_{\mathrm{desired}}$ and $N_{\mathrm{desired}}$ as ${10}^{-3}$ and 100, respectively.

The iterative optimization procedure is implemented in MATLAB (R2015a, 64 bits, The MathWorks, Inc., Natick, MA, USA).

2.3. 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³ was re-meshed as 2 mm³, 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].

2.4. 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 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_{\mathrm{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.

3. Results

3.1. 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. Reconstructed EPs of the phantom (based on simulated data).

SNR	Inner		Outer
	εr	σ(S/m)	εr	σ
Target	78.5	0.43	78.5	0.78
Noise-free	78	0.4309	78	0.7784
SNR = 15	77	0.4428	79	0.7644
SNR = 25	77	0.4212	78	0.7918
SNR = 35	78	0.4278	78	0.7843
SNR = 45	78	0.4326	78	0.7762
SNR = 55	78	0.4322	78	0.7824

SNR: signal-to-noise ratio. σ in S/m.

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

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.

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

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 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. Reconstructed EPs of the phantom (based on simulated data).

Compartment		Probe Measurement	Reconstructed
Inner	εr	78.5	78
	σ	0.78	0.7768
Outer	εr	78.5	78
	σ	0.43	0.4324

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 $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 $\left|B_z\right|$ was quite small. Close to the coil ends, the $\left|B_z\right|$ value was elevated. The availability of these magnetic field components should be helpful to improve the EPT solution and its applications in MRI.

4. Discussion

4.1. 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 $\left|B_1^{+}\right|$ 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 $\left|B_z\right|/\left|B_{x,y}\right|$ 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.

4.2. 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 $\left|B_1^{+}\right|$ 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.

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