Reconstruction of Conductivity Distribution of Brain Tissue from Two Components Magnetic Flux Density

In this paper the recent Magnetic resonance electrical impedance imaging (MREIT) technique is used to image non-invasively the three-dimensional continuous conductivity distribution of the head tissues. With the feasibility of the human head being rotated twice in the magnetic resonance imaging (MRI) system, a continuous conductivity reconstruction MREIT algorithm based on two components of the measured magnetic flux density is introduced. The reconstructed conductivity image could be obtained through solving iteratively a non-linear matrix equation. According to the present algorithm of using two magnetic flux density components, numerical simulations were performed on a concentric three-sphere and realistic human head model (consisting of the scalp, skull and brain) with the uniform and non-uniform isotropic target conductivity distributions. Based on the algorithm , the reconstruction of scalp and brain conductivity ratios could be figured out even under the condition that only one current is injected into the brain. The present results show that the three-dimensional continuous conductivity reconstruction method with two magnetic flux density components for the realistic head could get better results than the method with only one magnetic flux density component. Given the skull conductivity ratio, the relative errors of scalp and brain conductivity values were reduced to less than 1% with the uniform conductivity distribution and less than 6.5% with the non-uniform distribution for different noise levels. Furthermore, the algorithm also shows fast convergence and improved robustness against noise.


