Complexity reduction in human atrial modeling using extended Kalman filter

Abstract

Human atrial tissue electrophysiology is modeled upon biophysical details obtained from cellular level measurements. Data collected for this purpose typically represent a unique state of the tissue. As reproducing dynamic cases such as subject-varied and/or disorder-varied electrophysiological properties is in question, such complex models are typically hard to use. Hence, there is a need for simpler yet biophysically accurate and mathematically tractable models to be used for case-specific reproductions and simulations. In this study, a scheme for parameter estimation of a phenomenological cardiac model to match a targeted behavior generated from a complex model is used. Specifically, an algorithm incorporating extended Kalman filter (EKF) into the scheme is proposed. Its performance is then compared to that of particle swarm optimization (PSO) and sequential quadratic programming (SQP), algorithms that have been widely used for parameter optimization. Both robustness and adaptability performance of the algorithms are tested through various designs. For this, reproducing action potential (AP) waveforms of varying remodeling states of atrial fibrillation (AF) at different stimulus protocols was targeted. Also, randomly generated initial parameter sets are included in the tests. In addition, AP duration (APD) restitution curve (RC) is used for a multiscale evaluation of fitting performance. Finally, wavefront propagation on 2D of a selected AF remodeling state using parameter solutions from each of the algorithms is simulated for a qualitative evaluation. In general, PSO yielded superior performances than EKF and SQP with respect to fitting AP waveforms. Considering both AP and APD RC, however, EKF yielded the best accuracies. Also, more accurate spiral wave reentry is obtained with EKF. Overall, EKF algorithm yielded the best performance in robustness and adaptability.

Appendices

Appendix A: Minimal model equations

The states for minimal model is as follows:

$$ \dot{u}=I_{st}-(I_{\text{fi}}+I_{\text{so}}+I_{\text{si}}) $$

(29)

$$ \dot{v}=(1-H(u-\theta_{v}))(v_{\infty}-v)/\tau_{v}^{-}-H(u-\theta_{v})v/\tau_{v}^{+} $$

(30)

$$ \dot{w}=(1-H(u-\theta_{w}))(w_{\infty}-w)/\tau_{w}^{-}-H(u-\theta_{w})w/\tau_{w}^{+} $$

(31)

$$ \dot{s}=((1+tanh(k_{s}(u-u_{s})))/2-s)\tau_{s} $$

(32)

where u is membrane voltage, v, w, and s are gating variables. Here u is dimensionless.

Three currents are defined as follows:

$$ I_{\text{fi}}=-vH(u-\theta_{v})(u-\theta_{v})(u_{u}-u)/\tau_{\text{fi}} $$

(33)

$$ I_{\text{so}}=(u-u_{o})(1-H(u-\theta_{w}))/\tau_{o}+H(u-\theta_{w})/\tau_{\text{so}} $$

(34)

$$ I_{si}=-H(u-\theta_{w})ws/\tau_{\text{si}} $$

(35)

where I_fi is the fast inward, I_so is the slow outward, I_si is the slow inward current, and H(x) is the standard Heaviside function. The time constants are defined as follows:

$$ \tau_{v}^{-}=(1-H(u-\theta_{v}^{-}))\tau_{v1}^{-}+H(u-\theta_{v}^{-})\tau_{v2}^{-} \tau^{-}_{w}= $$

(36)

$$ \tau^{-}_{w}=\tau^{-}_{w1}+(\tau^{-}_{w2}-\tau^{-}_{w1})(1+\tanh(k^{-}_{w}(u-u^{-}_{w})))/2 $$

(37)

$$ \tau_{\text{so}}=\tau_{\text{so}1}+(\tau_{so2}-\tau_{\text{so}1})(1+\tanh(k_{\text{so}}(u-u_{\text{so}})))/2 $$

(38)

$$ \tau_{s}=(1-H(u-\theta_{w}))\tau_{s1}+H(u-\theta_{w})\tau_{s2} $$

(39)

$$ \tau_{o}=(1-H(u-\theta_{o}))\tau_{o1}+H(u-\theta_{o})\tau_{o2} $$

(40)

and the infinite values are defined as follows:

$$ v_{\infty} = \left\{\begin{array}{ll} 1, \quad u<\theta^{-}_{v} \\ 0, \quad u>=\theta^{-}_{v} \end{array}\right. $$

(41)

$$ v_{\infty} = (1-H(u-\theta_{o}))(1-u/\tau_{w\infty})+H(u-\theta_{o})w^{*}_{\infty} $$

(42)

Appendix B: Minimal model parameter regimes

Certain boundaries for each parameter is used in this study. Boundary for each parameter is determined in succession using following procedure. The parameter whose boundaries to be determined is varied as all other parameters are fixed in their nominal values. Through variation corresponding AP waveforms are generated and whenever a non-physiological AP is produced, the value before that is set as the boundary. The variation was towards the both directions, plus and minus infinity, from the nominal value. Resolution of variation differed accordingly with the magnitude of the nominal value of varied parameter. Endocardium parameter set provided in [8] is used as nominal parameter values.

Non-physiology of AP is detected based on certain criteria: (i) whenever APD is less than 90 ms or greater than 350 ms, (ii) whenever the value of u is greater than 1.8 during early repolarization phase, (iii) whenever the value of u is greater than 1.5 during the phases other than early repolarization.

Table 2 Upper and lower bounds of MM parameters

The correlations between the parameter groups are neglected in this study to simplify the assessment of boundary regime. Also, sensitivity of a certain parameter is expected to depend on a few number of other parameters since each phase of AP waveform in MM is governed by distinct clusters of parameters. However, parameter sensitivity of MM requires another independent study which is the beyond the scope of this one.