Current-density imaging using ultra-low-field MRI with zero-field encoding

Magnetic Resonance Imaging 32 (2014) 54–59

Contents lists available at ScienceDirect

Magnetic Resonance Imaging journal homepage: www.mrijournal.com

Current-density imaging using ultra-low-field MRI with adiabatic pulses Jaakko O. Nieminen a,⁎, Koos C.J. Zevenhoven a, Panu T. Vesanen a, Yi-Cheng Hsu a, b, Risto J. Ilmoniemi a a b

Department of Biomedical Engineering and Computational Science, Aalto University School of Science, P.O. Box 12200, FI-00076 AALTO, Finland Department of Mathematics, National Taiwan University, Taipei, TW-10617, Taiwan

a r t i c l e

i n f o

Article history: Received 14 February 2013 Revised 25 June 2013 Accepted 21 July 2013 Keywords: Ultra-low-field MRI Current-density imaging Adiabaticity Prepolarization Rotation-free

a b s t r a c t Magnetic resonance imaging (MRI) allows measurement of electric current density in an object. The measurement is based on observing how the magnetic field of the current density affects the associated spins. However, as high-field MRI is sensitive to static magnetic field variations of only the field component along the main field direction, object rotations are typically needed to image three-dimensional current densities. Ultra-low-field (ULF) MRI, on the other hand, with B0 on the order of 10–100 μT, allows novel MRI sequences. We present a rotation-free method for imaging static magnetic fields and current densities using ULF MRI. The method utilizes prepolarization pulses with adiabatic switch-off ramps. The technique is designed to reveal complete field and current-density information without the need to rotate the object. The method may find applications, e.g., in conductivity imaging. We present simulation results showing the feasibility of the sequence. © 2014 Elsevier Inc. All rights reserved.

1. Introduction Magnetic resonance imaging (MRI) allows noninvasive measurement of electric current density in an object [1–8]. The current density produces a magnetic field which affects the spin dynamics and causes changes to the MR phase images. Because the spatial origin of the signals can be encoded with high spatial accuracy using standard MRI encoding methods, the magnetic field inside the sample can be found without facing an ill-posed inverse problem, which occurs, e.g., when the current-generated magnetic field is measured only outside the sample [9]. Existing methods can produce images of static [1–3] or alternating [4–8] currents. The methods to image static currents or static magnetic fields rely on measuring changes in the magnetic field along the MRI main field B0. Thus, if all three components of the magnetic field are needed, the sample typically has to be rotated at least twice. The need for sample rotations can be relaxed with so-called radio-frequency current-density imaging methods, which are used for imaging time-varying current densities [4,5]. Despite the limitations, current-density imaging (CDI) can be used, e.g., to estimate the conductivity of an object by applying external currents through it [10–15]. Recently, it was demonstrated that the conductivity of an object can be imaged also without applying any external currents [16–18]. One example of MRI-based current-density imaging where no external current is applied is direct neuronal imaging (DNI). ⁎ Corresponding author. E-mail address: jaakko.nieminen@aalto.fi (J.O. Nieminen). 0730-725X/$ – see front matter © 2014 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.mri.2013.07.012

In DNI, the aim is to detect neuronal activity inside the brain by measuring its direct effect on MR signals [19–26]. In contrast to the indirect observation of brain activity in functional MRI (fMRI) [27], DNI is based on detecting the weak variations in the local magnetic fields inside the brain caused by neuronal activity. However, the neuronal magnetic fields are much weaker than the hemodynamic-related magnetic field changes measured by fMRI, or the fields related to externally applied currents, making DNI difficult. In this paper, we present a method to image static magnetic fields and current densities with prepolarized ultra-low-field (ULF) MRI with B0 on the order of 10–100 μT [28–31]. The approach, which we refer to as adiabatic current-density imaging (A-CDI), is based on adiabatic polarizing pulses and is designed to reveal complete field and current-density information without the need to rotate the object. In an accompanying paper [32], we present an alternative approach, called zero-field-encoded current-density imaging (ZCDI), for achieving the same goal.