INTRODUCTION
Knowledge of the electrical conductivity distribution in the human body is of great importance in many biomedical applications.The information of the electrical properties of head tissues is used for electroencephalographic source location and functional mapping of brain activity and functions.And it has been proved that information about the vivo tissue conductivity values improves the solution accuracy of bioelectrical field problems [1].
A non-invasive Magnetic Resonance Electrical Impedance Tomography (MREIT) imaging modality has been developed to reconstruct high-resolution conductivity distribution images for the biological issues.In this new imaging modality, electrical impedance tomography (EIT) is combined with magnetic resonancecurrent density imaging (MR-CDI) techniques to solve the well-known ill-posedness of the image reconstruction problem in traditional EIT.Zhang [2] proposed an image reconstruction algorithm using internal current density and peripheral voltage measurement to reconstruct static conductivity images which initiated the development on the theory of MREIT.
In MREIT, currents are injected into the subject through pairs of surface electrodes.A Magnetic Resonance Imaging (MRI) scanner is used to measure the induced magnetic flux density inside the subject and the current density distribution can be calculated according to the Ampere's law.The conductivity distribution images can be reconstructed based on the relationship between the conductivity and the measured magnetic flux density combined with the current density [3].
MREIT reconstruction algorithms fall into two categories: those utilizing internal current density [4][5][6][7] and those making use of only one component of measured magnetic flux density [3,[8][9][10][11][12][13][14][15].Considering the rotation problem of the object in the MRI system, the latter has the advantages over the former without the object rotation dilemma.MREIT algorithms that are based on current density require knowledge of the magnetic flux density vector B = (Bx, By, Bz).Due to the fact that only one component of the magnetic flux density which parallels to the direction of the main magnetic field of the MRI scanner can be measured once, the rotation of the object is required, which is impractical for MRI scanner.
Recently, several MREIT algorithms have been proposed which utilize only one component of magnetic flux density, such as the harmonic Bz algorithm [8,9], the gradient Bz decomposition algorithm [11], the algebraic reconstruction algorithm [12] and an anisotropic conductivity reconstruction algorithm [16].For the head tissue conductivity, the relatively novel and concise RBF (Radial Basis Function) and RSM (Response Surface Methodology) MREIT algorithms were proposed to focus on the piece-wise homogeneous head tissue conductivity reconstruction [13,14].However, the continuous conductivity estimation of the head tissue is more useful to obtain high resolution source localization and mapping results [17].The MREIT algorithms recently applied on the head homogeneous and inhomogeneous conductivities reconstruction [18][19][20] show better results.Therefore, in this paper a MREIT approach is developed to realize the estimation of continuous conductivity distribution in the human head tissues based on the three-layer realistic FEM head model.
For the head tissues, the low skull conductivity is a dilemma for the conductivity reconstruction, since much of the current is shunted through the scalp and does not enter the brain compartment.Nevertheless, the reconstruction accuracy of MREIT is the same through the imaging space and the inner part of it will not blur, which will occur in traditional EIT.In this paper, with the feasibility of the human head being rotated twice in the MRI system, a modified continuous conductivity reconstruction MREIT algorithm which is based on two components magnetic flux density as [14] did is developed for the head tissues.With one more rotation, more information could be gained and used for the reconstruction.Even under the situation that there is only one current injection used in the simulation, better results also can be achieved.Herein, the realistic FEM head model is utilized, as the original hexahedral element meshing is unavailable in the finite element modeling.The tetrahedral element meshing is an optimum choice considering the computational efficiency.Computer simulations are performed on a three-sphere head model and the results show that the three-dimensional head model continuous conductivity reconstruction method with two magnetic flux density components could get better results than the method with only one magnetic flux density component used.
For the present approach, some assumptions should be given.First, generally speaking, it is assumed that the head geometry is known, since this data can be obtained from the MRI scanner.Second, the skull conductivity of the head tissues is fixed to the constant value to focus on the continuous conductivity distribution reconstruction of the scalp and the brain.And in the following sections, the proposed algorithm based on two components magnetic flux density is explained and an iterative conductivity update method is derived.First, the formulation and numerical solution of the forward problem and the inverse problem are presented in Section 2.Then, the performance of the proposed approach is assessed by a series of computer simulations using a three-layer uniform and non-uniform conductivity distribution realistic head models in Section 3. In Section 4, discussions and conclusion about present results are given.

METHODS
In this paper, reconstruction of conductivity image by MREIT begins with a numerically relationship between the conductivity combined with the measured magnetic flux density due to the injected currents.In the forward problem, we calculate the peripheral voltage values and the magnetic flux density for a known conductivity distribution.And the inverse problem is to reconstruct conductivity σ from the measured magnetic field distribution B = (By, Bz) and J as well as physical laws of electromagnetic.
In MREIT, MR magnitude images provide excellent structural information, and then different regions (the scalp, the skull and the brain) of the three-layer realistic head can be determined.Within each region，the conductivity values are considered continuous distribution instead of piece-wise.Such continuous conductivity of head model could be more practical for forward or inverse problem solutions and source reconstruction in functional brain imaging using magnetoencephalography (MEG) and electroencephalography (EEG).

Forward Problem
Forward problem is defined as the calculation of magnetic field generated by the internal current pattern for a given boundary injected current profile and known conductivity distribution.Forward problem formulation is used to supply simulation data to test the proposed reconstruction algorithm and it is also used in formulating the inverse problem.The relation between the conductivity and the electrical potential U(r) induced by the injected current is given by Poisson's equation together with the Neumann boundary conditions as:     0, , on positive current electrode , negative current electrode 0 .elsewhere where σ(r) is the electrical conductivity and Ω is the imaging subject space.For complex conductivity distributions, analytical solution to the forward problem expressed in Eq.1 does not exist.Therefore, a numerical method must be applied.Finite element method (FEM) is used to calculate the electrical potential and corresponding magnetic flux density distribution for a given conductivity distribution and boundary conditions.After obtaining the electric potential distribution U(r) from solving Eq.1, the electric field E and the interior current density distribution J are given as: Then the magnetic flux density can be calculated using the Biot-Savat law: where B(r) is the magnetic flux density at the measurement point, and J(r') is the current density at the source point and μ 0 is the magnetic permeability of the free space.In order to avoid the singularity occurring when r = r', B(r) is treated as a node variable and J(r') is used at the centre of each finite element in (3) (Lee et al. 2003).The comparison between analysis solution and numerical solution by FEM method is performed by Nuo Gao [13] to indicate the feasibility of solving the forward problem using the FEM.

