Research Articles

Download PDF (831.89 KB)

TOMOGRAPHY, December 2015, Volume 1, Issue 2: 125-135
DOI: 10.18383/j.tom.2015.00142

Feasibility of Imaging Tissue Electrical Conductivity by Switching Field Gradients with MRI

Eric Gibbs1, Chunlei Liu2

1Brain Imaging and Analysis Center, and2Department of Radiology, Duke University School of Medicine, Durham, NC


Tissue conductivity is a biophysical marker of tissue structure and physiology. Present methods of measuring tissue conductivity are limited. Electrical impedance tomography and magnetic resonance electrical impedance tomography rely on passing an external current through the object being imaged, which prevents its use in most human imaging. More recently, tissue conductivity has been determined noninvasively from measurements of the radiofrequency (RF) field used in magnetic resonance imaging (MRI). This technique is promising, but conductivity at higher frequencies is less sensitive to tissue structure. Measuring tissue conductivity noninvasively at low frequencies remains elusive. It has been proposed that eddy currents generated during the rise and decay of gradient pulses could act as a current source to map low-frequency conductivity. This work centers on a gradient echo pulse sequence that uses large gradients before excitation to create eddy currents. The electric and magnetic fields during a gradient pulse are simulated by a finite-difference time-domain simulation. The sequence is also tested with a phantom and animal MRI scanner equipped with gradients of high gradient strengths and slew rates. The simulation demonstrates that eddy currents in materials with a conductivity similar to biological tissue decay with a half-life on the order of nanoseconds, and any eddy currents generated before excitation decay completely before influencing the RF signal. Gradient-induced eddy currents can influence phase accumulation after excitation, but the effect is too small to image. The animal scanner images show no measurable phase accumulation. Measuring low-frequency conductivity by gradient-induced eddy currents is presently unfeasible.


Electrical conductivity of biological tissue, the measure of how easily electrical current flows through tissue, is relevant to both biological and clinical studies. Tissue conductivity has been shown to be heterogeneous and anisotropic throughout the body, providing rich information related to tissue molecular structure and composition (1, 2). For example, the anisotropic ratio of conductivity when measured parallel/perpendicular to white matter tracts is 10:1 (3). Tissue conductivity is also of clinical interest for diagnosing tumor and diseased tissue (4, 5) and is required for accurate source location in electroencephalogram/magnetoencephalogram measurements (6).

Multiple biochemical factors influence bulk conductivity properties, including tissue structure, composition, ionic content, and temperature, each of which is reviewed in detail by Pethig (7). Tissue conductivity has a strong frequency dependence because different molecular and structural components dominate at different frequencies. In brief, at lower frequencies (<100 MHz), cell membrane structures contribute significantly to overall conductivity as mobile ions encounter membrane barriers; at higher frequencies (>100 MHz), ionic concentration becomes the dominant factor, and dielectric properties are independent of membrane structure, although the contents and thus dielectric properties of fluids differ both intra- and extracellularly.

One method for measuring low-frequency conductivity is electrical impedance tomography (EIT) (8). EIT relies on passing current from an external source through the object being imaged. Voltage electrode readings on the surface of the object are used to map the internal tissue conductivity. Determining a material's internal conductivity using surface voltage measurements is an ill-posed inverse problem and gives unreliable measurements in deeper tissue. Woo et al. (9) introduced magnetic resonance EIT (MREIT), a technique that uses magnetic resonance imaging (MRI) to measure magnetic field perturbations caused by injected current and—from these perturbations—calculate conductivity maps. Various techniques and algorithms have improved this process (10, 11), and MREIT has successfully created conductivity maps in phantoms (9), animals (12), and the human thigh (13), but the reliance on injected current limits MREIT's practical applications.

Alternatively, tissue conductivity can be determined noninvasively via B1 field mapping, albeit at the Larmor frequency, which is 128 MHz for 3T scanners (14). The B1 field is altered by the conductivity of the object being imaged and under certain assumptions can be used to determine the underlying dielectric properties. There has been rapid improvement in the reconstruction of conductivity maps from B1 maps (15, 16), but the technique is still resolution-limited. There is also the fundamental limitation that conductivity maps at the Larmor frequency are likely not sensitive to tissue microstructure.

To our knowledge, there is currently no effective method for measuring low-frequency tissue conductivity without injecting current. It was proposed that eddy currents from pulsed gradients could be used to map conductivity (1719). During an MRI experiment, the changing magnetic field of a gradient pulse induces an electric field that in conductive material drives eddy currents. The phase of the complex MRI signal is equal to γBzdt, where γ is the gyromagnetic ratio and Bz is the z-component magnetic field. The magnetic field created by eddy currents in gradient coils alters the overall magnetic field and requires compensating currents to avoid image artifacts (20). The magnetic field that arises from eddy currents in tissue also contributes to the overall magnetization phase but is less well studied. If magnetic field contributions from eddy currents to the overall phase signal can be isolated and measured, algorithms similar to those used in MREIT could be used to determine tissue conductivity.

