In silico optical control of pinned electrical vortices in an excitable biological medium

Vortices of excitation are generic to any complex excitable system. In the heart, they occur as rotors, spirals (2D) and scroll waves (3D) of electrical activity that are associated with rhythm disorders, known as arrhythmias. Lethal cardiac arrhythmias often result in sudden death, which is one of the leading causes of mortality in the industrialized world. Irrespective of the nature of the excitable medium, the rotation of a rotor is driven by its dynamics at the (vortex) core. In a recent study, Majumder et al (2018 eLife 7 e41076) demonstrated, using in silico and in vitro cardiac optogenetics, that light-guided manipulation of the core of free rotors can be used to establish real-time spatiotemporal control over the position, number and rotation of these rotors in cardiac tissue. Strategic application of this method, called ‘Attract-Anchor-Drag’ (AAD) can also be used to eliminate free rotors from the heart and stop cardiac arrhythmias. However, rotors in excitable systems, can pin (anchor) around local heterogeneities as well, thereby limiting their dynamics and possibility for spatial control. Here, we expand our results and numerically demonstrate, that AAD method can also detach anchored vortices from inhomogeneities and subsequently control their dynamics in excitable systems. Thus, overall we demonstrate that AAD control is one of the first universal methods that can be applied to both free and pinned vortices, to ensure their spatial control and removal from the heart and, possibly, other excitable systems.