Inverse Problem
The inverse problem is the reconstruction of σ(r) using the measured magnetic flux density combined with the calculated current density.Since static conditions are assumed, the equation below can be gained [21]: where S is defined as the natural logarithm of σ(r).
Based on the relation between current density and magnetic flux density, one can obtain: The equation is rearranged by the matrix form: With the two components magnetic flux density being used, the last two rows of the matrix equation can give: An iterative algorithm is used to solve the above nonlinear partial differential equation.The two components B z and B y used promises the equation a unique solution which is theoretically the same as the situation that two currents injection are applied.
Discretization equation is obtained by replacing derivatives with finite difference equivalents: 2 2 where Δx, Δy and Δz are the distances between adjacent pixels in the x, y and z directions, respectively.Operator  is denoted by the simple three-point difference scheme as follows: Then the equation is rearranged into the matrix form as: (10) where P is the (2 KM × M) coefficient matrix, M is the element number of the FEM head model and K denotes the injected current number.S is the M × 1 vector which is the each quadratic tetrahedral element value and Q is a 2 KM × 1 vector.Herein [ ] is the vector (KM × M) for internal current density distribution J = (J x , J y , J z ), and vector for 2  B  /μ 0 where  represents the x or y component, respectively.The matrix equation is being solved iteratively until the following stopping criterion is satisfied: where σ i+1 denotes the i + 1th iteration and σ i the ith iteration, ε is a given tolerance and  is the Euclidian norm.How to select the ε will affect the accuracy of the conductivity reconstruction algorithm.The smaller the selected ε is, the higher the accuracy of the reconstruction and the more time the convergence would be.
Then the steps of the iterative reconstruction algorithm are as follows: Step 1: According to the given target uniform and non-uniform conductivity distributions, the 'measured' two components magnetic flux density B z * and B y * for one and two currents injection are calculated using the FEM method.
Step 2: 2  B z * and 2  B y * are calculated with the three-point difference Eq.9, and the vectors Q  are obtained,  = x or y.
Step 3: An initial conductivity distribution σ m ini is assumed for m = 1, 2, … , M. This initial conductivity distribution is randomly selected to be values σ m ini ∈ (0, 15) for the first iteration.
Step 4: With σ m ini →σ m i , the current density distributions J = (J x , J y , J z ) are calculated using FEM and coefficient matrix P  is obtained for  = x or y.
Step 5: Rearrange P  and Q  as the Eq.10 and solve the combined system equation, then the solution S is found and the conductivity distribution σ m i+1 is achieved.
Step 6: Stop if the stopping criterion ( 11) is met.Othewise, set (i + 1)→i and go to Step 4 with σ m ini replaced by σ m i+1 .

RESULTS
In order to test the performance of the reconstruction algorithm using two magnetic flux density components, numerical simulations were performed on a concentric three-layer realistic human head model (consisting of the scalp, skull and brain) to estimate the continuous conductivity values σ = (σ scalp , σ brain ).

Preparation of Experiment Data
The three-layer realistic human head model has three compartments: scalp, skull and brain.Software ANSYS 10.0 was applied to the finite element modeling and meshing as well as the solving of the forward problem.A finite element meshing for the 3D three-layer realistic head model in Figure 1 with 57041 quadratic tetrahedral elements and 78893 nodes is used.And the conductivity in each element is assumed to be different values.
In order to get the current density image, a bipolar rectangular current (5mA) is injected into the human head model through a pair of opposite electrodes along the equator of the realistic head model (seeing Figure 2).A MRI scanner is used to measure the component of the induced magnetic flux density parallel to the main magnetic field of the MRI system (along the z direction), denoted as B z * .Then the head model is rotated once, the other component of the induced magnetic flux density denoted as B y * (along the y direction) is measured.The two components magnetic flux density is calculated with the given target human head conductivity distribution in the forward problem.To prove the noise tolerance of the algorithm, Gaussian White Noise (GWN) at different levels is added to the magnetic flux density to simulate the 'measured' noise-contaminated MR magnitude image B z * and B y * .The standard deviation of noise S B in 'measured' magnetic flux density is introduced by [22]: where γ = 26.75 × 10 7 radT -1 s -1 is the gyromagnetic ratio of hydrogen, T c is the duration of injection current pulse and SNR is the signal-to-noise ratio of the MR magnitude image.The SNR of the GWN is set to be infinite, 80, 60, 40, 20 and 10 with T c = 50ms, S B is obtained to be 2.9356 × 10 -9 T, 3.9142 × 10 -9 T, 5.8713 × 10 -9 T, 1.1742 × 10 -8 T, and 1.5657 × 10 -8 T, respectively.Given the skull conductivity, the continuous conductivity of the scalp and the brain is reconstructed with the proposed algorithm based on the measured two components magnetic flux density when one or two currents are injected, respectively.Two types of numerical simulations were performed on the three-layer realistic head model: uniform and non-uniform isotropic target continuous conductivity distribution reconstructions.The parameters of the head model are listed in Table 1, where σ m denotes the conductivity value of the mth element of the finite element realistic head model.
For uniform case, the conductivity distribution in each   element of the head FEM model is assumed to be identical in the scalp, skull and brain compartments, respectively.For non-uniform case, the conductivity distribution in each element is assumed to be different value in each compartment.In this study, the relative error (RE) between the estimated and the target conductivity distribution is used to quantitatively assess the performance of the MREIT reconstruction algorithm.The relative error is defined as follows: where σ * is the target conductivity distribution and σ is the estimation conductivity distribution.

