Modelling far field pacing for terminating spiral waves pinned to ischaemic heterogeneities in cardiac tissue

In cardiac tissue, electrical spiral waves pinned to a heterogeneity can be unpinned (and eventually terminated) using electric far field pulses and recruiting the heterogeneity as a virtual electrode. While for isotropic media the process of unpinning is much better understood, the case of an anisotropic substrate with different conductivities in different directions still needs intensive investigation. To study the impact of anisotropy on the unpinning process, we present numerical simulations based on the bidomain formulation of the phase I of the Luo and Rudy action potential model modified due to the occurrence of acute myocardial ischaemia. Simulating a rotating spiral wave pinned to an ischaemic heterogeneity, we compare the success of sequences of far field pulses in the isotropic and the anisotropic case for spirals still in transient or in steady rotation states. Our results clearly indicate that the range of pacing parameters resulting in successful termination of pinned spiral waves is larger in anisotropic tissue than in an isotropic medium. This article is part of the themed issue ‘Mathematical methods in medicine: neuroscience, cardiology and pathology’.

In cardiac tissue, electrical spiral waves pinned to a heterogeneity can be unpinned (and eventually terminated) using electric far field pulses and recruiting the heterogeneity as a virtual electrode. While for isotropic media the process of unpinning is much better understood, the case of an anisotropic substrate with different conductivities in different directions still needs intensive investigation. To study the impact of anisotropy on the unpinning process, we present numerical simulations based on the bidomain formulation of the phase I of the Luo and Rudy action potential model modified due to the occurrence of acute myocardial ischaemia. Simulating a rotating spiral wave pinned to an ischaemic heterogeneity, we compare the success of sequences of far field pulses in the isotropic and the anisotropic case for spirals still in transient or in steady rotation states. Our results 2017 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Introduction
Cardiac tissue is a nonlinear excitable medium whose spatio-temporal dynamics exhibits a complex behaviour ranging from regular to turbulent activity. The first is associated with plane waves, generated by specific pacemaker cells and travelling through the tissue triggering a coordinated mechanical contraction of the myocardium. Turbulent activity, instead, arises due to dysfunctions affecting the propagation of the electric pulse and seriously affects the pumping function of the heart, leading to the occurrence of life-threatening diseases, i.e. cardiac arrhythmias. Asynchronous electromechanical activity in the heart is known in the literature as fibrillation. During fibrillation, the regular heart contraction is altered by fast re-entrant waves (spiral or scroll waves) [1][2][3][4], which find a very critical substrate in the heterogeneous and anisotropic structure of the cardiac tissue [5]. Spiral and scroll waves can in fact get anchored (pinned) to localized heterogeneities, giving rise to abnormal high-frequency heart rhythms leading to the aforementioned arrhythmias. For this reason, if not followed by repinning, unpinning is an important step towards successful termination of re-entries.
In this scenario, one of the most challenging goals is controlling the complex dynamics of the heart under conditions of fibrillation. In clinical practice, the most often used and effective approach against re-entrant activity in the heart consists of delivering a high-energy electric shock that excites the entire medium at once and sets all the cells to their refractory state, thus terminating any electrical activity. Although this defibrillation technique is used routinely, the delivery of high-energy pulses is associated with undesired side effects that seriously affect the successful outcome of the therapy [6][7][8].
For this reason, there is a strong need for low-energy control methods. The approach presented in this paper is based on far field pacing (FFP), which recruits both the intrinsic (due to the complex geometry) and the induced (due to pathological conditions) multi-sized heterogeneities of the heart tissue to successfully terminate arrhythmias [9]. A heterogeneity represents an area of conduction changes, and when an electric field is applied to the whole heart, areas of depolarization and hyperpolarization (so-called Weidmann zones) [10] appear at the boundary of this 'obstacle' due to a redistribution of intracellular and extracellular currents. If the depolarization exceeds the excitation threshold, the heterogeneity acts as a source of excitation, i.e. a virtual electrode [11][12][13][14][15][16][17]. The intensity of the external electric field applied in FFP is much weaker than the defibrillation threshold [18]. It has been demonstrated both numerically and experimentally that waves emitted by virtual electrodes can be used to unpin or terminate re-entrant waves pinned to heterogeneities [19][20][21].
Detailed studies have been conducted in order to better understand these mechanisms and their impact on the success rate of FFP applied to isotropic media [12,19,22]. These investigations introduced concepts such as vulnerability and vulnerable window [19,23,24] and their relationship with the unpinning window, and refractory and excitable periods [19] on the obstacle's boundary to which the spiral is pinned.
Intensive investigation is still needed to translate these mechanisms in anisotropic media, where the role played by fibres and different conductivities in different directions has been demonstrated to be fundamental in order to obtain reliable simulation results [25,26]. This paper focuses on the effects exerted by anisotropy on the process of unpinning and termination of pinned spirals. We implement a bidomain formulation of phase I of the Luo  domains, separated by the cellular membrane. Therefore, it is the most realistic mathematical expression for macroscopic simulation of cardiac tissue [27]. The two spaces have different electrical conductivities and in each domain anisotropy is reproduced by assigning different conductivity values in the directions parallel and perpendicular to myocardial fibres. Termination of a spiral pinned to an ischaemic heterogeneity is attempted through the delivery of sequences of pulses in FFP at different pacing frequencies and at different angles of the pinned spiral rotating around the heterogeneity. The paper is organized as follows: §2 introduces the modelling used in this article; § §3 and 4 report the results of FFP applied in an isotropic and in an anisotropic medium, respectively; §5 sums up and compares the outcomes of the two previous sets of simulations.