Introduction
Vortices, in the form of self-organized rotors of spatial excitation, occur in various complex excitable systems [1]. These include, e.g. active galaxies [2], oscillatory chemical reactions [3,4], social amoebae Dictyostelium discoideum [5], Carbon monoxide oxidation on a Platinum surface [6], and cardiac tissue [7]. In the heart, such vortices are considered to be detrimental, as they underlie lethal cardiac arrhythmias.
The important features of the vortex, such as, the rotation frequency, spatial organization of emitted waves are determined by a core of a vortex, which is characterized by a phase singularity point, or a closed loop of phase singularities. It was shown that real-time manipulation of the core, could result in direct spatiotemporal control over these vortices in various excitable systems [8][9][10]. However, these studies focused on controlling vortices that existed freely in the medium, without attaching to a physical inhomogeneity. The validity of such control methods has not been tested in systems with anchored vortices, i.e. where the vortex attaches itself to a large physical inhomogeneity and starts circulating around it. Such pinned activity can occur in various excitable systems, see e.g. [11,12].
In the heart, both free and anchored vortices commonly occur. It is widely accepted that anchored vortices organize so called monomorphic tachycardias, characterized by periodic electrocardiograms (ECG), whereas free vortices, together with their complex spatial organizations, underlie polymorphic tachycardias and fibrillation, which are characterized by complex non-periodic ECGs. As these cardiac arrhythmias are serious medical conditions that, if left untreated, can lead to sudden death, it is very important to develop an understanding of the dynamics of these associated vortices and use the acquired knowledge to design methods for arrhythmia control. From a general point of view, anchored vortices are more difficult to control, as compared to free ones, because of the extra stability imposed upon these waves by the inhomogeneity, as it serves as the core organizing their dynamics. Therefore the methods developed to control free vortices [9,[13][14][15][16][17][18][19][20][21][22][23][24] can be ineffective for pinned vortices [25]. For pinned vortices, unpinning is effectuated by the application of external electrical field [26][27][28][29][30][31][32][33], or by a rotating electrical pulse [34]. One way to overcome this issue is to detach (unpin) the vortex from the inhomogeneity and then possibly use one of the mentioned methods. It was shown that a weak shock, with an amplitude of an order of magnitude less than the defibrillating shock, may unpin the vortices rotating around the obstacle, [35]. Other methods include regional cooling of the excitable medium [36], or application of periodic mechanical deformation [37]. However, in any of the afore-mentoned cases, the unpinning methdology by itself does not eliminates the vortex, while just unpinning that is not followed by attempts to eliminate the vortices, and can have dangerous consequences leading to more complex arrhythmia.
Recently Majumderet al [8] proposed a novel method to control spiral waves in cardiac tissue. They demonstrated realtime spatiotemporal control over spiral wave dynamics by pinning-based relocation of vortices in optogenetically-modified monolayers of neonatal rat cardiac cells. Their method (Attract-Anchor-Drag-control) was shown to be extremely powerful for the elimination of single or multiple spiral waves in twodimensional (2D), homogeneous cardiac tissue systems. Here we apply this method to pinned vortices and show that using 'Attract-Anchor-Drag' (AAD)-control method, it is possible to detach pinned spiral and scroll waves from conduction-type inhomogeneities over a wide range of sizes, and then remove these vortices from the tissue using the same approach. We investigate conditions necessary for the success of this phenomenon, study the dependence of its success rate on the characteristics of the pinned vortex and propose a modofication of the method which further increases its efficiency. Overall we show that AAD not only allows unpinning, but also removal of anchored vortices, making it one of the first universal control methods, that can be applied both for pinned, as well as free vortices in the heart and other excitable systems.

Results
We begin with an in silico demonstration of the success of AAD-control method in unpinning, and removing an anchored vortex. To this end, we first produced a stationary rotating vortex, pinned to a circular conductiontype inhomogeneity, in a simulated monolayer of neonatal rat cardiac atrial cells. Next, we created a circular light-induced inhomogeneity, close to the location of the anchored vortex, by applying a light spot to the in silico monolayer, as described in Majumder et al [8]. We observed unpinning of the anchored vortex from the conduction-type inhomogeneity, followed by its re-attachment to the light-induced inhomogeneity. Subsequent movement of the optical spot, away from the location of the obstacle, led to spatial movement of the light-induced inhomogeneity, along with the anchored vortex, as seen in case of free vortices with AAD-control. These findings are illustrated in figure 1.
In our attempts to use the AAD-control method to relocate a pinned vortex, we observed that detachment by the optical spot did not occur successfully, all the time and depended on timing and location of the light-induced inhomogeneity. This inspired us to investigate the factors responsible for detachment of the anchored vortex, and to use the acquired knowledge to improve upon the standard AAD control method, in order to maximize the probability of successful dragging. Therefore, as a first step, we studied the roles of sizes of the two types of inhomogeneities (conduction-type and light-induced) in the process of vortex unpinning and its re-attachment Figure 1. Termination of an anchored vortex by optical dragging (AAD method). (A1)-(A5) Show the representative states of the spiral wave at times 0 ms, 80 ms, 160 ms, 240 ms, and 320 ms, respectively. The vortex was initially anchored to a circular conductiontype (filled white) inhomogeneity of radius r C =0.78 mm, and it was dragged using a circular light spot (translucent blue) of the same size.
to the light-induced inhomogeneity. To this purpose, we initiated a spiral wave in 50 different realizations of the 2D realistic simulation domain (see Methods for details). In each case, we started with the spiral anchored to a conduction-type inhomogeneity of radius r C and placed a circular light-induced inhomogeneity of radius r L at some spatial separation from the conduction-type inhomogeneity, such that the distance between the centers of the two inhomogeneities was d. Then we allowed the system to evolve in time for 4 s and recorded the final state of the spiral. We calculated the relocation probability P relocation for the different combinations of r C , r L , and d. The results of our simulations are presented in figure 2. We observed that at small r C (1.10 mm), P relocation was significantly large, even at d=4 mm, with small r L . The maximum value of P relocation increased with increase in r L , for a given (r C , d), when there was no physical overlap between the inhomogeneities, i.e.
On the contrary, at r C >1.10 mm, detachment and relocation was possible only when there was physical overlap between the inhomogeneities, or the two inhomogeneities were so close that the cells occupying the domain space between them constituted part of the composite inhomogeneity (conduction+light-induced), due to the inexcitability zone created by depolarization gradients of the light-induced inhomogeneity [38]. At r C >1.10 mm and > + d r r C L ( ), relocation occurred very rarely (e.g. 5/50 cases for r C =1.23 mm, r L =1.35 mm, when d=4 mm). Thus, large inhomogeneities apparently exerted very strong pinning forces, which made it extremely difficult (almost impossible at d>4 mm), for a spiral to hop from one basin of attraction to the other.
AAD-control method exploits the use of light-induced inhomogeneities to relocate spiral waves. Thus, in order to ensure successful relocation of pinned waves from conduction-type inhomogeneities of all shapes and sizes, it was important to determine the critical distance at which the probability of successful relocation dropped to 1/e. In order to do so, we conducted a systematic exploration of the possible basins of attraction of circular conduction-type and light-induced inhomogeneities of different sizes, that were located at different distances from, and at different angular orientations relative to the tip of a spiral wave, located at the center of the simulaton domain at time t=0. For this study, we used 3 different sizes of the inhomogeneity (radius 0.78, 1.23, or 1.74 mm), 3 different distances (2, 4, or 6 mm) from the center of the domain, and 12 different orientational 360 in steps of 30°) relative to the instantaneous direction of movement of the tip of the spiral wave. We let the spiral evolve in time and recorded, at the end of 4 s, whether pinning occurred or not. We repeated these studies for 50 different realizations of the simulation domain and calculated the statistical likelihood of pinning (P anchor ) at each value of (r i , θ j ) (i=1, 2, 3, 1j12). The results of these studies are presented in figure 3. Note that pinning to light-induced inhomogeneities occurs faster (400 ms) than to conduction-type inhomogeneities (1000 ms). Our results confirmed that the typical basin of attraction of a circular inhomogeneity, for both conduction-type and light-induced cases, is circular, provided the spiral was exposed to the inhomogeneity for a reasonably long time (>3 s) (see figure 3. At short timescales, while the basin of attraction for a conduction-type inhomogeneity remained circular, with lower average probability of attraction at shorter observation times, the basin of attraction of a light-induced inhomogeneity developped a bias, as in figure 3 of Majumder et al [8] indicating that brief light pulses could be used for directed attraction-based movement, as is required for AAD. At long exposure times (4 s), we determined the radii (R BoA ) of the circular basins of attraction corresponding to different physical sizes of the conduction-type (r C ), and light-induced (r L ) inhomogeneities. If r be the minimum distance from the tip of the spiral wave at time t=0 to the surface of the inhomogeneity, then for small circular conduction-type inhomogeneities (r C 1.0 mm), P anchor was observed to vary linearly with r. At r C 1.0 mm, P anchor was highest (P max ) when the inhomogeneity shared the location of the spiral core. Even then, the absolute value of P max was rather low (see figure). P max increased with increase in r C . Somewhere between r C =0.95 and 1.10 mm, the dependence of P anchor on r changed from linear to inverse-sigmoidal (see figure 3(G)). We call this lengthscale, the critical radius () of the inhomogeneity. With a circular light-induced inhomogeneity, P anchor showed inverse-sigmoidal dependence of r, over the entire range of r L , investigated (i.e. from 0.78 to 1.35 mm)(figure 3(J)). The data obtained from our numerical experiments, was fitted using the following inverse-sigmoidal function: Here, a controlled the offset along the vertical axis. This was varied from 1.5 to 1.8 in order to obtain a perfect fit of the data. Parameters b and c were fixed at 1.8 and 1.9, respectively. These were obtained by optimizing f 1 (r) for the data corresponding to all 5 values of r L . We used the same function to fit the data for different r C . However, for the codunction-type inhomogeneity, the slopes of the inverse-sigmoidal curves were shallower than their light-induced counterparts. Hence b and c were reduced to 1.25 and 1.7, respectively. At r C =1.10 mm, the -P r anchor dependency curve underwent a transition from inverse-sigmoidal (r C >1.10 mm) to linear (r C <1.10 mm). Thus, the parameters in f 1 (r) needed to be adjusted further. We used b=1.15 mm and c=1.5 mm, with a normalization to 0.92, the value of P max , corresponding to r C =1.10 mm. The -P r anchor dependency curves, corresponding to r C =0.78 mm and 0.95 mm, in figure 3(G), were fit using the function ) values of (0.69, 0.44), and (0.52, 0.56), respectively. The radii of the different basins of attraction of our inhomogeneities were given by r BoA ). Figure 3(K) shows the dependence of r BoA on the radius of the inhomogeneity for both conduction-type and light-induced cases. Clearly, the basin of attraction of a light-induced inhomogeneity was almost -1.5 times larger than that of a conduction-type inhomogeneity of the same size (for =  r r r 1.10mm, 1.22 ). We also found that one of the primary causes for relocation failure was the tendency of the spiral to reanchor to the conduction-type inhomogeneity, even if detachment could be forced at initial times, by the presence of the light-induced inhomogeneity. Thus, successful relocation of a spiral wave initial detachment should be followed by rapid movement of the light-induced inhomogeneity that the core of the spiral wave, is brought away from the basin of attraction of the conduction-type inhomogeneity, in order to avoid reanchoring. However, it is not trivial to determine when such moving of the light spot should occur. One way to do it, could be via some feedback from spiral wave top position.
In figure 4, we exploited another possibility of unpinning by introducing a third light-induced inhomogeneity, in the form of a simple linear cut from the surface of the conduction-type inhomogeneity to the surface of the light-induced inhomogeneity. This cut provided the reentrant activity, a pathway to move from the conduction-type inhomogeneity to the light-induced inhomogeneity, by forming a composite inhomogeneity. Subsequent removal of the cut ensured successful relocation of the pinned spiral to the lightinduced inhomogeneity, which could then be dragged away from the conduction-type inhomogeneity, via the AAD-control method. These results are presented in figure 4, alongside the outcome of conventional RDP in the same system, whereby, one applies a train of electrical pulses of the same frequency as the anchored spiral wave, from a location close to the inhomogeneity. We found that AADwc method exhibited highest efficiency in controlling spiral wave dynamics, at the shortest-duration (T light ) light pulse, used for the generation of the lightinduced inhomogeneity (see figure 4(E)). Taking this into consideration, we found that the AADwc method clearly demonstrated a higher success rate in spiral wave relocation, than the conventional RDP or AAD methods (see figure 4(F)).
In the results presented so far, we had used idealized circular scars for the sake of simplicity. Next, we performed simulations with 3 different scar realizations that resemble anatomically accurate patterns of fibrosis [39]. We employed the AAD-control method with linear cut (AADwc) to unpin and relocate the spiral waves in each case. At T light =40 ms, we obtained successful relocation of the spiral from all 3 patterns, namely, compact-, diffuse-and patchy fibrotic inhomogeneities. These results are illustrated in figure 5.
Finally, we performed representative simulations to unpin anchored scroll waves from 3D slabs of neonatal rat atrial cardiac tissue, with a total fiber rotation of 60°across the thickness of the slab. We introduced a cylindrical conduction-type inhomogeneity of base radius r C =1.35 mm and height h=0.94 mm, 1.2 mm, and 1.56 mm, respectively, in a tissue slab of thickness 1.56 mm, which is the approximate thickness of the atrial wall in the neonatal rat heart. With a circular spot of light of r L =1.10 mm, applied to the surface of the slab which shared an end of the conduction-type inhomogeneity, and a linear connector of length 3 mm we detached and dragged a scroll wave to termination(see figure 6).

