### Theory

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

**J**_{e} 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

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

∂B∂t 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.

### Simulation

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 *t*_{n+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 *t*_{n}. After updating each electric field component, Faraday's law is then used to update each magnetic field component at *t*_{n+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.4 ps, where *c* is the speed of light and *dx*_{i} is the distance between simulation grid points in the *i*^{th} 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·10^{7} 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 ∂B∂t 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 ∂∂t→1dt+δ. 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∂ρfree∂t=−∇ ⋅ 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 ∂ρfree∂t 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.

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 Gibbs

^{1}, Chunlei Liu^{2}^{1}Brain Imaging and Analysis Center, and^{2}Department of Radiology, Duke University School of Medicine, Durham, NC## Abstract

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.

## Introduction

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

B_{1}field mapping, albeit at the Larmor frequency, which is 128 MHz for 3T scanners (14). TheB_{1}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 fromB_{1}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 (17–19). 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γ ∫ B z d t , where γ is the gyromagnetic ratio and

B_{z}is thez-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.

## Methods

## Theory

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:

## (1)

BandEare 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:J_{e}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## (2)

B. The secondary field affectsB. 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.## Simulation

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 pointt_{n+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 pointt_{n}. After updating each electric field component, Faraday's law is then used to update each magnetic field component att_{n+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,randr′are spatial vectors,Vis the voltage created by ρ_{free}, andE_{ρ}is the electric field created byV. 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 ind t ≈ Σ d x i 2 c = 4.4 ps , where

r, ϕ, andzwas 31 × 28 × 82, with a 1-cm resolution inrandzand π/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 meanscis the speed of light anddx_{i}is the distance between simulation grid points in thei^{th}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·10^{7}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∂ B ∂ t and not by the magnitude of

E_{ϕ}direction until reaching a constant final value. By directly increasingE_{ϕ}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 byB; 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 shortdt, 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∂ ∂ t → 1 d t + δ . 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∂ ρ f r e e ∂ t = − ∇ ⋅ J , where ρJ = σ E , allowing σ and ϵ to be inhomogeneous and assuming linear material, one can expand the conservation of charge to

_{free}is the free-charge density. By applying Ohm's law,## (3)

The charge creates a voltage, which is equal to

## (4)

randr′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 usedntime points later ifrisngrid points from a source atr′. 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∂ ρ f r e e ∂ t 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.

x,y, andz, anddt= 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 ofThe 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.

## Results

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

Each simulation with a dielectric phantom produces a magnetic fieldB E ∝ − ∂ B ND / ∂ t . The magnitude of σ p = 1 , 10, 100, and 1000 S/m, respectively. Note that

B_{D}. The magnetic field from the simulation with no dielectric isB_{ND}. The difference between the 2 represents the magnetic fields generated by eddy currents, denoted asB_{E}. The time course of each simulation at a fixed point in the phantom is shown in Figure 3A, and theB_{E}for each dielectric phantom is shown in Figure 3B. Maxwell's equations dictate thatB_{E}opposes changes in the magnetic field. This opposition creates a delay inB_{D}compared withB_{ND}. This delay increases in magnitude and duration as conductivity increases. For materials with a large conductivity and hence long-lasting eddy currents,B_{E}essentially opposes all changes and is nearly the opposite of the field driving the eddy currents, that is,B_{E}≈ −B_{ND}. At low conductivities, the eddy currents decay too quickly to accumulate over time, andB_{E}depends on geometry as explained in the following section, but the maximum value ofB_{E}at the point shown in Figure 3 was 0.05, 0.5, 3, and 6.5 μT forB_{E}initially almost completely opposes the maximum ofB_{ND}in the 1000-S/m conductivity phantom.## Figure 3.

Time courses of simulated magnetic fields under various conductivity values. Because∂ B ∂ t 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 ofB ′ E = ∂ I ∂ t ∗ A ⋅ e − b ⋅ t is fit to τ = ln ( 2 ) b and τ = 2, 11, and 119 ns for σϕ = γ ⋅ ∫ B z , E d t , where

B_{E}, it is useful to think ofB_{E}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 fieldB_{E}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_{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 byB_{z,E}is thez-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 monoexponentialA ⋅ e − b ⋅ t , which represents current decay in a single loop, is the impulse response function. After optimizing τ = ln ( 2 ) b = 2 ns . The fit is less acceptable as conductivity rises and is not reasonable for σ

Aandb, the predicted output closely follows the time course forB_{E}when σ_{p}= 1 S/m. The half-life of the mono exponential is_{p}= 1000 S/m.## Simulation of the Spatial Characteristics of the Eddy Current-Induced Magnetic Field

The largest determinant of

B_{E}is conductivity, but geometry plays a role as well. The magnitude ofB_{E}at a given point depends on the value ofB_{ND}at that point, butB_{E}/B_{ND}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 showsB_{D}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 ofB_{D}rises to values detectable by MRI, the boundaries between conductive regions are smoothed, and the percentage of variation drops. WhenB_{D}reaches 1 μT,B_{D}varies by only 1.5 nT in regions whose conductivity varies by 3.5 S/m. The ability to resolve variations inB_{D}for a given magnitudeB_{D}is independent of the underlying conductivity, althoughB_{D}rises more slowly in highly conductive material, increasing the duration of variations inB_{D}.## 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 ϕ atz= 4.5 cm andr= 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 andt= 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 andt= 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 ofB_{E}, the maximal current density was about 3.76 A/m^{2}, with an average magnitude of 0.87 A/m^{2}. 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∂ I ∂ t 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.

## Discussion

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 wavese i ( k · r − ω t ) , where

kis 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 vectorkis 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)## (5)

B_{E}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 byx 1 / 2 = ln ( 2 ) κ . Values of

x_{1/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,x_{1/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.## Conclusions

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.

## Notes

[1] Abbreviations:

TE

Echo time

EC

eddy current

EIT

electrical impedance tomography

FFT

fast Fourier transform

FOV

field of view

FDTD

finite-difference time-domain

MRI

magnetic resonance imaging

MREIT

magnetic resonance EIT

RF

radiofrequency

TR

repetition time

VSHARP

variable kernel sophisticated harmonic artifact reduction for phase data

## Acknowledgments

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.

## References

## Journal Information

Journal ID (nlm-ta): tom

Journal ID (publisher-id): TOMOG

Title: Tomography

Subtitle: A Journal for Imaging Research

Abbreviated Title: Tomography

ISSN (print): 2379-1381

ISSN (electronic): 2379-139X

Publisher: Grapho Publications, LLC (Ann Abor, Michigan)

## Article Information

Copyright statement: © 2015 The Authors. Published by Grapho Publications, LLC

Copyright: 2015, Grapho Publications, LLC

License (open-access, http://creativecommons.org/licenses/by-nc-nd/4.0/):

This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Publication date (print): December 2015

Volume: 1

Issue: 2

Pages: 125-135

Publisher ID: TOM-00142-15

DOI: 10.18383/j.tom.2015.00142

## PDF

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)