Open Access
30 January 2024 Simulation-based feasibility study of monitoring of intracerebral hemorrhages and detection of secondary hemorrhages using electrical impedance tomography
Jussi Toivanen, Antti Paldanius, Bachir Dekdouk, Valentina Candiani, Asko Hänninen, Tuomo Savolainen, Daniel Strbian, Nina Forss, Nuutti Hyvönen, Jari Hyttinen, Ville Kolehmainen
Author Affiliations +
Abstract

Purpose

We present a simulation-based feasibility study of electrical impedance tomography (EIT) for continuous bedside monitoring of intracerebral hemorrhages (ICH) and detection of secondary hemorrhages.

Approach

We simulated EIT measurements for six different hemorrhage sizes at two different hemorrhage locations using an anatomically detailed computational head model. Using this dataset, we test the ICH monitoring and detection performance of our tailor-made, patient-specific stroke-monitoring algorithm that utilizes a novel combination of nonlinear region-of-interest difference imaging, parallel level sets regularization and a prior-conditioned least squares algorithm. We compare the results of our algorithm to the results of two reference algorithms, a total variation regularized absolute imaging algorithm and a linear difference imaging algorithm.

Results

The tailor-made stroke-monitoring algorithm is capable of indicating smaller changes in the simulated hemorrhages than either of the reference algorithms, indicating better monitoring and detection performance.

Conclusions

Our simulation results from the anatomically detailed head model indicate that EIT equipped with a patient-specific stroke-monitoring algorithm is a promising technology for the unmet clinical need of having a technology for continuous bedside monitoring of brain status of acute stroke patients.

1.

Introduction

Stroke is the second leading cause of death and a major cause of disability worldwide.1 Intracerebral hemorrhage (ICH) accounts for about 20% of all strokes.1 ICH is associated with higher mortality than ischemic stroke especially at the early stage.2 ICH is typically caused by disruption of cerebral arteries, which leads to bleeding into brain tissue. ICH harms the brain via several different mechanisms,3 including the space-occupying mass effect in the acute phase, the toxic effect of the degradation products from the lysed red blood cells, and the inflammatory changes in the later stages. So far, there is no established outcome-modifying medical or surgical treatment for the index ICH. However, one of the major challenges in the acute stage (within 36 h from onset of ICH) is the risk of hematoma expansion and rebleeding.2 In such cases, intensive conservative treatment, including controlled mechanical ventilation and optimized blood pressure levels, and timely neurosurgical intervention are likely to be life-saving and even reduce disability. Such intensive care is not risk free, and it should be reserved for patients that show signs of hematoma growth and changes in the level of consciousness. In addition to patients with hemorrhagic stroke, patients with thrombolysis-treated acute ischemic stroke have an increased risk of cerebral hemorrhage as a complication of recanalization therapy. In both patient groups, the decision to advance to intensive care has to be done promptly, and therefore, hematoma expansion or rebleeding in ICH patients should be detected as soon as possible, optimally immediately.

The monitoring of hemorrhagic stroke in the intensive care unit (ICU) is based mainly on routine follow up of clinical signs and symptoms.2 However, evaluation of the vigilance of the patient, as well as observing subtle neurological signs can be very challenging, particularly if the patient is sedated and intubated. Currently, the most reliable way to monitor the progression of hemorrhagic stroke is repeated CT scanning2 that requires the patient to be repeatedly moved from ICU to CT and back. This is a demanding and time-consuming procedure and can even be dangerous for the patient. Rarely, in some individual stroke centers, moving the patient can be avoided if a bedside head CT scanner is available. Regardless, the selection of correct timing of the control CT imaging is itself already a difficult problem and CT images provide only occasional snapshots of the bleeding. A method for online bedside monitoring of hemorrhagic stroke at the ICU would most likely positively affect patients’ prognosis as there would be no delays in the detection of life-threatening complications and their care.

One potential bedside monitoring method is provided by electrical impedance tomography (EIT). In the EIT measurement setup, small alternating currents are fed through electrodes attached to the patient’s scalp, and the resulting voltages are measured on the electrodes. This safe and radiation-free measurement is then used to compute a three-dimensional image of the electrical conductivity of the brain, with hemorrhage showing as increased intensity with respect to normal brain tissue in the image due to the high conductivity of blood. This type of absolute imaging EIT has been investigated for early differentiation of stroke types, see e.g. Malone et al.,4 Goren et al.,5 Horesh et al.,6 Candiani et al.7 for absolute imaging and Agnelli et al.,8 McDermott et al.,9 Candiani et al.10 for machine learning-based approaches. However, EIT-based stroke differentiation is very challenging due to the instability of the absolute imaging EIT problem with respect to modeling errors, such as inexact electrode locations and head shape. Furthermore, for differentiation, a very high sensitivity is required because even very small quantities of blood should not be overlooked in clinical decision making. In a continuous bedside monitoring setup, the EIT measurement is repeated at a later time, and a three-dimensional image of the conductivity change of the brain between the measurements is computed. These types of difference imaging EIT approaches are promising for monitoring of hemorrhagic stroke and detection of secondary hemorrhages as they have good sensitivity, as suggested by simulation studies by Shi et al.,11 animal studies by Xu et al.,12 and human studies by Dai et al.13 and Yang et al.,14 and they are less prone to modeling errors.

