Tracking Arrhythmic Spiral Wave Tip Trajectories [Dec. 2015 – Mar. 2016]

This was my first ever research project (December 2015 – March 2016), facilitated via UC San Diego’s Faculty Mentor Program. My mentor was Dr. Wouter-Jan Rappel, who was extremely patient and kind to me. However, I failed him because I quit working on this endeavor after reaching my first milestone. I quit because, after the initial physics-related tackling, this project evolved into a problem of code optimization that I was not ready to commit to at the time. Besides, I did not want to waste Dr. Rappel’s time knowing that I was not mentally committed to this work.

While nothing substantial came out of this project, I am grateful for this opportunity as it was vital in my transition from imagining research to actually performing it. The amount of work required – from literature review, troubleshooting, setting-up, fighting the problem in my head alone, dabbling with theoretical study for the first time, transforming frustration from not-knowing to motivate learning, being lost for weeks, amongst other not-so-convenient-things – caught me off-guard. Regardless, this not-so-successful experience proved to be vital in subsequent research projects in several ways. For instance, I now try my best to: establish clear communication with collaborators; do more than what is expected; and ensure that I am true to myself so that I can commit myself to a problem completely instead of giving post-half-as*ery-guilt even the slightest chance.

Regardless, here is quick description of this project (a.k.a. TL;DR), followed by my original project proposal:

TL;DR: Made a MATLAB function to correct electronic spiral wave tip data from simulations on cardiac dynamics that investigated the effects of targeted ablation at sources of atrial fibrillation. Managed to track wave tip data on a grid with periodic boundary conditions and eliminate data points that violated theoretically imposed conditions.

Tracking Arrhythmic Spiral Wave Tip Trajectories Using Data from In Silico Investigations on the Effects of Targeted Ablation at Sources of Spiral Waves in Terminating Atrial Fibrillation

INTRODUCTION

Atrial Fibrillation (AF) is a common heart complication that affects over 30 million people worldwide. It increases the risk of heart failure, stroke, and other heart-related conditions1,2. Being a type of arrhythmia, AF involves fibrillation in the atria; i.e., the irregular and often-rapid contraction of the heart’s muscle fibers that results from a lack of synchrony between normal sinus rhythm and pulse. Despite its influence on the health of a significant portion of the population, there is but a limited understanding of the mechanisms of fibrillation – and this makes therapy for AF only partly effective. Thus, the study of the initiation, maintenance and control of arrhythmia comprises the work of several researchers who hope to use newfound knowledge of its mechanisms in bettering treatment for AF and other types of arrhythmia.1, 3, 4, 5, 9, 10

Normally, the heart’s muscular contraction is due to a single wave of electrical excitation that signals its cells to contract. During fibrillation, however, the sinus signal (that is generated by the sinoatrial node, and results in the heart’s normal/sinus rhythm) is interfered by other three-dimensional spiral waves of electrical activation that lead to out-of-phase localized contractions8. Regardless, there are two competing general theories: “stable spiral wave induced breakup” (SSWIB), which explains cardiac fibrillation in terms of stable periodic sources (such as motor rotors) that induce wave breakup depending on the heterogeneity of the tissue that it propagates in; and “multi-wavelet reentry” (MWR), which hypothesizes that multiple wavelets can interact with other wavelets, or create new rotors, to sustain AF independent of tissue structure.2, 4, 7, 10 

The first known studies of fibrillation and its mechanisms were performed only at the turn of the 20th century. Mayer showed that circulating electrical activity could be sustained in rings of muscle in jellyfish rings and turtle ventricular muscle indefinitely. Mines and Garrey then formed the concept of anatomic reentry using rings of canine ventricular muscle. Then, as a result of theoretical studies in the 1940s, Wiener and Rosenblueth suggested that reentry around an obstacle was necessary to sustain fibrillation. All these studies revolved around SSWIB. However, Moe then postulated the MWR theory; and this was experimentally proven by Allessie et al in 1985 after they used the results from experiments on rabbit atrial muscle to propose the leading circle hypothesis (according to which electrical waves rotate around an obstacle). However, the notion of rotors generating spiral waves reemerged at around the same time as a result of theoretical studies by Krinsky and Winfree (and was demonstrated experimentally by Davidenko in 1990). With the advent of sophisticated tools and computerized algorithms, much has been learnt about rotors and the mechanisms underlying spiral waves since then.5, 6, 9

