Accesso libero

Finite difference model of a four-electrode conductivity measurement system

INFORMAZIONI SU QUESTO ARTICOLO

Cita

Introduction

Four electrode measurement (4EM) systems have been used to measure conductivity in tissue where electrode polarization would cause unacceptable inaccuracy in measurement results and where directional measurement of conductivity is required

Less polarization than a two-electrode measurement system where the current is applied by the same electrodes that measures the potential and directional measurement as opposed to a coaxial measurement system that measures bulk conductivity but gives no information about direction.

(1, 2, 3). The system described here is intended for use in testing conductivity in brain tissues in vitro.

The response of the 4EM has been well characterized in homogeneous, isotropic media as (1, 3):

V=I4πσa$$\text V=\frac{I}{4\pi \sigma a}$$

Equation (1) assumes the four electrodes are linearly arranged with two high impedance potential measurement electrodes (V1,V2) flanked by two current injection electrodes (I1,I2). It also assumes the electrodes have uniform spacing (a) and are surrounded by a medium with conductivity σ. The injected current is from a point source (I1) and is collected by a point sink (I2) (Figure 1).

Figure 1

Electrode configuration of a 4EM system.

In addition to the response function (Equation 1), Robillard and Poussart (4) also characterized the field of view of the 4EM system as approximately 3a in homogeneous isotropic media. In other words, conductivity differences in media outside a radius of 3a from any of the electrodes should have little effect on the response.

Steendijk (5) found solutions for 4EM systems for several conductivity environments. Starting with the infinite space of isotropic material solution (1) (Equation 1), he found a solution for an infinite half space with electrodes placed on the air/media boundary:

V=I2πσa$$V=\frac{I}{2\pi \sigma a}$$

He also found equations for an infinite half space of anisotropic material, longitudinal (long) and transverse (trans) to the fiber direction:

Vsi,long=I2πaσtransVsi,trans=I2πa1σtrans1σlong$$\begin{align}& {{V}_{si,long}}=\frac{I}{2\pi a{{\sigma }_{trans}}} \\ & {{V}_{si,trans}}=\frac{I}{2\pi a}\sqrt{\frac{1}{{{\sigma }_{trans}}}\frac{1}{{{\sigma }_{long}}}} \\ \end{align}$$

Steendijk’s solutions only hold for specific geometries. A real world measurement in anisotropic or mixed tissue would necessarily require a much more complex solution, if such a solution is even possible.

This paper will demonstrate a finite difference numerical solution based on a three dimensional matrix of conductivity tensors to support any combination of included regions, limited only by the resolution and size of the grid. The tradeoff between either a finer grid size or increased spatial resolution is processing time.

Materials and methods

Dey and Morrison (6) developed a finite difference solution for regions of differing isotropic conductivity in geological media. This paper will extend their solution to allow differing anisotropic regions. The solution begins with a few basic equations:

J=σE$$\text J=\text{ }\!\!\sigma\!\!\text{ E}$$E=ϕ$$\text{E}=-\nabla\text{ϕ}$$J=ρtδ(xyz)$$\nabla⋅\text{J=}\frac{\partial\text{ρ}}{\partial\text{t}}\delta{(\text{xyz})}$$

Equation (4) is the differential form of Ohm’s law that relates current density (J) in a medium of conductivity σ to the electric field (E) in the medium. Equation (5) holds because the field is irrotational so E is the gradient of the potential (ϕ), and equation (6) relates the time derivative of a point charge (ρ) at location xyz to the divergence of J.

Combining these equations yields:

[σxyzϕxyz]=ρtδ(xyz)$$-\nabla \cdot \left[ {{\sigma }_{xyz}}\nabla {{\phi }_{xyz}} \right]=\frac{\partial \rho }{\partial t}\delta \left( xyz \right)$$

where the quantities σxyzandϕxyz${{\sigma }_{xyz}}\operatorname{and}{{\phi }_{xyz}}$ are conductivity and potential at point xyz.

The derivative of charge with respect to time is equivalent to the current within a volume element: ρt=IΔV.$\frac{\partial \rho }{\partial \text{t}}=\frac{\text{I}}{\Delta \text{V}}.$ Integrating both sides over the volume element (dV):

