Electrical wave propagation in an anisotropic model of the left ventricle based on analytical description of cardiac architecture.

We develop a numerical approach based on our recent analytical model of fiber structure in the left ventricle of the human heart. A special curvilinear coordinate system is proposed to analytically include realistic ventricular shape and myofiber directions. With this anatomical model, electrophysiological simulations can be performed on a rectangular coordinate grid. We apply our method to study the effect of fiber rotation and electrical anisotropy of cardiac tissue (i.e., the ratio of the conductivity coefficients along and across the myocardial fibers) on wave propagation using the ten Tusscher-Panfilov (2006) ionic model for human ventricular cells. We show that fiber rotation increases the speed of cardiac activation and attenuates the effects of anisotropy. Our results show that the fiber rotation in the heart is an important factor underlying cardiac excitation. We also study scroll wave dynamics in our model and show the drift of a scroll wave filament whose velocity depends non-monotonically on the fiber rotation angle; the period of scroll wave rotation decreases with an increase of the fiber rotation angle; an increase in anisotropy may cause the breakup of a scroll wave, similar to the mother rotor mechanism of ventricular fibrillation.

Introduction

The modeling of cardiac electrical function is a well-established area of research that began with early models of cardiac cells developed by D. Noble [1].

The importance of modeling in cardiology comes from the widespread prevalence of cardiac disease. For example, sudden cardiac death is the leading cause of death in the industrialized world, accounting for more than 300,000 victims annually in the US alone [2]. In most cases, sudden cardiac death is a result of cardiac arrhythmias that occur in the ventricles of the human heart [2].

When studying cardiac arrhythmias, it is important to understand that they often occur at the level of the whole organ and in these situations cannot be reproduced in single cells. Therefore, it is very important to model cardiac arrhythmias at the tissue level, preferably using an anatomically accurate representation of the heart. Compared to modeling at the single-cell level, anatomical modeling started much more recently [3], [4]. Using anatomical models, researchers have been able to obtain important results on the 3-D organization of cardiac arrhythmias in animal [5] and human [6] hearts. Moreover, the defibrillation process has been investigated [7], and the effects of mechano-electrical coupling on cardiac propagation have recently been modeled[8], [9]. Multi-scale anatomical cardiac modeling is becoming increasingly prominent in medical and pharmaceutical research [10].

In a broad sense, an anatomical model of the heart is a combination of models of cardiac cells and anatomical data. The development of models of the electrical and mechanical functions of cardiac cells is a well-established area of research, and many models have been developed, including models of the human cardiac cells [11]–[20]. The anatomical data necessary for cardiac modeling include not only the heart's geometry but also its anisotropic properties. Although the geometry of the heart can be obtained from routine clinical procedures such as MRI or CT scans [21], [22], anisotropy data are much more challenging to acquire. Currently, the acquisition can be done on explanted hearts only, using either direct histological measurements or time-demanding DT-MRI scans [23]–[26]. In addition to experimental noise, even perfect measurements will grant only the particular anisotropy of the imaged heart. Thus, to study the effects of anisotropy on wave propagation, one needs to vary the anisotropic properties and to separate the anisotropy effects from other factors. All these questions can be addressed with the development of models that account for the anisotropy of the heart using analytical or numerical tools.

In a previous article [27], we described an axisymmetric model of the left ventricle (LV) of the human heart. In the model, we represented the LV shape (including positions of cardiac fibers) as analytical functions of special curvilinear coordinates defined on a rectangular domain. Our model allowed the generation of not only a default architecture of anisotropy closest to the reality but also intermediate architectures that can be used to study the effects of any specific element of anisotropy on wave propagation in the heart.

In this paper, we build on our previous approach in two ways. First, we develop a numeral scheme for the integration of equations for wave propagation in our anatomical model of the LV, which is the best possible way to account for anisotropy. In particular, we develop a model on a rectangular domain and represent anisotropy and the LV shape by means of parameter changes. Second, we vary the geometry and anisotropy parameters to study how the rotation of the fiber orientation affects wave propagation and show that rotational anisotropy accelerates the spread of electrical excitation in the heart. We also study the behavior of scroll waves and their filaments. We show that the scroll waves drift and we calculate their drift velocity and period of rotation depending on the fiber rotation angle and the diffusion coefficients ratio.

Model Description

Geometrical model of the LV

In our model, the LV is represented as a body of revolution around the vertical axis Oz with the shape fitted to experimental data (for details, see [27]). A section of the LV is shown in Fig. 1. The rotation of the blue line delineates the epicardial surface, while the rotation of the red line yields the endocardial surface of the LV.

Figure 1

A radial section of the endocardial (solid red line) and epicardial (dashed blue line) surfaces of the LV model, from [27].

