SNR-Enhanced, Rapid Electrical Conductivity Mapping Using Echo-Shifted MRI

Magnetic resonance electrical impedance tomography (MREIT) permits high-spatial resolution electrical conductivity mapping of biological tissues, and its quantification accuracy hinges on the signal-to-noise ratio (SNR) of the current-induced magnetic flux density (Bz). The purpose of this work was to achieve Bz SNR-enhanced rapid conductivity imaging by developing an echo-shifted steady-state incoherent imaging-based MREIT technique. In the proposed pulse sequence, the free-induction-decay signal is shifted in time over multiple imaging slices, and as a result is exposed to a plurality of injecting current pulses before forming an echo. Thus, the proposed multi-slice echo-shifting strategy allows a high SNR for Bz for a given number of current injections. However, with increasing the time of echo formation, the Bz SNR will also be compromised by T2*-related signal loss. Hence, numerical simulations were performed to evaluate the relationship between the echo-shifting and the Bz SNR, and subsequently to determine the optimal imaging parameters. Experimental studies were conducted to evaluate the effectiveness of the proposed method over conventional spin-echo-based MREIT. Compared with the reference spin-echo MREIT, the proposed echo-shifting-based method improves the efficiency in both data acquisition and current injection while retaining the accuracy of conductivity quantification. The results suggest the feasibility of the proposed MREIT method as a practical means for conductivity mapping.


Introduction
Electrical conductivity inside the living body is determined by a number of factors, including its underlying cellular structure, ion mobility and concentration, and molecular composition [1,2]. Hence, reliable measurements of electrical conductivity would provide important insights into the physiology of biological tissues in health and disease. In fact, electrical conductivity information has been employed in a range of applications such as EEG and MEG neuronal source localization [3][4][5], quantitative monitoring of neuronal depolarization [6], ion mobility imaging [7][8][9], estimation of current distribution during therapeutic electrical stimulation [10][11][12][13], and evaluation of brain abnormalities [14].
Electrical impedance tomography (EIT) is a non-invasive imaging technique that permits the estimation of the electrical conductivity distribution inside an imaging object [15,16]. In EIT, current injection is applied to the imaging object through multiple surface electrodes, and induced voltages are then measured to reconstruct cross-sectional conductivity maps based on a nonlinear inverse solution. However, the boundary voltage is insensitive to internal conductivity variations, making the inverse problem highly ill-posed and resulting in inaccurate conductivity estimates in regions away from the measurement electrodes [16]. Additionally, it is difficult to achieve high spatial resolution conductivity mapping with a limited number of electrodes [16].
As an alternative, an MRI-based technique, typically referred to as magnetic resonance electrical impedance tomography (MREIT) [17], has been suggested to obtain electrical conductivity maps with sufficiently high spatial resolution and accuracy. In the method, an imaging object is subjected to a series of short current pulses in sync with an MR pulse sequence, and the induced magnetic flux density (B = [B x , B y , B z ]), which perturbs the local magnetic field, is estimated via phase analysis of acquired MR images. Once all three components of B are measured separately by rotating the object twice inside the MRI scanner, the current density (J) distribution is calculated using Ampere's law J = ∇ × B/µ (µ: magnetic permeability), leading to conductivity estimation [18,19]. More recently, a harmonic B z -based MREIT [20,21] has been introduced to avoid technical difficulties in object rotation in the J-based approaches. Given two sets of B z (rather than B) obtained with current injection in two orthogonal directions, the harmonic B z method reconstructs conductivity distribution (σ) by exploiting the relationship between ∇ 2 B z and ∇σ, thus obviating the need for object rotation in MREIT experiments. Despite the advantage, the method may suffer from noise amplifications in the numerical computation of ∇ 2 B z and thus requires a sufficiently high signal-to-noise ratio (SNR) of B z .
Previous studies analyzing B z SNR have shown that the standard deviation of B z is inversely proportional to the SNR of a magnitude image, duration of a current pulse, and further phase sensitivity of an MR pulse sequence to the injection current [22,23]. In this regard, the spin echo (SE) pulse sequence has been widely accepted in MREIT, as it provides a high SNR in magnitude images and allows a long time for current injection (TC) between an excitation RF pulse and signal readout. To further enhance the B z SNR without compromising imaging efficiency, multi-echo variants of SE imaging have also been explored [24,25]. Nevertheless, the utility of these SE-based methods has been limited by the impractically long scan times and low phase sensitivity. To address these issues in SE-based MREIT, alternating steady-state free precession (SSFP) MREIT has been suggested [23,26]. Compared with SE MREIT, the method enables rapid conductivity imaging while yielding high phase sensitivity resulting from the nonlinear behavior of SSFP signals in response to alternating current injection. However, in solving the corresponding nonlinear inverse problem, alternating SSFP MREIT requires knowledges of tissue relaxation times and transmit RF inhomogeneities, for which additional measurement scans need to be performed, making the B z estimation procedure rather complicated.
In this work, we aimed to overcome the above-mentioned limitations in current MREIT techniques by developing a multi-slice interleaving echo-shifted steady-state incoherent (ESSSI) imaging-based MREIT method so as to achieve current-efficient, high-speed conductivity imaging and direct extraction of B z values from acquired-image signals. Free induction decay (FID) signals are shifted in time over multiple imaging slices before forming an echo, thereby being exposed to a multitude of current pulses and accumulating current-induced phases successively, leading to an elevated B z SNR without an increase in TC. Number of echo-shifting (NES) values were optimized using numerical simulations. Experimental studies were performed in phantoms to evaluate the effectiveness of the proposed method in comparison to conventional SE MREIT.

