1.

Introduction

Electrical properties within the human body, primarily permittivity and conductivity, vary spatially due to natural contrasts created by fat, bone, muscle and various organs.¹ In addition, various disease states are known to cause changes in electrical properties. For example, Haemmerich et al.² showed that hepatic tumors exhibit elevated in vivo conductivity compared to normal tissue. Using skin electrodes to make impedance measurements, Songer³ demonstrated that muscle tissue impedance dropped at the linear rate of $0.148\mathrm{ohm}/\min$ during tourniquet “surgery,” intended to induce ischemia (see Ref. 3, Table 8.1, $0.8<R^2<0.99$). For these reasons, considerable effort has been made to devise strategies to image electrical property contrasts, especially through means not involving physical contact. Though magnetic induction tomography (MIT) is perhaps the most common no-contact approach currently under development,⁴ other no-contact methods have also attracted interest, such as noise tomography which attempts to image electrical conductivity by using MRI coils to detect and process the spatial variation of thermally induced electronic noise developed within conductive tissues.⁵

To help further motivate the need for MIT, we note that our interest here at Kimberly-Clark lies in a desire to have a low-cost, portable approach for imaging tissues in the immediate vicinity of an implanted medical device, such as an endotracheal tube or gastrostomy tube. An inexpensive imaging method with even moderate resolution could quickly ascertain proper positioning of the device, as well as the onset of possibly adverse effects of the device on tissues, such as the development of ischemic conditions due to excessive tissue pressure. A second application of MIT, perhaps more likely to see wide use, involves the low-resolution imaging of pressure-induced ulcers or bed sores which have been shown to result from pressure-induced ischemia in tissues.⁶ Assuming that such ischemic conditions cause an elevation of electrical conductivity, as suggested by the low-frequency results of Songer³ ($<1\mathrm{MHz}$), then MIT could offer a way to evaluate at-risk patients for early onset of pressure ulcers. In an effort to establish a connection between ischemia and elevated conductivity at frequencies exceeding 10 MHz, Feldkamp and Heller⁷ used a simpler version of the device described herein to show increased conductivity in human extremities when elevated for 60 s—meant to produce a mild level of ischemia. In particular, their results showed increased conductivity change in response to elevation change at the right mid-volar forearm that became more pronounced in subjects exhibiting low blood pressure (Figs. 3 and 4).

Most commonly, MIT imaging techniques involve placement of a large number of coils near the sample and an attempt to build an image based upon the measured mutual inductance of coil pairs within an array placed near the object.^8,9 Some authors have indicated a desire to reduce the number of coils, with one approach calling for the repositioning of either the coils in the array or the object itself.¹⁰ Whether MIT uses many or few coils, improved imaging calls for as many measurements as possible. With that as an incentive, research reported here aims to accomplish the task of MIT imaging with just one coil, functioning as both source and detector. A single coil is given RF excitation at a fixed frequency while the coil’s ohmic loss due to inductive coupling with the sample is measured. The intent is to allow the single coil to be sufficiently mobile that it can be rapidly placed into a very large number of locations and orientations, all while data collection proceeds. Inversion of data into a conductivity image would follow the completion of a physical “scan” which further helps to localize the target for meshing purposes.

In order to accomplish such a task, it is helpful to have a quantitative analytical model that establishes a clear relationship between an electrical property and a measured coil property. The only coil property we strive to compute from theory is ohmic loss resulting from inductive coupling between coil and object. To reduce the level of difficulty, we develop a solution of Maxwell’s equations for an arbitrary conductivity distribution, but we consider permittivity and magnetic permeability as spatially uniform. As shown by Harpen,¹¹ who considered the interaction of a conductive sphere with an EM field produced by a Helmholtz coil, permittivity will only have an effect on the imaginary part of coil impedance change, which is discarded here. Thus, even though our samples are expected to exhibit nonuniform permittivity, this is not expected to hamper efforts to measure ohmic loss. For example, Harpen^11,12 showed only small effects on coil self-resonant frequency as long as coil-sample contact is avoided. “Detuning,” resulting from capacitive coupling, was no more than 0.1 MHz at a 60 MHz excitation frequency.

Our coil geometry is kept as simple as possible—a collection of concentric circular loops, all lying within a common plane and connected in series, with the transient current considered to have the same value at all points along the loops. With a 4-cm diameter, the typical loop size is quite small compared with the wavelength (24 m) at 12.5 MHz, the intended operating frequency. A single current loop consisting of infinitely thin wire, placed parallel to a uniformly conductive half-space, has been treated analytically using a separation of variables approach by Zaman et al.¹³ Their solution was developed primarily for calculations involving metals and so was specialized further for the case of high conductivity. Here, the conductivity distribution is permitted to vary arbitrarily in space while a solution for the electric field is pursued in the limit of small conductivity ($<10\mathrm{S}/\mathrm{m}$), though higher conductivity is acceptable if the frequency is lowered. This leaves us with the task of solving a partial differential equation (PDE) for the electric field and then computing ohmic loss. With the formula in hand, comparison with experimental results is shown for the case of a coil consisting of five circular loops placed into each of two very closely spaced planes (0.5 mm plane separation), for a total of 10 loops placed on a two-layer printed circuit board (PCB) and connected in series. Model assumptions are realistic provided loops of the coil do not come into intimate contact with the dielectric, enabling capacitive coupling. Close contact to a dielectric and close spacing of loops within a plane are expected to lead to capacitive as well as inductive losses—only the latter is desired.