In our model, each point of the LV has three local coordinates (γ, ψ, ϕ). The coordinate γ (γ 0≤γ≤γ 1) gives us points between the endo- and epicardial surfaces in Fig. 1, i.e., it is a measure for transmural depth; the coordinate ψ (0≤ψ≤π/2) is explained in Fig. 1; and the coordinate ϕ (0≤ϕ<2π) is the rotation angle around the vertical axis Oz. The local coordinates are linked with the cylindrical coordinate system (CS) (ρ, φ, z) by the following formulae: with:

Below, we also useThe physical meaning of the parameters is as follows: r b is the LV cavity radius on the LV equator, z b is the LV cavity depth, l is the LV wall thickness on the LV equator, h is the LV wall thickness on the LV apex, and is a dimensionless parameter influencing the conicity-ellipticity characteristic of the LV shape. In this system, γ = γ 0 gives the epicardium and γ = γ 1 gives the endocardium; the value ψ = 0 describes the upper (basal) boundary of the LV.

Following Pettigrew's idea [28] about spiral surfaces and semicircle chords mapping on the surfaces, our LV model consists of spiral surfaces on which a set of curves is defined. The details are described in our previous work [27] and are summarized in Appendices A (spiral surfaces construction) and B (fiber equations). In Figs. 2 and 3, we show common views of spiral surfaces and the fibers inside them.

Figure 2

A spiral surface.

The lines on the surface have equations ρ = const and ϕ = const. Color corresponds to height (z coordinate).

Figure 3

A spiral surface viewed from the top (left panel) and side (right panel).

Two myofibers are displayed in red and black.

We demonstrated in [27] that the model approximates the real fiber orientation field in the LV reasonably well. A comparison of true fiber angle α and helix angle α 1 with MRI data showed that fiber architecture in the equatorial region of the heart was well reproduced in our anatomical model. In the middle (by height) and apical areas, the angles were reproduced both qualitatively and quantitatively well; the difference between the model and the experimental data was not more than 25°.

Overall, we can consider the anatomical model as a map from a rectangular domain γ 0≤γ ≤γ 1; 0≤ψ≤π/2; 0≤ϕ<2π to the shape of LV with anisotropy explicitly given by Eqs. (B.1)–(B.3). The total fiber rotation angle Δα 1 is defined as the difference between the epicardial and endocardial helix angles α 1 measured at the LV basal zone ψ = π/8 (see [29] for details). It can be varied by changing the values of the parameters γ 0 and γ 1:where p is the tangent vector to the parallel γ = const, ψ = const at the point γ = γ, ψ = ψ *, ϕ = 0 (the geometric model is axisymmetric; therefore the choice of ϕ is arbitrary); n is the normal vector to the epicardium passing through the same point; v is the fiber direction vector (it is defined by Eq. (B.1)–(B.3)); pr is the projection operator; and denotes the angle between the vectors u and w. We have chosen the value ψ = π/8 because in this case, the total fiber rotation angle changes uniformly enough depending on the values of γ 0,1 we use (see Table 1 and section “Parameter values” below). For short, below we will denote Δα 1(γ 0, γ 1) as α without argument. We use it in the present study to investigate the effect of fiber rotation on wave propagation in the heart.

Table 1

Dependence of the total fiber angle on the model parameters .

Model	γ 0	γ 1	The helix angle near the base at		The total fiber
			the epicardium	the endocardium	rotation angle α
1	0.3	0.55	−13°	3°	16°
2	0.2	0.7	−40°	29°	69°
3	0.1	0.85	−69°	64°	133°
4	0	1	−87°	87°	174°

Electrophysiological model

To describe the excitation of cardiac tissue, we use the detailed ionic model for human ventricular cells from [6], [11]. The model uses reaction-diffusion equations to describe the evolution of the transmembrane potential u = u(r, t):

Here, the intracellular processes are captured by I = I(r, t) which is the sum of the ionic transmembrane currents; C is the capacitance of the cell membrane. The locally varying diffusion matrix D accounts for myofiber anisotropy. As in [6], the diffusion matrix D = (D) was computed using the following formulawhere D 1 and D 2 are the diffusion coefficients along and across the fibers, v is the unit vector of fiber direction, i, j are Cartesian indices, and δ , is the Kronecker symbol.

Numerical integration scheme and boundary conditions

The aim in this paper is to present a numerical procedure that allows us to use the analytical representation of cardiac anatomy and anisotropy described in the previous section. In particular, as our anatomical model is just a map from a rectangular domain γ 0≤γ≤γ 1; 0≤ψ≤π/2; 0≤ϕ<2π to the shape of LV, we can formulate our approach in that rectangular domain. The shape of the heart, as well as anisotropy, will then be a curvilinear coordinate system (1)–(2) defined on that domain. We need to recalculate Eq. (4) with no-flux boundary conditions in those coordinates. A long computation presented in appendix C results in the following expression of the diffusion term:where ξ 1 = γ, ξ 2 = ψ, ξ 3 = ϕ, and p and q are coefficients given by the explicit analytical Eqs. (C.12) and (C.13) that depend only on the geometry of the LV and on the diffusion matrix. This representation is similar to that presented by Sridhar et al. in [30]. However, in our work, it was done for the 3-D case and for a special form of the diffusion matrix (see Eq. (6)).