Large crusher gradients maximize tissue eddy currents. However, the effect of the crusher rise and decay will affect the signal phase in an equal and opposite fashion, leaving no net effect. The crusher itself also has a much larger impact on the phase than the eddy currents. These problems can be mitigated by lengthening the crusher and placing it immediately before the radiofrequency (RF) pulse as shown in Figure 1. This arrangement allows the effects of the eddy currents from the crusher rise and the crusher itself to occur before excitation, but the effects of eddy currents from the crusher decay can potentially extend past excitation and be recorded in the phase signal.

Figure 1.

The proposed pulse sequence to map conductivity using eddy currents generated by an adjustable large crusher gradient. This gradient is positioned immediately before the excitation RF pulse with a flip angle of α. Eddy currents are generated during the rise and descent of the EC crusher. Using a 7T animal scanner, δ was 40 μs and reached a peak amplitude of 200 mT/m, giving a slew rate of 5000 T/m·s; Δ was 5 ms and τ was 100 μs. The intent is that the eddy currents created by the pulse rise will have decayed when the pulse decay occurs. If the eddy currents from the pulse decay last beyond the time period τ they will affect the MRI signal.


This strategy requires not only that the magnetic field to be strong enough to be detectable in MRI phase images but also that the eddy currents in tissue last beyond the RF excitation. In this work, we examine the feasibility of both of these points.



Maxwell's equations give the theory behind the use of secondary magnetic fields from eddy currents for conductivity mapping. Faraday's law of induction relates changes in the magnetic field, such as those during a switching gradient pulse, to an induced electric field as follows:

where B and E are the magnetic and electric fields, respectively. The curl operator dictates that the induced electric field circles in a plane perpendicular to the changing magnetic field. This electric field exerts a force on any charge in the material, and if the material is conductive this will create an eddy current according to Ohm's law: Je=σE, where Je is the current density of eddy currents and σ the material's conductivity. The secondary field generated by the eddy currents is described by Ampere's law. Assuming a linear electric and magnetic susceptibility Ampere's law states that
Et+Jeϵ= ×Bϵμ
where ϵ is the permittivity and μ is the permeability of the material. Either a changing electric field or eddy currents will create a magnetic field that circles perpendicular to the source. Because of the negative sign in Faraday's law, the net result is a secondary magnetic field that opposes the initial changes in B. The secondary field affects Bt in Faraday's law, and the coupled equations together describe the total electromagnetic field evolution. If ϵ and μ are known, the 2 equations allow one to express σ in terms of B. If magnetic field perturbations from eddy currents are detectable by MRI, algorithms similar to those used in MREIT could be used to map conductivity from the MRI phase signal.


A finite-difference time-domain (FDTD) simulation written in MATLAB version R2014a (MathWorks, Natick MA) was used to simulate the electric and magnetic fields during a gradient pulse. Each iteration of the FDTD simulation uses the discrete form of Maxwell's equations to advance the electromagnetic fields by a finite time dt, which allows dynamic tracking of the field strength and shape (21, 22). The blue arrows in the flow chart of Figure 2 outline the standard electromagnetic wave FDTD algorithm. In the first step, Ampere's law is discretized such that a component of the electric field at time point tn+1 can be written in terms of the discrete curl operator acting on the past magnetic field and the past electric field—both at time point tn. After updating each electric field component, Faraday's law is then used to update each magnetic field component at tn+1 in a similar manner but using the freshly determined electric field. The simulation then cycles back to updating the electric field for the next time point. As the curl acts nonlocally, each component of the electromagnetic fields has its own grid that covers the simulation volume but is slightly offset from the other components. This approach has been used in the past to calculate the strength of eddy currents in tissue, although the impact of the eddy currents on the magnetic field was not reported (23).

Figure 2.

A representation of the cylindrical phantom and gradient coils is shown on the top left, and an axial profile of the phantom is shown on the right. The cylinder itself (red portion) has conductivity σp, and the 4 inserts (blue regions) have conductivity that varied from σp by the shown amount (−0.5, +1, +3, and +2 S/m). Simulations were done with σp set at 1, 10, 100, and 1,000 S/m. To simulate the electromagnetic fields without conductive material, the entire phantom was uniformly set to the background conductivity of 0.001 S/m. The bottom flow chart gives the equations to create the simulation. The first simulation follows the blue arrows, whereas the second continues with the red arrows. In regard to the second simulation equations, ρfree is the free charge that accumulates at material interfaces, r and r′ are spatial vectors, V is the voltage created by ρfree, and Eρ is the electric field created by V. Details of discretizing these equations can be found in references 21 and 22.


The simulation volume modeled a magnetic bore partially loaded with dielectric material. The bore walls were modelled as perfect electric conductors and the open ends as perfectly matched layers (24). Differences between the simulation with and without a dielectric phantom were attributed to eddy currents. In addition to the standard FDTD simulation that assesses the effects of bulk conductivity, a second simulation that accounts for dielectric interfaces was done as well.