Several basic concepts pertaining to cardiac arrhythmias (such as Moe’s multi-wavelet hypothesis) were formulated using simple, generalized computerized models. Models are simpler than reality and therefore allow setups of varying parameters and initial conditions to be studied at different levels of complexity. In order to take into account all properties that influence arrhythmia (such as changes in ionic concentrations, and mechanisms that regulate these changes in the cell membrane), these simplified mathematical models take the form of differential equations that consider only specific elements of the more-complex reality. Mathematical analysis, even with the aid of computers, would otherwise be both time-consuming and very complex.5, 7, 11 However, even though much is learnt as a result of simplification through models, much more is left to be understood about the significance of the factors that are ignored during modeling. For instance, several early proposed mechanisms for spiral wave breakup in studies on the MWR theory used models that considered the domain as homogeneous tissue; when in reality, cardiac muscle tissue has inherent inhomogeneities4, 8. Even though several studies are being carried to investigate the effects of these factors – for instance: investigations on the critical mass hypothesis, which has to do with the tendency of fibrillation to terminate spontaneously in finite-sized tissue3 – much more remains to be done. 

Most of the models used in the study of cardiac dynamics are ionic models that attempt to reproduce cell membrane action potentials5. One such model is the Fenton-Karma model, which “is a simplification of complex ionic models of cardiac membrane that reproduces quantitatively many of the characteristics of heart cells”. The model uses ordinary differential equations to keep track of transmembrane voltage and three ionic currents; and behaves simply enough to be understood analytically.11 As with all models, the Fenton-Karma model has its limitations: it does not allow for derangements in specific ion channels and the influence of drugs on AF to be studied1.

Another reason for the use of models are the ethical and practical limitations to clinical/experimental research on human hearts. Thus, research is limited to: ECGs; body surface mappings; low risk invasive measurements such as catheter mapping; multi-electrode array mapping; and in vitro experiments on explanted human hearts.7 Indeed, a significant portion of the clinical data employed in the Rappel Laboratory’s project on cardiac dynamics is acquired using two 64-pole contact electrode catheters employed in a mapping technique called FIRM (focal impulse and rotor mapping)12.

A recent study1 employed both FIRM guided ablation (in clinical cases) and in silico experiments to show that limited ablation terminated spiral wave re-entry, and that relapse into AF after this treatment is less likely (in contrast to traditional ablation techniques). While these results support the SSWIB theory for cardiac fibrillation, it is difficult to use FIRM to investigate multiple wavelets (due to: the limited spatial resolution resulting from the 64 electrode catheter; specialized software being required to construct maps from these; and because wave tips and their trajectories should be identified from these maps through careful interpretation)2,10. An alternative to creating these spatiotemporal maps is to construct a binary synchronization network using the phase synchrony that is computed between distinct regions in tissue, as was done in another study2. The resulting analysis, when applied to clinical data (from patients via multi-electrode catheters) suggested that the SSWIB theory offers a more-general understanding of arrhythmia, compared to MWR.

