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

A method is presented for reconstructing images of fast neural evoked activity in rat cerebral cortex recorded with electrical impedance tomography (EIT) and a 6 × 5 mm2 epicortical planar 30 electrode array. A finite element model of the rat brain and inverse solution with Tikhonov regularization were optimized in order to improve spatial resolution and accuracy. The optimized FEM mesh had 7 M tetrahedral elements, with finer resolution (0.05 mm) near the electrodes. A novel noise-based image processing technique based on t-test significance improved depth localization accuracy from 0.5 to 0.1 mm. With the improvements, a simulated perturbation 0.5 mm in diameter could be localized in a region 4 × 5 mm2 under the centre of the array to a depth of 1.4 mm, thus covering all six layers of the cerebral cortex with an accuracy of <0.1 mm. Simulated deep brain hippocampal or thalamic activity could be localized with an accuracy of 0.5 mm with a 256 electrode array covering the brain. Parallel studies have achieved a temporal resolution of 2 ms for imaging fast neural activity by EIT during evoked activity; this encourages the view that fast neural EIT can now 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.

A method is presented for reconstructing images of fast neural evoked activity in rat cerebral cortex recorded with electrical impedance tomography (EIT) and a 6 × 5 mm 2 epicortical planar 30 electrode array. A finite element model of the rat brain and inverse solution with Tikhonov regularization were optimized in order to improve spatial resolution and accuracy. The optimized FEM mesh had 7 M tetrahedral elements, with finer resolution (0.05 mm) near the electrodes. A novel noise-based image processing technique based on t-test significance improved depth localization accuracy from 0.5 to 0.1 mm. With the improvements, a simulated perturbation 0.5 mm in diameter could be localized in a region 4 × 5 mm 2 under the centre of the array to a depth of 1.4 mm, thus covering all six layers of the cerebral cortex with an accuracy of <0.1 mm. Simulated deep brain hippocampal or thalamic activity could be localized with an accuracy of 0.5 mm with a 256 electrode array covering the brain. Parallel studies have achieved a temporal resolution of 2 ms for imaging fast neural activity by EIT during evoked activity; this encourages the view that fast neural EIT can now resolve the propagation of depolarization-related fast impedance Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Over the past two decades, there have been huge advances in the ability to image brain function. Techniques such as functional MRI and positron emission tomography have enabled imaging of changes in blood flow and metabolic changes, but these are only over seconds as they are secondary to brain activity. The 'holy grail' would be to image the primary substrate of brain function which is fast electrical activity over milliseconds. Electrical impedance tomography (EIT) is an emerging imaging technique which enables tomographic images of the electrical properties of a subject to be produced with rings of electroencephalography (EEG)-type electrodes. It has the unique potential to produce tomographic images of the small impedance changes which occur in the brain as ion channels open during neuronal depolarization. In 2011, our group published a method for measuring impedance changes due to fast neural activity in the somatosensory cortex using a planar epicortical electrode array during repeated sensory stimulation (Oh et al 2011). However, this method would only allow for imaging with temporal and spatial resolution of 8 ms and >500 μm. Since cortical activity has components up to 600 Hz, and it is desirable to be able to discriminate vertical propagation within cortical layers and columns which are about 250 μm across, this work was undertaken to improve the spatial and temporal resolution.
1.1. Background 1.1.1. Methods for imaging brain activity. Investigation of information processing in the cerebral cortex and wider brain is of great topical interest. The ideal method would enable imaging of all neurons in the brain with a temporal resolution of 1 ms or so. Currently no method accomplishes this. Established methods such as current source density analysis using linear microelectrodes (Vertes and Stackman 2011) or two photon imaging of intracellular calcium (Denk et al 1990) provide detailed information for a column or cube of cortex up to about 1 mm across, but do not provide information about lateral spread across the cortex. In principle, this information may be produced using microelectrode arrays such as the Utah (Maynard et al 1997) or Michigan (Langhals and Kipke 2009) arrays, but these are highly invasive penetrating arrays which may interfere with normal physiology. Electrical brain activity may be recorded with epicortical or scalp electrodes coupled to inverse source modelling of the EEG or magnetoencephalography (MEG), but these do not have a unique inverse solution and so require simplifying assumptions of uncertain validity especially for deep or widely distributed sources (Hämäläinen et al 1993, Baillet et al 2001. Other methods also have limitations: visible light optical mapping has limited penetration into the brain (Steinbrink et al 2005), MRI of neural currents has a low sensitivity and limited temporal resolution of about 30 ms due to echo-planar data acquisition (Parkes et al 2007), and the new method of MR-EIT, in which a dc current is applied to an object and distortions in the resulting magnetic field are imaged with MRI (Sadleir et al 2006) at present does not appear to have sufficient sensitivity.
1.1.2. Cortical anatomy and evoked responses in the rat. The depth of the cerebral cortex in the rat is 1.5-3 mm. It is organized into six horizontal layers in the granular cortex of the somatosensory and visual cortical areas which are 0.2-0.6 mm thick. The cortex is specialized into functional areas, such as somatosensory or visual, which are subdivided into cortical columns each about 0.5 mm in diameter which span the cortical layers. The function of the cerebral cortex may be investigated using EEG recording during a sensory stimulus, such as mechanical whisker displacement, or electrical stimulation of the forepaw. As the evoked potential is generally obscured by the spontaneous uncorrelated cortical activity, averaging following repeated stimuli is employed to separate out the evoked activity. Typical amplitude and bandwidth of cortical evoked potentials are 100-2000 μV and 0-1000 Hz respectively (Barth et al 1993).
1.1.3. Imaging fast neural activity with EIT. A method of recording impedance changes due to ion channel opening has been previously proposed (Oh et al 2011), in which a planar subdural electrode array (30 electrodes, 0.6 mm diameter each, 1.2 mm spacing, and 5 × 7 mm 2 total covering area) was placed directly on the exposed cortex over the primary somatosensory area. Alternating current (50 μA amplitude, 225 Hz) was passed through a single pair of electrodes at a time, and voltage was recorded on the remaining electrodes. One hundred twenty averages were collected over 1 min, and this was repeated with injection through 30 other electrode pairs to produce an imaging data set. The resulting voltage differences normalized over time were then demodulated around the carrier frequency and reconstructed into images using a generic FEM of the rat's head and truncated singular value decomposition inverse algorithm.
All inverse methods introduce prior assumptions about the conductivity distribution, such as the suspected size and shape, or spatial smoothness of the perturbation (Kaipio et al 1999). These assumptions are necessary for the stability of the solution; however, it is not currently clear what assumptions can be made in the application to neural activity, because it can have any size or form as it is a macroscopic representation of millions of microscopic events (synaptic activity, action potentials, etc), occurring in seemingly random places throughout the brain. At the same time, impedance changes of 1% within the cortical volume result in at most 0.1% changes on the boundary, and so require high precision in the measurement and reconstruction techniques.

Purpose
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.

Methods
Minimizing possible modelling errors in the forward problem is critical as the inverse problem, which uses the Jacobian matrix, is ill-posed and ill-conditioned. The most critical parameters of the solution are the quality of the mesh, which can be assessed by element quality and discretization errors, and quality of the applied boundary conditions, which are correlated. All related parameters involved in the forward solution were systematically addressed and optimized using mesh convergence analysis, which is the standard technique in FEM normally employed for mesh optimization (Schmidt et al 2009). The inverse problem was then optimized by improving the conditioning of the problem and avoiding the 'inverse crime' by using a detailed 5 M element FEM mesh for forward solutions and a reduced, coarsened regular hexahedral (cubical) mesh for inversion, a similar approach to the one used by (Borsic et al 2010). Zeroth-order Tikhonov regularization was applied (Tikhonov et al 1995) in order to minimize required priors to one parameter, which was cross-validated. The resulting imaging method showed non-uniform localization in depth, so a novel image post-processing method was developed and validated on simulations performed for all possible locations of neural activity underneath the array. These methods were developed for imaging just in the cerebral cortex. The applicability of the improved method to imaging throughout the rat brain was then tested by assessment of images of simulated fast neural activity in the hippocampus.

Simulation design
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 −1 ), white matter (0.15 S m −1 ), and CSF-ventricles (1.75 S m −1 ). 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 −1 ). It was placed in all possible locations underneath the electrode array, covering 7 × 9 × 2 mm 3 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.

