A method for reconstructing tomographic images of evoked neural activity with electrical impedance tomography using intracranial planar arrays

1.1.1. Methods for imaging brain activity

1.1.2. Cortical anatomy and evoked responses in the rat

1.1.3. Imaging fast neural activity with EIT

The purpose of the present study was to develop the imaging method for imaging fast neural activity in the cerebral cortex for the specific geometry of an epicortical planar electrode array. The desired specification was a spatial resolution of 0.2 mm to a depth of 1.5 mm, so as to be able to distinguish activity in any of the six cortical layers. The following sub-questions were investigated: what developments were needed to produce an optimal (a) forward, (b) inverse, and (c) image processing methods?

A related aim was to characterize the spatial accuracy of the applied methods throughout the rat brain, using electrode arrays placed over the whole cortical surface. The purpose was to assess the potential of fast neural EIT for imaging activity in deeper brain structures, such as the thalamus and hippocampus.

The scope of this work was to develop the imaging method with the required spatial resolution. Parallel studies have achieved a temporal resolution of 2 ms for imaging fast neural activity by EIT during evoked activity. Together, our results encourage the view that fast neural EIT can resolve the propagation of depolarization-related fast impedance changes in cerebral cortex and deeper in the brain with a resolution equal or greater to the dimension of a cortical column.

Simulation was performed in an anatomically realistic model of the rat brain obtained from an MRI image, with an electrode array placed on the cortex over the imaged volume (figure 1). The background conductivity was constant and isotropic throughout the grey matter (0.3 S m⁻¹), white matter (0.15 S m⁻¹), and CSF-ventricles (1.75 S m⁻¹). The anisotropy of white matter was not included. The modelled conductivity change during stimulation, the 'perturbation', was a sphere, 0.5 mm in diameter, with conductivity 1% higher than the background conductivity (0.303 S m⁻¹). It was placed in all possible locations underneath the electrode array, covering 7 × 9 × 2 mm³ in steps of 0.2 mm along each axis, which yielded 15 k simulated locations. Realistic additive Gaussian noise with zero mean and a standard deviation of 0.5 μV was added to all simulated voltage changes. Each perturbation was reconstructed using the proposed methods and its location and size compared to the actual perturbation (see section 2.5). A heuristic current injection protocol, based on the depth distinguishability (Kao et al 2003), concentrating current around the central part of the array using 30 linearly independent current injections was employed.

Zoom In Zoom Out Reset image size

2.2.1. Forward problem

The forward problem in EIT consists of the quasistatic Maxwell equations in space with the complete electrode model boundary conditions for pairs of electrodes (Somersalo et al 1992):

$$\begin{equation} \begin{array}{l} \nabla {\rm } \cdot {\boldsymbol \sigma }\nabla {\boldsymbol u} = 0, \\ {\boldsymbol \sigma }\displaystyle\frac{{\partial {\boldsymbol u}}}{{\partial {\boldsymbol n}}} = 0{\rm }\;{\rm off}\;{\rm }\bigcup\limits_{l = 1}^{\rm 2} {e_l ,} \\ {\boldsymbol u} + z_l {\boldsymbol \sigma }\displaystyle\frac{{\partial {\boldsymbol u}}}{{\partial {\boldsymbol n}}} = V_l ,\quad l = 1,2; \\ \end{array} \end{equation} \tag{ 1 }$$

where u is the electric potential, ${\boldsymbol \sigma }$ is the conductivity, n is the normal to the boundary, e_l is the lth electrode boundary, z_l is the contact impedance on lth electrode, and V_l the measured voltages. The forward problem is solved using the finite element method which employs a spatially discretized mesh of the modelled object, the brain in this case, and a numerical solution of equation (1) in matrix form for all elements. This process was repeated for all injecting electrode pair combinations with calculation of the voltages on all other electrodes. An assumption of linearity of the measured voltage changes with respect to the conductivity change in the volume was made. This allowed for the construction of a Jacobian matrix J, which mapped the changes in conductivity in each element to changes in voltages on the electrodes:

$$\begin{equation} {\rm \delta }{\bf u} = {\bf J}\delta {\boldsymbol \sigma }. \end{equation} \tag{ 2 }$$

The solution of equations (1) was performed with the finite element method using the open source finite element solver EIDORS (Adler and Lionheart 2006).

2.2.2. Meshing and mesh optimization