In the envisioned EIT-based stroke monitoring setup, the patient is admitted to the emergency department of the hospital because of symptoms of acute stroke. According to stroke treatment protocol, a CT scan is taken immediately for diagnosis and timely beginning of suitable treatment. The treatment of the patient continues in the ICU or in a stroke unit, where the EIT measurement device is attached using something akin to an EEG electrode cap. The EIT measurements, taken, e.g., every 2 to 10 min, and an imaging algorithm are then used to get snapshots of the conductivity changes in the brain, indicative of changes in the hemorrhage, which would indicate any immediate need for confirmatory CT imaging or medical intervention. The EIT-based stroke imaging could also be used in a similar manner for monitoring of occurrence of secondary hemorrhages in treatment of ischemic stroke patients. Furthermore, the CT image routinely taken at patient admission can be utilized for generation of a 3D computational model of the patient’s head and for regularization purposes in the image reconstruction.

In this paper, we study the feasibility of EIT for monitoring and detection of ICHs with simulated measurement data from an anatomically highly detailed head model. We compare the performance of several inverse estimation methods for the monitoring task. We use the algorithm first introduced by Toivanen et al.15 that has been constructed as a combination of the nonlinear region-of-interest difference imaging approach by Liu et al.16 and the parallel level sets regularization approach by Kolehmainen et al.17 In addition, it now also utilizes a prior-conditioned least squares algorithm, see Arridge et al.18 and Harhanen et al.,19 for computationally efficient solution of the lagged Gauss–Newton search direction for the minimization of the regularized nonlinear least squares functional. Furthermore, a more advanced initial estimation approach is used. We use an atomically detailed six-layered head model with intricate cerebrospinal fluid and brain geometry that was first introduced by Paldanius et al.20 to simulate EIT measurements with various ICH locations and volumes to test our stroke monitoring and detection algorithm. The results with our new monitoring algorithm are compared to the results of two reference algorithms, a total variation (TV) regularized absolute imaging algorithm and a linear difference (LD) imaging algorithm.

2.

Theory

2.1.

Modeling of EIT Measurements

We model the patient’s head in the EIT measurement setup as a domain ΩR3, and the L electrodes attached to its surface with circular surface patches e, l=1,2,,L. During the measurement, P patterns of currents, I(k)RL, k=1,2,,P, are consecutively injected through the electrodes, and the corresponding voltages U(k)RL are measured on all electrodes. Here I(k) and U(k) denote the applied current and measured voltage from the k’th current pattern on the ’th electrode for =1,2,,L. Based on the conservation of charge and our choice of electric potential ground, we have

Eq. (1)

=1LI(k)=0=1LU(k)=0.

The voltages U(k) are boundary measurements of the interior electromagnetic potential u(k)(x) that is modeled with the conductivity equation

Eq. (2)

·(σ(x)u(k)(x))=0,  xΩ
and the boundary conditions of the complete electrode model (CEM)21,22 for k=1,,P and =1,,L

Eq. (3)

u(k)(x)+zσ(x)u(k)(x)n=U(k),  xe,

Eq. (4)

eσ(x)u(k)(x)ndS=I(k),

Eq. (5)

σ(x)u(k)(x)n=0,  xΩ=1Le,
where z is the contact impedance between the electrode e and the body Ω, n denotes the outward unit normal vector on the boundary Ω, and the isotropic conductivity distribution σ is assumed to belong to L+(Ω){ςL(Ω)|essinfς>0}. The existence and uniqueness of the solution (u,U)H1(Ω)RL/R of the model [Eqs. (1)–(5)] were proven and its variational form derived by Somersalo et al.22

In this paper, for the imaging algorithms, the numerical solution of the model [Eqs. (1)–(5)] is based on the finite-element method (FEM), for details of the implementation see Vauhkonen et al.23 and Kaipio et al.24 In the following, we denote the FEM-based solution for a single current pattern I(k)RL by σU(σ;I(k))RL where we consider a discretized version of σ, such that σ=j=1Nσjφj, where σjR+ are the nodal coefficients, j=1,,N, and φjH1(Ω), j=1,,N, are the piecewise linear basis functions corresponding to an FE mesh of Ω. Similarly for other variables, we use identical notation for continuous and discretized forms, assuming the correct interpretation is clear from the context.

The measurement noise eRP·L is modeled as additive noise, leading to the measurement model

Eq. (6)

V=U(σ)+e,
where VRP·L is the vector of the measured noisy voltages for all applied current patterns and U(σ)=(U(σ;I(1)),,U(σ;I(P)))TRP·L.

2.2.

Algorithms

In EIT stroke monitoring, the occurrence of a secondary hemorrhage during treatment of ischemic stroke or the growth of a hemorrhagic stroke is detected as a change in conductivity

Eq. (7)

δσ=σ2σ1
between two consecutive measurement times t1 and t2. Because the hemorrhage-related changes are slow compared to the duration of a single EIT measurement set, it is reasonable to model the measurements at t1 and t2 with the stationary model Eq. (6) as

Eq. (8)

V1=U(σ1)+e1,V2=U(σ2)+e2.

With this setup, we test the feasibility of EIT for monitoring and detection of ICH by comparing the performance of three algorithms for obtaining δσ as follows.

  • (1) A total variation (TV) regularized absolute imaging algorithm that reconstructs σ1 and σ2 independently and obtains δσ by subtraction.

  • (2) A linear difference imaging algorithm (LD) that directly reconstructs δσ.

  • (3) A stroke-monitoring algorithm (MO) that reconstructs σ1 and δσ utilizing a novel combination of nonlinear region-of-interest difference imaging,16 parallel level sets regularization,17 and a prior-conditioned least squares algorithm for solution of the lagged Gauss–Newton search direction in the minimization of the regularized nonlinear least squares functional.18,19

2.2.1.

Absolute imaging-based algorithm

In the TV regularized absolute imaging algorithm, δσ is obtained by first solving two separate absolute imaging problems to obtain estimates of σ1 and σ2 and then computing δσ=σ2σ1. The absolute imaging problem is solved with one of the most popular reconstruction methods, the generalized Tikhonov regularization:

Eq. (9)

σ^t=argminσt>0{Le(VU(σt))2+pσ(σt)},
where t=1,2, Le is a Cholesky factor of the noise precision matrix, i.e., LeTLe=Γe1, and pσ(σt) is a regularization functional, which is in this algorithm chosen as smoothed TV regularization25

Eq. (10)

pσ(σt)=TV(σt)=αΩ(σt2+β2)1/2dx,
where α>0 is the regularization weight coefficient, σ is the gradient of the conductivity σ, and β>0 is a small smoothing parameter that ensures differentiability. The minimization problem [Eq. (9)] is solved with a lagged Gauss–Newton method equipped with a line search algorithm that also enforces the nonnegativity σ>0. To obtain a starting point σt(0) for the iteration, the best fitting constant conductivity can be estimated by solving a nonlinear least squares fitting problem.

Although the absolute imaging approach is generally very sensitive to geometric modeling errors, it can be expected to be a feasible choice in the stroke monitoring setup, where the patient-specific head geometry is available from the patient CT taken for diagnosis of the stroke and the geometry of the domain does not change during the monitoring.

2.2.2.

Linear difference imaging algorithm

In linear difference imaging, see Barber et al.26 and Bagshaw et al.,27 the aim is to reconstruct the change in conductivity between measurements, V1 and V2. In the linear difference imaging algorithm of this paper, the measurement models [Eq. (8)] are linearized at some conductivity σ0 using the first-order Taylor approximations

Eq. (11)

V1U(σ0)+J(σ1σ0)+e1,

Eq. (12)

V2U(σ0)+J(σ2σ0)+e2,
where the Jacobian matrix J is evaluated at σ0. A reasonable choice for the linearization point σ0 can be obtained by solving a nonlinear least squares fitting problem for the best fitting constant conductivity using the measurement data V1. With the linearized models Eqs. (11) and (12), the difference in measurements can be written as

Eq. (13)

δV=V2V1=(U(σ0)+J(σ2σ0)+e2)(U(σ0)+J(σ1σ0)+e1)=Jδσ+δe,
where δσ=σ2σ1 and δe=e2e1. The linear difference imaging problem is to reconstruct δσ based on the difference data δV and the solution can be obtained in the form of the generalized Tikhonov solution [Eq. (9)] as

Eq. (14)

δ^σ=argminδσ{Lδe(δVJδσ)2+pLD(δσ)},
where Lδe is the Cholesky factor of the noise precision matrix of δe so that LδeTLδe=Γδe1=(Γe1+Γe2)1. The linear regularization functional pLD(δσ) is in this paper a smoothness regularization implemented utilizing a distance-based correlation model, giving

Eq. (15)

pLD(δσ)=Lpσ2,
so that LpTLp=Γp1 where the covariance matrix Γp is constructed using distance-based correlation:28

Eq. (16)

Γp(i,j)=std(σ)2exp(xixj22a2),
where i,j=1,,N, and the parameter a controls the correlation length and can be solved by setting the distance xixj to a selected value d (e.g., half the radius of the target) and setting Γp(i,j) to the desired covariance for that distance (e.g., 1% of the variance).

One of the main reasons for the popularity of linear difference imaging is that in the difference measurements [Eq. (13)] usually at least part of the systematic modeling errors are cancelled out, making the approach more tolerant to modeling errors compared to most absolute imaging-based approaches. However, as the approach is based on a linear approximation of the nonlinear problem, it can lead to suboptimal results if the change between the states is large or, as usually in practice, the initial state is unknown, implying that σ0σ1 and the linearization is carried out in a nonoptimal point.

2.2.3.

Stroke-monitoring algorithm

The stroke-monitoring algorithm has been constructed as a novel combination of the nonlinear region-of-interest difference imaging approach by Liu et al.16 with the parallel level sets regularization approach by Kolehmainen et al.17 and it utilizes a prior-conditioned least squares algorithm similar to Arridge et al.18 and Harhanen et al.19 for computationally efficient solution of the lagged Gauss–Newton search direction in the minimization of the regularized nonlinear least squares problem.

The nonlinear difference imaging approach has been shown to be tolerant to modeling errors to at least the same extent as the linear difference imaging approach and it can produce quantitatively more accurate reconstructions of the conductivity change.29,30 This results from solving the full nonlinear EIT problem instead of the linearized one and from the parameterization of the problem with respect the initial conductivity state and the conductivity change. With this parameterization, the effects of the modeling errors that are invariant between the measurements V1 and V2 were found by Liu et al.30 to induce errors in the reconstruction of the initial state σ1, leaving the reconstruction of the conductivity change unaffected. Furthermore, in nonlinear difference imaging, it is possible to employ an a priori information-based region of interest (ROI) constraint for the conductivity change δσ so that

Eq. (17)

supp(δσ)=ΩROIΩ.

In monitoring of ICH a natural, not too constrictive ROI is the brain volume. If no ROI is desired, it is possible to use ΩROI=Ω, but restricting the ROI to a smaller subdomain has been shown to improve the reconstruction of δσ.16 Based on Eq. (17), the conductivity at the later measurement time t2 is modeled as

Eq. (18)

σ2=σ1+Kδσ,
where K is an extension mapping that extends the conductivity change from the ROI to the whole domain Ω, so that

Eq. (19)