Simulations
The impact of anisotropy on wave dynamics is investigated by means of a bidomain model of cardiac tissue [27].
The transmembrane potential V m and the extracellular potential V e are governed by where β is the ratio of the membrane surface area to the tissue volume, C m is the membrane capacitance per unit area, σ e and σ i are the diagonal extracellular and intracellular conductivity tensors, respectively (each one resolved into two components: one parallel (σ L ) and one perpendicular (σ T ) to the fibres), and I ion is the sum of the transmembrane ionic currents described by the local ionic model, being in our sets of simulations phase I of the Luo and Rudy action potential model [28].
Simulations are run in a two-dimensional sheet of cardiac tissue (4 cm × 4 cm) with a circular ischaemic heterogeneity (radius 0.6 cm) in the middle of it. In order to simulate acute ischaemia, the ionic model is modified accordingly. The ischaemic area is treated as a region with reduced conduction properties where conditions of hyperkalaemia and acidosis arising at the cellular level are modelled. The extracellular potassium concentration, [K + ] o , is incremented from 5.4 to 14 mmol l −1 [29] due to hyperkalaemia; the sodium and L-type calcium channel-specific conductances, g Na and g Ca , respectively (both chosen according to Luo & Rudy [28]), are both reduced by 25% [30] as an effect of acidosis. Changes in intracellular and extracellular conductivities are simulated according to the experimental data provided by Kléber & Riegger [31]. After the arrest of coronary flow, extracellular resistance is supposed to increase immediately by approximately 30%, due to a drop in perfusion pressure and a diminution of intravascular volume. After approximately 20 min, a second gradual increase up to 80% occurs. Regarding the intracellular resistance, it is supposed to remain constant in the first 10-12 min, followed by a steep increase due to the progression of cellular uncoupling and the breakdown of cellular homeostasis. Considering these phenomena and assuming that the dynamics of the system is simulated approximately 20 min after the onset of ischaemia, bidomain conductivities in the ischaemic heterogeneity are reduced by 40% in the extracellular domain and by 50% in the intracellular domain [31].
Each pulse is delivered for 5 ms in FFP by plane electrodes placed perpendicular to the longitudinal direction (x-axis = direction of the applied electric field). During the pulses, the boundary conditions are imposed according to Fishler [15]: for a longitudinally oriented shock, in both the intracellular and the extracellular domains potentials on the transverse boundary are governed by Neumann no-flux boundary conditions, while potentials along the longitudinal boundaries are governed by Dirichlet boundary conditions. In particular for Dirichlet boundary