Multi-Slice ESSSI Imaging with Current Injection
A timing diagram of the proposed multi-slice interleaving ESSSI imaging method with current injection is shown in Figure 1. Equidistant, spatially selective excitation RF pulses with constant flip angles (α•) are successively applied to corresponding, interleaved imaging slices. To achieve maximal incoherence between isochromats at the end of each time of repetition (TR), RF spoiling is applied with the quadratic phase cycling scheme [27]: where φ is the RF phase, n is the RF index for each slice, and ψ is the RF phase increment, typically set to 117 • so as to produce a signal intensity comparable to that achieved by ideal spoiling [27]. The time integral of gradient pulses between any two neighboring RF pulses is kept identical in all three directions (x, y, z) to avoid artifacts otherwise resulting from temporally varying spin dephasing. Further, a pair of spoiler gradients is inserted before and after each signal readout (see below for details) while a short, unipolar current pulse is applied between each pair of RF pulses. With the above configuration, a current-encoded incoherent steady state is established by increasing the number of RF pulses. Figure 1. A timing diagram of the proposed multi-slice interleaving echo-shifted MREIT pulse sequence for NES = 1 (GS: slice selection gradient; GP: phase-encoding gradient; GR: frequencyencoding gradient). Note that the FID signals are shifted in time over multiple imaging slices to form an echo at TE eff rather than TE. Note also that current pulses with a constant amplitude are applied during an interval between two neighboring RF pulses. Thus, the effective time of current injection (TC eff ) that each echo signal experiences is longer than the TC.
In the present pulse sequence, an FID signal produced in an imaging slice is shifted in time to the following blocks led by RF pulses for other imaging slices (Figure 1), forming an echo at an effective time of echo (TE eff ): where T r is the time interval between two neighboring RF pulses. With this echo-shifting scheme, each FID signal is exposed to more than one current pulse without an increase in TR, while each current pulse is effectively shared by NES slices. Hence, compared with non-echo-shifting multi-slice SSI imaging, the proposed method enables rapid and current-efficient conductivity imaging. A detailed analysis of the effect of echo-shifting on the B z SNR is provided in the next section.
To achieve echo-shifting over multiple slices, two spoiler gradient pulses (A and B in Figure 1) are positioned before and after each signal readout with the following design criteria: Here, γ is the gyromagnetic ratio; M A and M B are the zeroth moments of the spoilers A and B, respectively; ∆x is the voxel size; and m is the echo-shifting index ranging from 0 to NES-1. Equation (3) ensures that for any NES value, all FID signals experience spoilerinduced spin dephasing over 2π within each voxel so as to avoid their interferences within the window of signal readouts, thereby preventing resultant image artifacts. Given the above considerations, the minimum size of each spoiler gradient is calculated by: and: The total amount of gradient-induced spin dephasing during each TR is then given by 2π times the number of excitation slices, which is sufficiently large for effective RF spoiling [27,28].

