In-silico study of the cardiac arrhythmogenic potential of biomaterial injection therapy

Jul 31, 2020

Geometrical and morphological representation of swine hearts

This study considered one normal control heart (NC), one heart control with ischemic heart failure without biomaterial treatment (HFC), and one heart with ischemic HF and biomaterial injections (HFI). These subjects were selected as representative of larger cohorts considered in a previous morphological study of the effects of fiber remodeling under biogel treatment6. The three selected subjects included natural geometrical differences due to biological variance and remodeling to ischemia. This has been quantified in previous studies which showed local wall thinning and fiber reorientation due to ischemia, but that the global structure and morphology of the hearts were not significantly different6. The injection protocol for the hydrogel delivers a total of 12–14 intra-myocardial injections (0.3 mL each) in a circumferential pattern into the LV free wall during an open chest procedure. These are administered in roughly 1.5 cm apart and in two rows: one above and one below the mid-ventricular plane between the base and the apex. The hydrogel, which accounts for roughly 3% of the total LV wall volume, does not disrupt fiber orientation, but rather conforms to the native structure it is injected into, forming ellipsoidal shapes orientated with the local fiber structure6. The solidified injections mitigate the effects of adverse remodeling by anchoring and supporting the surrounding tissue in the nearby vicinity of the circumferential pattern. A deeper analysis of therapy efficacy and the mechanism of action in mechanical terms was recently provided in a recent study7.

Subject-specific accurate geometric representations of heart bi-ventricular structure, infarcted tissue, and biomaterial injections are used as the computational domain for the numerical simulations (Figs. 7 and 8). Imaging data originates from ex vivo segmentation of high-resolution MRI and DT-MRI of swine hearts. The experimental protocol, image acquisition, segmentation process and reconstruction methodology have been described previously8,55. In brief, myocardial infarction was induced by occluding the obtuse marginal branches of the left circumflex artery. Eight weeks after MI, animals underwent Algisyl-LVR injection, and hearts were excised eight weeks later. Anatomical MRI and DT-MRI were acquired using a readout-segmented diffusion-weighted spin-echo sequence with (1.0 times 1.0 times 1.0 ,hbox { mm}^{3}) spatial resolution. We discretized the heart domains using tetrahedral finite elements, identifying healthy, infarcted and hydrogel regions based on MRI observations (Fig. 7). Following previous works55, local properties for the infarcted and healthy tissues are modulated by the volume fraction of healthy tissue. This volume fraction is represented by the space-dependent scalar function (h=h(mathbf{x}) in [0,1]), where (mathbf{x}) represents the Cartesian coordinate vector. In particular, (h(mathbf{x})=0) defines properties of the infarcted zone (IZ), (h(mathbf{x})=1) identifies healthy tissue, and (0<h(mathbf{x})<1) defines the transition zone or gray zone (GZ), where mixed electrical properties of infarcted and healthy tissue are modeled according to the literature20,43. A gray zone ratio (GZR) could also be computed from this function to characterize the differences in GZ distribution between each heart under study. Myocardial fiber orientation, based on DT-MRI, was assigned to the mesh nodes and interpolated inside each finite element, delivering a continuous spatial vector field representation of the cardiac fiber orientation (mathbf {f}=mathbf {f}(mathbf {x})) for each heart analyzed55. To account for transmural dispersion of repolarization, we used a Laplace’s interpolation method56 to divide the heart wall into epicardial, mid-myocardial, and endocardial layers using a thickness ratio of 2:3:3, respectively. Figure 9 shows these transmural layers for the three hearts analyzed. The LVMM surface used in the study of APD distribution was constructed from the Laplace’s interpolation by creating a mesh at the mid-surface of the mid-myocardial region.

Numerical modeling of cardiac electrical activity

Given the cardiac domain (Omega in {mathbb {R}}^3) composed by conductive tissue, the propagation of electrical impulses within a time interval [0, T] was modeled using the monodomain electrophysiology equation derived from current balance57, which takes the form

begin{aligned} dfrac{partial V_m}{partial t}&= nabla cdot left( mathbf{D} , nabla V_m right) + I_s(V_m,mathbf{r}) quad text {in} quad {Omega times [0,T]} ,, end{aligned}

(1a)

begin{aligned} dfrac{dmathbf{r}}{d t}&= mathbf{g}(V_m,mathbf{r}),, end{aligned}

(1b)

where (V_m:Omega times [0,T] longrightarrow {mathbb {R}}) represents the transmembrane potential, defined as the difference between intracellular and extracellular voltages, (I_s) is the sum of the ionic currents depending on the transmembrane potential and (mathbf{r}:Omega times [0,T] longrightarrow {mathbb {R}}^{m}) is the set of state variables. This monodomain framework is constructed from a standard bidomain model58 by assuming that both intracellular and extracellular spaces have equal anisotropy ratios. In an averaged macroscopic sense, the conduction velocity is faster along the fiber direction and slower in a direction transverse to it. In particular, the conductivity tensor ({mathbf{D}} in {mathbb {R}}^{N times N}_{+}) is assumed to be transversely isotropic, which can be expressed as