2. Methods In prepolarized MRI, the sample is magnetized in a field higher than that used for spin precession [33]. This approach is commonly used to increase the signal strength in ULF MRI, where the prepolarizing field can be, e.g., 10–100 mT. Typically, the signals are detected with highly sensitive superconducting quantum interference devices (SQUIDs). To reduce external noise, ULF-MRI devices are often operated inside magnetically shielded rooms.

J.O. Nieminen et al. / Magnetic Resonance Imaging 32 (2014) 54–59

In the following, we will present the A-CDI method for measuring a current-density field J using ULF MRI. The method will concentrate on finding the magnetic field BJ due to J, which can be used to calculate J simply by 1 J ¼ ∇  BJ ; μ0

ð1Þ

where μ0 is the permeability of vacuum. Spins can be manipulated adiabatically or nonadiabatically [34]. In the adiabatic case, the total magnetic field changes slowly and the magnetization can follow its direction. In contrast, if the field changes abruptly, the magnetization is not able to follow it, and the change is nonadiabatic. Assume the magnetic field BJ is static and present during and after the initial polarization period. Now, we can switch off the polarizing field Bp adiabatically: when Bp reaches zero, the magnetization will point along BJ (if BJ = 0, the magnetization will stay along Bp). Finally, we can apply a measurement field, e.g., along Bp, and encoding gradients to image the transverse magnetization. The resulting image will contain information about BJ. 2.1. Waveform for the polarizing field In the following, we design the adiabatic decay waveform for the polarizing field. During the polarizing-field switch-off, the sample magnetization m stays parallel to the total magnetic field B(t) = Bp(t) + BJ, provided that the adiabatic condition, ωL ðt Þ≫jdα=dt j;

ð2Þ

holds at all times [34]. Here, ωL(t)=γB(t) is the Larmor angular frequency in B(t), γ is the gyromagnetic ratio, dα/dt defines the angular rate at which the field orientation changes, and t is time. First, let us consider a special case where Bp(t) and BJ are orthogonal. Now, we can find the angular velocity by noting that f(t) = Bp(t)/BJ = tan α(t); then, we get dα 1 df ðt Þ ¼ : dt 1 þ f 2 ðt Þ dt

ð3Þ

For the Larmor frequency, we obtain qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ωL ðt Þ ¼ γ B 2p ðt Þ þ B 2J ¼ ωJ f 2 ðt Þ þ 1;

ð4Þ

ð5Þ

In practice, we can find an adiabatic polarizing field waveform by setting R(t) = p, with some constant p ≫ 1. With the initial condition f(0) = f0, we get vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2 u u f 0p u q ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi −ω t u J u f 20 þ 1 f ðt Þ ¼ u !2 ≈ u u f 0p tp2 − qffiffiffiffiffiffiffiffiffiffiffiffiffi ffi−ωJ t f 20 þ 1

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 u  u p−ωJ t u u  2 ; t p2 − p−ωJ t

where we noted that f0 ≫ 1. From Eq. (6), we see that the polarizing field has been switched off when t = tend with t end ¼

ωJ

f 0p ffi ≈p=ωJ : qffiffiffiffiffiffiffiffiffiffiffiffiffi f 20 þ 1

ð7Þ

