On the Role of Ionic Modeling on the Signature of Cardiac Arrhythmias for Healthy and Diseased Hearts

Abstract

Computational cardiology is rapidly becoming the gold standard for innovative medical treatments and device development. Despite a worldwide effort in mathematical and computational modeling research, the complexity and intrinsic multiscale nature of the heart still limit our predictability power raising the question of the optimal modeling choice for large-scale whole-heart numerical investigations. We propose an extended numerical analysis among two different electrophysiological modeling approaches: a simplified phenomenological one and a detailed biophysical one. To achieve this, we considered three-dimensional healthy and infarcted swine heart geometries. Heterogeneous electrophysiological properties, fine-tuned DT-MRI -based anisotropy features, and non-conductive ischemic regions were included in a custom-built finite element code. We provide a quantitative comparison of the electrical behaviors during steady pacing and sustained ventricular fibrillation for healthy and diseased cases analyzing cardiac arrhythmias dynamics. Action potential duration (APD) restitution distributions, vortex filament counting, and pseudo-electrocardiography (ECG) signals were numerically quantified, introducing a novel statistical description of restitution patterns and ventricular fibrillation sustainability. Computational cost and scalability associated with the two modeling choices suggests that ventricular fibrillation signatures are mainly controlled by anatomy and structural parameters, rather than by regional restitution properties. Finally, we discuss limitations and translational perspectives of the different modeling approaches in view of large-scale whole-heart in silico studies.

1. Introduction

Recent developments in computational modeling of the heart’s bio-electrical activity have established the most highly detailed example of a virtual organ [1,2,3]. These advances are the result of the significant progress in cardiac cell modeling, corroborated with experimentation and clinical practice. Moreover, continuous increases in computational power have led to whole-heart electrical models that have shown promising translational outcomes, improving the general understanding of heart function, its pathologies, and therapies [4,5,6,7,8,9]. Still, critical modeling challenges arise when local heterogeneities at different spatio-temporal scales are taken into account [10,11,12,13,14,15,16].

Various approaches have been proposed to simulate complex electrical behavior of the heart, e.g., ventricular arrhythmias. First, detailed cell models, with highly accurate and validated biophysical relationships representing the ground truth, have been incorporated to improve the physiological relevance of in silico cardiac predictions. Nonetheless, these approaches may result computationally demanding, further requiring advanced optimization tools [17,18,19,20,21,22,23]. Alternatively, reduced phenomenological models were developed, preserving the main features of physiological approaches but still allowing for a robust in silico investigation. Depending on the application at hand, one of the two options is favorable, identifying the specific characters that minimize alterations of the final results [24]. This is especially true in the study of abnormal electrical waves, e.g., ventricular tachycardia (VT) and ventricular fibrillation (VF), in healthy and diseased conditions [25,26].

The computational cardiology literature contains a plethora of biophysical models that describe the electrical behavior of the human myocardium [27]. A detailed description of ion channels, pumps, and exchangers is usually included to explain experimental observations [28]. However, many parameters are highly difficult or nearly impossible to measure in vivo, and different formulations have been shown to replicate similar dynamics. For example, the O’hara Rudy dynamic (Ord) model (41 variables, 16 parameters) and the Ten Tusscher-Panfilov 2006 (TP06) model (19 variables, 48 parameters) produce similar electrical behaviors [29,30,31]. In addition, the multiple time scales involved in the formulation, from milliseconds to seconds, make the resulting dynamic system very stiff, which poses challenges in the numerical solution [32,33,34]. Alternatively, phenomenological models are derived from the overall description of cumulative inward and outward currents across the cell membrane and rely on a smaller set of variables. For instance, the Mitchell-Schaeffer model (2 variables, 9 parameters) [35] and its generalizations [36] or the Fenton-Karma (FK) model (3 variables, 14 parameters) and its later extensions (minimal model—4 variables, 28 parameters) [37,38] are well-known instances of simplified phenomenological descriptions able to reproduce several dynamic properties of the cardiac electrical activity: threshold of excitation, maximum upstroke velocity, action potential (AP) shape and morphology, restitution properties, action potential duration (APD) and alternans dynamics [39,40,41]. Increased numerical efficiency is enforced, lacking detailed ionic descriptions, although well-suited for large-scale whole-heart numerical studies.

