Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2841573

Formats

Article sections

Authors

Related links

Phys Med Biol. Author manuscript; available in PMC 2010 May 7.

Published in final edited form as:

Published online 2009 April 8. doi: 10.1088/0031-9155/54/9/005

PMCID: PMC2841573

NIHMSID: NIHMS169759

The publisher's final edited version of this article is available at Phys Med Biol

See other articles in PMC that cite the published article.

Magnetoacoustic Tomography with Magnetic Induction (MAT-MI) is a recently proposed imaging modality to image the electrical impedance of biological tissue. It combines the good contrast of electrical impedance tomography with the high spatial resolution of sonography. In this paper, three-dimensional MAT-MI forward problem was investigated using the finite element method (FEM). The corresponding FEM formulas describing the forward problem are introduced. In the finite element analysis, magnetic induction in an object with conductivity values close to biological tissues was first carried out. The stimulating magnetic field was simulated as that generated from a three-dimensional coil. The corresponding acoustic source and field were then simulated. Computer simulation studies were conducted using both concentric and eccentric spherical conductivity models with different geometric specifications. In addition, the grid size for finite element analysis was evaluated for model calibration and evaluation of the corresponding acoustic field.

Imaging electrical conductivity properties of biological tissues using non-invasive techniques has drawn much research interest during the last decades. These techniques include Electrical Impedance Tomography (EIT) (Paulson et al 1993, Metheral et al 1996), Magnetic Resonance Electrical Impedance Tomography (MREIT) (Kwon et al 2002, İder & Onart 2004, Gao et al 2005, Gao & He 2008) and Magnetic Induction Tomography (MIT) (Al-Zeibak et al 1995, Griffiths 2001). Among these methods, both EIT and MREIT need current injection into the sample from electrodes attached to its surface. In EIT electric potentials are measured on the surface of the sample to reconstruct its internal electrical impedance distribution. In MREIT magnetic flux density measurement data related to the injected current inside the sample are acquired from magnetic resonance phase images and these data are used to reconstruct the electrical conductivity distribution of the sample. In MIT, a dynamic magnetic field generated by an exciting coil is applied to induce eddy current in a conductive object. The magnetic field outside the object, perturbed by the eddy current, is detected to reconstruct the electromagnetic properties of the object. In addition, imaging modalities integrating electromagnetism and ultrasound technique were also introduced. Magnetoacoustic Tomography (MAT) was reported as a method to detect electrical current in biological tissue (Towe & Islam 1988). Another related technique, Hall Effect Imaging (HEI), was also reported to be able to image the electric properties of a sample (Wen et al 1998). In some operation models of MAT and HEI, electric current is injected into biological tissue that is placed in a static magnetic field and acoustic vibrations are generated inside tissue volume because of the Lorentz force. Acoustic signals are then measured around the object for inverse reconstruction (Roth et al 1994, Wen et al 1998).

Magnetoacoustic Tomography with Magnetic Induction (MAT-MI) (He 2005, Xu & He 2005, Li et al 2006, Xia et al 2007, Brinker & Roth 2008) is a new imaging modality recently proposed to image the electrical impedance of biological tissue. In this method, object is placed in a time-varying magnetic field and a static magnetic field. Eddy current is induced in the object by the time-varying magnetic field and as in MAT, acoustic vibrations are produced by the Lorentz force that comes from the interaction of the eddy current and the static magnetic field. Acoustic signal emitted from the object is collected to reconstruct the conductivity distribution within the sample object (Xu & He 2005). Compared to EIT, MREIT and MAT/HEI, the employment of magnetic excitation makes MAT-MI not affected by the “shielding effect” (Xu & He 2005) caused by low-conductivity tissue layers at the surface of the human body, such as the skull and the fat layer of the breast. In addition, other than MIT, MAT-MI utilizes the emitted megahertz ultrasound signal for conductivity reconstruction and is able to obtain high spatial resolution close to that of sonography.

In the present study, we solved the three-dimensional (3D) forward problem of MAT-MI with the aid of the finite element method (FEM). The corresponding FEM formulas describing the forward problem are introduced. A 3D coil model was used in the forward computation of acoustic fields. This kind of coil is easily implemented in practice and is related to a non-homogeneous magnetic stimulation. Magnetic induction in an object with conductivity values close to biological tissues was carried out using the FEM. The corresponding acoustic source, propagating field and pressure signal collected on the detecting surface were then calculated. Appropriate grid size for finite element analysis was determined according to the convergence in computing the acoustic field. Simulations with concentric and eccentric spherical conductive models with different geometric specifications were conducted to investigate the acoustic field and signal changes related to the conductivity distribution variations.

In MAT-MI, an excitation coil is loaded with electric current pulse, generating stimulating pulsed magnetic field **B**_{1} which induces eddy current **J** in a conductive object volume (Xu & He 2005). The induced eddy current in an existing static magnetic field **B**_{0} around the sample object subjects to the Lorentz force and consequently generates detectable acoustic vibrations. In principle, the signal generation mechanism of MAT-MI includes physical processes including magnetic induction, electromagnetic acoustic coupling through the Lorentz force and acoustic wave propagation. The governing equations and FEM formulas of these processes are described below.

The magnetic induction problem does not have an analytic solution unless with certain symmetric conductive geometry models (Li et al 2007) and needs finite element analysis in general. Here we used the state-space approach based on the magnetic vector potential and electric scalar potential (**A** − *ϕ*) method to derive the FEM solution (Mohammed & Üler 1992). Consider a whole volume *V* that consists of a conductive region *V*_{2} and an insulating region *V*_{1}, while the region *V*_{2} is surrounded by the region *V*_{1}. A current source **J*** _{i}* is applied inside