With the goal of demonstrating the feasibility of imaging using a single induction coil in a scanning mode, the agenda of this paper is as follows: a brief overview is provided of the analytic formula which facilitates mapping between electrical conductivity distribution and coil loss measurements—the full details are given in Appendix E; a description of coil loss measurements for the proposed coil geometry is given; a comparison is made of loss measurements with predictions of loss from the analytic formula while the coil is placed at known distances from a tank of aqueous potassium chloride; next, the procedure for solving the inverse problem, starting with a set of loss measurements and culminating in an electrical conductivity image, is presented; the inversion algorithm is then tested using a prescribed conductivity distribution that mimics size features expected for the thoracic spinal column; finally, clinical results are given to demonstrate preliminary “proof of concept” for the single-coil, scanning MIT imaging strategy by imaging a portion of the thoracic spinal column where natural electrical property contrasts exist between muscle and bone tissues.

2.

Arbitrary Conductivity Analytical Solution for Inductive Loss in a Multiloop Coil

For an RF-driven coil consisting of concentric, circular wire loops, all lying in the same plane, Maxwell’s equations may be solved in the quasistatic limit, leading to a solution that is valid for an arbitrary electrical conductivity distribution. Loops are connected in series, each carrying the same sinusoidal current. Assumptions required to build a solution for the electric field are that permittivity and permeability are uniform, electric charge is uniform, conductivity is isotropic and small ($\le 10\mathrm{S}/\mathrm{m}$), and loop wires have infinitesimal radii. The electric field can then be used to compute the EMF developed in the coil and associated impedance change due to coil-object inductive coupling. As developed in greater detail in Appendix E, coil loss (or the real component of impedance change) is given as a sum of integrals:

Arguments for the toroid (or ring) function $Q_{1/2}$ lie in the interval $1<\eta <\infty$, and are related to field position by

With the coordinate system origin located at the center of the set of concentric loops and $Z$-axis normal to the plane of the loops, other symbols in Eq. (1) are listed here

The ring function is readily evaluated using the hypergeometric series form as given in Appendix C, so that Eq. (1) provides a mapping between coil loss measurements and the conductivity distribution. As shown in Secs. 3 and 4, which focus on coil loss measurement, Eq. (1) gives a quantitative prediction of loss for a known conductivity distribution. As a result, Eq. (1) can be used to simulate measurements of coil loss for any coil position or orientation, providing a measurement set which can then serve as a starting point for the evaluation of inversion algorithms based upon Eq. (1) —an example is given in Secs. 5 and 6. The remaining sections then explore the feasibility of single-coil, scanning MIT imaging with a particular choice of inversion algorithm based upon Tikhonov regularization. The small clinical discussed in Sec. 7 shows preliminary feasibility of the approach by manually scanning the interscapulae region for six subjects, followed by inversion of data to produce images of the thoracic spinal region.

3.

Coil Geometry and Coil Loss Measurement

Coils built to test the theory consist of a total of 10 concentric, circular loops. Five loops are placed in an upper plane, stacked 0.5 mm above a lower plane also consisting of five identical loops. Circular loop radii are set at 4, 8, 12, 16, and 20 mm in either plane. Figure 1 shows the essential geometry of the stacked sets of loops, created with the help of the CAD module in FEKO (EMSS, Inc.) using “equivalent” wire segments (equivalent in the sense that segment diameters are assigned so that segment perimeter matches the actual PCB trace perimeter). The loops consist of 1 or 0.5 mm widths, 2 oz. copper traces on a multilayer PCB (1.5 mm thick), wired in series as shown in the figure. The smaller radial traces used to connect the loops in series are arranged so that current flow in the upper set of connectors is opposite to the current in the lower set of connectors, with the intent that radial connectors contribute nothing to the field. Ideally, the field produced by the stacked loops is only that due to the circular traces. Perfect circular loops were preferred over spirals since spirals do not match theory nearly as well as circular loops, especially as spacing grows.

Fig. 1

Geometry of two-layer, multiloop coil used for MIT imaging; loop radii—4, 8, 12, 16, and 20 mm; segment charges are shown for case with 1 V excitation at 12.5 MHz; 84 pF of capacitance and $20\mathrm{k}\mathrm{\Omega }$ of resistance were placed in parallel with the coil to resonate at 12.26 MHz; segment diameters were set at 0.6 mm to give a perimeter equivalent to that of a 1 mm, 2 oz. Cu PCB trace.

Furthermore, the inner loop of the lower layer (facing the specimen) is connected to the instrument common, whereas the inner loop of the upper layer is connected to the excitation. Coloration on loops (see Fig. 1) represents the amount of charge developed on coil wires, as computed from a simulation using the method of moments’ code known as FEKO (EMSS, Inc.). Inner loop charging is the highest, since the potential difference between the innermost loops of the two layers is greatest. Using a network analyzer, coil capacitance was measured to be about 14 pF (0.5 mm trace width) or 20 pF (1.0 mm trace width). Intrinsic coil capacitance causes no particular issues inasmuch as additional capacitance is added anyway in the form of a varactor to provide a means to tune the tank circuit. Capacitive coupling with the target, however, is another issue. Though the development of additional capacitance across the coil due to specimen proximity is acceptable, the appearance of losses associated with the additional capacitance is not.