Cardiac arrhythmias, such as VT and VF, have been studied in detail using physiologically-based computational approaches. Various cardiac diseases, such as long-QT syndrome or Brugada Syndrome, have been elucidated using advanced computational techniques [42,43,44]. Computational cardiology modeling has also been used to explain the relation between AP shape changes with the likelihood of a reentrant arrhythmia (see e.g., [24] and references therein). Estimates of membrane potential and ionic currents distributions in the core region of a sinoatrial node reentry was obtained through detailed computer simulations [45]. Subject-specific whole-heart computational models explored arrhythmia risk-stratification in patients with myocardial infarction (MI) [9,46], and complex VF dynamics have also been characterized in terms of scroll wave filaments through detailed ventricular models [42]. Likewise, the study of AP vortex dynamics has been conducted via phenomenological descriptions in both simplified [47,48,49] and whole-heart simulations [50,51,52,53]. These studies show that the in silico study of arrhythmias is currently being done either using phenomenological ionic models or biophysical ionic models, but how the choice of ionic model affects the arrhythmic signatures of VF in normal and failing hearts remains understudied. Furthermore, heart simulations based on phenomenological ionic models are expected to take less wall-clock time due to the simplified mathematical representation of simplified ionic models. This has been explored in previous works using the Eikonal equation to simulate action potential propagation [40]. However, using the complete monodomain or bidomain approach, a study of the differences in computing time for biophysical versus phenomenological ionic models in realistic heart geometries has not been fully elucidated.

The scientific question that guides the current study is the following: How does the choice of ionic model affect the dispersion of repolarization, arrrhythmic signatures of VF, and the overall computing time? To answer this question, we constructed subject-specific computational models of normal and failing hearts using two models of ionic transmembrane current, namely the FK and the TP06 models. To analyze the differences between these two modeling choices, we assessed the resulting dispersion of repolarization, VF signatures, and wall-clock time.

2. Materials and Methods

2.1. Subject-Specific Heart Models

We constructed three-dimensional bi-ventricular heart models representative of a healthy normal case (NC) and an ischemic heart failure case (HFC). See Figure 1 for a graphical representation of the cardiac fiber orientation distribution, transmural topology and cellular electrophysiological behavior. Layers from Endo to Epi showed in the figure refer to the parametrization used in the TP06 model. The geometry and structural information of each heart is taken from our previous publication [9].

Structural differences due to pathological remodeling processes appear in fibers orientation and ventricular wall thickness, though the overall organs are volumetrically and morphology similar [54]. Accurate geometry models are reconstructed from ex vivo segmentation of high-resolution swine hearts using MRI and DT-MRI datasets as described in [55,56]. For the HFC model, myocardial infarction was induced by occluding the obtuse marginal branches of the left circumflex artery. Healthy and infarcted regions were identified, and tetrahedral finite elements were used to discretize the bi-ventricular heart domains. In particular, a scalar function $h\left(\mathbf{x}\right)\in [0,1]$ was implemented to represent the volume fraction of healthy tissue, i.e.,: $h\left(\mathbf{x}\right)=1$ represents healthy tissue, $h\left(\mathbf{x}\right)=0$ represents the infarcted zone and $0<h\left(\mathbf{x}\right)<1$ represents the transition zone in which interpolated electrophysiological properties are imposed. A vector field, $\mathbf{f}\left(\mathbf{x}\right)$, representing fiber orientation, was related to the mesh nodes and interpolated to the finite elements.

In the case of the TP06 electrophysiological model we further accounted for the transmural dispersion of repolarization using a Laplace interpolation method [57]. Such a modeling choice allowed us to represent epicardial, mid-myocardial and endocardial layers by using a thickness ratio of 2:3:3 (see Figure 1c), and to modulate the AP morphology accordingly (see Figure 1d).

