A method for MREIT-based source imaging: simulation studies

This paper aims to provide a method for using magnetic resonance electrical impedance tomography (MREIT) to visualize local conductivity changes associated with evoked neuronal activities in the brain. MREIT is an MRI-based technique for conductivity mapping by probing the magnetic flux density induced by an externally injected current through surface electrodes. Since local conductivity changes resulting from evoked neural activities are very small (less than a few %), a major challenge is to acquire exogenous magnetic flux density data exceeding a certain noise level. Noting that the signal-to-noise ratio is proportional to the square root of the number of averages, it is important to reduce the data acquisition time to get more averages within a given total data collection time. The proposed method uses a sub-sampled k-space data set in the phase-encoding direction to significantly reduce the data acquisition time. Since the sub-sampled data violates the Nyquist criteria, we only get a nonlinearly wrapped version of the exogenous magnetic flux density data, which is insufficient for conductivity imaging. Taking advantage of the sparseness of the conductivity change, the proposed method detects local conductivity changes by estimating the time-change of the Laplacian of the nonlinearly wrapped data.


Introduction
Functional MRI, the most commonly used neuroimaging technique, is an indirect method of visualizing changes in cerebral blood flow associated with neuronal activities (Smith 2012). There have been numerous studies to detect endogenous magnetic field induced by evoked neuronal currents from MR magnitude and phase images (Bodurka andBandettini 2002, Luo et al 2011). However, the currently available MRI scanners seem to lack the sensitivity needed to produce any reproducible results (Huang 2014). Electroencephalogram (EEG) and magnetoencephalogram (MEG) source imaging methods entail the measurements of electric potential on the head and magnetic field outside the head, respectively, associated with neuronal current sources inside the head. The solutions of the corresponding inverse problems are not unique, and these neuroimaging methods have their own limitations (Wendel et al 2009).
The effective tissue conductivity inside the brain depends on molecular composition, amounts of intra-and extra-cellular fluids and cellular structure. Cell membranes are known to be highly resistive, hindering the DC current from passing through the cells. The opening of ion channels in cell membranes leads to an increase in the effective tissue conductivity since ion channels allow DC current to pass through the intracellular space. Animal studies have shown that neuronal depolarization leads to a net increase in the tissue conductance (Klivington and Galambos 1967, Velluti et al 1968, Oh et al 2011, Santos et al 2016. On the other hand, cell swelling accompanying epileptic seizures tend to decrease the effective conductivity. Santos et al (2016) characterized and visualized cortical conductivity changes during interictal and ictal activities in the anaesthetized rat. The decreased cerebral tissue conductivity during seizure activities can be attributed to the shrinkage of extracellular space by cell swelling that follows intensive neural activities (Van Harreveld and Schade 1962, Elazar et al 1966, Dietzel et al 1980, Olsson et al 2006.
There have been research efforts to detect conductivity changes associated with neuronal depolarization. If successful, this will provide a new direct neuroimaging method. Recent experimental studies show the feasibility of detecting fast neural evoked activities on rat cerebral cortex using an impedance imaging method employing a planar epicortical electrode array (Oh et al 2011, Aristovich et al 2014. This method may find numerous clinical applications where the invasive procedure to place epicortical electrodes is acceptable. Magnetic resonance electrical impedance tomography (MREIT) produces cross-sectional images of a conductivity distribution inside the human body by injecting currents and measuring induced magnetic flux densities using an MRI scanner. We may apply the technique to visualize local conductivity changes associated with evoked neuronal activities in the brain (Sadleir et al 2010, Woo 2011. The major technical challenge in this MREIT-based source imaging or so called functional MREIT is how to acquire the magnetic flux density data with an sufficient signal-to-noise ratio (SNR) to detect neural activities within an acceptable imaging time.
In MREIT, an almost constant current is injected into an imaging subject to produce a magnetic flux density , , x y z inside the body. The component B z is then found from the measured phase changes in the MR images. MREIT exploits the nonlinear relationship between the conductivity σ and B z via the Biot-Savart law, making a key observation that the Laplacian of the B z data probes changes in the logarithm of the conductivity distribution along any equipotential curve in each imaging slice , Woo and Seo 2008, Seo and Woo 2011, Seo et al 2013b. Given that the conductivity distribution changes with time due to neuronal activity, we consider σ and B z as functions of position ( ) = x y z r , , and time t. The relation between σ ∂ ∂t and ∂ ∂ B t z can then be expressed by the Biot-Savart law: for ∈ Ω r , (1) , Ω represents the imaging domain, and J and E are the current density and electric field, respectively, induced by the injection current. The inverse problem for MREIT source localization is to estimate location of the support of σ ∂ ∂t and its magnitude from knowledge of the measured data ∂ ∂ B t z and the relation (1). Noting that the magnitude of σ ∂ ∂t associated with neural activity is very small (less than a few %), the major challenge facing MREIT source localization is to obtain a sufficiently high SNR of ∂ ∂ B t z subject to a given amount of a given injection current and a fixed acquisition time. The SNR of ∂ ∂ B t z is proportional to that of the k-space data. In this paper, we use a subsampled k-space data set in the phase-encoding direction to increase SNR significantly for a fixed acquisition time. Unfortunately, it is very difficult to extract ∂ ∂ B t z with this sub-sampled k-space data set, because the B z data in the inverse Fourier transform of the sub-sampled k-space data are wrapped not only circularly but also nonlinearly. To overcome this difficulty, we take advantage of the sparsity of its Laplacian ∇ ∂ ∂ B t z 2 , which is closely related to σ ∂ ∂t being assumed to be sparse.
Numerical simulations demonstrate the performance of the proposed method and support its feasibility in future experimental studies. We will also address some technical issues for MREIT source imaging to advance it toward in vivo animal and human experiments.

Method
Assume that an imaging object occupies a three-dimensional domain Ω and that its conductivity distribution is σ. In MREIT, we inject a pulsed DC current of I mA into Ω using a pair of surface electrodes + E and − E to produce the time-independent current density J, electric field E and magnetic flux density B inside Ω as shown in figure 1. Since the conductivity ( ) σ t r, depends on position ( ) = x y z r , , and time t due to the conductivity change associated with neuronal activity, the static fields J E B , , also depend on time t due to the conductivity change σ ∂ ∂t . Throughout this section, we fix z = z 0 and consider all these quantities only on a slice { } Ω = Ω ∩ = z z z 0 0 of Ω. Assuming that the z-axis is the axial magnetization direction of the scanner, the MR signal contains only the information of the z-component of B (Joy et al 1989, Scott et al 1991, 1992. To be precise, the DC current of I mA is injected in pulses with the pulse width of T c between the 90° and 180° RF pulses and also between the 180° RF pulse and the read gradient, as shown in figure 1. The corresponding k-space MR signal due to the injection current can be described as , e e d d , where M is the MR magnitude image, γ = × ⋅ − 26.75 10 rad T s 7 1 the gyromagnetic ratio of hydrogen, and δ any systematic phase artifact. If k-space sampling is designed to meet the Nyquist criteria (Seo and Woo 2013a), the inverse discrete Fourier transformation (DFT) of the k-space MR signal ( ) S k k t , , x y provides the complex MR image , e e .
T B x y t xy i ,, i , . To be precise, we assume that ( ) M x y t , , is supported in the 2 and M is periodic in the sense that ( 1, , , . Let S N be a sub-sampled k-space by a factor N in the phase encoding direction (y-direction). From the Poisson summation formula (Seo and Woo 2013a), the inverse DFT of the subsampled k-space data S N produces the following folded image: , e x y i , . Figures 2(a) and (b) illustrate the real and imaginary parts of M N with N = 4 respectively.
In this paper, we use the sub-sampled k-space data set S N to increase its SNR significantly with a fixed acquisition time. The data acquisition time of k-space data S N is inversely proportional to the factor N due to the time consuming phase encoding, so that S N has higher averaging cycle with a fixed scan time.
Taking the time derivative on both sides of (4), we obtain To overcome this difficulty, we take advantage of the sparsity assumption of 2 z is related to σ ∂ ∂t through the following identity, which can be obtained by taking the Laplacian to the z-component of Biot-Savart law (1): Since the support (7) is contained in D, u can be viewed as a solution of the Poisson equation (7) with its source term supported in D. Hence, according to the standard PDE theory, = − ∇u E decays rapidly with the distance from D. Hence, we may assume that is known. This background image can be computed from the standard MREIT method. In this paper, we compute ψ = ∂ ∂ B t z by minimizing the following energy functional The first integral in (8), called the fidelity term, forces the residual ∑ + which is an approximation of (5), to be small. The second integral of (8), called the regularization term, enforces the sparsity of ψ small away from the small source region D. Here, λ is the regularization parameter which controls the tradeoff between the residual norm and the sparsity. For technical simplicity, instead of solving the misfit functional (8), we consider the following modified version: A simple calculation yields that a minimizer ψ satisfies the following partial differential equation approximately: Hence, it remains to select D out of the N-folded area D N . For the selection, we need to use the fact that the change ∂ ∂ in (1) can be approximated roughly by the following function ζ k : where u 0 is the potential corresponding to a reference conductivity σ 0 : and v k is a solution of Poisson's equation: Figure 4. Schematic diagram of the proposed MREIT-based source imaging. A head is placed in an MRI scanner. We use a sub-sampled k-space data set in the phase-encoding direction, in which a nonlinearly wrapped version of the exogenous magnetic signal is available.
Here, the reference conductivity σ 0 can be obtained by the standard MREIT algorithm, under the assumption that σ ∂ ∂t is small, and u 0 is computed by solving the forward problem (13). From (12)-(13), we observe that ζ k depends on the positions of electrodes + − E E , . Now, we use the fact that, if D = D k , the spatial distribution of ζ k is similar to that of 1 which minimizes the difference: Here, η is a scaling factor related to the percentage increase in conductivity. Figure 4 shows a schematic diagram of the proposed MREIT-based source imaging method.

Numerical simulation results
We carried out numerical simulations to test the feasibility of the proposed MREIT-based source localization method, which is described in the previous section. The procedure of the proposed method is summarized roughly as follows: (i) Using MRI scanner and the standard MREIT method, we get the MR magnitude image M and the reference information of B z (x, y, 0) and ( ) σ x y , , 0 on Ω z0 at time t = 0. (ii) Get the sub-sampled k-space data S N and the N-folded image of M N via inverse FFT for each time ⩾ t 0 . (iii) Solve the equation (17) iteratively with the initial guess ψ = 0 0 : For = n 1, 2, and each t, solve We construct a realistic three dimensional head model as shown in figure 1 (a). We place a spherical source domain D inside the subject, whose conductivity changes with time as σ = + t 0.09 0.0009 for ∈ t 3, 5, 8. The conductivity distribution when t = 5 of the slice Ω z0 is shown in figure 1(c). The domain of head was meshed into finite elements as in figure 1(b). Two pairs of surface electrodes are attached to reconstruct the background conductivity. We generate the simulated data B z using the forward solver of MREIT (Lee et al 2003). We use the standard MREIT algorithm, called the harmonic B z algorithm , to compute the reference conductivity distribution using two different reference B z data at t = 0, as shown in figure 6. Next, we consider the head model in figure 7, which is a bit different in terms of the size of electrodes. Figure 8 illustrates the temporal distribution of the conductivity. Through this pair of electrodes we inject a current with amplitude of 1 mA. Temporarily, we do not add any noise in the obtained data. Figure 9 illustrates the simulated data To test the noise tolerance of the algorithm, we add Gaussian random noise to the simulated data ( ) M x y t , , N . From the noise analysis in Sadleir et al (2005), Scott et al (1992), the noise standard deviation s B of B z is given by , where SNR is the signal-to-noise ratio of the MR magnitude image. Using the noisy B z (x, y, t), we can produce the sub-sampled k-space data and validate our algorithm. In this experiment, = SNR 2000 and = T 50 c ms.  Table 1 illustrates the process of selecting D from D 4 at t = 5 by evaluating the difference Err k defined by x y x y , , 5 , , 0 t , supported in a source region in Ω z0 . To enhance the image contrast, we multiply σ ∂ ∂t by 20.
x y x y x y min 1 , , 1 , , d d , where ν is a scaling constant given by ( . In this experiment, we choose η as 5% increase in conductivity: where ( ) σ x y , , 0 is the reference conductivity reconstructed from the standard MREIT algorithm . According to table 1, we have < < Err Err Err 1 2 0 in both cases of SNR = ∞ and 2000, and hence D = D 1 is selected.
Next, we performed numerical experiments for an F-shaped source region as shown in figure 12. Figures 13 and 14  We can see from numerical experiments in this section, the proposed algorithm is capable of visualizing the time changes of the conductivity less than 5% from the sub-sampled k-space data. The proposed method of detecting the support of D N with the undersampled k-space data works well with different simulation models. The performance does not change with a more realistic model with heterogeneous background. The performance of selecting D depends on the distance between D and electrode locations. Generally speaking, the nearer D is to the electrodes, the better the performance.

Phantom experiment
To experimentally test the proposed method, we built a saline-filled cylindrical phantom where an insulated carbon wire is immersed. Figure 15 shows the phantom with its diameter of 250 mm and height of 500 mm. The diameter of the carbon wire was 1 mm including insulation and 0.2 mm without insulation. We injected 1 mA imaging current through a pair of carbon electrodes attached on the boundary of the imaging plane { } Ω ∪ = z z 0 . The size of the carbon electrode was × 10 10 mm. We additionally passed 30 μA current through the insulated carbon wire to generate an equivalent current dipole, mimicking neural activity. The wire was partially exposed to saline in the imaging plane so that the dipole current could spread out inside the saline water. A 9.4T MRI scanner was used with a gradient echo-based MREIT pulse sequence (Seo et al 2013b). Imaging parameters were TR/TE = 300/3 ms, = T 27 c ms,   For the reader's convenience, we briefly explain the reason why the dipole current corresponds to an equivalent conductivity change. The 30 μA current spreads out in the imaging slice Ω z0 over a neighborhood of the uninsulated carbon wire segment. Now let us imagine that there occurred a conductivity change in the region D which is the neighborhood of the uninsulated wire segment. Denoting by ∂ ∂ u t the corresponding potential change, the time-difference

Conclusions and future works
Neuronal activities increase tissue conductivity values in a local region regardless of the directions of associated neuronal currents. Since this spontaneous relation is direct and instantaneous, neural activity can be in theory directly mapped by visualizing changes in the conductivity distribution. MREIT is capable of providing high-resolution conductivity images with a pixel size of about 1 mm, provided that the SNR of B z data is sufficiently high. For MREIT-based neuroimaging, the most challenging issue is to achieve a sufficiently high SNR of the acquired MR data to allow the detection of conductivity changes of less than 5%. Such an improvement of the SNR must occur under the constraint of a fixed total scan time. Assuming a certain MRI scanner with a given performance level, we may improve the SNR by increasing the number of averages. In order not to increase the total scan time, the proposed method uses sub-sampled k-space data. This is accomplished through expediting the data acquisition by skipping the time-consuming phase encoding lines. Using 20% k-space data and one injection current, for example, we can reduce the data acquisition time by about 90%, and hence the MR data can be averaged 10 more times.
Although the sub-sampled k-space data may be insufficient to obtain the complete B z data, it is capable of computing a folded version of ∇ ∂ ∂ B t z 2 , from which we can visualize local  perturbations of conductivity values. The results of numerical experiments show that 20% k-space data can successfully detect local conductivity changes associated with neural activities. The proposed method can either significantly reduce the scan time without deteriorating the quality of source images or improve the image quality without increasing the scan time. Seo