Cardiac tissue conductivity estimation using confirmatory factor analysis

Impaired electrical conduction has been shown to play an important role in the development of heart rhythm disorders. Being able to determine the conductivity is important to localize the arrhythmogenic substrate that causes abnormalities in atrial tissue. In this work, we present an algorithm to estimate the conductivity from epicardial electrograms (EGMs) using a high-resolution electrode array. With these arrays, it is possible to measure the propagation of the extracellular potential of the cardiac tissue at multiple positions simultaneously. Given this data, it is in principle possible to estimate the tissue conductivity. However, this is an ill-posed problem due to the large number of unknown parameters in the electrophysiological data model. In this paper, we make use of an effective method called confirmatory factor analysis (CFA), which we apply to the cross correlation matrix of the data to estimate the tissue conductivity. CFA comes with identifiability conditions that need to be satisfied to solve the problem, which is, in this case, estimation of the tissue conductivity. These identifiability conditions can be used to find the relationship between the desired resolution and the required amount of data. Numerical experiments on the simulated data demonstrate that the proposed method can localize the conduction blocks in the tissue and can also estimate the smoother variation in the conductivities. The conductivity values estimated from the clinical data are in line with the values reported in literature and the EGMs reconstructed based on the estimated parameters match well with the clinical EGMs


Introduction
The electrical conductivity of cardiac tissue plays an important role in the origin of heart rhythm disorders.It has been shown that impaired electrical conduction and damage of atrial tissue are related to the development and progression of cardiac tachyarrhythmias [1][2][3][4].For example, reduced conductivity facilitates development of atrial fibrillation, which comes with a high risk of strokes [1,2,5].Being able to infer the conductivity, is important to localize the arrhythmogenic substrate that causes atrial fibrillation.In this paper we focus in particular on estimating the conductivity for atrial tissue based on epicardiac electrograms (EGMs).
Tissue conductivity determines the ability of the tissue to transport electrical charge [6] and is related to the propagation of transmembrane current through the cardiac tissue.Previous research has proposed fully experimental methods or mathematical models to determine the conductivity from data, and the interested reader is referred to the reviews for more details [7][8][9].One of the classic approaches is based on the cable theory, which relates the changes in the transmembrane potential to changes in the total ionic current and the conductivity of the tissue [6].Based on this theory, bidomain models and monodomain models have been proposed to estimate tissue conductivity.The monodomain is a special case of the bidomain models which needs less model parameters and has been used extensively to simulate myocardiac cells [10][11][12] and nerve cells [13].
With the development of high resolution electrode arrays, it is possible to measure the transmembrane potential of the tissue at multiple locations in a parallel manner [14][15][16].With these high resolution spatial-temporal measurements, we can obtain a deeper understanding on the underlying structure of the tissue.In this work, we use the epicardial electrogram data measured using a high resolution mapping approach presented in [16].However, inferring the conductivity parameters from such large amount of high dimensional data by inverse techniques is highly challenging.It is due to many unknown model parameters that need to be estimated first [17] and also due to the high dimensionality, nonlinearity, and stochasticity of the ill-posed inverse problems [12].The work in [12] proposed a forward model using a monodomain approach and estimated tissue conductivity from measurements obtained by microelectrode arrays by solving a complex inverse problem.It is computationally intensive since ionic currents for all cells are needed to calculate the action potentials for conductivity estimation.To solve this problem, the authors in [18] proposed a compact model based on the work in [12,19] to reduce the cost for computing the ionic currents of all cells and used the model to estimate the tissue conductivity from epicardial electrograms measured by high resolution mapping arrays.It is based on the assumption that all cells generate the same stereotype action potential once activated.This assumption was also used in [20] to reduce the computational intensity in the ECG forward model.However, the work in [18] assumes that the anisotropy ratio is fixed everywhere in the tissue and known before estimating the conductivity parameters.Also, it does not explicitly make use of the spatial structure of the data.In this work, we do not use the fixed and known anisotropy ratio assumption, but estimate the anisotropy ratio and conductivity parameters for all cell positions jointly opposed to [18].We also make use of the spatial structure of the multi-electrode data to improve conductivity estimation.To do so, we exploit the cross power spectral density matrix (CPSDM) of the EGMs and estimate the conductivity parameters using the CPSDM in combination with confirmatory factor analysis [21] (CFA).
CFA is a type of structural equation modelling that deals with the relationships between observed measurements and latent variables or factors [21].With CFA, the observed measurements are intercorrelated and the variables or factors influence the correlations among the observations.In practice, CFA is usually confined to analyze variance-covariance structures and the parameters in the CFA model (such as factor loadings, factor variances and covariances), are estimated to reproduce the input variance-covariance matrix.Therefore, estimating the parameters from the CPSDM model (such as the power spectral densities of the sources) can be regarded as a CFA problem.For earlier works using CFA to estimate the parameters of a CPSDM model, one can refer to [22][23][24][25].
To estimate the parameters in the CPSDM model, the measurement model must first be identified.By analyzing the identifiability conditions in the CFA problem, we can verify whether the solution can be obtained with the given amount of data.In the current application this helps to analyze the relationship between the resolution that can be obtained for conductivity estimation and the required number of electrodes.
In order to easier satisfy the identifiability conditions in the CFA problem, we further propose simultaneous CFA across multiple temporal frequencies to estimate the tissue conductivity.This is based on the fact that the conductivity parameters are shared among multiple frequencies, which can increase the ratio between the number of equations and the number of unknown parameters.The proposed method also uses non-linear constraints to reduce the feasibility set of the parameter space with the prior knowledge of the parameters and thus increases the robustness of the estimation.The method also makes it possible to include constraints related to the mathematical properties of the problem at hand, e.g., to guarantee the positive semidefinite property of the CPSDMs.
The rest of the paper is organized as follows.In Sec.II, we introduce the notation used in this paper and introduce the EGM models.In Sec.III, we review the basic CFA theory and propose the simultaneous CFA for conductivity estimation and take some practical problems into account to increase the robustness of the proposed method.In Sec.IV, we conduct conductivity estimation experiments on the simulated data and discuss the results.In Sec.V we apply the estimation algorithm to the clinical data and evaluate its performance.We draw the conclusions in Sec.VI.