2.2. Numerical Modeling of Cardiac Electrical Activity

We model the electrical propagation within the cardiac domain $\mathsf{\Omega }$ by using the monodomain approach [58]. The time evolution of the transmembrane potential $V_m$: $\mathsf{\Omega }\times [0,T]\to \mathbb{R}$ is then represented by the balance of transmembrane and diffusive ionic currents, augmented by a set of local kinetics variables, i.e.,

$$\begin{array}{}\mathrm{(1a)}& \hfill \frac{\partial V_m}{\partial t}& =\frac{1}{\chi C_m}\left(\nabla ·\left(\mathbf{D}\nabla V_m\right)+\chi \left(I_{ion}(V_m,\mathbf{r})+I_{stim}\left(t\right)\right)\right)\mathrm{in}\mathsf{\Omega }\times [0,T],\hfill \mathrm{(1b)}& \hfill \frac{\partial \mathbf{r}}{\partial t}& =\mathbf{g}(V_m,\mathbf{r})\mathrm{in}\mathsf{\Omega }\times [0,T],\hfill \end{array}$$

where $I_{ion}$ is the sum of the ionic current that depend on the transmembrane potential and $\mathbf{r}$: $\mathsf{\Omega }\times [0,T]⟶{\mathbb{R}}^m$ is the set of state variables that characterize the selected cell model. $\chi$ is the surface-to-volume ratio and $I_{stim}$ represents the external electrical impulse. The conductivity tensor $\mathbf{D}\in {\mathbb{R}}_{+}^{N\times N}$ is assumed to be transversely isotropic and heterogeneous according to the identification of the infarcted regions

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathbf{D}=h\left(\mathbf{x}\right)\left[d_{\perp }\mathbf{I}+\left(d_{\perp }-d_{\Vert }\right)\mathbf{f}\otimes \mathbf{f}\right],\end{array}\end{array}$$

(2)

where $d_{\perp },d_{\Vert }$ represent the transversal and longitudinal conductivity, respectively, and $\mathbf{f}:\mathsf{\Omega }\to {\mathbb{R}}^3$ is the fiber direction vector field. The conductivity ratio $d_{\Vert }/d_{\perp }$ varies typically between 4 and 9 [37]. We selected a conductivity ratio of 4, establishing the longitudinal conductivity $d_{\Vert }=0.1{\mathrm{mm}}^2/\mathrm{ms}$, with a corresponding transversal conductivity $d_{\perp }=0.025{\mathrm{mm}}^2/\mathrm{ms}$. These conductivities were selected based on the standard values for the conduction velocity reported in the works of Fenton and Karma ($\mathrm{CV}\approx 0.42\mathrm{mm}/\mathrm{ms}$) [37] and Ten Tusscher and Panfilov ($\mathrm{CV}\approx 0.67\mathrm{mm}/\mathrm{ms}$) [30]. We selected these two models since they have been effectively used and validated in the study of VF formation and dynamics, alongside applications of anatomical-based models with transmural fiber rotation and ischemic heart disease [26,59,60].

2.2.1. The Fenton-Karma Model

The Fenton-Karma (FK) model is a 3-variable phenomenological model developed to study arrhythmogenesis in cardiac tissue that reproduce restitution properties and spiral wave dynamics by using a minimal set of field variables [37]. Model equations are formulated in dimensionless form, defining the nondimensional membrane potential $u=(V-V_0)/(V_{fi}-V_0)$, and scaling the ionic currents as $J_X=I_X/\left(C_m(V_{fi}-V_0)\right)$ (measured in 1/ ms). According to the monodomain approach, Equation (1) becomes