An equivalent circuit can be used to represent the coil in order to facilitate measurement of ohmic loss in the coil. Here, we represent the coil as a series combination of ideal inductance $L$ and ideal resistance $R$, which is then placed in parallel with ideal capacitance $C$—Appendix A gives a schematic of the equivalent circuit. Capacitance is mostly due to the interaction of loop segments with each other. However, capacitive coupling with the material can happen but would only constitute a problem if there were associated losses. An expression for circuit admittance $Y$ (reciprocal of impedance) can be easily written down for the proposed equivalent circuit. Since ohmic loss is expected to be very small, we specialize to the case that $R\ll \omega L$, where $\omega$ is angular frequency:

Only the real part of admittance is of interest, so that the loss $R$ can be measured as

Note that $R$ is the quantity that varies in response to eddy current generation within the space occupied by a conductive medium. A variety of ways exist to measure admittance. We use an AD8302 phase and gain detector from Analog Devices to simultaneously measure phase angle (between voltage and current) and voltage gain—the ratio of voltage across the coil to voltage across a precision sense resistor placed in series with the coil. The AD8302 measures gain by the use of a matched pair of demodulating log amplifiers having a 60 dB range. Appendix D shows a portion of the analogy circuitry needed to provide excitation to the AD8302. From gain and phase, computing the real part of admittance either at or off resonance is straightforward. Note that the real part of admittance is independent of capacitance, which is fortunate since some amount of capacitive coupling can be expected when the coil-specimen distance is small. Moreover, by placing a tunable varactor in parallel with the coil, the circuit can be tuned to any convenient point for measurement without affecting $Y_{\mathrm{re}}$.

We have measured $Y_{\mathrm{re}}$ both at resonance (minimum admittance) and at arbitrary points off resonance, obtaining the same result either way. Most commonly, off-resonance measurement is used here, at a phase angle of about $-60\mathrm{deg}$. Inasmuch as loss $R$ contains contributions intrinsic to the coil as well as from inductive coupling with a nearby conductive specimen, a free-space value for $R$ must first be determined and subtracted from all subsequent measurements. All experiments were performed using a fixed frequency, crystal controlled excitation at 12.5 MHz. Typically, admittance measurement precision is about $0.02\mu S$ while the span of admittance, from free space to within about 3 mm of human tissues, is about $12\mu S$—free space admittance is $112.55\mu S$. Though precision is very good, drift during 10 min of instrument warm-up was a nominal $0.05\mu S/\min$. To manage drift, appropriate time was allowed for instrument warm-up. Furthermore, free-space admittance values were obtained before and after a set of measurements and then used to supply a background correction based upon some estimated amount of time per measurement. A photo of electronics, coil, and enclosure is shown in Appendix A.

To implement Eq. (4), coil inductance is required. For loops of different radii, the mutual inductance contribution to overall coil inductance was computed from an equation given by Conway,¹⁴ rewritten here in a form that uses the zero order, half integer degree toroid function:

Loop radii are given by ${\rho }_j$ and ${\rho }_k$, while magnetic permeability is given by $1.2566\mu \mathrm{H}/\mathrm{m}$. The toroid function is available from tables or easily computed from a hypergeometric series. If the self-inductance of any individual loop $j$ is given by $L_{sj}$, then the mutual inductance between pairs of loops with the same radius, but in different layers, was just taken to be $L_{sj}$—in other words, the coupling constant is taken as unity. Self-inductance of individual loops was computed from simplified equations given by Terman,¹⁵ using a “wire” diameter having a circular perimeter equal to the trace perimeter. PCB traces are built at a width of either 1.0 or 0.5 mm using 2 oz. copper, equivalent to about 0.0694 mm thickness. Based upon equivalent perimeter, a 1.0 mm trace has an equivalent circular wire diameter equal to 0.68 mm and a 0.5-mm trace width has an equivalent diameter of 0.36 mm. The average of these two is close to 0.5 mm, which is the wire diameter used for the inductance calculation of either coil. According to the circular loop model, smaller diameter conductors will give a somewhat larger inductance, which agrees with network analyzer measurements. Overall, inductance of the 10-loop coil consists of

Inductance computed in this way was found to be $2.155\mu \mathrm{H}$. Inductance was also measured on an Agilent network analyzer by first obtaining the self-resonant frequency of the coil, and then a second time when placed in parallel with a precision capacitor ($\pm 5\%$). From the two resonant frequencies, $L$ was computed to be $2.132\mu \mathrm{H}$. The uncertainty in the ceramic capacitor, combined with the effects of added solder and capacitor leads, on both inductance and capacitance, makes the network analyzer result less reliable than the computed inductance—so the latter is used for both coils. However, the difference is only 1% which validates the inductance computation quite well.

4.

Comparison of Theory and Experiment for 10-Loop Planar Coil

Equation (1) may be used to compute coil loss when the planar coil is placed at various distances away from a plastic tank filled with aqueous potassium chloride. The plane of the coil was kept parallel to the liquid surface (within $\pm 0.5\mathrm{deg}$) and placed at various distances from the liquid surface using a precision $XYZ$ translator (Newport Universal Motion Controller/Driver Model ESP7000; equipped with M-MTM 150 CC1 motorized stage, capable of 150 mm of travel and $1.0\mu \mathrm{m}$ precision). In spite of superb translator precision, the measurement accuracy of coil-liquid-surface distance was about $\pm 0.5\mathrm{mm}$. Critical to imaging applications is the requirement that theory and experiment match nearly exactly or at least to within a known scale factor. This was tested for two variants of coil geometry; the only difference between these two is the trace width, 1 or 0.5 mm.