conditions,
where d is the edge length of the two-dimensional simulated tissue. Equation (2.3) ensures the application of a uniform electric field with intensity E 0 and direction parallel to the fibres. Using condition (2.4) avoids the near-field effects on V m at the tissue edge and thus allows the focus to be only on the impact of far field stimulations. During intervals between pulses, instead, Neumann no-flux conditions are applied at all boundaries. The set of parameters adopted in both the isotropic and the anisotropic computations are reported in table 1.

Isotropic medium
In this section, we present one simplifying assumption often used in bidomain models of cardiac tissue: the 'isotropic bidomain' approximation [32]. It affects the transversal and longitudinal conductivities that are set to be the same in both the intracellular and extracellular spaces, i.e. σ iT = σ iL = σ i and σ eT = σ eL = σ e . This assumption removes the presence of fibres from the model of the tissue.
Assigning conductivity values in bidomain model simulations is not straightforward: although these parameters have been measured experimentally [33][34][35][36][37], the data are inconsistent and therefore no consensus exists regarding their accurate values. In both isotropic and anisotropic media, our simulations refer to the experimental values provided by Clerc [33]. Figure 1 shows snapshots taken from a simulation computing the rotation of a spiral pinned to the ischaemic heterogeneity (highlighted by the white dashed circle). In the absence of external stimulation, the spiral is stably pinned to the obstacle and a rotation period (T sp ) of 380 ms is recorded. The positions of the spiral tip on the border of the ischaemic area are shown at 0 ms (figure 1a), 100 ms (figure 1b), 200 ms (figure 1c) and 340 ms (figure 1d), within one rotation period. The panels in figure 1 show the states of the spiral taken during an early phase of the pinning process, in which the spiral shows a stable behaviour in terms of constant speed and persistence of pinning in the absence of external perturbation, but has not yet spent a 'sufficiently long time' in this condition to be in its full steady rotation. The motivation behind the choice of a transient behaviour will be clarified further in the text. For the stimulation protocol, different times for delivering the first pulse (t 1 ) were chosen, corresponding to different angles of the spiral tip on the boundary of the ischaemic zone. This means that the positions of the spiral tip were chosen at an interval of 20 ms within one rotation period and each of them was set as the initial condition for a set of simulations. Figure 1 provides four examples of initial conditions at which the first pulse was delivered. In each set of simulations, the interval between pulses (T p ) was set taking into account two approaches of electrical stimulation: underdrive 1 and overdrive 2 pacing. Therefore, according to the first method, The number of delivered pulses ranged from one to five and the outcome of a stimulation (maintenance, unpinning or termination of the pinned spiral) was checked at each interval immediately before the potential application of another pulse. At this point, the pacing protocol was considered to be successful and thus interrupted only in the case of spiral termination.
The intensity of the electric field, E 0 , represents the minimum intensity which, under isotropic conditions, induces a depolarization that exceeds the excitation threshold at the boundary of the obstacle and creates a propagating wave. Figures 2 and 3 show two interesting phenomena arising in an isotropic tissue during FFP stimulation: unpinning and termination of a pinned spiral, respectively. Figure 2 shows snapshots taken from a simulation implementing overdrive pacing (T p = 0.9T sp ) with t 1 = 360 ms, after the delivery of the third pulse. In figure 2a, the application of the electric field induces areas of depolarization and hyperpolarization at the boundaries of the ischaemic zone. Since the excitation threshold is exceeded, a new wave (N) is generated. One 1 Electrical stimulation of the heart at a rate lower than that of an existing tachycardia. 2 Electrical stimulation of the heart at a rate faster than its intrinsic rate.    On the x-axis different time intervals between pulses (T p ) relative to the period of the spiral (T sp ) are shown, varying from overdrive to underdrive pacing. The y-axis reports several initial conditions at which the first pulse was delivered (t 1 ), each one corresponding to a specific position of the spiral on the border of the heterogeneity within one rotation period (cf. figure 1). The outcome is expressed in terms of the number of pulses that have to be applied in order to terminate the pinned spiral (as shown in figure 3f ). The maximum number of stimulations delivered during the pacing protocol was five, so the white area in the graph corresponds to sequences that were not successful within the five pulses. For this reason, no statement can be made about the possibility of terminating the spiral wave in that region.  The results obtained from simulations in the isotropic case are summarized in figure 4. The graph shows the success rate for sequences of far field pulses delivered at a certain interval between pulses (on the x-axis, varying from overdrive to underdrive pacing) and starting from a certain initial condition (on the y-axis, where each value corresponds to a specific position of the spiral on the border of the heterogeneity within one rotation period; cf. figure 1). The outcome is expressed in terms of the number of pulses that have to be applied in order to terminate the pinned spiral (as shown in figure 3f ). As already mentioned, the maximum number of stimulations delivered during the pacing protocol was five, so the white area in the graph corresponds to sequences that were not successful within the five pulses. For this reason, no statement can be made regarding the number of pulses required for that region. Regarding successful pulses, it is worth noting that at least three shocks have to be applied in order to terminate the pinned spiral. Moreover, the area of the graph where successful pulse sequences take place is restricted to the 'underdrive region' [22,38].