Kδσ={δσ,xΩROI0,xΩΩROI.

For the simultaneous estimation of σ1 and δσ, the measurement data V1 and V2 are combined into a single vector, leading to the measurement model

Eq. (20)

[V1V2]=[U(σ1)U(σ1+Kδσ)]+[e1e2].

This measurement model can be written as

Eq. (21)

V˜=U˜(σ˜)+e˜,
where

Eq. (22)

V˜=[V1V2],U˜=[U(σ1)U(σ1+Kδσ)],

Eq. (23)

σ˜=[σ1δσ],e˜=[e1e2],
and being now in a form similar to the measurement model Eq. (6), the joint reconstruction of σ1 and δ can be defined in the form of the generalized Tikhonov regularization [Eq. (9)] as

Eq. (24)

σ˜=argminσ˜{L˜e(V˜U˜(σ˜))2+p(σ˜)},
where the diagonal blocks of L˜e contain the Cholesky factors of the noise precision matrices of measurements V1 and V2, and the regularization functional

Eq. (25)

p(σ˜)=pσ1(σ1)+pδσ(δσ)
allows independent regularization models for δσ and σ1. The conductivity change caused by stroke expansion is expected to be localized and regular, and thus we use the smoothed TV regularization25

Eq. (26)

pδσ(δσ)=TV(δσ)=αδσΩ(δσ2+β2)1/2dx,
where αδσ>0 is the regularization weight coefficient, δσ is the gradient of the conductivity change, and β>0 is a small smoothing parameter that ensures differentiability. The initial conductivity σ1 is expected to correlate well with the structure of the patient CT that is always taken for diagnosis of the stroke at the time of patient admission to the hospital. This information can be utilized using a parallel level sets-based, spatially and directionally weighted TV regularization that promotes similar alignment of level sets in σ1 and the CT-based reference image,17 giving

Eq. (27)

pσ1(σ1)=WTV(σ1)=ασ1Ω(σ1B(κ)2+β2)1/2dx,
where ασ1>0 is the regularization weight coefficient, κ(x) is the reference image and the tensor B(κ) is chosen such that aligned edges (or level sets) in σ1(x) and the reference image κ(x) are preferred. One possible choice for the weighting tensor is17

Eq. (28)

B(κ)=I(1γ(x))ν^(x)ν^(x)T,
where

Eq. (29)

ν^(x)={0,if  κ(x)=0κ(x)/κ(x),otherwise,
and 0γ(x)1 is an edge indicator function such that γ(x)1 when κ(x)0. In this paper, γ(x) is chosen by thresholding the gradient of the reference image so that

Eq. (30)

γ(x)={γ1,if  κ(x)threshold1,otherwise,
where γ11.

The minimization problem [Eq. (24)] can be solved using a lagged Gauss–Newton method with the positivity constraints σ1>0 and σ1+Kδσ>0. However, solving of the Gauss–Newton search direction is very memory intensive and slow given the large number of unknowns and voltage measurements in the stroke-monitoring setup. In this paper, the search direction is solved in a less memory-intensive way using a prior-conditioned LSQR (MLSQR) iteration.18 This iterative method significantly reduces the required number of iterations by limiting its solution space based on the regularization used in the inverse problem. The formulation for prior conditioning also produces an alternative way for obtaining the gradients and Hessians of the regularization terms that are required for the lagged Gauss–Newton iteration (for more details, see Arridge et al.18 and Harhanen et al.19).

To obtain a good initialization for the optimization of Eq. (24), we utilize an anatomically guided initial estimate where the domain is divided to three subvolumes (scalp, skull, and brain) based on the reference image κ(x) and a three parameter nonlinear least squares estimation is carried out to obtain an initial conductivity with a constant conductivity value for each of the subvolumes.

3.

Methods

3.1.

Computational Head Model

The highly accurate head model used for simulations of the EIT measurements was published by Paldanius et al.20 and is based on a model from the population head model repository.31 The original head model consisted of surface mesh presentations of the scalp, skull, CSF, white matter, gray matter, and cerebellum. For 3D meshing, these surface meshes were imported into ScanIP Simpleware. 32 circular electrodes with a 10 mm diameter were assigned based on the modified 10-5 EEG electrode placement system designed for EIT measurements [Fig. 1(a)] by Goren et al.5 To simulate a hemorrhage growing over time, concentric spheres with diameters from 10 to 30 mm in 5 mm increments were manually placed in the volume of interest in the brain parenchyma. The 30 mm sphere has a volume of 14.14 ml, matching the 14 ml median size of ICH.32 For the meshing, the target maximum error from the original surface mesh was set to 0.15 mm to ensure preservation of thin details, such as the layer of CSF around the brain. The maximum element edge length on the electrodes was set to 1 mm to make the mesh denser in the regions near the electrodes as is required for sufficient numerical accuracy in EIT simulations. These settings resulted in a 3D mesh of 2.5 million tetrahedral elements, which was then exported to COMSOL for simulation of the measurement data.

Fig. 1

(a) Electrode locations chosen from the 10-5 system as described by Ref. 5. (b) The hemorrhage location close to the surface of the brain. (c) The hemorrhage location deep in the basal ganglia region.

JMI_11_1_014502_f001.png

Two locations of simulated hemorrhage were considered, one cortical hemorrhage close to the surface of the brain and one deep in the basal ganglia region of the brain. This provided varying level of challenge for the EIT algorithms, as the sensitivity of the EIT measurement is lower when the perturbation is deeper in the brain parenchyma.20 The basal ganglia region is also the most common location of ICH, as perforative small arteries are prone to rupture in patients with high blood pressure. The chosen hemorrhage locations are shown in Figs. 1(b) and 1(c).

Model setup for the simulation was performed in COMSOL 5.6. The tissue layers were assigned with corresponding dielectric material properties from Gabriel et al.33 (Table 1). The properties did not include values for the scalp as it is composed of multiple different tissues. Hence, values of muscle tissue which were within the range of values reported for the scalp in other literature34 were used. The dielectric values for cortical bone were used for the skull as it matches the recommended value for single-layer skull models.34 The material properties for 1 kHz frequency were used, as the measurements were simulated at that frequency.

Table 1

Material properties used in the simulation at 1 kHz frequency.33

TissueConductivity (S/m)Permittivity (F/m)
Scalp0.32434,932
Skull0.022702
CSF2.00109
White matter0.0669,810
Gray matter0.10164,062
Cerebellum0.12164,358
Blood0.705259

As COMSOL does not natively support the CEM, the CEM implementation tailor-made for COMSOL from Fouchard et al.35 was used. An iterative BiCGStab solver was selected for solving the FEM system.

The growth of the hemorrhage was simulated by changing the material properties of the concentric spheres representing the hemorrhage from white matter and gray matter to blood. In the first simulation, all the spheres were assigned material properties of white matter for the hemorrhage close to surface of the brain and gray matter for the hemorrhage deep in the basal ganglia region to represent a healthy brain. In the following simulations, the spheres were sequentially assigned material properties of blood with the hemorrhage finally being 30 mm in diameter in the last simulation.

3.2.

Simulated Measurement Setup

In the simulation, a total of 32 independent pairwise current injections at 1 kHz and 1 mA magnitude were used. The injection pattern used was 1-13, 2-14,…, 32-12, where the first number is the number of the current injecting electrode and the second number is the current sink electrode. In this current injection pattern, there is sufficient distance between the current injecting and current sink electrodes for some of the current to pass through the skull and the brain, instead of using adjacent electrodes resulting in current being mainly shunted along the well conducting scalp. For each injection, the noninjecting electrodes recorded the resulting potentials from the scalp. The simulated measurement data were saved as differential potentials between electrodes 1-2, 2-3,…, 32-1 and exported in .csv format.

Noisy realizations of the simulated measurement data were obtained by adding Gaussian zero mean random noise with a standard deviation of 1.84×105  V to the simulated noise-free measurements. The standard deviation of the noise was 0.067% of the maximum amplitude of the EIT data from a healthy brain and corresponds to the approximated relative noise level of the prototype stroke measurement device from Toivanen et al.15

3.3.

Computational Meshes and Regularization Parameters for the Imaging Algorithms

All imaging algorithms used a MATLAB-based solver for the FEM model [Eq. (6)] where the electric potential u(x) was approximated in a tetrahedral mesh of 116,235 nodes and 597,631 elements with refinement near the electrodes as is required for sufficient accuracy for the numerical solving of the model [Eqs. (1)–(5)]. Such refinement is not needed for the discretization of the conductivity and would unnecessarily increase the number of unknown conductivity values. Thus the conductivity was approximated in a coarser and uniform tetrahedral mesh of 38,433 nodes and 207,453 elements, leading to N=38,433 unknown conductivity values. These meshes used for the imaging algorithms are shown in Fig. 2 and they have a considerably smaller number of elements than the highly accurate, 2.5 million element mesh used in COMSOL for measurement data simulation from the anatomically detailed head model.

Fig. 2

Computational meshes used for the imaging algorithms. The computational mesh for (a) conductivity and (b) the electric potential with electrodes highlighted.

JMI_11_1_014502_f002.png

The regularization parameters for each algorithm were tuned manually to give the best visual quality of the reconstructions. For the TV regularized absolute imaging algorithm, the values α=0.01 and β=0.001 were used in Eq. (10). For the linear difference imaging algorithm, a standard deviation of conductivity std(σ)=2σ0 was used in Eq. (16) and the parameter a was calculated by setting a covariance of 1% of the variance at a correlation distance d=lΩ/4, where lΩ=20.4cm is the length of the domain Ω from the back of the head to the forehead. For the stroke-monitoring algorithm, the values αδσ=0.005 and β=0.001 were used in Eq. (26), the values ασ1=107 and β=0.001 were used in Eq. (27), and cross sections of the reference image κ(x), the indicator function γ(x) and the region of interest are shown in Fig. 3. The reference image was obtained by approximately mapping the true skull volume into the computational mesh used for discretization of the conductivity, in a similar manner as a CT image could be used in a clinical setup. The indicator function γ(x) corresponds to the boundaries in the reference image and the region of interest was chosen to approximately correspond to the brain volume. This would be a safe choice also when monitoring for secondary hemorrhages, where the hemorrhage can sometimes occur also outside of the brain volume weakened by the ischemic stroke.

Fig. 3

(a) Cross sections of the reference image κ(x), (b) the indicator function γ(x), and (c) the region of interest used for the stroke-monitoring algorithm.

JMI_11_1_014502_f003.png

4.

Results and Discussion

The simulated measurement data from all hemorrhage-growth and no-growth scenarios were used to compute estimates of δσ with all three algorithms, TV, LD, and MO. These estimates were computed for both hemorrhage locations, resulting in a total of 126 estimates of δσ. Results for a single hemorrhage progression chain (growth from 15 to 20 to 25 to 30 mm hemorrhage diameter, corresponding to volume increases of 2.42, 3.99, and 5.96 ml) and for a single no-growth case (20 to 20 mm hemorrhage diameter, corresponding to a volume change of 0 ml) are shown in Fig. 4 for the hemorrhage close to the surface of the brain and in Fig. 5 for the hemorrhage deep in the basal ganglia region.

Fig. 4

Slices of the computational head model at selected measurement times t corresponding to hemorrhage diameters Ø and hemorrhage volumes V, and estimates of the conductivity change δσ=σi+1σi between two consecutive measurement times ti and ti+1 for the hemorrhage close to the surface of the brain with the TV regularized absolute imaging algorithm (TV), the linear difference imaging algorithm (LD), and the stroke-monitoring algorithm (MO). (a) Single cross sections through the center of the simulated hemorrhage and (b) transparent 3D plots.

JMI_11_1_014502_f004.png

Fig. 5

Slices of the computational head model at selected measurement times t corresponding to hemorrhage diameters Ø and hemorrhage volumes V, and estimates of the conductivity change δσ=σi+1σi between two consecutive measurement times ti and ti+1 for the hemorrhage deep in the basal ganglia region with the TV regularized absolute imaging algorithm (TV), the linear difference imaging algorithm (LD), and the stroke-monitoring algorithm (MO). (a) Single cross sections through the center of the simulated hemorrhage and (b) transparent 3D plots.

JMI_11_1_014502_f005.png

For the cortical hemorrhage close to the surface of the brain, the cross sections in Fig. 4(a) show that all three algorithms are capable of indicating the largest and second largest volume change (fourth and third rows of estimates), but the smallest volume change (second row of estimates) is indicated only by the stroke-monitoring algorithm. Furthermore, an examination of the transparent 3D plots in Fig. 4(b) reveals that the estimates of the second largest volume change (third row of estimates) with the reference algorithms contain additional conductive inclusions that could be misinterpreted as additional hemorrhages. Overall, the estimates produced by the stroke-monitoring algorithm show a much more regular and better-localized conductivity change compared to the estimates of the reference algorithms.

The more challenging hemorrhage location deep in the basal ganglia region causes deterioration of estimate quality for all algorithms as is clearly visible in the cross sections in Fig. 5(a). For the reference algorithms, even the cross sections that seem to indicate the expansion somewhat correctly (second row of estimates), they show just one of many conductive artefacts, as revealed by a closer examination of the transparent 3D plots in Fig. 5(b). Also, the estimates from the stroke-monitoring algorithm are less regular and localized, and now more of the estimates of the conductivity change contain negative values, both indicative of less reliable results. Nevertheless, out of the three algorithms, the estimates of the stroke-monitoring algorithm best indicate the expansion of the hemorrhage also in this more challenging monitoring setup.

The results for all 126 estimation cases are shown as heat maps in Fig. 6. The first row shows for reference the true volume changes and the second row shows the norms of the noise-free difference data for the two hemorrhage locations. The third row shows the norms of the noisy difference data that show the magnitude of the change in the noisy measurement data from the hemorrhage growth, with the no-growth values showing only the level of the measurement noise. For comparison of the estimates, the integrals IDi,Dj=Ω(δσ)Di,Djdx were computed for all inclusion diameters Di and Dj and then normalized for each algorithm and simulated hemorrhage location separately. These were used to compute the adjusted normalized integrals shown in the heat maps

Eq. (31)

QDi,Dj=IDi,Djmax(abs(IDk,Dk)),
where IDk,Dk are the normalized integrals of the no-growth scenarios, and the maximum of their absolute values is here regarded as a rough approximation for the maximum deviation in IDi,Dj caused by measurement noise. Removing IDk,Dk adjusts the normalized integrals so that all resulting QDi,Dj<0 are expected to correspond to estimates that are comparable to those caused by noise alone and are thus not expected to be useful. This is further highlighted by the chosen color scheme where all QDi,Dj<0 are shown in white. In contrast, the closer the values are to 1, the higher the expected detectability of hemorrhage growth is.

Fig. 6

Comparison of all 126 cases. Diameter 1 and 2 are the diameters of the simulated hemorrhage at the time of the first and the second measurements, respectively. (a) First row: the true volume changes in ml; second row: the norms of the noise-free difference data for the two hemorrhage locations in mV; third row: norms of the noisy difference data in mV. (b) The adjusted normalized integrals QDi,Dj for the TV regularized absolute imaging algorithm (TV), the linear difference imaging algorithm (LD), and the stroke-monitoring algorithm (MO).

JMI_11_1_014502_f006.png

Overall, the heat maps of the stroke-monitoring algorithm have larger values than either of the reference algorithms, indicating better detectability of hemorrhage growth. The difference is clearest when comparing the second or the third columns in the heat maps. The values highlighted with black borders in the heat maps correspond to the estimates shown in Figs. 4 and 5. The highlighted norms of the noisy difference data show that for the hemorrhage close to the surface of the brain the two larger conductivity changes (20 to 25 and 25 to 30 mm) are clearly different from the no-growth cases, but the value corresponding to the smallest conductivity change (15 to 20 mm) is already very close to the no-growth cases. This small conductivity change was visible only in the estimate from the monitoring algorithm and this is at least partly because the algorithm benefits from utilizing both measurements and their connection, unlike the linear difference algorithm that uses difference data or the TV regularized absolute imaging algorithm that treats the measurements as completely separate.

5.

Conclusion

In this paper, we used simulated data from an anatomically detailed computational head model to study the feasibility of EIT for continuous bedside monitoring of ICHs and detection of secondary hemorrhages and got promising results using a tailor-made, patient-specific stroke-monitoring algorithm. The monitoring algorithm was referenced against a TV-regularized absolute imaging algorithm and a linear difference imaging algorithm, and it was shown that the proposed stroke-monitoring algorithm is capable of indicating smaller changes in the simulated hemorrhages than either of the reference algorithms. In the simulation tests of this paper, the smallest volume change in the simulated hemorrhage detected by the proposed stroke-monitoring algorithm was 2.42 ml for the cortical hemorrhage close to the surface of the brain and 3.99 ml for the hemorrhage deep in the basal ganglia region.

Further testing of the feasibility of EIT for monitoring of ICHs and detection of secondary hemorrhages is envisioned to include more complex simulation studies, laboratory phantom measurements that build upon the ones presented by Toivanen et al.,15 live animal measurements, and eventually measurements from human stroke patients. Future development presents a number of challenges, including the uncertainty of the electrode locations, electrode impedance variations, and effects of computational domain truncation. However, our results provide first simulation-based proof of concept, and we believe that EIT-based bedside stroke monitoring becomes a valuable tool in follow-up of ICH patients for an early warning of changes in the brain tissue of acute stroke patients that would prompt an immediate need for control CT scanning. In the most severe cases, this would lead to timely decision to advance to intensive care and to neurosurgical intervention.

Disclosures

Authors have no conflicts of interest to disclose.

Code and Data Availability

The codes and the data utilized in this study are not publicly available. The codes and the data may be requested from the authors, and permission for access will be granted or denied on a case-by-case basis.

Acknowledgments

This work was supported in part by the Academy of Finland (Project No. 336791, Finnish Centre of Excellence in Inverse Modeling and Imaging), the Jane and Aatos Erkko Foundation and Neurocenter Finland. Model geometry data were provided in part by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. The authors wish to acknowledge CSC – IT Center for Science, Finland, for computational resources.

References

1. 

M. Katan and A. Luft, “Global burden of stroke,” Semin. Neurol., 38 (02), 208 –211 https://doi.org/10.1055/s-0038-1649503 SEMNEP 0271-8235 (2018). Google Scholar

2. 

S. M. Greenberg et al., “2022 Guideline for the management of patients with spontaneous intracerebral hemorrhage: a guideline from the American Heart Association/American Stroke Association,” Stroke, 53 (7), e282 –e361 https://doi.org/10.1161/STR.0000000000000407 SJCCA7 0039-2499 (2022). Google Scholar

3. 

F. Schlunk and S. M. Greenberg, “The pathophysiology of intracerebral hemorrhage formation and expansion,” Transl. Stroke Res., 6 (4), 257 –263 https://doi.org/10.1007/s12975-015-0410-1 (2015). Google Scholar

4. 

E. Malone et al., “Stroke type differentiation using spectrally constrained multifrequency EIT: evaluation of feasibility in a realistic head model,” Physiol. Meas., 35 1051 –1066 https://doi.org/10.1088/0967-3334/35/6/1051 PMEAE3 0967-3334 (2014). Google Scholar

5. 

N. Goren et al., “Multi-frequency electrical impedance tomography and neuroimaging data in stroke patients,” Sci. Data, 5 (1), 180112 https://doi.org/10.1038/sdata.2018.112 (2018). Google Scholar

6. 

L. Horesh, Some Novel Approaches in Modelling and Image Reconstruction for Multi Frequency Electrical Impedance Tomography of the Human Brain, University of London( (2006). Google Scholar

7. 

V. Candiani et al., “Approximation error method for imaging the human head by electrical impedance tomography,” Inverse Prob., 37 (12), 125008 https://doi.org/10.1088/1361-6420/ac346a INPEEY 0266-5611 (2021). Google Scholar

8. 

J. P. Agnelli et al., “Classification of stroke using neural networks in electrical impedance tomography,” Inverse Prob., 36 (11), 115008 https://doi.org/10.1088/1361-6420/abbdcd INPEEY 0266-5611 (2020). Google Scholar

9. 

B. McDermott et al., “Multi-frequency symmetry difference electrical impedance tomography with machine learning for human stroke diagnosis,” Physiol. Meas., 41 (7), 075010 https://doi.org/10.1088/1361-6579/ab9e54 PMEAE3 0967-3334 (2020). Google Scholar

10. 

V. Candiani and M. Santacesaria, “Neural networks for classification of stroke in electrical impedance tomography on a 3D head model,” Math. Eng., 4 (4), 1 –22 https://doi.org/10.3934/mine.2022029 (2022). Google Scholar

11. 

Y. Shi et al., “Sparse image reconstruction of intracerebral hemorrhage with electrical impedance tomography,” J. Med. Imaging, 8 (1), 014501 https://doi.org/10.1117/1.JMI.8.1.014501 JMEIET 0920-5497 (2021). Google Scholar

12. 

C.-H. Xu et al., “Real-time imaging and detection of intracranial haemorrhage by electrical impedance tomography in a piglet model,” J. Int. Med. Res., 38 (5), 1596 –1604 https://doi.org/10.1177/147323001003800504 (2010). Google Scholar

13. 

M. Dai et al., “In vivo imaging of twist drill drainage for subdural hematoma: a clinical feasibility study on electrical impedance tomography for measuring intracranial bleeding in humans,” PloS One, 8 (1), e55020 https://doi.org/10.1371/journal.pone.0055020 POLNCL 1932-6203 (2013). Google Scholar

14. 

B. Yang et al., “Comparison of electrical impedance tomography and intracranial pressure during dehydration treatment of cerebral edema,” Neuroimage Clin., 23 101909 https://doi.org/10.1016/j.nicl.2019.101909 (2019). Google Scholar

15. 

J. Toivanen et al., “Monitoring hemorrhagic strokes using EIT,” Bioimpedance and Spectroscopy, 271 –298 Elsevier( (2021). Google Scholar

16. 

D. Liu et al., “Estimation of conductivity changes in a region of interest with electrical impedance tomography,” Inverse Prob. Imaging, 9 (1), 211 –229 https://doi.org/10.3934/ipi.2015.9.211 (2015). Google Scholar

17. 

V. Kolehmainen, M. J. Ehrhardt and S. R. Arridge, “Incorporating structural prior information and sparsity into EIT using parallel level sets,” Inverse Prob. Imaging, 13 (2), 285 –307 https://doi.org/10.3934/ipi.2019015 (2019). Google Scholar

18. 

S. Arridge, M. Betcke and L. Harhanen, “Iterated preconditioned LSQR method for inverse problems on unstructured grids,” Inverse Prob., 30 (7), 075009 https://doi.org/10.1088/0266-5611/30/7/075009 INPEEY 0266-5611 (2014). Google Scholar

19. 

L. Harhanen et al., “Edge-enhancing reconstruction algorithm for three-dimensional electrical impedance tomography,” SIAM J. Sci. Comput., 37 (1), B60 –B78 https://doi.org/10.1137/140971750 SJOCE3 1064-8275 (2015). Google Scholar

20. 

A. Paldanius et al., “Sensitivity analysis highlights the importance of accurate head models for electrical impedance tomography monitoring of intracerebral hemorrhagic stroke,” IEEE Trans. Biomed. Eng., 69 (4), 1491 –1501 https://doi.org/10.1109/TBME.2021.3120929 IEBEAX 0018-9294 (2021). Google Scholar

21. 

K. Cheng et al., “Electrode models for electric current computed tomography,” IEEE Trans. Biomed. Eng., 36 918 –924 https://doi.org/10.1109/10.35300 IEBEAX 0018-9294 (1989). Google Scholar

22. 

E. Somersalo, M. Cheney and D. Isaacson, “Existence and uniqueness for electrode models for electric current computed tomography,” SIAM J. Appl. Math., 52 (4), 1023 –1040 https://doi.org/10.1137/0152060 SMJMAP 0036-1399 (1992). Google Scholar

23. 

M. Vauhkonen et al., “Tikhonov regularization and prior information in electrical impedance tomography,” IEEE Trans. Med. Imaging, 17 (2), 285 –293 https://doi.org/10.1109/42.700740 ITMID4 0278-0062 (1998). Google Scholar

24. 

J. Kaipio et al., “Statistical inversion and Monte Carlo sampling methods in electrical impedance tomography,” Inverse Prob., 16 (5), 1487 –1522 https://doi.org/10.1088/0266-5611/16/5/321 INPEEY 0266-5611 (2000). Google Scholar

25. 

L. Rudin, S. Osher and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D, 60 (1–4), 259 –268 https://doi.org/10.1016/0167-2789(92)90242-F (1992). Google Scholar

26. 

D. C. Barber and A. D. Seagar, “Fast reconstruction of resistance images,” Clin. Phys. Physiol. Meas., 8 (4A), 47 https://doi.org/10.1088/0143-0815/8/4A/006 CPPMD5 0143-0815 (1987). Google Scholar

27. 

A. P. Bagshaw et al., “Electrical impedance tomography of human brain function using reconstruction algorithms based on the finite element method,” NeuroImage, 20 (2), 752 –764 https://doi.org/10.1016/S1053-8119(03)00301-X NEIMEF 1053-8119 (2003). Google Scholar

28. 

C. Lieberman, K. Willcox and O. Ghattas, “Parameter and state model reduction for large-scale statistical inverse problems,” SIAM J. Sci. Comput., 32 (5), 2523 –2542 https://doi.org/10.1137/090775622 SJOCE3 1064-8275 (2010). Google Scholar

29. 

D. Liu et al., “A nonlinear approach to difference imaging in EIT; assessment of the robustness in the presence of modelling errors,” Inverse Prob., 31 (3), 035012 https://doi.org/10.1088/0266-5611/31/3/035012 INPEEY 0266-5611 (2015). Google Scholar

30. 

D. Liu et al., “Nonlinear difference imaging approach to three-dimensional electrical impedance tomography in the presence of geometric modeling errors,” IEEE Trans. Biomed. Eng., 63 (9), 1956 –1965 https://doi.org/10.1109/TBME.2015.2509508 IEBEAX 0018-9294 (2016). Google Scholar

31. 

E. Lee et al., “Investigational effect of brain-scalp distance on the efficacy of transcranial magnetic stimulation treatment in depression,” IEEE Trans. Magn., 52 (7), 1 –4 https://doi.org/10.1109/TMAG.2015.2514158 IEMGAQ 0018-9464 (2016). Google Scholar

32. 

D. Robinson et al., “What is the median volume of intracerebral hemorrhage and is it changing?,” Int. J. Stroke, 17 (5), 576 –582 https://doi.org/10.1177/17474930211032594 (2022). Google Scholar

33. 

S. Gabriel, R. Lau and C. Gabriel, “The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues,” Phys. Med. Biol., 41 (11), 2271 https://doi.org/10.1088/0031-9155/41/11/003 PHMBA7 0031-9155 (1996). Google Scholar

34. 

H. McCann, G. Pisano and L. Beltrachini, “Variation in reported human head tissue electrical conductivity values,” Brain Topogr., 32 (5), 825 –858 https://doi.org/10.1007/s10548-019-00710-2 BRTOEZ 0896-0267 (2019). Google Scholar

35. 

A. Fouchard et al., “Flexible numerical platform for electrical impedance tomography,” in COMSOL Conf., (2015). Google Scholar

Biographies of the authors are not available.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 International License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Jussi Toivanen, Antti Paldanius, Bachir Dekdouk, Valentina Candiani, Asko Hänninen, Tuomo Savolainen, Daniel Strbian, Nina Forss, Nuutti Hyvönen, Jari Hyttinen, and Ville Kolehmainen "Simulation-based feasibility study of monitoring of intracerebral hemorrhages and detection of secondary hemorrhages using electrical impedance tomography," Journal of Medical Imaging 11(1), 014502 (30 January 2024). https://doi.org/10.1117/1.JMI.11.1.014502
Received: 11 September 2023; Accepted: 8 January 2024; Published: 30 January 2024
Advertisement
Advertisement
KEYWORDS
Brain

Electrodes

Head

Computer simulations

Detection and tracking algorithms

Biological imaging

Electrical conductivity

Back to Top