Electrode positioning to investigate the changes of the thoracic bioimpedance caused by aortic dissection – a simulation study

The volumetric expansion of the blood-filled aorta changes corresponding to the cardiac pulse wave. For the sake of simplicity, a spatial average time-dependent cross-sectional radius of the aorta has been used in the simulation model for two sections separately, the aortic arch and the descending aorta, see figure 2. The data are based on measurements provided in [11] from a young, healthy male volunteer at rest.

Fig. 2

The electrical properties of resting blood mainly depend on the volume fraction of red blood cells (RBCs) called haematocrit, the temperature, and cell shape. However, the electrical properties of flowing blood are found to be influenced by the flow rate [12].

A spatial average time-dependent velocity of the blood flowing inside the aorta, taken from the experimental data provided in [11], has been taken into account for the aortic arch and the descending aorta (figure 3). During the systolic phase of a cardiac cycle, the heart contracts to pump blood into the aorta, and in the diastolic phase, the heart relaxes after contraction. This pulsatile blood flow causes the variation of blood conductivity inside the aorta. The reason is the orientation and deformation of the RBCs in case of flowing blood. At higher velocities, the shear stress increases, which consequently deforms the RBCs in the layer with the highest stress close to the vessel wall and also aligns them throughout the vessel. Both effects lead to a higher conductivity than the resting blood (figure 4) [13,14].

Fig. 3

Fig. 4

The Maxwell–Fricke equation for the conductivity of blood reads [13]: (1) σbloodσplasma=1−H1+(C−1)H {{{\sigma _{blood}}} \over {{\sigma _{plasma}}}} = {{1 - H} \over {1 + (C - 1)H}} where σ_blood and σ_plasma are the conductivities of blood and plasma, respectively, H is the haematocrit expressed as the volume fraction of RBCs relative to the total blood volume, and C is a factor that depends on the geometry and orientation of the RBCs. Based on the formulation described in [13] (blood is considered as a Newtonian fluid with a steady flow), the blood conductivity changes as a function of reduced average velocity 〈v/R〉 (v is the spatial average blood velocity, and R is the spatial average cross-sectional radius of the aorta), for different haematocrit levels are shown in figure 5. From this figure, it is evident that more conductivity changes of blood exist with higher haematocrit levels. It is also evident from this figure that the conductivity changes mostly occur in the lower blood velocities and for higher velocities, the slope of the Δσ curve decreases significantly.

Fig. 5

In [15,16], it is shown that the blood conductivity during pulsatile blood flow is not the same at any given velocity during acceleration and deceleration. This disparity is a consequence of the RBCs inability to achieve complete randomization at end-systole, which leads to less but still considerable conductivity changes during the cardiac cycle. However, for simplicity, conductivity changes which are shown in figure 5 have been assumed in the simulation model.

A 3D numerical simulation model is used to investigate the changes in the electric potential and the impedance changes on the thorax surface. The model has been set up in COMSOL Multiphysics [17] for the underlying time-harmonic current flow problem. Since the duration of the cardiac cycle is much higher than the period of the injecting current, simulations can be performed in the frequency domain. The electric potential drop is evaluated between the measuring electrodes by solving the Laplace equation for the electric potential V: (2) ∇⋅([σ+jωε]∇V)=0. \nabla \cdot ([\sigma + j\omega \varepsilon ]\nabla V) = 0.

The model consists of a simplified geometry, as shown in figure 6. Three pairs of source (injection) electrodes are placed on the surface of the thorax (each pair in one vertical line) and inject an alternating current with a magnitude of 5 mA and a frequency of 100 kHz asynchronously. For each injection, the electric potential drop is evaluated between five measurement electrode pairs (each pair in one vertical line) which leads to the thoracic impedance (3) Z_=V_top−V_bottomI_. \underline Z = {{{{\underline V}_{top}} - {{\underline V}_{bottom}}} \over {\underline I}}.

Fig. 6

The boundary conditions are:

V = constrant, on the top source electrode;