$$\begin{array}{}\mathrm{(3a)}& \hfill {\displaystyle \frac{\partial u}{\partial t}}& =\nabla ·(\mathbf{D}\nabla u)-J_{fi}(u,v)-J_{so}\left(u\right)-J_{si}(u,w)+J_{stim},\hfill \mathrm{(3b)}& \hfill {\displaystyle \frac{dv}{dt}}& =H(u_c-u)(1-v)/{\tau }_v^{-}-H(u-u_c)v/{\tau }_v^{+},\hfill \mathrm{(3c)}& \hfill {\displaystyle \frac{dw}{dt}}& =H(u_c-u)(1-w)/{\tau }_w^{-}-H(u-u_c)w/{\tau }_v^{+},\hfill \end{array}$$

where

$$\begin{array}{}\mathrm{(4a)}& \hfill J_{fi}(u,v)& =-vH(u-u_c)(1-u)(u-u_c)/{\tau }_d,\hfill \mathrm{(4b)}& \hfill J_{so}\left(u\right)& =uH(u_c-u)+\frac{1}{{\tau }_r}H(u-u_c)/{\tau }_0,\hfill \mathrm{(4c)}& \hfill J_{si}(u,w)& =-w(1+tanh\left[k(u-u_c^{si})\right])/\left(2{\tau }_{si}\right).\hfill \end{array}$$

$J_{stim}=J_{stim}\left(t\right)$ represents the external electrical impulse, and $H\left(x\right)$ is the Heaviside step function defined as $H\left(x\right)=1$ for $x\ge 0$ and $H\left(x\right)=0$ for $x<0$. Different parametrizations of the FK model have been proposed [61,62]. In the current study, we refer to [37], with model parameters provided in

Table A3. Parameters used in the FK electrophysiologycal model.

$u_c$	$u_v$	${\overline{g}}_{fi}$	${\tau }_{v1}^{-}$	${\tau }_{v2}^{-}$	${\tau }_v^{+}$	${\tau }_0$	${\tau }_r$	k	${\tau }_{si}$	$u_{\mathit{csi}}$	${\tau }_w^{-}$	${\tau }_w^{+}$	$C_m$
0.13	0.04	4	1250	19.6	3.33	12.5	33.33	10	29	0.85	41	870	1

, Appendix A.

2.2.2. The TenTusscher-Panfilov 2006 Model

The TenTusscher-Panfilov 2006 (TP06) model is a 19-variable physiological model that includes an extensive description of the intracellular calcium dynamics and is good at reproducing cardiac action potentials of human epicardial, endocardial, and mid-myocardial cells [30]. In this case, the sum of ionic currents in Equation (1a) is defined as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill I_{ion}=I_{Na}+I_{K1}+I_{to}+I_{Kr}+I_{Ks}+I_{CaL}+I_{NaCa}+I_{NaK}+I_{pCa}+I_{pK}+I_{bCa}+I_{bNa}+I_{stim},\end{array}\end{array}$$

(5)

where $I_{NaCa}$ is the $N{a}^{+}/C{a}^{2+}$ exchanger current, $I_{NaK}$ the $N{a}^{+}/K^{+}$ pump current, $I_{pCa}$ and $I_{pK}$ are the $C{a}^{2+}$ and $K^{+}$ plateau currents, and $I_{bCa}$ and $I_{nNa}$ are the background $C{a}^{2+}$ and $N{a}^{+}$ currents. The description of each current and related gating variables is provided in Appendix A.1. Model parameters and initial conditions are in

Table A1. Parameter values used during VF simulations for the TP06 model. Parameters not included in this table take the same values reported in [30]. Modified parameters are the maximum conductance of the $I_{Kr}$,$I_{Ks}$, $I_{pCa}$ and $I_{pK}$ currents. The time constant of the f gate was also modified.

Section	$G_{Kr}$	$G_{Ks}$	$G_{pCa}$	$G_{pK}$	${\tau }_f$ Inactivation
Midmyocardium	0.172	0.0515	1.8545	0.00073	$\times 2$
Epicardium	0.172	0.2205	1.8545	0.00073	$\times 2$
Endocardium	0.172	0.2205	1.8545	0.00073	$\times 2$

and