The first simulation was done in a cylindrical coordinate system to more accurately portray the current density as explained in the next paragraph. The grid size in r, ϕ, and z was 31 × 28 × 82, with a 1-cm resolution in r and z and π/14 rad in ϕ. For the FDTD simulation to be stable, the fields cannot advance faster than the rate of propagation determined by the underlying equations, which in this case means dtΣdxi2c=4.4ps, where c is the speed of light and dxi is the distance between simulation grid points in the ith direction. The simulated phantom was a cylinder with a radius of 14 cm and height of 36 cm centered at the middle of the bore as shown in Figure 2. The phantom had a primary conductivity σp, with variations in the primary conductivity that were on the order of 1 S/m, as shown in the axial phantom slice. Four simulations were done with σp = 1, 10, 100, and 1000 S/m. To ease numerical discontinuities, regions outside the phantom were set to a background conductivity of 0.001 S/m. This background conductivity reduces noise by allowing nonphysical transients to decay without significantly altering the simulation values (24). For reference, σ is between 0.4 and 1.4 S/m in biological tissue and 5.96·107 in copper; ϵr was set at 80 throughout the phantom and 1 outside the phantom for all bulk conductivity simulations. Dielectric properties at interfaces were handled as described by Jewett et al. (25).

To stimulate the phantom with a gradient pulse, 2 rings (18-cm inner diameter, 19-cm outer diameter, separated by 33 cm, and centered about the isocenter) were designated as coil regions. The pulse was created by incrementally increasing the electric field in the coil region in the Eϕ direction until reaching a constant final value. By directly increasing Eϕ as opposed to modeling a voltage source, one bypasses the coil resistance and eliminates coil-induced eddy currents, whose magnetic field would mask those of tissue-induced eddy currents. Cylindrical coordinates were chosen to simulate the current distribution more realistically. In reality, current flows without any charge accumulation, or equivalently, the divergence of the current density is uniformly 0. Per Ohm's and Gauss's laws, a divergence in the current density creates a net charge. On a finite grid, cylindrical coordinates allow for a divergence-free current density, whereas Cartesian coordinates do not. The pulse rise was modeled as the first 8 terms of the Fourier series of a trapezoidal pulse and set constant when the first maxima was achieved. The result was a linear gradient along the bore axis, and the pulse parameters were empirically adjusted so that the rise time was 1 μs and the final gradient amplitude was 150 μT/m, giving a slew rate of 150 T/m·s. This large gradient, although not available on whole-body human scanners, is typical of small-bore animal scanners with small insert gradients as used in our laboratory. Eddy currents are driven by Bt and not by the magnitude of B; thus, for a symmetric pulse, the eddy currents caused by the pulse rise or decay differ only in sign. Because the FDTD simulation requires such a short dt, only the pulse rise was simulated.

Initial numerical transients at sharp dielectric property interfaces proved to be unstable when reflected off the conductive material. The simulations stabilized after adding a term to the discretized time derivative such that t1dt+δ. The effect of δ was to dampen dramatic initial swings as the fields grew from 0 magnitude. In simulations without conductive material, including δ increased field values by about 50%, but when the field was normalized by the maximum value, simulations with and without δ differed by less than 1%, indicating the field shape was conserved. Simulations with the same value for δ converged to the same steady state after the decay of eddy currents.

It has also been proposed that the magnetic field is influenced by a charge that accumulates at dielectric interfaces during gradient pulses (26). The FDTD simulation described previously was extended to include the effects of accumulated charges. This is shown with red arrows in Figure 2. At each time point after updating the electromagnetic fields, the electric field was used to determine the accumulated charge density.

The equations for updating charge density are derived using the conservation of charge equationρfreet=J, where ρfree is the free-charge density. By applying Ohm's law, J=σE, allowing σ and ϵ to be inhomogeneous and assuming linear material, one can expand the conservation of charge to

Equation 3 implies that existing charge decays with a characteristic decay time of ϵ/σ, whereas an electric field drives charge accumulation in regions where dielectric properties vary spatially. By discretizing this equation, one can advance the charge density from the calculated electric field.

The charge creates a voltage, which is equal to

where r and r′ are spatial vectors, and the integration is over the entire simulation volume, although only charge-accumulating interfaces have a nonzero charge density. This integral can be represented as a convolution and efficiently implemented using the fast Fourier transform (FFT) algorithm. To correctly propagate the voltage created from the accumulated charge, it is necessary to use the time-retarded potential. This can be done by solving for the voltage for the most updated charge distribution and saving the calculated voltage to be used n time points later if r is n grid points from a source at r′. The gradient of this voltage is an additional electric field contribution from an accumulated charge that is added to the electric field created by changes in the magnetic field. The simulation then restarts the cycle, updating the fields as described previously.

A smaller rectangular grid of 32 × 32 × 42 was used with a 1-cm resolution in x, y, and z, and dt = 19 ps. The grid size was reduced because the added complexity slowed the simulation and rectilinear coordinates were used so that convolutions could be done using the FFT. Regions where charges accumulated as a result of a nonzero divergence of the current density were set to 0. The phantom was adjusted to be a 16-cm cube with homogenous conductivity, which is sufficient for examining a charge that accumulates on the surface. The simulation is inherently unstable for reasons explained in the next paragraph, but under certain conditions the simulation does not diverge quickly. The driving term in Equation 3 is proportional to the gradient of σ and ϵ, and the sharp transitions from air to physiological values of σ and ϵ created large values of ρfreet that destabilized the simulation; consequently, the simulation diverged too quickly to be useful. To explore the relationship between σ, ϵ, and charge accumulation, simulations were done with lower values of σ and ϵ, specifically, σ=0.06,0.12, and 0.24 S/m, with ϵr=8 and ϵr=4,6, and 8, with σ = 0.12 S/m.

