|Home | About | Journals | Submit | Contact Us | Français|
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 B1 which induces eddy current J in a conductive object volume (Xu & He 2005). The induced eddy current in an existing static magnetic field B0 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 V2 and an insulating region V1, while the region V2 is surrounded by the region V1. A current source Ji is applied inside V1 to generate stimulating magnetic field and eddy current induction in V2. From Maxwell's equations, the curl-curl equation for the electric field E inside the volume V can be derived as in (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)
Incorporating the magnetic vector potential A and the electric scalar potential ϕ, while B1 = × A and E = −A/t − ϕ, together with the Coulomb Gauge condition · A = 0, we can derive (3) in the conductive region V2
In addition, in quasi-static condition, the eddy current has its solenoidal nature, which can be described as in (4)
On the other hand, in the nonconducting region V1, the governing equation can be represented through Ampere's law by
where Ni are vector weighted functions and Ni are scalar weighted functions.
Using vector identities · (M × F) = ( × M) · F − M · ( × F) and · (rF) = (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)
where S is the boundary between V1 and V2 and S0 is the outer surface of V1, 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 Ni and Ni, the field variables A and can be discretized in each finite element and by assembling the residuals shown in equations (7a), (7b) and (7c) in all the elements and setting the weighted residual integral associated to each nodes to zero, the matrix equations of the field variables can be written as
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 V1,V2. 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
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)
where eKt is a state transition matrix. Using time-domain discretization and letting t0 = nΔt, t = (n + 1)Δt, Eq. 10 can be written as
The state transfer matrix has the property as in (12)
where M is a nonsingular matrix and U is the identity matrix. Assuming Ji1(t) is piecewise homogeneous and using the property of the state transfer matrix in (12), Eq. 11 becomes (Mohammed & Üler 1992)
Expanding the state transition matrix eKΔt and item (eKΔt − U)K−1Z in series and taking the lower order components, Eq. (13) can be used to recursively solve the state variable Xξ2 (Mohammed & Üler 1992). With the solved state variable, we can then obtain the magnetic flux density B1 and electric field distribution E in the medium. The induced eddy current density J can be further calculated using Ohm's law.
The induced eddy current density J(r,t) in a static magnetic field B0 experiences the Lorentz force which can be described as J(r,t) × B0. 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)
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
where Green's function G(r | r′) = ejk|r − r′|/4π|r − r′|. Equations (15) and (17) describe the pressure field in the time and frequency domain, respectively. With these two equations, given the acoustic source term ·(J×B0), the acoustic pressure field can be solved by summing the product of the acoustic source and the corresponding Green's function for all the elements over the object volume.
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×B0) inside the conductive sample can be further simplified as (×J)·B0 taking the fact that the static magnetic field B0 is generated from sources outside the conductive object, for example, from a permanent magnet and ×B0 = 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 B1, the curl of J is not zero. However, the acoustic source term (×J)·B0 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):
For a first-order tetrahedral element, as shown in Fig. 1, where J1, J2, J3, J4 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)
where ae,be,ce,de are vector parameters as and so on. They can be determined by the coordinates of the four tetrahedral vertices and the corresponding nodal solutions on them. Substitute formula (19) into (18), we have
Where Jxj, Jyj, Jzj are the three Cartesian components of the nodal values of the eddy currents, are the coefficients which can be determined from expansion of the determinants of elemental interpolation (Jin 2002). With formula (20), we can compute the divergence of Lorentz force, i.e. the acoustic source. In a homogeneous media, the left side of formula (20) can also be written as , we can use the relation of electric conductivity, magnetic field stimulation and static magnetic field to calculate the sound source. With calculated acoustic source distribution and assuming acoustic speed in the medium to be constant, we can then compute acoustic pressure on a given spatial point using formula (17).
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 R1 = 30mm and R2 = 30mm, respectively. The conductivity values of the inner and outer layer are σ1=0.25S/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.
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 pai and pbi, where the index a and b correspond to two different sets of model grid size. The sum of absolute difference of the two data sets can be calculated as in (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.
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.
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.
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.
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.
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.
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−9s. 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.