In Eqs. (6) and (7), the approximations hold when f0 ≫ 1, i.e., the polarizing field before the switch-off is much stronger than the field to be imaged. The above derivation assumed Bp⋅ BJ = 0. Thus, if Bp and BJ are non-orthogonal, the goodness of the adiabaticity differs from that given by Eq. (5). However, in all cases, the performance of the ramp can be tuned by varying its duration. 2.2. From magnetization to image For imaging BJ, we suggest the following sequence (Fig. 1). First, a polarizing field is applied to magnetize the sample. Then, the polarizing field is reduced adiabatically to zero while BJ is on; m turns parallel to BJ. After this, BJ is switched off and a homogeneous measurement field B0 is switched on nonadiabatically. An excitation pulse with an appropriate phase is used to rotate the magnetization so that its desired components can be measured. Finally, the spatial origin of the signals can be encoded by conventional means. As we explain later, we use also an additional known magnetic field Bext for finding out the amplitude of BJ. Assume that we have obtained images Ix, Iy, and Iz by applying the described sequence to image the magnetization components perpendicular to the orthogonal directions ex, ey, and ez, respectively, and B0 = B0ez. The complex-valued images are 8 < I x ðrÞ ∝ −bz ðrÞ þ iby ðrÞ I ðrÞ ∝ b ðrÞ−ib ðrÞ ; : I y ðrÞ ∝ b x ðrÞ þ ibz ðrÞ z x y

ð8Þ

where r is the spatial position, i the imaginary unit, and h i BJ ðrÞ ¼ BJ ðrÞ bx ðrÞex þ by ðrÞey þ bz ðrÞez ;

ð9Þ

where bx, by, and bz define the direction of BJ. The solution of the system in Eq. (8) gives the direction of BJ:

where ωJ= γBJ. Now, we can write Eq. (2) as h i3=2 f 2 ðt Þ þ 1 ωL ðt Þ ¼ −ωJ ≫1: Rðt Þ ¼ df =dt jdα=dt j

55

ð6Þ

h i 8 > bx ðrÞ ¼ Re I y ðrÞ þ Re I z ðrÞ =2 > > < h i by ðrÞ ¼ Im I x ðrÞ þ Im I z ðrÞ =2 : > h i > > : b ðrÞ ¼ − Re I ðrÞ þ Im I ðrÞ =2 z x y

ð10Þ

In principle, we could obtain the same information by acquiring only two images. However, with three measurements, all components of BJ are measured twice, making the method more robust. The additional information could also be useful in estimating the accuracy of the measurement. With the above procedure, we get an image of the direction of BJ. However, without additional information, it is impossible to determine the amplitude of BJ. To obtain the amplitude of BJ, we repeat the sequence with the addition of applying a known magnetic field Bext to the object together with BJ. This yields an image of the direction of BJ + Bext. In principle, we could now calculate the amplitude of BJ simply by requiring that the two directions differ by the change caused by Bext. However, such a direct approach is highly sensitive to errors in the direction estimates.

56

J.O. Nieminen et al. / Magnetic Resonance Imaging 32 (2014) 54–59

470 μT/m). The time to echo was 50 ms and the polarizing time was 1 s. The sample in our simulations, shown in Fig. 2, was a cubic volume of 20×20×20 voxels positioned at the origin. The voxel volume was (1 mm) 3. We assumed the relaxation times to be T1 = T2 = 200 ms throughout the volume. The current density J was assumed to be that of three long 6-mm-thick tubes with circular cross sections passing through the origin along the x, y, and z axes, each carrying a 40-mA current. The simulated signals were sampled with a receiver having uniform sensitivity along the x direction. For the forward calculations, we added small random deviations in the positions of the dipoles representing the voxels: the dipoles were deviated uniformly in all three Cartesian directions with the maximum deviation being 5% of the voxel side length. After the signals were simulated, white Gaussian noise with standard deviation σ was added to them to obtain a signal-to-noise ratio (SNR) of 10. We used the definition

SNR ¼ Fig. 1. Schematic diagram of the A-CDI sequence for imaging BJ. First, a polarizing pulse is applied together with BJ. Then, the polarizing field is switched off adiabatically, after which BJ is switched off and a homogeneous measurement field B0 is switched on. An excitation pulse is used to rotate the magnetization so that its desired components can be measured. After these steps, gradient fields (Gx, Gy, and Gz) are used to encode the spatial origin of the signals. For determining the amplitude of BJ, the sequence is repeated by applying a known magnetic field Bext together with BJ.