Using an $XYZ$ translator to change the elevation of a mounted coil placed above a $30\times 30\times 13{\mathrm{cm}}^3$ deep tank of aqueous KCl, admittance change relative to free space was measured and then used to compute inductive loss according to Eq. (4). This was then compared to the same quantity computed from Eq. (1), which is equivalent to $R$ in Eq. (4). Results for both coil types are shown in Figs. 2 and 3. Loss from Eq. (1) was evaluated by using a finite element partition of the tank using pentahedral elements and nine-point integration. The tank was only approximately rectangular but was treated as such in building a finite element mesh.

Fig. 2

Comparing theoretical predictions from Eq. (1) and (infinite slab) with measurement on coil with 1 mm trace. No data fitting or scaling is used here; rather, the comparison is “head-to-head.”

Fig. 3

Comparing a theoretical prediction from Eq. (1) with measurement on coil with 0.5 mm trace. No data fitting or scaling is used here; rather, the comparison is “head-to-head.” Also included is the infinite slab result [derivation not included, but similar to Eq. (1)].

A conductance probe (sympHony™—VWR) was used to measure the electrical conductance required by Eq. (1) and the measured conductance was $2.86\mathrm{S}/\mathrm{m}$ for the 1.8% KCl solution (by weight). No data fitting or scaling is used to help improve fit in Figs. 2 and 3; rather, experimental and theoretical data are compared “as is.” Equation (1) calculations are only limited by the assumptions listed in Sec. 1 and numerical error. Assumptions behind Eq. (1) include small conductivity and uniform permittivity throughout all space. The latter assumption is not at all met by our system since permittivity abruptly changes at the air-solution interface. But as literature shows,^11,12 lowest order terms developed for the real component of impedance under nonuniform permittivity conditions [see his Eq. (21)] do not contain permittivity (or, equivalently, speed of light). Thus, a uniform permittivity approximation is not an issue if only conductivity is of interest. The situation is different when magnetic permeability varies in space, for then the theory behind Eq. (1) would not be useful for the accurate prediction of coil impedance. However, magnetic permeability of human body tissues and salt solutions is not significantly different from vacuum. Figures 2 and 3 also contain infinite slab results for a slab having a thickness equal to the finite slab—the infinite slab equation was derived in a manner similar to Eq. (1) though no mesh was needed for calculations.

Results for the narrow trace coil show somewhat better agreement with theory, which is not surprising since the theory was developed for infinitely thin conductors. The remaining disagreement between theory and experiment arises from several possible contributions: (1) coil not perfectly parallel with aqueous solution interface and vertical positioning accuracy; (2) computation/measurement of coil inductance; (3) admittance accuracy was judged to be about $\pm 0.2\mu S$, corresponding to a loss error of $\pm 0.006\mathrm{ohm}$, as computed from Eq. (4); (4) shape of container holding the KCl solution was not perfectly rectangular, and no attempt was made to work with the true container geometry, which was conceivable but time consuming; (5) numerical errors in computation of loss from Eq. (1)—accuracy could be improved by refining the mesh used to discretize the tank; and (6) short range losses due to capacitive coupling with target—no estimate was made of this. In either Fig. 2 or 3, coil loss is shown to decline substantially when coil-target separation reaches about one coil diameter. However, as Sec. 6 shows, the effective range of coil-specimen interaction can be up to five coil diameters when a specimen becomes very large.

5.

Inversion Algorithm for Single Coil, Scanning MIT

The intended application of Eq. (1) is the inversion of experimentally obtained coil loss data to form an image of the 3-D distribution of electrical conductivity for an arbitrary object—usually one of biological interest, such as human tissues. An inversion routine is developed using Eq. (1), first by discretizing some portion, if not all, of the target object on a finite element mesh. In some cases, discretization of the entire target may be impractical. Using six-node linear pentahedral elements, electrical conductivity is expanded in a linear basis

Substitution of Eq. (7) into Eq. (1), and writing the integral as a summation over all elements in the target geometry, yields for the $n$’th coil loss measurement:

Summation in Eq. (8) is over $M$ nodes and $E$ elements. Each of the measurements, $n$, involves a different position or orientation of the coil and thus affects how $\rho$, ${\eta }_j$, and ${\eta }_k$ are evaluated. Elemental integrals, after multiplication outside by the square root factor $\sqrt{{\rho }_j{\rho }_k}$, have dimensions of volume. Element integrations here are six-point, eight-point or nine-point, with only minor differences observed between those choices in most cases. Integrals are transformed to the local element coordinate system, i.e., in terms of area coordinates and a third coordinate spanning element height, $-1\le \xi \le 1$. In order to evaluate the functions under the integral, consisting of $\rho (\overrightarrow{r})$, $Q_{1/2}[{\eta }_j(\overrightarrow{r})]$ and the basis function ${\mathrm{\Phi }}_m$, it is necessary to determine the global coordinate system vector $\overrightarrow{r}$ that locates each integration point, and then determine from $\overrightarrow{r}$ the vector locating the same integration point from the point of view of the coil coordinate system. Axes in the coil coordinate system need not be aligned with the others.

Equation (8) can be written as a matrix equation—$N$ equations (or loss measurements) and $M$ unknowns (nodal electrical conductivities):

The matrix ${\tilde{T}}_{nm}$ has elements:

Note that the index “$n$” is attached to toroid function arguments, and ${\rho }_n$, to indicate that determination of these quantities is dependent on coil position and orientation. Though not ideal, the number of coil position-orientations selected for making coil-loss measurements will likely be less than the number of mesh nodes so that system (10) will be underdetermined. Regardless, the equations are solved in a least-squares sense while conductivities are forced to be non-negative. Furthermore, it is highly desired to penalize those solutions that involve very rapid changes in electrical conductivity since those solutions are not as likely. As a result, resolution will be limited by smoothing and element size. To penalize steep gradient solutions, a gradient matrix is defined from interpolation functions as

From the gradient matrix, conductivity gradient components within elements are readily computed:

To find a set of nodal electrical conductivities that best explain the observed coil-loss data, subjected to smoothing, we minimize the two-norm:

As the smoothing parameter, $\alpha$, is decreased, the first term becomes more dominant, thus “sharpening” our focus on any structure within the conductivity distribution. The approach is well-known and commonly referred to as Tikhonov regularization.¹⁶ Assignment of $\alpha$ was discussed by Rudnicki and Krawczyk-Stańdo,¹⁷ who surveyed a variety of regularization strategies, in particular the L-curve method. This approach constructs a log-log plot of the residual norm against the gradient norm,¹⁶ typically forming an L-shaped curve. Choosing $\alpha$ where a knee in the L-curve forms has been shown to provide a satisfactory compromise between achieving a small residual and adequate smoothing. The method we use here is a close variant of the L-curve approach and is described further in Sec. 7. Minimizing the 2-norm, or Euclidean norm, is equivalent to solving the following system:

Equation (14) is solved with the further constraint that solution components are non-negative. Ideally, Eq. (14) is sufficient for our needs. However, there are cases where instrumental offset related to drift and aging should be accounted for. This can be built into Eq. (14) by adding an additional unknown—the instrument offset, into row one of column vector ${\breve{\sigma }}_m$. This requires that a new first column of 1’s be added to the transfer, or translation, matrix $\hspace{0.17em}\tilde{T}$. Such a modification will then also require modification of $\hspace{0.17em}\tilde{G}$ to have a new first row and first column of just zeroes. System (14) is solved with the routine WNNLS (weighted, non-negative least squares), useful only for small to moderate-sized problems (Naval Surface Warfare Math Library). WNNLS is more general than is needed here since equality constraints could be included. Since equality constraints are not used, WNNLS is just an NNLS type problem. The strategy for enforcing non-negativity is the so-called active set method as discussed in Sec. 3.2 of Haskell and Hanson.¹⁸ Though there are variants to this method, active and passive sets of unknowns are defined—the former set consisting of mesh nodes which would violate non-negativity and the latter consisting of nodes with positive conductivity. Through an iterative process described in Sec. 3.2 of their paper, which they state will always converge, unknowns within the active set transfer to the passive set until a solution is eventually found with all nodes satisfying non-negativity.

Because few actual measurements were available per subject in the study described in Sec. 7, preliminary testing with “virtual measurements” was carried out to assess the feasibility of single-coil MIT imaging with fewer measurements than desired. The next section shows that even with 61 “simulated measurements” (132 per subject used in clinical), structures comparable in size with the spinal column and canal can be visualized.

6.

Imaging a Prescribed Conductivity Distribution from Simulated Measurements

Because the number of actual measurements per subject in the clinical study was limited, preliminary testing with simulated measurements was carried out to assess the feasibility of single-coil MIT imaging with fewer measurements than desired. Here, testing with 61 “simulated measurements” is done over a 2-cm thick, $9\times 9{\mathrm{cm}}^2$ slab with a conductivity distribution prescribed by

A simulated measurement consists of using Eq. (1) to compute coil loss for a particular coil position or orientation over the specimen. These simulated measurements are then used in the imaging algorithm as though they were actual measurements.

Note that the last term in Eq. (15) is centered on the slab, but the sine terms are not, which contributes to the asymmetrical features shown in the central region of Fig. 4. The prescribed conductivity distribution is meant to contain some asymmetry as well as structures comparable in size with the spinal column and canal in order to test inversion algorithm capability. Simulated measurements were made in just two planes above the 2-cm thick slab. At 2.0 mm above, 36 measurements: $x$, $y=2$, 3, 4, 5, 6, 7; at 6 mm above, 25 measurements: $x$, $y=2.5$, 3.5, 4.5, 5.5, 6.5. Staggering of the coil position from layer to layer is a good sampling practice. Clearly, the image inversion result shown in Fig. 5 indicates that such structures, especially the central “bump” and asymmetry, can be visualized even with fewer measurements than may be desired, and much less than the number used in the clinical study. Distortion is greatest toward the corners and edges of the image, where sampling frequency is less. Contour plots in Figs. 4 and 5 are obtained on a slice at a depth of 0.5 cm beneath the surface.

Fig. 4

Prescribed conductivity distribution used to create 61 virtual coil loss measurements using 4 cm diameter, 5-loop coil identical with that used in the study, over a 2-cm thick slab. Coil positions for sampling are noted in the text. Contours are shown at intervals of $0.2\mathrm{S}/\mathrm{m}$, with the contour at 0.4 given as reference.

Fig. 5

MIT image created on the basis of the 61 virtual measurements—compared with original in Fig. 4. Contours are shown at intervals of $0.2\mathrm{S}/\mathrm{m}$, with the contour at 0.4 given as reference.