The second simulation is unstable because of net charge growth. In reality, total charge is conserved, and positive charge created at one point is balanced by an equal negative charge elsewhere. During the simulation, the total charge is generally conserved, but any net imbalance in charge creates a self-feeding cycle that causes the charge density and hence electromagnetic fields to diverge. Unfortunately, because of the finite accuracy of numerical computing, this instability is an eventuality, but the simulation remains useful for a reasonable amount of time if carefully implemented. Unbalanced charges can be reduced, and the simulation diverges only after millions of iterations if each electric field component is interpolated onto the same grid, aliasing artifacts from the FFT are windowed, and the charge density calculated at each iteration is averaged with recent past values to prevent large fluctuations.

MRI Experiments

MRI experiments using the proposed sequence were conducted on a 7T small-bore animal system controlled by Agilent software version vnmr4 (Santa Clara, CA). The scanner is equipped with a gradient system that provides a maximum gradient strength of 620 mT/m and a maximum slew rate of 6200 T/m·s. A phantom was made of a plastic tube about 1.5 cm in diameter and 6.5 cm in length. The tube was filled with a solution whose Gd concentration was equivalent to a 5% ProHance solution (Bracco Imaging, Township, NJ). Inside the tube were 2 sealed straws with 1% and 2% NaCl by mass added to the Gd solution. For reference, the conductivity of gray matter tissue can be mimicked using an agar gel of 0.15% NaCl (27). The straws were held firmly against the tube walls to prevent vibrational motion artifacts.

The tube was placed in the magnetic bore perpendicular to the main magnetic field of the 7T animal scanner. The readout direction was aligned with the tube axis. A multigradient echo scan was modified to include an eddy current (EC)-generating crusher gradient. Scan parameters were TR = 82.2 ms, TE = 6.0 ms, TE2 = 10.0 ms, and α = 20°, with a 62.5-μm isotropic resolution. The EC crusher had a maximal amplitude of 200 mT/m with a rise and decay time of 40 μs, achieving a slew rate of 5000 T/m·s. RF excitation was done using a square pulse that lasted 100 μs. The gradient was placed on the readout axis, and 2 scans were done that differed only in the polarity of the EC crusher. For comparison, the scan was done without the EC crusher as well.

To isolate the effects of eddy currents, the phase signal from the positive EC crusher was subtracted from that of the negative EC crusher (28). The phase signal was then unwrapped using a Laplacian-based phase unwrapping technique (29). Next, eddy currents in gradient coils created a magnetic field that varied linearly in space, resulting in a phase term that varied in the same manner. This phase term was fit and removed from the phase data. Finally, phase accumulation resulting from magnetic sources outside the field of view (FOV) was removed via the VSHARP algorithm, which uses the spherical mean value theorem to separate phase contributions from sources inside and outside the FOV (30). Tissue eddy currents are within the FOV, and the phase originating from sources outside the FOV was considered background phase and removed. The principle behind VSHARP's ability to separate these 2 terms is that the phase originating from sources outside the FOV satisfies Laplace's equation, whereas those originating from within the FOV do not.


Simulation of the Temporal Characteristics of the EC-Induced Magnetic Field

Each simulation with a dielectric phantom produces a magnetic field BD. The magnetic field from the simulation with no dielectric is BND. The difference between the 2 represents the magnetic fields generated by eddy currents, denoted as BE. The time course of each simulation at a fixed point in the phantom is shown in Figure 3A, and the BE for each dielectric phantom is shown in Figure 3B. Maxwell's equations dictate that BE opposes changes in the magnetic field. This opposition creates a delay in BD compared with BND. This delay increases in magnitude and duration as conductivity increases. For materials with a large conductivity and hence long-lasting eddy currents, BE essentially opposes all changes and is nearly the opposite of the field driving the eddy currents, that is, BE ≈ −BND. At low conductivities, the eddy currents decay too quickly to accumulate over time, and BEBND/t. The magnitude of BE depends on geometry as explained in the following section, but the maximum value of BE at the point shown in Figure 3 was 0.05, 0.5, 3, and 6.5 μT for σp=1,10, 100, and 1000 S/m, respectively. Note that BE initially almost completely opposes the maximum of BND in the 1000-S/m conductivity phantom.

Figure 3.

Time courses of simulated magnetic fields under various conductivity values. BecauseBt differs only in sign between the pulse rise and decay of a symmetric pulse, only the pulse rise is simulated to reduce simulation time. (A) Time courses for the overall magnetic field. (B) Time courses for the secondary magnetic fields caused by eddy currents in the bulk eddy current simulation. The measurements were made at a point 4.5 cm above the isocenter on the bore axis. The eddy currents caused a delay in the growth of the overall magnetic field. Insets show a boxed region magnified. Different colors of the time course represent different conductivity values.