Electrogram model
In this section, we first introduce the notation used in this paper.Secondly, we introduce the background knowledge of electrical propagation in heart tissue and the EGM model that we use.After that, we introduce the time-domain impulse response model for the atrial EGM proposed in [18] and finally propose the CPSDM model, based on the impulse response model for conductivity estimation.

Notation
We use lower-case letters for scalars, bold-face lower-case letters for vectors, and bold-face upper-case letters for matrices.For example, a matrix A can be written as A = [a 1 , …, a m ], where a i represents its ith column and a ij represents an element of matrix A at position (i, j).The vector formed from the diagonal of the matrix A ∈ C m×m is denoted as Diag(A) = [a 11 , …, a mm ] T .A Hermitian positive semi-definite matrix is denoted as A⪰0.The Frobenius norm of a matrix is denoted as ‖ ⋅‖ F .

Electrical propagation & electrogram model
Pacemaker cells in heart tissue can initiate and conduct impulses to excite neighbouring cells.These atrial cells can conduct the impulses, generating a spread of excitation.When a cell is excited, its transmembrane potential increases until a threshold is reached and an action potential pulse is generated.
To describe the electrical propagation in a computational model, the tissue is discretized into a two-dimensional grid with N regular elements to model N cells.We model the electrical propagation in heart tissue using the cable theory and the monodomain approach, which assumes that the intracellular and extracellular domains have equal anisotropy ratios.Previous research has shown that the monodomain model can be a good approximation for its biodomain counterpart, even if the equal anisotropy assumption does not hold [26][27][28].Electrical propagation in monodomain cardiac tissue at time t is governed by the reaction-diffusion equation [6]: where v n (t) is the transmembrane potential of the nth cell, C is the membrane capacitance, I st,n (t) is the stimulus current, I ion,n (t, v(t)) is the ionic current, and I tm,n (t) is the cell-to-cell transmembrane current per unit area that accounts for the spatial evolution of the action potential.
We consider M electrodes indexed by m ∈{1, 2, …, M} that are positioned on the atrial area.An atrial EGM measures the changes in the action potential of a group of cells close to the electrode.We use r m,n to denote the distance between the mth electrode and the nth cell and let ] T , the atrial EGM at the mth electrode can be modelled as [18,19].
with a the area of each grid element, σ e the constant extracellular conductivity, S v the cellular surface to volume ratio, and D σ a double differentiation operator defined by [18].
Here, D x and D y are the first derivative operators.The conductivity vectors σ xx , σ yy , σ xy and σ yx stack the conductivity of all cells in the xdirection, the y-direction, and the two diagonal directions, respectively.
For example, If the tissue is homogeneous and isotropic, D σ is a scaled Laplacian operator.
Eq. ( 2) can then be simplified as To model the EGM with the model in (4), the action potentials of all M. Sun et al. cells are needed.The reaction-diffusion equation in (1) is often used to calculate the action potential of the cells, but it is complicated and timeconsuming to solve.To solve this problem, it was proposed in [18] to assume that all cells generate the same stereotype action potential once activated and used a simplified model -the impulse response model -to calculate the EGM.Next, we briefly introduce this model.

Impulse response model
Taking the cell with action potential v 0 (t) as reference, the action potential of the cell with index n is modelled as where * is the convolution operator, δ(⋅) is the Dirac delta function, and τ n is the time delay of the nth cell with respect to the reference cell.This implies that the action potential at cell n is modelled by a stereotype action potential that is only delayed compared to the reference cell.With τ 0 for the reference cell set to τ 0 = 0, the action potential of all cells can then be calculated as where Therefore, the atrial EGM measured by the mth electrode [ref.(4)] can be remodelled as This can be further simplified as with a m (t) = Qr T m D σ δ(t) the impulse response from all cells to the sensor at position m.We can see that the parameters of interest, i.e., the conductivities, are included in a m (t).
In addition to the atrial activity, the electrodes at the atrial area also measure the action potential of the ventricular cells and some noises.In this paper, we are only interested in the parameters of the atrial tissue and aim to estimate the parameters when atrial activity is present.To simplify the problem, the ventricular activity is not taken into account in the model.Uncorrelated sensor noise is always present and reflects the limited accuracy of a sensor.The sensor self noise of the mth electrode is denoted by u m (t).Finally, the atrial EGM measured by the mth electrode is modelled as Due to the convolution operator in (10), it would be more convenient to estimate the parameters in the frequency domain via the short-time Fourier transform (STFT).

