Optimization of magnetic flux density for fast MREIT conductivity imaging using multi-echo interleaved partial fourier acquisitions

Background Magnetic resonance electrical impedance tomography (MREIT) has been introduced as a non-invasive method for visualizing the internal conductivity and/or current density of an electrically conductive object by externally injected currents. The injected current through a pair of surface electrodes induces a magnetic flux density distribution inside the imaging object, which results in additional magnetic flux density. To measure the magnetic flux density signal in MREIT, the phase difference approach in an interleaved encoding scheme cancels out the systematic artifacts accumulated in phase signals and also reduces the random noise effect by doubling the measured magnetic flux density signal. For practical applications of in vivo MREIT, it is essential to reduce the scan duration maintaining spatial-resolution and sufficient contrast. In this paper, we optimize the magnetic flux density by using a fast gradient multi-echo MR pulse sequence. To recover the one component of magnetic flux density Bz, we use a coupled partial Fourier acquisitions in the interleaved sense. Methods To prove the proposed algorithm, we performed numerical simulations using a two-dimensional finite-element model. For a real experiment, we designed a phantom filled with a calibrated saline solution and located a rubber balloon inside the phantom. The rubber balloon was inflated by injecting the same saline solution during the MREIT imaging. We used the multi-echo fast low angle shot (FLASH) MR pulse sequence for MRI scan, which allows the reduction of measuring time without a substantial loss in image quality. Results Under the assumption of a priori phase artifact map from a reference scan, we rigorously investigated the convergence ratio of the proposed method, which was closely related with the number of measured phase encode set and the frequency range of the background field inhomogeneity. In the phantom experiment with a partial Fourier acquisition, the total scan time was less than 6 seconds to measure the magnetic flux density Bz data with 128×128 spacial matrix size, where it required 10.24 seconds to fill the complete k-space region. Conclusion Numerical simulation and experimental results demonstrated that the proposed method reduces the scanning time and provides the recovered Bz data comparable to what we obtained by measuring complete k-space data.