It has been shown in the literature [21, 22] that the blood flow is highly disturbed locally inside the aorta and changes to turbulent flow with strong recirculation. As depicted in figure 7, flow disturbances occur around the dissection, which inhibits the deformation and orientation of the RBCs. Thus, the flow shear rate and, consequently, the electrical properties of blood are altered. At the highest blood flow velocity and consequently the highest deformation and orientation rate of the RBCs, a remarkable difference in the electrical conductivity between the healthy (non-disturbed flow) and the aortic dissection conditions can be expected [6]. Since no experimental or simulation data exist regarding conductivity changes of blood in this kind of disturbed flow, it is assumed that with a radially growing false lumen also the blood flow disturbances increase and the conductivity changes decrease. To quantify this assumption, a damage factor DF has been introduced, which is the ratio of the volume of the dissection to the maximum volume of the false lumen: (4) DF=π⋅RFL⋅2hFLπ⋅max(RFL2)⋅hFL=RFL2max(RFL2) DF = {{\pi \cdot R_{FL \cdot}^2{h_{FL}}} \over {\pi \cdot \max \left({R_{FL}^2} \right) \cdot {h_{FL}}}} = {{R_{FL}^2} \over {\max \left({R_{FL}^2} \right)}} where R_FL and h_FL are the radius and the height of the false lumen, respectively. The maximum value for the radius of the false lumen in the simulation model is 1.5 cm. At this level, it is assumed that the flow disturbances in the descending aorta are at maximum (figure 8).

Fig. 7

Fig. 8

The damage factor DF is applied to the conductivity changes of blood in the descending aorta during the cardiac cycle to model the decrease in the conductivity changes of blood due to dissection (see the section on Implementation of global sensitivity analysis …, on the next page).

The aim of a Global Sensitivity Analysis (GSA) is to quantify the connection between the variance of the model output given the variability of its input. GSA is distinguished from a local sensitivity analysis since it investigates the whole input space of each random variable. This study aims to use a GSA technique to identify which electrode configuration has significant changes in the impedance cardiogram −|dZ/dt| given by the uncertainty on the developed status of AD. The discrepancies among obtained impedance cardiograms of the healthy and dissected models determine different states of the patient and are the criteria to simulate AD identification.

One of the most known techniques in GSA is the variance-based method. Here, the output variance is portioned in the sum of the contributions of each random variable. Consider a mathematical model output y as a function of an input random vector x of dimension n, the Sobol’ decomposition reads as [22] (5) y=f(x1,…xn)=f0+∑i=1nfi(xi)+∑1≤i≤j≤nfij(xi,xj)+⋯+f1,…,n(x1,….xn) \matrix{{y = f({x_1}, \ldots {x_n})} \hfill & {= {f_0} + \sum\limits_{i = 1}^n {{f_i}({x_i})}} \hfill \cr {} \hfill & {+ \sum\limits_{1 \le i \le j \le n} {{f_{ij}}({x_i},{x_j}) + \cdots}} \hfill \cr {} \hfill & {+ {f_{1, \ldots,n}}({x_1}, \ldots.{x_n})} \hfill \cr} where f₀ is a constant. The integral of each summand of (5) is equal to zero, which ensures the uniqueness of the decomposition given the integrity of the model over the random input space, i.e. ∫ f(x_i) f(x_j)dx_idx_j = 0. The terms of (5) can be then derived univocally from the conditional expectation E of the model output [23] (6) f0=E[y] {f_0} = E[y] (7) fi=E[y|xi]−E[y] {f_i} = E[y|{x_i}] - E[y] (8) fij=E[y|xi,xj]−fi−fj−E[y] {f_{ij}} = E[y|{x_i},{x_j}] - {f_i} - {f_j} - E[y]

A random variable x_i is considered to be influential (non-influential) to the model output if the conditional variance V(E[y|x_i]) is large (small) enough compared to the variance of the quantity of interest. Therefore, the first-order sensitivity index (or main Sobol’ index) can be derived as: (9) Si=V(E[y|xi])V(y) {S_i} = {{V(E[y|{x_i}])} \over {V(y)}} which represents the contribution of the random variable x_i to the change of model output, without taking into consideration the effect of its interaction with other input variables. Consequently, the total-order sensitivity index (or total Sobol’ index), which evaluates the total effect of such input parameter, has to account for the conditional variance of the output, but conditioning to all factors except one, x_~i: (10) STi=E[V(y|x∼i)]V(y)=1-V[E(y|x∼i)]V(y). {S_{Ti}} = {{E[V(y|{x_{\sim i}})]} \over {V(y)}} = 1 - {{V[E(y|{x_{\sim i}})]} \over {V(y)}}.

In this application, only the first-order and the total-order indices are considered. Any interaction between the input random variable can be derived by subtracting the first index to the total. Consequently, the difference will result in the amount of interaction present in the model.