For numerical integration of the model (B.1)–(B.3), (4), (5), we use the explicit finite difference method on a discrete grid in the (γ, ψ, ϕ) space. We initially use a uniform grid with the γ-indices of the nodes denoted as i = 0, 1, … n; the ψ-indices as j = 0, 1, … n; and the ϕ-indices as k = 0, 1, … n.

Although this grid is uniform in the (γ, ψ, ϕ) space, it is non-uniform in real space because distances between the grid points substantially decrease when ψ approaches π/2, which is similar to the situation at the pole in a polar CS. Note, however, that even the cubic Cartesian lattice in real space is not uniform with respect to anisotropic diffusion due to the curved space interpretation of anisotropy [31], [32]. To account for this problem, we exclude some points from our uniform grid in the following way. We first choose a threshold value of distance d min. Then, at γ = γ 1 (i.e., at the epicardial surface) and any given ψ = ψ, we calculate the distances between the node at ϕ = 0: (γ = γ 1, ψ = ψ, ϕ = 0) and node (γ = γ 1, ψ = ψ, ϕ = ϕ). We find minimal k satisfying two conditions: (1) the distance to the k-node from the node at ϕ = 0 is more than the threshold value d min; and (2) k is a divisor of n. We denote this number as K (as it depends on ψ). If ψ is far from π/2, then K = 1 and we use all nodes of our uniform (γ, ψ , ϕ) grid. When ψ approaches the value of π/2, K>1 and we drop all nodes between (γ = γ 1, ψ = ψ, ϕ = 0) and (γ = γ 1, ψ = ψ, ϕ = K), then the next node will be (γ = γ 1, ψ = ψ, ϕ = 2K) etc. Thus, only nodes with the ϕ-indices 0, K, 2K, … will be taken for evaluation of the Laplacian. After each time step, we compute values for all variables of our model in the omitted nodes using linear interpolation. With this approach, we reduce the number of grid elements in the ϕ-direction in the apical region; otherwise, we would have needed to lower the maximal time step in our explicit integration scheme, which would have slowed our computations significantly.

The no-flux boundary condition in our problem iswhere n is a normal to surface. We rewrite this equation in our special CS and obtain the following expression:where c γ, c, c are coefficients given by Eqs. (D.14) and (D.16) in Appendix D. In order to satisfy the boundary condition, we add nodes at the domain boundaries in the following way. In the special CS, the domain of integration is a rectangle (with periodic boundary in ϕ), and at ψ = π/2 we have a pole (i.e., also no boundary). Therefore, we have only three boundaries, namely at γ = γ 0 (i.e., i = 0), γ = γ 1 (i = n), and ψ = 0 (j = 0). Fictitious layers with (−1, j, k), (n +1, j, k), and (i, −1, k) are added. We then solve equation (9) on the three boundary surfaces to find values in the added nodes. Subsequently, we can compute Laplacian at all other nodes in the domain using, if necessary, the values in the additional nodes. Due to this procedure, all the nodes lying on the boundaries will satisfy the boundary conditions.

We have programmed this approach using C in the CodeBlocks IDE, a Mingw compiler. The calculations were performed in operating systems Windows 7 and Linux. The OpenMP and MPI technologies have been used for parallelization, and Paraview and Irfan have been used for visualisation. The formal parameters of our numerical scheme have been given in the methods section above. Such an approach allows us to compute various regimes of wave propagation in a model of LV with good representation of boundary conditions and to study various effects of anisotropy on wave propagation patterns.

We also studied the dynamics of scroll waves using the ten Tusscher–Noble–Noble–Panfilov (TNNP) model [11] and various anisotropy parameters. We initiated a scroll wave using the S1–S2 stimulation protocol and studied its dynamics 12 s. For drift velocity and rotation period calculations, we take into account only the last 8 seconds of the simulations to exclude transient processes. We calculated average periods in the section ϕ = 0, that is, x = ρ>0, y = 0. Filaments were analyzed as reported in [6], [33]. Finally, we computed their average velocity v in the Cartesian CS.

Parameter values

We used the following parameters from Table 2 in [29]: the LV equatorial radius mm, the equatorial wall thickness l = 10 mm, the LV cavity depth mm, the apical wall thickness h = 7 mm, the conicity-ellipticity parameter , and the spiral surface torsion angle (see Fig. 3c in [29]). The threshold distance between the adjacent nodes d min was set to 0.3 mm. Currently, the first four parameters are measured using modern experimental techniques such as MRI (see [34]–[36]). The values used in our paper are in agreement with these experimental data.