The proposed study will use MATLAB to analyze the spatiotemporal data generated by similar in silico studies that investigate the effect of targeted ablation on sources of AF. MATLAB was chosen as the programming tool for this project due to its data analysis capabilities, and more importantly due to its flexibility and accessibility to those who will work on this project in the future. A function that will identify and track the trajectories of arrhythmic spiral wave tips will be created. This function will account for possibly erroneous data points in the data generated by the simulations used in the in silico studies. Being able to determine the stability of spiral wave tips could allow for further analysis on certain mechanisms proposed for the MWR theory (that rely on the type of tip trajectory8). Additionally, the data that will be output by the function (number of arrhythmic spiral tip generations and destructions; their spatiotemporal information; average lifetime; and average Euclidean distance traveled) could be used to possibly contribute to reaction-diffusion equations10 that could allow for the swift prediction of the locations of motor rotors given a new system (a patient, for instance), and ultimately, more-precise targeted ablation.

METHODOLOGY

Data analysis will be performed on a dataset output by in silico experiments on cardiac dynamics. This output will contain spatiotemporal data of spiral wave tips in a periodic two-dimensional domain. In general, the input data will be expected to have the following format:

MATLAB will be used to identify and track the trajectories of spiral wave tips that were generated and destroyed in data pertaining to a given timeframe. The function created will attempt to account for erroneous coordinates that may have resulted from assumptions and approximations made in models used by the computerized simulations.

Two matrices containing x-coordinates and y-coordinates respectively will be created using this data, with each matrix column containing all coordinates pertaining to each consecutive time step. Next, the Euclidean distance between each coordinate pair in a certain time step, (x1, y1), will be calculated with respect to all coordinate pairs in the next time step using the following formula:

            where (x2, y2) is an arbitrary coordinate in the time step that directly follows the time step being considered.

Since, in reality, a spiral wave tip can travel at most a certain maximum Euclidean distance in a given time interval, the calculated Euclidean distances allow for the isolation of coordinate pairs (in the next time step) that are most likely to belong to the trajectory of the spiral tip with the coordinate pair in the time step under consideration, i.e. (x1, y1).

To account for the periodic boundary conditions of the domain, the Euclidean distance for each set of coordinates is calculated for four cases: without shifting any coordinates by the boundary values of the domain; shifting only x-coordinates; shifting only y-coordinates; and shifting both x and y-coordinates. If a coordinate pair that meets the criteria (of having traveled at most the set threshold value for the Euclidean distance) is found, the coordinates corresponding to the two time steps involved are flagged (by setting their values as NaN; “Not a Number”) in both x and y-coordinate matrices. This is to ensure that the same coordinate pair is not considered in the trajectories of two different spiral tips. If more than one match is found, only the case that will lead to the minimum value for the Euclidean distance will be considered.

On the other hand, if no match is found, the function will search for matches in the time step that follows, and the one after that. This is to account for possible errors in data points that may have resulted from approximations/assumptions made by the models used in the simulations that generated the data being analyzed. Finding no match even then implies that the spiral tip has been destroyed, and the spatiotemporal data pertaining to this tip will be recorded in a separate matrix if the tip’s time step is not the dataset’s first time step; because only the spiral tips specifically generated and destroyed within the time period of the data given strongly correspond to spiral waves that are linked to arrhythmia (the others correspond to the heart’s normal electrical activity).

Certain checks will be done to ensure that the final dataset is a valid and accurate representation of the arrhythmic wave tips in the input data. Since spiral wave tips should be generated and destroyed in pairs1, 4, the function will check: that the number of spiral tips is even at all times; and that the spiral tips generated/destroyed are in relatively close proximity to each other. Additionally, the user will be alerted of a possible discrepancy (in the output data) if the function gives a different dataset when fed with input data whose x-coordinates are in reverse order. This is done because the order in which the data is processed might affect the trajectories determined for two spiral tips that have similar trajectories. Since the algorithm chooses the smallest Euclidean distance from all four cases (pertaining to periodic boundaries), and because it allows a range for the Euclidean distances for each trajectory, a coordinate pair that in fact belongs to another tip’s trajectory could be flagged, resulting in the other tip being marked as destroyed (if it will not have another qualifying coordinate pair in its trajectory). This error could propagate to more trajectories, and in such scenarios, the function will ultimately output misrepresentative trajectories for both (and possibly more) spiral tips.