B z Estimation and SNR Analysis
With the integration of injection currents into the ESSSI pulse sequence, the effective TC (TC eff ), the time for which an FID signal encodes the current-induced magnetic field (B z ), is defined as: TC e f f = TE + NES·TC As a result, the phase of a forming echo at TE eff is given by: while its magnitude signal is determined by a transverse relaxation time constant, T * 2,e f f , expressed as: Here, φ b and φ c are phases due to the background magnetic field (∆B 0 ) and B z , respectively. In Equation (8), it is assumed that ∆B 0 lends itself to a Lorentzian spectral distribution while B z is piecewise constant.
To selectively extract φ c values from acquired signals, in this work two separate ESSSI MREIT scans were performed, in which current pulses with positive (+I) and negative (−I) polarities were applied, respectively, with the same amplitudes and durations. The corresponding steady-state signals in a voxel, S + and S − , can be written as: where S M is the signal magnitude, M 0 is the magnetization in the thermal equilibrium state, and E 1 = e −TR/T 1 . The current-induced B z value in each voxel is then calculated by: Here, r is the voxel position and arg(·) is the operator that yields the phase of its argument.
Noise analysis in MREIT [22,23] reveals that the standard deviation of B z (SD B z ) can be , where Y M is the SNR of magnitude image and ξ is the phase sensitivity of an MR pulse sequence to the injection current. In ESSSI MREIT, both TC and Y M are functions of NES (Equations (6) and (9)), while the image phase accumulates linearly with ξ = 1. Hence, SD B z in the proposed method can be written as: The above equations imply that the two factors, TC eff and Y M,ESSSI , conflict with each other in maximizing B z SNR. Furthermore, with an Ernst flip angle employed, the signal intensity (S M in Equation (9)) is determined by TE eff and TR. Given these considerations, for a given number of imaging slices, the two imaging parameters, NES and T r , are critical determinants of B z SNR and are optimized using numerical simulations in the next section.

Numerical Simulations
Numerical simulations were performed in the proposed ESSSI MREIT process to determine an optimal combination of NES and T r that yields a maximal B z SNR for a range of tissue relaxation times. To this end, the gain of B z SNR (η) achievable with the echo-shifting scheme relative to non-echo-shifting (i.e., NES = 0), was defined as: With increasing NES values from 0 to 7 and varying T r times from 5 ms to 15 ms, contour plots of η were generated for three different tissues with T * 2 values of 40, 70, and 100 ms, respectively. In the ESSSI pulse sequence, as NES increases the size of the two spoiler gradients needs to be enlarged, thereby increasing the minimum possible T r under the limit of a maximum gradient amplitude (G max ). This systematic lower bound of T r was calculated for each NES value and was indicated in the contour plots. The simulation parameters were: TE = 5 ms, number of slices = 8, current-induced B z = 20 nT, RF pulse duration = 1 ms, and T 1 /M 0 = 500 ms/1.0. Both TR and TC were adjusted with T r . ∆x = 1.5 mm and G max = 28 mT/m being assumed.
To evaluate the proposed method's performance over the conventional SE MREIT approach in terms of B z SNR per unit scan time, simulations were performed using the following equation: where SNR Bz is the SNR of the B z values calculated in each imaging method. Hence, ξ represents the B z SNR efficiency of the proposed ESSSI technique relative to SE imaging.
Here, SD Bz,ESSSI was obtained using Equation (11), while SD Bz,SE was derived using For T 2 and B z values in the ranges of 30-150 ms and 0-50 nT, ξ was computed by varying NES from 0 to 5, leading to contour plots of ξ with respect to T 2 and B z . The simulation parameters in the ESSSI pulse sequence were kept identical to those above, while those in SE imaging were: TR = 1000 ms, TE = 40 ms, and TC = 32 ms.