Of interest is the range of interaction expected between a coil and nearby conductive sample. Coil-loss results are shown in Fig. 6 for a square slab, 12-cm thick, and spanning a range of side dimensions up to 50 cm. Conductivity is set at a nominal $1.0\mathrm{S}/\mathrm{m}$ within the slabs. The coil consists of five circular, concentric loops having radii 0.5, 1.0, 1.5, 2.0 and 2.5 cm and is placed parallel to the specimen at 1.0 mm above the slab center. Also shown is the result from a model of a 12-cm thick slab having an infinite extent along $X$ and $Y$ axes. It is noted here that the computed loss only becomes asymptotic to the infinite slab result when the slab dimension reaches about five coil diameters. Thus, a coil has significant “reach.”

Fig. 6

Impedance change in 5-loop coil (5 cm diameter) when positioned 1.0 mm above a 12-cm thick square slab of variable side dimension. The results are compared with the infinite slab value computed from an analytical formula. The finite element mesh used for computing impedance change from Eq. (1) has resolution of about 0.5 cm for slabs up to 30 cm along a side, thereafter 1.0 cm. Frequency is set at 12.5 MHz.

7.

Imaging the Thoracic Spine—Clinical Testing

A number of interrelated reasons conspire to make actual imaging much more difficult. For one, we now need to contend with measurement noise and inaccuracies—which includes coil position and orientation, not just coil ohmic loss. The other is that the target is generally much larger than what we may wish to mesh—those nodes that are remote are expected to contribute far less than those that are near, and the task of accurately apportioning those contributions to measured loss can be daunting. However, a large number of distant or remote nodes might be expected to produce a collective response in the coil that is non-negligible but still very difficult to correctly apportion to those nodes. Thus, in view of the simplicity of the proposed inversion algorithm, the extent of target meshing will be limited herein. Because of the natural variability of electrical conductivity in the vicinity of the thoracic spine, due to interspersed bone, muscle, fat and connective tissue, our imaging effort centered on this location. Another reason is the relative flatness that exists in the interscapulae region.

Six male subjects, all in good health, were recruited for single-coil MIT testing under study number 13–356 as approved by the New England Institutional Review Board on October 8, 2013. Each subject was required to sign an informed consent document prior to participation. In order to acquire coil-loss measurements at a reasonably large number of positions, a very simple locator-template was constructed to control positioning, as shown in Fig. 7. The template was sketched onto two low density polyethylene (LDPE) sheets, each having about 1.60 mm thickness, one serving as the bottom sheet, contacting the back, the other as a top sheet. To build distance from the spine, ethylene vinyl acetate (EVA) foam sheets were progressively added in between the LDPE sheets, each sheet is about 2.0 mm. A total of nine EVA spacer sheets, together with the two LDPE sheets, allowed us to span from 1.6 mm up to about 21 mm, excluding the 1 mm buildup on top of the outermost coil trace. Template dimensions are $10.2\times 12.2{\mathrm{cm}}^2$. The template consists of 12 locations where the coil is centered for coil-loss measurement. Together with a total of 11 sheets, a total of 132 coil loss measurements were made per subject. While making the measurements, maintaining coil axis orientation consistent from one location to the next was impossible. As a result, distortion within recovered images was expected. By visual inspection, it is easy to see that the centers of blocks 11 and 12 are the most heavily sampled—these locations are just above and below the template center. Thus, the recovered images beneath those regions are expected to be the most reliable. Spatial sampling frequency becomes less nearer the outer boundaries of the template so that images are expected to be less reliable there (also noted in virtual tests described in Sec. 6).

Fig. 7

$102\times 122{\mathrm{mm}}^2$ template used for guiding placement of the MIT single coil—the PCB supporting the coil has square dimensions that allow it to just fit within any one of the squares.

The mesh used to discretize the space beneath the template (in the body) is shown in Fig. 8—note that it extends somewhat beyond the sampling boundary of the template, given by the middle solid border. The mesh of Fig. 8 was extruded into either 10 (mesh 1) or 9 (mesh 2) layers, with mesh 2 better resolving the region nearer the template and mesh 1 emphasizing greater depths in an effort to facilitate the resolution of the spinal canal. The $Z$-locations of mesh horizons for the two mesh choices are also shown in Fig. 8.

Fig. 8

Details of finite element meshes used for inversion of MIT clinical data—left hand graph shows $Z$-elevations of layer horizons, with 0 cm at greatest body depth and 7 cm at lower boundary of template.

In order to retrieve the greatest amount of information available from the 132 coil loss values obtained from subjects, the smoothing parameter $\alpha$ must be adjusted in a consistent and meaningful way. The so-called L-curve method¹⁶ is one approach for finding an optimal value for $\alpha$, which involves finding the knee in the log-log plot of residual norm versus gradient norm. Here, a similar, but simpler approach involves plotting the residual norm against the log of $\alpha$. Starting from 0.01, $\alpha$ was decreased until the root mean square error [left term of Eq. (13)] showed little further decrease, i.e., where the knee is formed in the proposed plot. For five of the subjects, the results of this exercise are shown in Fig. 9, indicating that $\alpha \approx 0.002$ should be the least value one could use and still represent an appropriate tradeoff between smoothing and a small residual on mesh #1, while 0.003 was judged best for mesh #2. For all images shown here from mesh #2, 0.003 is used, while 0.002 is applied for all images based upon mesh #1. Figure 10 shows results for mesh #1, which are somewhat different than how $\alpha$ affects mesh #2. In either case, the “knee” of the curve is viewed as the best location for choosing $\alpha$, as proposed in the L-curve approach. Also, results appeared to be best with subject #2, so images are shown for subject #2 only. Typically, inversion required about 1 h of processing on a Xeon processor and about 1.3 GB of physical memory.

