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 (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 , 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.
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:
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).
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 , 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 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 . 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, where ρfree is the free-charge density. By applying Ohm's law, , allowing σ and ϵ to be inhomogeneous and assuming linear material, one can expand the conservation of charge to
The charge creates a voltage, which is equal to
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 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, and 0.24 S/m, with and 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 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 . 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 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.
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 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 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 , 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.
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.
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 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 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.
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.
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 , 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)
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 . 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.
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.