Table A2. Initial conditions used for the TP06 model.

V	${\mathit{X}\mathit{r}}_{\mathbf{1}}$	${\mathit{X}\mathit{r}}_{\mathbf{2}}$	$\mathit{X}\mathit{s}$	m	h	j	d
$-85.23$	$0.00621$	$0.4712$	$0.0095$	$0.00172$	$0.7444$	$0.7045$	$3.373\times {10}^5$
$\mathit{f}$	${\mathit{f}}_{\mathbf{2}}$	${\mathit{f}}_{\mathit{cass}}$	$\mathit{s}$	$\mathit{r}$	${\mathit{Ca}}_{\mathit{i}}$	$\overline{\mathit{R}}$	${\mathit{Ca}}_{\mathit{SR}}$	${\mathit{Ca}}_{\mathit{ss}}$	${\mathit{Na}}_{\mathit{i}}$	${\mathit{K}}_{\mathit{i}}$
$0.7888$	$0.9755$	$0.9953$	$0.999998$	$2.42\times {10}^{-8}$	$0.000126$	$0.9073$	$3.64$	$0.00036$	$8.604$	$136.89$

.

Layer-specific electrophysiological properties are included in this model through the transient outward and slow delayed rectifier currents. In particular, epicardial cells show higher $I_{to}$ current (via $G_{to}$) and higher inactivation recovery (via ${\tau }_s$). In constrast, mid-myocardial cells are characterized by a lower $I_{\mathrm{K}s}$ current.

2.3. Numerical Discretization

The monodomain model presented in Equation (1a) along with the evolution equations for the gating variables were discretized in space using a standard Galerkin finite element (FE) approach, as detailed in [10]. Linear tetrahedral elements were used in the FE discretization. The resulting semi-discrete dynamical system was integrated using a Strang operator-splitting method, as detailed in [63]. Time integration of the gating equations in all simulations was performed using the Rush-Larsen method. All computational implementations were developed in Python, using the Fenics library for the discretization and solution of the monodomain system. Given conduction velocity dependency to the mesh size, a convergence test was carried out, finding that a characteristic mesh size of ≈0.8 mm provided a good approximation to the original conduction velocity at a normal pacing (≈0.42 mm/ms for the FK model and ≈0.67 mm/ms for the TP06 model). This mesh size implied a set of nonlinear equations of about 7 million degrees of freedom.

2.4. Stimulation Protocols and Post-Processing

We computed APD restitution distribution for both ionic models by applying repetitive electrical stimulations at the apex of the bi-ventricular heart models (see S1 location in Figure 2a) with varying pacing frequency. The total stimulation time corresponded to $17\mathrm{s}$ of physical time (32 electrical stimulations) for each case. At every pacing cycle length (CL), the APD was computed locally to obtain a regional distribution of the restitution curves. An empirical probability distribution function (PDF) was estimated using APD data extracted from selected nodes of the finite element discretization highlighted in Figure 2b. Activation and repolarization times were obtained from unrefined mesh data and interpolated over the finite element discretization.

We further implemented an S1–S2 stimulation protocol to induce a sustainable VF within each bi-ventricular heart geometry. The S1 stimulation site was located on the LV endocardium, and the S2 stimulus was delivered on the posterior epicardium layer in both models. We explored various S1–S2 stimulation locations, performing a preliminary sensitivity analysis combined with stimulus strength, finding no differences in VF dynamics when stimuli locations were varied. The results analyzed in the following are referred to the cases requiring minimum stimulation current to induce sustained fibrillation, i.e., the most high-risk scenarios. This is supported by previous research on the effects of the initial site on VF formation in dog and human heart geometries concluding that, although the stimulation site may affect VF onset, the number of filaments (e.g., VF signatures) reaches the same steady-state value [26,50] such the underlying physics is preserved [64].