ΔV[σxyzϕxyz]dV=ΔVIΔVδ(xyz)dV$$-\iiint_{\Delta V}{\nabla \cdot \left[ {{\sigma }_{xyz}}\nabla {{\phi }_{xyz}} \right]}\,dV=\iiint_{\Delta V}{\frac{I}{\Delta V}}\,\delta \left( xyz \right)dV$$

The left hand side of Equation (8) can be converted to a surface integral using Green’s theorem:

ΔV[σxyzϕxyz]dV=SσϕndA=sσϕndA$$\begin{align}& -\iiint_{\Delta V}{\nabla \cdot \left[ {{\sigma }_{\text{xyz}}}\nabla {{\phi }_{\text{xyz}}} \right]}\,\text{dV} \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\text{=}-\iint_{\text S}{\sigma \,\nabla \phi }\cdot \text{n}\,\text{dA} \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\text{=}-\iint_{\text s}{\sigma \frac{\partial \phi }{\partial \text{n}}\text{dA}} \\ \end{align}$$

where dA is the surface area element.

And the right hand side of Equation (8) is simply equivalent to the current at the point xyz, Ixyz. The resulting equation can be discretized:

sσϕndA=Ixyz$$-\iint_{s}{\sigma \frac{\partial \phi }{\partial \text{n}}}\text{dA=}{{\text{I}}_{\text{xyz}}}$$

To simplify the equations, we utilize a uniform spatial grid throughout the matrix where the voxels have the same dimension in each direction, Δx = Δy = Δz = h. The left-hand side of Equation (10) is integrated by breaking into separate equations for each face and using the finite difference approximation for the normal derivative.

Dey and Morrison estimated the conductivity at each voxel face by averaging the conductivities of four adjacent voxels. In order to account for anisotropic conductivities, we used conductivity of the adjacent cell in a direction normal to the face being integrated (σnorm). This should introduce minimal error as long as the change in conductivity is gradual local to the voxel. Care should be taken in interpreting the potentials immediately adjacent to a boundary between two conductivities.

The discretized solutions for each face are summed to comprise the surface integral:

faceσnorm(ϕadjϕselfh)h2$$-\sum\nolimits_{\text{face}}{{{\sigma }_{\text{norm}}}}\left( \frac{{{\phi }_{\text{adj}}}-{{\phi }_{\text{self}}}}{\text{h}} \right){{\text h}^{2}}$$

where:

sdA=h2$$\iint_{\text s}{\text{dA=}{{\text{h}}^{2}}}$$

and the values of the conductivity (σadj) and potential (ϕadj) of the adjacent voxel are defined as in Figure 2. Coupling

Figure 2

Orientation, conductivity, and potential of each face of a cube.

constants as defined by Dey and Morrison (6) are developed for each face and for the object voxel, refer to Figure (2) for definition of conductivity for each face:

C1=σ1*hC2=σ2*hC3=σ3*hC4=σ4*hC5=σ5*hC6=σ6*hCself=(C1+C2+C3+C4+C5+C6)$$\begin{align}& \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{C}_{1}}={{\sigma }_{1}}*h \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{C}_{2}}={{\sigma }_{2}}*h \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{C}_{3}}={{\sigma }_{3}}*h \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{C}_{4}}={{\sigma }_{4}}*h \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{C}_{5}}={{\sigma }_{5}}*h \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{C}_{6}}={{\sigma }_{6}}*h \\ & {{C}_{self}}=-\left( {{C}_{1}}+{{C}_{2}}+{{C}_{3}}+{{C}_{4}}+{{C}_{5}}+{{C}_{6}} \right) \\ \end{align}$$

All the analysis to this point applies to voxels in the interior of the volume of a conductive region. Because the size of the region is necessarily limited, some attention needs to be paid to the behavior at the limits of the computational volume. Voxels at the boundary of the volume require a boundary condition to allow accurate prediction of potentials across the boundary.

The simplest boundary condition is to assume the potential at the boundaries is a constant (usually zero), as if the boundaries are at infinite distance from the probe (Dirichlet boundary condition). Implementation of the Dirichlet boundary condition resulted in an overprediction of the measured potentials by 25% for anisotropic cases, presumably due to interactions of fields with finite boundaries. Although Robillard and Poussart (4) found a field of view of 3a in an isotropic medium, an anisotropic medium may alter the field of view by extending the range of the currents in the low conductivity direction.

A mixed boundary condition:

αϕ+βϕn=f(x,y,z)$$\alpha \phi +\beta \frac{\partial \phi }{\partial \text{n}}=\text{f}\left( \text{x,}\,\text{y,}\,\text{z} \right)$$

may yield a more accurate prediction. The field at a distance (r) from the point current source behaves as ϕ1r,$\phi \propto \frac{1}{\text{r}},$ so the boundary condition becomes:

ϕn1r2e^rn^$$\frac{\partial\phi}{\partial\text{n}}\propto-\frac1{\text{r}^2}{\widehat{\text{e}}}_\text{r}⋅\widehat{\text{n}}$$

where e^r${\widehat{\text{e}}}_\text{r}$ is the unit vector in the direction from the point current source to the boundary location and n^$\widehat{\text{n}}$ is the unit vector normal to the face at the boundary. e^r^=rcos(θ)$\widehat{\text{e}}⋅\widehat{\text{r}}\text{ }\text{=}\text{ }\text{r}\text{ }\text{cos}{(\theta)}$ so the boundary condition is:

αϕ+βϕcos(θ)r=f(x, y, z)$$\alpha \phi +\beta \frac{\phi \cos \left( \theta \right)}{\text{r}}=\text{f}\left( \text{x, y, z} \right)$$

where cos (θ) is the cosine of the angle between the vector from the point current source to the boundary voxel and the normal direction of the face on the boundary. And r is the distance from the source to the voxel. This definition greatly complicates the execution of the model because creation of a new finite difference matrix is necessary for each source location. To simplify model creation, we can assume the source is always at the middle of the geometric model volume. This causes minimal error when the source is near the center of the volume (6).

As a result of the mixed boundary condition, the potential for faces on the boundary is assumed to be zero (for instance, C1 = 0) and Cself (Equation 16) gets an extra term. (Note: Face 6 is opposite Face 1).

The subscripts ctr and bnd refer to the center of the geometric model volume and the voxel on the boundary, respectively. Because the Cxbnd terms are zero except for the face on the boundary, the inclusion of the extra term only slightly increases the complexity of the model.

Cself=(C1+C2+C3+C4+C5+C6)(C1bnd+C2bnd+C3bnd)+C4bnd+C5bnd+C6bnd)$$\begin{align}& {{\text C}_{\text{self}}}=-\left( {{\text C}_{1}}+{{\text C}_{2}}+{{\text C}_{3}}+{{\text C}_{4}}+{{\text C}_{5}}+{{\text C}_{6}} \right) \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-\left( {{\text C}_{1\text{bnd}}}+{{\text C}_{2\text{bnd}}}+{{\text C}_{3\text{bnd}}} \right) \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+{{\text C}_{4\text{bnd}}}+{{\text C}_{5\text{bnd}}}+{{\text C}_{6\text{bnd}}}) \\ \end{align}$$

Where

