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.
2.1
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 = D u
Iion ;
Cm 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
In the following, we use Cartesian coordinates xi, i 2 f1; 2; 3g, x1 = x, x2 = y, x3 = z, to describe points in space. On the surface, we use toroidal surface coordinates A, A 2 f1; 2g. 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 b cos ; a y = a + b sin b sin ; a z = b cos b where 2 [0; 2 b] and 2 [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: where gABgBC = CA, 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 To write the diffusion term, it is convenient to introduce M AB = pggAB, i.e.
M 11( ) = h2( );
M 22( ) = 1=h2( );
;
(
K =
1 R1R2 =
sin( =b) 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 = 2k1k2 is also negative there. Note that R reaches local maxima at =b = ( 1=2 + 2`) and local minima at =b = (+1=2 + 2`) , where ` 2 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
X_ A = q1gAB@BR where = X(t), = Y (t) are the surface coordinates of the spiral rotation center and the pair (X; Y ) is denoted as XA. Furthermore, q1, q2 are coefficients which depend on the chosen reaction-diffusion model. Therefore we compute
From Eq. (
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 2 [ L=2; L=2], 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
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
X Discretisation using central finite differences on a rectangular grid for ; delivers with space step d = d = dr then 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:
C0;0( ) = C 1;0( ) =
C0; 1( ) = C 1; 1( ) =0: V ( ) =
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 = p0 + p1t + A cos(p2t + p3); = p0 + p1 t + A sin(p2 t + p3 ) 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 p2 = p2 . The drift velocity has components ~v = (v ; v ) = (p1; p1 ).
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 .
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. 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. 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.
We performed our simulations on the clusters of Krasovskii Institute of Mathematics and Mechanics (¾URAN¿) and Ural Federal University.