CPSDM model
To estimate the conductivity, we first deconvolve the time domain model in the STFT domain.In the STFT domain, the EGM model from (10) can be expressed as where l is the time frame index, k is the frequency bin index, and ỹm (l,k), ãm (l,k), ṽ0 (l,k), and ũm (l, k) are the STFT coefficients of y m (t), a m (t), v 0 (t), and u m (t), respectively.Stacking all electrodes in a vector, the EGM model in the STFT domain using vector notation is given by where and where τ n (l) represents the time delay of the nth cell with respect to the reference cell in the lth frame.
The CPSDM of the EGM in the lth frame and the kth frequency is then given by where is the PSD of the reference cell, and the PSD of the sensor self noise of the mth sensor, which is assumed to be uncorrelated across sensors, across time and across frequency bins.
The unknown parameters of interest, i.e., conductivity vectors, are included in the spatial operator D σ in ã(l, k) [ref.(3)].We assume that The ratio of the conductivity in the y-direction and the conductivity in the x-direction is called the anisotropy ratio.We denote the anisotropy ratio at the nth cell position as α n and stack the anisotropy ratios at all cell positions as a vector α = [α 0 , α 1 , …, α n− 1 ] T .Let σ xx = σ, then the operator D σ now can be calculated as For isotropic tissue, the anisotropy ratio is fixed to α = 1.

Problem formulation
The goal of this paper is to estimate the parameters σ and α, using only the noisy CPSDMs Φ y (l, k) estimated from the EGMs and possibly having an initial estimates of φ(l, k), the activation time τ n , and/or Φ u (l, k).
To solve the problem, the number of knowns should be larger than the number of unknowns.The number of knowns depends on the number of electrodes.Given a certain number of electrodes, the number of unknowns depends on the chosen spatial resolution.In addition to estimating the conductivity, we also like to analyze the relation between the number of electrodes and the resolution that we can obtain.This can be achieved with the help of CFA, which we will introduce in the next section.To simplify the notation, we will neglect the indices of time samples, frequency bins and time-frames wherever possible.

Confirmatory factor analysis
In this section, we first review the CFA problem for parameter estimation and the identifiability conditions.Then, we adapt the CFA formulation to our problem and introduce the parameter estimation problem in the CPSDM model.To easier satisfy the identifiability conditions, we further propose simultaneous CFA, which estimates the parameters across multiple temporal frequencies.Practical considerations are also discussed at the end of this section.

The general CFA problem
CFA methods have been proposed to estimate the parameters of the following model [21][22][23][24]: where Φ y is an M × M variance-covariance matrix of the measurements, A is an M × r matrix of unknown factor loadings, Φ is an r × r variancecovariance matrix of the r common factors, and Φ u is an M × M variancecovariance matrix of the residuals.The factors and the residues are assumed to be uncorrelated.The residues are also assumed to be mutually uncorrelated, i.e., Φ u = Diag([q 1 , …, q M ] T ).As an example, in [25], CFA is used to jointly estimate the parameters in the multi-microphone signal model, and the matrices Φ y , Φ, and Φ u are interpreted as the CPSDMs of the noisy signals of M microphones, the r sources, and the noises, respectively, while A was used in [25] to model the early relative acoustic transfer functions of the sources with respect to the microphones.In the current work, we assume a single source, i.e., the action potential v 0 of the reference cell as introduced in (5), and therefore, r = 1 and Φ = φ [ref.(17)].The matrices Φ y and Φ u are the CPSDMs of the EGMs and the sensor self noises, respectively, and A models the transfer function of the source, which describes the propagation from a source (reference) cell to the electrode.CFA can be used to calculate the parameters in (20) using the following constrained optimization problem  where F( Φy , Φ y ) is a cost function such as the maximum likelihood, least square, or generalized least square, and where the constraints are used to incorporate model knowledge, for example the a priori knowledge from (20) and incorporating the knowledge that Φ is positive definite, and Φ u is diagonal and positive definite.In addition, constraints can be put on the elements of A and Φ.This is expressed by the last two constraints in (21).In these constraints  and  are the sets of the selected row-column index-pairs of the matrices A and Φ, respectively, with their elements a ij and φ i ′ j ′ fixed to some known constants âij and φi ′ j ′ .Note that the problem in ( 21) is not convex and may have multiple local minima.
There are two necessary conditions for the parameters of the CPSDM model to be uniquely identifiable.The first identifiability condition is that the number of equations should be larger than the number of unknowns.Therefore, some of the elements in A and Φ are often fixed to make the remaining variables uniquely identifiable.
For guaranteeing unique identifiability, the first identifiability condition is not enough and the second identifiability condition is needed.For any arbitrary non-singular matrix T ∈ C r×r , we have Therefore, there are infinite optimal solutions to the CFA problem in (21).Because there are r 2 variables in T, the second identifiability condition of the CPSDM model in (20) is that we need to fix at least r 2 parameters in A and Φ, i.e.
This condition is necessary but not sufficient.Because we need to fix the proper parameters and not just any r 2 parameters such that the only possible matrix T is T = I.By analyzing the identifiability conditions in the CFA problem, we can determine the resolution that we can obtain with a given number of electrodes.In the following we will formulate the CFA problem for conductivity estimation and propose the simultaneous CFA for the robust estimation of the conductivity with some practical considerations.