The Sobol’ indices are computed from a Polynomial Chaos Expansion (PCE) of the model [24], which also represents a valid mathematical metamodel. PCE consists of the sum of orthogonal, multivariate polynomials ψ_α of increasing order up to some maximal polynomial order p. The polynomials are multiplied by expansion coefficients y_α, which can be estimated with different methods [25]. A PCE is written as: (11) y(x)≈fPCE(x)=∑α∈Ayαψα(x) y({\bf{\it x}}) \approx {f_{PCE}}({\bf{\it x}}) = \sum\limits_{\alpha \in A} {{y_\alpha}{\psi _\alpha}({\bf{\it x}})} where α is a multi-index that refers to the degree of each polynomial and each input parameter, the multivariate polynomials ψ_α are defined as the product of univariate polynomials of order α_i, i.e. ψ_αi. The univariate polynomials are generated following the Askey scheme for the composition of polynomials [26]. Finally, from the PCE, it is possible to estimate the two sensitivity indices as the ratio between the PCE coefficients. Since the case study is evaluated in time, the implementation of time-dependent indices is implemented following [27] as: (12) Si≈∑α∈Aiyα2(t)∑α∈A;α≠1yα2(t)=∫0t∑α∈Aiyα2(t′)dt′∫0t∑α∈A;α≠1yα2(t′)dt′ {S_i} \approx {{\sum\nolimits_{\alpha \in {A_i}} {y_\alpha ^2(t)}} \over {\sum\nolimits_{\alpha \in A;\alpha \ne 1} {y_\alpha ^2(t)}}} = {{\int_0^t {\sum\nolimits_{\alpha \in {A_i}} {y_\alpha ^2(t')dt'}}} \over {\int_0^t {\sum\nolimits_{\alpha \in A;\alpha \ne 1} {y_\alpha ^2(t')dt'}}}} (13) STi≈∑α∈AiTyα2(t)∑α∈A;α≠1yα2(t)=∫0t∑α∈AiTyα2(t′)dt′∫0t∑α∈A;α≠1yα2(t′)dt′. {S_{Ti}} \approx {{\sum\nolimits_{\alpha \in A_i^T} {y_\alpha ^2(t)}} \over {\sum\nolimits_{\alpha \in A;\alpha \ne 1} {y_\alpha ^2(t)}}} = {{\int_0^t {\sum\nolimits_{\alpha \in A_i^T} {y_\alpha ^2(t')dt'}}} \over {\int_0^t {\sum\nolimits_{\alpha \in A;\alpha \ne 1} {y_\alpha ^2(t')dt'}}}}.

The PCE is computed through the UQLab toolbox for Matlab [28,29]. However, the time-dependent indices are developed manually from the extrapolation of the PCE coefficients.

To assess the sensitivity analysis, the input and output spaces together with the numerical models, have to be established. Since the aim of the study is to catch the difference between different health conditions, two numerical models are set. The first one refers to the healthy condition and the second one to the dissected condition. As described in section 2–2, the latter differs from the first one in the presence of the false lumen and another blood flow profile. Therefore, different input spaces are produced for each model. Besides, introducing variability in the input space of the models will guarantee the realization of as many patient-specific cases as the dimension of the input sample. Thus, a deeper understanding of the impedance cardiography for a human thorax can be revealed. From the models’ evaluations, the PCEs for the healthy and dissected conditions are constructed, and analyzing the differences between them will guide toward the choice of the best electrode configuration.

The input space of the healthy case is composed of only two random variables, namely the maximum radius of the true lumen R_TL and the blood conductivity coefficient θ_H. R_TL is considered uniformly distributed between 1.35 and 1.95 cm, according to the study [30]. As it was shown in figure 2, the average radius of the aorta changes in time due to pressure changes over a cardiac cycle. Since the aorta in the simulation model is considered as a blood-filled lumen, different values of the R_TL emulates different blood volumes dilating the aorta, in other words, different stroke volume values. It has to be mentioned that changes in the stroke volume will also vary the peak velocity of the blood passing the aorta. Since for higher velocities, almost all the RBCs are entirely aligned and deformed, the differences among the blood conductivity changes for different stroke volumes are not significant. Therefore, here, for simplicity, it is neglected. Based on the distribution's moments of R_TL, the stroke volume changes approximately between 62 and 140 ml in the simulation model.

The blood conductivity changes as a function of the reduced average velocity 〈v/R〉 for five different haematocrit levels are shown in figure 5. For the parameterization of the haematocrit-dependent blood conductivity changes corresponding to other haematocrit values, the coefficient θ_H has been introduced. It emulates a scaling factor for Δσ at a haematocrit level of 35% up to 55% at maximum. Given h as the index that represents different haematocrit levels, the parameter θ_H,h is set as: (14) θH,h={1if h=11+ahif h=2,…,5 {\theta _{H,h}} = \left\{{\matrix{1 \hfill & {if\,h = 1} \hfill \cr {1 + {a_h}} \hfill & {if\,h = 2, \ldots,5} \hfill \cr}} \right. where h = 1 to h = 5 refer to H = 35%. H = 40%, H = 45%, H = 50% and H = 55% respectively.

ah=∑k=2h(ΔσH,k−ΔσH,k−1)/ΔσH,k−1 {a_h} = \sum {_{k = 2}^h} \left({\Delta {\sigma _{H,k}} - \Delta {\sigma _{H,k - 1}}} \right)/\Delta {\sigma _{H,k - 1}} computes the difference in conductivity changes between two successive haematocrit levels, Δσ_H,k are the conductivity changes of the blood for a given haematocrit level H and index k. In conclusion, since the variable H is uniformly distributed between 35% and 55%, the parameter θ_H,h is also uniformly distributed between 1 and 1.1. Therefore, the blood conductivity changes for the healthy case are computed as: (15) Δσ(v(t),H)=θH,h ⋅ Δσ(v(t))|H=35% \Delta \sigma (v(t),H) = {\theta _{H,h\, \cdot \,}}\Delta \sigma (v(t)){|_{H = 35\%}} where v(t) is the velocity of blood and Δσ(v(t))|_H=35% represents the velocity-based conductivity changes for H = 35%, as shown in figure 5.

The dissected condition includes both the random variables of the healthy condition plus the radius of the false lumen R_FL and the radial position of the false lumen to the true lumen α_FL (representing different possible positions of the false lumen). The two new parameters have been considered uniformly distributed since knowledge regarding the dimension and position of the false lumen is not available prior to the measurement, and the uniform distribution better represents the lack of knowledge regarding a model variable. The description of the input space for both case studies is given in Table 1. Furthermore, since in the case of aortic dissection the damage factor DF also affects the conductivity changes, the blood conductivity results in: (16) σ(v(t),H)=σ0(H)+Δσ(v(t),H)⋅(1−DF) \sigma (v(t),H) = {\sigma _0}(H) + \Delta \sigma (v(t),H) \cdot (1 - DF) where σ₀(H) is the initial conductivity of blood for a given haematocrit level. In equation (16), DF = 0 corresponds to the conductivity of the flowing blood in the aorta for the healthy condition.

The two models produce a measurement of the impedance cardiograms for each source electrode pairs N_inj and for each measuring electrode pairs N_meas at each time step N_t, as described earlier. The time interval is limited to the first half of the cardiac cycle. After setting the input and output space characteristics, the next step is to set up the options for the PCE. It is essential to notice that a minimum number of point evaluations N_s is needed to have an accurate expansion. N_s is defined to satisfy (17) Ns≥2⋅(M+p)!(M!p!) {N_s} \ge 2 \cdot {{(M + p)!} \over {(M!p!)}} where M is the input space dimension, and p is the order of the polynomial. The healthy condition model has N_s equal to 100, while for the dissected condition N_s is assumed 150. To ensure proper construction of the metamodel and an accurate representation of the interaction between the random input variables, polynomial order of 4 has been chosen for both models. Such a degree guarantees the presence of enough interaction terms in the expansion.

Two PCE functions HC^PCE(t) {\widehat {HC}^{PCE}}(t) and DC^PCE(t) {\widehat {DC}^{PCE}}(t) are introduced. They express the impedance cardiograms of the healthy condition (HC) and the dissected condition (DC), respectively: (18) HC^n,mPCE(t)=−|dZ_dt|HC \widehat {HC}_{n,m}^{PCE}(t) = - {\left| {{{d\underline Z} \over {dt}}} \right|_{HC}} (19) DC^n,mPCE(t)=−|dZ_dt|DC \widehat {DC}_{n,m}^{PCE}(t) = - {\left| {{{d\underline Z} \over {dt}}} \right|_{DC}} where n and m are the indices related to the source and the measurement electrode pairs, respectively. The difference Y^n,mPCE(t) \widehat Y_{n,m}^{PCE}(t) between (22) and (23) is then calculated as (24) (20) Y^n,mPCE(t)=HC^n,mPCE(t)−DC^n,mPCE(t) \widehat Y_{n,m}^{PCE}(t) = \widehat {HC}_{n,m}^{PCE}(t) - \widehat {DC}_{n,m}^{PCE}(t) which expresses the PCE output of the difference between the healthy and dissected metamodels for each electrode configuration and time step.