Spiral wave drift on the surface of a torus Sergei F. Pravdin Krasovskii Institute of Mathematics and Mechanics (Yekaterinburg, Russia) Ural Federal University (Yekaterinburg, Russia) Hans Dierckx Alexander V. Panfilov Ghent University Ghent University (Ghent, Belgium) (Ghent, Belgium) Abstract Theory of drift of spiral waves has been developed for isotropic and anisotropic plane and isotropic sphere. Here, we proceed to parts of the isotropic torus surface with no-flux boundary conditions and show how velocity components and drift coefficients depend on core position, grid size and local surface curvature. 1 Introduction The theory of drift of spiral waves has been created for isotropic and anisotropic plane and isotropic sphere. The case of isotropical plane was considered in several works devoted to the stationary rotation, meandering and drift of spirals (see for example [3, 2, 12]). Spirals on a spherical surface drift only when the sphere is small enough to cause pairwise interaction [7] or when its radius is modulated in time [1, 8]. Drift of spiral waves on non-uniformly curved surfaces with isotropic diffusion was examined in [1, 8]. Chemical spiral waves on the isotropic sphere were considered in [15, 14]. Symmetry of spirals on the spherical surface was examined in [17]. The anisotropic plane case has been studied motivated by cardiac applications [6, 16, 13, 5, 10]. Finally, a unified theory producing drift laws for generic isotropic and anisotropic curved surfaces were described in [10]. Recently, a theory was also created for a mechanically deformable medium [9]. Here, we study drift of spiral waves on a torus surface with isotropic diffusion. Since the local curvature depends only on one torus coordinate, it is a simple test case for the asymptotic theory developed in [10]. The aim of the present paper is to study drift of spiral waves on toroidal surfaces numerically. We will use the Panfilov–ten Tusscher 2006 model of human ventricular tissue [18] and investigate how drift velocity depends on torus and numerical grid parameters. 2 Methods 2.1 Tissue model We use Panfilov–ten Tusscher model of human ventricular tissue TP06 [18]. The monodomain reaction-diffusion equations describe an homogeneous isotropic active medium: ∂u Iion = D∆u − , ∂t Cm Copyright c by the paper’s authors. Copying permitted for private and academic purposes. In: A.A. Makhnev, S.F. Pravdin (eds.): Proceedings of the International Youth School-conference «SoProMat-2017», Yekaterinburg, Russia, 06-Feb-2017, published at http://ceur-ws.org 260 Iion = IKr + IKs + IK1 + Ito + IN a + IbN a + ICaL + +IbCa + IN aK + IN aCa + IpCa + IpK , where u = u(~r, t) is transmembrane potential of a cell at the point ~r in time t, D is diffusion coefficient and Cm is the electrical capacitance of cell membrane. Intracellular processes are captured by the sum of ionic currents Iion = Iion (~r, t) which depend on 17 additional non-diffusing state variables that describe local ion concentrations and membrane channel states. Parameter values are taken from [18] and correspond a normal cardiomyocyte. 2.2 Torus geometry In the following, we use Cartesian coordinates xi , i ∈ {1, 2, 3}, x1 = x, x2 = y, x3 = z, to describe points in space. On the surface, we use toroidal surface coordinates ξ A , A ∈ {1, 2}. The theoretical analysis then starts from the mapping xi = ri (ξ A ) which defines the chosen geometry in space [10]. For the case of a torus with major radius a and minor radius b < a, we can take ξ 1 = ξ, ξ 2 = η to find the explicit functions ri (ξ A ):     ξ η ξ η ξ x = a + b sin cos , y = a + b sin sin , z = b cos (1) b a b a b where ξ ∈ [0, 2πb] and η ∈ [0, 2πa]. From this relation, one can define the metric tensor gAB = ∂A~r · ∂B ~r, that is used in the Laplacian formula on a curved surface: ∆V = g −1/2 ∂A (g 1/2 g AB ∂B V ) (2) where g AB gBC = δC A , g = det(gAB ). Summation over repeated upper and lower indices is assumed here and below. It then follows that the components of the metric tensor g (i.e. the first fundamental form) of the surface are given by g11 = h21 = 1 g12 = 0 (3) √ g22 = h22 (ξ) g = h1 h2 = h2 (ξ) (4) in terms of the scale factors h1 (ξ) = 1, h2 (ξ) = (1 + (b/a) sin(ξ/b)). (5) We chose the torus example to make h1 , h2 independent of η. Note that near ξ = 0, the chosen coordinates are locally Euclidean. Since the metric is block diagonal, one has g 11 = h−2 1 (ξ) = 1, g 22 = h−2 2 (ξ), g 12 = 0. (6) √ To write the diffusion term, it is convenient to introduce M AB = gg AB , i.e. M 11 (ξ) = h2 (ξ), M 22 (ξ) = 1/h2 (ξ), M 12 (ξ) = 0. (7) The Gaussian curvature K of the torus equals the product of the principal inverse radii of curvature, i.e. 1 sin(ξ/b) K= = . (8) R1 R2 b(a + b sin(ξ/b)) By Gauss’ Theorema Egregium, K also equals half of the Ricci curvature R of the surface. Thus, the Ricci curvature is positive only if sin(ξ/b) > 0, i.e. on the outer half of the torus. On the inner half, R < 0 which can be easily understood since all points there are saddle points which have principal curvatures k1 , k2 of different sign. Hence, R = 2k1 k2 is also negative there. Note that R reaches local maxima at ξ/b = (−1/2 + 2`)π and local minima at ξ/b = (+1/2 + 2`)π, where ` ∈ Z. Gradients of R are known to cause drift of the spiral wave. In [10] it was computed in leading order of a curvature expansion that Ẋ A = −q1 g AB ∂B R − q2 g −1/2 BA ∂B R. (9) where ξ = X(t), η = Y (t) are the surface coordinates of the spiral rotation center and the pair (X, Y ) is denoted as X A . Furthermore, q1 , q2 are coefficients which depend on the chosen reaction-diffusion model. Therefore we compute a cos(ξ/b) ∂ξ R = , ∂η R = 0. (10) b2 (a + b sin(ξ/b))2 261 2.3 Numerical methods From Eq. (2) and the block-diagonal structure of the metric, we infer that 1 1 ∆V = √ ∂ξ M 11 ∂ξ V + √ ∂η M 22 ∂η V .   (11) g g Discretisation using central finite differences on a rectangular grid for ξ, η with space step dξ = dη = dr then delivers X ∆V (ξ, η) ≈ Cm,n (ξ)V (ξ + m dr, η + n dr). (12) m,n∈{−1,0,1} This expression replaces the usual 5-point Laplacian stencil. The coefficients Cm,n (ξ) are only computed at the start of the simulation for every grid point and then stored. For the torus case under consideration, they are: h2 (ξ)[h2 (ξ + dr/2) + h2 (ξ − dr/2)] + 2 C0,0 (ξ) = − , (dr h2 (ξ))2 h2 (ξ ± dr/2) C±1,0 (ξ) = , dr2 h2 (ξ) 1 C0,±1 (ξ) = , [dr h2 (ξ)]2 C±1,±1 (ξ) =0. (13) We used the explicit Euler method for integration over time t. For the simulation domain, we took a part of the toroidal surface that corresponds to a rectangular domain in the surface coordinates ξ, η. At the edges of the domain ξ ∈ [−L/2, L/2], η ∈ [−L/2, L/2], no-flux boundary conditions were imposed. We took L = 200 mm close to the average length of circumference of human left ventricle at its equator. Spiral waves are made using a S1S2 protocol. Calculations are performed in a C program on cluster URAN (IMM UB RAS). We use OpenMP parallelisation technology and GNU C Compiler «gcc». Parameters that are constant among all simulations are the following. Diffusion coefficient D = 0.154 mm2 /ms. In S1S2 stimulation protocol, S1 is applied to the area x < L/5, S2 is applied to the area y < L/2. The spiral wave period was 249 ms on a planar domain. 2.4 Analysis of tip motion The dynamics of spiral waves are usually studied based on trajectories of the points (tips) around which the spiral rotates. To find such points, a level of transmembrane potential u∗ is set, which describes the union of the wave front (which moves in space to more negative u∗ ) and wave back (moving to more positive u∗ ). In two dimensions, the spiral wave tip is the point where the wave front meets the wave back. In [11], it was shown that the tip coordinates ~rtip can be numerically found from simulation results by solving for the intersection of two isolines: u(~rtip , t) = u∗ , u(~rtip , t + ∆t) = u∗ , where u∗ and ∆t have model-specific values. We take u∗ = 0 mV and ∆t = 8 ms (for dt = 0.02 ms, see Table 2) or ∆t = 5 ms (for dt = 0.005 ms). Tip trajectories enable us to find an average drift velocity and the type of spiral wave dynamics. Hereto, we time-average of the spiral tip trajectory to find a curve without self-intersections, which approximates the trajectory of the spiral wave’s rotation center ξ(t), η(t) in surface coordinates. Then, we can locally find the slope of these functions, which define the velocities Vξ , Vη . In components, we have that Vξ (ξ) = −q1 g 11 ∂ξ R(ξ) = −q1 ∂ξ R(ξ) (14) 1 Vη (ξ) = −q2 √ ∂ξ R(ξ). (15) g(ξ) 262 By measuring spiral drift velocity at a given ξ, one can thus try to measure q1 and q2 . Due to discrete sampling in time with interval ∆t, the tip trajectory is found as a sequence of points (ξj , ηj ) at times tj = t0 + j∆t. We cut a starting part of the sequence (at t = t0 ) to omit all initial transient fluctuations and then approximate the trajectory by the functions ξ = pξ0 + pξ1 t + Aξ cos(pξ2 t + pξ3 ), η = pη0 + pη1 t + Aη sin(pη2 t + pη3 ) using MathCAD genfit (least squares method) function. The core’s semiaxes A are estimated as halves of the differences between maximum and minimum of ξi , ηi for 50∆t (250 or 400 ms). We obtain two sets of coefficients pξ and pη where pξ2 = pη2 . The drift velocity has components ~v = (vξ , vη ) = (pξ1 , pη1 ). We study how drift velocity depends on torus parameters, on grid size and on spiral core coordinates. For each set of torus and grid parameters, we use 5 moments of time to apply S2 stimulus (see Table 1), to create spiral waves located at different ξ. Table 1: Moments of applying S2 stimulus Simulation S2 time number 1 310 2 345 3 380 4 415 5 450 Torus and grid parameters are given in Table 2. We used spatial step 0.1 mm, which is 2.5 times less than the usually used step for this model. We hope that such a small step, together with two other steps we try, will help us to get more precise results. Table 2: Grid and torus parameters in simulation series Simulation Grid size Spatial step Timestep Simulation Torus parameters series n dr, mm dt, ms time, s a, mm b, mm A 500 0.4 0.02 40 127.4 63.7 B 500 0.4 0.02 40 254.8 127.4 C 500 0.4 0.02 40 509.6 254.8 D 1000 0.2 0.02 40 254.8 127.4 E 1000 0.2 0.02 40 127.4 63.7 F 2000 0.1 0.005 20 254.8 127.4 G 2000 0.1 0.005 20 127.4 63.7 3 Results First of all we checked spiral wave dynamics on a plane isotropic square. We observed an almost rigid rotation on all three values of gridsize dr (the velocity components were less than 0.005 mm/sec). For the torus surfaces, velocity components are shown in Fig. 1. That Figure shows that vξ is proportional to vη , confirming that drift occurs under a fixed angle with the gradient (here, the ξ-direction). Figure 2 displays a typical trajectory of the spiral wave tip. We see an oval-shaped curve slowly shifting in some direction. Our task was to isolate the slow low-frequency drift and to exclude the high-frequency core rotation. Figure 3 shows plots of the averaged velocity components on [t0 , t] where t is the abscissa. We still can see some high-frequency oscillations and a limit for each velocity component. Figure 4 shows how velocity component Vξ depends on ∂ξ R. We see that Vξ is positive and decreases with a growth of positive curvature ∂ξ R. Figure 5 shows how velocity component Vη depends on ∂ξ R. There is a weak change of positive Vη when ∂ξ R grows. Finally, Figure 6 displays law of change of coefficient q1 on ∂ξ R as nearly hyperbolic. 263 Figure 1: Velocity of the spiral wave drift in all simulation series. Logarithmic scale. Figure 2: Simulation E1. Simulation time 40 sec. Plot of the tip trajectory after removing the initial part. Plot coordinates vary from 0 to n dr and they are shifted by n dr/2 relatively to ξ and η, which vary from −n dr/2 to n dr/2. 264 Figure 3: Simulation E1. Plot of the velocity components vξ (on the left), vη (on the right) against time, sec. The limits limt→+∞ vξ and limt→+∞ vη are approximately taken as v at the full simulation time. Figure 4: Simulation series F and G. Plot of velocity Figure 5: Simulation series F and G. Plot of velocity component Vξ against ∂ξ R. component Vη against ∂ξ R. 265 Figure 6: Simulation series F and G. Plot of coefficient q1 against ∂ξ R. 4 Discussion and conclusions We present the first simulations on drift of spiral waves in an ionic model of cardiac tissue on a surface with with a non-zero gradient of the Ricci curvature scalar R. The linear theory of the spiral wave drift predicts that both velocity components of drift are directly proportional to the gradient of R and should become zero at ∂ξ R = 0. This is not the case in our simulations. We see more complex dependency of the drift velocity on the R gradient. There are several possible reasons for such discrepancy. First, it potentially can be due to numerical issues. The most crucial parameter in the integration of ionic models is spatial discretisation step. We used a spatial step of 0.1 mm which is widely accepted as sufficient, as typical simulations employ a spatial resolution of between 0.1 and 0.2 mm [4]. Moreover, as cardiac cells have a typical length of 0.1 mm, lower spatial resolutions are not relevant for studies of macroscopic wave propagation in the heart. Thus, it is unlikely that the observed discrepancy can be due to numerical arguments. Next, the deviations can be due to the fact that linear theory may be insufficient. This can also be caused by several factors. First of all, on the considered surface, the gradient of Ricci curvature scalar is not constant and thus higher order terms can be essential. In addition, the linear theory is based on a theory of response functions [10] which in turn relies on stationary rotation of spiral waves. In our case, we have a strong front-tail interaction during spiral wave rotation causing non-stationary dynamics and thus the response functions theory should be adjusted for such case. Therefore additional studies are necessary to find the limits of application of linear theory. All our simulations show that the instantaneous value of coefficient q1 is negative. It indicates that spiral drifts opposite to the gradient of R. In case of an isotropic surface, it means that a spiral wave in TP06 model will drift towards a larger Gauss curvature of the surface. Note, however, that in some cases, we observed that the spiral waves stabilises at points which are not maxima of Gauss curvature. That once again may indicate a possible non-linear interaction of curvature of the surface with spiral wave drift dynamics. Acknowledgements We performed our simulations on the clusters of Krasovskii Institute of Mathematics and Mechanics («URAN») and Ural Federal University. 266 References [1] A.Yu. Abramychev, V.A. Davydov, and V.S. Zykov. Drift of spiral waves on nonuniformly curved surfaces. Zh. Eksp. Teor. Fiz., 97:1188–1197, 1990. [2] V. N. Biktashev, A. V. Holden, and E. V. Nikolaev. Spiral wave meander and symmetry of the plane. International Journal of Bifurcation and Chaos, 06(12a):2433–2440, 1996. [3] V.N. Biktashev and A.V. Holden. Resonant drift of autowave vortices in two dimensions and the effects of boundaries and inhomogeneities. Chaos, Solitons and Fractals, 5(3):575 – 622, 1995. [4] R.H. Clayton, O. Bernus, E.M. Cherry, H. Dierckx, F.H. Fenton, L. Mirabella, A.V. Panfilov, F.B. Sachse, G. Seemann, and H. Zhang. Models of cardiac tissue electrophysiology: Progress, challenges and open questions. Progress in Biophysics and Molecular Biology, 104(1–3):22–48, 2011. [5] V.A. Davydov, V.G. Morozov, N.V. Davydov, and T. Yamaguchi. Propagation of autowaves in excitable media with chiral anisotropy. Physics Letters A, 325(5-6):334–339, 2004. [6] V.A. Davydov and V.S. Zykov. Spiral waves in anisotropic excitable media. Zh. Eksp. Teor. Fiz., 95:139–148, 1988. [7] V.A. Davydov and V.S. Zykov. Kinematics of spiral waves on nonuniformly curved surfaces. Physica D, 49:71–74, 1991. [8] V.A. Davydov, V.S. Zykov, and T. Yamaguchi. Drift of spiral waves on nonuniformly curved surfaces. Macromol. Symp., 160:99–106, 2000. [9] Hans Dierckx, Sander Arens, Bing-Wei Li, Louis D Weise, and Alexander V Panfilov. A theory for spiral wave drift in reaction-diffusion-mechanics systems. New Journal of Physics, 17(4):043055, 2015. [10] Hans Dierckx, Evelien Brisard, Henri Verschelde, and Alexander V. Panfilov. Drift laws for spiral waves on curved anisotropic surfaces. Phys. Rev. E, 88:012908, Jul 2013. [11] Flavio Fenton and Alain Karma. Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation. Chaos: An Interdisciplinary Journal of Nonlinear Science, 8(1):20–47, 1998. [12] James P. Keener. A geometrical theory for spiral waves in excitable media. SIAM Journal on Applied Mathematics, 46(6):1039–1056, 1986. [13] Victor G LeBlanc. Rotational symmetry breaking for spiral waves. Nonlinearity, 15(4):1179, 2002. [14] Jerzy Maselko. Symmetrical double rotor spiral waves on spherical surfaces. J. Chem. Soc., Faraday Trans., 94:2343–2345, 1998. [15] Paul McQuillan and Jagannathan Gomatam. Rotating chemical waves on the sphere. The Journal of Physical Chemistry, 100(13):5157–5159, 1996. [16] J.M. Rogers and A.D. McCulloch. Nonuniform muscle fiber orientation causes spiral wave drift in a finite element model of cardiac ation potential propagation. J Cardiovasc Electrophysiol, 5:496–509, 1994. [17] Rachel Sigrist and Paul Matthews. Symmetric spiral patterns on spheres. SIAM Journal on Applied Dynamical Systems, 10(3):1177–1211, 2011. [18] K.H. ten Tusscher and A.V. Panfilov. Alternans and spiral breakup in a human ventricular tissue model. Am. J. Physiol. Heart Circ. Physiol., 291:H1088–1100, 2006. 267