CFA problem for conductivity estimation
Estimating parameters from the EGM CPSDM model in (17) from Section II-D can be regarded as a special case of the general CFA problem.To simplify the problem, we fix the PSD of the reference cell φ by a known constant φ, based on the action potential model of a cell.For the reference cell we set τ 0 = 0. We estimate the activation times τ n of the cells using the steepest descent approach [29].With this approach, the activation time of the cells under the electrodes are determined using the steepest descent method followed by interpolation of the activation time for other cells.With this prior knowledge, we propose the following CFA problem formulation for the conductivity estimation at a particular frequency band (For ease of notation the frequency index is left out): The loss function in (24) is defined using the least square cost function as where ‖ ⋅‖ F represents the Frobenius norm.
Prior to solving the problem in (24), we need to analyze the number of knowns and unknowns in (24).We start the analysis by taking one frequency bin into account.Since the estimated noisy CPSDM Φy is Hermitian symmetric, there are M(M+1) 2 known values, while there are 2 N unknowns due to σ and α in D σ and M unknowns due to Φ u .Therefore, the first identifiability condition is given by This means that the number of cells N for which we can determine the conductivity is constrained by M as N ≤ M 4 (M − 1).In our problem, the number of sources r = 1.According to the second identifiability condition, we need to impose r 2 = 1 independent restrictions on ã and φ.Since we already fixed φ, the second identifiability condition is satisfied.

Simultaneous CFA in multiple frequencies
In the previous section we considered estimation of the conductivity parameters σ and the anisotropy ratio parameters α observing the EGMs in only a single frequency band.However, note that we can assume that σ and α are roughly constant across different frequencies within a certain range.By using multiple frequency bands, we can increase the ratio of known-to-unknown variables, when estimating the conductivity using the CPSDM.This allows to estimate the conductivity of more cells using the same number of electrodes.This can be done by solving the

Table 2
Comparison of estimation error with different methods.

Tissue type
With/Without noise CMM (α is known) SCFA (α is known) SCFA (α is estimated) where and the total number of unknowns is 2 The first identifiability condition is thus given by Without sensor self noise, the first identifiability condition becomes Comparing (28) to (26), we can see that the SCFA problem needs less sensors than the CFA problem to satisfy the first identifiability condition when |S f | > 1.The second identifiability condition in the SCFA problem is also always satisfied since r = 1 and we have fixed φ(k) for all k.

Practical considerations 3.4.1. Over-determined & model-mismatch problem
Increasing the ratio of the number of equations over the number of unknowns makes the model fit better to the measurements when assuming the CPSDM model is accurate enough.The number of knowns are determined by the number of measurement positions, i.e., the sensors, while the number of unknowns is determined by the resolution at which we would like to estimate the conductivity.At the highest resolution this means the number of unknowns is determined by the number of cells.One way to increase this ratio is by reducing the resolution at which we estimate the conductivity.
Since the number of cells is large, estimating the conductivity at the finest resolution results in an unsolvable problem, i.e., the first identifiability condition in Eqs. ( 26) or ( 28) is unsatisfied.Therefore, instead of considering individual cells, we consider groups of cells and use the position of the center cell of the group to denote the position of the group.From now on we use the symbol N to denote the number of cell groups.The array that is used to measure the EGMs has electrode dimensions M x × M y .The inter-electrode distance is constant and is denoted by d.The total number of electrodes is M = M x × M y .The tissue is first discretized into M x × M y equal square areas and each electrode is assumed to be at the center of the square area.Then each area is discretized into Z x × Z y regular grid points.One grid point models one cell group.Let N x = Z x × M x and N y = Z y × M y , which means the total number of cell groups is N = N x × N y = Z x × M x × Z y × M y .We define the spatial resolution of the cells as the number of cell groups in a given tissue area.Then the spatial resolution of the cells for which the conductivity is to be obtained, i.e., the number of unknowns to be estimated, is equal to N x × N y .Let R denote the resolution scale compared to the number of electrodes, then R = Z x × Z y .Fig. 1 shows an example.
Rewriting the first identifiability criterion now in terms of R, the first identifiability condition in (28) becomes Decreasing R decreases the resolution N (increases the number of cells per group) and reduces the number of unknowns.However, it also reduces the model accuracy and increases the model-mismatch error compared to the conductivity of the actual cells.
We can also increase the ratio of knowns-to-unknowns by exploiting the fact that the conductivity is constant over frequency and estimate the conductivity simultaneously over multiple frequencies as proposed in Section III-C.As the impact of model-mismatch problem can be different at different frequencies, we propose to use the frequencies that are less influenced by the decrease of the resolution.The measurements of the EGM can be regarded as the output of a low-pass filter on the heart tissue.Therefore, the high-frequency components of the EGM are less accurate than the low-frequency components and it is better to use low frequency components to estimate the parameters.This will be further investigated in Section V.
The presented framework depends on the estimated noisy CPSDM Φy , which we obtain using the maximum likelihood estimate Φy = ỹỹ H .Note that as this is an estimate, it also increases the model-mismatch error.