We calculated the vulnerability window for each electrophysiological model by measuring the elapsed time between the S1 and S2 pulses for cases that resulted in sustained ventricular fibrillation (i.e., VF duration > 2 s). To further characterize VF dynamics, we computed evolution of scroll wave filaments during $10\mathrm{s}$ of physical time with a time step of 10 ms. After identification of the intersection between an isopotential line ($-70\mathrm{mV}$) and the constraint $d{V}_m/dt=0$, each intersection was related with a finite element and filaments were labeled and counted using a density-based clustering algorithm.

To assess the global behavior of cardiac simulations, pseudo-ECGs were computed to study VF Signatures such as the fundamental frequency for all of the ionic models analyzed. To this end, surface potential was estimated during numerical simulations using the classical approximation

$$\begin{array}{c}\hfill \begin{array}{c}\hfill V_e=-{\int }_{\mathsf{\Omega }}\nabla V_m·\nabla \frac{1}{\left|\right|\rho \left|\right|}d\mathsf{\Omega },\mathrm{with}\rho =|{\mathbf{x}}_e-\mathbf{x}|.\end{array}\end{array}$$

(6)

where $\rho \in {\mathbb{R}}^3$ represents the distance from each point of the heart domain to the virtual external electrode, located $2\mathrm{cm}$ away from the left ventricular wall. This setting is commonly used to mimic precordial leads [65,66]. After the pseudo-ECG was constructed, the fundamental frequency was identified by building up the power spectra of the signal [42].

3. Results

3.1. Dispersion of APD Restitution

APD restitution curves were computed for the two electrophysiolgical models over four selected cardiac surfaces (see Figure 3a): epicardium (EPI), left ventricular endocardium (LV), right ventricular endocardium (RV), left ventricular mid-myocardium (LVMM). Accordingly, a numerical probability density function (PDF) of the action potential duration was derived for each anatomical region both for NC and HFC geometry models. The eight distributions of restitution curves computed are shown in Figure 3b,c. For comparison purposes, APDs were normalized according to the formula ${\mathrm{APD}}_n=(\mathrm{APD}-{\mathrm{APD}}_{min})/({\mathrm{APD}}_{max}-{\mathrm{APD}}_{min})$, where, ${\mathrm{APD}}_{min}$ and ${\mathrm{APD}}_{max}$ are provided in

Table 1. Normalization values for FK and TP06 models.

	FK	TP06
${\mathrm{APD}}_{min}$	250 $\mathrm{ms}$	250 $\mathrm{ms}$
${\mathrm{APD}}_{max}^{EPI}$	273 $\mathrm{ms}$	335 $\mathrm{ms}$
${\mathrm{APD}}_{max}^{LV}$	273 $\mathrm{ms}$	350 $\mathrm{ms}$
${\mathrm{APD}}_{max}^{RV}$	273 $\mathrm{ms}$	350 $\mathrm{ms}$
${\mathrm{APD}}_{max}^{LVMM}$	273 $\mathrm{ms}$	390 $\mathrm{ms}$

for each anatomical surface.

In the FK case, ${\mathrm{APD}}_n$ distributions stabilize around $0.78\pm 0.07$ for $\mathrm{CL}\ge 500\mathrm{ms}$, whereas, in the TP06 case, PDF restitutions show different average values and a non-uniform variance. In addition, a steeper and narrower distribution is observed for the FK model for $\mathrm{CL}<400\mathrm{ms}$, whereas a gradual reduction is noted for TP06. Interestingly, the FK distribution display peaks at the center of the distributions that remain stable during the whole restitution curve. The four anatomical surfaces show similar trends in both healthy and diseased conditions. Unlike FK, TP06 presents a strong dependence of the distribution peaks on the selected region. In particular, LV and RV surfaces show severe differences among NC and HFC models. Moreover, the left ventricle shows the highest dispersion of repolarization, i.e., higher variance, in the pathological case. Median and standard deviations of the distributions are provided in

Table 2. Medians and standard deviations for FK and TP06 model, corresponding to Figure 4 ($\mathrm{CL}=400\mathrm{ms}$).

	NC		HFC
	Median	SD	Median	SD