Fig. 9

Root mean square error [first term, Eq. (13)] as a function of smoothing parameter $\alpha$, indicating the best choice should be about 0.003; similar results were found for mesh #1.

Fig. 10

Root mean square error [first term, Eq. (13)] as a function of smoothing parameter $\alpha$, indicating the best choice should be about 0.002; similar results were found for mesh #2.

Inversion of data from any subject, on either mesh, shared the common and very important feature that nearly all conductivity values across FE mesh nodes fall within the range, from 0 to $2\mathrm{S}/\mathrm{m}$. This result indicates successful inversion, at least in an overall sense, since values outside this range would be unexpected. An exception occurs only along the bottom surface of the rectangular mesh, at points having deepest locations within the body, or most remote from the coil. This is most likely a result of not extending the mesh over the entire tissue volume, so that the collective effect of remote points on coil loss is to cause some amount of overshoot there, essentially allowing those points to “over contribute.” This is illustrated in Fig. 11, which shows a sagittal slice located 2.5 cm offset from and parallel to the spinal column center line. To capture the distribution of conductivity, 0 to $1\mathrm{S}/\mathrm{m}$ transitions are shown from black to orange, 1 to $2\mathrm{S}/\mathrm{m}$ transitions from orange to yellow, 2 to $3\mathrm{S}/\mathrm{m}$ transitions from yellow to white, and finally, greater than $3\mathrm{S}/\mathrm{m}$ is all white—to isolate regions of anomalously high conductivity. First, we note that white only appears at the deepest locations ($Z=0$), suggesting that there is indeed a tendency for very remote points in the body to collectively contribute excessively high values to the deepest mesh nodes.

Fig. 11