Conductivity Image Reconstruction
The MREIT reconstruction algorithm based on two com-ponents magnetic flux density for the two different target conductivity distribution cases were tested with six noise levels, and some compared results are given below.With the tolerance ε = 0.01, it took seven iterations to reconstruct these continuous conductivity distributions.In each situation (uniform conductivity or non-uniform conductivity) with one current injection, simulation results are listed at Table 2.But with one current injection, meaningful reconstructions are not achieved even for noise-free simulations based on the only one component magnetic flux density algorithm.And given the two currents injection, the compared results between the algorithms based on two components and only one component magnetic flux density are given at Table 2.
For uniform simulation in Table 2, with different noise levels, all the REs are less than 1% for one current injection and two currents injection simulations, which are reduced considerably compared with the algorithm based on only one magnetic flux density.It also suggests that the RE is insensitive to the added noise.The data also show that with the noise level of SNR = 10, one current injection results are more influenced than two injection ones.The case is same to the non-uniform simulation.
For non-uniform simulation, the REs are less than 6.5% for one current injection and two currents injection simulations, which are not much less compared with the only one magnetic flux density algorithm.The reason is that the non-uniform distribution influences the performance of the algorithm.However, the non-uniform range of the real head conductivity distribution in one tissue changes not so severely, the feasibility of the algorithm on human head is confirmed.
The reconstruction conductivity images are shown in Figure 3 and  present simulation results demonstrate the excellence performance of the algorithm for conductivity reconstruction of human head tissue.
In the procedure of the computation of 2  B, threepoint difference method is utilized which will arise errors and vulnerable to measurement noise.So the means of finite difference modeling correction to the finite element modeling is a good choice to lessen the errors and improve the robustness against noise.A more sophisticated head model such as real geometry head model would be utilized to enhance the accuracy and resolution of reconstruction image at the cost of computation time.It is also found that increasing the number of injection currents does not obviously improve the solution of the reconstruction conductivity image when three or four currents are injected in the simulation.
In the present simulation, a current of 5 mA was used, which is thought to be the upper safe limit for human beings (IEC criterion).And for human head, it is a little higher.So it would be desirable to utilize a better MRI scanner, some denoising techniques and improved approaches.
In summary, we proposed MREIT approach based on two components measured magnetic flux density for noninvasive imaging of the three-dimensional continuous conductivity distribution of the head tissues.A series of computer simulations demonstrate the feasibility of the algorithm and fast convergence ability combined with improved robustness against noise.In our future studies, researches should focus on real geometry head model study and the experimental validation the algorithm on the human head phantom experiment, as well as ways to reduce the amount of the injection current down to less than 1 mA for human security consideration in clinical experiment and setting.

Figure 1 .
Figure 1.Three-layer realistic finite element head model: (a) a tetrahedral element is used to mesh the three-layer realistic head model (b) meshed model of the finite element realistic head model.

Figure 2 .
Figure 2. Current injection and electrodes locations of the realistic finite element head model.

Figure 4 ,
respectively.Figures 3(a) and (b) demonstrate the reconstruction conductivity images under the condition that without adding any noise, SNR = Infinite, SNR = 40 and SNR = 10 for uniform distribution.Figure 4 is for the non-uniform case, respectively.

Figure 4 .
Figure 4. Non-uniform conductivity reconstruction images with different noise levels: (a) corresponding to one current injection, (b) corresponding to two currents injection.(a) current injection k = 1; (b) current injection k = 2.

Table 1 .
Uniform and non-uniform conductivity head model parameters.