FK Model
EPI	0.41	±0.07	0.39	±0.08
LV	0.41	±0.08	0.41	±0.06
RV	0.39	±0.07	0.41	±0.07
LVMM	0.41	±0.05	0.41	±0.06
TP06 Model
EPI	0.27	± 0.08	0.26	±0.07
LV	0.29	±0.08	0.36	±0.11
RV	0.37	±0.10	0.23	±0.06
LVMM	0.41	±0.05	0.43	±0.06

for $\mathrm{CL}=400\mathrm{ms}$. Using 62 CPUs within a parallel computing platform, the restitution protocol simulation took around 55 h of wall-clock to compute the FK model and 82 h for the TP06 model.

To better quantify the previous observations, a detailed overview of the computed numerical PDFs is shown in Figure 4 and Table 2 at CL = 400 ms. Remarkable discrepancies in density distributions between the two electrophysiological models can be seen among the NC and HFC geometries. In particular, the FK model displays larger median values and narrower distribution profiles in both healthy and disease cases (Figure 4a). Medians and standard deviations do not vary significantly at different myocardial layers. The TP06 model shows bimodal distributions on the two endocardial surfaces (LV and RV) and the epicardial region, resulting in smaller median values and larger variance (Figure 4b). On the LVMM surface, the PDF median values are similar among the two ionic models, and negligible differences arise between healthy and pathological cases for the TP06 model. In particular, the expected median value is ≈0.41, and the standard deviation is $\pm 0.05$. In Appendix A.3 we provide the same statistical analysis at $\mathrm{CL}=500\mathrm{ms}$ and $\mathrm{CL}=600\mathrm{ms}$ further confirming these features.

3.2. VF Sustainability

We studied cardiac arrhythmia signatures by implementing an S1–S2 stimulation protocol to induce sustained VF, i.e., multiple scroll wave formations. Figure 5 shows a qualitative visual assessment of the scroll wave evolution for the two ionic models within the HFC geometry. Three selected times of transmembrane activation, $V_m>-75$ mV, are compared where the intramural transparency of the excitation highlights both the complexity of the dynamics and the differences among the two ionic models (Figure 5a for FK and Figure 5b for TP06). The ratio of the amplitude for the two S1–S2 stimuli differed. The FK model had an equal 1:1 ratio with a vulnerable window of $345\mathrm{ms}$ in both NC and HFC subjects. In contrast, the TP06 model required a much higher stimulation amplitude ratio of 11:2 for the NC case and 15:2 for the HFC case, with a higher vulnerable window of $408\mathrm{ms}$.

To provide a quantitative indication of VF dynamics and sustainability, we assessed the time evolution of scroll waves by identifying the total number of 3D filaments during 10 s of physical time (see Figure 6b). Employing the same computational settings adopted for the restitution protocol, the VF wall-clock time was 34 h for the FK and 48 h for the TP06 model.

The two ionic models are similar in terms of rolling averages computed with a time window of $500\mathrm{ms}$. In all cases, the number of filaments stabilizes after about 2.5 s from arrhythmia induction (see Figure 6c). Mean value of filaments, mean and variance, and maximum and minimum values are compared in

Table 3. Rolling averages of the number of filaments: mean and variance, maximum and minimum values after 2.5 s from arrhythmia induction. Comparison among ionic (FK and TP06) and geometry (NC, HFC) models.

	NC			HFC
	Mean	Max	Min	Mean	Max	Min
FK	25 ± 5	44	10	33 ± 6	52	15
TP06	27 ± 6	46	11	34 ± 8	60	14

. We computed the $L_1$ error between the number of filaments since it is less sensitive to extreme values compared with the $L_2$ error. The $L_1$ error between the mean number of filaments in the FK and TP06 models was 14.2% and 11.4% for NC and HFC, respectively.

Finally, the quantification of pseudo-ECG provides the time course of the computed signal for both models in healthy (Figure 7a) and diseased (Figure 7b) heart geometries. The resulting fundamental frequencies were about 11 Hz for FK and 5 Hz for TP06 independently of the anatomical model considered.

