Drift of Scroll Waves in a Mathematical Model of a Heterogeneous Human Heart Left Ventricle

: Rotating spiral waves of electrical excitation underlie many dangerous cardiac arrhythmias. The heterogeneity of myocardium is one of the factors that affects the dynamics of such waves. In this paper, we present results of our simulations for scroll wave dynamics in a heterogeneous model of the human left ventricle with analytical anatomically based representation of the geometry and anisotropy. We used a set of 18 coupled differential equations developed by ten Tusscher and Panﬁlov (TP06 model) which describes human ventricular cells based on their measured biophysical properties. We found that apicobasal heterogeneity dramatically changes the scroll wave dynamics. In the homogeneous model, the scroll wave annihilates at the base, but the moderate heterogeneity causes the wave to move to the apex and then continuously rotates around it. The rotation speed increased with the degree of the heterogeneity. However, for large heterogeneity, we observed formation of additional wavebreaks and the onset of complex spatio-temporal patterns. Transmural heterogeneity did not change the dynamics and decreased the lifetime of the scroll wave with an increase in heterogeneity. Results of our numerical experiments show that the apex may be a preferable location of the scroll wave, which may be important for development of clinical interventions.


Introduction
Vortices in excitable medium, which are called spiral waves in 2D or scroll waves in 3D, were found in many physical, chemical, and biological systems [1]. They play an important role in dynamics of these systems. For example, formation of such vortices in the heart causes cardiac arrhythmias, which remain the largest cause of death in the industrialized countries [2]. The dynamics of such vortex in the heart determines the type of cardiac arrhythmia. For example, the drift of a scroll wave in the ventricles of the heart may result in the onset of an arrhythmia called polymorphic ventricular tachycardia [3] and a breakup of the vortex into complex spatio-temporal patterns results in the onset of ventricular fibrillation, which causes cardiac arrest and sudden cardiac death [4]. Therefore, the factors responsible for drift of scroll waves are of great interest.
In a very general sense, the factors which can cause drift of a scroll wave include geometrical factors, anisotropy, and tissue heterogeneity. Previous research on generic models showed that, in homogeneous isotropic medium with high excitability, scroll waves drift to the regions where the