http://www.biomedical-engineering-online.com/content/12/1/82 Background Magnetic resonance electrical impedance tomography (MREIT) utilizes a magnetic resonance imaging (MRI) scanner to measure magnetic flux density B z data inside an imaging object induced by the externally injected current. The internal current density distribution has been studied in magnetic resonance current density imaging (MRCDI) by measuring the whole magnetic flux density data B = (B x , B y , B z ) [1,2]. Combining MRCDI and electrical impedance tomography (EIT) technique, MREIT provides the cross-sectional conductivity images of the object with high spatial resolution [3][4][5][6][7][8][9][10]. Since an MRI scanner measures only one component B z of B without rotating the imaging object, most MREIT algorithms assumed that the internal conductivity is isotropic and focused on visualizing its distribution by using one component of the magnetic flux density data B z of B [11][12][13][14][15][16][17][18].
Recent MREIT imaging techniques have been developed with respect to both the capacity of measurement techniques and the numerical reconstruction algorithms. Experimental results from in vivo animal and human have been reported [19,20] in MREIT. As an innovation of current MREIT, a fast MREIT imaging technique referring to the continuous monitoring of objects includes various wide application areas [21]. Recently, one of challenging problem in MREIT is to implement a new imaging technique with a very short acquisition time for the imaging of neural activities of brain related to the conductivity change.
Since current MREIT experiments suffer from poor SNR of the measured B z data under the typical data acquisition durations and a small amount of injected current, it is important to reduce the scan time, while maintaining the spatial-resolution and sufficient contrast, for practical implementations of in vivo MREIT. Recently, to reduce the scan time in MREIT, Hamamura et al [22] reconstructed the interior conductivity using a single-shot spin-echo echo planar imaging (SS-SEPI) pulse sequence and Muftuler et al [23] used a SENSE-accelerated imaging technique to acquire phase signal by the injected current. One of basic approaches for maintaining the spatial resolution is to reduce the number of phase encoding steps because each phase encoding step requires a certain amount of time for execution. Since the MREIT techniques use an interleaved phase encoding acquisition scheme to double the B z signal, Park et al [24] reconstructed the phase signal B z by filling the skipped k-space region using the interleaved measurement property.
To obtain the static conductivity image in MREIT [19,20], a spin-echo based MREIT pulse sequence has been predominantly used to reduce the background artifact and to increase the imaging quality. In real situations, it is difficult to employ the fast conventional MR pulse sequences because the noise standard deviation of B z is inversely proportional to the width of injection current and the intensity of MR magnitude, simultaneously. An MREIT pulse sequence should be devised to enhance changes in MR phase images for given current amplitudes.
In this paper, we used the multi-echo fast low angle shot (FLASH) MR pulse sequence which allows the reduction of imaging time without any substantial loss in image quality. In addition, the multi-echo FLASH sequence maximizes the width of injection current extending the duration of injection current until the end of a readout gradient in MREIT. To reconstruct the internal conductivity distribution, most algorithms require at least two independent injection currents in an interleaved sense, which require relatively a http://www.biomedical-engineering-online.com/content/12/1/82 long scanning duration [7,11]. To reduce the scanning time, we adopt a partial phase encoding acquisition scheme using the multi-echo FLASH MREIT pulse sequence and rigorously investigate the relationship between the convergence ratio of the algorithm and the background field inhomogeneity [24]. We consider the discrete 2 -norm to evaluate the convergence ratio.
To show the feasibility of the proposed algorithm, we performed numerical simulations and compared the performance to the simulated true B z data. We designed a cylindrical acrylic phantom filled with a calibrated saline solution and located a rubber balloon inside the phantom. The rubber balloon was inflated by injecting the same saline solution during the scan. The phantom was designed to provide a homogeneous magnitude image, but distinguishable signals of measured B z between the inside and outside the balloon. The phantom experiment demonstrated that the proposed method reduces the scanning time and recovers the reasonable resolution of B z , which is comparable to the recovered B z using the complete k-space data.

k-space signal and B z data
In a conventional spin echo MREIT pulse sequence, both positive and negative currents of the same amplitude and duration are injected with reverse polarity. These injection currents with the pulse width of T c accumulate extra phases. Corresponding k-space MR signals can be described as where ρ is the T 2 weighted spin density, δ is any systematic phase artifact, and is a fieldof-view (FOV). Here, the superscript of S ± (k x , k y ) denotes a brief notation for S + (k x , k y ) and S − (k x , k y ). For the standard coverage of k-space, we set where γ = 26.75 × 10 7 rad/T · s is the gyromagnetic ratio of hydrogen, t is the time between samplings, G x is the frequency encoding gradient strength, T E is the echo time, G y is the phase encoding step, and T pe is the phase encoding time. The induced magnetic flux density ±B z is generated by the positive and negative injection currents I ± . Applying the inverse Fourier transform to the measured k-space data sets in (1), we can compute the magnetic flux density B z as where α and β are the imaginary and real part of ρe iδ e iγ B z T c / ρe iδ e −iγ B z T c , respectively [2]. The current MREIT method is based on the electromagnetic information embedded in the measured B z data in order to visualize the conductivity(or current density) based on Biot-Savart law: where μ 0 = 4π10 −7 Tm/A is the magnetic permeability of the free space. The current density J in is given for the isotropic conductivity σ in J(r) = −σ ∇u(r) (5) and satisfies the following elliptic equation where ν is the outward unit normal vector on ∂ and g is an applied current density on the surface.