$$\nabla \times \left(\frac{1}{\mu}\nabla \times \mathbf{\text{E}}(\mathbf{\text{r}},t)\right)+\epsilon \frac{{\partial}^{2}\mathbf{\text{E}}(\mathbf{\text{r}},t)}{\partial {t}^{2}}+\sigma \frac{\partial \mathbf{\text{E}}(\mathbf{\text{r}},t)}{\partial t}=-\frac{\partial {\mathbf{\text{J}}}_{i}(\mathbf{\text{r}},t)}{\partial t}\phantom{\rule{2em}{0ex}}\mathbf{\text{r}}\in V$$

(1)

where *μ* is the permeability, *ε* is the permittivity and *σ* is the conductivity of the whole volume. In MAT-MI, as we are considering magnetic induction in biological tissue with mega-hertz frequency the induced displacement current is much smaller than the conduction current (Xu & He 2005). Consequently, the magnetic induction problem here can be considered quasi-static and the displacement current item in (1) can be neglected as shown in (2)

$$\nabla \times \left(\frac{1}{\mu}\nabla \times \mathbf{\text{E}}(\mathbf{\text{r}},t)\right)+\sigma \frac{\partial \mathbf{\text{E}}(\mathbf{\text{r}},t)}{\partial t}=-\frac{\partial {\mathbf{\text{J}}}_{i}(\mathbf{\text{r}},t)}{\partial t}\phantom{\rule{4.5em}{0ex}}\mathbf{\text{r}}\in V$$

(2)