begin{aligned} mathbf{D} = h(mathbf{x}) left[ d_{perp } mathbf{I} + left( d_{perp }-d_{parallel }right) mathbf{f} otimes mathbf{f}right] ,, end{aligned}

(2)

where (d_{perp },d_{parallel }) represent the transversal and longitudinal conductivities, respectively, and (mathbf{f}=mathbf{f}(mathbf{x})) is the fiber direction vector field. The conductivity ratio (d_{parallel }/d_{perp }) varies typically between 4 and 959. We selected a conductivity ratio of 4, establishing the longitudinal conductivity as (d_{parallel }=0.1;mathrm{mm}^2/mathrm{ms}), with a corresponding transversal conductivity of (d_{perp }=0.025;mathrm{mm}^2/mathrm{ms}). To complement equations (1a) and (1b), we assume the no-flux boundary condition by denoting the boundary normal by (mathbf{n}) and enforcing that

begin{aligned} nabla V_m cdot mathbf{n} = 0 quad text {on} quad partial Omega , end{aligned}

(3)

which includes the internal interfaces between cardiac tissue and gel injection, as the gel is initially assumed non-conductive. For the simulations where the biogel is assumed to be electrically conductive, we consider (Omega _{gel}) as the domain of injectates, and represent the electrical conduction in that region by the Laplace equation

begin{aligned} nabla cdot left{ mathbf{D}_{gel}nabla V_m right} = 0 quad text {in} quad Omega _{gel}. end{aligned}

(4)

In writing Eq. (4) we have considered (-V_m) as the voltage field in the biogel. This consideration can be justified by assuming that the biogel domain constitutes an extracellular domain, and that the voltage in the intracellular domain is zero in (Omega _{gel}), allowing the use of the transmembrane potential field in the simulation of electrical conduction in the biogel domain.

The transmembrane ionic current is modeled according to the model proposed by ten Tusscher and Panfilov (TP06)11 for ventricular human cells. This model is fitted to reproduce APD restitution curves measured in humans, and have an extensive description of the intracellular calcium dynamics. It has been broadly used for the study of reentrant arrhythmia and electrical instability at cellular and tissue levels. In this model the ionic current density, defined as the sum of all transmembrane current densities, is given by

begin{aligned} I_s = I_{Na}+I_{K1}+I_{to}+I_{Kr}+I_{Ks}+I_{CaL}+I_{NaCa}+I_{NaK} +I_{pCa}+I_{pK}+I_{bCa}+I_{bNa}, end{aligned}

(5)

where (I_{NaCa}) is (Na^+/Ca^{2+}) exchanger current, (I_{NaK}) is (Na^+/K^+) pump current, (I_{pCa}) and (I_{pK}) are plateau (Ca^{2+}) and (K^+) currents, and (I_{bCa}) and (I_{nNa}) are background (Ca^{2+}) and (Na^+) currents. One advantage of using this model is that, by changing the parameter values it allows for the representation of transmural heterogeneity observed in myocardial tissue. To this end, the cardiac wall was divided into epicardial, mid-myocardial and endocardial layers using the distribution ratio 2:3:3 as adopted in previous works60, see Fig. 9. The set of parameters corresponding to each layer, which were used in this work for the simulation of the restitution protocol, have been reported elsewhere11. For the simulation of VF, the parameter values reported in21 were employed, which are included in Table 1. Initial values for the gating and internal variables are included in Table 2.

The governing equations described in Eq. (1a) are discretized in space using a standard Galerkin finite-element scheme61, where linear tetrahedral elements are employed to approximate the transmembrane potential field. Time integration was performed using an operator splitting method58. All numerical implementations were developed in Python, using the FENICS Project software and the Cbcbeat software collection in an in-house parallel computing platform. Given the well-known dependency of the conduction velocity to the mesh size, conduction velocity convergence was carried out for linear tetrahedral elements (Fig. 10). In particular, we found that a characteristic mesh size of (approx 0.8;mathrm{mm}) provides a good approximation of the sought conduction velocity ((approx 0.67,mathrm{mm}/mathrm{ms})) at a normal pacing. This implies a set of nonlinear equations with over 7 million degrees of freedom. Time integration of gating evolution equations at quadrature points is performed using the explicit Rush-Larsen method with a time step of (dt=0.1;mathrm{ms}). Overall, 1 second of simulation required about 4.8 hours of computation using 62 CPU’s within a parallel computing platform AMD Opteron Processor 6378. For the APD restitution estimation, this meant a total of 81.6 hours for each heart, whereas VF simulations involved around 48.0 hours of computation for each heart.