An alternative to this precaution would be to track the trajectories of all spiral tips for all cases that have a match (instead of choosing only the case that led to the minimum Euclidean distance) and choose only the final data set that has the least amount of spiral tip creations and destructions (or vice versa, depending on which case is more representative of arrhythmia in the model used). However, this will come at the cost of processing time and increased function complexity.

The proposed method used will have other limitations too. For one, it will be time-consuming regardless, with a runtime of slightly over five minutes for a dataset with 20,000 time steps and between 1 to 22 coordinate pairs in each time step. This is very likely because the function will use iterations (loops) and calculate Euclidean distances within these loops. Thus, a simple remedy would be to pre-calculate all Euclidean distances and store them in separate matrices for all cases prior to iteration. This is expected to be more efficient because the basic data element in MATLAB is the matrix.

The fact that the threshold permitted for the Euclidean distance is a constant may also be a limitation, as it may depend on a case-by-case basis, and also on the model used in the in silico experiment. All other theoretical limitations of the models used in the simulation studies will apply to this function too.

REFERENCES

  1. Rappel, Wouter-Jan, Junaid A.B. Zaman, and Sanjiv M. Narayan. “Mechanisms for the Termination of Atrial Fibrillation by Localized Ablation.” Circ Arrhythm Electrophysiol (2015). Print.
  2. Vidmar, David, Sanjiv M. Narayan, and Wouter-Jan Rappel. “Phase Synchrony Reveals Organization in Human Atrial Fibrillation.” Articles in PresS. Am J Physiol Heart Circ Physiol (2015). Print.
  3. Zhilin Qu. “Critical Mass Hypothesis Revisited: Role of Dynamical Wave Stability in Spontaneous Termination of Cardiac Fibrillation.” Am J Physiol Heart Circ Physiol (2006): 290:H255–H263. Print.
  4. Karma, Alain. “Physics of Cardiac Arrhythmogenesis.” Annu. Rev. Condens.Matter Phys (2013): 4:313–37. Print.
  5. Fenton, F. H. “Numerical Simulations of Cardiac Dynamics. What Can We Learn from Simple and Complex Models?” IEEE Computers in Cardiology (2000): 27:251-54. Print.
  6. Pandit, Sandeep V., and José Jalife. “Rotors and the Dynamics of Cardiac Fibrillation.” Circulation Research (2013): 849-63. Print.
  7. Tusscher, Kirsten Ten. “Spiral Wave Dynamics and Ventricular Arrhythmias.” Thesis. Utrecht University, 2004. Print.
  8. Fenton, Flavio H., Elizabeth M. Cherry, Harold M. Hastings, and Steven J. Evans. “Multiple Mechanisms of Spiral Wave Breakup in a Model of Cardiac Electrical Activity.” Chaos (2002): 852-92. Print.
  9. Jalife, José. “Déjà Vu in the Theories of Atrial Fibrillation Dynamics.”Cardiovascular Research (2011): 766-75. Print.
  10. Hwang, Minki, Jun-Seop Song, Young-Seon Lee, Changyong Li, Eun Bo Shim, and Hui-Nam Pak. “Electrophysiological Rotor Ablation in In-Silico Modeling of Atrial Fibrillation: Comparisons with Dominant Frequency, Shannon Entropy, and Phase Singularity.” PLOS ONE PLoS ONE 11.2 (2016). Print.
  11. Tolkacheva, E. G., D. G. Schaeffer, D. J. Gauthier, and C. C. Mitchell. “Analysis of the Fenton–Karma Model through an Approximation by a One-dimensional Map.” Chaos Chaos: An Interdisciplinary Journal of Nonlinear Science 12.4 (2002): 1034. Print.
  12. “Rappel Lab Webpage.” Rappel Laboratory, UCSD. Web. 19 Mar. 2016. <http://physics.ucsd.edu/~rappel/>.

Leave a Comment