Box constraints on the parameters
To overcome CPSDM estimation errors and model-mismatch errors, we propose additional constraints on the parameters.Extra robustness can be achieved if the conductivities and the anisotropy ratios are boxconstrained as σ low ≤ σ ≤ σ upper and 0 < α ≤ 1, where σ low and σ upper are real valued upper and lower bounds of conductivity, respectively.Since the conductivity is non-negative, we can simply set σ low = 0.For σ upper , we select proper values based on the previous research on the conductivity of atrial tissue.

Solver
The CFA problem that we propose is a non-convex optimization problem, which can be solved with various solvers.We used the MAT-LAB optimization toolbox and the 'interior-point' algorithm to solve the problem, which is based on a combination of the methods in [30][31][32].

Numerical experiments on the simulated data
To verify the performance of the proposed method, we need the ground truth on the conductivity.As for clinical data, the true values of the physiological parameters are unknown, we perform numerical experiments on simulated data in this section, after which we perform experiments with clinical data in Sec.V.

Data generation
In order to generate simulation data, we need to model individual cells, as well as cells integrated in larger tissue.In this work, we use the Courtemanche model proposed in [33] to model the atrial cells and calculate the total ionic current.Once the cell model is implemented, we can couple the cells together to form the tissue.We define the tissue to be two-dimensional and discretized into 90 × 90 regular grid points to model cells with 0.02 cm cell-to-cell distance.Notice that these cells are actually larger than real cells, which is done to reduce the computational complexity and memory requirements.To demonstrate the proposed algorithm under various conditions, we generate several isotropic and anisotropic tissues with different settings for conduction blocks, including: (i) homogeneous and isotropic tissue without areas of conduction block (conductivity is σ n = 1.1 mS/cm for all n and α = 1); (ii) inhomogeneous tissue with one area of conduction block and the area outside the conduction block is homogeneous and isotropic as in (i); (iii) inhomogeneous tissue with two areas of conduction block and the area outside the conduction block is homogeneous and isotropic as in (i); (iv) inhomogeneous tissue with one area of conduction block and the area outside the conduction block has a smooth variation in conductivity, varying from 0.77 mS/cm to 1.1 mS/cm and has a constant anisotropic ratio α n = 0.4; (v) inhomogeneous tissue with two areas of conduction block and the area outside the conduction block has a higher variation than (iv) in conductivity which varies from 0.75 mS/cm to 1.3 mS/cm and its anisotropic ratio α n changes from 0.7 to 0.5 from left to right.
Before the tissue is activated by a stimulus, all the cells are at the initial conditions specified by the Courtemanche model.Then, we apply an external stimulus current I st = 200 μA to some source cells at the upper left corner of the tissue with duration of 0.5 ms, which triggers the propagation of the action potential through the tissue.By solving the reaction-diffusion equation in (1) using the forward Euler's method with a time step of 0.05 ms, we can calculate the action potential at every point of the tissue.The moment when a cell's potential reaches the threshold of − 40 mV is considered as the cell's activation time.Since the electrical wave should not propagate further than the tissue boundaries, we take here into account the no-flux boundary conditions.The values of the other parameters that we used for the experiments are summarized in Table 1.
Fig. 2(a) plots the central area of 45 × 45 cells and a 5 × 5 electrode array with inner-electrode distance 0.18 cm on top of this area.The EGMs are generated using (2).In the simulation experiments, we start with the case without noise added to the atrial EGM.Then, we generate noisy EGMs by adding Gaussian noise at 50 dB signal-to-noise ratio (SNR) to simulate sensor self noise.We also estimate the conductivity at different resolutions.As a reference, the proposed method is compared to the conductivity estimation method from [18].

Model testing
To reduce the number of unknowns in the estimation problem, we consider groups of cells as we introduced in Section III-D.For example, we take a block of 3 × 3 initial cells as a cell (group) in our model, as shown in Fig. 2(b).Then we use the conductivity of the central cell of the block to replace the conductivity of the cell group in the low-resolution model.The spatial resolution of the cells to be estimated reduces from As we discussed in Section III-D, a reduced spatial resolution will increase the model-mismatch error.For comparison, we further generate the atrial EGMs using the impulse response model with different spatial resolutions and compare the generated EGMs.Fig. 3(a) shows the conductivity map of the tissue, Fig. 3(b) shows the activation time map of the tissue, and Fig. 3(c) gives examples of the cell action potential v(t) and Fig. 3(d) plots the spectral amplitude of the reference cell, respectively.The spectral amplitude of the reference cell is used to calculate parameter φ in the CPSDM model.Fig. 3(e) shows the EGMs generated using different models and different resolutions, including 'E1': the EGM generated based on the reaction diffusion model in (2) using the original resolution (45 × 45), 'E2': the EGM generated based on the impulse response model in (8) using the original resolution (45 × 45), and 'E3': the EGM generated based on the impulse response model in (8) using the lower resolution (15 × 15), respectively.Fig. 3(f) shows the spectral amplitudes of the EGMs in Fig. 3(e).We observed that the generated EGMs have small differences in the morphology and decreasing the resolution increases the difference, as shown in Fig. 3(e).From Fig. 3(f) we also see that the model-mismatch error increases when the spatial resolution decreases, and the high frequency components are more affected than the low frequency components.Motivated by this, we use the lower frequency components rather than the higher frequency components as long as the identifiability conditions are satisfied.This can reduce the negative effect of the model-mismatch problem, which is unavoidable in practice.Some experiments are done as a function of the number of frequency components involved in the estimation algorithm.In that case, the number of frequency bands involved is increased by adding frequency bands towards the higher frequencies.
From Fig. 3(f), we observe that the model-mismatch error is relatively large at 0 Hz.This is due to the assumption made in the impulse response model that all cells generate the same stereotype action potential once activated, which ignores the very small change in the waveform from cell-to-cell.In the experiments, we do not use the 0 Hz frequency component when estimating the parameters.