To obtain a more robust reconstruction, we also require one of Maxwell’s equations, ∇⋅B ¼ 0;

ð11Þ

to hold, where B = B1 = BJ in the first case and B = B2 = BJ + Bext in the second case. By a voxel-wise discretization of Eq. (11) presented in Appendix A, we can reformulate the problem of finding the field amplitudes to 0 @

A∇⋅B1 0 −AB1

1 0 1   0 a1 A @ A∇⋅B2 ¼ 0 A; a2 AB2 d21 0

ð12Þ

where A∇⋅Bi (i = 1,2) is a matrix operator implementing Eq. (11) in a discretized sense for Bi, ABi is a matrix that gives the Cartesian components of Bi, ai is a vector of the discretized magnetic field amplitudes for Bi, and d21 is a vector corresponding to B2− B1, i.e., a vector containing the components of Bext at each voxel. By solving Eq. (12) in the least-squares sense, we get a1 and a2, i.e., the amplitudes of BJ and BJ + Bext at each voxel. We can reduce the effect of noise in the estimate of BJ by combining the two results: from the solved B2 field we subtract Bext and average the results, obtaining BJ. Finally, we can calculate J using a discretized version of Eq. (1) (see Appendix B).

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u∫jsðt Þj2 dt t ∫σ 2 dt

¼

1 σ

sffiffiffiffiffiffiffi sT s ; Ns

ð13Þ

where s(t) is the signal as a function of time, s a vector containing the sampled data points, and Ns the number of acquired samples. In the simulations, the polarizing field was switched off adiabatically in 250 ms using a switch-off ramp designed using Eq. (6) with BJ = 5 μT. The current was applied between the start of the polarizing pulse and the end of the adiabatic ramp. The transverse magnetization after the adiabatic switching was solved for the three repetitions by constructing a reconstruction matrix for the two transverse magnetization components, see Ref. [35], and by solving the corresponding equation numerically. Then, the underlying direction of the magnetic field was solved from Eq. (10), and Eq. (12) was applied to find the amplitude of the magnetic field. Finally, J was calculated using a discretized version of Eq. (1). We also studied how the root-mean-square error of the reconstructed current density depends on the SNR by performing the reconstruction for different amounts of added noise. The reconstruction error was calculated as E¼

1 J0

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   1  ∑ ^J ðri Þ−Jðri Þ2 ; N i

ð14Þ