4. Discussion

In this work, we investigated how the choice of ionic models affects the repolarization properties and VF signatures of computational models of normal and failing hearts. Numerical analyses allowed us to build up probability distribution functions of APD restitution curves for four selected anatomical surfaces. Accordingly, we identified left and right ventricular endocardium as the anatomical regions characterized by a higher degree of dispersion of repolarization (multimodal PDFs) for the myocardial infarction case. On the other hand, strong electrotonic effects within the mid-myocardial layer homogenize the resulting electrical behavior, thus providing coherent PDFs among the two ionic models for healthy and diseased cases, although non-negligible differences appear for epicardium, left and right endocardium.

This result, consistent with the literature and specific kinetics described by the TP06 ionic model, was further confirmed by long-run simulations conducted for ventricular fibrillation. The statistical examination of the evolution and the number of vortex filaments provided an enhanced arrhythmogenicity in the diseased case. Constitutive differences among the two electrophysiological models appear in terms of restitution dynamics at large cycle lengths, stimulus intensity, and vulnerable windows. Such a discrepancy is due to the intrinsic complexity of the ionic models in delivering an external stimulation current signal among the different gating variables. The membrane capacitance and the surface-to-volume coefficients regulate the electrical current’s spreading, with a predominant role in ionic modeling cases. Therefore, we highlight the necessity of fine-tuning these parameters when using a physiological approach other than a phenomenological one. However, dissimalirities were negligible for ventricular fibrillation signature. Specifically, irrespective of the healthy or pathological geometry model, the fundamental frequency of the pseudo-ECG during sustained VF is constant. Moreover, the number of filaments in the HFC case was statistically equal among FK and TP06 models.

Accordingly, our results suggest that VF signatures are mainly controlled by anatomical and structural features rather than regional restitution properties. Therefore, given the large-scale computational models and clinical translation, the choice of a simplified phenomenological description seems favorable. From a computational point of view, phenomenological models allow the implementation of highly efficient and simple integration schemes, whereas biophysical-based models, such as the TP06, require specialized and adequate time integration schemes that demand a complex implementation. Several studies have tackled the development of time-integration strategies to cope with the stiff nature of these models [34,67]. Besides the high computational demand for solving the electrophysiology equations, the used programming language (Python) is limited by its lack of performance and scalability. Moreover, more efficient numerical approaches, like spatio-temporal adaptivity or strong scalable solvers, are also needed to cope with the computational burden imposed by high-resolution electrophysiological problems. Future research may use a programming language better suited in terms of scalability and performance, such as C++, which may be used in conjunction with Python and highly advanced numerical techniques [60,68]. Despite this, we solved a highly computational demanding problem in a relatively modest amount of time. Moreover, several parameters required for postprocessing (i.e., pseudo-ECG, APD distribution, and filaments) were computed during simulation time, saving time for analysis.

The current work can be extended in several directions. First, a general comparison between phenomenological and biophysical models requires the investigation of the bidomain theoretical framework to emphasize the role of structural material heterogeneities during defibrillation procedures [69,70,71]. In this direction, a proper mathematical and computational modeling of torso will be critical for obtaining clinically relevant ECG signals, i.e., acceleration of forward ECG numerical modeling with the use of personalized torso geometries [72]. Second, comparing linear, nonlinear, and fractional diffusion formulations would significantly help to identify the best computational approach to apply in subject-specific cardiac studies [10,73,74,75]. Third, characterization of the border zone of the scar and gray zones in the infarcted myocardium is related to small-scale components [9,76] that are accurately modeled by detailed biophysical constitutive models reproducing modified ion channel dynamics (electrical remodeling). However, reduced-order models have been shown to provide a viable alternative also in this context [77]. This is key in clinical-related applications such as drug delivery and stem cell delivery in which disease characterization and treatment depend on the electrochemical dynamics of the cell membrane and the associated local heterogeneities.