Magnetic resonance-based reconstruction method of conductivity and permittivity distributions at the Larmor frequency

Magnetic resonance electrical property tomography is a recent medical imaging modality for visualizing the electrical tissue properties of the human body using radio-frequency magnetic fields. It uses the fact that in magnetic resonance imaging systems the eddy currents induced by the radio-frequency magnetic fields reflect the conductivity ($\sigma$) and permittivity ($\epsilon$) distributions inside the tissues through Maxwell's equations. The corresponding inverse problem consists of reconstructing the admittivity distribution ($\gamma=\sigma+i\omega\epsilon$) at the Larmor frequency ($\omega/2\pi=$128 MHz for a 3 tesla MRI machine) from the positive circularly polarized component of the magnetic field ${\bf H}=(H_x,H_y,H_z)$. Previous methods are usually based on an assumption of local homogeneity ($\nabla\gamma\approx 0$) which simplifies the governing equation. However, previous methods that include the assumption of homogeneity are prone to artifacts in the region where $\gamma$ varies. Hence, recent work has sought a reconstruction method that does not assume local-homogeneity. This paper presents a new magnetic resonance electrical property tomography reconstruction method which does not require any local homogeneity assumption on $\gamma$. We find that $\gamma$ is a solution of a semi-elliptic partial differential equation with its coefficients depending only on the measured data $H^+$, which enable us to compute a blurred version of $\gamma$. To improve the resolution of the reconstructed image, we developed a new optimization algorithm that minimizes the mismatch between the data and the model data as a highly nonlinear function of $\gamma$. Numerical simulations are presented to illustrate the potential of the proposed reconstruction method.

1. Introduction. Magnetic resonance imaging (MRI) system can visualize both the conductivity, σ, and permittivity, , of biological tissues at the Larmor frequency, which is approximately 128 MHz for a 3 tesla MRI machine. Magnetic resonance electrical property tomography (MREPT) uses a time-harmonic magnetic field inside an imaging object. The standard radio-frequency coil of the magnetic resonance scanner produces the field by feeding in a sinusoidal current at the Larmor frequency. The time-harmonic magnetic field, denoted by H = (H x , H y , H z ), reflects both the conductivity σ and permittivity of human tissues through the following arrangement of time-harmonic Maxwell's equations: in Ω, (1.1) occupying an imaging object. Here, we use the fact that the magnetic permeability of the human body is approximately equal to µ 0 . Clinical MRI scanners measure the positive rotating magnetic field, H + := (H x + iH y )/2, which is the component of the magnetic field H in the direction (1, i, 0)/2. This is because the MR signal, denoted by S, contains partial information about the time-harmonic magnetic field H = (H x , H y , H z ) in the following way S τ (r) ∝ M (r)H − (r)H + (r) sin(ατ |H + (r)|) |H + (r)| for r = (x, y, z) ∈ Ω, (1.2) where H − = (H x −iH y )/2 is the negative rotating magnetic field, M (r) is the standard MR magnitude image at position r, and α is a constant. Here, τ is the duration of the radio-frequency pulse that controls the intensity of the signal S τ . Acquiring two MR signals S τ1 and S τ2 with suitably chosen τ 1 and τ 2 , we can extract the H + data through (1.2) with the assumption that H + /|H + | ≈ H − /|H − |. This data acquisition technique is called B1 mapping, and was first suggested by Haacke et al . [10] in the early nineties. For details on the B1 mapping technique measuring H + , we refer to numerous published works in the literature [1,4,19,26,29]. The inverse problem of MREPT consists of reconstructing distributions of σ and from H + . To solve the inverse problem, we need to represent the distributions of σ and with respect to the data H + . Under the assumption of the local homogeneity, ∇(σ + iω ) = 0, the governing partial differential equation (1.1) directly gives the following simple relation between σ + iω and H + ; The most widely used MREPT reconstruction methods [30,12,13,14,15] are based on (1.3) as it gives the direct representation formula for σ + iω with respect to H + , in Ω. (1.4) However, when ∇(σ + iω ) is not small, the direct formula (1.4) produces serious reconstruction errors [21]. The local homogeneity assumption neglects the contribution of ∇ ln γ × (∇ × H) in (1.1). Such reconstruction errors are rigorously analyzed in [21]. We need to remove the local homogeneity assumption to develop a reconstruction method. Recently, a reconstruction method [22] removing the assumption of ( ∂ ∂x , ∂ ∂y )(σ + iω ) = 0 has been developed, although it still requires the assumption of ∂ ∂z (σ + iω ) = 0. The method is based on the finding that, under the assumption of longitudinal homogeneity, σ + iω is a solution of a semilinear elliptic PDE with coefficients that only depend on H + [22].
In this paper, with no assumption of local homogeneity for σ + iω , we develop a new reconstruction method. We find that σ and satisfy the elliptic partial differential equation, To improve the spatial resolution of the reconstructed image, we develop an optimal control method for the parameters σ and . In the proposed adjoint-based optimization method, the conductivity and permittivity distributions are updated iteratively by a nonlinear optimization algorithm which minimizes the discrepancy function describing the L 2 -mismatch between the forward model and the observed data. We compute the Fréchet derivatives of the discrepancy function with respect to σ and . This optimal control method requires a very good initial guess. Fortunately, we can obtain a good initial guess using (1.5). Several numerical simulations are carried out to show the validity of the proposed reconstruction method.
2. Governing equation for the admittivity reconstruction. We assume that an imaging object occupying a three-dimensional domain Ω with its boundary ∂Ω being of class C 2 . Let γ = σ + iω denote the admittivity of the subject at the MR Larmor frequency. For simplicity, we assume that γ is a constant near the boundary; that is, γ = γ 0 in the region Ω d := {x ∈ Ω | dist(x, ∂Ω) < d} for some d > 0, where γ 0 = σ 0 + iω 0 with σ 0 and 0 being known reference quantities.
Let H s (Ω) denote the standard Sobolev space of order s. We assume that the admittivity distribution γ = σ + iω belongs to the following admissible set A: The inverse problem is to invert the map γ → H + where H + represents the measured data extracted from the MR signal in (1.2) and the relation between H and γ is given in (1.1). Noting that the component H z is known to be relatively small with a regular birdcage coil of MRI scanner [27], we assume H z = 0.
To solve the inverse problem, we need to express σ + iω in terms of H + only using the governing equation (1.1). Taking the inner product of both sides of equation (1.1) with the vector a = (1, i, 0)/2, we have It follows from the result of [22] that the contribution of H − in (2.2) can be eliminated from the identity Equation (2.2) with the above identity gives the following lemma [22]. Lemma 2.1. The γ in (1.1) satisfies the following first-order partial differential equation