Anisotropic medium
In this section, we present the results for unpinning and termination of a pinned spiral in an anisotropic medium. In this case no simplifying assumption was made for the bidomain model, but four different conductivities were introduced, as shown in table 1. Myocardial fibres are modelled as straight, uniform and directed along the x-direction. Figure 5 shows snapshots taken from a simulation implementing the rotation of the spiral pinned to the ischaemic heterogeneity without any far field stimulations. Under these conditions, in which anisotropy ratios in the intracellular and extracellular domains are approximately 9 and 2.6 [33], respectively, the rotation period (T sp ) of the spiral changes to 350 ms. The positions of the spiral tip on the boundary of the ischaemic area are recorded at 80 ms (figure 5a), 140 ms (figure 5b), 250 ms (figure 5c) and 300 ms (figure 5d) within one rotation period.
As already mentioned for the isotropic domain, these snapshots refer to the spiral during its transient behaviour. The spiral has a constant period and a 'sufficiently strong' pinning to the heterogeneity, requiring external perturbations in order to be detached from it. The pacing protocol adopted in this section was the same as the one performed for isotropic computations and described in §3. Each set of simulations had initial conditions chosen at different positions of the spiral tip on the heterogeneity, which correspond to the time at which the first pulse was delivered (t 1 ). Again, the difference between the t 1 values was set to 20 ms, in order to cover a wide range of spiral tip positions in the application of the first stimulation. The interval between pulses ranged from overdrive to underdrive pacing.
The intensity of the electric field E 0 is different from the isotropic simulations, since this time a smaller value was already sufficient to induce a propagating wave out of the obstacle.
In contrast with the results shown in §3, anisotropic simulations resulted either in the maintenance or in the termination of the spiral, but no episodes of unpinning were observed.
Snapshots from a simulation resulting in successful termination of the pinned spiral after the second stimulus are shown in figure 6. The overdrive pacing protocol was applied, with (T p = 0.9T sp ) and t 1 = 300 ms. In figure 6a, the applied far field pulse is strong enough to exceed the excitation threshold and to induce a new wave nucleating from the heterogeneity. Since part of the boundary is still refractory due to the tail of the original spiral, only a small portion of the new wave is able to propagate. Within this part, one end disconnects from the obstacle because its propagation is inhibited by the refractory tail of the spiral, as shown in figure 6b. The other end forms a pinned wave and collides with the tip of the original spiral ( figure 6c,d). An overview of the results obtained from the anisotropic simulations with respect to the time interval between pulses and the time of the first pulse t 1 within one spiral rotation is shown in figure 7. The plot reports the success rate of the pacing protocols in terms of the number of pulses (ranging from 1 to 5) that it is necessary to deliver in order to achieve the termination of the pinned On the x-axis different time intervals between pulses (T p ) relative to the period of the spiral (T sp ) are shown, varying from overdrive to underdrive pacing. The y-axis reports several initial conditions at which the first pulse was delivered (t 1 ), each one corresponding to a specific position of the spiral on the border of the heterogeneity within one rotation period; compare with figure 5. The outcome is expressed in terms of the number of pulses that have to be applied in order to terminate the pinned spiral (as shown in figure 6f ). The maximum number of stimulations delivered during the pacing protocol was five, so the white area in the graph corresponds to sequences that were not successful within the five pulses. For this reason, no statement can be made about the possibility of terminating the spiral wave in that region. (Online version in colour.) spiral. As in the previous case, on the x-axis the time interval spans both the overdrive and the underdrive pacing protocols; on the y-axis a set of 18 initial conditions covers a wide range of positions taken by the pinned spiral tip during one rotation (cf. figure 5). In this case, successful termination is already obtained after the delivery of only one pulse for most positions of the spiral tip on the boundary of the obstacle. In addition, most of the combinations of parameters result in spiral termination within the delivery of five pulses (the white region, in which no statement can be made, is much smaller than that in figure 4). The most successful episodes in the upper graph, i.e. terminations requiring a smaller number of pulses, tend to gather in the underdrive region.
All the results shown so far for the anisotropic case refer to the transient behaviour of the pinned spiral. Further computations have been carried out in order to investigate the impact that the different (transient or asymptotic) states of the spiral dynamics can have on the outcome of the pulses. For this reason, the same pacing protocol previously described was also applied during the asymptotic regime of the spiral and the success rate was computed. The success rate of the protocol in the new dynamic state is shown in figure 9. As in figures 4 and 7, the results summarized in figure 9 report the number of pulses (ranging from one to a maximum of five) that enable termination of the spiral. The x-axis is partitioned into overdrive   figure 9 where the first pulse is already sufficient to terminate the pinned spiral. The main difference consists in the extension of this area, which in figure 7 spreads from the underdrive to the overdrive region and covers a wide range of parameter combinations, while in figure 9 it 'shrinks' to a much smaller area in the underdrive region. Note that, in both figure 7 and figure 9, the most successful cases occur for underdrive pacing. In addition, figure 9 shows that in most cases five pulses (or fewer) are sufficient to terminate a steadily rotating spiral.