Baseline Homogeneous Model. Numerical Approach and Software
The numerical computation of electrophysiological activity requires models at three levels: cell, tissue, and organ.
We used the LV anatomical model described in [24]. This model is axisymmetric and represents an average ventricle of a healthy adult human. The LV model was represented as a body of revolution with the shape fitted to experimental data [30]. The rotation axis is vertical axis Oz. An example of the LV model is shown in Figure 1C. The epicardium is the colored surface, and the endocardium is the mesh surface. Figure 1. Distribution of APD 90 (shown by color, from blue, 353 ms, to red, 248 ms) in the considered heterogeneous models of the LV. (A) apicobasal heterogeneity is set by a decrease of IKs at the apex (shades of blue) or ICaL at the base (shades of red). APD90 increases with increase of Ψ from 0 (base) to π/2 (apex); (B) Transmural heterogeneity is set by decrease of IKs at the subendocardium. γ = 0.1 means subepicardium cells, γ = 0.9, subendocardial; (C) points with maximal difference of APD90 in cases of apicobasal (Ψ = 0, Ψ = π/2) and transmural (γ = 0.1, γ = 0.9) heterogeneity. Color shows the z-coordinate; (D) plot of the action potential against time in LV points with maximal heterogeneity.
The model has the following list of parameters. R b is the LV outer radius on the LV equator; Z b is the LV height including the apical wall thickness; L is the LV wall thickness at the LV equator; h is the LV wall thickness at the LV apex; and ∈ [0, 1] is a dimensionless parameter influencing the conicity-ellipticity characteristic of the LV shape.
The local coordinates are linked with the cylindrical coordinates (ρ, ϕ, z) by the following formulae: In every LV point, a vector of myofibre direction is defined as described in [24]. The parameters of the anatomical model were kept constant between all our simulations: outer (epicardial) radius at the base R b = 33 mm, basal thickness L = 12 mm, full height Z b = 60 mm, and ellipticity coefficient ε = 0.85.
A set of parameters was taken as reference. The reference ventricle had anisotropy coefficients in ratio 9:1, transmural fibre rotation angle (FRA) of 147 • , and apical LV thickness h = 6 mm.
The LV model allows us to easily change LV geometrical parameters-for example, to vary its thickness at the base and apex. These parameters are important for studying scroll waves drift which goes, accordingly with scroll waves theory, to a thin place if filament tension is positive.
For the tissue level, we used a monodomain approach [31] and anisotropic medium. Reaction-diffusion equations describe processes in the modelled body and have the following form: I ion = I Kr + I Ks + I K1 + I to + I Na + I bNa + I CaL + I bCa + I NaK + I NaCa + I pCa + I pK .
Here, u = u(r, t) is the transmembrane potential; the intracellular processes are captured by I ion = I ion (r, t), which is the sum of the ionic transmembrane currents. C m is the capacitance of the cell membrane. In a general sense, this equation shows that the potential difference V m over the cell membrane changes due to currents (I) which flow through the membrane. These currents are conveyed by different ions (Na, K, Ca) and via various biophysical processes. In modern computational cardiac electrophysiology, the properties of all such currents are fitted to their experimentally measured values and their dynamics is fitted using additional differential equations. Each of the currents typically depends on V m and time, and the time-dependency is normally given by an exponential relaxation equation for the gating variable(s) g. For example, consider a hypothetical current I * which conveys ion ' * ': It has a maximal conductivity of G * = const and is zero for V * , the so-called Nernst potential for ion ' * '. This Nernst potential can be easily computed from the concentrations of specific ions outside and inside the cell membrane. The time dynamics of this current is given by two gating variables g * , g • to the power α, β. The variables g * , g • approach their voltage-dependent steady state values g ∞ i (V m ) with characteristic time τ i (V m ) (4). All parameters and functions are chosen here to fit experimentally measured properties of the specific ionic current. Most of ionic currents have one or two gating variables, with α = β = 1. In our simulations we use equations derived and fitted in [32], which has 18 state variables. Note that a proper description of ionic currents is very important as most cardiac drugs affect the maximal conductivities of specific ion channels. Furthermore, cardiac tissue is heterogeneous with respect to expression levels of these channels. Thus, using Equations (2)-(4), one can model the effect of pharmacological preparations on cardiac tissue, and, as we do in this paper, study effects of heterogeneity of the heart on cardiac arrhythmias. One of the standard ways to do so is to make the maximal conductivity of a specific current dependent on space in a manner based on experimental data. As described in the next subsection, the most important currents accounting for the heterogeneity are I CaL and I Ks .
As in [33], the diffusion matrix D = (D ij ) was computed from the unit vectors in fibre direction v using the formula where D 1 and D 2 are the diffusion coefficients along and across the fibres and δ i,j is the Kronecker symbol. Ionic model parameters were taken from [32]; they correspond to physiological characteristics of cardiomyocytes.
A scroll wave was initiated using a standard protocol S1S2 [29] and rotated counterclockwise. We used a standard boundary condition -zero flux of potential through the LV surface. A uniform grid with spatial step dr = 0.28 mm was set in Cartesian coordinates.
To solve the differential Equations (1)-(3), we used a finite difference approach. To approximate the diffusion term, we used a stencil of 18 points for 3D using the following equation: where l is an index over the 18 neighbors of the point (i, j, k) including itself, and w l is the weight of the voltage of a particular neighboring grid point which accounts for its contribution to the Laplacian. The weights w l were computed based on the conduction tensor D ij and location of the point with respect to the local boundary. All points outside the heart geometry have weights w l = 0. The weights were precomputed for each geometry used and allowed an efficient evaluation of the Laplacian during the simulations. The gating variables in the TP06 model were integrated using the Rush and Larsen approach [34]. To speed up calculations, voltage-dependent functions in gating variables equations were precomputed and placed into look-up tables. Additional details on the integration procedure can be found, for example, in [33]. We studied scroll wave dynamics depending on LV wall thickness at the apex and degree of heterogeneity. In our simulations, LV basal thickness was L = 12 mm, while apical thickness was varied to h = 6, 12, 18 mm. We considered the physiological value of the anisotropy ratio D 1 /D 2 = 9.
Our procedure of postprocessing filament coordinates was thoroughly described in [29]. All calculations were performed on a C program on clusters "URAN" (IMM, Ural Branch of RAS) and "UrFU" of Ural Federal University (Ekaterinburg). The program uses CUDA for GPU parallelization and was compiled with a Nvidia C Compiler "nvcc". Computational nodes have graphical cards Tesla K40m0.