To better understand the temporal evolution of BE, it is useful to think of BE as the output of a linear time invariant signal system and the changing current in gradient coil current as the input. Eddy currents have been successfully modeled by exponential decay (20). Because the eddy currents are small, it is reasonable to assume that the immediate surrounding material forms a single current loop that perturbs the magnetic field, from which one expects the current and hence magnetic field to decay monoexponentially. To test this premise, a monoexponential function was used as the impulse response function to the gradient coil current. An estimated eddy current field BE=ItAebt is fit to BE for each value of σp. Reasonable fits are achieved for σp= 1, 10, and 100 but not for 1000 S/m. As shown in Figure 4, the fit is best for σp = 1 S/m. This supports our premise that eddy currents act as a single current loop at low conductivity. Because an excellent fit was achieved for the physiologically relevant 1 S/m phantom, more complicated models were not explored. The impulse response half-life is defined by τ=ln(2)b and τ = 2, 11, and 119 ns for σp = 1, 10, and 100 S/m, respectively, and although a good fit is not achieved, σp = 1000 S/m seems to decay over about 10 μs. The rate of decay varies depending on where the field is measured, but only eddy currents produced from the phantom with σp = 1000 S/m lasted longer than a few microseconds. Because the eddy currents decay before RF excitation, there is no phase accumulation, but if the phase were acquired during this time, the net phase accumulation given by ϕ=γBz,Edt, where Bz,E is the z-component of the eddy current-induced magnetic field, would be 1.7·10−4, 1.6·10−3, 2.6·10−2, and 1.2 rad for σp= 1, 10, 100, and 1000 S/m, respectively.

Figure 4.

The magnetic field resulting from eddy currents can be modeled using a linear system. Rapidly decaying eddy currents are assumed to behave as a single current loop from immediately neighboring areas. The time derivative of the coil current is the input function and a monoexponential Aebt, which represents current decay in a single loop, is the impulse response function. After optimizing A and b, the predicted output closely follows the time course for BE when σp= 1 S/m. The half-life of the mono exponential is τ=ln(2)b=2ns. The fit is less acceptable as conductivity rises and is not reasonable for σp= 1000 S/m.


Simulation of the Spatial Characteristics of the Eddy Current-Induced Magnetic Field

The largest determinant of BE is conductivity, but geometry plays a role as well. The magnitude of BE at a given point depends on the value of BND at that point, but BE/BND is largest in the areas most embedded in conductive material, as shown in Figure 5A. To determine whether variations in conductivity can be resolved with the overall magnetic field, Figure 5B, C shows BD along the inner circumference of the phantom at a depth of 7 cm. The phantom depth is constant over this line, but conductivity varies from σp by −0.5 to 3 S/m. Different conductivity regions are initially distinct, with sharp variations within 1% of the field amplitude. However, as the amplitude of BD rises to values detectable by MRI, the boundaries between conductive regions are smoothed, and the percentage of variation drops. When BD reaches 1 μT, BD varies by only 1.5 nT in regions whose conductivity varies by 3.5 S/m. The ability to resolve variations in BD for a given magnitude BD is independent of the underlying conductivity, although BD rises more slowly in highly conductive material, increasing the duration of variations in BD.

Figure 5.