Evaluation
We estimate the conductivity at a lower resolution than what was used for data generation.In this way we can study the effect of the model-mismatch as well.We consider the resolution is R times the number of electrodes and set R = 1.8, R = 2.25, R = 3, and R = 4.5 by taking 5 × 5 cells, 4 × 4 cells, 3 × 3 cells, and 2 × 2 cells per block, respectively.The larger R represents the higher resolution, which implies a smaller model-mismatch.However, having a higher resolution also makes it harder to satisfy the identifiability conditions.
We calculate the normalized mean square error (NMSE) to quantify the estimation performance.When the tissue is isotropic, that is, the longitudinal conductivity is equal to the transverse conductivity, the estimation error is equal in the two directions.In this case, the NMSE is calculated as Fig. 6.Examples of real and estimated activation time maps and conductivity maps.In the first row, the subfigures (iii-1)/(iii-2) show the real estimated activation time maps with original resolution 45 × 45 and the estimated activation time map with resolution 15 × 15 using the steepest descent method, respectively; the subfigure (iii-3) shows the estimated conduction velocity (CV) map; the subfigures (iii-4)/(iii-5)/(iii-6)/(iii-7) show the real conductivity map, the estimated conductivity map using the CMM method when the anisotropy ratio α is known, the estimated conductivity map using the proposed SCFA method when α is known, and the estimated conductivity map using the proposed SCFA method when α is unknown for tissue type (iii), respectively.The subfigures (iv-1)/(iv-2) and (v-1)/(v-2) show the real and estimated activation time maps for tissue type (iv) and (v), respectively.The subfigures (iv-3) and (v-3) show the estimated CV maps for tissue type (iv) and (v), respectively.The subfigures (iv-4-L)/(iv-4-T), (iv-5-L)/(iv-5-T), (iv-6-L)/(iv-6-T), and (iv-7-L)/(iv-7-T) are the real and estimated longitudinal/transversal conductivity maps for tissue type (iv), and (v-4-L)/(v-4-T), (v-5-L)/(v-5-T), (v-6-L)/(v-6-T), and (v-7-L)/(v-7-T) are the real and estimated longitudinal/ transversal conductivity maps for tissue type (v), respectively.
where N p represents the number of unknowns, σ represents the real conductivity and σ represents its estimation.When the tissue is anisotropic, the longitudinal conductivity is unequal to the transverse conductivity.In this case, we first calculate the errors in the two different directions using (31) and then calculate the averaged error to evaluate the performance.
To focus on the conductivity estimation, we start the evaluation with the case where the activation time is known and there is no sensor self noise.Then we evaluate the algorithm having inaccuracies in activation time.

Assuming known activation time
In this evaluation we first assume the anisotropy ratio is known, and only estimate the conductivity, and next consider the case when both α and σ are unknown and estimated.We use the proposed method to estimate the parameters given different frequency bands.We denote the frequency band for estimation by [1, f max ] Hz and calculate the estimation errors when f max = 1, 10, …, 100 for different resolutions and tissue types, respectively.That is, increasing |S f | [ref.(30)] from 1 to 100.The iterations in the algorithm stops when the change of the objective function is small enough.Fig. 4(a) and Fig. 4(b) show the estimation errors for the estimation where the anisotropy ratio is known and for the situation where the anisotropy ratio is unknown and has to be estimated, respectively.Comparing Fig. 4(a) and (b), we see that the algorithm still works well without knowing the anisotropy ratio, although its performance degrades a little bit.We also see that the curves in Fig. 4(a) and (b) show a similar trend.From both figures (a) and (b), we see that for all tissue types the estimation error first decreases and then tends to stabilize with the increase of f max in most cases.When R increases from 1.8 to 3, the resolution increases and the estimation error decreases, because the model-mismatch error gets smaller.In most of the cases, we see the performance degrades when R increases from 3 to 4.5.Although the model is more accurate when R = 4.5, it has more parameters to fit and this estimation problem therefore is more challenging.We also see that the estimation error increases when increasing the number of areas with conduction block, since the underlying structure of the tissue is more inhomogeneous and reducing the resolution of the tissue increases the error in the model.Moreover, when increasing f max , the estimation error decreases faster for low resolution compared to a high resolution.
Fig. 5 shows an example of the true simulated conductivity maps and the estimated conductivity maps given different number of frequency components and different resolutions for the five tissue examples (i)-(v) specified earlier in Section IV-A.Fig. 5(a) shows the results for tissue (i), (ii), and (iii), where the tissue areas outside the blocks are isotropic, i.e., α = 1.Fig. 5(b) shows the results for tissue (iv) and (v), where the tissue areas outside the blocks are anisotropic, i.e., α ∕ = 1.The anisotropic ratio is constant in example (iv) and is spatially-varying in example (v).Therefore, we show in Fig. 5(b) both the longitudinal and transverse conductivity maps.In Fig. 5, the first column shows the real conductivity maps, the second to the sixth columns show the estimated conductivity maps for the increasing resolution with the assumption of knowing the anisotropy ratio, and the last column shows the conductivity maps estimated without knowing the anisotropy ratio.We can observe from the second to the sixth columns that the performance is obviously improved when increasing the resolution from R = 1.8 to a higher resolution and increasing f max from f max = 1 to f max = 100.We also find that the area of the conduction block is more accurate in the high resolution map.From the estimated conductivity maps of tissue type (iv) in Fig. 5 (b), we find that the algorithm also works well when there is smooth variation in the conductivity maps.However, the estimated conductivity maps seem less accurate for more heterogeneous simulations (iv-v), mainly in case of non-constant anisotropy ratios and when the ratio is unknown.This could be a consequence of the fixed stereotype waveform for all cells.