An optimized high-quality mesh generator based on the CGAL geometrical processing engine (Gärtner et al 2007) was used to produce tetrahedral meshes with desired quality measures directly from MRI images of the subject. Since element sensitivity decreases exponentially away from the planar array used in the present investigation, the generator was set to produce non-uniform meshes with a positive linear gradient of element size starting from the position of the array (top of the brain, z = 0) down to the opposing surface of the brain (z_p) along the depth z (figure 2). The area of the finest elements with size s_i was preserved directly underneath the array down to depth z_i. In addition, the area near the electrodes was refined independently in order to investigate optimal element size s_e for complete electrode boundary condition application. Mesh parameters (overall finest element size s_i, fraction of finest elements with respect to total depth z_i/z_p, and element size near the electrodes s_e) were optimized using mesh convergence analysis.

Zoom In Zoom Out Reset image size

The principle of mesh convergence entails refining the element size until the convergence error (difference between the solution with given element size and solution with previous element size) reaches the required precision. In the present study, the mesh was generated from an MRI image with 0.2 mm resolution. Such generated meshes with the same element size have natural variability, so the desired precision was calculated as the variance of the solution computed on the meshes produced with the same element size (ten meshes were used for each element size). Hence mesh convergence criteria were formulated with following:

$$\begin{equation} E_c = \mathop {\max }\limits_i \frac{{\overline {V_i^n } - \overline {V_i^{n - 1} } }}{{\overline {V_i^{n - 1} } }} < E_v = \mathop {\max }\limits_i \frac{{V_i^n - \overline {V_i^n } }}{{\overline {V_i^n } }}. \end{equation} \tag{ 3 }$$

Here E_c is the mesh convergence error, E_v is the variance of the solution obtained on nth refinement step, and $V_i^n$ is the absolute voltage obtained on ith electrodes with nth mesh refinement step for a particular parameter.

Obtained errors were assessed from plots with respect to the optimized parameters. All considered meshes were quality-checked with the Joe–Liu quality measure (Liu and Joe 1994):

$$\begin{equation} q = 12\frac{{( {3V})^{\frac{2}{3}} }}{{\mathop \sum _{0 \le i < j \le 3} l_{ij}^2 }}, \end{equation} \tag{ 4 }$$

where q is quality, V is the element volume, and l is the edge length between ith and jth element vertices. All meshes were insured to have average q > 0.9 and min  q > 0.1.

The inverse solution in EIT for time difference mode comprises inversion of J in order to obtain a linear map between electrode voltages and conductivity change in volume (J⁻¹). As J is ill-conditioned, regularization is required to perform this operation. Zeroth-order Tikhonov regularization minimizes the l2-norm magnitude of the changes:

$$\begin{equation} \widehat{{\rm \delta }{\boldsymbol \sigma }} = {\rm arg}\;{\rm min}_{{\rm \delta \sigma }} \;\Vert{\rm \delta }{\bf v} - {\bf J}{\rm \delta }{\boldsymbol \sigma }\Vert^2 + {\lambda} \Vert\delta {\boldsymbol \sigma }\Vert^2 , \end{equation} \tag{ 5 }$$

$$\begin{equation} {\bf J}_\lambda ^{ - 1} = ({{\bf J}^T {\bf J} + \lambda {\bf I}})^{ - 1} {\bf J}^T . \end{equation} \tag{ 6 }$$

$\widehat{{\rm \delta }{\boldsymbol \sigma }}$ is the estimated conductivity change, δv is the boundary voltage change, and is a constant hyperparameter which was set using cross-validation.

The optimal value of λ (the value that minimized the cross-validated error) was established in the range from 10⁻²⁰ to 1 in 3000 logarithmically spaced steps, and the cross-validated error was estimated using the leave-one-out cross-validation method (Picard and Cook 1984), which used a single observation from the original sample as the validation data and the remaining observations as the training data.

This was repeated such that each observation in the sample was used once as the validation data. This was the same as a K-fold cross-validation, with K being equal to the number of observations in the original sampling. After the optimal value for λ, was obtained, the inverted Jacobian matrix was computed using equation (6).

In order to improve the conditioning of the inverse problem, and reduce the computational time, a hexahedral mesh was produced as a subdomain of the main sensitivity matrix as follows (figure 3): (a) a region of interest volume under the most sensitive area of the array, 2.5 mm deep, was extracted from the tetrahedral mesh, (b) depth slicing was performed and tetrahedra were labelled according their position with respect to the outer surface in 0.1 mm steps, (c) each layer was unwrapped to a planar cuboid and tetrahedra were labelled with their x and y positions with a step of 0.1 mm. The resulting mesh represented the unwrapped region of the cortex 2.5 mm underneath the array, where all elements were cubes with 0.1 mm side length, and the total number of elements was 60 000. The Jacobian matrix was then mapped onto the new hexahedral mesh:

$$\begin{equation} J_{{\rm hex}}^i = \mathop \sum _{k\epsilon P_i } J_{{\rm tetra}}^k , \end{equation} \tag{ 7 }$$

where P_i is the subset of the elements from the tetrahedral forward mesh which belong to the ith hexahedron of the inverse mesh.

Zoom In Zoom Out Reset image size

In order to validate the coarsened hexahedral mesh and sensitivity matrix, ten reconstructions were performed using the original tetrahedral sensitivity matrix and compared to those obtained with the hexahedral sensitivity matrix. The spatial difference between solutions was less than 0.1 mm, which is the size of one hexahedron.

EIT with planar arrays produces distortion, especially with respect to the distance from the electrode plane, due to a poor general sensitivity field for planar arrays (Kao et al 2003). Several methods exist for performing regularization which account for the extreme inhomogeneity of sensitivity in this case (Kaipio et al 1999); however, all of these require significant additional priors regarding the assumed perturbation and therefore are unsuitable for fast neural EIT as the shape and properties of the impedance change are uncertain in advance.

Here we present and validate a novel technique which is based on noise-based, t-score masking of the conductivity changes obtained with Tikhonov regularization. The method comprises use of the inverted Jacobian matrix $J_\lambda ^{ - 1}$ in order to compute the standard deviation of the estimated conductivity change in each hexahedron due to random Gaussian noise in the voltage measurements. k = 10 000 samples of the computer-generated Gaussian normal distribution were taken to compute deviation correctly:

$$\begin{equation} \begin{array}{l} \delta {\boldsymbol \sigma }_n = {\bf J}_\lambda ^{ - 1} \delta {\bf v}_n , \\ \delta {\bf v}_n \sim N ({0,0.5 \cdot 10^{ - 6} } ). \\ \end{array} \end{equation} \tag{ 8 }$$

$\delta {\boldsymbol \sigma }_n$ is the conductivity change matrix of size m by k (m is the number of elements) caused by noise in the voltage change δv_n on the electrodes. The t-score for the actual conductivity change in ith element can then be computed:

$$\begin{equation} t_i = \frac{\delta \sigma ^i}{\big\langle {\delta \sigma _n^i - \overline {\delta \sigma _n^i }} \big\rangle }. \end{equation} \tag{ 9 }$$

Then the image was masked by the t-score (all voxels with t-score less than 3 are considered non-significant: .δσ|_{t < 3} = 0), and spatially smoothed with a Gaussian filter (with 0.3 mm kernel).

The accuracy and resolution of image reconstruction were determined by computing the centre of mass and standard deviation of the spatial conductivity distribution for each reconstructed perturbation:

$$\begin{equation} \begin{array}{l} {\bf r}_c = \displaystyle\frac{{\mathop \sum _{i = 1}^m ( {\delta \sigma ^i \cdot r_i })}}{{\mathop \sum _{i = 1}^m \delta \sigma ^i }} \\ \varepsilon _r = \displaystyle\frac{{\langle {\delta {\boldsymbol \sigma } \cdot {\bf r} - \overline {\delta {\boldsymbol \sigma } \cdot {\bf r}} } \rangle }}{{\mathop \sum _{i = 1}^m \delta \sigma ^i }}. \end{array} \end{equation} \tag{ 10 }$$