Table 2

The initial excitation areas.

Series	γ-indices	ψ-indices	ϕ-indices	area
A	i≥nγ−4	j≥nψ−4	all	apical epicardium
B	i≥nγ−4	|j−(nψ/2)|≤2	k≤4	central epicardium
C	i≤4	|j−(nψ/2)|≤2	k≤4	central endocardium

Our mesh had a distance of 0.2 to 0.3 mm between the nodes and before the deletion of the nodes described above had n = 40, n = 300, n = 800. The diffusion coefficient along the fibers was D 1 = 0.3 mm2/ms. The diffusion coefficient across the fibers D 2 was varied between different experiments depending on whether we modeled isotropy or anisotropy.

We applied our approach to study the effect of anisotropy on the spread of excitation in the heart. In particular, we initiated a wave at several locations and studied how the wave arrival time depends on the two main features of anisotropy. Our first anisotropy parameter was the ratio of the diffusion coefficients D 1:D 2 along and across the fibers. Also, we independently varied the total rotation angle of fibers through the myocardial wall by adjusting γ 0, γ 1 and keeping wall thickness constant. We also compared our results with the spread of excitation in an isotropic model of the LV where D 1 = D 2.

The dependency of wave velocity on the direction of propagation in the heart was measured in [37]–[39]. Experimental data show that the ratio of longitudinal to transverse conduction velocities ranges between 3 [38] and 2.1 [39]. Since D ∝ c 2, where c is wave velocity (see, e.g., [40]), we used the ratios D 1∶D 2 = 1∶0.111 and D 1∶D 2 = 1∶0.25, which correspond to the experimental data. These anisotropy ratios were also used in the modeling studies [33], [41]–[44].

For point stimulation, we increased the value of the variable u from the resting potential of −86.2 mV to u = 0 mV at the first time step in small regions located at three different sites. In the A series, it was a small epicardial region at the apex; in the B series, at the centre of LV epicardium; in the C series, at the centre LV endocardium (see Table 2).

In this paper, we study the effect of the fiber rotation on the spread of excitation. With this purpose, we generated a series of LV models that differ in the total fiber rotation angle through the myocardial wall. The parameters of the model are listed in Table 1. Note that although the values of γ 0 and γ 1 differ between the models, they affect only fiber rotation, and the LV geometry is exactly the same for all models due to the rescaling procedure described in Sec. 1.4 in our previous article [27]. Also note that the change in fiber rotation results in the change in fiber angle at the epicardial surface (see column “α” of Table 1).

The time step was equal to 0.005 ms for the isotropic cases and 0.01 ms for the anisotropic ones.

We considered that a wave came to a node when the potential in the node was more than −80 mV the first time.

Numerical Results

Activation maps

Figure 4 shows the wave arrival time after stimulation of the small apical epicardial zone A for ratios D 1∶D 2 equal to 1∶0.111 or 1∶0.25 (shown at the top of the figure) and four different rotation angles of the fibers, which are displayed in the left column. We see that in this first example, all the figures are axisymmetrical, which is a consequence of the axisymmetric properties of our model and the initial conditions.

Figure 4

Arrival times, in ms, of the waves after point stimulation at the apex for various values of anisotropy and fiber rotation.

The values of anisotropy are shown at the top of the figure and the values of the fiber rotation are shown in the left column. For details, see Table 1. Arrival times are color-coded in ms.

We observe that for the low rotation angle (the upper row), the speed of the upward wave propagation for the diffusion coefficient ratio of 1∶0.111 is substantially smaller than that for the ratio of 1∶0.25. However, if the fiber rotation angle increases (the lower rows), the difference in the speed between the two anisotropy ratios decreases. For the rotation angle of 174° (the lower row), the excitation patterns for both anisotropy ratios become close to each other. Thus, we observe that fiber rotation increases the velocity of the spread of excitation and also decreases the effect of anisotropy.

Let us now consider the case of lateral stimulation for a given ratio of the diffusion coefficients of 1∶0.111; the results for epi- and endocardial stimulation are presented in Figs. 5 and 6. After epicardial stimulation, the wave initially follows the fiber direction. In Fig. 5, we note a displacement of the early activation zones (red) due to the change in fiber direction at the epicardium in our anatomical model (see column 4 of Table 1). However, for endocardial stimulation (Fig. 6, the first column), the shift of the early activation zone (red) on the epicardium is attenuated by fiber rotation. As in Fig. 5, we see that an increase in the rotation angle causes a decrease in the arrival time. In addition, in the second column in Fig. 5, the excitation patterns have a clear V shape on the surface opposite the stimulation site, which is a mere consequence of the shape of the heart (see [3], [4]).