Recovery of complex T * 2 weighted spin density using interleaved partial Fourier acquisition
To simplify, we develop a theory using a conventional cartesian k-space that has three zones within the k y (phase-encode) domain; the central region (P 0 ), the positive region (P + ) and the negative region (P − ): Here, N c denotes the number of partial Fourier over-sampling phase-encodes. The interleaved k-space data S + p (k x , k y ) and S − p (k x , k y ) as partially acquired for the phase k y = γ 2π m G y T pe can be expressed where S ± p represents S + p and S − p simultaneously. We propose an algorithm to determine the T * 2 weighted spin density ρ : using the partially scanned k-space data S ± p . The systematic phase artifact e iδ(x,y) , unavoidable artifacts due to the main field inhomogeneity and the mismatch between the center of data acquisition interval and echo formation, arises from a low frequency field, which mainly belongs to the central region P 0 . Including a small perturbed phase artifact e iδ(x,y) , we start with the initial guess S ± p and design an alternating procedure by updating the skipped k-space regions.
The recovered ρ ± p = FT −1 (S ± p ) can be formally expressed by the equation where FT −1 P − (S ± p (k x , k y )) denotes the inverse Fourier transform by zero-filling the k-space except in the region P − and the remainder term E ± p is Let us define a support region D δ for a low spatially varying magnetic field due to background field inhomogeneities: The proof of observation 1 is provided in the Appendix A. By using the observation 1, The proof of observation 2 is provided in the Appendix B. The observation 2 shows that the skipped region in the measured k-space S + can be recovered by using the interleaved acquired S − and the estimated background sensitivity map.

Convergence characteristics
When the common measured P 0 does not cover the support region D δ , N δ > N c , for the phase-encode (k x , k y ) ∈ P + , the Fourier transform of ρ − p e −2iδ can be written by following the observation 1: From the relation (14), for the phase-encode (k x , k y ) ∈ P + , we have We set where R is a subregion of the k-space and FT(ψ)| R denotes the restriction of FT(ψ) to the region R. The discrete 2 where the constant C is independent to the function ψ and the region R. http://www.biomedical-engineering-online.com/content/12/1/82 When the skipped k-space region, P − , are filled the previously updated as in (13), we have For the k-space region (k x , k y ) ∈ P + , we have the following identity where Since the updated complex density Thus, the discrete 2 -norm of the difference between the true and the iteratively updated T * 2 weighted spin density can be estimated Detailed estimates of 2 -norm calculation are presented in the Appendix C.
Define an estimator for the convergence of the proposed algorithm Using the same procedures of (17) and (19), we have The relations (17), (19) and (21) show that the convergence of the proposed method depends on the estimator Z 2δ,P ± :

Optimization of B z using gradient multi-echo data
Since the noise standard deviation s B z of the measured B z is inversely proportional to injection current duration T c and the SNR of MR magnitude image ϒ M [25,26] as To reduce the noise level of B z and the imaging time, we applied the proposed method to the gradient multi-echo pulse sequence as a fast MR imaging technique. Subsequently, http://www.biomedical-engineering-online.com/content/12/1/82 by using the gradient multi-echo, it is possible to inject the current for a long duration to maximize the off-resonance phase.
When T E 1 denotes the first echo time and T E is the echo spacing, the m-th echo time is where N E is the echo number. The phase artifact δ T Em depends on the echo time T E m . Using a priori estimation of the phase artifact map δ T Em from a reference scan, for a time varying functional MREIT technique, the measured k-space region including P 0 and P + can be determined by taking account of the imaging multi-echo times T E m , the echo number N E and the repetition time T R .
The recovered multiple T * 2 -weighted complex densities ρ ± T Em , m = 1, · · · , N E , using the proposed algorithm in each echo time T E m can be optimized to generate a representative measured B m z data where α m and β m are the imaginary and real parts of ρ + T Em /ρ − T Em , respectively. The recovered multiple B m z data include pixel-by-pixel different noise level depending on the different imaging time T E m and the width of injection current. The multiple measured B m z data are optimally combined to reduce the noise level of B z [27]: where the point-wise weighting factor ω m (r) is given as Now, we setup an algorithm to reconstruct the magnetic flux density B z data using the partially measured k-space data S ± p and involving the following steps: 1. Take an initial guess ρ ± 0 = FT −1 (S ± | P 0 ∪P + ). 2. Transform the n-th updated ρ ± n e −2iδ to k -space by taking the Fourier transform. 3. Update the measured k -space data S ± n+1 by filling the skipped region in S ± p with the transformed FT(ρ ± n e −2iδ ) data. 4. Update ρ ± n+1 by taking the two-dimensional inverse Fourier transform.