Discussion
Spatiotemporal control of spiral (scroll) wave dynamics is an important problem in nonlinear sciences and complex systems; in particular, the heart, where their existence is associated with lethal caridiac arrhythmias. Dynamical control can include regulated initiation of new spiral (scroll) waves, suppression of pre-existing waves, anchoring of spiral (scroll waves) and detachment or unpinning of these waves. For decades, researchers have worked at developing control mechanisms for these reentrant waves in the heart. However, only recently, after the advent of cardiac optogenetics, a method to exert real-time spatiotemporal control over spiral wave dynamics by targeting an area close to the core of the vortex, was proposed [8]. Arrhythmias driven by pinned spiral (scroll) waves are extremely important in cardiology. They often occur in patients after myocardial infarction, or as a result of structural remodeling of heart which occur in many forms of heart disease. In this article, we propose a gentle method which can be used to manage such arrhythmias, which is further improvement of the AAD-control method, suggested by Majumder et al [8]. In this method we first unpin and then relocate spiral (scroll) waves that are anchored to conduction-type inhomogeneities. Important that our approach uses the same methodology for both unpinning and relocation processes, thus they can be realized on the same platform.
Although unpinning works in non-modified AAD method, we showed that its efficiency can be improved by producing of temporary simple linear cut. Note, that this is only one possibility and efficiency of unpinning can  potentially be increased by other approaches. For example one can try to use scanning of the tissue, or changing of radius of the illuminated spot or other dynamical protocols etc. This issue needs to be further investigated and the method should be further tuned for each possible particular application of the method.
Our 2D studies predict the existence of a critical size () of the inhomogeneity, below which, the inhomogeneity fails to attract a spiral wave sufficiently in order to pin it. Above , the attraction profile of a conduction-type inhomogeneity resembles that of a light-induced inhomogeneity, which is an established attractor, and therefore makes detachment of spiral waves by conventional methods, rarely possible. Our studies demonstrate that with the suggested improvement, AAD-control method has a close-to-100% likelihood for success in relocating a pinned vortex, which is far better than what can be achieved with other methods such as unmodified AAD or RDP.
In previous studies, a perturbation theory based on properties of response functions, was developed for spiral wave dynamics. These functions are eigen modes of the adjoint operator, linearized on the spiral wave solution [40,41]. Interestingly, for rigidly rotating, as well as meandering, spirals the response functions were numerically found to be localized in the space around the spiral core. This localization area was identified as the key determinant of the region where external perturbation could affect spiral wave dynamics [42][43][44]. It would be interesting to apply response function theory to the problem of anchoring. One could, for example, estimate the distance between the core of the spiral wave and the location of an obstacle, such that, the mere presence of the obstacle can affect the dynamical properties of the spiral wave and possibly predict the likelihood of its anchoring. In addition, it is possible to study the rotation of a spiral wave close to an obstacle and observe more complex dynamics, such as precession around a localized inhomogeneity [45]. However, in order to apply the response function theory to the results obtained in our paper, few questions still need to be addressed. First, one needs to find effective tools for calculation of response functions for complex ionic models such as one used in our study. Second, perturbation theory normally deals with small-amplitude, smooth perturbations, whereas, processes such as unpinning change the topology of the system. These additional challenges need to be addressed in order to reach a full quantitative understanding of the pinning and unpinning processes.
The AADwc method can provide an efficient and optimal strategy in the design of control schemes for the removal of any type of anchored reentry in the heart, as it clearly demonstrates that real-time spatiotemporal control of scroll waves can be achieved, even with sub-transmural illumination of cardiac tissue. Although in practice, the application of optogenetics to cardiology, is, in itself, a debatable topic, recent advances in cardiac optogenetics [46][47][48][49][50][51][52] promise a brighter future. One of the biggest challenges faced by existing scroll wave control methods is that the control pulses (electrical or optical), can only be applied to the surface of the heart. Blue light does not have sufficient penetration depth into human cardac tissue. Currently, the optogenetic tools used to study such control methods respond only to light of short wavelengths in the visible spectrum. Thus, although it would be ideal to have a optogenetic tool, responding to red light, which has deepest penetration into cardiac tissue, such a tool has not been realized so far [53]. This leaves us to consider improvements to our current methods to achieve real-time control over spiral (scroll) wave dynamics, as the only option, aside from drug interventions or ablation (in case of AF). The success of AADwc method in controlling scroll wave dynamics in a 3D slab of cardiac tissue, has opened new doors in cardiac optogenetics, to the design of arrhythmia control schemes by sub-transmural illumination of the heart wall. Detailed investigations into the influence of fiber orientation, increased tissue thicknesses and realistic geometries remain warranted for the future.