3)
where L is the linear differential operator given by According to Lemma 2.1, the inverse problem is reduced to solve the first-order partial differential equation (2.3) for γ. Unfortunately, it may not be possible to solve the first-order partial differential equation (2.3). As the direction vector field of LH + is not a real-valued function, the method of characteristics can not be applied. Indeed, Hörmander [11] and Lewy [17] provided non-existence results for the first order partial differential equation with complex-valued coefficients. To be precise, the governing equation (2.3) for H + can be rewritten in the standard form F · ∇u = f (·, u), where F = LH + , u = log γ, and f (·, u) = iωµ 0 e u H + − ∆H + . According to the Cauchy-Kowalevski theorem [16], the equation F · ∇u = f (·, u) with suitable initial data, can be locally solvable only when f is analytic. On the other hand, this local solvability can not be guaranteed for general f ∈ C ∞ from Lewy's example [17]. This is why we do not use the model (2.3) to compute γ.

Elliptic equation for the admittivity.
In this subsection, we prove that σ and satisfy the elliptic partial differential equation (1.5) which is one of our main results in this paper. This key observation follows from long and careful computations.
Theorem 2.2. The distributions of σ and satisfy the following equation: in Ω, (2.5) where A[H + ] is a positive semi-definite matrix given by (2.14) Proof. We first separate the governing equation (2.3) into its real and imaginary parts. Let H + r and H + i be the real and imaginary parts of H + , i.e., H + = H + r + iH + i . Let γ i be the imaginary part of the admittivity.
Equation (2.3) can be expressed as (2.15) The real and imaginary parts of equation (2.15) are given, respectively, by The real part (2.16) can be written as Similarly, the imaginary part (2.17) can be written as where F 1 is given in (2.8). Similarly, subtracting (2.20) from (2.23) yields with F 2 being given by (2.9). A direct computation shows that equation (2.24) can be expressed as Since P x = −Q y and P y = Q x , P x P y + Q x Q y = 0 and therefore, equation (2.26) can be rewritten as Similarly, (2.25) gives where det denotes the determinant. Hence, all the eigenvalues of the matrix A are non-negative.