C1bnd=C6*|zctrzbnd|(xctrxbnd)2+(yctrybnd)2+(zctrzbr$${{\text C}_{1\text{bnd}}}={{\text C}_{6}}*\frac{\left| {{\text{z}}_{\text{ctr}}}-{{\text{z}}_{\text{bnd}}} \right|}{\sqrt{{{\left( {{\text{x}}_{\text{ctr}}}-{{\text{x}}_{\text{bnd}}} \right)}^{2}}+{{\left( {{\text{y}}_{\text{ctr}}}-{{\text{y}}_{\text{bnd}}} \right)}^{2}}+\text{(} {{\text{z}}_{\text{ctr}}}-{{\text{z}}_{\text{br}}} }}$$

So the discretization of Equation (10) is:

C1ϕ1+C2ϕ2+C3ϕ3+C4ϕ4+C5ϕ5+C6ϕ6+Cselfϕ(i,j,k)=Ixyz$$\begin{align}& {{\text{C}}_{1}}{{\phi }_{1}}+{{\text{C}}_{2}}{{\phi }_{2}}+{{\text{C}}_{3}}{{\phi }_{3}}+{{\text{C}}_{4}}{{\phi }_{4}}+{{\text{C}}_{5}}{{\phi }_{5}}+{{\text{C}}_{6}}{{\phi }_{6}} \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+{{\text{C}}_{\text{self}}}{{\phi }_{\left( \text i,\text j,\text k \right)}}={{\text{I}}_{\text{xyz}}} \\ \end{align}$$

This equation can be written in matrix form as:

[C][ϕ]=[S]$$\left[ \text{C} \right]\left[ \phi \right]=\left[ \text{S} \right]$$

Where [C] is the matrix of coupling coefficients, [ϕ] is the vector of unknown potentials and [S] is the source vector which is zero everywhere except at the location of the point current sources.

In order to construct the finite difference matrix, the 3-dimensional physical matrix must be mapped to a 2-dimensional matrix where each column and each row represent one of the unknown potentials. For instance, a 10×10×10 cell region (1000 cells or voxels) would map to a 2-dimensional matrix with 1000×1000 elements representing 1000 equations with 1000 unknowns that must be solved simultaneously to determine the potential in each cell. Each row of the matrix is the equation for one unknown potential and each column element represents one unknown of the equation.

Figure 3 shows the top right-hand corner of the finite difference matrix representing the system of n equations and n unknowns where the Cs are the coupling coefficients with the adjacent cell. Because the potential in a cell only depends on the six adjacent potentials and its own potential, the finite difference matrix is sparse and can be manipulated using iterative techniques in Matlab (MathWorks Inc., Natick, MA) to solve for the potential in each cell. The full Matlab code is http://dx.doi.org/10.5617/jeb.2641available (see bottom of page).

Figure 3

Top right had corner of finite difference matrix. C is the multiplier of the unknown potential at each adjacent face.

Several aspects of the Matlab code deserve mention. The input to the Matlab code is a geometric model consisting of a 3 dimensional array of diagonalized conductivity tensors. The implementation assumes the grid spacing to be the same in all directions (Δx = Δy = Δz = h). The largest practical size of the grid is 100×100×100 voxels on a quad-core MacBook Pro. With a grid spacing of 0.1 mm, this provides a space of 10 mm × 10 mm × 10 mm that is large enough to contain the 4EM system electrodes with a boundary of approximately 3a. The grid is built by a custom Matlab function (GM_builder.m) using a mathematical description of the medium. The medium can have isotropic and anisotropic regions and could be imported from a diffusion tensor image to calculate the coupling coefficients if the diffusion tensor is assumed to be proportional to the conductivity tensor.

The most computationally expensive part of the code is construction of the finite difference matrix (Matlab function SetUpFDMatrix_impbound.m). The matrix is built one line at a time using logical indexing and the command sub2ind that converts the three dimensional coordinates (i,j,k) to a serial number (1…N). Matlab references matrix indices in row, column, slice order so a large flat matrix can be built fairly efficiently. On a quad-core MacBook Pro this step can take from 20 minutes to one hour to complete. Fortunately, this step only needs to be accomplished once for each geometry, as it is not dependent on the location of the sources or measurement points.

The final step in predicting the response of the 4EM system is to specify the location of the point current sources (one positive and one negative) and the location of the measurement points (as shown in Figure 1). The Matlab functions SinglePoint.m and MultiPoints_angulareror.m compute the potential for a single probe location or an array of probe locations within the medium. The functions take the center of the probe location and the direction of the probe axis and compute the location of each current point source and measurement location to maintain an appropriate geometry. In a physical system there could be error in the positioning of the probes, so the function MultiPoints_angularerror.m allows specification of three angular errors, θx, θy, and θz. The function can show the impact of having the probes misaligned with the interface directions on the response of the 4EM.

The solution of the finite difference equation cannot be accomplished using the straightforward Matlab backslash operator ([ϕ] = [S]\[C]) due to memory limitations of the computer. Matlab offers several methods of iterative solutions for sparse matrices that make the solution possible. The General Method of Residuals (gmres) (7) is a robust method that is successful at solving the system regardless of the content of the geometric model. However, it can take up to an hour to solve the equation. The Symmetric LQ (symmLQ) (8) method is much faster for some geometries, but fails to converge if the geometries are asymmetric such as an infinite half-space. The method we selected is the Quasi-minimal residual (qmr) (9) method because it is very fast (20-40 seconds) and converges for all geometries attempted. Note that all the methods give very similar results for conditions where they converge.

Results

The finite difference model can be compared to closed form solutions like Steendijk (5). Equation (1) provides a solution for an infinite space with isotropic conductivity. Equation (4) provides solutions for electrodes on the surface of an anisotropic halfspace. Equation (4) can be modified to predict the response in an infinite anisotropic space by multiplying by a factor of 1/2. As long as the distance between the outer electrodes and the boundary are more than approximately three times the electrode spacing (3×a) the boundary will have minimal impact on the measurement.

Finite difference models are validated based on convergence to the known exact value for the solution with decreasing voxel size (h). Figure (4) shows a graph of predicted potential in an isotropic infinite space with decreasing ‘h’ as compared to the exact solution from Equation (1). Figures (5) and (6) show graphs of predicted potential versus ‘h’ for longitudinal and transverse probe orientations in an infinite anisotropic space (σ1 = 0.5, σt = 0.1). These plots show fairly close convergence to the analytical solution with decreasing voxel size (Approximately 8% difference for the anisotropic transverse cases, approximately 2% for the other cases).

Figure 4

Convergence plot for an isotropic infinite space.

Figure 5

Convergence plot for an anisotropic infinite space, transverse direction

Figure 6

Convergence diagram for an anisotropic infinite space, longitudinal direction

Models can be further validated by comparing them to theory over a range of values. Figures (7), (8), and (9) compare model predictions at increasing depth with theoretical predictions (Equation (4)) at the surface of an infinite half space and at depth (essentially an infinite space). Note that the finite difference model slightly underpredicts (by 8%) the value at the interface. This may be inaccuracy caused by the simplistic formulation (kernel) of the finite difference model. However, at depth, the model predicts within approximately 5% of the theoretical value.

Figure 7

Potential on the surface and at depth in an isotropic infinite halfspace

Figure 8

Potential on the surface and at depth in an anisotropic infinite halfspace, transverse direction

Figure 9

Potential on the surface and at depth in an anisotropic infinite halfspace, longitudinal direction

Accepting that the finite difference model is a fairly good predictor of the response of the 4EM, a number of configurations can be examined. In use, the 4EM is advanced down through tissue through a layer that is thought to be anisotropic. In order to understand the response as the probes traverse the anisotropic layer, a geometry model was built with an anisotropic layer, 2mm in thickness, flanked by infinite isotropic layers. Additionally, an isotropic infinite space, and an isotropic infinite halfspace will be analyzed (Figure (10)).

Figure 10

Geometries for analysis of 4EM system response. Top figure: infinite space. Central figure: infinite halfspace with probes on surface. Bottom figure: layer inclusion in infinite space.

Figure (11) shows the predicted response of the 4EM system as the probe is traversed from 3 to 7mm in depth (the z direction) in the center (x and y directions) of the space. The isotropic media in these analyses have an electric conductivity of σ = 0.2 mS/cm and the anisotropic media have conductivities longitudinal σ1 = 0.5 mS/cm, and transverse σt = 0.1 mS/cm.

Figure 11

Predicted response of the 4EM system at depths from 3 to 7 mm

The response of the 4EM system probes in the anisotropic layer is reduced by interactions with the surrounding medium as compared to the infinite case. In the transverse direction, the response is reduced by 30.2% and in the longitudinal direction the response is reduced by 20.2%.

Another use of the finite difference model is to explore the impact of probe positioning error on the 4EM response. Figure (12) shows the probe configuration with a 5-degree positioning error. Note that the probes enter the layer at different probe ‘z’ depths. The effect of this can clearly be seen in the graphs in Figure (13). The error causes the response of the 4EM system to be reduced by 4% in the longitudinal direction and 3.3% in the transverse direction.

Figure 12

Probe with 5 degree positioning error

Figure 13

Effect of 5 degree positioning error on 4EM system response

Discussion

The FD model converges to the closed form result and agrees with the closed form solutions across a range of values within 10% in all cases. Furthermore, the predicted response of the 4EM system is reasonable as the probes traverse down through an anisotropic layer. Future work on this effort could include measurement of anisotropic media of known conductivity to further validate the model.

The FD model will be extremely useful to those who need to directly measure conductivity with depth in tissue if the nature of the tissue structure is known well enough to form a geometric model. It is also useful to analyze a range of tissue and geometric properties (e.g. layer thickness, anisotropy ratio) for comparison to measured data and to explain 4EM system response in anisotropic heterogeneous media. The demonstrated errors caused by a 5-degree error in probe positioning as shown in Figure 13, help provide confidence in the usefulness of the 4EM system in a practical measurement situation.

Matlab code:

https://github.com/joeb-files/2016_Montgomery