2.3. Simulation methods The sequence depicted in Fig. 1 was implemented to simulate ACDI. For simplicity, we decided to implement the spatial encoding by a 3D gradient-echo sequence. We assumed Bp = Bpex, B0 = B0ez, Bext = Bextey, and the gradient fields linearly varying and purely zdirectional. Frequency encoding was applied in the z direction with a gradient strength of 470 μT/m. The field amplitudes were Bp = 100 mT, B0 = 50 μT, and Bext = 1 μT. In the phase-encoding directions, 20 k-space lines were acquired with a 10-kHz sampling frequency (maximum phase-encoding-gradient strengths were

Fig. 2. Simulation geometry. The orange tubes represent the region of the simulated current flow. The simulation results are plotted along the planes indicated in green.

J.O. Nieminen et al. / Magnetic Resonance Imaging 32 (2014) 54–59

57

Fig. 3. Reconstructed current density at planes indicated in Fig. 2 using A-CDI. The dashed lines indicate the regions of the current flow. Without any reconstruction error, inside these regions, the current density normalized to the simulated true current density should be 1, which corresponds to 1415 A/m2, and outside them 0. Positive and negative values correspond to current densities along the positive and negative axes, respectively.

where the summation is performed over the N positions where the current density was reconstructed, J and ^J are the simulated and reconstructed current density, respectively, and J0 = 1415 A/m 2 is the current density along the tubes. 3. Results Fig. 3 shows the current density reconstructed using A-CDI. While the current density is reproduced well, some weak artefacts appear mainly because of measurement noise and imperfect adiabaticity. In Fig. 4, we show how the goodness of the adiabatic passage varies across the sample. For voxels where BJ N 1 μT the ramp is highly adiabatic. On the other hand, for most of the voxels in which BJ b 0.5 μT the deviation from perfect adiabaticity is over 20° and there are some voxels for which the passage is clearly non-adiabatic. The difference in the directional error between voxels with equal BJ is due to the different orientation of BJ in them. Fig. 5 illustrates how the reconstruction error of Eq. (14) depends on the SNR. As can be seen, the reconstruction quality improves with the SNR. However, some reconstruction error remains even without added noise. This is mainly because the adiabatic ramping of the fields is not perfect. Assuming perfectly adiabatic ramps and an error-free decoding of the magnetization orientation gave an error E = 0.06.

formalism, only the total field Brem + BJ can be measured. If Brem is of the same order as BJ or larger, it significantly affects the results. Thus, the method can be directly used only to image fields stronger than Brem. On the other hand, if Brem is static, the problem can also be dealt with by first imaging Brem without applying BJ. By subsequently measuring also the total field Brem + BJ, the field of interest, BJ, can be recovered by subtraction. Alternatively, the first measurement can be done by applying −BJ followed by a measurement with +BJ. Although we demonstrated an A-CDI sequence, which uses a known external magnetic field to obtain information about the strength of BJ, it is also possible to do the reconstruction without Bext.

4. Discussion In this study, we assumed that BJ was the only unknown field present in the experiments. In reality, however, a remnant field Brem exists even inside a magnetically shielded room. With the presented

Fig. 4. The angle between m and BJ at the end of the adiabatic passage for each voxel as a function of BJ in the voxels.

J.O. Nieminen et al. / Magnetic Resonance Imaging 32 (2014) 54–59

Reconstruction error

58

0.40 0.35 0.30 0.25 0.20 0

5

10

15

20

Signal-to-noise ratio

Fig. 5. Effect of the SNR on the reconstruction quality. The reconstruction error was calculated using Eq. (14). The dashed line indicates the reconstruction error without added noise.

In this case, one can use Eq. (1) at the object boundary to set a boundary condition for J. In practice, this can be done, e.g., by dividing the boundary into triangles and integrating Eq. (1) over their surfaces. The advantage of this approach is that the imaging time can be reduced by 50%. However, it is then necessary to know where the current enters and leaves the sample. In some cases, it may be advantageous to interrupt the adiabatic field ramps abruptly. For example, the use of a long field ramp would guarantee good adiabaticity throughout the field decay but result in a loss in the magnetization amplitude. By interrupting such a ramp before its natural end, the direction of the magnetization would still point in the direction of the total magnetic field before the interruption. By measuring the magnetization component perpendicular to the direction of the polarizing field and by applying the reconstruction algorithm, an image of BJ could be obtained. This type of sequence would resemble the use of low flip angles in traditional MRI. We note also that a range of adiabatic switch-off ramps different from Eq. (6) may be designed, e.g., by letting R(t) in Eq. (5) be a function of time. Because A-CDI is based on adiabatic field ramping, imaging weak current densities requires long polarizing-field ramps, during which the sample magnetization will decay. Thus, imaging current densities corresponding to magnetic fields below 1 μT will require ramps longer than ~200 ms. This means that highquality imaging of such weak current densities would require either a strong initial sample magnetization or an object with long relaxation time constants. On the other hand, because the

method is not susceptible to phase wrapping, A-CDI is wellsuited for imaging strong current densities (or magnetic fields) which can be turned off during the acquisition period. Thus, A-CDI may find applications, e.g., in mapping the magnetic field of coils which can be immersed in a liquid and driven with a suitable current. A comparison of A-CDI with zero-field-encoded current-density imaging (Z-CDI), presented in the accompanying paper [32], reveals several differences. First, the methods are based on completely different physical principles: A-CDI utilizes adiabatic magnetic field ramps for the encoding, whereas Z-CDI is based on spin precession in the unknown magnetic field without a B0 field present. Second, because Z-CDI uses a voxel-wise reconstruction of the current density, the noise characteristics of the Z-CDI reconstruction differ from those of A-CDI. 5. Conclusion We have presented a sequence to image magnetic fields and current densities using prepolarized MRI. The A-CDI method utilizes an adiabatic magnetic field ramp of the polarizing field to rotate the magnetization of an object parallel to the unknown magnetic field. It gives complete 3D magnetic field and current-density information without the need to physically rotate the object. The method may be useful, e.g., for conductivity imaging or mapping the magnetic field of coils. Acknowledgment The research leading to these results has received funding from the European Community’s Seventh Framework Programme (FP7/ 2007–2013) under grant agreement No. 200859, from the Instrumentarium Science Foundation, the Emil Aaltonen Foundation, the International Doctoral Programme in Biomedical Engineering and Medical Physics (iBioMEP), the Finnish Cultural Foundation, and the Academy of Finland. The authors thank Jukka Sarvas for useful discussions.

Appendix A This appendix outlines the Cartesian discretization of Eq. (11). We assume that B = B/B is known in a rectangular grid defined by the MR image, i.e., at the voxel centers. By joining the centers of a group of eight neighboring voxels, we can form a cube (Fig. 6A). Inside this cube, we can place two regular tetrahedra whose vertices

Fig. 6. An illustration of the notation used for the Cartesian discretization of Eq. (11). (A) The numbers show the positions of the voxel centers for a group of eight neighboring voxels in a uniform Cartesian grid. They are also used to refer to the voxels. In (B) and (C), we show two regular tetrahedra that are formed by joining the voxel centers. These tetrahedra define the integration surfaces of Eq. (A.1).

J.O. Nieminen et al. / Magnetic Resonance Imaging 32 (2014) 54–59

are the voxel centers (Fig. 6B and C). These tetrahedra define volumes over which we can integrate Eq. (11): ∮∂V Bb⋅dS ¼ 0:

ðA:1Þ

Assuming the magnetic field is a linear function of position on each face of a tetrahedron, we obtain for the tetrahedron of Fig. 6B ∮1764 Bb⋅dS ¼

1 ½B b ⋅ðv −v6 Þ  ðv7 −v6 Þ þ B 4 b4 ⋅ðv6 −v1 Þ 6 1 1 4  ðv7 −v1 Þ þ B 6 b6 ⋅ðv7 −v1 Þ  ðv4 −v1 Þ þ B 7 b7 ⋅ðv1 −v6 Þ  ðv4 −v6 Þ;

ðA:2Þ

where the subscripts 1, 4, 6, and 7 refer to the vertices of the tetrahedron and vi is the position of the ith vertex. By doing the same for all groups of eight neighboring voxels, we can construct the matrix A1764 whose rows correspond to the tetrahedra in the respective cubes and the columns to the voxels. When multiplied with a vector containing the voxel-wise magnetic field amplitudes, this matrix implements Eq. (A.1). We can also do equivalent calculations for the tetrahedron shown in Fig. 6C to obtain matrix A2583. Finally, we can join the matrices to obtain a more accurate discretization of Eq. (11). We used the combination  A∇⋅B ¼

 A1764 þ θA2583 ; A2583 þ θA1764

ðA:3Þ

with parameter θ = −0.5, because it produced better results than θ = 0. In order to balance the matrix in Eq. (12), we further multiplied A∇⋅B by a weighting factor of 5⋅10 6. This operation made the maximum element of (− AB1 AB2) and the maximum element of Eq. (A.3) approximately equal. Appendix B This appendix presents the Cartesian discretization of Eq. (1). We assume that BJ is known in a rectangular grid at the voxel centers. First, we integrate Eq. (1) over a surface: ∫S J⋅dS ¼

1 ∮ B ⋅dl: μ 0 ∂S J

ðB:1Þ

Let us consider a square surface limited by four neighboring voxels (their centers) at its corners. The positions of the voxels, in counter-clockwise order, are given by vi, i = 1,…,4. Assuming the magnetic field is a linear function on this surface, we obtain i 1 1 h ∮1234 BJ ⋅dl ¼ BJ ðv1 Þ−BJ ðv3 Þ ⋅ðv2 −v4 Þ μ0 2μ 0 h i þ BJ ðv2 Þ−BJ ðv4 Þ ⋅ðv3 −v1 Þ :

f

ðB:2Þ

g

Thus, the normal component of J at the center of the square can be estimated by dividing Eq. (B.2) by the area of the square. By considering separately all groups of four neighboring voxels, the three components of J can be calculated in the volume of interest. References [1] Joy MLG, Scott GC, Henkelman RM. In-vivo detection of applied electric currents by magnetic resonance imaging. Magn Reson Imaging 1989;7:89–94. [2] Pešikan P, Joy MLG, Scott GC, Henkelman RM. Two-dimensional current density imaging. IEEE Trans Instr Meas 1990;39:1048–53. [3] Scott GC, Joy MLG, Armstrong RL, Henkelman RM. Measurement of nonuniform current density by magnetic resonance. IEEE Trans Med Imaging 1991;10: 362–74. [4] Scott GC, Joy MLG, Armstrong RL, Henkelman RM. Electromagnetic considerations for RF current density imaging. IEEE Trans Med Imaging 1995;14:515–24.

59

[5] Scott GC, Joy MLG, Armstrong RL, Henkelman RM. Rotating frame RF current density imaging. Magn Reson Med 1995;33:355–69. [6] Ider YZ, Muftuler LT. Measurement of AC magnetic field distribution using magnetic resonance imaging. IEEE Trans Med Imaging 1997;16:617–22. [7] Mikac U, Demšar F, Beravs K, Serša I. Magnetic resonance imaging of alternating electric currents. Magn Reson Imaging 2001;19:845–56. [8] Wang D, DeMonte TP, Ma W, Joy MLG, Nachman AI. Multislice radio-frequency current density imaging. IEEE Trans Med Imaging 2009;28:1083–92. [9] Hämäläinen M, Hari R, Ilmoniemi RJ, Knuutila J, Lounasmaa OV. Magnetoencephalography—theory, instrumentation, and applications to noninvasive studies of the working human brain. Rev Mod Phys 1993;65:413–97. [10] Oh SH, Lee BI, Woo EJ, Lee SY, Cho MH, Kwon O, et al. Conductivity and current density image reconstruction using harmonic B z algorithm in magnetic resonance electrical impedance tomography. Phys Med Biol 2003;48: 3101–16. [11] Seo JK, Yoon J-R, Woo EJ, Kwon O. Reconstruction of conductivity and current density images using only one component of magnetic field measurements. IEEE Trans Biomed Eng 2003;50:1121–4. [12] Lee J-Y. A reconstruction formula and uniqueness of conductivity in MREIT using two internal current distributions. Inv Prob 2004;20:847–58. [13] Nachman A, Tamasan A, Timonov A. Conductivity imaging with a single measurement of boundary and interior data. Inv Prob 2007;23:2551–63. [14] Woo EJ, Seo JK. Magnetic resonance electrical impedance tomography (MREIT) for high-resolution conductivity imaging. Physiol Meas 2008;29:R1-26. [15] Hasanov KF, Ma AW, Nachman AI, Joy MLG. Current density impedance imaging. IEEE Trans Med Imaging 2004;27:1301–9. [16] Voigt T, Katscher U, Doessel O. Quantitative conductivity and permittivity imaging of the human brain using electric properties tomography. Magn Reson Med 2011;66:456–66. [17] Kim DH, Choi N, Gho SM, Shin J, Liu C. Simultaneous imaging of in vivo conductivity and susceptibility. Magn Reson Med, In Press, http: //dx.doi.org/10.1002/mrm.24759. [18] De Geeter N, Crevecoeur G, Dupré L. A numerical study on conductivity estimation of the human head in the low frequency domain using induced current MR phase imaging EIT with multiple gradients. IEEE Trans Magn 2013;49:5004–10. [19] Bodurka J, Jesmanowicz A, Hyde JS, Xu H, Estkowski L, Li S-J. Current-induced magnetic resonance phase imaging. J Magn Reson 1999;137:265–71. [20] Bodurka J, Bandettini PA. Toward direct mapping of neuronal activity: MRI detection of ultraweak, transient magnetic field changes. Magn Reson Med 2002;47:1052–8. [21] Kraus Jr RH, Espy MA, Volegov PL, Matlachov AN, Mosher JC, Urbaitis AV, et al. Toward SQUID-based direct measurement of neural currents by nuclear magnetic resonance. IEEE Trans Appl Supercond 2007;17:854–7. [22] Kraus Jr RH, Volegov P, Matlachov A, Espy M. Toward direct neural current imaging by resonant mechanisms at ultra-low field. Neuroimage 2008;39:310–7. [23] Burghoff M, Albrecht HH, Hartwig S, Hilschenz I, Körber R, Sander Thömmes T, et al. SQUID system for MEG and low field magnetic resonance. Metrol Meas Syst 2009;16:371–5. [24] Burghoff M, Hartwig S, Cassará A, Trahms L. Approaches to detect neuronal currents using a DC-mechanism by means of low field magnetic resonance. World Congress on Medical Physics and Biomedical Engineering, September 7 – 12; 2009. p. 31–3 Munich, Germany. [25] Burghoff M, Albrecht HH, Hartwig S, Hilschenz I, Körber R, Höfner N, et al. On the feasibility of neurocurrent imaging by low-field nuclear magnetic resonance. Appl Phys Lett 2010;96:233701. [26] Höfner N, Albrecht H-H, Cassará AM, Curio G, Hartwig S, Haueisen J, et al. Are brain currents detectable by means of low-field NMR? A phantom study. Magn Reson Imaging 2011;29:1365–73. [27] Belliveau J, Kennedy Jr DN, McKinstry RC, Buchbinder BR, Weisskoff RM, Cohen MS, et al. Functional mapping of the human visual cortex by magnetic resonance imaging. Science 1991;254:716–9. [28] McDermott R, Lee SK, ten Haken B, Trabesinger AH, Pines A, Clarke J. Microtesla MRI with a superconducting quantum interference device. Proc Natl Acad Sci USA 2004;101:7857–61. [29] Zotev VS, Matlachov AN, Volegov PL, Sandin HJ, Espy MA, Mosher JC, et al. Multichannel SQUID system for MEG and ultra-low-field MRI. IEEE Trans Appl Supercond 2007;17:839–42. [30] Vesanen PT, Nieminen JO, Zevenhoven KCJ, Dabek J, Parkkonen LT, Zhdanov AV, et al. Hybrid ultra-low-field MRI and MEG based on a commercial whole-head neuromagnetometer. Magn Reson Med 2013;69:1795–804. [31] Hilschenz I, Körber R, Scheer H-J, Fedele T, Albrecht H-H, Cassará AM, et al. Magnetic resonance imaging at frequencies below 1 kHz. Magn Reson Imaging 2013;31:171–7. [32] Vesanen PT, Nieminen JO, Zevenhoven KCJ, Hsu Y-C, Ilmoniemi RJ. Currentdensity imaging using ultra-low-field MRI with zero-field encoding. Magn Reson Imaging [E-pub ahead of print]. [33] Packard M, Varian R. Free nuclear induction in the Earth’s magnetic field. Phys Rev 1954;93:941. [34] Abragam A. The Principles of Nuclear Magnetism. Oxford, UK: Clarendon Press; 1961. [35] Nieminen JO, Ilmoniemi RJ. Solving the problem of concomitant gradients in ultra-low-field MRI. J Magn Reson 2010;207:213–9.