Having inaccuracies on the activation time
Next, the conductivities are estimated using the activation time estimated from the measured EGMs instead of using the true activation time.They are estimated using the steepest descent approach [29], followed by linear interpolation to obtain a higher resolution activation time map for all modelled cells.We now use the proposed method (fixing f max = 100) to estimate the conductivity and use the CMM method from [18] as a reference.With the CMM method, the activation time of the cells is also estimated by the steepest descent method.Note that the CMM method assumes the anisotropy ratio is known and constant, while this is not required in the proposed SCFA method.We here calculate the estimation error for the CMM method taking known anisotropy ratio into account and calculate the estimation error using the proposed method taking both known and unknown anisotropy ratio into account.We also provide the estimated local conduction velocity (CV) map [34], which is commonly used in the literature for analysis of conduction block in tissue.Note that CV does not measure the conductivity parameters but only provides an estimation of the local velocity of the wave propagation in tissue.
Three tissue examples (iii)-(v) specified in Section IV-A are used to test the estimation performance when there is no sensor self noise and with simulated white Gaussian sensor self noise by setting the SNR at 50 dB.Table 2 lists the estimation error for different types of tissue obtained by the proposed SCFA method and the CMM method, respectively.It can be seen that SCFA performs better than CMM in all scenarios even when CMM uses the true α and SCFA estimates α, which implies that using spatial information of the multi-electrode data helps to improve estimation of the conductivity.The performance of the methods degrades a little bit when the EGMs are disturbed by the noise.Fig. 6 shows the real activation time maps, the estimated activation time maps estimated by the steepest descent approach followed by linear interpolation, the estimated conduction map, the real conductivity maps and the conductivity maps estimated by the CMM method with a known α and estimated by the SCFA method with known and estimated α when there is no noise.Compared to CMM, the proposed method localizes the conduction block more accurate and estimates the conductivities of the tissue regions outside the conduction blocks better.We also find that the higher resolution captures the underlaying structure better, which further results in better estimation results.The CV maps only indicate the fast and slow conduction areas, which are not accurate compared to the other two methods and cannot estimate tissue conductivity and infer electropathology in tissue.

Experiments on clinical data
In this section we evaluate the performance of the proposed method on clinical data.The data originate from epicardial EGMs measured from human atria during open-heart surgery using a high-resolution mapping approach [16].The study protocol was approved in February 2010 by the Medical Ethics Committee (2010-054) in the Erasmus Medical Center, Rotterdam, The Netherlands.A mapping array of 8 × 24 electrodes with an inter-electrode distance of 2 mm is used for data collection.The array measures epicardial EGMs during sinus rhythm and during induced atrial fibrillation at each site.The acquired EGMs are amplified, filtered (bandwidth 0.5-400 Hz), sampled (1 kHz), and analogue to digital converted (16 bits).For more details on the mapping approach and the electrode array we refer to [16].The far field artefacts and ventricular activities in the EGMs are cancelled using the method in [35].
Before estimating the tissue conductivity, we need to discretize tissue into groups of cells as introduced in Section III-D.More concretely, the tissue under the mapping array is first discretized into M x × M y = 8 × 24 equal square areas of 2 mm by 2 mm each and each area is discretized into Z x × Z y = 3 × 3 cell groups.Then the spatial resolution of the tissue is N x × N y = 24 × 72.Each electrode is assumed to be at the center of the square area.To reduce the computation cost, we took a smaller tissue area with 5 × 5 electrodes in each estimation and the cell groups to be estimated is N x × N y = 15 × 15.
The activation times of the cells are estimated using the steepest descent approach with interpolation, similar as done in the simulated data.After computing the activation times of the cells, we use the proposed method to estimate the tissue conductivity.Since there is no ground truth of the tissue conductivity, we compare the estimated conductivity parameters to the experimentally determined conductivity values reported in the literature [7,36].However, the experimental results in the literature are inconsistent so far.To further demonstrate the performance of the proposed method, we also compare the real EGMs to the EGMs reconstructed using the estimated conductivity.The reconstruction error is calculated by where y m and ŷm respectively represent the real and the estimated atrial EGM segments during atrial activity, and ‖ ⋅‖ 2 represents the l 2 -norm of a vector.
We compare the performance of the proposed method with the CMM method proposed in [18].Since the anisotropy ratio is unknown and cannot be estimated by the CMM method, we just fix it to one as which has been done in [18].Table 3 shows the reconstructed errors of the EGMs of five patients.We can see that the proposed method achieves better performance than the CMM method.We show an example of the estimation on the clinical data in Fig. 7.The estimated activation time map is plotted in Fig. 7 (a).The estimated conductivity map in the longitudinal and the transverse directions which are obtained by the proposed method are plotted in Fig. 7 (b) and Fig. 7 (c), and the estimated conductivity map obtained by the CMM method is plotted in Fig. 7 (d).Examples of the real and the reconstructed EGMs obtained by the proposed method and the CMM method are shown in Fig. 7 (e).In Fig. 7 (b), (c), and (d), the estimated conductivities outside the conduction block vary from 0.84 mS/cm to 1.66 mS/cm, 0.38 mS/cm to 0.95 mS/cm, and 0.48 mS/cm to 0.88 mS/cm, respectively, which are in line with the values reported in the literature [7,36].From Fig. 7 (e) we see that the reconstructed EGMs are relatively similar to the original EGMs, in particular when only single deflections are present.Compared to the CMM method, the reconstructed EGMs fit the original EGMs better, which further validates the performance of the proposed method.