Figure 5

Arrival times, in ms, of the waves after point stimulation at the epicardial surface for a large anisotropy ratio D 1∶D 2 = 1∶0.111.

The notation is the same as in Fig. 4.

Figure 6

Arrival times, in ms, of the waves after point stimulation at the endocardial surface for a large anisotropy ratio D 1∶D 2 = 1∶0.111.

The notation is the same as in Fig. 4.

Figs. 7 and 8 show the arrival time after lateral stimulation in a case of decreased anisotropy, i.e., for a diffusion coefficient ratio of 1∶0.25. The excitation patterns resemble those from Fig. 5 and Fig. 6 respectively, with similar effects of fiber rotation on the epicardial stimulation and V-shaped patterns. Here, we also observe that an increase in the rotation angle decreases the overall excitation time. A compensating effect of fiber rotation on the degree of anisotropy can be noted. The difference between the corresponding panels in Fig. 5 and Fig. 7 (and also Fig. 6 and Fig. 8) is more pronounced for a lower fiber rotation rate (rotation angle 16°).

Figure 7

Arrival times, in ms, of the waves after point stimulation at the epicardial surface for an intermediate anisotropy ratio D 1∶D 2 = 1∶0.25.

The notation is the same as in Fig. 4.

Figure 8

Arrival times, in ms, of the waves after point stimulation at the endocardial surface for an intermediate anisotropy ratio D 1∶D 2 = 1∶0.25.

The notation is the same as in Fig. 4.

Average speed of excitation

Now let us quantify the effects of rotation and anisotropy on wave propagation. In order to do this, we use the following procedure. We group all points of the heart to bins differing by their “distance” from the stimulation point. We define the distance as the arrival time from the stimulation to a given point. To calculate the distances, we perform simulations in which we initiate a wave at the same locations as in Figs. 4–8. However, for the isotropic medium, we use a diffusion coefficient of D 1 = D 2 = 0.3 mm2/ms. We generate 40 groups in which points differ in the arrival time by 2 ms. Then we determine the average arrival time for each of these groups for various anisotropic conditions and compare these arrival times to the arrival time in the isotropic model. As in the isotropic model, the velocity of wave propagation in all directions is the same; this dependence gives us the dependence of the wave arrival time on the distance from the stimulation point.

In Fig. 9, the red lines correspond to α = 174° and the black lines correspond to α = 16° fiber rotation angle in the LV wall. The fact that the red lines are always located below the black lines shows that the increase of the fiber rotation angle results in faster wave propagation.

Figure 9

Arrival times, in ms, as a function of the distance from the stimulation point for the apical (A), epicardial (B), and endocardial (C) stimulation.

The distance on the horizontal axis is measured in ms as the arrival time of the wave in the isotropic model (see text for more details). The red lines represent numerical experiments for total rotation angle α = 174°; the black lines represent for α = 16°; and the blue lines represent isotropy. The solid lines correspond to the case D 1∶D 2 = 1∶0.111, while the dashed lines correspond to the case D 1∶D 2 = 1∶0.25. The vertical segments display minimal and maximal arrival times in each group of nodes. The average, min, and max arrival times are displayed on the leftmost panels for D 1∶D 2 = 1∶0.111 and in the middle column for D 1∶D 2 = 1∶0.25. The right column compares the average arrival times.

Also, all solid lines correspond to D 1∶D 2 = 1∶0.111, while dashed lines correspond to D 1∶D 2 = 1∶0.25. If we now compare the solid and dashed lines of the same color (third column in Fig. 9), we see that the red lines are closer to each other than the black lines. This indicates that in presence of higher fiber rotation (the red lines) the decrease of D 2 (i.e., solid vs. dashed) has a smaller effect on the arrival time. This once again illustrates that anisotropy is compensated for by the rotation of the fibers.

In addition, we see in Fig. 9 that the red lines always have a less steep slope than the corresponding black lines. As going from black to red shows an increase in fiber rotation, we can conclude that the increase in rotation makes propagation faster in all cases.

Scroll wave dynamics

We have also studied scroll wave dynamics for the same values of anisotropy and fiber rotation. We generated a single scroll wave located approximately at the middle between the apex and the base of the ventricle and studied its behavior for different model parameters γ 0, γ 1 and D 1∶D 2. We found that anisotropy substantially affects the dynamics of scroll waves. In all cases, the increase in the fiber rotation angle results in a decrease in the period of scroll wave rotation (Fig. 10). We see that for D 1∶D 2 = 1∶0.25, when the fiber rotation angle is increased from 16° to 174°, the period drops significantly, from 277 to 257 ms. For D 1∶D 2 = 1∶0.111, we see similar dependency. In addition, the period for the same rotation angle for D 1∶D 2 = 1∶0.111 was slightly longer than for D 1∶D 2 = 1∶0.25.