Experimental Studies
Experimental studies were performed at 3 T (Siemens Trio, Erlangen, Germany) in three different custom-made cylindrical conductivity phantoms (hereafter referred to as phantoms A, B, and C). A 32-element head coil was used for signal reception. Phantoms A and B consisted of a homogeneous agar gel object with σ = 2.1 S/m (CuSO 4 , 0.5 g/L; NaCl, 10 g/L) but differed in their agar concentrations (A: 25 g/L; B: 10 g/L); thus, they simulated tissues with different T * 2 values. Phantom C was composed of two cylindrical agar gels surrounded by saline solutions, in which the agar objects with larger and smaller diameters and the background solutions presented σ = 2.79 S/m (CuSO 4 , 1.25 g/L; NaCl, 12.5 g/L; agar, 25 g/L), σ = 1.14 S/m (CuSO 4 , 1.25 g/L; NaCl, 3.75 g/L; agar, 25 g/L), and σ = 0.2 S/m (CuSO 4 , 1 g/L; NaCl, 1.5 g/L), respectively. Prior to conductivity imaging, a multi-echo gradient echo scan (TR = 500 ms, flip angle = 60 • , TEs = 2.2/4.0/6.2/8.4/10/6/15/20/30/50/80 ms) was performed to measure T * 2 via a linear regression. The obtained T * 2 maps served as a reference when experimentally validating the relationship between NES and B z SNR simulated above.
Two sets of MREIT data were acquired with opposite polarities of the injection current, yielding a B z map via Equation (10). In each measurement, a single line of the projection signal was collected without current injection, and then its phase was demodulated from current-encoded signals so as to correct for any systematic global phase offset. The above procedure was repeated with the direction of current injection rotated by 90 • . Finally, given the two sets of B z estimates (B z,1 , B z,2 ) in the orthogonal directions, a conductivity map was constructed using the CoReHa software package [29], in which the harmonic B z algorithm is implemented in a semi-automatic manner.
To investigate the effect of echo-shifting on SD B z and for conductivity estimates, data were acquired in phantoms A and B using the proposed multi-slice interleaving ESSSI pulse sequence with increasing NES from 0 to 5 (increment: 1), leading to four sets of imagesthe magnitudes of S + (|S + |), B z,1 , B z,2 , and σ. The magnitude SNR and SD B z in seven regions-of-interest (ROIs) of the |S + | and B z,1 images, were estimated using the NEMA-I method [30]. To this end, two sets of images were obtained independently using the same imaging parameters, and the standard deviation of their difference in each ROI was computed, serving as a noise statistic. Given the SD B z measurements, η was then calculated using Equation (12). Furthermore, standard deviations of σ (SD σ ) in seven ROIs were measured and compared across the examined NES values. The imaging parameters were as follows: field-of-view = 180 × 180 mm 2 , slice thickness = 4 mm, number of slices = 6, in-plane matrix size = 128 × 128, readout bandwidth = 500 Hz/pix, TE = 5.5 ms, flip angle = 25 • , number of signal averages = 2, and I = +10/−10 mA. For a given NES value, both TR (and thus T r ) and TC were adjusted to minimum and maximum possible values, respectively, which were achievable with the above scan parameters: TR/TC = 50/6 ms (NES = 0, 1), 54/7 ms (NES = 2), 58/7.5 ms (NES = 3), 62/8 ms (NES = 4), and 68/9 ms (NES = 5).
Data were collected in phantom C using the proposed ESSSI MREIT method with increasing NES values from 0 to 3 (increment: 1), along with the conventional SE MREIT method for comparison. With both methods, three sets of images, |S + |, B z,1 , and σ, were presented, respectively. Scanning parameters common to both imaging techniques were: field-of-view = 150 × 150 mm 2 , slice thickness = 4 mm, number of slices = 8, in-plane matrix size = 128 × 128, and I = +10/−10 mA. Parameters specific to the proposed method were: TR = 72 ms, TE = 4 ms, TC = 6 ms, flip angle = 15 • , readout bandwidth = 500 Hz/pix, number of signal averages = 10, and scan time = approximately 6.5 min. The parameters specific to conventional SE MREIT were: TR = 1000 ms, TE = 40 ms, TC = 32 ms, number of signal averages = 3, and scan time = approximately 26 min. The experiments performed in this study and the relevant imaging parameters are summarized in Table 1.     (Figure 3f). For all of the parameter combinations examined, ξ is larger than 1, thus suggesting that compared with the SE technique, the proposed ESSSI pulse sequence presents a substantially higher B z SNR efficiency. When a relatively smaller value of NES is employed (NES = 0-2), ξ decreases more rapidly with increasing T 2 more than it does with increasing B z (Figure 3a-c). In contrast, with moderately large NES values (NES = 3, 4), ξ changes predominantly toward the direction of B z (Figure 3d,e). Finally, when NES = 5, ξ increases with increasing T 2 and decreasing B z (Figure 3f).   (Figure 6d-f). As the NES values increase from 0 to 5, the magnitude image SNR drops approximately exponentially in both phantoms, while its decreasing rate in phantom A is relatively higher than that in phantom B (Figure 6a vs. Figure 6d). With increasing NES, η rises and falls at the NES values of 3 and 4 in phantom A (Figure 6b) and phantom B (Figure 6e), respectively, which experimentally validates the results of the numerical simulations (Figure 2a,b). As expected, SD σ presents a reversed pattern in comparison to the variations of η against NES. Figure 7 compares three sets of images, |S + |, B z,1 , and σ, in phantom C obtained using the ESSSI MREIT method with NES = 0-3 and conventional SE MREIT. As NES values increase from 0 to 3, the proposed method results in decreased signal intensity for the magnitude images but yields a gradual reduction in noise for the B z estimates, leading to an elevated conductivity contrast. Furthermore, compared to the reference, namely the conventional SE MREIT method, the conductivity map obtained using the proposed method with NES = 3 depicts a similar level of conductivity contrast. Although the results for NES > 3 in the proposed method are not shown, the above results from numerical simulations ( Figure 2) and homogeneous phantom experiments (Figures 4-6) suggest that B z SNR and conductivity reconstruction would be degraded after a certain NES threshold.