Approximate solution.
Using the elliptic partial differential equation (2.5) in Theorem 2.2, we can compute a fairly good approximation of the true admittivity. Since the matrix A in (2.5) is degenerate, we need a regularization strategy. By adding a regularization term ρe T 3 e 3 to the matrix A, we can compute viscosity solution U ρ = (σ ρ , ω ρ ) T of the elliptic partial differential equation (2.5): with the Dirichlet boundary condition U 0 = (σ 0 , ω 0 ) T on ∂Ω, where e 3 = (0, 0, 1), superposed T denotes the transpose, and ρ is a small positive constant. Note that the matrix A + ρe T 3 e 3 is positive definite, since the eigenvalues of the matrix A + ρe T 3 e 3 are By solving equation (2.30), we can get the blurred admittivity image of the true distribution.
3. Adjoint-based optimization method. This section presents adjoint-based optimization method for finding admittivity distribution. A Newton iteration is used to find optimal solution, hence, a fairly good initial guess is required. The approximated solution in section 2.2 is used as an initial guess of the Newton iteration. Let H + m ∈ H 1 (Ω) be the measured data corresponding to the true admittivity γ * ∈ A; hence H + m satisfies For γ ∈ A, let H + [γ] be a solution of the Dirichlet problem: (3.1) The equation (3.1) has a unique solution for properly chosen c 1 in the definition of A in (2.1). From now on, we assume that c 1 is chosen so that (3.1) has a unique solution. Then, the map is well-defined. We define the misfit function J[γ] of the variable γ = σ + iω by the L 2 -norm of the difference between H + [γ] in (3.1) and the measured data H + m : In this minimization problem, we need to determine the Fréchet derivative of the misfit function J with respect to the control variable γ. Let A be defined by The following theorem proves the Fréchet differentiability of H + [γ] under the assumption that c 1 < 1. Theorem 3.1. Let c 1 < 1. For γ ∈ A, the map γ → H + is Fréchet differentiable. Let δ ∈ A be such that γ + δ ∈ A. The Fréchet derivative DH + [γ](δ) at δ is given by the solution u of the following equation in Ω, Proof. First, remember that for w ∈ H 1 (Ω), Then, defining in Ω, (3.6) Therefore, we have . (3.7) By Hölder's inequality and Sobolev embedding theorem (see (3.5)), the first term of the right-hand side of (3.7) can be estimated by Combining (3.7) and (3.8), we have . (3.9) By Hölder's inequality and Sobolev embedding theorem, (3.9) can be estimated by if c 1 < 1. Since the data difference w δ satisfies (3.6) and u is the solution of equation (3.4), the difference w δ − u ∈ H 1 (Ω) satisfies .

(3.12)
Again, by Hölder's inequality and Sobolev embedding theorem, the first term of the right-hand side of (3.12) can be estimated by Combining (3.12) and (3.13), we have . (3.14) 8 By Hölder's inequality and Sobolev embedding theorem, (3.14) can be estimated by if c 1 < 1. By inequalities (3.10) and (3.15), it follows that Thus, The following theorem expresses the Fréchet derivative of J[γ].

17)
where p is the solution of the adjoint problem: By (3.10), Thus, On integrating by parts, it follows that Moreover, Hence, Note that w δ satisfies the following identity: Therefore, which completes the proof.
It is worth mentioning that the smallness assumption on the bound c 1 defined in (2.1) ensures the well-posedness of (3.17) with homogeneous Dirichlet boundary condition.
In the next lemma, we rewrite the adjoint problem (3.18) as a second-order elliptic partial differential equation. Lemma 3.3. For γ = σ + iω , the adjoint problem (3.18) can be rewritten as with the Dirichlet boundary condition p = 0 on ∂Ω, where Proof. Denote by v := log γ. Since the linear operator L is given by (3.23) Hence, if we substitute (3.23) into the adjoint problem (3.18), we get the second-order elliptic partial differential equation (3.22) and the proof is complete.
4. Newton-type reconstruction algorithm. In the section 3, the Fréchet differentiability of the discrepancy functional J is proven. To find γ ∈ A such that J[γ] = 0, we apply the Newton method. The Newton method starts from the linearization of the functional J: where h ∈ A and γ, γ + h ∈ A. Let γ n be the n-th iteration. Given γ n , Newton's method seeks to find h n such that If we update the iteration as γ n+1 = γ n + h n , J[γ n+1 ] ≈ 0 by (4.1) and (4.2). Note that DJ[γ](kh) = kDJ[γ](h) for any real number k by (3.17). The following lemma shows how to find h n . Hence, for any f n ∈ A such that DJ[γ n ](f n ) = 0, h n = − J[γn] DJ[γn](fn) f n is the solution of (4.2).
Lemma 4.1 proves that we can make the step size h n of (4.2) if we choose the function f n ∈ A such that DJ[γ n ](f n ) = 0. Equation (3.17) shows that DJ[γ n ](f n ) can be represented by L 2 inner product, DJ[γ n ](f n ) = < f n , g n >, where g n := 1 γn ∇ · (p n LH + [γ n ]) + iωµ 0 H + [γ n ]p n . Note that g n is computed from given γ n and the solution of the adjoint problem (3.18) p n . If we choose f n = g n , then DJ[γ n ](f n ) = ||g n || 2 2 = ||g n || 2 2 . In that case, DJ[γ n ](f n ) = 0 unless g n = 0 and the step size h n becomes So, the Newton iteration algorithm is given by where g n = 1 γn ∇ · (p n LH + [γ n ]) + iωµ 0 H + [γ n ]p n . It is worth emphasizing that the Newton method guarantees the convergence only when the initial guess γ 0 is close enough to the true solution.
To find a good initial guess for γ, we use an iteration scheme to solve (2.30) with small regularization parameter ρ: .