Figure 10

Scroll wave rotation period T, ms, as a function of total fiber rotation angle α, deg.

The scroll wave dynamics was also substantially affected by the anisotropy. In all cases, we observed a drift of the filament (Fig. 11). The drift always had two components, both in the vertical (ψ) and circumferential (ϕ) directions. The total velocity of drift (Fig. 12A) was very small, about 1 mm/s, which is about 0.2 mm per rotation; however, the drift was monotonic and persistent. The value of velocity had no clear relationship with the rotation angle; for D 1∶D 2 = 1∶0.25, we see some tendency for velocity decrease with an increase in fiber rotation, while for D 1∶D 2 = 1∶0.111, the dependency is strongly non-monotonic and it is maximal for the intermediate values of the fiber rotation angle. The drift direction can be seen from the sign of the vertical and horizontal components of the velocity (Figs. 12B, C). Here again, the direction is affected by the rotation angle; however, we also did not find any clear tendency for either drift to the apex or base of the heart depending on the rotation angle.

Figure 11

Potential, mV, on the LV surface during scroll wave rotation (left) and tip trajectory for D 1∶D 2 = 1∶0.111 (red line) and for D 1∶D 2 = 1∶0.25 (black line) (right).

The results are shown for model 2 (γ 0 = 0.2, γ 1 = 0.7, see text and Table 1 for details).

Figure 12

Velocity of scroll wave filament drift for the simulation of 8 s.

Average filament velocity V, mm/s (A). Velocity components multiplied by 1000, per second: latitudinal component v (B) and longitudinal component v (C).

For D 1∶D 2 = 1∶0.25, the initial scroll wave was always stable and did not break down to multiple scroll waves. For D 1∶D 2 = 1∶0.111, we did observe formation of the additional sources of excitation. However, in most cases, they appeared simultaneously at a substantial distance from the initial filament and not as a result of filament buckling and breakup due to rotational anisotropy in the way it was reported in [33]. The onset of new sources had a clear correlation with the fiber rotation angle (Fig. 13). We did not observe any instabilities for small and big α (models 1 and 4 in Table 1, Figs. 13A, D), however, for intermediate and large values of α (Figs. 13B, C), we observed new sources, and their number increased with the increase of the fiber rotation angle (compare panels B and C). Note that cases presented in panels B and C correspond to the largest drift velocities of a scroll wave; thus, the onset of secondary filaments may be related to the filament drift.

Figure 13

Scroll wave filaments in the LV model.

The anisotropy ratio is D 1∶D 2 = 1∶0.111. Panels A, B, C, D: models 1, 2, 3, 4 (see Table 1), fiber rotation angle in the LV wall increases from panel A to panel D. The epicarduim (semitransparent colored surface; color denotes height from the red base to the purple apex), the endocardium (opaque white meshy surface), and filaments (black lines and dots).

We have also studied change in filament shape over time. For small values of total fiber angle α, the filament remained transmural, nearly straight, and stable. This case is shown in Fig. 13A. For intermediate α equal to 69° and 133°, the filament not only drifted faster but also deformed to a transmural S or L-shape (Fig. 13B). For larger values of α, the filament again had a nearly straight shape (Fig. 13D).

Discussion

In this paper, we used our recent anatomical model of the LV of the human heart using a special CS (γ, ψ, ϕ) which gives an explicit analytical map from a rectangular domain to the heart shape and fiber orientation field. This allowed us to represent the heart's geometry on a rectangular grid and explicitly write expressions for boundary conditions. This approach may be helpful for studies of any phenomena in which boundary effects are of great importance.

One important feature of our model is the possibility to change the properties of anisotropy. The most significant characteristic of LV anisotropy is the rotation of myocardial fibers through the myocardial wall. As shown in [45], relatively simple rule-based global models of LV myofiber directions yield no worse results than complicated image-based locally optimized models. Since any locally optimized model needs to be regular enough and the regularity requires smoothing and interpolation, no image-based model can be an untouched copy of a real heart.

We can change the degree of the fiber rotation in a consistent way and study its effect on normal and abnormal wave propagation in the heart. In this paper, we investigated two main features of wave propagation using the detailed human ventricular electrophysiological model: the distribution of effective excitation speed and the dynamics of transmural scroll waves.