Parasagittal slice (mesh #2) with subject in prone position, parallel to spinal column, about 2.5 cm to the right of column center. $Z=0$ is farthest from the coil and lies at the deepest point in the body for either mesh.

Articulations with rib bones nominally begin about 1.75 cm to either side of the spinal column center line (based upon a 3-D model of the human skeleton). The sagittal slice of Fig. 11 is 2.5 cm right of the spine center, so that either transverse spinal processes or ribs should be distinguishable, as indeed they are. From the slice, transverse processes repeat about every 3.75 cm which is the correct spacing as taken from anatomical models using a caliper. Furthermore, transverse processes should lie at a depth of about 4 cm beneath the template, which is the case from the figure ($Z=3$ on axis). Another sagittal slice taken at the same distance to the left of spinal column center is very similar to that shown in Fig. 11. In order to facilitate comparison of Fig. 11 structures to standard reference images, Fig. 12 shows a similar parasagittal slice, offset by about 2.5 cm from the spine center, taken from full body CT scans of a human male cadaver (Visible Human, Research Systems, Inc., version 2.0) placed in a supine position, though shown in the figure as standing. The dimensions of the superimposed red box give a best estimate of the region scanned and meshed during MIT imaging, based upon scan dimensions given in the Visible Human documentation (Table 5.3 of documentation). Note that the estimated region includes four transverse processes or ribs joining to a thoracic vertebral body.

Fig. 12

Parasagittal CT slice through the upper torso, offset from the spinal axis by about 2.5 cm. The red box gives an estimate (see text) of the region probed and meshed for the MIT scan.

If a sagittal slice is taken as close to the spinal column centerline as possible, shown in Fig. 13, then the transverse processes disappear as they should, and just a darkened, low conductivity region is evident in the space of the spinal column. However, the spinal canal is not visible with this mesh, which offers 1 cm resolution, at best—note that the canal diameter is typically 1.5 cm and should start about 5 cm beneath the template (or $Z=2$ on the axis). If the smoothing parameter is reduced to 0.002, the canal is just barely visible, but not convincing. Thus, we proceed to consider results from mesh #1, which uses 10 layers of pentahedral elements and places more refinement at a greater depth.

Fig. 13

Sagittal slice (mesh #2) with subject in prone position, parallel to and centered along the spinal column—canal not visible.

To improve prospects for resolving the spinal canal with mesh #1, which is easily accomplished in MRI or CT scans, the orange threshold was moved from $1\mathrm{S}/\mathrm{m}$ down to $0.5\mathrm{S}/\mathrm{m}$. Other thresholds were left the same as before—black at 0.0, yellow at 2.0, and white for all values above $3.0\mathrm{S}/\mathrm{m}$, which are considered anomalous as stated before. Figure 14 shows that the canal is resolved in a transverse slice offset by about 2 cm from the template center when using mesh #1. The image shows that it is located at a depth of about 5.5 cm ($Z=1.5$) and has a diameter of about 1.0 cm. Directly above the canal, there is the hint of an outward spinal process on a vertebra. Also visible are the rib articulations coming into the column at an angle of about 45 deg. The absence of any resolved structure above the spinal column, closer to the skin surface is discussed in the next section. Again, to facilitate comparison of the structure visible in Fig. 14 with known standards, a transverse slice through a thoracic vertebral body is shown in Fig. 15, which is taken from full body CT scans of a human male cadaver (Visible Human, Research Systems, Inc., version 2.0) placed in a supine position, though shown here as prone.

Fig. 14

Transverse slice (mesh #1) with subject in prone position, about 2 cm below template center—spinal canal is visible at a depth of about 5.5 cm (1.5 cm on $Z$-axis), and centered horizontally.

Fig. 15

Prone-oriented transverse CT slice through a human male thoracic vertebral body; the red box is drawn to approximate the region that was scanned and meshed by the single-coil MIT method. Rib articulations with transverse processes are clearly visible here.

The final MIT image, shown in Fig. 16, is a sagittal slice parallel to and centered on the spinal column as nearly as possible, with the intent of resolving the canal over the full length of the space that was sampled. Clearly, in the regions beneath the template that are most visited by the coil (just offset from the template center), the canal is readily visible ($Z=1.5$), but is not fully visible along its entire length. This is understandable at the far right and left of the image where sampling frequency is lowest in the template, but better canal resolution was expected near the template center. Once improved methods are available for accurately positioning the coil relative to target, resolving structures in and around the spinal column should become more feasible. Note that at both the lower right and lower left of Fig. 16 the spinal column fades out, replaced by higher conductivity material. This is clearly in error due to poor sampling of these outlying regions.

Fig. 16

Sagittal slice (mesh #1) with subject in prone position, parallel to and centered along the spinal column; canal is visible along a portion of the spinal column along the line $Z=1.5\mathrm{cm}$—in particular, beneath those regions of the template that are better sampled.

8.

Discussion of Imaging Results

Spinal column images were essentially the same across all subjects. In all cases, ribs or transverse spinal processes were visualized in parasagittal slices on either side of the spinal column regardless of which mesh was used for inversion. Similarly, the spinal canal could be readily seen for all subjects, but was only seen clearly with mesh #1. Some variations across subjects were found, in the form of distortion—presumably the result of our inability to control and know coil rotation and position with sufficient accuracy. As far as we know, these are the first in vivo MIT images of electrical conductivity ever reported that visualize the interior of any portion of the human body. Electrical conductivity within the body (lungs) has been imaged using electrical impedance tomography, but resolution is limited to about 12% of the thoracic diameter,¹⁹ or about 40 mm. There are published results of MIT-imaged phantoms, composed of containers filled with various concentrations of saline. However, the body is much more challenging due to its greater spatial extent compared with the rather small phantoms commonly used.

A puzzling result associated with our set of spinal images is that while there is an ability to resolve objects at a significant depth beneath the template and skin, resolvable structures nearer the skin where the coil was actually placed seem absent. A part of the reason may be less conductivity contrast in the soft tissues around the spinal column than surmised, although it would seem that fat and muscle tissue should be distinguishable—fat conductivity is less than about $0.1\mathrm{S}/\mathrm{m}$ while muscle conductivity is nearly $1\mathrm{S}/\mathrm{m}$.¹ Given that virtual data were imaged successfully when conductivity varied by a factor of nearly 2 (Sec. 6), then surely fat and muscle tissues create sufficient contrast to permit their resolution. However, since all subjects had BMI $<25$, fat content between scapulae may be far too low, leaving essentially only muscle and bone. Given our crude scanning technique, it is also likely that structural resolution is deficient just beneath the skin since position error is much more likely to impact these locations. Revisiting Fig. 3 should make plain that a position error will be much more consequential at coil positions located just a few millimeters above the skin than when distances are around 15 mm. Observing the nurse while relocating the sensor indicated that precise positioning could never be assured with our scanning method since, among other reasons, the amount of pressure applied to place the device could not be kept the same, and even if possible, variations in upper back curvature would prevent knowing coil location and orientation with accuracy. This interpretation is consistent with the results of Figs. 9 and 10 which show that we actually achieved a more satisfactory optimal RMS fit error (lower) when most mesh nodes are farther from the coil locations—as with mesh #1.

In view of the impact of mesh design here, in particular mesh spatial refinement at greater depth, it would seem important to further investigate the value of spatially variant regularization,^20,21 which was shown effective for a related inverse problem solved in diffuse optical tomography (DOT). Permitting the reguarlization parameter to vary with increased distance between mesh nodes and sensor could be useful for penalizing solutions at nodes that are expected to contribute little to the loss signal measured at the coil. But it may be wiser to penalize solutions at nodes nearer the surface, where relative position of coil and surface is harder to establish at present (again, note position sensitivity as revealed in Fig. 3). At this point, whether to assign the regularization parameter to a higher value for nodes nearer or farther from the single sensor is not clear, but some experimentation with the approach is entirely appropriate. A difference between DOT and single coil MIT is that here we have but one sensor, while DOT uses at least two. Thus, spatially variant regularization may be simpler to implement and have greater effect for single coil MIT.

Clearly, sampling is vital to improving our prospects for single-coil MIT imaging. Ideally, additional internal (or external) instrumentation needs to capture pitch and roll rotations together with coil position, all while the instrument is moved within the space near the target. We envision that during collection of coil loss data, coil position and orientation are automatically collected and logged. A further benefit of position and orientation logging is the automatic determination of locations in space that are strictly outside of the body being scanned. Then, coil position and orientation data not only permit us to know where the coil is when coil loss is measured, such data also establish the boundary surface of the scanned object—since the coil only closely approaches the surface of the object, never pushing into it. Tracking coil position and orientation during a physical scan provides a powerful benefit to the inversion step—only the region of space not visited by the coil would ever need to be meshed, or discretized by grid construction. In addition to making more measurements, independent measurement sets could be created from sensors having different coil diameters. There is nothing to stop us from combining contributions from coils of different diameters in Eq. (8). This would certainly be helpful when there is a desire to probe structures at greater depth.