Magnetic fields vary as a function of both spatial position and conductivity value. (A) Ratio of the eddy current-induced magnetic fields to the fields without material along the z-axis of the cylinder phantom. The ratio is greatest in the center of the phantom. (B, C) Overall field variations as a function of ϕ at z = 4.5 cm and r = 6 cm, as indicated on the phantom schematic. The variations are distinct when the amplitude of the overall magnetic field is very small in (B) (corresponding to σp = 1000 S/m and t = 0.13 microseconds), but they become blurred, and a smaller fraction of the overall field grows in magnitude in (C) (corresponding to σp = 1000 S/m and t = 1.35 microseconds. When the overall field is 1 μT in (C), the variation between regions of 4 S/m is about 1 nT.


Along with the magnetic fields, the simulation also reported the electric fields and hence the current density. The current density magnitude was largest on the edges of the phantom close to the coil and decreased deeper in the phantom. Regions of higher conductivity had higher current densities as well. For the σp = 1 S/m phantom at the peak amplitude of BE, the maximal current density was about 3.76 A/m2, with an average magnitude of 0.87 A/m2. This is comparable to the eddy currents reported by Liu et al. in their FDTD simulation of the human body (21), although the results here are slightly larger because the coils were simulated close to the phantom and because 1 S/m is more conductive than most human tissue.

Simulation of Surface Charge

The simulation of charge accumulation at interfaces tracks the charge density at dielectric interfaces. The circular electric field induced by the coils creates a charge density gradient that goes from positive to negative across the surfaces of the phantom. Although the charge seems to be balanced, a net charge accumulates, and the simulation diverges during the pulse plateau. The maximum charge density is tracked through the pulse rise. This curve is fit using It as input and a monoexponential impulse response function as described previously. The maximal charge density and fit are shown in Figure 6, and the impulse response has a half-life of 4.9 ns when ϵr=8 and σ = 1.2 S/m. Increasing the conductivity increases the rate of response, but as shown in Figure 6B, the accumulated charge density does not change. Conversely, as shown in Figure 6C, materials with higher permittivity have a higher charge density that is reached in roughly the same amount time for each material. The results suggest the effects of accumulated charge dissipate faster than those of eddy currents.

Figure 6.

Maximal charge density from the charge interface simulation. (A) The charge essentially follows the slew rate with a delay of about 10 nanoseconds. This delay is more clearly seen in (B) and (C), which show the first 50 nanosecond of (A) (boxed region) for different values of σ and ϵ. (B) The charge equilibrates more quickly in more conductive materials but to a constant value. (C) Larger permittivity values increase the induced charge density.


MRI Experiments

Images from the 7T animal scanner are presented in Figure 7. The magnitude image showed no substantial difference between the straws and background water. Despite the straw ends being snugly fit to plastic anchors against the tube wall to reduce movements, there were faint artifacts that resulted from the vibrations of the straws that were even fainter in phase data. Figure 7 shows the phase image from the third echo at each postprocessing step. Low-signal regions are masked out and set in black. Once the effects of the linear gradient and background field were removed, there was no significant phase contrast between the upper tube of 1% NaCl, the lower tube of 2% NaCl, and the background, which has no NaCl. There was also no phase accumulation when the EC crusher was played on either phase-encoding axis.

Figure 7.

Phase maps of a tube and plastic straw phantom from the proposed pulse sequence. The phantom consists of 2 straws, 1 filled with 1% NaCl and the other with 2% NaCl. The background is filled with distilled water. Two images were acquired with opposing EC crusher polarities. The phase images of the 2 acquisitions were subtracted. In principle, this cancels phase contributions unrelated to eddy currents. The grayscale bars are in units of radians. There is a substantial phase contrast across the phantom in the subtracted phase and the linear phase-filtered phase. However, when the VSHARP algorithm is used to remove background sources, there is no significant phase contrast. Regions of low signal were masked out for filtering purposes.



Using quasi-static and FDTD approaches, we investigated the feasibility of imaging tissue conductivity at a low or near direct current frequency using eddy currents induced by the switching of field gradients. The resulting simulated values are comparable to past work (19, 23, 31). In addition, we also calculated the effect of charge accumulation at dielectric interfaces and found that its effect decays even faster than eddy currents. This section will focus on motivating the spatiotemporal evolution observed in the simulations and discuss the implications for imaging experiments.

As required by Maxwell's equations, all electromagnetic changes are propagated as waves, and the electromagnetic field of a gradient pulse can be represented in free space as a packet of plane waves ei(k·rωt), where k is the spatial wave vector. In material, ignoring the effects of accumulated charge, Maxwell's equations still allow for a plane wave solution, except the wave vector k is complex. This leads to the wave being attenuated as it penetrates the material. To simplify the analysis, consider only waves travelling perpendicular to the surface; in this case, the exponential attenuation coefficient is (32)

Note that the decay is a function of space that explains the relationship between BE and phantom depth. κ increases monotonically from 0 as a function of both frequency and conductivity, creating a low-pass filter that is more selective as σ increases.

The low-pass filter accounts for the temporal delay of the magnetic fields in conductive material. The distance at which the plane wave is half of its original amplitude is given by x1/2=ln(2)κ. Values of x1/2 can be tabulated as a function of ω. Figure 8 is a 2-dimensional schematic that shows half-decay lines as a function of phantom depth and frequency. Note that at larger conductivities most frequencies hardly penetrate the material, giving rise to skin effects. In relation to the simulation, the largest frequency component in the pulse rise is 500 kHz. At this frequency, for the 1000 S/m phantom, x1/2 is only a few centimeters (meaning it is well attenuated), whereas the 10 S/m phantom needs to have a depth of about 40 cm to have a similar attenuation. To translate this to imaging experiments, the experiments shown in Figure 7 have an excitation pulse of 100 μs. Inverted, this gives a frequency of 10 kHz. To attenuate this component by half, a tissue with a conductivity of 1 S/m needs to extend to a depth of 2.5 m. The effect to which the total field is modified depends on the spectra of the pulse. Gradients with rapid changes will have a broader spectra, and attenuation of higher frequencies will have a greater impact on the overall field evolution.

Figure 8.

Two-dimensional plots showing electromagnetic wave decay as a function of frequency and material depth. Lines show the frequencies and depths at which a wave is 50% attenuated for a given conductivity. The larger the conductivity, the more frequencies are attenuated at a given depth. In conductive material, the attenuation of high frequency components creates a low-pass filter that accounts for the temporal delay seen in simulations. Inset shows the boxed magnified.


The analysis described previously is simplified in that it assumes a flat surface and direct incidence. Accurate analytical analysis is impractical because it requires complete knowledge of the frequency composition of the fields and proper accounting of oblique incidence and internal reflection. As a surrogate, the impulse response function approach used in this work allowed one to extrapolate field evolution. Each point has its own impulse response that is a function of both the phantom depth and surrounding material but the magnetic field evolution at each point qualitative resembles the field evolution shown in Figure 3. The monoexponential model of eddy currents worked best for the lowest-conductivity phantom. The exponential has a half-life of about 2 ns. Eddy currents in this phantom have undergone 50,000 half-lives by the time the RF pulse has finished. This decay is so rapid that one cannot expect to have any tissue-related eddy currents at the time of acquisition. Even if the eddy currents occurred during acquisition, however, the accumulated phase at 1 S/m (1.7·10−4 rad) is currently too small to be resolved in phase MRIs.

The results of the second simulation also demonstrate that an induced charge dissipates even faster than the bulk eddy currents, and the rate of dissipation increases as a function of conductivity. This seems paradoxical because eddy currents last longer at higher conductivities but is sensible when considering the actual physical dynamics, where a freely flowing charge can rapidly balance any net charge and that current can flow in electrically neutral material.

Other strategies have been proposed that use gradients after RF excitation and a spin echo to generate the phase (18, 31). If the gradient is on continuously during the spin echo across the refocusing pulse, the phase accumulated during the rise and decay do not cancel but add in this case as a result of the spin echo reversal. This allows for net phase accumulation from tissue-induced eddy currents after RF excitation. Such images, as with our unfiltered images, have a phase accumulation that is much higher than the expected accumulation from tissue-induced eddy currents (19, 31). Mandija et al. (19) explained this phase accumulation as a result of geometric distortions resulting from coil eddy currents that shift the positive and negative gradient polarity images from one to another. This is observed as well in our images because the eddy currents in the gradient coils endure until image acquisition. The effects of the RF-induced phase are not perfectly cancelled when the images are subtracted as a result of this shift. This results in RF leakage in the phase signal, which is influenced by the material conductivity at the Larmor frequency. The theoretical phase acquisition resulting from eddy currents in tissue during a slice-selection spin-echo experiment is reported to be between 10−5 and 10−3 rad, which is similar to the phase accumulated in our work. These small values are completely obscured by RF leakage and noise (19, 31). Even if the phase resulting from the eddy currents could be isolated, this work suggests that resolving variations in conductivity similar to those found in tissue would be another challenge.


We have presented simulations that explore the time evolution of eddy currents induced in tissue as a result of large crusher gradients. These results and methods may be useful for materials of intermediate conductivity, but at higher conductivities the FDTD simulation must take surface current into account because the electromagnetic fields are attenuated just beneath the surface (33). The proposal to use eddy currents before excitation is completely unfeasible because of the rapid decay of tissue eddy currents. Using the simulation half-life of 2 ns, the amplitude of the eddy currents undergoes 50,000 half-lives before RF excitation in materials comparable in size and electrical properties to that of a human head. The phase that accumulates from a switching gradient pulse after RF excitation is small, but given a 100-fold increase it could be detected. This is still a large hurdle to cross. The spin echo approach allows one to acquire phase from tissue-induced eddy currents after RF excitation, avoiding the problem of rapid decay, but the phase that accumulates from a switching gradient is still about a 100 times smaller than the expected noise and RF leakage. Although not yet feasible, low-frequency conductivity mapping remains an interesting opportunity that merits an eventual solution.


[1] Abbreviations:


Echo time


eddy current


electrical impedance tomography


fast Fourier transform


field of view


finite-difference time-domain


magnetic resonance imaging


magnetic resonance EIT




repetition time


variable kernel sophisticated harmonic artifact reduction for phase data


We would like to acknowledge the help and support of former lab member Ioannis Argyridis as well as Gary Cofer and John Nouls at the Center for In Vivo Microscopy at Duke University. We also thank Craig Henriquez for his discussions.


    Gabriel C, Gabriel S, Corthout E. The dielectric properties of biological tissues: I. literature survey. Phys Med Biol. 1996;41(11):2231.
    Geddes L, Baker L. The specific resistance of biological material—a compendium of data for the biomedical engineer and physiologist. Med Biol Eng. 1967;5(3):271–293.
    Nicholson PW. Specific impedance of cerebral white matter. Exper Neurol. 1965;13(4):386–401.
    Hancu I, Roberts JC, Bulumulla S, Lee SK. On conductivity, permittivity, apparent diffusion coefficient, and their usefulness as cancer markers at MRI frequencies. Magn Res Med. 2015;73(5):2025–2029.
    McEwan A, Romsauerova A, Yerworth R, Horesh L, Bayford R, Holder D. Design and calibration of a compact multi-frequency EIT system for acute stroke imaging. Physiol Meas. 2006;27(5):S199.
    Awada K, Jackson DR, Baumann SB, Williams JT, Wilton DR, Fink PW, Prasky BR. Effect of conductivity uncertainties and modeling errors on EEG source localization using a 2-D model. IEEE Trans Biomed Eng. 1998;45(9):1135–1145.
    Pethig R. Dielectric and electrical properties of biological materials. J Biolectr. 1985;4(2):vii–ix.
    Cheney M, Isaacson D, Newell JC. Electrical impedance tomography. SIAM Rev. 1999;41(1):85–101.
    Woo EJ, Lee SY, Mun CW. Impedance tomography using internal current density distribution measured by nuclear magnetic resonance. In: International Symposium on Optics, Imaging, and Instrumentation. Bellingham, WA: International Society for Optics and Photonics; 1994:377-385.
    Oh SH, Lee BI, Woo EJ, Lee SY, Cho MH, Kwon O, et al: Conductivity and current density image reconstruction using harmonic Bz algorithm in magnetic resonance electrical impedance tomography. Phys Med Biol. 2003;48(19):3101.
    Jeong WC, Chauhan M, Sajib SZ, Kim HJ, Serša I, Kwon OI, Woo EJ. Optimization of magnetic flux density measurement using multiple RF receiver coils and multi-echo in MREIT. Phys Med Biol. 2014;59(17):4827.
    Oh TI, Jeong WC, McEwan A, Park HM, Kim HJ, Kwon OI, Woo EJ. Feasibility of magnetic resonance electrical impedance tomography (MREIT) conductivity imaging to evaluate brain abscess lesion: in vivo canine model. J Magn Res Imaging. 2013;38(1):189–197.
    Kim HJ, Kim YT, Minhas AS, Jeong WC, Woo EJ, Seo JK, Kwon OJ. In vivo high-resolution conductivity imaging of the human leg using MREIT: the first human experiment. IEEE Trans Med Imaging. 2009;28(11):1681–1687.
    Katscher U, Voigt T, Findeklee C, Vernickel P, Nehrke K, Dossel O. Determination of electric conductivity and local SAR via B1 mapping. IEEE Trans Med Imaging. 2009;28(9):1365–1374.
    Liu J, Zhang X, Van de Moortele P-F, Schmitter S, He B. Determining electrical properties based on B1 fields measured in an MR scanner using a multi-channel transmit/receive coil: a general approach. Phys Med Biol. 2013;58(13):4395.
    Liu J, Zhang X, Schmitter S, de Moortele V, He B. Gradient-based electrical properties tomography (gEPT): a robust method for mapping electrical properties of biological tissues in vivo using magnetic resonance imaging. Magn Res Med. 2014;74(3):634–646.
    Liu C, Li W, Argyridis I. Imaging electric conductivity and conductivity anisotropy via eddy currents induced by pulsed field gradients. In: Proceedings of the Annual Meeting of the International Society for Magnetic Resonance in Medicine; May 10-16, 2014; Milan, Italy. Abstract 0638.
    van Lier ALHMW, ATC, van den Berg CAT, Katscher U. Measuring electrical conductivity at low frequency using the eddy currents induced by the imaging gradients. In: Proceedings of the Annual Meeting of the International Society for Magnetic Resonance in Medicine; May 5-11, 2012; Melbourne, Australia. Abstract 5955.
    Mandija S, van Lier AL, Katscher U, Petrov PI, Neggers SF, Luijten PR, van den Berg CA. A geometrical shift results in erroneous appearance of low frequency tissue eddy current induced phase maps. Magn Res Med. In press.
    Liu Q, Hughes D, Allen P. Quantitative characterization of the eddy current fields in a 40-cm bore superconducting magnet. Magn Res Med. 1994;31(1):73–76.
    Spencer RL, Ware M. Computational physics 430: partial differential equations; 2012. Accessed August 12, 2015.
    Dib N, Weller T. Two-dimensional finite difference time domain analysis of cylindrical transmission lines. Int J Electronics. 2000;87(9):1065–1081.
    Liu F, Crozier S, Zhao H, Lawrence B. Finite-difference time-domain-based studies of MRI pulsed field gradient-induced eddy currents inside the human body. Concepts Magn Res. 2002;15(1):26–36.
    Schneider JB. Understanding the finite-difference time-domain method; 2010.∼schneidj/ufdtd. Updated December 6, 2013. Accessed August 12, 2015.
    Witwer JG, Trezek GJ, Jewett DL. The effect of media inhomogeneities upon intracranial electrical fields. IEEE Trans Biomed Eng. 1972(5):352–362.
    Su J, Zheng B, Li SFY, Huang SY. Further study of the effects of a time-varying gradient fields on phase maps—theory and experiments. In: Proceedings of the Annual Meeting of the International Society for Magnetic Resonance in Medicine; May 30–June 5, 2015; Toronto, Canada. Abstract 3291.
    Kato H, Kuroda M, Yoshimura K, Yoshida A, Hanamoto K, Kawasaki S, Shibuya K, Kanazawa S. Composition of MRI phantom equivalent to human tissues. Med Phys. 2005;32(10):3199–3208.
    Scott G, Joy M, Armstrong R, Henkelman R. Sensitivity of magnetic-resonance current-density imaging. J Magn Res. 1992;97(2):235–254.
    Schofield MA, Zhu Y. Fast phase unwrapping algorithm for interferometric applications. Opt Lett. 2003;28(14):1194–1196.
    Wu B, Li W, Guidon A, Liu C. Whole brain susceptibility mapping using compressed sensing. Magn Res Med. 2012;67(1):137–147.
    Oran OF, Gurler N, Ider YZ. Feasibility of conductivity imaging based on slice selection and readout gradient induced eddy-currents. In: Proceedings of the Annual Meeting of the International Society for Magnetic Resonance in Medicine; May 30–June 5, 2015; Toronto, Canada. Abstract 0930.
    Griffiths DJ, College R. Introduction to Electrodynamics: Upper Saddle River, NJ: Prentice Hall; 1999.
    Harris C, Haw D, Handler W, Chronik B. Application and experimental validation of an integral method for simulation of gradient-induced eddy currents on conducting surfaces during magnetic resonance imaging. Phys Med Biol. 2013;58(12):4367.


Download the article PDF (831.89 KB)

Download the full issue PDF (12.92 MB)

Mobile-ready Flipbook

View the full issue as a flipbook (Desktop and Mobile-ready)