r_c is the centre of mass, ε_r is the standard deviation of the conductivity distribution in three dimensions, which represents the radius of the fussy reconstructed object, r_i is the coordinates of ith element, and δσⁱ is the conductivity value in ith element. Localization error (difference between the true and reconstructed position of the perturbation: $\left| {{\bf r}_c - \widehat{{\bf r}_c }} \right|$ was computed and displayed, as well as shape error $\left| {\varepsilon _r - \widehat{\varepsilon _r }} \right|$, where $\widehat{{\bf r}_c }$ is the true position of the centre of perturbation, and $\widehat{\varepsilon _r }$ = 0.25 is the radius of perturbation.

After evaluation, the developed methods were employed to study the feasibility of imaging deep brain neural activity with a 3D arrangement of electrodes. A 15 M element optimized forward mesh, 100 k inverse hexahedral mesh, and 256 electrodes placed around the brain, with 64 of them being injecting, were employed (figure 4). A 0.5 mm spherical perturbation with a 10% conductivity change, simulating neuronal activity, was placed in the hippocampal area CA1 and the improved image reconstruction procedures described above were utilized.

Zoom In Zoom Out Reset image size

The mesh convergence analysis indicated readily identifiable optimal parameters, which satisfied the convergence criteria. s_i and $z_i /z\_p$ were set at the minimum value where the convergence error reached mesh variability (figures 5(a) and (b)). On the other hand, the error did not vary across the search space for parameter s_e (figure 5(c)), and therefore s_e was set at the minimum value. The search space for s_e commenced at s_i/2 for optimal s_i, because electrodes had to be refined for the mesh where convergence for overall problem has been reached (table 1).

Zoom In Zoom Out Reset image size

Total 15 k of simulations has been performed for all possible perturbation locations, of which the single example, in which the perturbation was placed in the middle of the array at 0.8 mm depth, is shown on figure 6. It demonstrates the improvement of the noise-based, t-score masking in comparison to raw conductivity change.

Zoom In Zoom Out Reset image size

The localization error was assessed at depths of 0.4 and 0.8 mm throughout the whole area of the array. At a depth of 0.4 mm, corrected and uncorrected reconstruction produced similar localization errors throughout the volume under the array (figure 7). Uncorrected or corrected reconstruction produced an average error of 0.16 or 0.11 mm respectively.

Zoom In Zoom Out Reset image size

At 0.8 mm, the region where depth accuracy for imaging evoked sensory neural activity is arguably most important, as this corresponded to layer IV which receives the initial cortical input (Armstrong-James et al 1991), uncorrected reconstruction produced a mean error of 0.5 mm. Corrected reconstruction matched the expected depth to within 0.08 mm (figure 8).

Zoom In Zoom Out Reset image size

Further assessment of the localization error at depths 0.2–2 mm throughout the volume under the array indicated that under the centre of the array (about 4 × 5 mm²), there was a uniform accuracy of 0.4 mm for uncorrected, and 0.1 mm for corrected reconstruction. The average shape error $\left| {\varepsilon _r - \widehat{\varepsilon _r }} \right|$ was also homogenous throughout this area, being 0.12 and 0.01 mm for uncorrected and corrected reconstruction respectively.

The accuracy of localization for a spherical perturbation 0.5 mm ($\hat \varepsilon _r$ = 0.25 mm) in diameter with respect to depth in the centre of the array to 2 mm deep was for uncorrected reconstruction was 0.24 mm (no noise) and 0.4 mm (with noise). With noise-based correction it was 0.05 mm for both cases (figure 9).

Zoom In Zoom Out Reset image size

The average shape error $\left| {\varepsilon _r - \widehat{\varepsilon _r }} \right|$ was 0.05 mm for uncorrected without noise, 0.06 mm for uncorrected with noise, and 0.02 mm for corrected (both with and without noise) image reconstruction.

Overall, with noise-based correction, the total region of acceptable accuracy and resolution of <0.1 mm was an area 4 mm × 5 mm around the centre of the array down to 1.4 mm deep.

The 0.5 mm diameter perturbation reconstruction for the deep brain neural activity produced a localization error $\left| {r_c - \widehat{r_c }} \right|$ = 0.5 mm with the shape error being $|\varepsilon _r - \hat \varepsilon _r |$ = 0.35 mm (figure 10).

Zoom In Zoom Out Reset image size

The results show that the proposed chain of methods and image reconstruction methodology yields an accuracy within the desired parameters of the experimental setup (SNR = 5, accuracy <0.2 mm). The forward and inverse problems were optimized for the application of EIT to imaging neural activity with intracranial planar arrays. Mesh optimization resulted in the identification of optimal parameters for this application. A novel noise-based image post-processing technique was proposed and tested, resulting in significant improvements in accuracy of activity localization, in comparison to a standard reconstruction technique based solely on Tikhonov regularization.

Overall, the volume of acceptable accuracy with localization error <0.1 mm, was a region 4 mm × 5 mm around the centre of the array, 1.4 mm deep. This indicates that the developed methods may be expected to be reliably applied for imaging neural activity with planar arrays.

In addition to this, feasibility simulations have been performed which showed that deep brain hippocampal or thalamic activity imaging with the existing experimental parameters is possible with an accuracy of 0.5 mm.

During the study, the current injection protocol was selected heuristically; however, superior results may be obtained with a more optimized current injection pattern. This is a substantial project which is work in progress and is outside of the scope of this paper. In addition, the feasibility of deep brain neural activity has been performed only for the example of the hippocampus; further feasibility studies regarding deep cerebral structures should be guided by experimental hypotheses to be addressed.

The developed methods for planar array imaging are currently being applied to experimental data collected in vivo. In these experiments, simultaneous recording is being performed using intrinsic optical imaging and local field potentials as independent yardsticks. These results display the potential of EIT for deep brain imaging and imaging with humans. It is planned to assess the resolution throughout the entire brain further and more methodically using 256 electrodes wrapped around the brain. In order to achieve this aim, new current injection protocols for such high density recordings are also under development.