We developed in-house routines to post-process the simulation results to compute parameters such as activation (AT) and repolarization times (RT), during numerical simulations with a time step (dt^*=5dt) over selected nodes in the finite element mesh (Fig. 10). The nodes were selected using information of unrefined meshes. This setup reduced computational costs, especially for post-processing purposes such as the APD restitution curve. Given the complex geometry of each model, we also performed region labeling for the different heart surfaces by using a Laplace interpolation approach. In particular, we computed APD restitution curves in the LV endocardium (LV), RV endocardium (RV), left ventricle mid-myocardium (LVMM), and epicardium (EPI) on a simple subdivision of the heart geometric features directly during simulations. Accordingly, we produced a probability density function (PDF) using a Gaussian kernel density approach to analyze the distribution of the APD restitution in these regions.

Simulation protocols and analysis

To compute the APD restitution distribution, the three heart models were paced at the apex location (see star in Fig. 1) 32 times with varying pacing frequency. The total simulation time reached 17 s of physical time for each geometry model, with stimulus value of (40, mathrm{mV}/mathrm{ms}) and duration of (2,)ms. At each pacing cycle length (CL), the APD was computed locally to obtain a regional distribution of the restitution curves. In particular, APD was computed for each selected node (see Fig. 10) of the finite-element discretization and the PDF was estimated from this data (see Fig. 2). For comparison purposes, the PDF was examined within the EPI, LV, RV and LVMM surfaces.

An S1–S2 protocol was implemented to induce VF in each heart model. The S1 stimulation site was the septum LV endocardium, while the S2 stimulus was delivered at the posterior zone of the epicardium, near the tail of the S1 wave. The timing between S1 and S2 was fine-tuned to obtain excitation wave-break, formation of a scroll-wave and evolution into VF. A vulnerable window was measured as the elapsed time between the S1 and S2 pulses for which a sustained VF was achieved ((> 2) s). Depending on the heart model used, different values for the S2 stimulus were needed. The S2 stimulus value was (300, mathrm{mV}/mathrm{ms}) (7.5 times the S1 stimulus value) for the NC heart, while for both HFC and HFI hearts the S2 stimulus value was around (220, mathrm{mV}/mathrm{ms}) (5.5 times the S1 stimulus value). The vulnerable window was 408 ms, 406 ms and 423 ms for the NC, HFC and HFI heart, respectively. The complexity of activation patterns developed during VF dynamics was quantified by computing the number of rotors in time. These rotors are 3D scroll waves that rotate around a filament24,62. Scroll wave filaments were distinguished using the algorithm proposed by Fenton and Karma59. In particular, a singular point is found by computing the intersection of an isopotential line (-70 mV) with the condition (dV_m/dt =0). In the computational model here implemented, each singular point is related to a single finite element. Afterwards, scroll wave singular line, e.g. filament, was identified and labeled by using a density-based spatial clustering algorithm (DBSCAN)63. This method allows one to group elements that are closely packed together forming a specific filament. Each group (filament) formed is classified, updated and counted every (10,)ms of physical time directly during simulations.

The electrical properties of alginate hydrogel implants have not been reported to date, thus in the previous experiments we have assumed that its conductivity is zero. To assess the effect of gel conductivity on VF sustainability in the treated heart, we performed a sensitivity analysis where the biogel is assumed to behave as a passive conductor. VF simulations were performed assuming an isotropic conductivity in the gel of the form

begin{aligned} mathbf{D}_{gel} = c d_{perp } quad text {in} quad Omega _{gel}, end{aligned}

(6)

where c is a conductivity ratio that modulates the conductivity of the gel based on the value used for transversely-isotropic propagation in normal cardiac tissue. Accordingly, for values of the conductivity ratio (c<1) the biomaterial region is less conductive than the normal tissue and for values of conductivity ratio (c>1) the biomaterial region is more conductive than the normal tissue. We performed simulations for (c = {0.0, 0.5, 1.0, 1.5}), and computed the number of rotors in each case to evaluate the spatiotemporal dynamics during VF. pseudo-electrocardiograms (EGCs) were computed by estimating the surface potential using the approximation64

begin{aligned} V_e= – int _{Omega } nabla V_m cdot nabla frac{1}{||rho ||} dOmega ,quad text {with} quad rho =| mathbf{x}_e- mathbf{x} |. end{aligned}

(7)

Here, (rho in {mathbb {R}}^3) defines the distance from each point of the heart to a point placed at some distance from the ventricles (i. e., electrode location). This position was established at 2 cm away from the left ventricular wall, as is commonly defined to mimic pre-cordial leads to study T waves and QT intervals65. ECG signals were analyzed using Fourier Transform, from which power spectra were constructed, and the fundamental frequency was identified.