Stroke type differentiation using spectrally constrained multifrequency EIT: evaluation of feasibility in a realistic head model

We investigate the application of multifrequency electrical impedance tomography (MFEIT) to imaging the brain in stroke patients. The use of MFEIT could enable early diagnosis and thrombolysis of ischaemic stroke, and therefore improve the outcome of treatment. Recent advances in the imaging methodology suggest that the use of spectral constraints could allow for the reconstruction of a one-shot image. We performed a simulation study to investigate the feasibility of imaging stroke in a head model with realistic conductivities. We introduced increasing levels of modelling errors to test the robustness of the method to the most common sources of artefact. We considered the case of errors in the electrode placement, spectral constraints, and contact impedance. The results indicate that errors in the position and shape of the electrodes can affect image quality, although our imaging method was successful in identifying tissues with sufficiently distinct spectra.


Background
MFEIT (multifrequency electrical impedance tomography) is a method for imaging biological tissues with frequency-dependent conductivity. Electrodes placed on the skin are employed to inject alternating current through the body and measure the resulting boundary voltage differences. Tissues are differentiated on the basis of their conductivity spectra by varying the modulation frequency of the current. The advantage of MFEIT over standard EIT, which usually involves the observation of a change in the voltages which occurs over time, is that baseline measurements are not required. Therefore MFEIT can be used as a diagnostic tool for a multitude of applications (Malich et al 2003, Hampshire et al 1995, Brown et al 1995, Shi et al 2008, including stroke type differentiation (Horesh et al 2005, Romsauerova et al 2006, Packham et al 2012. While cerebral haemorrhagic stroke is caused by bleeding in the brain and requires surgery for treatment, ischaemic stroke is an interruption of blood flow in a region of the brain caused by an embolism and can be treated with a thrombolytic drug (rt-PA). The current diagnostic procedure is to take a CT image and results in only 2.5%-5% of the 80% of ischaemic strokes to be treated in time (Power 2004). In the case of ischaemic stroke, cell swelling caused by energy failure results in an impedance increase (Hansen andOlsen 1980, Holder 1992). In the case of haemorrhagic stroke, increased blood volume results in higher conductivity. Although EIT cannot compete directly with CT in terms of image quality, low cost and portability could make EIT scanners immediately available in the ambulance. Once a method for fast application and localization of the electrodes in emergency situations is developed, application of EIT to stroke imaging could result in fast administration of thrombolytic drugs and significantly improve the outcome of treatment.
Application of EIT to brain imaging is complicated by the presence of the skull, which is highly resistive and limits the current flowing into the centre of the head. This results in a low signal-to-noise ratio because the areas with the highest current density contribute most to the measurements. The amplitude of the injected current is limited by medical safety regulations, therefore the obtainable signal amplitude is limited. Electrode positions must be measured accurately and skin-to-electrode contact impedance is an unpredictable source of noise. MEIT imaging would allow for the subtraction of such modelling errors and, given the knowledge of the spectral properties of blood and ischaemic tissue, is suitable for stroke classification. However, a nonlinear, large-scale inversion framework is required (Horesh et al 2004(Horesh et al , 2007. A novel method for performing MFEIT using spectral constraints was recently proposed by our research group (Malone et al 2013). The inverse problem was modified by substituting conductivity with an auxiliary variable, which depends on the conductivity and the spectral information. The new variable, the volume fraction of each tissue in each voxel, describes the physical distribution of tissues in the domain, and is independent of frequency. Therefore data acquired at different frequencies can be used simultaneously and frequency independent modelling errors can be significantly reduced.
The effect of modelling errors has recently been investigated in the case of 2D timedifference EIT imaging (Boyle and Adler 2011). The results indicate that errors in the shape of the electrodes and boundary and in the contact impedance can produce artefacts in the reconstructed images. These effects are likely to be more severe in 3D multifrequency imaging. In this paper, the first application of the fraction reconstruction algorithm to a realistic threedimensional model of the head with skull and scalp is demonstrated. The influence of imprecise modelling is evaluated in three different cases: electrode positions, electrode contact impedance and tissue conductivity spectra. This simulation study illustrates that our new imaging method might be used to differentiate between different types of stroke in real measurements and shows where changes in the measurement paradigm are necessary.

Forward and inverse problem
The forward problem consists in determining the expected boundary voltage data for a given object with known conductivity and is solved using the finite element method (FEM). The commonly used complete electrode model (CEM) accounts for two observed effects, electrodes shunting current due to their high conductivity and a voltage drop at the interface of the electrodes and the skin, which is due to an electrochemical effect (Somersalo et al 1992).
The inverse problem consists in determining the conductivity of an object σ from knowledge of the boundary voltages v. The problem is solved by minimizing a regularized objective function, which depends on the modulation frequency of the injected current ω. Assuming the noise is Gaussian, the objective function assumes the form where is a regularizing function, and τ is the regularization parameter.

Purpose
The purpose of the study in this paper is to evaluate the robustness of our new current imaging method to various sources of error. This reproduces the conditions of a real experiment, and helps us assess if the method is suitable for application to human subjects. Specifically, we are interested in answering the following questions.
• What is the effect of using a coarse mesh for the reconstruction?
• What is the effect of adding errors to the position of the electrodes (thereby also changing the area and shape)? • What is the effect of adding errors to the assumed spectral information?
• What is the effect of adding errors to the contact impedance of the electrodes?

Experimental design
We construct a numerical head phantom with homogeneous layers for the brain, skull and scalp. The meshes and surfaces are obtained from a CT scan of a human head. For simplicity, in this first feasibility study the model does not include the cerebro-spinal fluid, a simplification which is commonly made in EIT research. Furthermore, the electrodes are placed in the same positions used to acquire EEG measurements on the scalp. The advantage of this setup is that electrode caps and other equipment intended for EEG applications can be used. Realistic conductivities for all tissues are taken from the literature for a range of frequencies. Two tetrahedral finite element meshes with different resolution, one coarse and one very fine, are generated in order to avoid the inverse crime (Lionheart 2004). The fine mesh is used to simulate the boundary voltage data, and the coarse mesh to reconstruct images. The injecting pairs of electrodes are chosen to maximize the distance between the electrodes, while acquiring the maximum number of independent measurements. This is achieved by finding the maximum spanning tree of the electrodes, weighted by the distance between the electrodes. In order to simulate the most common sources of artefact in an experimental setup, errors are added to the model. For each case an EIT image is reconstructed using our spectrally constrained fraction reconstruction method. The images are evaluated and compared using an objective image quality quantification method.
We consider the modelling errors that are the most common sources of artefacts in EIT images obtained experimentally (McEwan et al 2007). The error levels are chosen such that they can be achieved by a modern experimental setup.
• Instrumentation noise is chosen to match that of measurements acquired using the KHU Mark 2.5 EIT system (Oh et al 2011) in a saline filled tank and averaged over 64 frames. This noise level is achievable with most EIT measurement systems and can be reduced by use of better instrumentation. The standard deviation of the proportional noise was ς p = 0.02% and the standard deviation of the additive noise was ς a = 5 μV. These values illustrate that the additive noise dominates EIT measurements. • Electrode positions can currently be measured to around 1 mm precision using photogrammetry (Qian and Sheng 2011). Other technologies, such as the commercial MicroScribe, laser 3D scanners, or electrode helmets, can achieve an even higher precision in electrode localization. To demonstrate the importance of using such localization technologies, we simulated electrode position errors of around 1 mm and 2 mm. These relatively small errors result in a remarkable degree of image degradation. Given that the electrodes are represented on a discrete mesh, the shape and size of the electrodes also change when an error is added to the position of the centre. This could be accounted for by refining the mesh, however we decided not to do this in that a coarse representation of the electrodes constitutes an unpredictable source of errors, and thus provides a greater similarity between simulation and experiment (Kolehmainen et al 1997, Boyle andAdler 2011). Work is in progress to investigate the separate contributions of the errors caused by the shape, size and position of the electrodes. Errors are added to the (x, y, z) positions of all electrodes before simulating the data. Deviations of up to three times the standard deviation ς of the error are expected in the majority of cases. Therefore the overall shift of the centre of each electrode will normally be less than or equal to where (x,ỹ,z) is the position of the shifted electrode. For the errors chosen, we have -√ 3(3 × 0.25) ≈ 1.3 mm for a standard deviation of 0.25 mm; -√ 3(3 × 0.5) ≈ 2.6 mm for a standard deviation of 0.5 mm.
• Prior spectral information is imprecise because tissue conductivity is influenced by anisotropy, inhomogeneity and temperature. Because the resulting effect of these factors is difficult to predict, we simulate errors based on the literature that roughly represent their frequency-dependent contribution (Edd et al 2005). To test the limitations of our reconstruction method, we consider a reasonable and a 'worst-case' level of error: 1% and 5%. The errors are added independently to each frequency and each tissue type.
It is important to note that the multifrequency reconstruction algorithm used in this study is insensitive to conductivity changes with a flat frequency-spectrum. Therefore we only have to consider frequency-dependent errors, which constitute a small fraction of the above mentioned error sources. • Contact impedance errors are chosen on the basis of our experience of conducting experiments. It is assumed that all electrodes have sufficiently low contact impedance.
In an experimental setup, this is equivalent to discarding any electrodes with near-infinite impedance, which may have detached from the head, or any broken measurement channels. If the variance of the contact impedance across electrodes is approximately 20%, the setup is considered suboptimal. If the variance of the contact impedance is larger than 50%, then the electrodes have to be re-applied. We thus choose these two levels of error to illustrate that the contact impedance mostly influences the instrumentation rather than the modelling and image reconstruction. This is because a high transfer impedance introduces large noise levels in the instrumentation.
We acknowledge that the number of test cases used in this paper does not constitute a complete test over all possible permutations. The intention of this paper is to show a proof of principle for MFEIT stroke type detection. The number of averages chosen for each simulated error source is sufficient to illustrate which effects have strongest contribution to image artefacts.

Model and tissue impedance spectra
A three-dimensional model of a human head was used to simulate EIT data. The model comprised three layers, corresponding to the scalp, skull, and brain. A fine mesh, with 5 million elements, was used to simulate the boundary voltages, and a coarse mesh, with 180 thousand elements, was used to reconstruct the images. A spherical perturbation of diameter 3 cm was placed in two different positions inside the brain: lateral and posterior (figures 1(b) and (c)). The conductivity spectra of the tissues (scalp, skull, brain, ischaemic brain, blood) were obtained from the literature (Romsauerova et al 2006, Horesh et al 2005 ( figure 1(a)). In order to simulate an ischaemic stroke, the conductivity of the perturbation was set to the conductivity of ischaemic brain approximately one hour after onset, and in order to simulate a haemorrhagic stroke, the conductivity of the perturbation was set to the conductivity of blood. Twelve frequencies were chosen in the range 5 Hz-5 kHz based on the observation that the slopes of the tissues are most different in this region (figure 1(a)), and the boundary voltages were simulated for each frequency.

Data simulation
Thirty-two electrodes of diameter 10 mm were placed on the surface of the model. The electrodes were modelled using the CEM, and the contact impedance was set to 1 k for all electrodes. The peak-to-peak amplitude of the current was set to 140 μA. Voltage measurements were made for each injection on all adjacent pairs not involved in delivering the current. The total number of measurements acquired for each frequency was 869.
The boundary voltages were computed with the recently developed parallel EIT solver (PEITS) (Jehl et al 2014) written in C++ using the DUNE-FEM package (Dedner et al 2010). PEITS is a finite element solver specifically written for forward simulation in EIT using the CEM. It was designed to have a very good performance on multi-core machines and clusters using MPI and thus reduce the computation time for forward simulations. PEITS partitioned the mesh into largely independent parts on which the weak formulation of the CEM was assembled: We applied a Dirichlet boundary condition on one surface node of the mesh by setting it to 0 V to make the system uniquely solvable. The system was then solved iteratively with a conjugate gradient solver which was preconditioned by the algebraic multigrid implementation ML (Gee et al 2006). Using PEITS on all 16 cores of a workstation with two 2.4 GHz Intel Xeon CPUs with eight cores and 20MB cache each, reduced the computation time for 31 forward solutions on the fine 5 million element mesh to less than 2 min.

Image reconstruction
Images were reconstructed using the fraction reconstruction method (Malone et al 2013). For all frequencies {ω i ; i = 1, . . . , M} and tissues t j ; ∀ j = 1, . . . , T , the conductivity of a tissue sample i j = σ t j (ω i ) was assumed to be known. The conductivity of an object at a certain frequency σ n (ω i ) was modelled as the linear combination of the conductivities of individual tissues i j ; j = 1, . . . , T : where n runs over the elements of the mesh, 0 f n j 1, and T j=1 f n j = 1. The weightings f n j of the linear combination are called fractions, and are defined for each tissue and element. The fractions represent the physical distribution of the tissues, and are independent of the frequency at which the data is acquired.
The objective function (1) was expressed in terms of the fractions by substituting the model for the conductivity (4): where f j = f n j ; n = 1, . . . N ∈ R N×1 and F = f j ; j = 1, . . . M ∈ R N×M . Given that the reconstructed parameter is frequency independent, all multifrequency measurements can be considered simultaneously. Furthermore, changes in boundary voltages across frequencies can be taken, rather than the absolute data, in order to suppress frequency independent modelling errors. The resulting objective function is A Markov random field regularization term of the form was chosen, where l(n) runs over all neighbours of the nth voxel. The fractions were recovered simultaneously for all tissues and elements by minimizing the objective function (F). The latter is differentiable in the fractions, and the gradient was obtained via the chain rule. The reconstruction of the fractions was constrained to the closed interval [0, 1], and the constraint T j=1 f n j = 1 was enforced by substituting f 1 = 1− T j=2 f j in the objective function. The minimization was performed by alternating steps of gradient projection and damped Gauss-Newton algorithms (Nocedal and Wright 1999).

Numerical validation
In order to validate the method, images were reconstructed from simulated data without the addition of modelling errors except those due to mesh discretization and noise. The data were simulated using the fine mesh and the images were reconstructed using the coarse mesh. In order to simulate instrumentation error, both proportional and additive noise were added to the data: where rand(ς ) indicates a random number drawn from a Gaussian distribution with zero mean and standard deviation ς . The standard deviation of the proportional noise was ς p = 0.02% and the standard deviation of the dominating additive noise was ς a = 5 μV. The skull and scalp were fixed in place, and it was assumed that the area inside the skull was occupied by either the brain, or the stroke. The initial guess was the healthy brain. The optimal regularization parameter was approximated by computing the L-curve for one step of Gauss-Newton descent. The corner of the L-curve was selected for the first step of the reconstruction, and the value was divided at each step by 2 for the ischaemic stroke and, given that the contrast was lower, by 3 for the haemorrhagic stroke (Viklands and Gulliksson 2001). The automatic selection of the regularization parameter was repeated in all the following cases, and the maximum number of steps was fixed to ten.

Error simulation
Modelling errors were simulated by altering the model used to simulate the voltages. Three types of errors were considered: erroneous electrode positions, erroneous tissue spectra, and erroneous contact impedance values. The position errors were added to the (x, y, z) coordinates of the electrodes separately. The error in conductivity spectra was added to each tissue at each frequency individually. The error in the contact impedance was added to each electrode separately.
In order to evaluate the response of the fraction reconstruction method to different errors, the study was repeated for normally distributed errors with two different levels of variance. The following cases were considered, • electrode positions: standard deviation 0.25 mm and 0.5 mm on (x, y, z) coordinates separately; • tissue conductivities: standard deviation 1% and 5% on i j , where i indicates the frequency and j the tissue; • contact impedance: standard deviation 20% and 50% from 1 k .
Additionally, proportional and additive noise was added to each data set.

Image quantification
An image quantification method was devised to evaluate the images objectively. The quality of an image was assessed in terms of the ability to distinguish an anomaly (the stroke) from a background (the brain). For this reason the fraction f s corresponding to the tissue that makes up the anomaly was assessed. The area P corresponding to the reconstructed perturbation was identified as the largest connected cluster of voxels with values larger than 50% of the maximum of the image ( (i) Image noise: inverse of the contrast-to-noise ratio between the perturbation P and the background B wheref P s andf B s are the mean intensities of the perturbation and background and std is the standard deviation. (ii) Localization error: ratio between the norm of the x-y-z displacement of the centre of mass of the reconstructed perturbation P from the actual position (x, y, z) true , and the norm of the dimensions of the mesh (d x , d y , d z ) where (x n , y n , z n ) is the position of the centre of the nth tetrahedron. (iii) Shape error: mean ratio of the difference between the dimensions of the simulated and reconstructed perturbations, respectively (l x , l y , l z ) and (l x , l y , l z ), and the dimensions of the mesh 1 3

Numerical validation
The data were simulated on the fine mesh, noise was added to the data, and the images were reconstructed on the coarse mesh. For comparison, the process was repeated using data simulated on the coarse mesh (figure 2). The discretization errors introduced differences in the area of each electrode between the fine (5 million elements) and coarse (180 thousand elements) meshes. The difference in the area of the electrodes between the fine and coarse mesh was, in average 5.6 ×10 −7 m 2 , over an average electrode area of 7.7 ×10 −5 m 2 , i.e. about 1.4% (figure 4). Image quantification measures were computed for each of the reconstructed images ( figure 3). The images obtained from data simulated and reconstructed on the same mesh are superior in terms of the previously defined measure of image quality (section 2.6) to those obtained from data simulated on the fine mesh. The contrast recovered in the images relating to ischaemic stroke is greater than of the images relating to haemorrhagic stroke. The images obtained for the posterior position are in most cases better than for the lateral position.

Erroneous electrode positions
The images were reconstructed assuming that the electrodes were fixed in the original positions ( figure 5) and image quality measures were computed for each image (figure 6).  The perturbation was recovered only in the case of 0.25 mm standard deviation error added to ischaemic stroke data. In all other cases the images quality is deteriorated to the point that the imaging must be considered unsuccessful.

Erroneous tissue spectra
Errors were added to the conductivity values of the model to simulate a deviation from the literature values. The images were reconstructed using the original values for the conductivities of the brain and stroke (figure 7), and image quality measures were computed for each image (figure 8). The perturbation was recovered successfully in all cases for 1% error, but in the case of 5% error only the ischaemic stroke was identified correctly. Figures 9(a) and (b) display the frequency-difference spectra for brain, ischaemic brain and blood, with the associated error bars. The variance of the error on the relative spectra is given by the sum of the variance of the errors added to the absolute values, and the error bars indicate the minimum and maximum limit within the majority of the errors are drawn, given by ±3 times the standard deviation.

Erroneous electrode impedances
Images were reconstructed assuming a value of 1 k for the contact impedance of all electrodes (figure 10), and image quality measures were computed for each image ( figure 11)   are nearly unchanged by the introduction of 20% errors on the contact impedance, and image quality is slightly diminished for 50% errors.

Numerical validation
The contrast obtained in the images of ischaemic stroke was higher than in the images of haemorrhagic stroke. This may be attributed to the difference in the impedance spectra of ischaemic brain and blood. The difference between the slope of the conductivity spectrum for ischaemic brain and healthy brain ( figure 1(a)) is greater than the difference between blood and healthy brain. Therefore, given that the method uses data referred to the lowest frequency,   Figure 11. Image quantification results for images reconstructed with errors of 20% and 50% standard deviation added to the electrode contact impedances: (a) ischaemic stroke, (b) haemorrhagic stroke. the signal given by ischaemic stroke is greater than that of a haemorrhagic stroke of the same size and in the same location.
The reduction in image quality between the case of data simulated on the fine and coarse mesh was primarily caused by the modelling of the electrodes and skull. Given the different resolutions, the shape and size of the electrodes and the skull differ between the two meshes. The purpose of using a fine mesh to simulate the data, and a coarse mesh to reconstruct the image, is to make the simulation more realistic (figure 3). In the case of imaging a real human head, the size and thickness of the skull will not be known exactly. Furthermore, the discrete representation of the electrodes and skull on the mesh does not reflect the reality precisely. Therefore it was necessary to consider these discrepancies in a realistic simulation. If the same mesh is used to simulate and reconstruct the data the problem is over-simplified with respect to the real-case scenario, and conclusions drawn from simulation results may not be applicable in practice (figure 2).

Erroneous electrode positions
Errors added to the electrode positions severely affect image quality. This highlights the importance of registering the position of the electrodes accurately. Furthermore, in order to preserve the shape and size of the electrodes, the mesh must be sufficiently refined on the boundary (figure 4).

Erroneous tissue spectra
The fraction reconstruction method requires knowledge of the conductivity spectra of all tissues, and these are assumed to be fixed and exact. The performance of the algorithm is therefore diminished if the assumed spectral constraints are incorrect ( figure 8). Furthermore, the confidence with which the reconstruction algorithm distinguishes between different tissues depends on the difference between the conductivity spectra of the tissues. Specifically, given that frequency difference data is used, the tissues are distinguished on the basis of the respective change in conductivity between the lowest and the other frequencies. In order to determine if it is possible to distinguish between healthy and ischaemic brain, or between brain and blood, we investigated the effect of adding uncertainty to the spectra. If a random error is added to the absolute spectrum, then the error on the difference in the spectrum with respect to the lowest frequency is given by the sum of the absolute errors. For 1% error, all the spectra are distinct, but for 5% error, the spectra overlap for some or all frequencies (figures 9(a) and (b)). For this reason it was not possible to locate the haemorrhagic stroke in the case of 5% error added to the conductivities (figure 7(c)), and the ischaemic stroke was only identified in the lateral position ( figure 7(a)). In the case of haemorrhagic stroke the addition of a proportional 5% error caused a large degree of uncertainty because the absolute value of the conductivity of blood is large. In the case of ischaemic stroke, the uncertainty was caused by the similarity in the spectra of healthy and ischaemic brain.

Erroneous electrode impedances
The effect of the errors added to the contact impedance is very limited. A change in the contact impedance will cause a change in the current distribution around the electrode. However, given that the conductivity of the electrode is very high relative to the conductivity of the object, changes in the electrode impedance have a small effect on the current flow inside the object. For this reason the images obtained after adding errors to the contact impedance are similar to the original images without modelling errors.

Technical remarks
Ideally, several images would have been created for each noise level in order to characterize the effect of modelling errors over a very large number of samples. The computational expense of multiple repetitions was prohibitive, in that reconstruction of a single image took 5-6 h. However, given that the electrode specific errors (contact impedance and position) were sampled on the 32 electrodes individually, this provides a sufficiently large number of samples to give a reasonable characterization of the influence of the noise. Likewise, in the case of errors added to the tissue spectra, we added the noise independently to each tissue and at each frequency, and this allows us to describe the effect of the spectral error reasonably well. Thus, the conclusions derived from this relatively small number of images appear to be valid in principle. In future, examination of more permutations in simulation and tank studies may allow us to define quantitative limits to the acceptable variation of each parameter.

Conclusion
We have applied the fraction reconstruction method using spectral constraints to a numerical head phantom with realistic conductivities. We have added noise and modelling errors to investigate the robustness of our method. The results show a varying degree of sensitivity of the method to different sources of error. We conclude that • the method is highly sensitive to errors in the position and shape of the electrodes, and these must be modelled with the highest achievable level of accuracy; • the fraction reconstruction method allows to distinguish between tissues if the respective spectra are sufficiently distinct; • the method is highly robust to errors in the assumed contact impedance.