In our study of the effect of the fiber rotation and anisotropy on wave propagation, the initial stimulation area was located at the apex and on the lateral epi- and endocardium of the LV. We found that the rotation of myocardial fibers accelerates the spread of excitation waves in the heart, which was explicitly demonstrated using models with different fiber rotation angles. This acceleration of wave propagation was discussed in [32]; it occurs because the wave can propagate with maximal speed in more directions with a larger rotation angle, which results in an overall faster wave propagation. Note that if the rotation angle is 2π or more, propagation with maximal speed will be possible in any direction, and in the limit of a large medium, the arrival time will be the same as in an isotropic medium with the velocity determined by the velocity along the fiber [32]. We were able to demonstrate this in an anatomical setup in which we explicitly changed the rotation of the fibers, while in [32] such an assessment was made for a single anisotropy configuration. While Young and Panfilov represented the tissue as a simple rectangular 3-D slab with plane–parallel fibers [32], we adopt a fully 3-D fiber architecture together with an LV shape. In this paper, we particularly considered a more realistic ventricular architecture and morphology that includes the following features:

A more realistic method for the fiber rotation angle around the axes [27], so called “Japanese-fan arrangement” [29]; and

A realistic change of the fiber rotation angle values due to the displacement of the transmural axis from the LV apex to the base [27].

As in [32], our model also shows the faster excitation propagation with an increase in the fiber rotation angle and a decrease in the anisotropy.

A second set of results in this paper concerns the dynamics of transmural scroll waves. The negative correlation found for the rotation period versus fiber rotation angle in Fig. 10 is different from the observations in [46] made by Qu et al., who observed an increasing period with a faster fiber rotation rate. Note, however, that their simulations used in-plane fiber rotation, while we work with a 3-D ventricular geometry. Both complementary cases can be qualitatively explained on geometrical grounds. In [46], the period was only affected by fiber rotation, which causes negative intrinsic curvature of the associated curved space [32]. By its definition, negative geometrical curvature represents a saddle-like space with positive angular deficit. More specifically, in a space with curvature , the circumference C of a circle (ball) with small radius r amounts to [47] As the spiral tip in each cross-section needs to travel along a larger closed path of length C before completing a period, negative is expected to increase the rotation period. Therefore, the trend found in [46] can be expected if the fiber rotation is the main determinant of the rotation period. In our present model, however, a transmural filament consists of spiral waves' tips in different layers of constant depth γ. These γ-surfaces are sphere-like and therefore have positive . Thus, under normal excitability, the sphericity of the LV will decrease the rotation period [48], [49]. We conclude that in general, the rotation period of scroll waves in the heart may decrease or increase with increasing fiber rotation angle α, depending on the relative strengths of the competing effects of the fiber rotation rate and the extrinsic curvature (sphericity) of the LV cavity.

Regarding the non-monotonic dependence of filament drift velocity versus total rotation angle, we first note that the creation of secondary filaments in the regime of intermediate α is consistent with the clinical or experimental picture of a “mother rotor” during cardiac arrhythmias [50], [51]. In such a scenario, the primary filament (mother rotor) remains stable and creates secondary sources that further disturb the electrical excitation of the heart, leading to cardiac arrest. In our simulations, we also see a similar situation with a stable mother rotor, which was always sustained until the end of the simulation time, and secondary excitation sources induced at some distance from it.

In addition to the direct simulation results detailed and discussed above, our anatomical model [27] coupled to detailed electrophysiological model [11] may prove useful in future studies for the following reasons. First, our model can be used to investigate the possible contribution of the LV geometry to the propagation of the excitation waves. In particular, in certain heart diseases (e.g. dilated and hypertrophic cardiomyopathy, eccentric and concentric cardiac hypertrophy, etc., see chapter 8 in [52]), the shape of the ventricle becomes more spherical and the thickness of the wall also increases. Such changes in geometry can easily be accommodated in our model. In general, the study of the effects of the LV geometry on excitation seems to be of great importance because many cardiac pathologies tightly correlate with changes in the LV geometrical characteristics. The LV becomes more dilated near the apex and thicker near the base during stress-induced (“Takotsubo”) cardiomyopathy, or transient apical ballooning syndrome (see chapter 8 in [52]). Such remodeling of the LV geometry might be mimicked in similar mathematical models via the fitting of the geometric parameters values to account for the LV shape of a particular pathology.

The results of our simulations can be verified by direct measurements of wave propagation on the whole heart preparations, such as [53]. Note, however, that another factor important for overall excitation of the heart is the Purkinje conduction system. In order to compare experiments with our simulations, such measurements should be performed after chemical ablation of the Purkinje system [54].

Another potentially interesting application of this approach may be its application in studies on the defibrillation of cardiac tissue, in which case tissue texture and boundary conditions are also of key importance [55]. However, a bi-domain representation of cardiac tissue must be used for defibrillation [55]–[57], which does not fall under our present scope. Nonetheless, the extension of our approach for such cases is straightforward. The formulae for the diffusion term (7) and for the boundary conditions (9) can be directly used for representation of the bi-domain equations in the special CS. Then, the finite difference problem can be formulated in the same way as in our case and can be solved using any existing method (see [58]).