Heterogeneity Representation
The heterogeneity in the ventricular myocardium is a result of changes of balance between inward currents (mainly ICaL) and outward potassium currents. Although in reality these currents are changing in combination, the exact degree of modifiication is not fully quantified. In order to compare effects of the heterogeneities in the same conditions, we decided to study a separate effect of ICaL and one of the most important potassium currents IKs.
We considered LV models with the following types of apicobasal heterogeneity: 1. Apicobasal heterogeneity caused by decrease of ICaL current with APD base < APD apex . We reduced ICaL current at the LV base to 75%, 50%, and 25% of its original value, which resulted in the gradients of 14, 28, and 38 ms. We denote these cases as ICaL-75, ICaL-50 and ICaL-25.

2.
Apicobasal heterogeneity caused by decrease of IKs current with APD base < APD apex . We reduced IKs at the apex to 75%, 50%, and 25% of its original value, which resulted in the APD gradients of 14, 34, and 68 ms. We denote these cases as IKs-75, IKs-50, and IKs-25.
To represent transmural (TM) heterogeneity, we decreased the IK current because IK is the main current affecting APD, and it is responsible for transmural heterogeneity [35]. We followed the work [36] and made models with APD epi < APD endo by multiplying IKs conductivity by 75%, 50%, and 25%, which resulted in the APD gradients of 12, 31, and 58 ms. We denote these cases as IKs-75-TM, IKs-50-TM, and IKs-25-TM.
The effects of the changes of the currents on AP shape are shown in Figure 1D. Figure 1A shows distribution of APD from apex to base for each type of apicobasal heterogeneity and Figure 1B shows the three types of transmural heterogeneity.

Results
Below, we present results of our simulations for scroll wave dynamics in ventricles of different geometries and various types of heterogeneity.
We considered three geometries with the same anisotropy but different in the thickness of the apex h = 6, h = 12, and h = 18 mm. In all these geometries in homogeneous ventricles, the scroll wave drifted to the base and disappeared. We studied the dynamics of the scroll waves with the same initial conditions, but in the presence of the heterogeneity.
We considered two types of regional heterogeneity: apicobasal and transmural. Figure 2A,B show a typical trajectory of the scroll wave drift in the absence of heterogeneity. The scroll was initially located between apex and base, then drifted upwards and disappeared after collision with the boundary at the base. We found that the presence of apicobasal heterogeneity substantially changed the scroll wave dynamics. Figure 2C shows an example of scroll wave drift from the same initial conditions, but for a small apicobasal heterogeneity of 14 ms, IKs-75. We see that the scroll wave drifts to the apex, stabilizes at 8 mm from it, and continues its circumferential rotation with a speed of 0.19 mm/s. We performed 18 computations for six types of the heterogeneity in three LVs of different geometry. The results are presented in Figure 3, where we mark the type of observed dynamics (e.g., annihillaiton at the base and breakup) and in case of stabilisation at attractor we show the coordinates of the attractor, and the drift speed along it.

Apicobasal Heterogeneity
Our findings were the following. From the 18 cases studied, in only one case (IKs-75 and h = 18 mm), we observed the same dynamics as in the homogeneous case: annihilation at the base.
In 12 out of 18 cases, we observed a dramatic change of dynamics: instead of annihilation at the base, the scroll wave stabilized at the apex and then continued moving around it along the circle. One example of such dynamics was already shown in Figure 2C, where the filament stabilized at 8 mm from the apex and rotated with a speed of 0.19 mm/s. The apex-base heterogeneity in that case was 14 ms (the corresponding grey columns in Figure 3A). For the same geometry of h = 6 mm and a larger heterogeneity (IKs-25, DR AB = 68 ms), the scroll wave filament stabilized also at approximately 8 mm from the apex (ψ = 1.24) and was moving along the circle with speed 2.27 mm/s. The speed of the motion along the attractor increased with the increase in the heterogeneity. Note that the location of attractor in most of the cases did not depend on the degree of heterogeneity. However, for h = 18 mm, we saw that the location of the attractor for ICaL-50 was closer to the base than for ICaL-75. In the five other cases, all with large heterogeneity, we observed a more complex dynamics. Here, the initial scroll wave generated waves with a high frequency. Such waves were propagating toward the apex where the APD and refractory period were the longest and tissue was not able to recover. As a result, we observed formation of the wave breaks. Figure 4 shows an example of such process. Here, the scroll wave, whose tip is located at the other side of the LV, generates two wavebreaks ( Figure 4B) which evolve to a complex spatiotemporal pattern ( Figure 4C). This occurred in cases IKs-25 and ICaL-25. Interestingly, long filaments with both ends at the epicardial surface sometimes emerged ( Figure 5). Such filament may have a complex behavior which can lead to annihilation of the pattern or formation of new filaments. However, they never stabilized to attractors as in the previous 12 cases. Overall, we can say that, while a moderate increase in heterogeneity resulted in stabilization of filament towards the region of longer APD, a larger degree of heterogeneity resulted in the formation of additional wavebreaks and complex spatio-temporal patterns.