(4.4)
Based on the above iteration scheme, we develop the following reconstruction algorithm.
Step 4. Compute H + [γ n ] from the given γ n , for n ≥ 0. From t(3.1), the following equation for H + can be obtained: with H + = H + m on ∂Ω and G[γ n ] being the vector field in (3.22).

Numerical simulations.
In this section, we will present numerical simulation results from two models to validate the proposed algorithm. In the first model, we set the domain Ω to be a cylindrical model where the admittivity distribution does not change along the z-direction.   In subsection 2.1, we proved that the solution of (2.5) is the blurred approximation of true admittivity. However, (2.5) is degenerate since the diffusion matrix A[H + ] is singular. So, we modified (2.5) to (4.4) by adding the regularization term ρe T 3 e 3 . Figure 5.3 shows the determinant of A and A + ρe T 3 e 3 , where ρ is 5% of the maximum of A 33 , P 2 z + Q 2 z in (2.5). Figure 5.3 explains that the regularized semi-elliptic PDE is also degenerate near l = {(0, 0, z) | z ∈ R}. To avoid this, we segmented subdomain D near l, as shown in Figure 5.4. In the subdomain Ω\D, we applied the iteration method (4.4). Figure 5.4 illustrates solutions of (4.4), U k = (σ k , ω k ), in Ω\D with various iteration numbers k. We set the initial values to be constant: σ 0 = 1 and ω 0 = 0. In order to check the convergence and the accuracy of the proposed algorithm (4.4), we plotted ||γ k − γ * || 2 and ||γ k − γ k−1 || 2 with k = 1, 2, · · · , 10 in Figure 5.5. Figure 5.5 shows that the iteration method (4.4) converges as k increases and the error between true admittivity and reconstructed admittivity decreases. We choose U 3 to be the solution  of the iterative algorithm. We used direct method (1.4) for the admittivity value γ in the segmented subdomain D. So, we let U 3 with the value obtained from the direct method in D to be the initial guess of the Newton method (4.7). Figure 5.6 illustrates the reconstructed conductivity and relative permittivity distribution by the Newton iteration (4.7).  images by the direct formula (1.4) and the proposed method. Figure 5.10 compares between the two methods, the direct formula (1.4) and the proposed method, for imaging small anomalies.
In the second numerical model, we simulate the model with admittivity changing along z-direction, i.e., ∂γ ∂z = 0. The domain Ω is decomposed into two parts, Ω − = Ω ∩ {z < 0} and Ω + = Ω ∩ {z ≥ 0}. In Ω − , the admittivity distribution is the same as in Model 1. However, the admittivity distribution in Ω + is different from Model 1 and is such that ∂γ ∂z = 0 in Ω 0 . Figure 5.11 shows the second configuration model, the conductivity and the relative permittivity values in the domain. Figure 5.12 shows the reconstruction results using the direct method and the proposed method. Figure  5.13 shows the accuracy of the proposed method applied to the second model. Figure  5.14 presents the reconstructed conductivity distribution of the second model in the slice of Ω − using the proposed method. Figure 5.13 and Figure 5.14 demonstrate that the proposed method works well in the case of ∂γ ∂z = 0.
6. Concluding remarks. In this paper, we have developed an iterative novel scheme for reconstructing electrical tissue properties at the Larmor frequency from measurements of the positive rotating magnetic field. We first suggest the elliptic partial differential equation (2.30) which provides a blurred reconstructed image. By considering the blurred reconstructed image as an initial guess of the Newton iteration, the Newton iteration for finding the minimizer of the functional J in (3.3) finds the final reconstruction admittivity. Note that our scheme does not require a local homogeneity assumptions on γ and allows to reconstruct inhomogeneous distributions accurately.