Discussion
This work introduces a new MREIT method based on an echo-shifted steady-state incoherent imaging pulse sequence for rapid and current-efficient conductivity mapping without an apparent loss of measurement accuracy. In the proposed method, FID signals are shifted over one or more imaging slices and experience an accordingly increased number of injection current pulses. Therefore, the effective current duration is lengthened, leading to enhanced B z SNR values relative to the non-echo-shifted counterpart. Nonetheless, the number of echo-shifts is limited by the T * 2 -related signal loss. Hence, the optimal NES value varies with the intrinsic relaxation of the tissues, as shown via numerical simulations (Figures 2 and 3) and further validated by phantom experiments (Figures 4-6).
In the original implementation of echo-shifted imaging [31] and its variants thereafter [32,33], FIDs were shifted over the direction of phase-encoding lines for the same imaging slice. Hence, in these early techniques, steady-state signals were weighted by (cos α/2) 2NES and remained at a relatively low level if a high NES value was employed. In contrast, since echo-shifting in the present pulse sequence is integrated into the multi-slice interleaving data acquisition process, the FID signal for a given slice is unaffected by RF pulses for other imaging slices. Furthermore, the multi-slice interleaving configuration employs a long TR to accommodate multiple RF pulses, leading to an enhanced level of steady-state signals. Hence, compared with the original echo-shifted imaging technique, the proposed method substantially elevates the SNR of magnitude images and accordingly the B z SNR.
As an alternative to the T r -periodic spoiler gradient pulses in the proposed method, echo-shifting with (NES + 1)T r -periodic spoilers can be considered [32], such that the size of spoiler A (M A ; Figure 1) varies across the pulse train with a period of (NES + 1)T r , while that of spoiler B is set to zero (i.e., M B = 0). For example, if NES = 1, M A changes its sign alternately over T r with its absolute moment held constant. Compared with the present implementation of the spoilers, the (NES + 1)T r -periodic approach shortens the minimum possible T r and potentially enhances the scan efficiency. Nonetheless, with reduced T r , the duration for the current pulse needs to be decreased accordingly, leading to reduced B z SNR. Additionally, in the presence of any uncompensated residual eddy currents, (NES + 1)T rperiodic spoilers make the signal phase vary periodically along the direction of phase encoding, potentially resulting in ghosting artifacts in the obtained images.
In the proposed MREIT method, two measurements need to be performed separately with opposite polarities of injecting current pulses to eliminate the background magnetic field in the estimation of current-induced B z . Instead, alternating current injection may be applied by continuously switching the polarities of current pulses across the entire pulse train to make the B z estimation relatively more immune to global phase offset over scans, as well as potential magnetic field drifts during each measurement. However, alternating current injection in sync with an echo-shifted pulse sequence would be undesirable because current-induced phases are successively cancelled out along the pulse train if the echoshifting mode is turned on.

Conclusions
In conclusion, it was demonstrated that a new echo-shifted steady-state incoherent imaging-based MREIT method enables rapid, high-resolution conductivity mapping. The echo-shifting strategy in combination with the multi-slice interleaving data acquisition approach allows efficient encoding of the current-induced magnetic field, thereby enhancing the SNR of B z without an explicit increase in the number of current pulses or the current duration. It is expected that the proposed method will provide a novel and highly efficient way to measure conductivity information in MREIT studies.  Data Availability Statement: Data will be made available on request to the corresponding author.