Forward solution
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): (1) where u is the electric potential, σ 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: The solution of equations (1) was performed with the finite element method using the open source finite element solver EIDORS (Adler and Lionheart 2006).

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.
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: Here E c is the mesh convergence error, E v is the variance of the solution obtained on nth refinement step, and V n i 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): 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.

Inverse solution
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 −1 . As J is ill-conditioned, regularization is required to perform this operation. Zeroth-order Tikhonov regularization minimizes the l2-norm magnitude of the changes: δσ = arg min δσ δv − Jδσ 2 + λ δσ 2 , δσ 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 −20 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: where P i is the subset of the elements from the tetrahedral forward mesh which belong to the ith hexahedron of the inverse mesh. 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.

Image post-processing: a novel noise-based correction approach
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 −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: δσ n = J −1 λ δv n , δv n ∼ N(0, 0.5 · 10 −6 ).
δσ 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: 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).

Accuracy and resolution
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: 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 δσ i is the conductivity value in ith element. Localization error (difference between the true and reconstructed position of the perturbation: |r c − r c | was computed and displayed, as well as shape error |ε r − ε r |, where r c is the true position of the centre of perturbation, and ε r = 0.25 is the radius of perturbation.

Feasibility of deep brain neural activity imaging
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.

Mesh convergence
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).

Simulation example
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.

Resolution and accuracy analysis
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. 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). Figure 6. Improvement of depth localization using the novel t-score masking procedure: true perturbation (top left), raw reconstructed conductivity change (bottom left), t-score mask, applied for the correction (top right), and resulting conductivity change reconstructed using novel approach (bottom right).
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 2 ), there was a uniform accuracy of 0.4 mm for uncorrected, and 0.1 mm for corrected reconstruction. The average shape error |ε r − ε r | 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 (ε 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).
The average shape error |ε r − ε r | 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.

Deep brain neural activity imaging
The 0.5 mm diameter perturbation reconstruction for the deep brain neural activity produced a localization error |r c − r c | = 0.5 mm with the shape error being |ε r −ε r | = 0.35 mm (figure 10).

Summary
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.

Technical issues
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.

Future work
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.