Numerical simulation setup
To validate the proposed algorithm, we performed numerical simulations with the twodimensional finite-element model of a object 20 × 20 cm 2 with 256 × 256 rectangular elements and with the origin at its bottom-left, shown in Figure 1. We added the different complex field inhomogeneity artifacts to the simulated spin density image in Figure 1(a). The target magnetic flux density B z in Figure 1(b) was generated by solving the elliptic equation (6) and by using the Biot-Savart law given by (4). Since the magnetic flux density http://www.biomedical-engineering-online.com/content/12/1/82 B z in Figure 1(b) is continuous and has no abrupt changes, we used |∇B z | image displayed in Figure 1(c) to enhance image for the magnetic flux density. Figure 1(d) shows a partially measured k-space data corresponding to ρ + .
The target conductivity distribution σ had different anomalies with different conductivity values and the amount of injection current was 10 mA. Set an applied current density g on the surface as

Phantom imaging experimental setup
For the practical application of the proposed method as a fast MREIT imaging, we designed a cylindrical phantom filled with the saline solution of conductivity 1 S/m (shown in Figure 2(a)), including a rubber balloon for the visualization of isotropic conductivity excluding other artifacts by any concentration gradient in the phantom. The inside of balloon was filled with the same saline solution and the volume of balloon was controlled by injecting saline solution during the imaging experiment. After positioning the phantom inside the bore of 3T MR scanner (Achieva TX, Philips Medical Systems, Best, The Netherlands) with 8 channel RF coil, we collected k-space data using the gradient multi-echo injection current nonlinear encoding (ICNE) pulse sequence which was originated from FALSH (Figure 2(b)). To obtain the MR magnitude and magnetic flux density (B z ) images, it extends throughout the duration of injection current until the end of a readout gradient [28]. Since the multi-echo ICNE pulse sequence was synchronized to the injection currents with alternating polarity, it enabled to maximize the width of the injection currents and minimize the noise standard deviation of the measured B z data. The maximum amplitude of injection current was 5 mA and the total imaging time was 10.24 second to fill the k-space in the interleaved sense. Since the total imaging time was corresponding to the whole k-space scan, the actual imaging time would be reduced to 5.52 and 5.92 seconds for the number of partial region N (P 0 ∪ P + ) = 69 and 74, respectively. The imaging parameters were followings: slice thickness 5 mm, number of imaging slices one, repetition time T R = 40 ms, echo spacing T E = 6 ms, flip angle 40 degree, and multi-echo time T E m = 6 + (m − 1) × 6 ms for N E = 4. The FOV was 160 ×160 mm 2 with a matrix size of 128 × 128. The duration of current injection T c m was almost same to the multi-echo time T E m = 6 + (m − 1) × 6, m = 1, 2, 3, 4, because the current was continuously injected until the end of the readout gradient. Figure 3(a) shows the multiply acquired magnitude images |ρ + T Em |, m = 1, · · · , 4, where ρ + T Em was the m-th measured T * 2 weighted complex spin density, Figure 3