In this paper, we studied wave propagation due to point stimulation and scroll waves. Other important wave propagation regimes include various types of scroll waves and turbulent patterns [19], [59]–[61]. It was shown that heterogeneity [62], [63] and anisotropy [31], [64], [65] of the tissue are significant factors determining the dynamics of these sources. The effect of heterogeneity on the dynamics of spiral waves was also studied in a series of papers by Shajahan, Sinha and co-authors [30], [41], [55], [67]. In particular, in [66] they showed that some changes in the position, size, and shape of a conduction inhomogeneity can transform a single rotating spiral to spiral breakup or vice versa. Since our model provides tools for changing anisotropic properties and allows one to add heterogeneity, these effects can also be studied using our approach.

Appendices

A Construction of spiral surfaces

We model myocardial sheets as surfaces filling the LV. The filling was obtained by rotation of one surface around the vertical LV symmetry axis. We call these surfaces “spiral” (see Fig. 2).

We introduce the special CS (γ, ψ, ϕ), which is linked with the cylindrical CS (see (1) and (2)). In this CS, the equation of the spiral surface iswhere different values of ϕ 0 allow us to obtain different spiral surfaces and ϕ max is a constant of the model.

The parametrical equation of a spiral surface in cylindrical CS is

B The myofibers' equations

The equations are (see Fig. 3):where different values of the parameter Y ϵ (0, 1) correspond to different myofibers, is the explicit equation of a spiral surface in the cylindrical CS.

At each point (γ, ψ, ϕ), 0≤ψ≤π/2, a fiber orientation vector v = (x′, y′, z′) is given by the formulae: where ρ is determined by Eq. (1),

where ϕ max is a dimensionless parameter affecting fiber twist.

C Laplacian in implicit curvilinear coordinates

The Laplacian is an important term in the reaction-diffusion equation as it is responsible for the modeling of electrical wave spreading. It can be written as div(D grad f) where f is the transmembrane potential and D is an anisotropic local diffusion matrix.

Below, we calculate the Laplacian for anisotropic diffusion in the Cartesian and the special CS.

C.1 The Cartesian CS

For an arbitrary diffusion matrix D = D In consideration of Eq. (6), one can write:

C.2 The curvilinear CS

Here we deduce from Eq. (C.2) the proper form of the Laplacian in the curvilinear coordinates ξ 0, ξ 1, ξ 2. Let us consider v and f as functions of ξ 0, ξ 1, ξ 2 We calculate three types of derivatives. First, one haswhereSecondly, we evaluateand thirdly,where

The difficulty now lies in the fact that the functions x(ξ) define the ξ only implicitly. To evaluate the necessary derivatives, we need the following matrices:

The matrices are linked between themselves with the following relations: We substitute (C.3), (C.5), (C.7) to (C.2) and get:where

q are elements of matrix Q:

D Boundary conditions

Let n be a normal vector to one of the LV boundary surfaces. For outer domain boundaries, we use a no-flux condition on the current.

D.1 Isotropic case, cylindrical CS

One can write the boundary condition as

By the definition of directional derivative and since the normal vector to the LV boundary lies in our problem always in the corresponding meridional half-plane,where n and n are the normal vector components in the meridional half-plane.

D.1.1 The equator

On the equator, n = 0, so (D.2) reduces to

D.1.2 The epicardium

On the epicardium, n = (z+h) cos ψ, , so one can write (D.2) as

Everywhere on the epicardium, except the apex,

On the apex, cos ψ = 0, , so the boundary condition looks like .

D.1.3 The endocardium

On the endocardium, n = z cos ψ, , so (D.2) can be written as

Everywhere on the endocardium, except the apex,

On the apex of the endocardium, like the epicardium, the boundary condition is .

D.2 Isotropic case, special CS

In the special CS, we use and . Let us take into account thatand

D.2.1 The equator

Formula (D.3) can be rewritten as

D.2.2 The epicardium

Formula (D.4) becomes

D.2.3 The endocardium

Formula (D.6) yields

Note that if , the two last formulae cannot have zero values in the denominators and can be directly applied to the LV apex points.

D.3 Anisotropic case, special CS

The boundary condition iswith n, the normal to the LV surface.

D.3.1 The epi- and endocardium

Let us write (D.13) in detail: Here, the J , are columns of derivatives of these special variables by Cartesian coordinates (see (C.8)):and so on.

D.3.2 The equator

Let us write (D.13) detailed as the following:

Vector n is collinear to the Oz axis, so n = n = 0. The equation (D.15) goes over

Writing down the matrix product:we substitute (C.5) to this equation:Let us express :We can calculate the derivatives of γ, ϕ, ψ by x, y, z: here, r = r+(1−γ)l, z = z+(1−γ)h.

So