Conclusion
In this paper, we demonstrate the crucial role of anisotropy in the process of termination of a spiral pinned to an ischaemic heterogeneity during far field stimulation. A comparison between results obtained under isotropic and anisotropic conditions highlights the fundamental differences that must be taken into account when comparing simulations with experimental results. Simplifying assumptions in the conductivities of the bidomain model (as the 'isotropic bidomain' approximation) drastically reduce the computational load of simulations, but have been demonstrated both analytically and numerically to be in disagreement with experimental observations [26]. Therefore, unequal anisotropy ratios must be used in order to achieve reliable results [26,39]. T p /T sp Figure 9. Results obtained from simulations in the anisotropic case during an asymptotic regime. On the x-axis different time intervals between pulses (T p ) relative to the period of the spiral (T sp ) are shown, varying from overdrive to underdrive pacing. The y-axis reports several initial conditions at which the first pulse was delivered (t 1 ), each one corresponding to a specific position of the spiral on the border of the heterogeneity within one rotation period (cf. figure 8). The maximum number of stimulations delivered during the pacing protocol was five, so the white area in the graph corresponds to sequences that were not successful within the five pulses. For this reason, no statement can be made about the possibility of terminating the spiral wave in that region. (Online version in colour.) Furthermore, we compared the impact of the external field pulse on spirals at an early stage when they were still in the transient regime (figures 4 and 7) with spirals that had already reached the asymptotic state of steady rotations (figure 9). This comparison was motivated by clinical practice, when it is important to have rapid and immediate intervention, because waiting too long before defibrillating the heart can result in the death of the patient. Therefore, with the aim of investigating the critical role played by the 'waiting time' and the dynamic state of the heart tissue during defibrillation, figure 9 has been computed to measure the success rate of defibrillation attempts during steady rotation of the pinned spiral.
In our simulations, the main differences between isotropic and anisotropic simulations lie in the amount of energy required to recruit the heterogeneity as a virtual electrode and in the probability of successful termination assigned to a specific sequence of pulses.
In more detail, the intensity of the electric field applied at the virtual electrodes required to exceed the excitation threshold and induce a propagating wave at the boundary of the ischaemic area was five times stronger in the isotropic case than in the one used in the anisotropic case due to different conductivity values. Interestingly, we found that, although a higher amount of energy was released in the isotropic regime, the probability of successful termination was lower. This statement is supported by the comparison between figures 4, 7 and 9.
During the transient regime, in the isotropic case, a higher number of pulses have to be delivered in order to induce termination of the pinned spiral. As shown in figure 4, a sequence of pulses starts to be successful only after the delivery of the third stimulation, while the outcome from anisotropic simulations shows a region where the first pulse is already sufficient to achieve the goal ( figure 7). A smaller intensity of the electric field in anisotropic tissues also ensures spiral termination within the five pulses for a larger range of parameter combinations, as shown by the  figure 7. Checking the outcome of a stimulation (maintenance, unpinning or termination of the pinned spiral) at each interval immediately before the potential application of another pulse allowed the protocol to be stopped as soon as termination was detected. In this way, it was possible to avoid the delivery of unnecessary pulses. In this last case, one would expect that, in the case of completely quiescent tissue, the delivery of a new pulse and the recruitment of the heterogeneity as a virtual electrode should not alter the success rate, since the new wave would spread in the tissue without hitting any further obstacles. The results could be very different in the case of existing residual waves propagating in the tissue at the moment when a new far field pulse is applied and a new wave originates from the heterogeneity. In these circumstances, the collision between the waves might lead to a more complex dynamics that might affect the success rate of the defibrillation attempts. When translating this pacing protocol to clinical practice (in experiments or in modern pacemakers, for example), one may think of regulating the electrical behaviour of the heart by sensing and stimulating activities using specific electrodes implanted at different points in the heart. When abnormal activity due to incoming arrhythmias is detected by the sensing electrodes, a variable number of external pulses could be delivered by the stimulating electrodes until the physiological sinus rhythm is restored and the defibrillation attempt is considered successful.
The results from our simulations show that successful termination in anisotropic media is not restricted to underdrive pacing, as in the case of isotropic substrates [22], but it also extends to overdrive pacing. The exact mechanism that underlies this behaviour in the anisotropic computations still has to be completely understood. A preliminary guess could be made by comparing figures 7 and 9, i.e. transient and steady rotation of the pinned spiral. Our results show that the transient regime seems somehow to support the process of spiral termination, since the spiral dynamics is less stable than in its steady state and, for this reason, perturbations by external far field pulses are more likely to result in spiral termination. This 'vulnerability' of the spiral, together with the higher conductivity values in the longitudinal direction, i.e. the direction of the fibres, might induce a major drift of the spirals towards the boundary of the tissue and cause their annihilation. However, as previously stated, this interpretation is only preliminary and would require additional simulations to explore a broader range of possible scenarios.
Another interesting observation that we made was the occurrence of unpinning and repinning of the spiral only in isotropic simulations. These phenomena have been observed in previous isotropic investigations, where the concepts of unpinning and vulnerable window have been studied [18,19]. Further research is needed to properly translate the mechanisms of re-entrant dynamics into anisotropic computations. This includes the impact of the relative angle between the direction of the electric field and the orientation of the myocardial fibres. In this study, as a first approach, an electric field parallel to myocardial fibres has been modelled, but we expect the relative angle to be a 'critical parameter'. To summarize, our results provide clear evidence of the strong impact exerted by anisotropy in the tissue and of misleading results that can be obtained when introducing any simplifying approximations for the conductivity values in bidomain modelling.