Simulation results
We compared the reconstructed magnetic flux density B z using the complete k-space data to the B z achieved using the partial k-space data. To evaluate the convergence characteristics of the proposed algorithm, we define the relative 2 -errors: where ψ and ψ n are the recovered images with the complete k-space data and the n-th updated data using the partially measured k-space data, respectively, and · denotes the 2 -norm.
To investigate the estimator Z 2δ,P ± for the convergence of the proposed iterative algorithm to fill the sipped k-space region P − , we fixed the number of partially measured k-space region as N (P 0 ∪ P + ) = 138, i.e., N (P 0 ) = 20 and N (P + ) = 118, and changed the frequency range of background field inhomogeneity. We generated several background field inhomogeneities changing the phase frequency range in k-space region by taking N δ = 5, 10, 20, 30 where the background field inhomogeneity δ satisfies FT(e iδ )(k x , k y ) = 0 for |k y | > N δ . Figure 4(a)-(d) show the background field inhomogeneities used in the reconstruction procedure and Table 1 shows the estimated Z 2δ,P ± in each background field inhomogeneity for N δ = 5, 10, 20, and 30, respectively. Figure 5(a) shows the reconstructed |∇B z | images using the fixed background field inhomogeneity with N δ = 5. From the top-left to the bottom-right, each image was corresponding to the j-th iterative updated |∇B z |, j = 0, 1, · · · , 10. The recovered B z using the partially measured k-space data included a large amount of artifacts (the top-left image in Figure 5(a)). However, since the value of Z 2δ,P ± was small, the first updated magnetic flux density almost recovered the true B z . Figure 5(b) shows reconstructed |∇B z | images using the fixed background field inhomogeneity with N δ = 30 corresponding to Figure 5(a). Since the value of Z 2δ,P ± was 0.7738, the convergence ratio of ρ ± was relatively slow comparing to the field inhomogeneity with N δ = 5. Table 2 shows the relative discrete 2 -errors of updated complex spin density for each iteration number and the convergence ratio of ρ ± was depending on the number of N δ . Table 3 shows the relative 2 -errors of reconstructed ∇B z for each iteration number depending on the number of N δ . The decay rates of the relative 2 -error were very fast as the number of N δ was small, but we needed relatively many iterations to approach the  Table 1 Calculated estimator Z 2δ,P ± in (20) for N δ = 5, 10, 20, 30 and fixed measured k-space data N (P 0 ∪ P + ) = 138 Z 2δ,P ± 0.0013 0.0749 0.4927 0.7738 required accuracy as the number of N δ was increase, even though the update procedure was rapidly computed by use of the fast Fourier transform.

Phantom experimental results
For the phantom experiment, we changed N c = 5, · · · , 10 for the set P 0 to investigate the convergence behavior with respect to a given background field inhomogeneity.
Using the collected k-space data with 8 channel RF coil and the gradient multi-echo by alternating readout gradient, we measured the T * 2 weighted complex densities ρ ±n m , n = 1, · · · , N CH , m = 1, · · · , N E , where N CH = 8 denotes the coil number and N E = 4 is the echo number. Figure 6 shows the measured background field inhomogeneities by displaying the real part of e 2iδ mn corresponding to the n-th coil and the m-th echo image. According to the increase of echo number, the accumulated background field inhomogeneity also increased. Figure 5 Reconstructed |∇B z | images using the fixed background field inhomogeneity with N δ = 5 and N δ = 30. a) reconstructed |∇B z | images using the fixed background field inhomogeneity with N δ = 5. b) reconstructed |∇B z | images using the fixed background field inhomogeneity with N δ = 30. From top-left to bottom-right, each image corresponds to the j-th iterative updated |∇B z |, j = 0, 1, · · · , 9. http://www.biomedical-engineering-online.com/content/12/1/82  Figure 7(a) and (c) show the measured T * 2 weighted magnitude and magnetic flux density B z images at the time T E m , m = 1, 2, 3, 4, using partially acquired k-space region P 0 ∪ P + with N c = 5 by a transversally injected current. Although the amount of accumulated phase signal by the injected current increased as the echo time varied from T E 1 to T E 4 , the magnitude image at the 4-th echo was more deteriorated comparing to the 1st echo case. Figure 7 Comparing to the measured images in Figure 7(a)-(d), in contrast to the recovery of low phase frequency information corresponding to P 0 , the increased background field inhomogeneity caused relatively high frequency artifacts. Figure 8(a)-(d) shows iteratively updated T * 2 weighted magnitude and magnetic flux density B z images using P 0 ∪P + with N c = 5. We fixed the update iteration number as 20 for all experiments. When we fixed N c = 5, the 1-st and 2-nd recovered T * 2 weighted complex densities in Figure 8(a) and (c) were relatively close to the recovered ones using the complete k-space data. However, as the phase artifact increased, the 3-rd and 4-th recovered T * 2 weighted complex densities were deficient in reflecting full information of B z signal. Especially, the 4-th updated magnetic flux density B z image shows some defective region due to the insufficient recovery of T * 2 weighted complex density. Figure 8(e)-(f ) shows iteratively updated T * 2 weighted magnitude and magnetic flux density B z images corresponding to Figure 8(a)-(d) using P 0 ∪P + with N c = 10. When we used N c = 10, the updated T * 2 weighted complex densities almost recovered the magnetic flux density B z data comparing to those using the complete k-space data.