Conclusions and discussions
In this work we proposed to estimate the tissue conductivity using the cross power spectrum of the multi-electrode EGMs in combination with the confirmatory factor analysis (CFA).Based on the fact that the conductivity parameters are shared among multiple frequencies, we proposed the simultaneous CFA (SCFA) to estimate the parameters of interest in the cross power spectral density matrix model.With SCFA, we are able to determine the resolution that we can obtain with a given number of electrodes.Compared to traditional CFA, SCFA needs less sensors to obtain the same resolution.Experiments on simulated data and clinical data demonstrate that SCFA obtains good estimation of the tissue conductivity.
In this work, several assumptions are made to simplify the problem, which comes with some limitations.For example, we assume a fixed stereotype waveform for all cells although the action potential waveforms are possibly heterogeneous.For example, in the areas of compromised conductivity the upstroke may be slow.Also, we considered two-dimensional tissue which is more appropriate for thin-walled tissues, and assumed the conductivity tensors are aligned along the axes, which cannot model the curved paths of the fibre and the diseased atria with potential dissociation of endo-/epicardial layer.These assumptions reduce the accuracy of the EGM model, but they help to simplify the problem and reduce the computational cost for the estimation.The experiments that account for both isotropic and anisotropic propagation also demonstrate the robustness of the proposed method.We would like to improve the model in the near future, taking into account more realistic conditions, while reducing the computational cost for estimating the parameters.
Previous research has shown that alterations of atrial conduction is related to multiple cardiovascular diseases, such as atrial fibrillation, atrial fibrosis, atrial enlargement, etc.The potential applications of this model can be automated detection of the areas of impaired conduction that possibly need to be ablated in an ablation operation to deal with atrial fibrillation, or staging of the abnormalities of the arrhythmias.

Fig. 1 .
Fig. 1.Discretization of tissue.There are M x × M y electrodes placed on the tissue.Dividing the tissue into M x × M y areas with Z x × Z y cell groups for each area.Each electrode is assumed to be at the center of the area.

Fig. 3 .
Fig. 3. (a) Conductivity map; (b) Activation time map; (c) Examples of cell action potential; (d) Spectral amplitude of the reference cell; (e) Generated EGMs; (f) Spectral amplitudes of the EGMs.In (e), E1 represents the EGM generated based on the reactiondiffusion model.E2 and E3 represent the EGM generated based on the impulse response model using the original resolution (45 × 45) and the low resolution (15 × 15), respectively.In (f), A1, A2 and A3 represent the spectral amplitude of E1, E2, and E3, respectively.

Fig. 4 .
Fig. 4. Estimation errors of conductivity parameters when f max = 1, 10, …, 100 for different resolutions for different types of tissue, respectively.(a) The anisotropy ratio α is known.(b) The anisotropy ratio α is unknown and both σ and α are estimated.For the anisotropic tissue examples (iii), (iv), and (v), we calculate the average estimation errors in the longitudinal and the transverse directions and then plot the average values in the figures.

Fig. 5 .
Fig. 5.The subfigures in each row from the first to the sixth represent the real conductivity map, the estimated conductivity map when R = 1.8 with f max = 1, R = 1.8 with f max = 100, R = 2.25 with f max = 100, R = 3 with f max = 100, R = 4.5 with f max = 100 assuming the anisotropy ratio is known, and the last subfigure shows the conductivity map estimated without knowing the anisotropy ratio when R = 4.5 with f max = 100, respectively.(a) Conductivity maps for tissue types (i)/(ii)/(iii).(b) Longitudinal (L) and transversal (T) conductivity maps for tissue types (iv)/(v).

Fig. 7 .
Fig. 7.An example of the estimation on clinical data.(a) Estimated activation time map.Five blue points represent five electrodes; (b) and (c) Estimated conductivity maps in the longitudinal and the transverse direction by the proposed SCFA method, respectively; (d) Estimated conductivity map by the CMM method; (e) Examples of the real (black) and the reconstructed (red/blue) EGMs at the five electrode positions marked in (a).

Table 1
Summary of parameters used in the experiments.
is the set of the frequency indices for conductivity estimation.Using the noisy CPSDM in |S f | frequency bands, the total number of knowns is

Table 3
Reconstruction errors of clinical EGMs using different methods.