Methods
In order to simulate the electrophysiological properties of neonatal rat atrial cardiomyocytes, we followed the model proposed by Majumder et al [54], whereby, the temporal evolution of the single-cell transmembrane potential V is described as follows: I ion is composed of 10 ionic currents, namely, the fast Na + current (I Na ), the L-type Ca 2+ current (I CaL ), the inward rectifier K + current (I K1 ), the transient outward K + current (I to ), the sustained outward K + current (I Ksus ), the background K + (I Kb ), the background Na + (I Nab ) and Ca 2+ (I Cab ) currents, the hyperpolarizationactivated funny current (I f ), and the acetylcholine-mediated K + current (I K,ACh ).
In extended media, the cardiomyocytes were coupled such that the spatiotemporal evolution of V obeyed the following reaction-diffusion-type equation: Here,  is the diffusion tensor. In 2D, the system is isotropic. Thus,  reduced to a scalar =  1 , respectively. For our studies in 2D, we used 50 different realistic in silico monolayers of neonatal rat atrial cardiomyocytes. By realistic, we mean that (a) the simulation domains were circular, with a physical diameter of 1.56 cm (same as in a typical 24-well plate), (b) the cardiomyocytes composing the simulation domains expressed a natural cellular heterogeneity as described in Majumder et al [54], and (c), ≈17% of the cells constituting the simulation domains, were cardiac myofibroblasts [8,54,56]. The cardiac myofibroblasts were modeled as passive circuit elements, according to the MacCannell formulation [57], which were electrically coupled to the cardiomyocytes. These cells were distributed randomly throughout the domain. For our studies in 3D, we used homogeneous slabs of neonatal rat atrial tissue (no cardiac myofibroblasts. We used rectangular parallelepipedal domains of physical size 1.6 cm×1.6 cm×0.156 cm, with fibers rotated through 60°from the base of the slab to the top. We used a model of Chlamydomonas reinhardtii channelrhodopsin-2 mutant H134R, as the optogenetic tool, following the studies of Boyle et al [58], with the exact same parameter values for the opto-channel, as in [8,58]. In general, light has poor penetration inside cardiac tissue. Thus, a computationally efficient and appropriate method to model light attenuation effects, is to solve the photon diffusion equation with a zero source term inside the tissue. In the simplest case, we assume homogeneous absorption and isotropic scattering of light photons. Thus the irradiance can be approximated by an exponential function, as explained in detail in [59]. For our 3D studies, irradiance E e (r) inside the tissue was modeled using the simple distribution function [59] [8], and δ=0.588 mm, the penetration deption for blue light in cardiac tissue [60].
In our studies, we used 2 kinds of inhomogeneity: conduction-type and light-induced. The conduction-type inhomogeneties were modeled as regions of zero diffusion within the simulation domain, with no flux exchange across the boundaries. The light-induced inhomogeneities were produced by affecting the kinetics of the simulated opto-channels in localized regions, thereby forcing most of the opto-channel subunits to move to the open state, and ensuring depolarization of the cell membrane.
For the simulations using realistic inhomogeneities in 2D, we adopted the different fibrosis patterns from [39] to scale, and placed then in the middle of large simulation domains of physical size 3.2 cm×3.2 cm. We initiated spiral waves in these domains by the standard cross-field S1-S2 protocol. Scroll filaments were tracked in 3D, by the method described by Qu et al [55].