Incorporating the magnetic vector potential **A** and the electric scalar potential *ϕ*, while **B**_{1} = × **A** and **E** = −**A**/*t* − *ϕ*, together with the Coulomb Gauge condition · **A** = 0, we can derive (3) in the conductive region *V*_{2}

$$\nabla \times \left(\frac{1}{\mu}\nabla \times \mathbf{\text{A}}(\mathbf{\text{r}},t)\right)+\sigma \left(\frac{\partial \mathbf{\text{A}}(\mathbf{\text{r}},t)}{\partial t}+\nabla \phi \right)=0\phantom{\rule{4.5em}{0ex}}\mathbf{\text{r}}\in {V}_{2}$$

(3)

In addition, in quasi-static condition, the eddy current has its solenoidal nature, which can be described as in (4)

$$\nabla \cdot \sigma \left(\frac{\partial \mathbf{\text{A}}(\mathbf{\text{r}},t)}{\partial t}+\nabla \phi \right)=0\phantom{\rule{4.5em}{0ex}}\mathbf{\text{r}}\in {V}_{2}$$

(4)

On the other hand, in the nonconducting region *V*_{1}, the governing equation can be represented through Ampere's law by

$$\nabla \times \left(\frac{1}{\mu}\nabla \times \mathbf{\text{A}}(\mathbf{\text{r}},t)\right)={\mathbf{\text{J}}}_{i}\phantom{\rule{4.5em}{0ex}}\mathbf{\text{r}}\in {V}_{1}$$

(5)

Using Galerkin's method for finite element analysis, the corresponding weighted residuals derived from (3), (4) and (5) are given by (Silvester and Ferrari 1996)

$${\mathit{\text{Res}}}_{2a}={\iiint}_{{V}_{2}}{\mathbf{\text{N}}}_{i}\cdot \nabla \times \left(\frac{1}{\mu}\nabla \times \mathbf{\text{A}}(\mathbf{\text{r}},t)\right)dV+{\iiint}_{{V}_{2}}\left\{{\mathbf{\text{N}}}_{i}\cdot \sigma \left(\frac{\partial \mathbf{\text{A}}(\mathbf{\text{r}},t)}{\partial t}+\nabla \phi \right)\right\}dV$$

(6a)

$${\mathit{\text{Res}}}_{2b}=-{\iiint}_{{V}_{2}}{N}_{i}\left\{\nabla \cdot \sigma \left(\frac{\partial \mathbf{\text{A}}(\mathbf{\text{r}},t)}{\partial t}+\nabla \phi \right)\right\}dV$$

(6b)

$${\mathit{\text{Res}}}_{1}={\iiint}_{{V}_{1}}\left\{{\mathbf{\text{N}}}_{i}\cdot \nabla \times \left(\frac{1}{\mu}\nabla \times \mathbf{\text{A}}(\mathbf{\text{r}},t)\right)-{\mathbf{\text{N}}}_{i}\cdot {\mathbf{\text{J}}}_{i}\right\}dV$$

(6c)

where **N*** _{i}* are vector weighted functions and

Using vector identities · (**M** × **F**) = ( × **M**) · **F** − **M** · ( × **F**) and · (*r***F**) = (*r*) · **F** + *r* · **F** together with the Gauss divergence theorem, the weak-form finite element formulation for (6a), (6b) and (6c) are given by (Silvester & Ferrari 1996)

$$\begin{array}{l}{\mathit{\text{Res}}}_{2a}={\iiint}_{{V}_{2}}\frac{1}{\mu}(\nabla \times {\mathbf{\text{N}}}_{i})\cdot (\nabla \times \mathbf{\text{A}}(\mathbf{\text{r}},t))dV-{\u222f}_{S}\frac{1}{\mu}{\mathbf{\text{N}}}_{i}\times (\nabla \times \mathbf{\text{A}}(\mathbf{\text{r}},t))\cdot \widehat{\mathbf{\text{n}}}dS\\ \phantom{\rule{4em}{0ex}}+{\iiint}_{{V}_{2}}{\mathbf{\text{N}}}_{i}\cdot \sigma \left(\frac{\partial \mathbf{\text{A}}(\mathbf{\text{r}},t)}{\partial t}+\nabla \phi \right)dV\end{array}$$

(7a)

$${\mathit{\text{Res}}}_{2b}={\iiint}_{{V}_{2}}\left\{\nabla {N}_{i}\xb7\sigma \left(\frac{\partial \mathbf{\text{A}}(\mathbf{\text{r}},t)}{\partial t}+\nabla \phi \right)\right\}dV-{\u222f}_{S}{N}_{i}\sigma \left(\frac{\partial \mathbf{\text{A}}(\mathbf{\text{r}},t)}{\partial t}+\nabla \phi \right)\cdot \widehat{\mathbf{\text{n}}}dS$$

(7b)

$$\begin{array}{l}{\mathit{\text{Res}}}_{1}={\iiint}_{{V}_{1}}\left\{\frac{1}{\mu}(\nabla \times {\mathbf{\text{N}}}_{i})\cdot (\nabla \times \mathbf{\text{A}}(\mathbf{\text{r}},t))-{\mathbf{\text{N}}}_{i}\cdot {\mathbf{\text{J}}}_{i}\right\}dV-{\u222f}_{S}\frac{1}{\mu}{\mathbf{\text{N}}}_{i}\times (\nabla \times \mathbf{\text{A}}(\mathbf{\text{r}},t))\cdot \widehat{\mathbf{\text{n}}}dS\\ \phantom{\rule{4em}{0ex}}-{\u222f}_{{S}_{0}}\frac{1}{\mu}{\mathbf{\text{N}}}_{i}\times (\nabla \times \mathbf{\text{A}}(\mathbf{\text{r}},t))\cdot \widehat{\mathbf{\text{n}}}dS\end{array}$$

(7c)

where *S* is the boundary between *V*_{1} and *V*_{2} and *S*_{0} is the outer surface of *V*_{1}, represents the vector normal to the boundary surfaces. In Equation (7a) and (7b) we also employed a transform *ϕ* = /*t* to make the later introduced coefficient matrices of the discrete algebraic equations symmetric.

Using the basis functions similar as the weighted functions **N*** _{i}* and

$$\left[\begin{array}{cc}{\gamma}_{11}& {\gamma}_{12}\\ {\gamma}_{21}& {\gamma}_{22}\end{array}\right]\phantom{\rule{0.2em}{0ex}}\left[\begin{array}{c}{\mathbf{\text{A}}}_{\xi 1}\\ {\left[{\mathbf{\text{A}}}_{\xi 2},\nabla \varphi \right]}^{\text{T}}\end{array}\right]=\left[\begin{array}{c}0\\ {\left[{\dot{\mathbf{\text{A}}}}_{\xi 2},\nabla \dot{\varphi}\right]}^{\text{T}}\end{array}\right]-\left[\begin{array}{c}{\mathbf{\text{J}}}_{i1}\\ 0\end{array}\right]$$

(8)

where *γ*_{11}, *γ*_{12} *γ*_{21}, *γ*_{22} are the coefficient matrices that can be determined by the residuals in all the elements and the selected weighted and basis functions. Note here the first row and the second row of (8) do not have the same units and the units of *γ*_{11}, *γ*_{12} are different from those of *γ*_{21}, *γ*_{22}. The subscript *ξ* denotes the *x*, *y*, *z* component of magnetic vector potential **A**, the subscripts 1, 2 denote the variants and coefficients for *V*_{1},*V*_{2}. and are the first-order derivatives of **A** and in the time domain.

Using elimination method and let **X**_{ξ2} = [**A**_{ξ2}, ]^{T}, we have

$${\dot{\mathbf{\text{X}}}}_{\xi 2}=\mathbf{\text{K}}{\mathbf{\text{X}}}_{\xi 2}+\mathbf{\text{Z}}{\mathbf{\text{J}}}_{i1}$$

(9)

where **K** and **Z** are determined by the coefficient matrices in (8). If the electric and magnetic parameters are not changing over time, the matrices **K** and **Z** become constant. In that case, Eq. (9) can be considered as a linear time invariant state equation. The corresponding solution to (9) is (Skelton 1988)

$${\mathbf{\text{X}}}_{\xi 2}(t)={e}^{\mathbf{\text{K}}(t-{t}_{0})}{\mathbf{\text{X}}}_{\xi 2}({t}_{0})+{\int}_{{t}_{0}}^{t}{e}^{\mathbf{\text{K}}(t-\tau )}\mathbf{\text{Z}}{\mathbf{\text{J}}}_{i1}(\tau )d\tau $$

(10)

where *e*^{K}* ^{t}* is a state transition matrix. Using time-domain discretization and letting

$${\mathbf{\text{X}}}_{\xi 2}((n+1)\mathrm{\Delta}t)={e}^{\mathbf{\text{K}}\mathrm{\Delta}t}{\mathbf{\text{X}}}_{\xi 2}(n\mathrm{\Delta}t)+{\int}_{n\mathrm{\Delta}t}^{(n+1)\mathrm{\Delta}t}{e}^{\mathbf{\text{K}}((n+1)\mathrm{\Delta}t-\tau )}\mathbf{\text{Z}}{\mathbf{\text{J}}}_{i1}(\tau )d\tau $$

(11)

The state transfer matrix has the property as in (12)

$${\int}_{0}^{t}{e}^{\mathbf{\text{M}}\tau}d\tau ={\mathbf{\text{M}}}^{-1}\left({e}^{\mathbf{\text{M}}t}-\mathbf{\text{U}}\right)=\left({e}^{\mathbf{\text{M}}t}-\mathbf{\text{U}}\right){\mathbf{\text{M}}}^{-1}$$

(12)

where **M** is a nonsingular matrix and **U** is the identity matrix. Assuming **J**_{i}_{1}(*t*) is piecewise homogeneous and using the property of the state transfer matrix in (12), Eq. 11 becomes (Mohammed & Üler 1992)

$${\mathbf{\text{X}}}_{\xi 2}((n+1)\mathrm{\Delta}t)={e}^{\mathbf{\text{K}}\mathrm{\Delta}t}{\mathbf{\text{X}}}_{\xi 2}(n\mathrm{\Delta}t)+({e}^{\mathbf{\text{K}}\mathrm{\Delta}t}-\mathbf{\text{U}}){\mathbf{\text{K}}}^{-1}\mathbf{\text{Z}}{\mathbf{\text{J}}}_{i1}(n\mathrm{\Delta}t)$$

(13)

Expanding the state transition matrix *e*^{K}^{Δ}* ^{t}* and item (

The induced eddy current density **J**(**r**,*t*) in a static magnetic field **B**_{0} experiences the Lorentz force which can be described as **J**(**r**,*t*) × **B**_{0}. The time-varying Lorentz force induces acoustic vibrations in the object and the wave equation for acoustic pressure *p*(**r**,*t*) is given in (14) (Roth et al 1994)

$${\nabla}^{2}p(\mathbf{\text{r}},t)-\frac{1}{{c}_{s}^{2}}\frac{{\partial}^{2}p(\mathbf{\text{r}},t)}{\partial {t}^{2}}=\nabla \cdot (\mathbf{\text{J}}(\mathbf{\text{r}},t)\times {\mathbf{\text{B}}}_{0})$$

(14)

where *c _{s}* is the acoustic speed in the media. Using Green's function, the solution of (14) can be written as (Morse & Ingard 1968)

$$p(\mathbf{\text{r}},t)=-\frac{1}{4\pi}{\iiint}_{V}d\mathbf{\text{r}}\prime {\int}_{-\infty}^{+\infty}\nabla \cdot (\mathbf{\text{J}}(\mathbf{\text{r}},t\prime )\times {\mathbf{\text{B}}}_{0})\frac{\delta \left(t-t\prime -\frac{R}{{c}_{s}}\right)}{R}dt\prime $$

(15)

where R = |**r** − **r**′| and the integration volume *V* covers the entire conductive object.

Taking Fourier transform on both hand sides of equation (14), the wave equation in frequency-domain is then given as

$${\nabla}^{2}p(\mathbf{\text{r}},\omega )+{k}^{2}p(\mathbf{\text{r}},\omega )=\nabla \cdot (\mathbf{\text{J}}(\mathbf{\text{r}},\omega )\times {\mathbf{\text{B}}}_{0})$$

(16)

where k = *ω*/*c _{s}* is the wave number and