Discussion
We used the gradient multi-echo MREIT pulse sequence to reduce the imaging time and to maximize injection current duration. Since the MREIT techniques utilize accumulated phase signal by the injected current, it requires enough repetition time T R to accumulate the phase signal. In this sense, the gradient multi-echo MREIT pulse sequence seems practical approach for the improvement of B z quality as well as reducing the imaging time.  In this paper, we used a partially acquired k-space data in the phantom experiment by filling the k-space as much as 74 line by line, results in 5.92 second to image the resolution of 128 × 128. Experimental results show that the proposed interleaved partial Fourier strategy for MREIT has a potential to reduce scan times and maintain the information of B z data comparable to what is obtained with complete k-space data.
The convergence ratio of the iteratively updated phase signal heavily depends on the frequency of the background filed inhomogeneity and the number of half-Fourier oversampling phase-encodes P 0 . Instead of the gradient multi-echo, if we use the spin multiecho pulse sequence, the proposed iterative algorithm would rapidly recover T 2 -weighted complex spin density due to a small amount of background field inhomogeneity. However, in spite of some advantages of the spin multi-echo MREIT pulse sequence, for a realtime MREIT imaging, MR pulse sequence should be carefully investigate by taking into account of the width of injection current, the scan duration and the low SNR of measured B z signal.
In this paper, we assumed a priori background field inhomogeneity which is typically used in the sensitivity encoding (SENSE) as a fast MRI measurement technique. Since the MREIT techniques typically used interleaved acquisition by injecting alternative currents, http://www.biomedical-engineering-online.com/content/12/1/82 it may be possible to extract background field inhomogeneity information under a low frequency range assumption and by cancelation of B z information: ρ + (x, y)ρ − (x, y) = ρ 2 (x, y)e 2iδ(x,y) Several studies reported for the feasibility of MREIT to detect neural activities in the brain, directly [29,30]. Functional MREIT technique is suggested to image brain activity via conductivity change related to neural activity through the fast MREIT pulse sequence. Our future study will focus on applying the proposed method to produce functional conductivity images of animal and/or human brain to pursue rapidly changing conductivity associated with neural activities.

Conclusion
In MREIT, the inherent challenges are to reduce the scan time and maintain current injection duration to make it feasible for the clinical applications. We developed an iterative method to optimize the measured magnetic flux density B z using the multi-echo interleaved partial Fourier acquisitions for fast imaging in MREIT. The proposed method used a fast gradient multi-echo MR pulse sequence to reduce the scan time and to maximize the phase signal by injection current. Under the assumption of a priori background field inhomogeneity map, we rigorously investigated the convergence ratio of the proposed method using the discrete 2 -norm, which was closely related with the number of measured phase encode set and the frequency range of the background field inhomogeneity. To evaluate the proposed method, a specially designed conductivity phantom was used to provide a homogeneous magnitude, but it yielded distinguishable B z signal inside and outside the anomaly. For the phantom experiment, total imaging time was 10.24 seconds to fill the complete k-space region in the interleaved sense and it was less than 6 seconds to fill the partial k-space region to implement the proposed method. The proposed interleaved partial Fourier strategy for the fast MREIT has a potential to reduce scan times and maintain the information of B z data comparable to what is obtained with the complete k-space data.