Transmural Heterogeneity
To investigate the effect of transmural heterogeneity, we performed nine simulations: three LVs with different geometry with three degrees of heterogeneity, where we changed IKs. This is because IK is the main ionic current responsible for the transmural heterogeneity. The results of these simulations are shown in Figure 6.  Figure 1C).
In seven out of nine cases, we observed the same type of dynamics as in the absence of the heterogeneity. The scroll wave drifted to base and annihilated there. We saw, however, that the lifetime of the scroll wave decreased with increase in heterogeneity. In two cases (IKs-25-TM; apex thickness 6 and 12 mm), the scroll wave stabilized at the apex and either rotated with a velocity of 0.19 mm/s (h = 6 mm) or did not move at all (h = 12 mm).

Discussion
In this paper, we study the effects of heterogeneity on scroll wave dynamics in a model of the left ventricles of the human heart. We used a detailed ionic model for human ventricular cell and analytical representation of the geometry and anisotropy. We considered two types of heterogeneity: apicobasal and transmural, and we represent them by changing ICaL and IK currents.
The apicobasal heterogeneity was found to change the dynamics of the scroll wave. Instead of annihilation at the base, the scroll wave in most of the cases stabilized at the apex. This result can be explained by the following known dynamics of the spiral wave in cardiac tissue. In [20,21], it was shown that spiral waves drift to the region of longer APD. Because in our case APD was longest at the apex, we can assume that stabilization at the apex is a result of such dynamics. Note that for h = 6 mm and 12 mm the location of the attractor did not change with the degree of the heterogeneity. However, for h = 18 mm, we had a gradual change of the location of the attractor Ψ = 0.72 for ICaL-75, while Ψ = 1.2 for ICaL-50. In addition, for h = 18 mm and IKs-75, the scroll wave still drifted to the base and annihilated there. The mechanisms of the observed phenomenon can be the following. In case h = 18 mm, the apex is substantially thicker than the base. Because in the absence of heterogeneities the scroll wave drifts toward the region of minimal thickness, such gradient attempts to push the scroll wave towards the base. On the other hand, the apicobasal heterogeneity thrusts the scroll wave toward the apex. As a result of interaction of these opposite processes, we observe a stabilization of the scroll closer to the apex for larger heterogeneity in ICaL. In addition, the small heterogeneity IKs-75 was not sufficient to stabilize the scroll and it annihilated at the base.
Another important type of the dynamics is a breakup. Basically in all cases of large heterogeneity, some form of breakup was observed. These results can also be explained in the following way. The formation of wavebreak at the regions with a longer refractory period (APD) is one of the most classical mechanisms of generation of spiral waves in 2D [37]. As in our case, the breakup occurred via formation of wavebreak at the apex, where the APD is longest, we think that the mechanism [37] is also responsible for the onset of breakup in our case.
For the transmural heterogeneity, we did not observe substantial effects on the filament dynamics. It did not change the location of the attractor: in almost all cases as in the control case, the scroll wave drifted towards the base and annihilated there. We saw a slight decrease in the life time of the scroll wave, but the effect was small and not significant. In a few cases of largest heterogeneity, we saw stabilization at the apex. Unfortunately, the mechanisms of such behaviour are not clear. In [22,23], it was shown that transmural heterogeneity affects the period of rotation of a scroll wave, and in extreme cases can cause its breakup. However, studies there were perfomed in a simplified model of cardiac tissue. In [38], a sproing (elongation) of the filament was also observed in a simplified model of cardiac tissue. However, we did not observe such phenomena and found effects that require additional study.
It should be noted that filaments have complex dynamics; however, not all complexity of the dynamics comes from the heterogeneity. As pointed out by Papadimitriou [39], the complexity does not necessarily require heterogeneity and that is definitely true for filaments. For example, even in completely homogeneous cardiac tissue, for certain parameter values, one can have a breakup of filaments due to negative filament tension [40], or dynamic instabilities [41]. In addition, anisotropy by itself can produce complex filament shapes [42]. It would be interesting to further characterize complexity of filament dynamics in the heart in a wider range of parameters and study if heterogeneity will be a determinant for their dynamics as it was in our study.
In a previous paper [26], we studied scroll wave drift on the Aliev-Panfilov model and symmetrical LV model with isotropy and anisotropy 9:1 in homogeneous and heterogeneous myocardium. In that paper, we also found a tendency of scroll wave to move towards the apex and an increase of the rotation speed of waves with an increase in the heterogeneity. The drift speed in the AP model was smaller than one found in this paper, which is a result of using of a more accurate model for cardiac excitation.
Overall, these results show that the scroll wave drifts and stabilizes at the apex of the left ventricle in the presence of moderate apicobasal heterogeneity. Thus, independently of their initial location, the scroll waves are after some transient processes more likely to be observed at the apex of the heart. Therefore, the apex may be a preferable location of the scroll wave. Identification of the location of a source of cardiac arrhythmias is an important question in practical cardiology. One of the most effective ways to treat cardiac arrhythmia is a cardiac ablation procedure, during which flexible catheters are passed into the heart through a blood vessel. Once located near the target region, the catheter is activated to locally destroy or modify the arrhythmia substrate by delivering radio-frequency energy or extreme cold. Thus, locating the arrhythmia source is a problem of great importance. In this paper we show that the source of cardiac arrhythmia is more likely to be located at the apex of the heart. Hence, the modification of the apex by some sort of ablation, which makes rotation of scroll wave impossible, can eliminate the cardiac arrhythmia. Recently, a non-invasive method of ablation was proposed, which enables modifying the tissue not only on the surface but also inside the myocardial wall by noninvasive delivery of radiation with stereotactic body radiation [43]. For such methodology, the location of the arrhythmia source is again crucial. Note, however, that ablation in 3D ventricle to eliminate a scroll wave is a non-trivial issue and needs to be further studied, e.g., using a modeling methodology applied in our paper.
A limitation of our approach is that we considered only spiral waves with one chirality. The chirality may substantially affect the dynamics of scroll waves [29]. In this paper, we chose that particular case for counterclockwise rotation as without heterogeneity we had always one regime (annihilation at the base). It was thus straightforward to study the effect of heterogeneity in that situation. The effect of chirality in the presence of heterogeneity should be investigated in a subsequent study.
We studied transmural and apicobasal heterogeneity separately. It would be also interesting to study their combined effect on scroll wave dynamics.
In this paper, we performed our studies only in a limited parameter range. This is because such modeling is extremely challenging, as the studied effects can be observed only in anatomical models of the heart on a long time scale (of the order of minutes) while typical simulations in anatomical models of the heart are normally performed for a few seconds (see e.g., [44]). It would be interesting in the future to extend studies for a larger parameter range and also model direct clinical interventions which can stop arrhythmias caused by scroll waves, such as catheter ablation.

Conclusions
We performed a study of filament dynamics in the left ventricle of the heart using a state-of-art ionic model for cardiac tissue and combining all the most important factors which can be responsible for the filament drift: the shape of the ventricle, variation in thickness, anisotropy, and heterogeneity. Under those conditions, we find that, in spite of the importance of all factors, the apicobasal heterogeneity was the most important in determining the final position of the filament. For moderate heterogeneity, the filament was stabilized at the apex of the ventricle, while for a large heterogeneity a break-up into a complex spatio-temporal pattern was observed. The transmural heterogeneity did not substantially affect the filament dynamics. As predicting of location of the arrhythmia sources is important for clinical interventions, our results can help in developing clinical procedures to remove arrhythmias organized by filaments.