$$p(\mathbf{\text{r}},\omega )=-{\iiint}_{V}\nabla \cdot (\mathbf{\text{J}}(\mathbf{\text{r}}\text{'},\omega )\times {\mathbf{\text{B}}}_{0})G(\mathbf{\text{r}}|\mathbf{\text{r}}\text{'})d{\mathbf{\text{r}}}^{\prime}$$

(17)

where Green's function *G*(**r** | **r**′) = *e ^{jk}*

From the finite element analysis of the magnetic induction problem mentioned above, we can evaluate the eddy current **J** on each element inside the conductive object volume. The acoustic source term ·(**J**×**B**_{0}) inside the conductive sample can be further simplified as (×**J**)·**B**_{0} taking the fact that the static magnetic field **B**_{0} is generated from sources outside the conductive object, for example, from a permanent magnet and ×**B**_{0} = 0 inside the object volume (Xu & He, 2005). Note here as the eddy current **J** inside the conductive volume is induced by the time varying magnetic field **B**_{1}, the curl of **J** is not zero. However, the acoustic source term (×**J**)·**B**_{0} contains spatial derivative of the eddy current and thus in order to solve the acoustic source in each element, we need make further treatment on the FEM calculated eddy current data. Actually this term can be expanded inside the object volume as in (18):

$$\nabla \cdot (\mathbf{\text{J}}\times {\mathbf{\text{B}}}_{0})={\text{B}}_{0z}\frac{\partial {\text{J}}_{y}}{\partial x}-{\text{B}}_{0y}\frac{\partial {\text{J}}_{z}}{\partial x}+{\text{B}}_{0x}\frac{\partial {\text{J}}_{z}}{\partial y}-{\text{B}}_{0z}\frac{\partial {\text{J}}_{x}}{\partial y}+{\text{B}}_{0y}\frac{\partial {\text{J}}_{x}}{\partial z}-{\text{B}}_{0x}\frac{\partial {\text{J}}_{y}}{\partial z}$$

(18)

For a first-order tetrahedral element, as shown in Fig. 1, where **J**_{1}, **J**_{2}, **J**_{3}, **J**_{4} are the nodal solutions of induced eddy current density, using finite element interpolation, the eddy current density **J**(*x*, *y*, *z*) on a point in the element can be written as in (19) (Jin 2002)

$$\mathbf{\text{J}}(x,y,z)={\mathbf{\text{a}}}^{e}+{\mathbf{\text{b}}}^{e}x+{\mathbf{\text{c}}}^{e}y+{\mathbf{\text{d}}}^{e}z$$

(19)

where **a*** ^{e}*,

$$\begin{array}{l}\nabla \cdot (\mathbf{\text{J}}\times {\mathbf{\text{B}}}_{0})={\text{B}}_{0z}{\text{b}}_{y}^{e}-{\text{B}}_{0y}{\text{b}}_{z}^{e}+{\text{B}}_{0x}{\text{c}}_{z}^{e}-{\text{B}}_{0z}{\text{c}}_{x}^{e}+{\text{B}}_{0y}{\text{d}}_{x}^{e}-{\text{B}}_{0x}{\text{d}}_{y}^{e}\\ \phantom{\rule{5em}{0ex}}=\frac{1}{6{V}^{e}}({\text{B}}_{0z}\sum _{j=1}^{4}{\text{b}}_{j}^{e}{J}_{yj}-{\text{B}}_{0y}\sum _{j=1}^{4}{\text{b}}_{j}^{e}{J}_{zj}+{\text{B}}_{0x}\sum _{j=1}^{4}{\text{c}}_{j}^{e}{J}_{zj}\\ \phantom{\rule{7em}{0ex}}-{\text{B}}_{0z}\sum _{j=1}^{4}{\text{c}}_{j}^{e}{J}_{xj}+{\text{B}}_{0y}\sum _{j=1}^{4}{\text{d}}_{j}^{e}{J}_{xj}-{\text{B}}_{0x}\sum _{j=1}^{4}{\text{d}}_{j}^{e}{J}_{yj})\end{array}$$

(20)

Where *J _{xj}*,

In previous studies of MAT-MI, numerical simulation was conducted on symmetrical conductivity models with pulsed and spatially homogeneous magnetic stimulation (Xu & He 2005, Li et al 2006). In the present study, using finite element analysis, we employed a real shape excitation coil model and conducted simulations on two-layer concentric and eccentric spherical conductivity models. The excitation coil produces a more realistic but on the other hand more complex magnetic field than the homogeneous one, allowing us to simulate more sophisticated eddy current distribution in the object volume.

In the present simulation study, induced eddy current within the conductivity model was computed based on the electromagnetic finite element analysis and implemented using the FEM software ANSYS. The acoustic source is calculated in every tetrahedral element according to the interpolation formula (20). Pressure signals on the detection surface were calculated using the integral formula of (17).

The present simulation model includes a two-layer conductive sphere, a real shape coil, and the surrounding water and air media models. The water media surrounding the conductive sphere has a cylindrical shape with radius of 300mm and height of 600mm. The water media is surrounded by air media with radius of 450mm and height of 800mm. The real-shape coil model with a concentric and an eccentric two-layer spherical conductivity model are shown in Fig. 2. Fig 2(a) shows the concentric spherical model, together with the coil model. The coil has a height of 10 mm, with inner radius and outer radius to be 35mm and 45mm, respectively. The axis of the coil is in Z direction, while the center of the coil is 110mm away from the center of the conductive sphere, which resides at the origin. The coil is excited by current pulses to produce magnetic stimulation, with 20 *μs* pulse width and 10 *kHz* repetition frequency. Simulation was done with 200 *kHz* sampling frequency. The strength of the magnetic stimulation at the edge of the coil is 0.1T. The static magnetic field is presumed to be homogeneous with magnetic flux density to be 1T. For the conductivity models, the radii of the inner and outer layer sphere are *R*_{1} = 30*mm* and *R*2 = 30*mm*, respectively. The conductivity values of the inner and outer layer are *σ*_{1}=0.25*S*/*m* and *σ*_{2}=0.04 *S*/*m*, respectively. In all the simulations, these parameters of spherical radii and conductivity values were fixed. Fig 2(b) shows the eccentric spherical model together with the surrounding spherical detecting surface where acoustic signals were simulated. The detecting surface includes those potential positions where we may put transducers on in experiment system. The radius of the detection surface was set to be 140 *mm*. We built the solid models with ANSYS and meshed them with tetrahedrons. The meshed grids of the coil model and a two-layer spherical conductivity model are shown in Fig. 3. The mesh of water media and air media are fixed in all the simulation. The numbers of nodes and elements for water media are 883 and 5732, the numbers of nodes and elements for air media are 3969 and 22565, respectively.

(a) A concentric spherical conductivity model together with the realistic geometry coil model. (b) An eccentric spherical model

In order to determine the appropriate finite element size employed in the conductive model for the forward calculation, we tested its convergence property. The convergence test is also an indirect validation of the forward calculation. Specifically, with the grid size of the coil and the driving current parameters unchanged, while refining the grid size of the two-layer conductive sphere, we compared the corresponding acoustic signal sets calculated on the detection surface surrounding the object. When the calculated acoustic signal set converges with the grid size becoming smaller, we consider the grid size to be appropriate with which the FEM model produces reliable and precise enough solution.

With N sampling positions on the detection surface surrounding the object, the acoustic pressure signal on the *i* th sampling position at certain time point can be calculated as *p _{ai}* and

$${D}_{ab}=\sum _{i=1}^{N}\left|{p}_{ai}-{p}_{bi}\right|$$

(21)

If this difference decreases as we refining the model grids, we consider the forward calculation converging. In addition, by this way, we can find out the right FEM model with appropriate grid size that produces small enough error and reliable solution.

- Perform finite element analysis on the magnetic induction problem and calculate the induced eddy current in the conductive model with
*M*tetrahedral elements. - Calculate acoustic source in each finite element based on Eq. (20).
- Perform Fourier transform to acoustic source signals on every element with
*N*discrete sampling time points. - For every discrete frequency, use (22), the discrete form of the integral formula (17), to calculate the acoustic pressure at a given spatial point
**r**$$p(\mathbf{\text{r}},{\omega}_{i})=-\sum _{j=1}^{M}(\nabla \cdot ({\mathbf{\text{J}}}_{j}({\mathbf{\text{r}}}_{j}\text{'},{\omega}_{i})\times {\mathbf{\text{B}}}_{0}))G(\mathbf{\text{r}}|{\mathbf{\text{r}}}_{j}^{\prime})$$(22)where*ω*_{i}is the*i*th discrete frequency,**r**′ is the geometrical center of the_{j}*j*th tetrahedral element. Here, we have*N*complex numbers in frequency domain for each spatial point. - Perform inverse Fourier transform on the sequence in frequency domain calculated in step (4). This gives a temporal acoustic pressure sequence with
*N*discrete points.

In the present study, we conducted simulations in one concentric spherical model and five eccentric spherical models. The center coordinates of the inner and the outer layer spheres for both concentric and eccentric models are shown in Table 1.

With the aforementioned coil model, we conducted the finite element analysis of magnetic induction on both concentric and eccentric spherical models. A simulation example is shown in Fig. 4. The left part is the result using the concentric spherical model while the right part is the result using the X-29.5 eccentric model. In both cases, the magnitudes of magnetic flux density and induced eddy current distribution are illustrated in Z=0 and X=0 plane.

Magnitudes of magnetic flux density and eddy current density distributions on Z=0 and X=0 plane calculated with concentric and eccentric spherical models.

We used three sets of different grid size in the two-layer concentric model for calibration and determination of the appropriate grid size. The element and node numbers in the inner and outer layers corresponding to these three settings are listed in Table 2.

Number of elements and nodes in inner and outer layer spheres with three different grid size settings

With these three settings acoustic pressure signals were calculated on 5298 points uniformly distributed on the spherical detection surface. Using formula (21), the sum of absolute difference between each two acoustic pressure data sets were calculated and shown in Table 3.

The sum of absolute difference between each two pressure data sets corresponding to the three grid size settings

As shown in Table 3, the difference between the pressure data sets *b* and *c* is much smaller than those between *a*,*c* and *a*,*b*, indicating the trend of convergence of the forward calculation relative to the fineness of the grid. The grid size *b* is then considered to be fine enough for acoustic computation for MAT-MI method.

With the determined appropriate grid size for forward calculation, we meshed the concentric and eccentric models with almost the same grid size. The number of elements and nodes in the inner and outer spheres for all the models used in the present study are shown in Table 4.

In order to study the acoustic field and pressure signal variation under different conductivity distributions, we compared the simulation results using the concentric and eccentric models. The acoustic field was calculated over the three dimensional spatial domain at different time points, with the field at 5, 25, 35, 45, 55 *μs* after magnetic stimulation illustrated.

Fig. 5 shows the acoustic pressure distribution on planes Z=50mm, 25mm, 5mm, -5mm, -25mm, -50mm calculated with the concentric model. As in Fig. 5, the acoustic pressure distribution in every slice perpendicular to Z axis is centrosymmetric. This is because of the symmetry of the concentric conductivity model and the cylindrical symmetry of the magnetic stimulation generated by the coil. From Fig. 5 the acoustic propagation can also be clearly seen as the acoustic field evolves over time. In addition, it is clear that after magnetic stimulation, the conductive volume closer to the coil has acoustic source with larger amplitude than the volume further away from the coil.

Acoustic pressure distribution in units of Pascal (Pa) on plane z=50mm, z=25mm, z=5mm, z=-5mm, z=-25mm, z=-50mm at time 5, 25, 35, 45, 55 *μs* after magnetic stimulation with concentric spherical model.

The acoustic pressure distribution on Z=0 plane obtained from all six models listed in Table 1 are shown in Fig. 6. From Fig. 6, we can see that in MAT-MI, conductive volumes with higher conductivity values and located closer to the simulating coil produce acoustic source with higher strength and play dominant roles in the acoustic signal generation. In addition, the farther the eccentric sphere is away from the symmetric axis, the lower strength of the acoustic source is produced.

Fig. 7 shows the pressure distribution on Y=0 plane for concentric model, Z+10 and Z-10 models. The acoustic pressure is not centrosymmetric on Y=0 plane even for the concentric model. This is because the magnetic stimulation generated by the real shape coil becomes smaller as distance from the coil increases. In consequence, as shown in Fig. 7, the closer the inner sphere to the coil, the larger the acoustic source can be induced.

Acoustic pressure in units of Pascal (Pa) on Y=0 plane at time 5, 25, 35, 45, 55 *μs* after magnetic stimulation for concentric model and eccentric models Z+10 and Z-10 when the homogeneous static magnetic field is in Z direction.

As in Fig. 5 to Fig. 7, we note that there are some inhomogeneous distributions shown in the calculated acoustic fields 5 *μs* after the stimulation. This is because in the simulation, we assume the sound source in each element is a point source located at the element center. The amplitude of the point source in each element is set to be proportional to the element volume. Within elements having different sizes, the corresponding point source amplitudes are then different. Consequently, the calculated pressure is only accurate when the distance R=|**r**-**r**′| is large enough compared to the element size. That is the reason the acoustic field calculated shortly after the stimulation has those inhomogeneous appearance. This problem, however, can be solved by using smaller element.

Finally, we compared the acoustic signals simulated on six cross points between the detection surface and X, Y and Z axis. The coordinates of these six points are (-140, 0, 0), (0, -140, 0), (0, 0, -140), (140, 0, 0), (0, 140, 0) and (0, 0, 140) mm, respectively. The time sequences of the acoustic signal on the six points obtained from six different models (concentric, eccentric X-5, X-15, X-29.5, Z+10, Z-10) are shown in Fig. 8. For concentric model, due to its symmetry in conductivity distribution and the cylindrical symmetric magnetic stimulation, the acoustic signal detected at points on X and Y axis are similar as shown in Fig. 8 (a), (b), (c) and (d). In addition, as shown in Fig. 8 (a) and (b), the acoustic signals collected on X axis are sensitive to the conductivity model changes that happens on X axis. As the inner sphere getting closer to the detecting point, the peak time of the curve appears earlier, and vice versa. However, this kind of signal variation is insensitive to those model changes happening on Z axis. In a similar sense, from Fig. 8 (c) and (d), we can see that the acoustic signals detected at points on Y axis are similar for all six models because of their symmetry in Y axis. Furthermore, as shown in Fig. 8 (e) and (f), the peak time of acoustic signal curve for concentric and eccentric models X-5, X-15, X-29.5 are similar, but for eccentric models Z+10 and Z-10, the peak time of the curve are different. For models with inner sphere closer to the detecting point, the peak time appears earlier.

In the present study, we have conducted 3D numerical simulation of MAT-MI forward problem using the finite element method. Compared to the homogeneous spherical model (Xu & He 2005) and concentric spherical model (Li et al 2007) used in pervious simulation studies, the FEM method developed here makes us capable of doing MAT-MI forward simulation with arbitrary geometry model. Simulations using both concentric and eccentric conductivity models have been conducted and the effects of conductivity asymmetric distribution on the generated acoustic field and wave propagation have been investigated.

The aim of MAT-MI is to image electrical conductivity distribution of biological tissue. It is important to analyze the electromagnetic field and induced acoustic filed in MAT-MI by solving the MAT-MI forward problem. Some researchers have studied the induced electric field in biological tissue with Transcranial Magnetic Stimulation (TMS) (Miranda et al 2003) and Magnetic Induction Tomography (MIT) (Griffiths 2001). Here we performed finite element electromagnetic analysis with a real shape coil and spherical conductivity models. MAT-MI acoustic sources were evaluated based on this analysis. Acoustic calibration on the calculated field with different grid sizes was used to make sure that the acoustic field solution is convergent outside the object. In the present FEM analysis, we used first-order tetrahedral finite element and the distribution of the acoustic signal were interpolated and set constant inside each finite element. In order to have high accuracy of finite element solution, higher-order finite elements can be introduced to MAT-MI electromagnetic analysis.

In the derivation of the FEM formula, we considered magnetic diffusion in general as shown in the first item in Eq. (3). However, as in MAT-MI we are considering imaging biological tissue with around Megahertz stimulation, the magnetic diffusion is actually ignorable. For example, as for a conductive sphere with radius to be 60 mm and conductivity value 0f 0.25 *S* / *m*, the resistive diffusion time is calculated to be 4.52×10^{−9}*s*. Compared to the time step corresponding to 200 *kHz* sampling frequency, which is 5 *μs*, the magnetic diffusion time can be ignored. From the skin effect point of view, we note that the size of our conductivity model is much smaller than the magnetic skin depth and the skin effect can be ignored. This is also confirmed by our simulation result that the magnetic field caused by the induced eddy current is 10E-5 times less than the main stimulating magnetic field.

As shown in the simulation results, in MAT-MI, the tissue with higher conductivity value produces higher strength of sound source and so plays a bigger role in the acoustic signal generation. This property suggests that MAT-MI has potential to locate those tissue or organs, which have lower impedance value than their surrounding tissue. As there are both in vivo and in vitro studies showing that the conductivity values of malignant tumor tissue is significantly larger than the surrounding normal tissues, (Surowiec et al 1988, Campbell & Land 1992, Jossinet 1996, Morimoto et al 1993), MAT-MI has the potential to become a noninvasive method for tumor screening and detection. From the present simulation results, we can also see that the acoustic pressure distribution in spatial domain has a close relation to the conductivity distribution, and the acoustic pressure sequence measured on detection surface is sensitive to the change of conductivity distribution inside the object. This indicates that the detected acoustic signals contain information for inverse conductivity reconstruction.

The present simulation is performed under the acoustic homogeneity assumption, while the acoustic speed is set to be constant (1500m/s). This assumption certainly simplifies the simulation algorithm and limits the application of current algorithm to imaging soft biological tissue. In fact, the acoustic heterogeneity in soft tissue is within 10%, and in MAT-MI the wave reflection and scattering is negligible, the distortion caused by acoustic heterogeneity is small (Xu & He 2005).

Also note that relatively simple volume conductor models (concentric and eccentric spheres) were used in the present simulation, for the purpose of evaluating the developed FEM algorithm. To our knowledge, the developed FEM algorithm, for the first time, solves the forward problem of MAT-MI in a conductive object with arbitrary geometry together with realistic geometry excitation coil. Such rigorous forward solution is important for accurately evaluating MAT-MI imaging and assessing parameters of importance in MAT-MI. The eccentric spherical models allow us to simulate the acoustic fields in a well-controlled computational setting and provide evaluation of the developed FEM algorithm. Also, the present simulation provides detailed results of simulated acoustic fields in the MAT-MI, which allows us to correlate them with experimental recordings.

In summary, we have developed a forward solver for MAT-MI by using the finite element method, and conducted computer simulations for MAT-MI forward problem with a realistic geometry excitation coil. We have simulated the acoustic sources and fields using concentric and eccentric spherical conductivity models, and examined the relations between acoustic pressure and conductivity distributions. The MAT-MI forward solutions enable us to quantitatively evaluate the sensitivity of acoustic signals on the detection surface to the change of electrical conductivity within an object.

The authors are grateful to the anonymous reviewers for the constructive comments to the initial version of manuscript. This work was supported in part by NIH R21EB006070 and NSF BES-0602957 to BH, and partly supported by the Supercomputing Institute of the University of Minnesota.

- Al-Zeibak S, Goss D, Lyon G, Yu ZZ, Peyton AJ, Beck MS. A feasibility study of electromagnetic inductance tomography. Proc. 9th Int. Conf. on Electrical Bio-Impedance (Heidelberg, 1995); 1995. pp. 426–9.
- Brinker K, Roth BJ. The effect of electrical anisotropy during magnetoacoustic tomography with magnetic induction. IEEE Trans Biomed Eng. 2008;55:1637–9. [PubMed]
- Campbell AM, Land DV. Dielectric properties of female breast tissue measured in vitro at 3.2 GHz. Phys Med Biol. 1992;37:193–210. [PubMed]
- Gao N, Zhu S, He B. Estimation of electrical conductivity distribution within the human head from magnetic flux density measurement. Phys Med Biol. 2005;50:2675–87. [PubMed]
- Gao N, He B. Noninvasive Imaging of Bioimpedance Distribution by means of Current Reconstruction Magnetic Resonance Electrical Impedance Tomography. IEEE Trans Biomed Eng. 2008;55:1530–8. [PubMed]
- Griffiths H. Magnetic induction tomography. Meas Sci Technol. 2001;12:1126–31.
- He B. High-resolution Functional Source and Impedance Imaging. Proc. Of Annual Int. Conf. of IEEE-EMBS; 2005. pp. 4178–82. [PubMed]
- İder YZ, Onart S. Algebraic reconstruction for 3D magnetic resonance-electrical impedance tomography (MREIT) using one component of magnetic flux density. Phys Meas. 2004;25:281–94. [PubMed]
- Jin J. The Finite Element Method in Electromagnetics. 2rd. New York: John Wiley & Sons; 2002.
- Jossinet J. Variability of impedivity in normal and pathological breast tissue. Med Biol Eng Comput. 1996;34:346–50. [PubMed]
- Kwon O, Woo EJ, Yoon JR, Seo JK. Magnetic resonance electrical impedance tomography(MREIT): Simulation Study of J-Substitution Algorithm. IEEE Trans Biomed Eng. 2002;49:160–7. [PubMed]
- Li X, Xu Y, He B. Magnetoacoustic tomography with magnetic induction for imaging electrical impedance of biological tissue. J Appl Phys. 2006;99:066112. [PMC free article] [PubMed]
- Li X, Xu Y, He B. Imaging Electrical Impedance From Acoustic Measurements by Means of Magnetoacoustic Tomography With Magnetic Induction (MAT-MI) IEEE Trans Biomed Eng. 2007;54:323–30. [PubMed]
- Metheral P, Barber DC, Smallwood RH, Brown BH. Three dimensional electrical impedance tomography. Nature. 1996;380:509–12. [PubMed]
- Miranda PG, Hallett M, Basser PJ. The Electric Field induced in the Brain by Magnetic Stimulation: A 3-D Finite-Element Analysis of the Effect of Tissue Heterogeneity and Anisotropy. IEEE Trans Biomed Eng. 2003;50:1074–84. [PubMed]
- Mohammed OA, Üler FG. A State Space Approach and Formulation for the Solution of Nonlinear 3-D Transient Eddy Current Problems. IEEE Trans Magnetics. 1992;28:1111–4.
- Morimoto T, Kimura S, Konishi Y, Komaki K, Uyama T, Monden Y, Kinouchi DY, Iritani DT. A study of the electrical bioimpedance of tumors. J Invest Surg. 1993;6:25–32. [PubMed]
- Morse PM, Ingard KU. Theoretical Acoustics. McGraw-Hill; 1968.
- Paulson K, Lionheart W, Pidecock M. Optimal experiments in electrical-impedance tomography. IEEE Trans Med Imag. 1993;12:681–6. [PubMed]
- Roth BJ, Basser PJ, Wikswo JP., Jr A Theoretical Model for Magneto-Acoustic Imaging of Bioelectric Currents. IEEE Trans Biomed Eng. 1994;41:723–8. [PubMed]
- Silvester PP, Ferrari RL. Finite elements for electrical engineers. 3rd. Cambridge: Cambridge University Press; 1996.
- Skelton RE. Dynamic systems control: linear systems analysis and synthesis. John Wiley & Sons; 1988.
- Surowiec AJ, Stuchly SS, Barr JB, Swarup A. Dielectric properties of breast carcinoma and the surrounding tissues. IEEE Trans Biomed Eng. 1988;35:257–63. [PubMed]
- Towe BC, Islam MR. A Magneto-Acoustic Method for the Noninvasive Measurement of Bioelectric Currents. IEEE Trans Biomed Eng. 1988;35:892–4. [PubMed]
- Wen H, Shah J, Balaban RS. Hall Effect Imaging. IEEE Trans Biomed Eng. 1998;45:119–24. [PMC free article] [PubMed]
- Williams EG. Fourier Acoustics: Sound Radiation and Nearfield Acoustical Holography. London: Academica press; 1999.
- Xia R, Li X, He B. Magnetoacoustic tomographic imaging of electrical impedance with magnetic induction. Appl Phys Lett. 2007;91:083903. [PMC free article] [PubMed]
- Xu Y, He B. Magnetoacoustic tomography with magnetic induction (MAT-MI) Phys Med Biol. 2005;50:5175–87. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |