Pseudospectral Solutions of Reaction-Diffusion Equations that Model Excitable Media: Convergence of Solutions and Applications. by Daniel Olmos Liceaga B.Sc, University of Sonora, 1999 M.Sc, CIMAT, 2002 A THESIS SUBMITTED IN PARTIAL F U L F I L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF Doctor of Philosophy in The Faculty of Graduate Studies (Mathematics) The University Of British Columbia July, 2007 © Daniel Olmos Liceaga 2007 Abstract In this thesis, I develop accurate and efficient pseudospectral methods to solve Fisher's, the Fitzhugh-Nagumo and the Beeler-Renter equations. Based on these methods, I present a study of spiral waves and their interaction with a boundary. The solutions of Fisher's equation are characterized by propagating fronts with a shock-like wave behavior when large values of the reaction rate coefficient is taken. The pseudospectral method employed for its solution is based on Chebyshev-Gauss-Lobatto quadrature points. I compare results for a single domain as well as for a subdivision of the main domain into subintervals. Instabilities that occur in the numerical solution for a single domain, analogous to those found by others, are attributed to round-off errors arising from numerical features of the discrete second derivative matrix operator. However, accurate stable solutions of Fisher's equation are obtained with a multidomain pseudospectral method. A detailed comparison of the present approach with the use of the sine interpolation is also carried out. Also, I present a study of the convergence of different numerical schemes in the solution of the Fitzhugh-Nagumo equations. These equations, have spatial and temporal dynamics in two different scales and the solutions exhibit shock-like waves. The numerical schemes employed are Chebyshev multidomain, Fourier pseudospectral, finite difference methods and in particular a method developed by Barkley. I consider two different models of the local dynamics. I present results for plane wave propagation in one dimension and spiral waves for two dimensions. I use an operator splitting method with the Chebyshev multidomain approach in order to reduce the computational time. I conclude this thesis by presenting a study of the interaction of a meandering spiral wave with a boundary, where the Beeler-Reuter model is considered. The phenomenon of annihilation or reflection of a spiral at the boundaries ii Abstract of the domain is studied, when the trajectory of the tip of a spiral wave is essentially linear. This phenomenon is analyzed in terms of the variable j, which controls the reactivation of the sodium channel in the Beeler-Reuter model. Table of Contents Abstract ii Table of Contents iv List of Tables vii List of Figures ix List of Abbreviations xii Acknowledgements xiii Dedication xiv Statement of C o - A u t h o r s h i p xv 1 1 3 5 8 10 12 12 Introduction 1.1 Preliminaries 1.2 Excitability and the Hodgkin Huxley Equations 1.3 Propagation of an action potential 1.4 Pseudospectral Methods 1.5 Thesis objectives 1.5.1 A pseudospectral solution of the Fisher's Equation. . 1.5.2 A pseudospectral solution for the Fitzhugh Nagumo Equations 1.5.3 Annihilation and reflection of spiral waves at a boundary for the Beeler-Reuter model 15 17 Bibliography 24 2 34 34 37 A Pseudospectral M e t h o d of Solution of Fisher's E q u a t i o n . 2.1 Introduction 2.2 Chebyshev-Lobatto Spectral Approach to Fisher's Equation iv Table of Contents 2.3 2.4 2.2.1 The modified F E ; scaling the p dependence 2.2.2 Pseudospectral solution of the modified F E 2.2.3 Numerical results Chebyshev-Lobatto Multidomain Spectral Method Discrete Singular Convolution; Whittaker's Sine Interpolation 2.4.1 Analysis of the round-off error for D , and D f 2.5 37 40 41 44 47 50 Summary 53 Bibliography 63 3 68 Pseudospectral method of solution 3.1 3.2 3.3 3.4 3.5 Numerical methods 3.1.1 Chebyshev Pseudospectral Multidomain Method. 3.1.2 The Fourier method 3.1.3 Barkley's method 3.1.4 Initial and boundary conditions Numerical studies with kinetic model 1 3.2.1 Plane I D waves 3.2.2 Spiral waves Numerical studies for kinetic model II 3.3.1 Plane waves 3.3.2 Spiral waves Operator splitting method for C M D Summary. . . 72 73 75 76 77 78 78 80 84 84 85 87 89 Bibliography 107 4 112 115 117 A n n i h i l a t i o n a n d r e f l e c t i o n of s p i r a l waves 4.1 Membrane models and numerical methods 4.2 Numerical Results 4.3 Boundary effects on the rotation period in compound rotation 4.4 Annihilation and reflection as a function of the incident angle. 4.5 The role of excitability in controlling reflection and annihilation at the boundary 4.6 A rationalization of the fraction of trajectories annihilated, F (6i), Fig. 4.9 A 4.7 Summary ' 120 122 125 130 132 v Table of Contents Bibliography 149 5 Conclusions 152 5.1 Future directions 155 5.1.1 The role of obstacles and application to cardiac arrhythmias 155 5.1.2 Quantitative factors explaining the annihilation-reflection phenomenon 156 5.1.3 The Bidomain model 157 Bibliography 158 Appendices A Beeler-Reuter equations B Initial conditions used to generate spiral waves • 161 162 vi List of Tables 2.1 2.2 2.3 Analysis of the multidomain method for p = 10 , A t = 10~ , total time of integration T = 0.003; M is the number of subintervals and N is the number of Chebyshev points per subinterval Analysis of the multidomain method for p = 10 , At = 10~ , T — 0.001; M is the number of subintervals and N is the number of Chebyshev points per subinterval Maximum value of the elements in the sum, E q . (2.41), with C equal to IV ) and • When reduced bandwidth (RB) is considered, the index in the sum, E q . (2.41), runs from — W to W, W = 32 Comparison between multidomain method and the Sine ap4 8 5 8 55 55 2 2.4 56 (2) proach with (W=32) for different values of p. T is the total time of integration, M is the number of subintervals and N is the number of Chebyshev points per subinterval 3.1 Convergence of the wavefront speed versus ./V, for kinetic model I with the Fourier method and Dirichlet boundary conditions, e = 0.005, a = 0.3, b = 0.01, x = - 3 0 , x = 30. . . . Convergence of the wavefront speed versus N, for kinetic model I with C M D and no flux boundary conditions, e = 0.005, a = 0.3, b = 0.01, x = - 3 0 , x = 30. N = (N h — 2)N + 2 is the total number of grid points, with N = 40 subdomains and N h is the number of collocation points per subdomain Convergence of the wavefront speed versus A^, for kinetic model I with C M D and no flux boundary conditions, e = 0.005, a = 0.3, b = 0.01, x = - 3 0 , x = 30. N = {N h — %)N + 2 is the total number of grid points, with N subdomains and N h = 13 collocation points per subdomain. L 3.2 L c 91 R S s 3.3 R 56 c L c S 91 R s c 91 vii List of Tables 3.4 Convergence of the wavefront speed for kinetic model I with F D and no flux boundary conditions, e = 0.005, a = 0.3, b = 0.01, x = - 3 0 , x = 30 92 Convergence of the wavefront speed for kinetic model I with the Barkley method and no flux boundary conditions, e = 0.005, a = 0.3, b = 0.01, x = - 3 0 , ^ = 30 92 Convergence of the wavefront speed versus N, for kinetic model II with C M D and no flux boundary conditions, e = 0.001, a = 0.1, b = 0.5, A = 1, x = - 2 0 , x = 20. N = (N h — 2)N + 2 is the total number of grid points, with N subdomains and N h = 13 collocation points per subdomain. 93 Convergence of the wavefront speed versus N, for kinetic model II with F D and no flux boundary conditions, e = 0.001, a = 0.1, A = 1, b = 0.5, x = - 2 0 , XR = 20 93 L 3.5 R L 3.6 L c s 3.7 R S c L 4.1 Measure of the total length Lo in millimeters, and period of rotation To in milliseconds of a unit of trajectory for three cases shown in Figs. 4.6(B-D) Notice that for case D corresponding to F i g . 4.6D, there is no information available for the dashed unit of trajectory, as the trajectory was annihilated at the boundary. 134 List of Figures 1.1 1.2 1.3 1.4 1.5 Phase portrait and solutions of (A) Lotka-Volterra equations (Eq. 1.4) and (B) Perturbed Lotka-Volterra (Eq. 1.5) Two examples of action potential Electrical circuit considered to model the membrane of an excitable cell Propagation of an A P along a nerve axon Corresponding circuit of the membrane of a nerve cell. The R C circuits are given by the local circuit in F i g . (1.3). Current in the x direction (axial current) propagates in a medium with intra and extracellular resistance per unit of length rj and r . It is the transmembrane current e 2.1 2.2 2.3 2.4 Time dependent profiles versus x for p and ./V The relative error versus x in the application of the second derivative matrix operator to the solution at t = 0 for different values of N with (A) p = 10,000 and (B) p = 15,000 The relative error RE versus N in the application of the second derivative matrix operator to the solution at t = 0 at the boundaries, x = XR (upper curve) x = x_ (lower curve) with p equal to (A) 10000 and (B) 15000 (A) Variation of the logarithm of the second derivative versus p, at XR and x_ for the analytic solution at t = 0. (B) \og (D ) versus j Partition of the interval Time dependent profiles U(x, t) versus x for p = 10 with 140 subdomains and 20 points per domain 2 10 2.5 2.6 2.7 2.8 Nj 21 22 22 23 23 57 58 58 59 59 6 Logarithmic value of the i row of (A) D^) (B) D^, i = 1,16,32,48 and 64 Oscillations of U(x,t) after the first time step th 2 60 for 61 61 ix List of Figures 2.9 (A) Instability of U{x,t) at x = x , p = 10 , A t = 1 x 10~ ... (B) Time dependent profiles versus x for p = 10 and times t = 0.001, 0.0015, 0.002, 0.0025 and 0.003 with the 4 L 7 4 (2) (N + 2W + 1) x (N + 1) second derivative operator Dy . . . Local dynamics of the spatially homogeneous case of E q . (3.2) for kinetic models I and II 3.2 Numerical solution for kinetic model I, at a specific time, in a square domain of length 60 x 60 with (A) Dirichlet boundary conditions; (B) Periodic boundary conditions 3.3 Time dependent profiles u(x, t*) versus x for t* = 0,1.66, 3.33, 5, for the propagation of a I D pulse with different numerical schemes for kinetic model 1 3.4 Error, EN(X), (Eq. 3.22) versus x for the I D solution of E q . (3.5) for kinetic model I and a = 0.3, b = 0.001, e = 0.005 at t* = 5 3.5 Convergence of the numerical wavefront speed to the asymptotic speed c as e —> 0, with a = 0.3, b = 0.01, D\ = 1 3.6 Contour plots of the solution u(x, y, t*) of E q . (3.5) for kinetic model I and a = 0.3, b = 0.001, e = 0.005, u = 0.5 at t* = 8. . 3.7 Graphs on the left: I D profiles from cuts of the spirals at y = —15, as given by the different numerical schemes at time t* = 8 for kinetic model I; Graphs on the right: Analysis of EN(X) error (Eq. 3.22) versus x for the profiles shown on the left 3.8 Comparison of contour plots (u(x,y,t*) = 0.5) of the numerical solution of kinetic model I (Eq. (3.5) with (3.3)) for different numerical schemes at t* = 100 3.9 Time dependent profiles u(x, t*) versus x at t* = 0,0.6,1.2,1.8, for the propagation of a I D pulse with different numerical schemes for kinetic model II and a = 0.1, b = 0.5, e = 0.001, A = l 3.10 Error, Efj(x), (Eq. 3.22) versus x for the I D solution of E q . (3.5) for kinetic model II and a = 0.1, b = 0.5, e = 0.001, A = 1 at t* = 1.8 3.11 Convergence of the numerical wavefront speed to the asymptotic speed c as e —> 0 with a = 0.1, b — 0.5, D\ = 1 and A = l 62 3.1 a 94 95 96 97 98 99 100 101 102 103 a 104 x List of Figures 3.12 Convergence analysis for the 2D solution u(x,y,t*) of E q . (3.5) with (3.4) (kinetic model II) with e = 0.001, a = 0.1, A = 1 and b = 0.5 and t* = 5 105 3.13 Comparison of contour plots (u = 0.5) of the numerical solution of (A) kinetic model I for different numerical schemes at t* = 8 106 4.1 4.2 Compound rotation or flower like trajectories, resemble two different curves (A) Hypotrochoid for outward petal flowers (B) Epitrochoid for inward petal flowers Spiral tip trajectories obtained with the B R model with respect to two parameters, gjq and E^ Spiral tip trajectories with the B R model with g^ = 2.38 m S c m and gs = 0.03 m S c m Spiral tip trajectories with the B R model with gs = 0.03 m S c m " and (A,B) g = 2.41 m S c m " ; (C,D) g = 2.34 mScm" Quantities associated with reflections at the boundary. Dashed lines indicate the direction of the trajectory of the tip (A) 6i and 6 are the incident and reflected angle is 95° A ) A unit of a spiral tip trajectory far from the boundary for the RQO case. (B-D) Effects of the boundary on the spiral tip units Spiral interactions with the boundary for an incident angle of 8i = 120°, and different values of t T i p trajectories with 6\ = 70°, and different values of t . . . . Fraction of the trajectories, .FA(0J), that were absorbed by the boundary as a function of the incidence angle 0$ (A) Propagating pulse for the B R model in I D A ) A plane front at a fixed time, with a free end (point c) that propagates from left to right Reflection of a spiral Annihilation of a spiral • • • • Particular examples of trajectories of spirals that (A) annihilate at the boundary. (B) Are reflected at the boundary. . . . a 4.3 - 2 2 Na r 4.6 4.7 g 4.8 4.9 4.10 4.11 4.12 4.13 4.14 B.l 137 2 Na 2 4.5 136 a - 2 4.4 a 135 g Scheme used to generate a tip trajectory with incident angle (A) 6i = 63°, and (B) 0; = 120° 138 139 140 141 142 143 144 145 146 147 148 164 xi List of Abbreviations AP AV BR CMD DSC FD FE FFT FHN HH LV N N ODE OS RD SA SSW RB ch s Action Potential. Atrioventricular. Beeler-Reuter. Chebyshev Multidomain. Discrete Singular Convolution. Finite difference. Fisher's Equation. Fast Fourier Transform. Fitzhugh-Nagumo. Hodgkin-Huxley. Lotka-Volterra. Number of Chebyshev Points. Number of Subdomains. Ordinary Differential Equation. Operator Splitting. Reaction Diffusion. Sinoatrial. Super Speed Wave. Reduced Bandwidth. Acknowledgements I would like to thank G o d for giving me the opportunity to come to Canada to study at such a great University, and also for being my source of wisdom, patience and strength to let me complete this program. I am indebted to Professor Bernie Shizgal, for his patience, time, encouragement and guidance throughout my studies as my research supervisor. I am honoured to have been his student. Special thanks to Eric Cytrynbaum ( U B C ) and Jae-Hun Jung (UMass) who always had the time for answering my questions in mathematical modelling and spectral methods. Special thanks also to Roman Baranowski ( U B C Westgrid) for helping me with all my questions about debugging my codes and guiding me to learn parallel computing. I am very grateful to the math biology people at U B C for letting me be part of their group as well as for all the feedback they gave me on my research. Thanks also to the people in the I A M for the good time in these years. Thanks also to my parents and sisters for having the advice I always needed to hear. Yuliana, thank you for your support and love despite of the distance. Thanks to all of my friends in Vancouver for making these five years in my life, worth living. xiii Dedication To my parents Daniel and Martha. To my sisters Tania, Elizabeth and Graciela. To Yuliana. To duraznito, my peach tree at home. A mi rancho. Nogales, Sonora. xiv Statement of Co-Authorship In 2002, I came to U B C with the idea of pursuing a project related to nonlinear dynamics and synchronization. Professor Shizgal suggested to work on reaction diffusion systems and the modelling of cardiac arrhythmias, a new research area for both of us. We agreed to consider the application of pseudospectral methods to the solution of Fisher's equation, a simple I D reaction diffusion problem. The long term goal was the study of the physics of cardiac arrhythmias. We thus agreed to consider the numerical solution of the Fitzhugh-Nagumo equations, a reduced model of the Hodgkin-Huxley equations; Finally, we decided to work with the Beeler-Reuter model, where I choose to work with the Roo case and the phenomenon of annihilationreflection of a spiral wave at a boundary. M y research started by writting codes initially in Matlab and later in Fortran 90. The codes in this thesis were written and validated by me, except the differentiation matrices, which were written in M A T L A B by J . Weidemman and S. Reddy in their paper " A Matlab Differentiation matrix suite" A C M Transactions on mathematical software, V o l 26 (2000) 465-519. A l l the simulations and results in this thesis are my own work and I am entirely responsible for any errors and inaccuracies. In this thesis, I did write the first draft of each chapter, and subsequent drafts were done together with Professor Shizgal. xv Chapter 1 Introduction The electrical activity of the heart and its role in the activation of cardiac contraction, is an important field of study. A periodic electrical wave of excitation, called an action potential wave, is initiated at the sinoatrial (SA) node, the natural pacemaker of the heart. This wave propagates throughout the atria and arrives at the atrioventricular (AV) node, where after some time delay it propagates to the ventricles via the Purkinje fibers [118]. In normal conditions, this process is repeated approximately 70 to 100 times each minute and is commonly referred as the heartbeat. A n arrhythmia is an abnormal heart rhythm due to atypical generation or propagation of excitation waves, during the process described above. Each year, more than 20,000 Canadians die as a consequence of developing fatal ventricular arrhythmias [96, 122]. One of the proposed mechanisms involved in the development of arrhythmias are spiral waves, a particular form of functional reentry [34, 106]. Spiral waves are self sustained waves of excitation that rotate freely or around an obstacle, reactivating the same area of tissue at a higher frequency than the normal S A node would do, increasing the normal heartbeat rate. In the worst scenario, a spiral wave might break up into smaller spiral waves giving uncoordinated contractions of the heart . When this phenomenon occurs in the ventricles, the heart quivers and loses its ability to pump blood to the body leading to immediate cardiac arrest [34, 118]. 1 Therefore, the development of techniques that help in the understanding of the control and annihilation of spiral waves is an important endeavor. In order to understand spiral wave dynamics, it is necessary to consider excitability of cardiac tissue, which also occurs in nerve cells [57] and to some extent in chemical reactions such as the Belouzov-Zhabotinsky reaction [102]. Computer simulations of models in excitable media have been of great utility to understand and support results obtained experimentally x T h i s condition is referred as fibrillation [34] 1 Chapter 1. Introduction. [24, 76], as well as for proposing new ideas about different regimes in cardiac wave propagation [31, 83, 99, 117]. Computer simulations of excitable media usually involve solving partial differential reaction diffusion (RD) equations, where the reaction term models physiological ion kinetics [16, 67, 68]. A very important task is to find fast and reliable numerical methods to solve these R D equations [25]. Therefore, this thesis has two objectives. The first objective is the development of fast and accurate numerical methods to simulate models for excitable media, with an emphasis on cardiac wave propagation and spiral waves. The numerical schemes considered in this thesis are based on spectral methods, which are known to provide exponential convergence of the solution with respect to the number of collocation points used in the spatial discretization [19, 87]. The second objective corresponds to the understanding of the process of annihilation and reflection of a spiral wave when it interacts with a non excitable boundary; a physical region that blocks the propagation of waves of excitation. As discussed in the previous paragraphs, annihilation of a spiral wave in the context of cardiac wave propagation will terminate the reentry behaviour and inhibit the development of atrial or ventricular fibrillation [34]. In order to get a better understanding of the work in this thesis, we devote the rest of this introductory section to give an overview of the concepts and terminology relevant to spiral waves in excitable media. Therefore, some general concepts about nonlinear dynamics for ordinary differential equations are discussed in Section 1.1. In Section 1.2, we provide a general description of excitability in nerve cells, where also, we present in detail the Hodgkin-Huxley model for the propagation of an action potential in nerve cells. The Hodgkin-Huxley model is the basis for physiological ionic models used for modelling cardiac wave propagation, in particular the one given by Beeler and Reuter [16], that we will consider in Chapter 4. In Section 1.3, we describe the propagation of an action potential along a nerve axon, where also, the property of multiple temporal and spatial scales associated with this phenomenon is discussed. In Section 1.4, we provide a general description of Chevyshev pseudospectral methods which will serve as a basis for developing further methods to solve equations in excitable media. We conclude this introductory chapter by giving an overview of the proposed thesis to be discussed in the following chapters (Section 1.5). 2 Chapter 1. 1.1 Introduction. Preliminaries The well known harmonic oscillator is one of the simplest oscillators. A n example of harmonic motion is a mass attached to a perfect spring. In this case there is no dissipation, giving as a result periodic motion. However, friction or other external perturbations play an important role i n the dynamics of such physical system. Therefore, i n the presence of damping, the periodic motion is lost and the mass tends to equilibrium. The equation for an oscillator with damping is dx dx 2 m—jrz dt + c— at .„ „. + kx = 0 (1.1) where for the perfect spring problem, x is the position, m , c and k are the mass, damping coefficient and spring constant respectively. E q . (1.1) is an example of a linear differential equation. The general form of a linear differential equation is given by f = A» (1.2, where x is a vector and A is a matrix. Its solution is well known and is given i n [62]. However, most physical processes are described by nonlinear P D E s . Simple examples include the carrying capacity for the dynamics of a population or the modelling of an autocatalytic reaction [62]. For these examples, it is not possible to find an equation of the form (1.2) able to simulate such processes. Therefore, an extension of E q . (1.2) is considered and is given by — =Ax + f{2,t) (1.3) where f(x,t) is a nonlinear function, e.g. a polynomial of degree greater than one, a trigonometric function, etc. E q . (1.3), is termed a nonlinear differential equation. For this thesis, we will consider nonlinear partial differential equations. A second example of periodic behaviour is the Lotka-Volterra (LV) model [62]. Lotka [32], developed a theoretical model for chemical oscillations which was also considered by Volterra [32] to model the interaction between two species, a predator and a prey. The L V equations are given by ft = ax-bxy | = -cy + dxy [ • > 3 Chapter 1. Introduction. where x and y represent prey and predator populations, respectively. Solutions to Eqs. (1.4) are given by the pair of time dependent functions (x(t),y(t)) satisfying (1.4) and particular initial conditions x(0) — XQ, y(0) = y . In F i g . 1.1A, solutions (x(t),y(t)) of E q . (1.4) for different initial conditions are shown in the so called phase plane. The blue lines, called trajectories, represent two different solutions for different initial values. The arrows, which are given by the derivative of the solution with respect to time, are tangent to a solution curve passing through a point (x*,y*) in the plane and point in the direction of positive time. In F i g 1.1 A , it is observed that the amount of individuals in each population not only oscillates in time but also do so in a periodic way. In general, the family of solutions of E q . (1.4) are closed orbits that never intersect each other (Fig 1.1A). Each closed orbit has assigned a period and a waveform. However, if the solution is perturbed, the new solution will remain periodic but the period and waveform will change. 0 A small variation of Eq. (1.4) can be considered by adding carrying capacity effects for the prey population. We obtain d J = = ^ A _ h x y = a x _ ^ _ h x y ~cy + dxy If K » 1, Eq.(1.5) is a perturbed system of E q . (1.4). A trajectory solution for E q . (1.5) in the phase plane is shown in F i g 1.1B, where the solution oscillates but it converges to a steady state. Therefore, a small perturbation in E q . (1.4) breaks the periodic behaviour. In real situations, most physical, chemical or biological systems are always subject to perturbations. Therefore, it is desirable to preserve the waveform and period of an oscillatory solution under small perturbations in such systems. This is obtained when a stable limit cycle occurs. A stable limit cycle is a periodic solution such that for any sufficiently small perturbation, the original period and waveform are recovered. In F i g 1.1C the concept of a stable limit cycle is shown. The blue solid curve is the stable limit cycle and the red long dashed trajectories are attracted to the blue one. Stable limit cycles is a property shared by any nonlinear O D E attempting to simulate periodic behaviour. Although the concept of stable limit cycles is very useful in modelling periodic behaviour, there are usually other features in oscillatory phenomena 4 Chapter 1. Introduction. that are important to consider. The amplitude, period and time scale of a process provide important information when trying to describe or predict some phenomenon. A n interesting problem is given when a process occurs on two different time scales. This is the case for some electronic circuits [75, 104], autocatalytic chemical reactions [86, 102] or biological processes [57, 65]. A different class of problems than those having oscillatory behaviour, is given by a phenomenon called excitability. The phenomenon of excitability, which is one of the main themes in this thesis, arises in different systems such as cardiac and nerve cells [57, 103] as well as in some autocatalytic reactions [86], and was introduced in 1952 by A . Hodgkin and A . Huxley [57]. In a series of five papers [57], Hodgkin and Huxley explained the mechanisms by which the membrane potential U in a nerve cell changes when the cell is under activity. Therefore, we present an overview of the model proposed by Hodgkin and Huxley, and discuss the concept of excitability. 1.2 Excitability and the Hodgkin Huxley Equations. The concept of excitability is described by the transmembrane potential U of a neuron axon, when the cell is under activity [57]. A t rest, U has the value of about — 60mV. If a short time pulse of current is applied in such a way that the new potential is —56mV, it will be observed that the value of U will return to its resting initial value U = — 60mV immediately (dashed curve in F i g . 1.2A). However, if the potential is raised until U = — 50mV, the transmembrane potential will go to a large excursion raising its value approximately to 4.0mV and then returning back to the value of —60mV (solid curve in F i g . 1.2A). This phenomenon is called an action potential (AP) and cells with this property are called excitable cells. The value of U above which-an A P is elicited is called the threshold value U h (blue straight line in Fig. 1.2A). t The change of the potential in the cell membrane is due mainly to the passage of ions (Na ,K ) via ion channels through the cell membrane. The ion channels are membrane proteins which allow the passage of specific type of ions. A n essential part of the work by Hodgkin and Huxley was to establish that the Na and K channels can be opened or closed, and that state depends + + 5 Chapter 1. Introduction. on the membrane potential at a given time. The time and voltage dependent conductance of the channels are discussed i n more detail later in this section. The main purpose of the Hodgkin-Huxley (HH) equations is to describe the change of membrane voltage U during an A P . To better understand the dynamics of U during an A P , it is useful to consider the equivalent electrical circuit. Because the membrane area covered with ion channels, which allow the passage of ions, is 100 times smaller than the membrane area that acts as an insulator, the cell membrane can be considered as a leaking capacitor. Due to this property, it is possible to formulate a parallel R C circuit where the insulator membrane can be considered as a capacitor and each type of ionic channel as a variable resistor. The circuit is presented i n F i g 1.3. In the circuit, the resistance satisfies Ohm's Law Is = -gs(U - E ) (1.6) s where S refers to either K, Na or CI ion; Is represents the current of the ion S through the membrane, gg represents the conductance of the membrane for the ion S, U is the membrane potential and Es, referred to as the reversible potential and is given by the Nernst potential. W i t h Kirchoff's law for the circuit i n F i g 1.3, one gets the following equation dU lion (1.7) Cr, dt where U is the membrane potential, C is the membrane capacitance and lion = 9Na(U - U ) + g (U - U ) + 9L(U - U ) with I = g (U - U ) representing a leakage current carried mostly by sodium and chloride ions. However, Hodgkin and Huxley found that the gjy and gx conductances, were voltage as well as time dependent. The use of the so called voltage clamp technique [57], combined with the use of channel blockers, provide an empirical fitting of the dependence of the Na and K conductances on voltage and time. Therefore, g^ and g^ given by m Na K K L L L L a + a r + e a 9K = g~K gNa = g~Na nA m3fl . (I- ) 8 where g , g represent the maximum conductance and the variables m, n, h are voltage and time dependent functions that take values between 0 and 1. The variable n satisfies the equation K Na d?7j — = a ( l - n) - f3 n n n (1.9) 6 Chapter 1. Introduction. where a and f3 are functions of voltage. The variables TO and h satisfy equations similar to E q . (1.9), with parameters a ,a^,,f3 ,PhThe variables m,n, h are known as gating variables; n and m h can be interpreted as the fraction of gates that are open at a time t. Finally, the variables TO, n, h and their respective a's and /3's were calculated empirically i n order to fit the data obtained with the voltage clamp technique [57]. n n m m 4 3 To get a better idea of the role of the gating variables, we focus on the conductances and describe the process of an A P . Initially the Na and K voltage dependent ion channels are closed. A t rest, U « — 60mV as seen in F i g . 1.2A. If U is raised to — 50mV by adding positive ions to the cytoplasm, the Na channels will open in a fast time scale, generating a flux of Na ions to the interior of the cell, so that U increases to almost AOmV (See Fig. 1.2A). During this fast transition in U, two other mechanisms initiate on a slower time scale: the K channels start to open allowing the passage of K outside of the cell, and the Na channels start to close. For this reason, the membrane potential U returns to its original state U = — 60mV. When the Na channels close, they enter into a stage called refractoriness. In a sense, they remain closed as they do not allow the passage of ions, but in refractory state they are prevented from opening again until after a certain period of time, called the refractory period. This refractoriness prevents the cell from being excited before returning to its initial stage. + W i t h the information provided above, the use of the variables m,n and h is clear. The variable n is associated with the opening of K channels, when U changes from —60 to 40mV. In the case of Na channels, they open very fast but then close when there is a change of U from —60 to 40mV\ That is why there are two variables m and h associated with the conductance of Na . The variable TO controls the opening of the channels and h controls the closing of them. The final system of equations provided by Hodgkin and Huxley is given by + Cmf = % -9 nHU-U )-g m h(U-U )-g (U-U )+I 3 K = K Na Na L L app a {l-n)-!3 n=^n n n § = (l - h) - p h = a°°=± ah h (1.10) where I is the applied stimulus current, p^ = °+p and r = \p^ for p = m,h and n. The notation of the gating variables has been changed in order to understand the concept of the multiple time scales. The fast app a p a 7 Chapter 1. Introduction. opening of Na channels, and the slow opening and closing of K and Na channels, respectively, are quantified by the variables r , r „ and Th in E q . (1.10), with the relations r « Th and r << T , giving a multiple time scales problem. If the process of opening the Na and K channels is on the same time scale, it is not possible to generate an A P . m m m n The H H equations have been used as a basis of different models for propagation of potentials, mostly in heart cells [16, 23, 67]. Some examples of such models are the Beeler-Reuter (1977) and the Luo-Rudy (1991,1994) equations for ventricular cells [16, 67, 68], the Noble model (1984) for the SA node [78] and the Courtemanche et al (1998) and Nygren et al (1998) models for atrial cells [23, 79]. The multiple time scale problem for an A P is shown in Figs. 1.2A and 1.2B, for typical A P for a nerve cell and myocardial cell, respectively. In the myocardial A P , a plateau at U « 0 is observed and is due to the balance between the influx and efflux of Ca and K ions, respectively, to the cardiac cell. The Ca ions are involved in the contraction mechanisms of cardiac cells. 1.3 Propagation of an action potential. In the previous section, we described the behaviour of U for the spatially homogeneous case. However, an A P travels spatially along the neural or the heart cells. The propagation of an A P can be explained from F i g . 1.4. In the Figure, the cylindrical structure at the lower part represents an axon in a nerve cell. For simplicity we assume propagation of the A P in the one dimensional longitudinal direction x, where propagation goes from left to right as shown by the arrow. We focus on the point in space along the cell where an excitation occurs as shown in F i g . 1.4. As discussed in the previous section, an A P is elicited when the membrane potential reaches a value above the threshold of excitation. A t the beginning of an excitation, sodium ions get into the cell very fast giving the abrupt change in the A P in a phase called depolarization. The sodium ions that enter the cell during the depolarization, move in the interior of the cell following its concentration gradient to parts where the cell is ready to accept an A P . In these neighboring parts the movement of sodium ions cause the transmembrane potential to increase until it reaches the threshold of excitation, giving rise to a new A P . This process is repeated and an A P propagates along a cell. 8 Chapter 1. Introduction. The process of a propagating A P can be associated with the circuit shown in F i g . 1.5. In this case, local R C circuits such as the one shown in F i g . 1.3 are interconnected in parallel, each one representing a patch in the cell membrane. The flow of positive ions from left to right i n F i g . 1.5 define the axial current which has two components, an intracellular, It, and an extracellular, I . These currents are proportional to the changes in voltage with respect to x, and this relationship is given by Ohm's law as Ii = 7 ^ 7 ^ and I = 7^7^?, where r-j and r are the intra and extracellular resistances per unit of length of the medium. The changes of the intra and extracellular currents at a point x along the cell are due to the transmembrane current per unit of length I , which means that I = Qj^ = §£ . B y taking U = Ui — U , the transmembrane potential and, IT = I% + I to be the total axial current, we obtain that e e e = t L t e e dx \ri + r e dx where IT was considered constant. Because I is given as the sum of the ionic and capacitive currents, we get t where p is the perimeter of and ionic current per unit respectively. B y taking the membrane potential U equation the axon, and C and Ii are the capacitance of area of membrane, as given in E q . (1.7), and r independent of space, we observe that satisfies the one dimensional reaction diffusion m on e *L- #!!L_± ; D I dt - dx* f l l 3 ) cj v t m { 1 A 6 ) where D = „ / ,—r is the diffusion coefficient. The typical value for D for a myocardial cell is about 0 . 1 m m / c m [34]. Equation (1.13) is the simplest P D E considered for modelling the propagation of an A P in nerve tissue. 2 2 The profile of a propagating A P in space is similar to the profiles shown in F i g . 1.2 where U is plotted versus time. Therefore, U also changes i n two different spatial scales, where this phenomenon is emphasized in F i g . 1.2B for the A P in a myocardial cell. To solve numerically equations of the reaction diffusion type for nerve or cardiac wave propagation (Eq. 1.13), is a very difficult task due to the different spatiotemporal scales as discussed 9 Chapter 1. Introduction. by Cherry et al [25]. The main problem arises because of the fast transitions in U versus time as shown in F i g . 1.2B. In order to solve accurately the dynamics of U during this transition, it is necessary to consider a very small time step which might not be needed for other regions of the A P , as shown in F i g . 1.2B for the time interval t e [50, 300]. In the same way, the propagating pulses of excitation, which have a fast changing transition at the front, require refinement in the mesh in order to obtain an accurate solution. In most of the cases, simulations of cardiac tissue need to be run for a total time two or three orders of magnitude larger than the time at which the fast transitions in time occur. Solving accurately the fast time process implies a small time step, increasing the computation time of the overall process in an unnecessary way. This issue, together with the stiffness of Eqs. (1.10), presents a very difficult computational problem to solve. 1.4 Pseudospectral Methods. In order to simulate the dynamics given by E q . (1.13), we develop numerical schemes based on spectral methods. In this section, we present one of the most important spectral methods in the literature, the Chebyshev pseudospectral method. This method is based on the expansion of a solution in Chebyshev polynomials. Consider the Chebyshev-Gauss polynomials, Tk(x), which are orthogonal with respect the weight function w(x) = (1 — x )~ l on the interval [—1,1], gives, 2 - f w{x)T (x)Ti{x)dx = k l 2 \*6 ktl (1.14) where c = 1 for all k except for cn = 2. The Lobatto quadrature points and weights associated with the Chebyshev-Gauss polynomials [19] are given by Xi = — cos ( ^ ) and the weights are Wi = jj for all i except WQ = WN = 577 [18, 19, 87]. These points and weights provide the approximate quadrature, k 1 / N w(x)f(x)dx~Y^ if( J w x ( ) L15 where N is the number of points. Since any piecewise continuous function, / G L^fO, 1] can be expanded in a Chebyshev polynomial series that is 10 Chapter 1. convergent in the mean of the Introduction. norm, we have N /(z)«/AKx)=^a f c T f c (x) (1.16) k=0 where 2 a* = J w{x)f{x)T {x)dx l k (1.17) CfcTT W i t h Eqs. (1.15) to (1.17) we obtain the interpolation algorithm JV where the interpolating polynomials, Ij(x), are given by 2 N Ij(x) = ^-^2 T ( )T (x) Vk k Xj k (1.19) k=0 where = v?j = 1/2 and ^ = 1 if k ^ 0, N. The nt/i derivative of f(x) at the quadrature points is then given approximately by f P i x ^ ^ l f X x , ) ! ^ ) (1.20) j=0 If we denote by f, the vector of the function evaluated at the ChebyshevLobatto points, E q . (1.20) can be rewritten as f (n) = D (n) . ( f L 2 1 ) and thus, the second derivative matrix is then identified with Bj? - (1.22) This is the basis for the Chebyshev pseudospectral method. The application of the formalism discussed in this section to the solution of E q (1.13) in one dimension reduces the problem to the solution of a set of N + 1 coupled nonlinear system of O D E s , given by d ^ OX = >mW-J- l AT (1.23) I n O m 11 Chapter 1. Introduction. where, U and I\ are U and I given in E q . (1.13), respectively, evaluated at the collocation point Xi. In this case, A = — ® includes the diffusion coefficient D from E q . (1.13), and a constant that arises from the transformation of [X^^XR] to [—1,1], where [XL^XR] defines the range of the position inside the cell in one dimension. l on ion 4 1.5 Thesis objectives. The objectives discussed at the beginning of this chapter, are explained in this section in greater detail, where also more concrete questions are addressed. We start by discussing the first objective, which consists in the development of numerical schemes for P D E s of the form of E q . (1.13), based on spectral methods (Sections 1.5.1 and 1.5.2). In Section 1.5.3, a numerical study of the annihilation and reflection of a spiral wave at a boundary in excitable media is discussed. In Section 1.3, we show that the solutions of E q . (1.13) are propagating pulses that change over multiple spatial scales, simulating the process of a propagating shock. In order to solve E q . (1.13) numerically, we consider a spectral method approach to be discussed in Chapters 2 and 3. Spectral methods have long been known to provide very accurate and rapidly convergent solutions of partial differential equations with smooth solutions [18, 19, 87]. These methods generally provide an exponential convergence of the solution versus the number of collocation points. In recent years, spectral methods have also been used for the solution of differential equations with solutions that resemble shock waves or fronts typical of hyperbolic partial differential equations [41]. There is an ongoing interest to further adapt spectral methods to differential equations like E q . (1.13), which have rapidly varying and propagating solutions. Therefore, in this thesis we test the performance of spectral methods in two specific examples, named the Fisher's equation and the Fitzhug-Nagumo equations, which we describe now, and discuss at length in Chapters 2 and 3, respectively. 1.5.1 A pseudospectral solution of the Fisher's Equation. Fisher's equation ( F E ) , which was originally proposed by Fisher [35] as a model for the spatial and temporal propagation of a virile gene in an infi12 Chapter 1. Introduction. nite medium, represents one of the easiest models of an equation of the R D type, and is given by E q . (2.1). The solutions for this equation over an infinite domain with initial and boundary conditions given by E q . (2.2) are propagating fronts as shown in F i g . 2.1. A closely related problem of E q . (2.1) is the modified F E given in E q . (2.5). In this case, the propagating fronts are very steep having a shock-like wave behaviour as shown in Fig. 2.6. Although F E is not related to the phenomenon of excitable media, it is an excellent example to test the performance of pseudospectral methods for problems that present propagation of fronts in two different spatial scales. This same feature is exhibited by the solutions of equations for an excitable cell (Eq. 1.13) as discussed in the previous section. This common feature makes F E a good starting point in the effort of developing a numerical scheme for problems with propagating shock-like wave solutions, like those arising in excitable media. Different numerical schemes have been considered to study F E . One of the first numerical solutions of Fisher's equation was presented by Gazdag and Canosa [38] with a pseudo-spectral approach. Implicit and explicit finite difference algorithms have been reported by different authors such as Parekh and Puri [84] and Twizell et al [100]. The works of Mickens [71, 72] considered time stepping aspects for finite difference algorithms. The work by Hagstrom and Keller [48], where the main goal was to develop asymptotic boundary conditions, considered a centered finite difference algorithm. Rizwan-Uddin [89], compared a nodal method with non standard finite differences scheme. A Galerkin finite element method was used by Tang and Weber [97] whereas Carey and Shen [20] employed a least-squares finite element method. A collocation approach based on Whittaker's sine interpolation function [18, 94] was also considered by Al-Khaled [1], and Zhao and Wei [120]. The work by Gourley [44] considered a nonlocal form of F E . The study by Zou [123] was concerned with another modified form of F E including time delay and the work by Roessler and Hussner [90], considered finite elements for a two dimensional F E . More recent studies have considered comparisons of the numerical and exact solutions of E q . (2.5), where the exact solution given by E q . (2.6). Solutions of E q . (2.5) with large p have been referred to by Zhao and Wei [120] as super speed wave (SSW) types. W i t h an increase in p, the propagating front steepens and this presents a challenging numerical problem to both resolve and track the front. This rescaled version of F E was considered by L i et 13 Chapter 1. Introduction. al [66] in their study of moving mesh strategies in finite difference methods of solution of partial differential equations. They commented that moving mesh methods are not recommended for such reaction diffusion problems for which the diffusion term is much smaller than the reactive term. Subsequently, Q i u and Sloan [88] carried out a detailed comparison of different moving mesh strategies [49, 50] and concluded that these methods are not easily adapted for equations analogous to F E with steep fronts. The authors have also applied these methods to Burger's equation [14, 15, 74] with similar shock like solutions. In the course of these studies on F E , numerous workers have also reported a sensitivity of the numerical solutions to perturbations owing to numerical noise in the solutions or round-off errors resulting in instabilities. Therefore, one of the objectives of this thesis to be discussed in length i n Chapter 2, is to apply spectral methods to the solution of F E . A second objective is to establish a relationship between the problems given by E q . (2.1) and E q . (2.5), which we demonstrate lead to equivalent numerical problems. In the first instance, we consider a single fixed domain suitably truncated to the interval [XL,XR] and apply a spectral method based on Chebyshev-Lobatto points. This approach gives good results but only for relatively small values of p. The difficulty encountered for larger values of p is traced to round off errors in the application of the second derivative spectral matrix operator to the solution. The occurrence of round-off errors for these matrix derivative operators have been well documented [5, 6, 7]. We examine in detail the origin of these round-off errors i n the application of spectral methods to F E . A multidomain approach developed by Shizgal and coworkers [69, 115, 116] for the solution of Burgers equation and the Navier Stokes equation was then employed for the present thesis. We obtain accurate stable solutions of F E for relatively large values of p, with the appropriate division of the domain [XL,XR] into subdomains. Also in Chapter 2, we present a comparison of the multidomain approach with the one based on Whittaker's sine interpolation [18, 21, 109] employed by Zhao and Wei [120] and by Wei [108] that they refer to as the discrete singular convolution (DSC). 14 Chapter 1. Introduction. 1.5.2 A pseudospectral solution for the Fitzhugh Nagumo Equations. A second example that we consider for developing fast and accurate algorithms for solving R D equations with excitable dynamics is given by the Fitzhugh-Nagumo ( F H N ) equations (Eq. 3.2). The F H N equations, which were developed by Fitzhugh [36] and Nagumo et al [75], are a reduction of the set of four ordinary differential equations given by the H H equations (Eq. 1.10) to a two dimensional system of O D E s . The variables u and v in E q . (3.2), represent the fast (U, m) and slow (n, h) changing variables in E q . (1.10). However, the F H N equations have been applied to numerous other problems for over four decades [8, 11, 54, 56, 57, 81, 102]. For example, the F H N equations are employed to describe the C O oxidation on Pt(110) [8], the oscillation of calcium concentration in cells [57] and the study of reentry in heart tissue [34, 61]. The algorithm developed for the F H N equations, to be discussed at length in Chapter 3, can be applied to study the propagation of waves in excitable media including spiral waves [10, 54, 55, 59, 83, 111]; to analyze spiral tip trajectories and their transition to complex patterns [10, 54, 111]; the same algorithm can be used to study spiral breakup [46, 82], drift of a spiral tip due to the interaction between spiral waves or via wave trains [43, 113], and competition between spirals and periodic circular pacemakers [63]. In this thesis we consider two examples of the F H N equations (Eq. 3.2) with different kinetics, given by Eqs. (3.3) and (3.4), which we refer as kinetic model I and II, respectively. Kinetic model I was proposed by Barkley [11] whereas kinetic model II is the classic cubic F H N local dynamics [57]. In both models, the variable u(t) has fast transitions as shown in F i g . 3.1. These fast transitions are observed to occur also in space with the solutions of the R D equation (Eq. 3.2) with the reactive term given by kinetic models I and II, giving for each kinetic model, a two temporal and spatial scales problem. For kinetic model I, Barkley [11] developed a numerical scheme based on a finite difference method, which made use of the particular feature that u(t), the local behaviour of u, is very small during the segment d-e shown in F i g . 3.IB. In his numerical scheme, Barkley sets u(x,t) to zero if it is less than some particular small value and thus, avoids the computation of u(x, t) during this time interval. As a result the computational time is reduced considerably. 15 Chapter 1. Introduction. In order to solve the spatial dynamics of F H N equations in two dimensions, different methods have been considered. K a r m a [54] and Mitkov et al [73] used a Fourier pseudospectral method; X i e [113] and Amdjadi [2] used an explicit finite difference scheme; Gottwald [43] and Diks [28] used a method based on finite differences reported by Barkley [11]; Jones and O'Brien [53] solved the Gray-Scott equations with a Fourier pseudospectral method. The F H N and the Gray-Scott equations belong to the family of excitable systems [86] and exhibit the problem of multiple time and spatial scales. However, despite all the numerical experiments realized with the F H N equations to date, it remains important to develop accurate numerical methods for the solution of E q . (3.2), which occurs in two different spatial and temporal scales. In order to obtain an accurate numerical solution for the F H N equations (Eq. 3.2), it is necessary to increase the number of sample points in the spatial grid in regions where fast transitions take place. Therefore, the problem of obtaining a reliable solution when solving systems such as the F H N equations involving different time and space scales represents a challenging endeavor. Therefore, the study of the F H N equations has two main objectives. The first is to continue the ongoing effort to develop pseudospectral methods for an accurate and efficient solution of reaction diffusion equations. This study follows on the work dealing with Fisher's equation in Chapter 2, for which the solutions develop propagating fronts with short spatial variation. The second objective is to compare the solutions obtained with pseudospectral methods with those obtained by Barkley [11] for E q . (3.2) with (3.3), and the finite difference method, in accuracy and computational time. There have been many studies of reaction diffusion systems that are based on k i netic model I proposed by Barkley [3, 8, 43, 64, 70, 73, 121]. Pseudospectral methods are generally considered useful for solving smooth problems and to provide exponential convergence of the solution with respect to the number of collocation points used [19, 87]. However, it has been demonstrated recently that pseudospectral methods can provide a significant improvement over finite difference methods for non smooth problems that develop shocks and steep fronts [42, 80], features shared by the solutions of the F H N equations. Pseudospectral methods provide solutions to partial differential equations defined on a grid of collocation points which are the quadrature points for a given polynomial basis set [18, 19, 37, 40, 87]. The main advantage of pseudospectral (and spectral) methods is the "ex16 Chapter 1. Introduction. ponential" convergence of the solutions versus the number of collocation points. There have been only a few studies of the applicability of pseudospectral methods to the numerical solution of reaction diffusion equations [9, 30, 53, 55]. Another aspect of reactive diffusion equations of which the F H N equations are a subset, is that the time scale for the u variable is short compared to the time scale for the v variable. The discretized equations that result from a pseudospectral method are stiff and generally a very small time step is required for their integration. In Chapter 3, we provide the details of the various numerical methods for the solution of the F H N equations that are compared. The first method is the Chebyshev-multidomain, based on Chebyshev polynomials, which is used also for the the study of fronts in one dimension for Fisher's equation. Also, the Fourier pseudospectral (Fourier) method is considered. Then, we compare the solutions for kinetic model I obtained with pseudospectral methods with those with finite difference techniques and the method presented by Barkley [11]. We validate their efficiency in the solution of the F H N equations and the description of spiral waves. We also generalize the Fourier method to the non periodic case which has not been considered previously. Analogous results for kinetic model II are also presented. The use of an operator splitting method for the time integration to decrease the computational time for the Chebyshev-multidomain is discussed. 1.5.3 Annihilation and reflection of spiral waves at a boundary for the Beeler-Reuter model. The numerical methods developed in Chapters 2 and 3 for the F E and F H N equations are considered to solve the Beeler-Reuter (BR) model, to study the phenomenon of reflection and annihilation of spiral waves at a boundary. The B R model, developed by G . W . Beeler and H . Reuter i n 1977, was the first model based on the Hodgkin-Huxley formalism (Section 1.2) to study excitability in ventricular cells [16, 57]. The equations of the B R model are provided in the Appendix of Chapter 4. At the beginning of this chapter, we discuss the importance of spiral waves in cardiology. Spiral waves are thought to be responsible for the development of certain arrhythmias [34, 60], which may lead to sudden cardiac death [34]. Therefore, an understanding of their control and annihilation is a very important task [34, 52]. 17 Chapter 1. Introduction. A very important feature of spiral waves is the behaviour of their tip, as discussed by Fenton et al [34] and Comtois et al [22]. For a particular set of parameters characterizing the medium, the spiral tip executes different trajectories which can be circular or more complicated patterns [12, 54, 111]. This phenomena is referred to as spiral meandering and has been studied with different kinetic models [13, 31, 34, 111]. A second feature of the spiral tip, called spiral wave drift, is the response of the spiral wave to an external perturbation. Some examples of external perturbations include interaction between two spirals [114], interaction of spirals with a boundary [39, 117] and the drift induced by light as in the Belousov-Zhabotinsky reaction [47]. In Chapter 4, we study numerically a spiral wave in the meandering regime, and the drift effects of the boundary of the domain on the tip trajectory of the spiral. In the simulations, we observe two different phenomena due to the interaction of the spiral wave with a boundary, which are reflection and annihilation of the spiral wave. Experiments in isolated cardiac tissue [27, 51, 52, 85] have been carried out to study the interaction of spiral waves with obstacles and the boundaries of the tissue. Davidenko et. al. [27], Pertsov et. al. [85], and Ikeda et. al. [51] showed annihilation of spiral waves at the boundary. Also, Ikeda et. al. [52] observed attachment o f meandering spirals to an obstacle of some minimum size. Therefore, the phenomena of annihilation and reflection observed in the computer simulations in Chapter 4 suggest that the interaction of spiral-boundary and spiral-obstacle do not necessarily end in annihilation at the boundary or attachment to the obstacle as observed experimentally. This is a very important issue in cardiology, and a proper understanding of the conditions that lead to the annihilation of the spiral at some inexcitable obstacle can help in the development of techniques on how to eliminate arrhythmias generated by the spiral behaviour [31, 51]. Spiral drift due to boundary effects has been considered previously. It has been found experimentally [39] and studied numerically with no-flux boundary conditions [117]. Gomez-Gesteira et al. [39] considered the BelousovZhabotinsky reaction and found that the boundary affected the trajectory of the spiral tip. The trajectory moved along the boundary, whereas in other cases the spiral was annihilated at the boundary [39]. Yermakova and Pertsov [117] analyzed the effects of the boundary on the trajectory of the spiral tip that followed a circular path. B y considering no flux boundary conditions, they showed that the period of the spiral increases when the core 18 Chapter 1. Introduction. of the spiral is close to the boundary. Also, they showed that the center of the circular trajectory, drifted at constant speed along the boundary giving as a result a trajectory resembling the shape of a trochoid. However, the case when the trajectory of the spiral tip meanders and traces a more complex pattern other than a circle, has not been considered. The phenomena of meandering, also referred to as compound rotation, was first noted by Winfree [110]. Figure 4.1 illustrates the two types of trajectories of the tip of a spiral wave constructed by considering the motion of the tip of the arrow attached to the small circle of radius h that rotates either on the inside or outside of the large circle of radius R. When meandering occurs, the trajectory of the tip executes a flower like pattern [12, 111], where the petals lie on a circle of radius R. When the petals lie outside the circle, the trajectory resembles a curve called a hypotrochoid (Fig. 4.1A), whereas in the case the petals lie inside the circle, the trajectory resembles an epitrochoid (Fig. 4.IB). B y considering different parameter values in a particular model with excitable kinetics [12, 31], it is possible to take the limit R —> oo. In this case, which we refer to as the limiting Roc case with R > > 1, the flower has almost an infinite radius and the petals lie essentially on a straight line. When the spiral meanders and is close to the boundary (with no flux boundary conditions), it is observed that the trajectory is annihilated or reflected by the boundary. In the case where the trajectory is reflected, the angle of reflection is not necessarily equal to the angle of incidence. Also, the reflection angle is very sensitive to the position along the spiral tip trajectory, at which the trajectory hits the boundary. Therefore, the main question we address in Chapter 4 is to find under which conditions the spiral is annihilated at the boundary. In order to analyze the effects of the boundary on a meandering spiral, we focus the attention on the degenerate limiting i?oo case. The infinite radius regime is just a transition from the outward petal to the inward petal flower tip trajectory and therefore is not a generic behaviour [12]. The analysis of the R^, case is considered due to its simplicity compared to the case of finite R. Near the boundary, the behaviour of the tip of a spiral can be approximated by the Roo case, and therefore, the results obtained for this limit may also provide an understanding for the case when R is finite. Therefore, in Chapter 4 we provide the equations of the model used for the simulations followed by a description of the numerical methods employed 19 Chapter 1. Introduction. in their solution. We then present the results of numerical simulations for different values of R, including the limiting case. Then, we concentrate on the i?oo case where we analyze some of the properties of annihilation. A n argument based on the reactivation variable j for the sodium channels, is discussed to explain the phenomenon of annihilation and reflection of a spiral at a boundary. It is important to mention that for the present thesis, we only provide a qualitative description of the phenomenon of annihilation and reflection at a boundary. The interaction of the trajectories in the R^ case with a boundary has not been considered before. In Chapter 5, we summarize the results obtained i n the thesis and provide some potential applications of the research findings as well as the future directions of the research area. 20 Chapter 1. Introduction. 11111111111111111111111111111 0.5 1 1.5 2 2.5 (C)~ Figure 1.1: Phase portrait and solutions of (A) Lotka-Volterra equations (Eq. 1.4) and (B) Perturbed Lotka-Volterra (Eq. 1.5). The arrows are tangent to the solution curves. Note that in (B) the periodic behaviour is broken and the solution converge to a stable steady state.(C) Stable limit cycle (solid blue). Close solutions to the limit cycle (long dashed red), are attracted to the stable limit cycle. 21 Chapter 1. Introduction. 400 Figure 1.2: Two examples of action potential. Membrane potential U versus time (A) In a nerve axon. A n A P (solid curve) is elicited when U is above. threshold U h (blue straight line). No A P is elicited if U is below U h (dashed curve). (B) In a myocardial cell. Note the fast and slow changes in U in time. t t Q- g 7 Na + + + c u r - ir ^ u oFigure 1.3: Electrical circuit considered to model the membrane of an excitable cell. The capacitor represents the cell membrane and the resistors the ion channels. gNa,9K d 9L represent the conductance, where g = ^, i.e. the inverse of the resistance. a n 22 Chapter 1. Introduction. Direction of propagation Figure 1.4: Propagation of an A P along a nerve axon. The propagation is considered in one dimension. U h is the threshold of excitation. t Outside cell r AW W\r J RC e RC RC -n— —vw1 r Inside cell Figure 1.5: Corresponding circuit R C circuits are given by the local direction (axial current) propagates resistance per unit of length r-j and ; of the membrane of a nerve cell. The circuit in F i g . (1.3). Current in the x in a medium with intra and extracellular r . It is the transmembrane current. e 23 Bibliography [1] K . Al-Khaled, Numerical study of Fisher's reaction-diffusion equation by the sine collocation method, J . Comput. Appl. Math. 137 (2001) 245-255. [2] F . Amdjadi, J . Gomatam, Spiral waves on static and moving spherical domains, J . Comp. Applied Math., 182 (2005) 472-486. [3] I. Aranson, H . Levine, L.Tsimiring, Controlling spatiotemporal chaos, Phys. Rev. Lett. 72 (1994) 2561-2566. [4] E . M . Azene, N . A . Trayanova, E . Warman, Wavefront-Obstacle interactions in cardiac tissue: A computational study, A n n . Biomed. Eng., 29 (2001) 35-46. [5] R. Baltensperger, J . P. Berrut, The errors in calculating the pseudospectral differentiation matrices for Chebyshev-Gauss-Lobatto points, Comp. and Math. Appl. 37 (1999) 41-48. [6] R. Baltensperger, Improving the accuracy of the matrix differentiation method for arbitrary collocation points, Applied N u m . Math. 33 (2000) 143-149. [7] R. Baltensperger, M . R . Trummer, Spectral differencing with a twist, S I A M J . Sci. Comput. 24 (2003) 1465-1487. [8] M . Bar, N . Gottschalk, M . Eiswirth and G . E r t l , Spiral waves in a surface reaction: model calculations, J . Chem. Phys. 100 (1994) 12021214.' [9] E . Barillot, J . Boissonade, Asymptotic pseudospectral method for reaction-diffusion systems, J . Phys. Chem. 97 (1993) 1566-1570. [10] D . Barkley, M . Kness and L . Tuckerman, Spiral-wave dynamics in a simple model of excitable media: The transition from simple to compound rotation, Phys. Rev. A 42 (1990) 2489-2493. 24 Bibliography [11] D . Barkley, A model for fast computer simulation of excitable media, Physica D , 49 (1991) 61-70. [12] D . Barkley, Spiral meandering, Chemical and wave patterns, (Kluwer academic publishers, Netherelands, 1995). [13] J . Beaumont, N . Davidenko, J . M . Davidenko, J . Jalife, Spiral waves in two-dimensional models of ventricular muscle: formation of a stationary core, Biophys. J . 75 (1998) 1-14. [14] G . Beckett, J . A . Mackenzie, A . Ramage, D . M . Sloan, O n the numerical solution of one-dimensional P D E s using adaptive methods based on equidistribution, J . Computat. Phys. 167 (2001) 372-292. [15] G . Beckett, J . A . Mackenzie, A . Ramage, D . M . Sloan, Computationalsolution of two-dimensional unsteady P D E s using moving mesh methods, J . Computat. Phys. 182 (2002) 478-495. [16] G . W . Beeler and H . Reuter, Reconstruction of the action potential of ventricular myocardial fibres, J . Physiol. 268 (1977) 177-210. [17] A . Beuter, L . Glass, M . C . Mackey, M . S . Titcombe, Nonlinear dynamics in Physiology and medicine, Int. Appl. M a t h . [18] J . P. Boyd, Chebyshev and Fourier Spectral Methods (Dover, New York, 2000). [19] C . Canuto, M . Y . Hussaini, A . Quarteroni and T . A . Zang, Spectral methods in fluid dynamics, Springer Series in Computational Physics, (Springer, New York, 1988). [20] G . F . Carey, Y . Shen, Least-squares finite element approximation of Fisher's reaction-diffusion equation, Numer. Methods Partial Differential Equations, 11 (1995) 175-186. [21] D . T . Colbert, W . H . Miller, A novel discrete variable representation for quantum mechanical reactive scattering v i a the S-matrix K o h n method, J. Chem. Phys. 196 (1992) 1982-1991. [22] P. Comtois, J . Kneller, S. Nattel, Of circles and spirals: Bridging the gap between the leading circle and spiral wave concepts of cardiac reentry, Europace 7 (2005) S10-S20. 25 Bibliography [23] M . Courtemanche, R . J . Ramirez, S. Nattel, Ionic mechanisms underlying human atrial action potential properties: Insights from a mathematical model, A m . J . Physiol, 275 (1998) H301-H321 [24] T . R. Chay, Proarrhytmic and antiarrhythmic actions of ion channel blockers on arrhythmias in the heart: model study, A m . J . Physiol. 271 (1996) H329-H356. [25] E . M . Cherry, H . S. Greenside, C . S. Henriquez, A space-time adaptive method for simulating complex cardiac dynamics, Phys. Rev. Lett. 84 (2000) 1343-1346. [26] D . J . Christini, L . Glass, Introduction: Mapping and control of complex cardiac arrhythmias (Focus issue), Chaos 12 (2002) 732-739. [27] J . Davidenko, A . V . Pertsov, R. Salomonsz, W . Baxter, J . Jalife, Stationary and drifting spiral waves of excitation in isolated cardiac muscle, Nature, 355 (1992) 349-351. [28] C . Diks, B . Hoekstra, J . Degoede, Spiral wave dynamics, Chaos Sol. Frac., 5 (1995) 645-660. [29] D . R . Durran, Numerical methods for wave equations in Geophysical fluid dynamics, Springer. [30] J . C . Eilbeck, A collocation approach to the numerical calculation of simple gradients in reaction-diffusion systems solutions, J. M a t h . Biol. 16 (1983) 233-249. [31] I. G . Efimov, V . Krinsky, J . Jalife, Dynamics of rotating vortices in the Beeler-Reuter model of cardiac tissue, Chaos, Sol. Frac. 5 (1995) 513-526. [32] I. R. Epstein, J . A . Pojman, A n introduction to nonlinear chemical dynamics, (Oxford University Press, New York ,1998). [33] F . Fenton, A . Karma, Vortex dynamics in three dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation, Chaos 8 (1998) 20-47. [34] F . Fenton, E . Cherry, H . M . Hastings, S. J . Evans, Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity, Chaos 12 (2002) 852-892. 26 Bibliography [35] R . A . Fisher, The wave of advance of advantageous genes, A n n . Eugenics 7 (1937) 355-369. [36] R. Fitzhugh, Impulses and physiological states in theoretical models of nerve membrane, Biophy. J . , 1 (1961) 445-466. [37] B . Fornberg, A Practical Guide to Pseudospectral Methods (Cambridge University Press, Cambridge, 1996). [38] J . Gazdag, J . Canosa, Numerical solution of Fisher's equation, J . A p p l . Prob. 11 (1974) 445-457. [39] M . Gomez-Gesteira, A . P . Munuzuri, V . Perez-Munuzuri, V . PerezVillar, Boundary imposed spiral drift, Phys. Rev. E , 53 (1996) 54805483. [40] D . Gottlieb and S. A . Orszag, Numerical Analysis of Spectral Methods ( S I A M , Philadelphia, 1977). [41] D . Gottlieb, J . S. Hesthaven, Spectral methods for hyperbolic problems, J. Comput. Appl. Math., 128 (2001) 83-131. • [42] D . Gottlieb, J.S. Hestheaven, Spectral methods for hyperbolic problems, J . Comput. Appl. Math. 128 (2001) 83-131. [43] G . Gottwald, A . Pumir, V . Krinsky, Spiral wave drift induced by stimulating wave trains, Chaos, 11 (2001) 487-494. [44] S. A . Gourley, Travelling front solutions of a nonlocal Fisher equation, J . Math. Biol. 41 (2000) 272-284. [45] R. A . Gray, J . Jalife, Ventricular fibrillation and atrial fibrillation are two different beasts, Chaos 8 (1998) 65-77. [46] R. Gray, Termination of spiral wave breakup in a Fitzhugh-Nagumo model via short and long duration stimuli, Chaos 12 (2002) 941-951. [47] S. G r i l l , V . S. Zykov, S. C . Miiller, Spiral wave dynamics under pulsatory modulation of excitability, J . Phys. Chem. 100 (1996) 1908219088. [48] T . Hagstrom, H . B . Keller, The numerical calculation of travelling wave solutions of nonlinear parabolic equations, S I A M J . Sci. Stat. Comput. 7 (1986) 978-988. 27 Bibliography [49] W . Huang, Y . Ren, R. D . Russell, Moving mesh partial differential equations ( M M P D E s ) based on the equidistribution principle, S I A M J . Numer. A n a l . 31 (1994) 709-730. [50] W . Huang, R. D . Russell, A moving collocation method for solving time dependent partial differential equations, Appl. N u m . M a t h 20 (1996) 101-116. [51] T . Ikeda, T . Uchida, D . Hough, J. J . Lee, M . C . Fishbein, W . J . Mandel, P. S. Chen, H . S. Karagueuzian, Mechanisms of spontaneous termination of functional reentry in isolated canine right atrium, Circulation, 94 (1996) 1962-1973. [52] T . Ikeda, M . Yashima, T. Uchida, D . Hough, M . C . Fishbein, W . J . Mandel, P. S. Chen, H . Karagueuzian, Attachment of meandering reentrant wavefronts to anatomic obstacles in the atrium, Circ. Res. 81 (1997) 753-764. [53] W . B . Jones, J . J . O'Brien, Pseudo-spectral methods and linear instabilities i n reaction-diffusion fronts, Chaos 6 (1996) 219-228. [54] A . K a r m a , Meandering transition in two-dimensional excitable media, Phys. Rev. Lett., 65 (1990) 2824-2828. " [55] A . Karma, Scaling regime of spiral wave propagation in single-diffusive media, Phys. Rev. Lett., 68 (1992) 397-400. [56] J . P. Keener, A geometrical theory for spiral waves in excitable media, S I A M J . A p p l . Math. 46 (1986) 1039-1056. [57] J . P. Keener, J . Sneyd, Mathemathical Physiology, Interdisciplinary Applied Mathematics (Springer, New York, 1998) [58] J . P. Keener, K . Bogar, A numerical method for the solution of the bidomain equations in cardiac tissue, Chaos 8 (1998) 234-241. [59] D . A . Kessler, H . Levine and W . N . Reynolds, Theory of the spiral core in excitable media, Physica D 70 (1994) 115-139. [60] A . G . Kleber, Y . Rudy, Basic mechanisms of cardiac impulse propagation and associated arrhythmias, Physiol Rev, 84 (2004) 431-488. [61] V . Krinsky, A . Pumir, Models of defibrillation of cardiac tissue, Chaos 8 (1998) 188-203. 28 Bibliography [62] L . Edelstein-Keshet, Mathematical models in biology, Classics in A p p . Math., S I A M . [63] K . Lee, Wave pattern selection in an excitable system, Phys. Rev. Lett. 79 (1997) 2907-2910. [64] A . Z. Lei, et al, Spiral wave maintained by noise, Fluctuation and Noise Letters, 4 (2004) 447-452. [65] Y . X . L i , J . Rinzel, Equations for InsP receptor mediated [ C a ] i oscillations derived from a detailed kinetic model: a Hodgkin-Huxley like formalism, J . Theor. B i o l , 166 (1994) 461-473. 2+ 3 [66] S. L i , L . Petzold, Y . Ren, Stability of moving mesh systems of partial differential equations, S I A M J . Sci. Comput. 20 (1998) 719-738. [67] C h . Luo, Y . Rudy, A model of the ventricular cardiac action potential, Circ. Res. 68 (1991) 1501-1526. [68] C h . Luo, Y . Rudy, A dynamic model of the cardiac ventricular action potential, Circ. Res. 74 (1994) 1071-1096. [69] G . Mansell, W . Merryfield, B . Shizgal, U . Weinert, A comparison of differential quadrature methods for the solution of partial differential equations, Comp. Meth. A p p l . Mech. Engrg. 104 (1993) 295-316. [70] D . Margerit, D . Barkley, Selection of twisted scroll waves in threedimensional excitable media, Phys. Rev. Lett. 86 (2001) 175-178. [71] R . E . Mickens, A best finite-difference scheme for Fisher's equation, Numer. Meth. Part. D . E . 10 (1994) 581-585. [72] R . E . Mickens, Relation between the time and space step-sizes in nonstandard finite difference schemes for the Fisher equation, Numer. Meth. Part. D . E . 13 (1997) 51-55. [73] I. Mitkov, I. Aranson, D . Kessler, Meandering instability of a spiral interface i n the free boundary limit, Phys. Rev. E , 54 (1996) 6065-6069. [74] L . S. Mulholland, Y . Qiu, D . M . Sloan, Solution of evolutionary partial differential equations using adaptive finite differences with pseudospectral post-processing, J . Comput. Applied Math. 131 (1997) 280-298. [75] J . Nagumo, S. Arimoto, S. Yoshizawa, A n active pulse transmision line simulating nerve axon, Proc I R E , 50 (1962) 2061-2070. 29 Bibliographylid} S. Nattel, J . Kneller, R. Zou, J . Leon, Mechanisms of termination of atrial fibrillation by class I antiarrhythmic drugs: Evidence from clinical, experimental, and mathematical modeling studies, J . Cardiovasc. Electrophysiol. 14 (2003) S133-S139. [77] K . T . Ng, R. Yan, Three-diemnsional pseudospectral modelling of cardiac propagation in an inhomogeneous anisotropic tissue, Med. Biol. Eng. Comput. 41 (2003) 618-624. [78] D . Noble, S. J . Noble, A model of S. A . node electrical activity using a modification of the DiFrancesco-Noble(1984) equations, Proc. Roy. Soc. Lond. 222 (1984) 295-304. [79] A . Nygren, C . Fiset, L . Firek, J . W . Clark, D . S. Lindblad, R . B . Clark, W . R . Giles, Mathematical model of an adult human atrial cell: the role of K currents i n repolarization, Circ. Res. 82 (1998) 63-81. + [80] D . Olmos, B . Shizgal, A Spectral Method of Solution of Fisher's Equation, J . Comp. Appl. Math. 193 (2006) 219-242. [81] N . Otani, A primary mechanism for spiral wave meandering, Chaos 12 (2002) 829-842. [82] A . Panfilov, P. Hogeweg, Spiral breakup in a modified FitzhughNagumo model, Phys. Lett. A 176 (1993) 295-299. [83] A . V . Panfilov, J . P. Keener, Effects of high frequency stimulation on i cardiac tissue with an inexcitable obstacle, J . theor. Biol, 163 (1993) 439-448. [84] N . Parekh, S. Puri, A new numerical scheme for the Fisher equation, J. Phys. A : Math. Gen. 23 (1990) L1085-L1091. [85] A . M . Pertsov, J . M . Davidenko, R . Salomonsz, W . T . Baxter, J . Jalife, Spiral waves of excitation underlie reentrant activity in isolated cardiac muscle, Circ. Res. 72 (1993) 631-650. [86] V . Petrov, S.K. Scott, K . Showalter, Excitability, wave reflection and wave splitting in a cubic autocatalysis reaction diffusion system, P h i l . Trans. R. Soc. Lond. A 347 (1994) 631-642. [87] R . Peyret, Spectral Methods for Incompressible Viscous Flow, Springer, New York, 2002. 30 Bibliography Y. Qiu, D . M . Sloan, Numerical solution of Fisher's equation using a moving mesh method, J . Comput. Phys. 146 (1998) 726-746. Rizwan-Uddin, Comparision of the nodal integral method and non standard finite-difference schemes for the Fisher equation, S I A M . J . Sci. Comput. 22 (2001) 1926-1942. J. Roessler, H . Hiissner, Numerical'solution of the 1 + 2 dimensional Fisher's equation by finite elements and the Galerkin method, M a t h . Comput. Modelling, 25 (1997) 57-67. B . J . Roth, A r t Winfree and the bidomain model of cardiac tissue, J . Theo. Biol., 230 (2004) 445-449. B . D . Shizgal, H . Chen, The quadrature discretization method ( Q D M ) in the solutions of the Schrodinger equation with nonclasical basis functions, J . Chem. Phys., (1995) 4137-4150. G. Skinner, H . L . Swinney, Periodic to quasiperiodic transition of chemical spiral rotation, Phys. D , 48 (1991) 1-16. F . Stenger, Numerical methods based on sine and analytic functions, Springer Series i n Comp. Math., 20 (1993) 91-96. S. Takagi, et al, A physical approach to remove anatomical reentries: a bidomain study, J . Theo. Biol., 230 (2004) 489-497. A . S. Tang, H . Ross, C . S. Simpson, L . B . Mitchell, P. Dorian, R . Goeree, B . Hoffmaster, M . Arnold, M . Talajic, Canadian cardiovascular society / Canadian heart rhythm society position paper on implantable cardioverter defibrillator use in Canada, The Canadian J . Cardiol. 21 (2005) 11A-18A. S. Tang, R. O. Weber, Numerical study of Fisher's equation by a Petrov Galerkin finite element method, J . Aust. Math. Soc. B 33 (1991) 27-38. J. A . Trangenstein, C. K i m , Operator splitting and adaptive mesh refinement for the Luo-Rudy I model, J . Comp. Phys. 196 (2004) 645-679 K . H . W . J . ten Tusscher, A . V . Panfilov, Influence of nonexitable cells on spiral breakup i n two-dimensional and three-dimensional excitable media, Phys. Rev. E , 68 (2003) 062902(1-4) [100] E . H . Twizell, Y . Wang, W . G . Price, Chaos free numerical solutions of reaction-diffusion equations, Proc. R. Soc. Lond. A . 430 (1990) 541-576. 31 Bibliography [101] J . Tyson, J . Keener, Spiral waves in a model of myocardium, Physica D, 29 (1987) 215-222. [102] J . Tyson, What everyone should know about the belousovZhabotinsky reaction, Lecture notes in biomathematics, Springer, V o l . 100 (1994) 569-587 [103] Vander, Sherman, Luciano, Human Physiology, McGraw H i l l (2001) [104] B . Van der Pol, O n relaxation oscillations, P h i l . Mag., 2 (1926). 978992. [105] B . Van der Pol, J . Van der Mark, The heartbeat as a relaxation oscillation and an electrical model of the heart, (1926) [106] G . D . Veenhuyzen, C . S. Simpson, H . Abdollah, A t r i a l C M A J 171 (2004) 755-760. fibrillation, [107] E . J . Vigmond et. al., Computational Techniques for solving the bidomain equations in three dimensions, I E E E Trans. Biomed. E n g . 49 (2002) 1260-1269. [108] G . W . Wei, Solving quantum eigenvalue problems by discrete singular convolution, J . Phys. B : A t . Moi. Opt. Phys. 33 (2000) 343-352. [109] E . T . Whittaker, O n the functions which are represented by the expansions of the interpolation theory, Proc. Roy. Soc. Edinburgh 35 (1915) 181-194. [110] A . T . Winfree, When time breaks down, Princeton University, 1987. [Ill] A . T . Winfree, Varieties of spiral wave behavior: A n experimentalist approach to the theory of excitable media, Chaos, 1 (1991) 303-334. [112] A . T . Winfree, Evolving perspectives during 12 years of electrical turbulence (Focus issue), Chaos 8 (1998) 1-19. [113] F . X i e , Z. Q u , J . Weiss, A . Garfinkel, Coexistence of multiple spiral waves with independent frequencies in a heterogeneous excitable medium, Phys. Rev. E , 63 (2001) 031905 1-4. [114] F . X i e , Z. Qu, J . N . Weiss, A . Garfinkel, Interactions between spiral waves with different frequencies in cardiac tissue, Phys. Rev. E . 59 (1999) 2203-2205. 32 Bibliography [115] H . Yang, B . Shizgal. Chebyshev pseudospectral multidomain technique for viscous flow calculation, Comput. Methods A p p l . Mech. Eng. 118 (1994) 47-61. [116] H . Yang, B . Seymour, B . Shizgal, A Chebyshev pseudospectral multidomain method for steady flow past a cylinder, up to Re=150, Comput. fluids 23 (1994) 829-851. [117] Ye. A . Yermakova, A . M . Pertsov, Interaction of rotating spiral waves with a boundary, Biophys., 31 (1986) 932-940. [118] B . L . Zaret, M . Moser, L . S. Cohen, Yale University school of medicine heart book, (Hearst books, New York, 1992). [119] Z. Zhan, K . T . Ng, Two-dimensional Chebyshev pseudospectral modelling of cardiac propagation, Med. Biol. Eng. Comput. 38 (2000) 311318. [120] S. Zhao and G . W . Wei, Comparison of the discrete singular convolution and three other numerical schemes for solving Fisher's equation, S I A M J . Sci. Comput. 25 (2003) 127-147. [121] Z. Zhu, H . X i n , Dynamics of spiral waves under the modulation of noise pulses, Phys. Rev. E , 64 (2001) 056124(1-5). [122] D . P. Zipes, Epidemiology and mechanisms of sudden cardiac death, The Canadian J . Cardiol. 21 (2005), 37A-40A. [123] X . Zou, Delay induced traveling wave fronts in reaction diffusion equations of K P P - F i s h e r type, J . Comput. A p p l . M a t h . 146 (2002) 309-321. 33 Chapter 2 A Pseudospectral Method of Solution of Fisher's Equation. 2.1 Introduction Spectral methods have long been known to provide very accurate and rapidly convergent solutions of partial differential equations with smooth solutions [132, 135, 164]. These methods generally provide an exponential convergence of the solution versus the number of collocation points. In recent years, spectral methods have also been used for the solution of differential equations with solutions that resemble shock waves or fronts typical of hyperbolic partial differential equations [146]. There is an ongoing interest to further adapt spectral methods to differential equations with rapidly varying and propagating solutions. The purpose of the present chapter is to apply a spectral method to the solution of Fisher's equation ( F E ) , which was originally proposed by Fisher [143] as a model for the spatial and temporal propagation of a virile gene in an infinite medium. It is a one-dimensional reaction diffusion model for the evolution of the infected population, U(x',t'), with a quadratic reactive term corresponding to logistic growth. The equation is defined by 2 where t' is the time and x' € ( — 0 0 , 00) is the position. The diffusion and reactive processes are parameterized by a diffusion coefficient, D, and a reactive rate coefficient, k, respectively. We consider solutions to Eq.(2.1) subject to the following initial and boundary conditions, lira U{x\t') = 0 (2.2) x'—*oo A version of this chapter has been published. D. Olmos, B. Shizgal, A Pseudospectral Method of Solution of Fisher's Equation, J. Comp. Appl. Math. 193 (2006) 219-242. 2 34 Chapter 2. A Pseudospectral Method of Solution of Fisher's lim U(x',t') = 1 = U (x') Equation. x'—*—OO U(x',0) 0 It has been shown that with the appropriate boundary conditions F E will support travelling waves of the form U(x' — c't') moving in the positive a;-direction, provided that the speed c! is greater than the critical value c' = 2\/kD. Equation (2.1) is the simplest reaction diffusion equation employed to model many problems in mathematical biology [160]. W i t h the change of variables, min / k\ t = kt' 1 / 2 x = x'(^r) (2.3) E q . (2.1) becomes and travelling wave solutions exist for dimensionless speeds c > 2, [152]. The mathematical properties of F E have been studied extensively and there have been numerous discussions in the literature. Excellent summaries have been provided in [133, 153, 160]. One of the first numerical solutions was presented by Gazdag and Canosa [145] with a pseudo-spectral approach. Implicit and explicit finite differences algorithms have been reported by different authors such as Parekh and Puri [163] and Twizell et al [173]. The works of Mickens [157, 158] considered time stepping aspects for finite difference algorithms. The work by Hagstrom and Keller [149], where the main goal was to develop asymptotic boundary conditions, considered a centered finite difference algorithm. Rizwan-Uddin [167], compared a nodal method with a non standard finite difference scheme. A Galerkin finite element method was used by Tang and Weber [172] whereas Carey and Shen [136] employed a least-squares finite element method. A collocation approach based on Whittaker's sine interpolation function [132, 171] was also considered by Al-Khaled [125], and Zhao and Wei [182]. The work by Gourley [147] considered a nonlocal form of F E . The study by Zou [183] was concerned with another modified form of F E including time delay and the work by Roessler and Hiissner [168], considered finite elements for a two dimensional FE. Solutions of F E exhibiting propagating fronts thus possess features similar to those of shock waves that arise with hyperbolic equations. There are also special interesting features of the solution in terms of the relation between 35 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. the speed of the wavefront and the behavior of the solution at infinity [148, 153]. Larson [153] and Hagan [148], proved that for any initial condition of Eq. (2.4) such that U (x) 0 Ae~ Px then, U(x, t) evolves as a wave front with speed given by 0</3< 1 2 P>1 There is also an interest concerning the instability of the solution to small perturbations in the solution particularly when U(x,t) w 0 as discussed by different authors [137, 145, 148, 165]. Canosa [137] proved stability of the solution to perturbations of compact support, whereas instability occurred when the perturbation vanished at infinity. This property plays a fundamental role when E q . (2.4) is solved numerically [145]. A closely related problem is to consider a modified form of F E introduced by L i et al [154], for which the nonlinear reactive term is made arbitrarily larger than the diffusion term for the purpose of testing algorithms. This modified F E is given by, with initial and boundary conditions similar to E q . (2.2), and the reaction rate coefficient is generally chosen so that p >> 1. A particular solution of Eq. (2.5) considered by L i et al [154], was found by Ablowitz and Zepetella [124]. It has the form of a travelling wave front, and is given by U(x,t) = (l + 1 2 (2.6) e x p ^ - ^ t ) ) which travels with constant speed c = 5 ^ / | . The initial condition of E q . (2.5) is clearly given by U (x) = U(x,0). Solutions of E q . (2.5) with large p have been referred to by Zhao and Wei [182] as super speed wave (SSW) types. W i t h an increase in p, the propagating front steepens and this presents a challenging numerical problem to both resolve and track the front. This rescaled version of F E was considered by L i et al [154] in their study of moving mesh strategies in finite difference methods of solution of partial differential equations. They commented that moving mesh methods are not recommended for such reaction diffusion problems for which 0 36 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. the diffusion term is much smaller than the reactive term. Subsequently, Qiu and Sloan [165] carried out a detailed comparison of different moving mesh strategies [150, 151] and concluded that these methods are not easily adapted for equations analogous to F E with steep fronts. The authors have also applied these methods to Burger's equation [130, 131, 159] with similar shock like solutions. In the course of these studies on F E , numerous workers have also reported a sensitivity of the numerical solutions to perturbations owing to numerical noise i n the solutions or round-off errors resulting i n instabilities. One of the objectives of the present chapter is to apply spectral methods to the solution of F E . A second objective is to establish a relationship between the problems given by E q . (2.4) and E q . (2.5), which we demonstrate lead to equivalent numerical problems. In the first instance, we consider a single fixed domain suitably truncated to the interval [xi,Xft] and apply a spectral method based on Chebyshev-Lobatto points. This approach gives good results but only for relatively small values of p. The difficulty encountered for larger values of p is traced to round off errors.in the application of the second derivative spectral matrix operator to the solution. The occurrence of round-off errors for these matrix derivative operators have been well documented [126, 127, 128]. We examine i n detail the origin of these round-off errors i n the application of spectral methods to F E . A multidomain approach developed by Shizgal and coworkers [155, 180, 181] for the solution of Burgers equation and the Navier Stokes equation was then employed for the present work. We obtain accurate stable solutions of F E for relatively large values of p, with the appropriate division of the domain [XL,XR] into subdomains. In Section 2.2, we outline the spectral method applied to F E in one domain and analyze the problem of round-off errors. The multidomain approach is presented i n Section 2.3. In Section 2.4, we present a comparison of the present approach with the one based on Whittaker's sine interpolation [132, 140, 179] employed by Zhao and Wei [182] and by Wei [174] that they refer to as the discrete singular convolution (DSC). 2.2 2.2.1 Chebyshev-Lobatto Spectral Approach to Fisher's Equation The modified F E ; scaling the p dependence. Several groups ( L i et al [154], Qiu and Sloan [165] and Zhao and Wei [182]), have employed the modified form of F E (Eq. 2.5) for the S S W situation. 37 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. The exact solution of the modified F E exhibits a shock-like front for large p and speed c = 5^/p/6. For the infinite spatial domain, the rapidly varying shock front is considered to be stiff with the stiffness depending on p . In the next section, we propose a numerical solution of Eq. (2.5) which involves the expansion of the solution in Chebyshev polynomials orthogonal on [—1,1]. This requires that the boundary conditions be applied on a truncated domain [XI,XR\. AS a consequence of the use of a finite domain in the numerical solution, there is an important dependence of the solution on p and the width of the interval considered. O n this truncated interval, we consider the differential equation for V(x,t), V = V T XX + V(1 - V) (2.7) P with the initial condition given by V(x,0) = V (x) x £ [ x , x ] 0 L (2.8) R and boundary conditions at the ends of the truncated interval, that is V{x ,t) = L V{x ,t) = R 1 . 0, i € [ 0 , T ] ^' . y j where VQ{X) is given by U(x,0) i n E q . (2.6) but over [XL,X ]. From E q . (2.6), the dependence on the parameter p can be removed with a scaling of the space and time variables, i.e. R z = ^/px T= p t (2.10) and E q . (2.7) can be written with the dependence of p occurring only at the end points of the computational domain, that is V = V T ZZ + V(l-V) (2.11) with V(z,0) = V(^px ,r) = V{Jpx ,T) = L R V (z) z£ 1 0, r€[0,pT] 0 \yjp~XL, yfpxR] (2.12) Due to the boundary condition V ( ^ / p x , T ) = 0, and i n order to preserve a good approximation of the solution. T is taken as a time before the wavefront hits the right boundary at x . B y virtue of the transformations Eq. (2.10), the solution of Eq. (2.7) for p = p i on [ x i , x ] is equivalent to E q . (2.11) on [y/pixi, ^ / p i X \ . Similarly, the R R R R 38 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. solution for p = p on [x' , x' ] is equivalent to E q . (2.11) on [^/p2x' , s/p2x' ]. The solutions for different p are equivalent to each other provided that 2 L R L R (2.13) are satisfied. If we consider the problem given by E q . (2.7) with p = p on [x' , x' ], and a numerically equivalent problem given by pi on [xi, XR] (i.e. satisfying (2.13)), then, if pi is modified to p,, the problem with p = p on [ r E j r ^ y and the modified problem with p = ~p on [XL, XR], are no longer equivalent unless XL and XR are modified in order to satisfy condition (2.13). In other words, relation (2.13) establishes a property between p and the length of the interval L = XR — XL, where modifying p or L and fixing the other one, represents the same problem as long as relation (2.13) is satisfied. It follows that increasing the value of p leads to a steeper front only if we consider a fixed numerical domain. This aspect of F E appears not to have been mentioned previously in the literature. The case of S S W studied in [154, 165, 182] has been considered previously by different authors with p = 1 and different interval lengths. As an example, Qiu and Sloan [165] and Zhao and Wei [182] considered the computational domain with end points given by XL = —0.2 and XR = 0.8 with p = 10 whereas Gazdag and Canosa [145] considered p — 1, XL = 0, XR = 140. According to (2.10), the problem considered by Gazdag and Canosa [145], is equivalent to p = (140) = 19,600, XL = 0 and XR = 1. Another example is the problem studied by Parekh and Puri [163] where they chose p = 1, XL = 0 and XR = 300 which is equivalent to p = (300) = 90,000 with XL = 0 and XR = 1. However, these previous works did not consider the case of solving E q . (2.7), with the initial and boundary conditions given by Eqs. (2.8) and (2.9). Instead, properties such as the numerical stability of algorithms, or features of reaction diffusion processes, were considered. A s an example of previous numerical studies, we mention the works in [125, 136, 167, 172] that analyzed the phenomena of reaction and diffusion on F E for different initial conditions. Gazdag and Canosa [145] based on Canosa's previous work [137], considered initial conditions of E q . (2.4) with an asymptotic behavior at infinity such that the speed is faster than the minimum, c = 2. However, due to their proposed numerical scheme, they obtained oscillations in their numerical solution near the right boundary and at the onset of the front, giving an unstable solution. The instability problem was solved considering V(x,t) = 0 whenever | V ( x , t)| < e, where e is some small quantity. However, this assumption leads to the loss of the initial theoretical speed, and to the convergence to the minimum speed front c = 2. 2 L R 2 l 4 2 2 39 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. The variable p in the numerical problem (2.11) with (2.12), does not play the role of a reaction rate coefficient but a scaling over the numerical domain as considered in [125, 136, 145, 163, 167, 172]. However, as the work presented in this chapter follows from the work by L i et al [154], Q i u and Sloan [165] and Zhao and Wei [182], we will consider the numerical problem given by (2.7) with (2.8). ' 2.2.2 Pseudospectral solution of the modified F E . Expanding U(x,t) or V(x,t) i n the Chebyshev Gauss polynomials, T (x), which are orthogonal with respect the weight function w(x) = (1 — i ) / on the interval [—1,1], gives, k 2 cfc 7_i w^TkixW^dx = ±w5 kil - 1 2 (2.14) where c = 1 for all k except for CQ = 2. The Lobatto quadrature points and weights associated with the Chebyshev polynomials are given by Xi = — cos ( ^ ) and the weights are Wi = fr for all i except WQ = WN = ^ [132, 135, 164]. These points and weights provide the approximate quadrature, k N 1 / w(x)f(x)dx~Yl *f( i) w x (- ) 2 15 i=0 where N+l is the number of points. Since any piecewise continuous function, / € L [ 0 , 1 ] can be expanded i n a Chebyshev polynomial series that is convergent in the norm, we have 2 JV f(x)^f (x) N = J2 kTk(x) a (2.16) fc=0 where a = — [ w(x)f(x)T (x)dx k k (2.17) CfcTT 7 _ i W i t h Eqs. (2.15) to (2.17) we obtain the interpolation algorithm f (x)^^Ij^)f(^) (- ) 2 N 18 j=0 where the interpolating polynomials, Ij(x), are given by 2 ^) = N E ^Tk{xj)T (x) k fc=0 (2.19) 40 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. where VQ = — 1/2 and v = 1 ii k ^ 0, N. The nth derivative of f(x) at the quadrature points is then given approximately by k f%\x )^lf(x )f(x ) k k (2.20) 3 3=0 If we denote by f, the vector of the function evaluated at the ChebyshevLobatto points, E q . (2.20) can be rewritten as f(") = D ( n ) •f (2.21) and thus, the second derivative matrix is then identified with This is the basis for the Chebyshev pseudospectral method. It is now straightforward to apply the pseudospectral method to the modified F E , by approximating the second derivative operator as in E q . (2.22). From Eq. (2.7) and the linear transformation from [XL,XR] to [-1,1], the following system of ordinary differential equations is obtained. ^=Aj2D^V V (l-V ) j+P l l • (2.23) where A = 4/ (XR — XL) , Vi = V(xi,t) and the time dependence is considered implicitly in Vj. The spectral derivative matrices were calculated with the M A T L A B suite of programs developed by Weideman and Reddy [177]. The set of ordinary differential equations was integrated with a Runge-Kutta integrator i n M A T L A B subject to the boundary conditions, V(XL) = 1 and V(xR) = 0 for all t. 2 2.2.3 Numerical results. The formalism i n Sections 2.2.1 and 2.2.2 was applied to a study of the behavior of the numerical solution of F E versus p. The behavior for a fixed interval for different values of p is considered. In the present study, Uo(x) will refer to the initial condition depending on p, that is E q . (2.6) with t — 0, unless it is otherwise indicated. We choose XL = —0.2, XR — 0.8 as done by Q i u and Sloan [165] and by Zhao and Wei [182] and vary p. In F i g . 41 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. 2.1, we compare the analytic solution U(x,t) (Solid curves) evaluated at the Chebyshev collocation points with the numerical solution Vi (symbols), at different t. The values of p are (A) 2,000, (B) 5,000 and (C,D) 10,000, with the observation that p = 10,000 was considered i n the previous works [154, 165, 182]. The steepening of the front on increasing p for a fixed interval is clearly seen by comparing Figs. 2.1C and 2.ID with Figs. 2.1A and 2.IB. In Fig. 2.1A, the good agreement between the numerical results and the analytic solution with p — 2,000 and 40 Chebyshev points is shown. As the initial profile (Eq. (2.6), t = 0) behaves asymptotically as given by Eq. (2.29), the wave speed is constant at c = 5-y//o/6. The wave speeds calculated from the numerical solutions agree with the theoretical speed to 4-6 significant figures depending on the value of p. In order to get a stable and accurate solution for p = 5,000 shown in Fig. 2.IB, the number of points had to be increased from N = 40 to N = 64. F i g . 2.1C for p = 10 and iV = 64, shows an instability (negative part) that appears at x ~ 0.7 for t = 0.002. When N is increased to 150, a stable solution is obtained as shown in F i g . 2.ID. To further validate the numerical method used, we have considered the solution of the diffusion equation (p = 0) but with an exponentially decaying initial profile (UQ{X),P = 10 ). We find that instabilities at small times are quickly damped by diffusion. However, for F E with p — 1.5 x 10 , a stable solution could not be obtained even with N as large as 250. 4 4 4 It is of considerable interest to understand the origin of this instability, which we attribute to round-off error associated with the application of the second derivative operator on the solution at some time step. This is a problem common to many applications of spectral methods as recently discussed [126, 127, 128]. We therefore study the error in the numerical computation of the second derivative of the initial condition Un(x). In F i g . 2.2, we show the relative error defined by d U (x) dx 2 0 RE{xk) = l o g (2.24) 2 10 d U (x) dx 2 0 2 I I x=x k versus x , for p equal to 10 and 1.5 x 10 and N = 64, 128 and 200. These results clearly demonstrate that the error is larger at the right boundary than at the left boundary. Furthermore, for x > 0.1 where U(x) « 0, the error increases exponentially with respect to x, reaching its largest value at the right boundary. In Fig. 2.3, the relative errors at the right and left boundaries versus N are compared. As in Fig. 2.2, the error is bigger at XR than at xi. It decreases 4 4 k 42 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. with N until about N = 200 for p = 10 (iV = 230 for p = 1.5 x 10 ) and then increases slowly with a further increase in iV. Thus, it is clear that an increase in ./V decreases the error and as a consequence a stable solution can be obtained, F i g . (2.ID), iV = 150 whereas an instability occurs with a smaller N = 64 for Fig. (2.1C). For p = 10 , the improvement of the accuracy can be done up to iV = 200, where for larger values of N round-off error begins to become significant in the calculations. When p = 1.5 x 10 , N = 250 is not large enough to obtain a stable solution. From F i g . 2.3B, it is noticed that values of N greater than 250 do not provide an improvement of the accuracy due to the effects of round off error that have become significant. It is important to mention that round off errors are present even for relatively small values of p. However, since the amplitude of the error is small throughout the time integration considered, the accumulated error does not play a crucial role, and a stable solution is obtained. The main contribution of the round off error, is due to (i) the alternation in sign and the magnitude of the elements of the second derivative operator D( ), and (ii) the small values of d Uo(x )/dx relative to Uo(xi). Uo(x) varies approximately exponentially with x for x « x and UQ(X) w 0 in the region near the right boundary. Consequently, d Uo(x) / dx is very small. The behavior of D UQ(X) / dx versus p at the boundary points is shown in Fig. 2.4A. Whereas the value of the second derivative at x^ remains almost constant, the value of d Uo(x )/dx is small and decreases rapidly with an increase in p. The numerical approximation of the second derivative at the right boundary, d Uo(x )/dx involves the summation of the form, 4 4 4 4 2 2 2 R R 2 2 2 2 2 2 R 2 2 R N S = Y1 NI V( I<) D N U (- ) X 2 25 k=0 where the subindex A'" refers to the Nth row and D ' are the elements in the Nth row of the second derivative matrix operator (2.22). It is well-known that (i) the size of the largest element of increases as A ^ , whereas the smallest increases as N and (ii) the elements alternate in sign [164]. These features of the elements of D( ) are confirmed in F i g . (2.4B) for N = 250. As a result, the sum (2.25) consists of adding alternately positive and negative terms, due to the alternation of the signs of the elements in D^^. For p = 15,000 and N = 250, d U (x )/dx is of the order of 1 0 ~ (Fig. 2.4A), whereas the terms in (2.25) take values from 1 0 ~ to 10 i.e., the terms in (2.25) become relatively large compared to d Uo(x )/dx . K N t 4 2 2 2 2 2 0 31 R 27 4 2 2 R 43 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. The key point is observed i n the last two elements of the sum, E q . (2.25). It is the difference between two numbers of the order of 10 , that should approximate a quantity of the order of 1CT . This subtraction leads to a poor approximation of the second derivative. In the same way, the approximated second derivative over the rest of the points near XR, presents the same problem. It has to be mentioned that this problem is less severe for 4 31 d U (x )/dx , 2 2 0 L as \og (d U {x )/dx ) 2 w 2 0 L « 0. From the previous analysis, the smallness of d Uo(xR)/dx , is the main source of round off error at XR. A s a larger domain is considered, an increase in XR and/or decrease in xi, d Uo(xR)/dx decreases, giving a larger error at XR. Then, it is clear that the effects of round off error are greater for a larger domain XR — XL- In order to reduce the effects of round off error a smaller domain has to be considered. As we will see i n the next section, we partition the main domain into smaller domains, contributing to a reduction in the round off error. From relation (2.13), it follows that increasing the length of the interval and fixing p, is equivalent to an increase of p and fixing the length of the interval. Then,,the problem of round off error will be present with the same magnitude whether p or the length of the interval are varied, according to relation (2.13). The main consequence due to the round off error problem is the development of unwanted oscillations at some particular time step in the integration of E q . (2.23). This type of oscillation is similar to the one reported by Zhao and Wei [182], for the Fourier pseudo-spectral method. The difference between the oscillations reported by Zhao and Wei [182] and the instabilities in this chapter with the Chebyshev-Lobatto collocation, is the location of the oscillations. Whereas in [182] the higher amplitude oscillations are at the foot of the wave front, in this chapter they are located i n a neighborhood of XR as discussed previously. 2 2 2.3 2 2 Chebyshev-Lobatto Multidomain Spectral Method In order to overcome the numerical round-off errors discussed in the previous section, we employ a multidomain approach used previously [141, 142, 155, 162, 180, 181]. This involves splitting the domain [XL,XR], into K subdomains denoted by 44 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. where each subinterval has length L and is discretized with M+l Chebyshev points as shown i n F i g . 2.5. The first two quadrature points of the interval Ifi+i coincide with the last two points of the interval 1^. The overlap of these points of neighboring domains is possible provided that each subinterval is of the same length L: K+ (l-jV,)(l-cos(|.)) 2 and the same number of Chebyshev collocation points is used i n each interval. Thus, the new grid of points over the whole interval can be represented by, { k} = { 0> •••) \d-l x x x 0i M = x X = \i M-2> M-1 x X X> = 0 x + > ••• M) (2.26) X The overlap of the subdomains as described is very important for the correct construction of the derivative operators. In each sub-domain 1^, we proceed as before and have an equation similar to E q . (2.1), given by lir = ~dx r + p U " { 1 ~ U ^[a&a&L n te[0,t ] f (2.27) where is the solution over the fi interval, and XQ and x^ are the left and right boundaries of the u. interval, respectively. The discretized form of E q . (2.27) for each subinterval 7 is given by th M th M m jrD(fu? 4 dt (M X o) x + pU?(l-Un (2.28) j=0 where t/f = U (x^,t) and are the Chebyshev-Lobatto points over the interval 1^. The column vector U^(i), consists of joining the column vectors L/f according to E q . (2.26). >x 45 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. The first derivative matrix operator is then defined by, n n 1 -^1,0 ' -^M-1,0 ••• ••• 1 l,M-l u 1 ^l,M -^M-l.M-l n ^M-l.M D -^M-1,0 -^M-1,1 2 D= n D 2 ••• 2 ^ M - l . M - l D Ni ntii D \ n Ni 2 1,M ^M-1,M D Ni N D t D Ni and the remaining components are zero. The size of the derivative operator is (JV + 1) x (N + 1), with N = M+(M - 1)(K - 1). The multidomain approach is expected to provide more stable and accurate results. The round-off errors that occur in the application of D^ ) to the solution vector U ( i ) are smaller as the main interval has been subdivided. We demonstrate the success of the multidomain approach i n F i g . 2.6 where the successive profiles are shown for p = 10 . In this calculation, 140 subdomains with 20 collocation points in each subdomain were used. The time step was A t = 8 x 1 0 with an integration up to T = 3.3 x 1 0 such that Loo = 3.12 x 1 0 . The agreement between the numerical and analytic solutions is excellent. The value of p used here is significantly larger than the value used by numerous other researchers whose works have been cited. To further benchmark this method, we choose p = 10 for which we obtained accurate solutions with the single domain and N = 150. W i t h the multidomain method we study the effect of varying the number of domains and the number of points in each domain. The results are summarized in Table 2.1. It is clear from these results that there is a considerable improvement with the multidomain method. The best results for the case in Table 2.1 is with 15 subdomains and 10 points per domain. The error is an order of magnitude smaller than the one domain approach with the same total number of points. In Table 2.2, a similar analysis is carried out for p = 10 with N = 400. The single domain calculation leads to an instability in this case. The multidomain approach provides an excellent result with 20 domains and 20 points in each domain. 2 6 - 8 - 4 - 1 2 4 5 46 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. As discussed in the previous section, UQ(X) given by E q . (2.6) with t = 0, behaves as a negative exponential where UQ(X) « 0. More specifically, U (x) ~ 0 e~ V$ 2 x U (x) w 0 for (2.29) 0 Moreover, for values of x where UQ(X) « 0, the second derivative d Uo(x) / dx behaves like UQ(X), up to a factor depending on p, i.e., 2 ^ g ^ f t M x ) for 2 C/ (x)«0 (2.30) 0 O n the other hand, we know that the largest value of the sum in E q . (2.25) that is used also in the multidomain approach for each subdomain, is of the order of R = L 0 G L ° + {(XM-XO) ) 2 l o g l o ( 7 V 2 ) " °Zio(U (xo))\ 0 where from Eq. (2.30), the smallest value to be approximated is ^C/O(IM), which is of the order of E = log (2.31) 1l 10 (jfj - | log (U (x ))\ 10 0 M d Uo(xM)/dx 2 2 (2.32) In order to get a good approximation for E, the difference R — E given by should be reduced. For a fixed value of p and N, the first term of E q . (2.33) is a decreasing function of the length of the interval whereas the second is an increasing function. Then, for a fixed p and iV there is an optimal interval length to consider in the multidomain approach. For this last analysis, the region considered is where UQ(X) « 0. The main reason is due to the importance of having a good approximation for UQ[X) W 0, as the solution U(x) = 0 for F E , is unstable. 2.4 Discrete Singular Convolution; Whittaker's Sine Interpolation. In this section, we compare the results i n Sections 2.2 and 2.3 with the results obtained with the discrete singular convolution (DSC) method employed by 47 i Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. Zhao and Wei [182], to solve F E . We also provide a detailed analysis of the numerical aspects of their method. This method is based on the generic Cardinal function due to Whittaker [179] and discussed previously by others [132, 140, 156, 169, 170], and defined in terms of the sine function, C k { X ) ~ l(x-x ) ( 2 k ' 3 4 ) A uniform grid of JV +1 points x — xi + hk, is defined for the finite interval [XL,XR\ where the grid spacing is h = y . This Cardinal function satisfies the interpolation requirement C (xj) = Sj , and the second derivative of f(x) is approximated by k l k < k N f%\x)^Y, k\x)f{x ). C k (2.35) The use of E q . (2.35) on a finite interval instead of the infinite interval results in truncation errors unless the function of interest is well localized in the selected interval [132, 156, 169, 170]. From the explicit differentiation of E q . (2.34), the representation of the second derivative operator in this scheme is, D<$ = I (2.36) We distinguish the sine second derivative matrix operator with a tilde overbar. This scheme is similar to the one used by Wei [174] and by other researchers [132, 140, 156, 169, 170]. These methods employ a uniform grid based on the sine cardinal function as the interpolation function, although they are presented from different perspectives. The time dependent solution is then determined with E q . (2.23), where A = 1 and the derivative matrix operator given by E q . (2.36). Although Wei [175, 176] presents his methodology i n the language of wavelets and signal analysis, it is useful to recognize that it is simply an alternate interpolation scheme analogous to other interpolations schemes such as Lagrange [129, 156, 166] as well as the interpolation defined in terms of orthogonal polynomials in Section 2.2.2. We present an analysis of the numerical aspects of the D S C method with concern to the round off errors associated with the application of the second derivative matrix operator, D ( ) , to UQ(X) and to F E analogous to the analysis for D^ ) in Section 2.3. Zhao and Wei [182] report excellent results for the solution of F E . It is useful to consider a detailed analysis of the second 2 2 48 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. derivative matrix operator that they used analogous to the study in Section 2.2.3 of the Chebyshev-Lobatto spectral method. As the sine function decays slowly as 1/x, Zhao and Wei [182] multiplied the sine function in E q . (2.34) with a Gaussian window function of the form R (x) a = exp ^—j such that Thus, the second derivative matrix operator matrix E)( ) is thus changed to D ' ' with elements, 2 2 f f , ( 38) ~3 + 2^) 3 ~ k where j = As R (x) is a Gaussian and R (0) = 1, only the off-diagonal elements of D^ ) are modified. Zhao and Wei [182] chose the optimum value a = 3.5h. Chen and Shizgal [138, 139] employed a similar weighted interpolation scheme for the solution of Sturm-Liouville problems and the Poisson equation. This was originally proposed by Weideman [178]. a a 2 ~ (2) (2) We show the matrix elements of the fcth row of log (Dj ) and l o g ( i ? ^ j ) in Figs. (2.7A) and (2.7B), respectively, versus j for fixed k =1, 16, 32,'48 and 64, and N = 64. From the results i n Fig. (2.7A), it can be seen that the matrix elements of E)( ) are of the order 10 along the diagonal and decay quickly away from the diagonal. This should be compared with the matrix elements of D ^ whose elements range from 10 on the diagonal to as small as about 1 0 ~ for the extreme off-diagonal elements. The addition of the Gaussian window function has resulted in a rapid decrease of the derivative matrix operator away from the diagonal as expected, and thus the derivative matrix is banded about the diagonal. We consider an additional modification of the differentiation algorithm, E q . (2.35), such that 10 2 k 10 fc 4 4 a 65 f%\x -)* 3 E (2.39) %kf(*k) D where 0 -W j <W j>W _ j W +j m a x ~ \ N j < N- W j>N-W (2.40) 49 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. and W is the number of elements on either side of the diagonal that are retained with 2W + 1 < N. The elements of decay exponentially as we move away from the diagonal due to the introduction of the Gaussian window function in E q . (2.37). From F i g . (2.7B) we see that the terms neglected vary from about 1 0 ~ to 1 0 ~ . This should not influence the final results as these elements of D</) are small. This procedure reduces the bandwidth (RB) from N to 2W + 1 and the derivative algorithm is now more local than global. In the next subsection, we report an analysis of the two second derivative matrix operators, E)( ), D ' ' given by (2.36) and (2.38), respectively, with and without R B . While the polynomial interpolation, E q . (2.20), and the associated differentiation algorithm, E q . (2.22), can give excellent results on a finite interval, the interpolation based on the Sine functions, E q . (2.35), and the differentiation algorithm, E q . (2.39), involve truncation errors especially at the boundaries of the finite interval [XL,XR]. W i t h this in mind, we consider an alternative interpretation of E q . (2.39) whereby the differentiation algorithm is always centred on the point x with W points to the left and W points to the right of the point of interest. Thus 16 67 2 2 2 ff k (2) we define a second derivative matrix operator D which is of dimension (JV + 2W +1) x (JV +1) and operates on solution vectors of F E of dimension (JV + 2fV" + l ) of the form (1,1, 1, U(x ) = U(x ) = 1, U( ),U(x ) = U(x ) = 0,0 0,0) ff 0 L Xl N R with the first W components set equal to unity and the last W components set equal to zero consistent with the boundary conditions. The (JV + 1) components U(xi) are computed from the solution of F E within the computational domain [XL,XR\. 2A.l Analysis of the round-off error for D ^ , D ^ ^ and 2 CT . We consider a solution of F E with the same values used in Section (2.2.3), that is, x = - 0 . 2 , x = 0.8, p = 10 and JV = 64. In F i g . (2.8), we show the computed profile after the first time step for E)() without R B in Fig. (2.8A) and with R B in F i g . (2.8B). In both cases, small amplitude oscillations develop at the foot of the wave front. These oscillations are responsible for the unstable numerical solutions that are obtained. O n the other hand, when the window function R is considered, i.e. the second 2 derivative operator matrix is given by D,^ ), there are no oscillations. As a result, a stable solution is obtained. To understand the origin of the oscillations when the operator I)( ) with and without R B is considered, we study the details in the sum involved in 4 L R 2 a 2 50 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. the calculation of the second derivative of UQ(X), given by JV Sj = Y,£jkU {x ) 0 (2.41) k fc=0 where the matrix operator C can be £)( ) or with or without R B . For the present analysis, two test points x = 0.3 and x = 0.8, which correspond to j = N/2 and j = N respectively, were considered to approximate the second derivative. The elements in the sum i n E q . (2.41), alternate in sign along a row, as shown i n Eqs. (2.36) and (2.38). T h e calculation of the second derivative for j = N/2 and j = iV, is subject to round-off errors analogous to the discussion i n Section 2.3. The values d Uo(xpj/ )/dx — 2.93 x l O " and d U (x )/dx = 2.85 x 1 0 " (N = 64, N/2 = 32) are approximated by a difference of two numbers whose order will depend on the second derivative operator matrix L. A summary of the order of magnitude of the elements, whose difference will approximate d Uo(x3 ) / dx and d Uo(xQ )/dx for different second derivative operators C, is presented i n Table 2.3. For the case of LV ) without R B , both d U (x )/dx = 2.93 x 10~ and d U (x )/dx = 2.85 x 1 0 ~ are approximated by the difference of two numbers of the order of 10 . This round-off error problem is similar to the problems encountered with Chebyshev Lobatto collocation discussed previously. The result of such round-off error is shown i n F i g . (2.8A). Similar instabilities are obtained when E)() with R B is considered. However, as seen i n F i g . (2.8B), there is a notable improvement after the first time step integration for x = x^4- This improvement is a consequence of the approximation of d Uo(x$4)/dx , that is calculated as a difference of two numbers of the order of 10~ (Table 2.3). 2 2 2 2 7 2 2 0 2 2 5 N 2 2 2 2 4 2 2 2 0 7 2 32 2 0 25 64 1 2 2 2 9 (2) On the other hand, when given by (2.38) with and without R B is considered, the solution is free of oscillations and as a consequence, a stable solution is obtained. The results in Table 2.3 show the improvement of the approximation of the second derivative at x = x and x = XQ . It is important to mention that the approximation of d Uo(xQi)/dx — 2.85 x 1 0 ~ , with D^, is numerically inaccurate even though the elements of the sum E q . (2.41) range from 1 0 to 1 0 . The inaccurate approximation is due i n part to the use of the window function R , that modifies the structure of the elements of jLV ). Since the error i n the calculation of d Uo(x64)/dx is at most of the order of 1 0 ~ , the error generated is small for the time integration considered. This leads to an excellent solution for F E with large values of p. ; 4 32 2 2 25 - 2 6 - 2 1 a 2 2 2 21 51 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. In F i g . (2.9A) we show the solution to F E with D^ \ We find a small growing instability near x = XL after several time steps and which after some time ceases to grow. We find that the amplitude of this instability increases with N. The origin of the instability near x = XL in Fig. (2.9A), can be understood in terms of the second derivative operator D„' ' applied to UQ(X). From E q . (2.39), S « O(10 ) is the approximation of d U (x )/dx = 0(10°). This poor approximation is due to the truncation of the physical domain to [XL,XR]. The effect of the truncation on the approximation of d Uo(xk)/dx depends directly on the combination of the values of Uo(xk) and D^ in 2 2 4 2 0 2 0 0 2 2 k E q . (2.39). For x = xr,, the largest value of both Un(xk) and D^\ is for k = 0, and thus the main contribution of the sum, E q . (2.39), is from the first element with magnitude 10 . Therefore, if the vector U denotes the numerical solution at the n time step, we find that U (x2) —> rj ^ 1 as n —> oo, as shown in Fig. (2.9A). In this case, n satisfies the relation k 4 n th n lim D^U \ = pq(l - n) n x=X2 (2.42) n—*oo From E q . (2.42) , it is shown that for any time t = t the poor approximation of d Uo(xo)/dx at x = x\ is always present and So « O(10 ) for t = t . However, when x = XR, the truncation of the domain does not have the same effect as the one shown in F i g . (2.9) for x = XL- The quantity to be approximated is of the order of 1 0 ~ whereas the largest element of E q . (2.39) for x = XR is also of the order of 1 0 . However, the instability that occurs near x = XL, does not destabilize the 2 rest of the solution. This is due to the distribution of the elements of D f f ' ' as shown in F i g . (2.7B), and the magnitude of elements of UQ(X). In order to approximate d Uo(x)/dx near x — XR, the values of Un[x) near x = XL do not play any role in the calculation, even if R B is not considered. This is an important difference between Chebyshev-Lobatto collocation and the weighted D S C (Eq. 2.38) differentiation matrices, applied to F E . Whereas in Chebyshev-Lobatto collocation, the second derivative at some point x = x has a strong dependence on all the Chebyshev-Lobatto points, that is, is global, the D S C with the window function becomes a more local method as only a small number of points play a major role in the calculations. It is important to mention that for a fixed value of p, XL and XR, it is not possible to eliminate the instability at x = XL in Fig. (2.9A). One possible alternative to get a good solution on XL and XR, is to increase the physical domain, keeping only the interval in x where the solution has a good precision. However, this alternative will break the structure given by p, XL and n 2 2 4 n 25 - 2 5 2 2 k 52 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. xn in relation (2.13), giving as a result a different numerical problem. However, when the solution vector is "padded" with W points to the left of XL and W points to the right of XR such that rr/ X f 1 if ~W<k<-\ , n and of dimension (N + 2W + 1) x (N + 1) is used, we obtain the result shown i n Fig. (2.9B) without the instability at x = XL shown in F i g . (2.9A). In this calculation with p = 10 and x G [-0.2,0.8], we verified that both d Uo(x) / dx \ and its numerical approximation So given by E q . (2.41) are of order one 0(1). Again the approximation of d U~o(x)/dx \ = is numerically inaccurate, due to the use of the window function R and because U(x) lies outside the space of bandlimited functions as described in [170]. However, the numerical solution that is obtained is accurate as now n = 1 in E q . (2.42). Therefore, we obtain 4 2 2 x=Xl 2 2 x xi a lim DMlJ \ = n x Xl = pn(l - 77) = 0 (2.44) n—>oo giving as a result what is shown in Fig. (2.9B). Table 2.4 shows a comparison between the results with the sine method with D and the multidomain method presented in Section 3. The multidomain approach provides a more accurate result for the largest value of p. a 2.5 Summary The objective of the present chapter was to develop an accurate and efficient pseudospectral solution of the F E , a prototypical reaction diffusion equation. The collocation method used the Chebyshev-Gauss-Lobatto quadrature points. The solutions of F E are characterized by propagating fronts that can be steep depending on the value of the reaction rate coefficient, p. We compared results for a single domain as well as for a subdivision of the main domain into subintervals. On a single domain the integration of the F E lead to instabilities at some time step. These instabilities were a consequence of the numerical round-off errors arising from the numerical form of the discrete second order derivative matrix operator. We have demonstrated the importance of constructing the differential matrix operator accurately. From a detailed numerical analysis, we have also identified the nature of the round-off errors that occur in the use of the differential matrix operator in 53 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. the numerical solution of F E . This complements the work by Baltensperger and Trummer [128]. The exponentially small values of the solution, U(x,t), combined with the size of the elements in D^) which oscillate i n sign along a row, was the main source of round off error, when a single domain was used. In order to reduce the effects of round off error for the Chebyshev collocation, the main domain was subdivided into smaller subintervals as proposed by Shizgal and coworkers [155, 180, 181]. The multidomain method provided stable and accurate solutions of F E for values of p as large as 10 . We also compared the present numerical treatment with the D S C approach of Zhao and Wei [182] who employed an interpolation based on the sine function. They added a window function R to the sine interpolation function and also limited the number of terms in the differentiation algorithm to a small number about the diagonal elements of the derivative matrix operator. We refer to this procedure as reduced bandwidth R B . We also studied the occurrence of numerical round-off with their method. We found that for fixed values of p, xi and XR, the results with the multidomain method did not present any problem at the left boundary at x = XL as did the D S C method. The instability at the left boundary shown in Fig. (2.9A) obtained with D S C , is attributed to the truncation of the spatial domain as discussed in Section 2.4.1. Another important result of this chapter is the demonstration of the equivalence of the numerical problems denned by Eqs. (2.7-2.9) and Eqs. (2.11, 2.12). However, for Eqs. (2.7-2.9), the parameter p plays an important role and for large p we have the S S W situation, whereas for Eqs. (2.11,2.12) the main parameter is the length of the numerical domain. The results in the present chapter for Fisher's equation, a prototypical reaction diffusion equation, will play an important role in the application of pseudospectral methods to more complex physical systems such as the Fitzhugh-Nagumo equation [134, 144, 161]. 2 6 a 54 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. Table 2.1: Analysis of the multidomain method for p = 10 , At = 10~ , total time of integration T — 0.003; M is the number of subintervals and N is the number of Chebyshev points per subinterval. 4 Scheme Chebyshev Multidomain M(sub) 1 2 6 10 15 25 N points 150 75 25 15 10 6 8 Loo error 0.332(-2) 0.305(-2) 0.456(-3) 0.407(-3) 0.357(-3) 0.724(-2) Table 2.2: Analysis of the multidomain method for p = 10 , At = 1 0 " , T = 0.001; M is the number of subintervals and A'' is the number of Chebyshev points per subinterval. 5 Scheme Chebyshev Multidomain M(sub) 1 10 16 20 25 40 A'' points 400 40 25 20 16 10 8 Loo error NAN 0.132(-3) 0.134(-4) 0.511(-5) 0.773(-4) 0.471(-2) 55 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. Table 2.3: Maximum value of the elements in the sum, E q . (2.41), with C equal to E>() and D£^. When reduced bandwidth (RB) is considered, the index in the sum, E q . (2.41), runs from — W to W, W = 32. 2 Operator C d' U (x ) = 2.93 x 1CT7 Max* ax T)W (No R B ) 10 10 10- 2 0 32 8' Uo(x64) c)xz = 2.85 x 2 Max f e (RB) D^ (No R B ) D^ (RB) 10 1 io- IO" 2) 25 f c 1 1 4 2) 10" io- 4 1 0 9 2 1 -21 Table 2.4: Comparison between multidomain method and the Sine approach with (W=32) for different values of p. T is the total time of integration, M is the number of subintervals and N is the number of Chebyshev points per subinterval. Scheme Sine P 10 10 10 10 10 10 4 5 6 Multidomain 4 5 6 M(sub) 1 1 1 3 20 140 N points 75 400 2800 25 20 20 1X 1X 1X 1X 1X 8X At io-' IO" IO" IO" IO" IO" 8 8 7 8 8 T 0.0032 0.001 0.00033 0.0032 0.001 0.00033 Loo error 4.976(-4) 1.395(-10) 4.966(-8) 4.839(-3) 5.115(-6) 3.125(-12) 56 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. -0.2 0 0.2 0.4 0.6 0.8 -0.2 0 0.2 x C -0.2 . 0 0.2 0.4 0.6 0.8 0.4 0.6 0.8 x 0.4 D 0.6 0.8 -0.2 0 x 0.2 x Figure 2.1: Time dependent profiles versus x for p and N equal to (A) 2 x 10 , 40; (B) 5 x 10 , 64; (C) 10 , 64 and (D)10 , 150. The successive profiles in each graph are for times (A) t = 0.002,0.003,0.004,0.005,0.006; (B) t = 0.0005,0.0015,0.0025,0.0035,0.0045; (C) t = 0.0004,0.0008, 0.0012,0.0016,0.002; (D) t = 0.001,0.0015,0.002,0.0025,0.003. The solid line represents the analytic solution, U(x,t), and the symbols, Vi, the numerical solution evaluated at the Chebyshev points. 3 3 4 4 57 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. Figure 2.2: The relative error versus x in the application of the second derivative matrix operator to the solution at t = 0 for different values of N with (A) p = 10,000 and (B) p = 15,000 Figure 2.3: The relative error RE versus N in the application of the second derivative matrix operator to the solution at t = 0 at the boundaries, x = XR (upper curve) x — XL (lower curve) with p equal to (A) 10000 and (B) 15000. 58 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. 0 4000 8000 12000 16000 20000 0 • 50 100 150 200 j Figure 2.4: (A) Variation of the logarithm of the second derivative versus p, at XR and xi for the analytic solution at t = 0. (B) l o g ( D ) versus j. The elements alternate in sign. Negative terms (o), Positive terms (•) N = 250. 2 1 0 V j Subdomain 1^ X 0 X M-1 X M X 2M-1 Figure 2.5: Partition of the interval [ X ^ J X R ] into Ni subdomains. 59 250 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. 1.2 1 0.8 xT 0.6 0.4 0.2 0 -0.2 -0.2 0 0.4 0.2 0.6 0.8 x Figure 2.6: Time dependent profiles U(x,t) versus x for p = 10 with 140 subdomains and 20 points per domain; The profiles are for times t = 8 x I O , 8 x I O , 1 . 6 x 1 0 , 2 . 4 x 10~ ,3.2 x 1 0 ; the total integration time is T = 3.3 x I O " . 6 - 9 -5 -4 4 - 4 4 60 Chapter 2. A Pseudospectral Method of Solution of Fisher's -0.2 0 0.2 0.4 X 0.6 0.8 -0.2 0 Equation. 0.2 0.4 0.6 X Figure 2.8: Oscillations of U(x,t) after the first time step. D^) is given by (2.36). (A) Without reduced bandwidth ( R B ) , (B) W i t h R B . p = 10 , TV = 64. 2 4 61 0.8 Chapter 2. A Pseudospectral Method of Solution of Fisher's Equation. -0.2 0 0.2 0.4 0.6 0.8 -0.2 X 0 0.2 0.4 X Figure 2.9: (A) Instability of U(x,t) at x = x , p = 10 , At = 1 x I O . The (N + l) x (N + l) second derivative operator considered is given by the algorithm, Eqs. (38) and (39); N = 64; total integration time T = 0.004. (B) Time dependent profiles versus x for p = 10 and times t = 0.001, 0.0015, 0.002, 0.0025 and 0.003 with the (N + 2W + 1) x (N + 1) second derivative 4 - 7 L 4 operator D£°; N = 64, W = 32. 62 0.6 0.8 Bibliography [124] M . J . Ablowitz, A . Zepetella, Explicit solutions of Fisher's equation for a special wave speed, Bull, of Math. Bio. 41 (1979) 835-840. [125] K . Al-Khaled, Numerical study of Fisher's reaction-diffusion equation by the sine collocation method, J . Comput. Appl. M a t h . 137 (2001) 245-255. [126] R . Baltensperger, J . P. Berrut, The errors in calculating the pseudospectral differentiation matrices for Chebyshev-Gauss-Lobatto points, Comp. and Math. Appl. 37 (1999) 41-48. [127] R . Baltensperger, Improving the accuracy of the matrix differentiation method for arbitrary collocation points, Applied N u m . M a t h . 33 (2000) 143-149. [128] R . Baltensperger, M . R. Trammer, Spectral differencing with a twist, S I A M J . Sci. Comput. 24 (2003) 1465-1487. [129] D . Baye, M . Hesse, M . Vincke, The unexplained accuracy of Lagrange mesh-method, Phys. Rev. E 65 (2002) 026701 1-8. [130] G . Beckett, J. A . Mackenzie, A . Ramage, D . M . Sloan, On the numerical solution of one-dimensional P D E s using adaptive methods based on equidistribution, J . Computat. Phys. 167 (2001) 372-292. [131] G . Beckett, J . A . Mackenzie, A . Ramage, D . M . Sloan, Computational solution of two-dimensional unsteady P D E s using moving mesh methods, J . Computat. Phys. 182 (2002) 478-495. [132] J . P. Boyd, Chebyshev and Fourier Spectral Methods (Dover, New York, 2000). [133] P. K . Brazhnik, J . J . Tyson, On travelling wave solutions of Fisher's equation in two spatial dimensions, S I A M J . A p p l . Math. 60 (1999) 371-391. 63 Bibliography [134] A . Bueno Orovio, V . Perez Garcia, F . Fenton, Spectral methods for partial differential equations in irregular domains: The spectral smoothed boundary method, S I A M J . Sci. Comp. 28 (2006) 886-900. [135] C . Canuto, M . Y . Hussaini, A . Quarteroni and T . A . Zang, Spectral Methods in F l u i d Dynamics, Springer Series in Computational Physics, Springer, New York, 1988. [136] G . F . Carey, Y . Shen, Least-squares finite element approximation of Fisher's reaction-diffusion equation, Numer. Methods Partial Differential Equations, 11 (1995) 175-186. [137] J . Canosa, On a nonlinear diffusion equation describing population growth, I B M J . Res. Develop. 17 (1973) 307-313. [138] H . Chen, Y . Su, B . D . Shizgal, A direct spectral collocation Poisson solver in polar and cylindrical coordinates, J . Computat. Phys. 160 (2000) 453-469. [139] H . Chen, B . D . Shizgal, A spectral solution of the Sturm-Liouville equation; comparison of classical and nonclassical polynomials, J . Computat. Phys. 136 (2001) 17-35. [140] D . T . Colbert, W . H . Miller, A novel discrete variable representation for quantum mechanical reactive scattering via the S-matrix K o h n method, J . Chem. Phys. 196 (1992) 1982-1991. [141] P . J . Diamessis, J . A . Domaradzki, J . S. Hesthaven, A spectral multidomain penalty method model for the simulation of high Reynolds number localized incompressible stratified turbulence, J . Computat. Phys. 202 (2005) 298-322. [142] W . S. Don, D . Gottlieb, J . - H Jung, A multidomain spectral method for supersonic reactive flows, J . Comput. Phys. 192 (2003) 325-354. [143] R . A . Fisher, The wave of advance of advantageous genes, A n n . E u genics 7 (1937) 355-369. [144] R . Fitzhugh, Impulses and physiological states in theoretical models of nerve membrane, Biophys. J . 1 (1961) 445-466 [145] J . Gazdag, J . Canosa, Numerical solution of Fisher's equation, J . A p p l . Prob. 11 (1974) 445-457. 64 Bibliography [146] D . Gottlieb, J . S. Hesthaven, Spectral methods for hyperbolic problems, J . Comput. Appl. Math., 128 (2001) 83-131. [147] S. A . Gourley, Travelling front solutions of a nonlocal Fisher equation, J . M a t h . Biol. 41 (2000) 272-284. [148] P. Hagan, Travelling wave and multiple travelling wave solutions of parabolic equations, S I A M J . Math. Anal. 13 (1982) 717-738. [149] T. Hagstrom, H . B . Keller, The numerical calculation of travelling wave solutions of nonlinear parabolic equations, S I A M J . Sci. Stat. Comput. 7 (1986) 978-988. [150] W . Huang, Y . Ren, R. D . Russell, Moving mesh partial differential equations ( M M P D E s ) based on the equidistribution principle, S I A M J . Numer. Anal. 31 (1994) 709-730. [151] W . Huang, R. D . Russell, A moving collocation method for solving time dependent partial differential equations, A p p l . N u m . M a t h 20 (1996) 101-116. [152] A . Kolmogorov, I. Petrovsky, N.Piscounov, Etude de Pequation de la diffusion avec croisance de la quantite de matiere et son application a un probleme biologique, Bull. Univ. Etat. Moscu, A l (1937) 1-25. [153] D . A . Larson, Transient bounds and time-asymptotic behavior of solutions to nonlinear equations of Fisher type, S I A M J . Appl. Math., 34 (1978) 93-103. [154] S. L i , L . Petzold, Y . Ren, Stability of moving mesh systems of partial differential equations, S I A M J . Sci. Comput. 20 (1998) 719-738. [155] G . Mansell, W . Merryfield, B . Shizgal, U . Weinert, A comparison of differential quadrature methods for the solution of partial differential equations, Comp. Meth. Appl. Mech. Engrg. 104 (1993) 295-316. [156] E . H . W . Meijering, W . J . Niessen, Viergeveer, The Sine approximating kernels of classical polynomial interpolation, I E E E Conf. on Image Proc. 3 (1999) 652-656. [157] R. E . Mickens, A best finite-difference scheme for Fisher's equation, Numer. Meth. Part. D . E . 10 (1994) 581-585. 65 Bibliography [158] R. E . Mickens, Relation between the time and space step-sizes in nonstandard finite difference schemes for the Fisher equation, Numer. Meth. Part. D . E . 13 (1997) 51-55. [159] L . S. Mulholland, Y . Qiu, D . M . Sloan, Solution of evolutionary partial differential equations using adaptive finite differences with pseudospectral post-processing, J. Comput. Applied M a t h . 131 (1997) 280-298. [160] J . D . Murray, Mathematical Biology I: A n introduction, Interdisciplinary Applied Mathematics, 19 (2001) 437-459. [161] J . Nagumo, S. Arimoto, S. Yoshizawa, A n active pulse transmision line simulating nerve axon, Proc I R E , 50 (1962) 2061-2070. [162] S. A . Nielsen, J . S. Heasthaven, A multidomain Chebyshev collocation method for predicting ultrasonic field parameters in complex material geometries, Ultrasonic 40 (2002) 177-180. [163] N . Parekh, S. Puri, A new numerical scheme for the Fisher equation, J. Phys. A : M a t h . Gen. 23 (1990) L1085-L1091. [164] R. Peyret, Spectral Methods (Springer, New York, 2002). for Incompressible Viscous Flow [165] Y . Qiu, D . M . Sloan, Numerical solution of Fisher's equation using a moving mesh method, J . Comput. Phys. 146 (1998) 726-746. [166] R. E . Robson, K . F . Ness, G . . E . Sneddon, L . A . Viehland, Comment on the discrete ordinate method in the kinetic-theory of gases , J . Computat. Phys. 92 (1991) 213-229. [167] Rizwan-Uddin, Comparision of the nodal integral method and non standard finite-difference schemes for the Fisher equation, S I A M . J . Sci. Comput. 22 (2001) 1926-1942. [168] J . Roessler, H . Hiissner, Numerical solution of the 1 + 2 dimensional Fisher's equation by finite elements and the Galerkin method, M a t h . Comput. Modelling, 25 (1997) 57-67. [169] C . Schwartz, High accuracy approximation techniques for analytic functions, J . Math. Phys. 26 (1985) 411-415. [170] C . E . Shannon, Communication in the presence of noise (Classic paper), Proc. of the I E E E , 86 (1998) 447-457 66 Bibliography [171] F . Stenger, Numerical methods based on sine and analytic functions, Springer Series in Comp. Math., 20 (1993) 91-96. [172] S. Tang, R. O. Weber, Numerical study of Fisher's equation by a Petrov Galerkin finite element method, J . Aust. Math. Soc. B 33 (1991) 27-38. [173] E . H . Twizell, Y . Wang, W . G . Price, Chaos free numerical solutions of reaction-diffusion equations, Proc. R. Soc. Lond. A . 430 (1990) 541-576. [174] G . W . Wei, Solving quantum eigenvalue problems by discrete singular convolution, J . Phys. B : A t . M o l . Opt. Phys. 33 (2000) 343-352. [175] G . W . Wei, A unified approach for the solution of the Fokker-Planck equation, J . Phys. A : Math. Gen. 33 (2000) 4935-4953. [176] G . W . Wei, A new algorithm for solving some mechanical problems, Comput. Methods Appl. Mech. Engrg. 190 (2001) 2017-2030. [177] J . A . C. Weideman, S. C . Reddy, A M A T L A B matrix differentiation suite, A C M Transactions on Mathematical Software, 26 (2000) 465-519. [178] J . A . C . Weideman, Spectral methods based on nonclassical orthogonal polynomials, Int. Series of Numer. Math. 131 (1999) 238-251. [179] E . T . Whittaker, On the functions which are represented by the expansions of the interpolation theory, Proc. Roy. Soc. Edinburgh 35 (1915) 181-194. [180] H . Yang, B . Shizgal. Chebyshev pseudospectral multidomain technique for viscous flow calculation, Comput. Methods A p p l . Mech. Eng. 118 (1994) 47-61. [181] H . Yang, B . Seymour, B . Shizgal, A Chebyshev pseudospectral multidomain method for steady flow past a cylinder, up to Re=150, Comput. fluids 23 (1994) 829-851. [182] S. Zhao and G . W . Wei, Comparison of the discrete singular convolution and three other numerical schemes for solving Fisher's equation, S I A M J . Sci. Comput. 25 (2003) 127-147. [183] X . Zou, Delay induced traveling wave fronts in reaction diffusion equations of K P P - F i s h e r type, J . Comput. Appl. M a t h . 146 (2002) 309-321. 67 Chapter 3 Pseudospectral method of solution of the Fitzhugh-Nagumo equations. Reaction diffusion systems have been a very active area of research for many decades. They arise in areas such as population dynamics and epidemiology [194], physiology [215] and biology [224]. Classic examples are the modelling of animal coat patterns [224], the Belousov-Zhabotinsky reaction [202, 218, 224, 239], the Hodgkin-Huxley model of the propagation of the action potential along nerve cells [215], and models of the propagation of a disease in an ecosystem [194]. The general form of a reaction diffusion system is given by ^ = D - V 2 V + F ( V ) (3.1) where V = (Vi, V ,Vjv) with Vi = V ; ( x , f ) , (i — 1,JV) can represent the concentrations of N chemical species in a chemical reaction [202], or the number of susceptible, infected and recovered individuals at position x and time t for a model in epidemiology [194]. The function F is called the reaction term and models the local dynamics due to the interactions of the V£. The spatial variation of V j ( x , t) is modelled with the diffusion term D • V V , where the matrix D is the diffusion coefficient matrix. 2 2 The present chapter is devoted to reaction diffusion equations which model the propagation of action potentials in cardiac muscle cells, analogous to the Hodgkin-Huxley model [215]. Fitzhugh [205] and Nagumo et. al. [225] simplified the local dynamics of the Hodgkin-Huxley model. The set of four ordinary differential equations to describe the change in the potential across the membrane of a nerve cell in the giant axon of the squid of the Hodgkin-Huxley model is reduced to a two dimensional system of O D E s called the Fitzhugh-Nagumo ( F H N ) equations. The F N H equations have 3 A version of this chapter has. been submitted for publication. D. Olmos, B . Shizgal, Pseudospectral method of solution of the Fitzhugh-Nagumo equations. J. Comp. Phys. 3 68 Chapter 3. Pseudospectral method of solution ... been applied to numerous other problems for over four decades [187, 190, 212, 214, 215, 227, 239]. For example, the F H N equations are employed to describe the C O oxidation on Pt(110) [187], the study of C a waves on Xenopus oocytes [215] and Medaka eggs [224], and the study of reentry in heart tissue [203, 217]. + 2 The F H N equations preserve the most important feature exhibited by the Hodgkin-Huxley equations, that is excitability. The F H N equations [191, 214, 240] with diffusion can be written as du f = = D V u + \f(u,v) D S7 v + g(u,v) 2 x 2 ( 3 2 ' 2 ) where from E q . (3.1), V = (u,v) , F = (\f(u,v),g(u,v)) and D is the diagonal matrix with elements D\,D , known as the diffusion coefficients for u and v, respectively. The symbol T denotes the transpose of a vector or a matrix. T T 2 In the present chapter, two different examples of F H N equations in the form of E q . (3.2) are considered. The first kinetic model (I) was proposed by Barkley [190] for which fi{u,v) gi(u, v) = = u(l - u) (u - u ) u —v / th 3 N 3 where u h = 2 ^ and a and b are dimensionless parameters. From E q . (3.3), it is clear that the nullclines for the spatially homogeneous case, fi(u, v) = 0, and g\(u,v) = 0 are u = 0, u = 1, u = u h(v) and u = v. A qualitative description of the local dynamics was provided by Barkley [189, 190]. The second kinetic model studied (II) is the classic cubic F H N local dynamics [215], f {u,v) = Au(l - u)(u - a) - v , g (u, v) = u — bv t t 2 3 4 2 where A, a and b are dimensionless constants. The nullclines f (u,v) = 0 and g (u,v) = 0 are v = Au(\ — u)(u — a) and u = bv (see Fig. 4.15 in [215]). 2 2 The phase portraits for these two kinetic models are shown for particular parameter values in Figs. 3.1A and 3.1C. The explicit time variations of u(t) and v(t), where ^ = ^f(u,v) and | ^ = g(u,v), are shown in Figs. 3.IB and 3.ID, respectively. The segments a-b and c-d in all figures occur 69 Chapter 3. Pseudospectral method of solution on a fast time scale determined by e whereas the segments b-c and d-e occur on a slow time scale. It is useful to notice that u for kinetic model II goes negative whereas for kinetic model I, it is flat and close to zero. Also, the segments b-c are different. The point (0,0) on the phase portraits is a stable excitable fixed point. For the particular initial condition, v — 0 and u > u h, u has a large excursion (a-e) as shown in Figs. 3.1A and 3.1C. The quantity u h is referred as the excitation threshold value. When u and v are located at the point d in F i g . 3.1, it is not possible to generate a large excursion when u is taken above any u h < 1, the upper bound of u. In this case, u is said to be in refractory state. The variable u is generally referred to as the excitation (activator) variable and v is the recovery (inhibitory) variable. The fast and slow changes in the F H N equations are characterized by the parameter e [205]. t t t The present chapter has two main objectives. The first is to continue the ongoing effort to develop pseudospectral methods for an accurate and efficient solution of reaction diffusion equations. This chapter follows on the previous work dealing with Fisher's equation in Chapter 2, for which the solutions develop propagating fronts with short spatial variation. The second objective is to compare the solutions obtained with pseudospectral methods with those obtained by Barkley [190] for E q . (3.2) with (3.3). There have been many studies of reaction diffusion systems that are based on kinetic model I proposed by Barkley [185, 187, 209, 220, 222, 223, 242]. Pseudospectral methods are generally considered useful for solving smooth problems and to provide exponential convergence of the solution with respect to the number of collocation points used [195, 231]. However, it has been demonstrated recently that pseudospectral methods can provide a significant improvement over finite difference methods for non smooth problems that develop shocks and steep fronts [207, 226], features shared by the solutions of the F H N equations. Pseudospectral methods provide solutions to partial differential equations defined on a grid of collocation points which are the quadrature points for a given polynomial basis set [193, 195, 204, 208, 231]. The main advantage of pseudospectral (and spectral) methods is the "exponential" convergence of the solutions versus the number of collocation points. There have been only a few studies of the applicability of pseudospectral methods to the numerical solution of reaction diffusion equations [188, 201, 211, 213]. Another aspect of reactive diffusion equations of which the F H N equations are a subset, is that the time scale for the u variable is short compared to the time scale for the v variable. The discretized equations that result from a pseudospectral method (Section 3.1) are stiff and generally a very small 70 Chapter 3. Pseudospectral method of solution time step is required for their integration. Barkley [190] developed a finite difference method for the solution of kinetic model I, which made use of the particular feature that u{t) is very small during the segment d-e. He developed a computer code that sets u(x, t) to zero if it is less than some particular small value and thus, avoids the computation of u(x, t) .during this time interval. For the spatially inhomogeneous F H N equations with nonzero diffusion coefficients, the solutions w(x, t) and v(x,t), subject to no flux, periodic or Dirichlet boundary conditions, are dependent on both position and time. In order to solve the spatial dynamics of F H N equations in two dimensions, different methods have been considered. K a r m a [212] and Mitkov et al [223] used a Fourier pseudospectral method; X i e [241] and Amdjadi [184] used an explicit finite difference scheme; Gottwald [209] and Diks [200] used a method based on finite differences reported by Barkley [190]; Jones and O'Brien [211] solved the Gray-Scott equations with a Fourier pseudospectral method. The F H N and the Gray-Scott equations belong to the family of excitable systems [230] and exhibit the problem of multiple time and spatial scales. However, the reaction term in E q . (3.1) for the Gray-Scott equations leads to a different class of physical problems and will not be considered in this work. As a result of adding diffusion, pulses of excitation propagate through the medium giving rise to different kinds of spatiotemporal behaviors, such as plane waves, target patterns, spiral waves and scroll waves [215, 238]. The algorithm tested in the present chapter will have applications to the study of the propagation of waves in excitable media including the study of spiral waves, [189, 212, 213, 216, 229, 240], the analysis of spiral tip trajectories and their transition to complex patterns [189, 212, 240]. Also, spiral breakup [210, 228], drift of a spiral tip due to the interaction between spiral waves or via wave trains [209, 241], and competition between spirals and periodic circular pacemakers [219], has been studied with F H N type equations. Studies in one dimension [186, 199, 206] and extension to spirals in three dimensions, called scroll waves [198, 222, 238], have been considered. Additional studies have been performed with variants of E q . (3.2), such as the generation of spiral waves due to obstacles [229] and spirals in a medium with random obstacles [237]. In particular, kinetic model I, has been used to simulate the drift of spirals under periodic stimulation [209], to model chemical reactions such as the isothermal C O oxidation on Pt(110) by Bar et al[187], and the study of the effect of noise on spiral dynamics [220, 242]. 71 Chapter 3. Pseudospectral method of solution Kinetic model II (the classic cubic F H N model) has been used to understand the effects of curvature on the speed of a two dimensional propagating front [214], as well as for the comprehension of cardiac arrhythmias [206]. However, despite all the numerical experiments realized with the F H N equations to date, it remains important to develop accurate numerical methods for the solution of E q . (3.2). It is well known [190, 196] that an accurate numerical solution of the dynamics of excitation, which occur on a fast time scale, requires a small time step. Moreover, when considering the reaction diffusion equation associated with the F H N dynamics, E q . (3.2), the problem becomes more complex as the fast transitions observed locally in time are observed to occur in the spatial variables as well. Therefore, for an accurate numerical solution, it is necessary to increase the number of sample points in the spatial grid in regions where fast transitions take place, giving as a result a multiple spatial scales problem as well. Therefore, the problem of obtaining a reliable solution when solving systems such as the F H N equations involving different time and space scales represents a challenging endeavor. This chapter is organized as follows. In Section 3.1, we provide the details of the various numerical methods for the solution of the F H N equations that are compared. The first method is the Chebyshev-multidomain ( C M D ) , based on Chebyshev polynomials, which was shown to be useful in the study of fronts in one dimension [226]. Also, the Fourier pseudospectral (Fourier) method is considered. In Section 3, we compare the solutions for kinetic model I obtained with pseudospectral methods with those with finite difference techniques and the method presented by Barkley [190]. We validate their efficiency in the solution of the F H N equations and the description of spiral waves. We also generalize the Fourier method to the non periodic case which has not been considered previously. In Section 4, we present analogous results for kinetic model II. The use of an operator splitting method for the time integration to decrease the computational time is discussed in Section 5. A summary of the results is presented in Section 6. 3.1 Numerical methods. To solve E q . (3.2) numerically, we consider a two dimensional truncated domain given by Q, — [XL,XR] X [VL^VR]- For the present chapter, D = 0 2 72 Chapter 3. Pseudospectral method of solution as considered by Krinsky and Pumir [217], Starobin and Starmer [235] and Xie et. al. [241] as our purpose is the modelling of propagation of action potentials through excitable tissue for which v represents the non-diffusive gating variable. We consider the system of differential equations for u and v given by f = D^u + \f{u,v) I = 9(u,v) t - (3 5) The initial condition is given by f u(x,y,0) \ v(x,y,0) = u (x,y) = v (x,y) . 0 \ iy) ^ x 0 . W and no flux boundary conditions or Dirichlet type conditions are employed as discussed in Section 3.1.4. 3.1.1 Chebyshev Pseudospectral Multidomain Method. We present in this section the basis for the pseudospectral method of solution of the F H N equation based on Chebyshev polynomials. The Chebyshev polynomials T (z) are orthogonal with respect the weight function w(z) = (1 — z ) / on the interval [—1,1], that is, k 2 - 1 2 - f w(z)T (z)Ti(z)dz = \ix5 c 7-i 2 1 k kjl (3.7) k where c = 1 for all k except for CQ = 2. The Lobatto quadrature points and weights associated with the Chebyshev-Gauss polynomials are given by Xi = — cos (^) and the weights are Wi = jj for all i except = wpf = igfl [193, 195, 231]. These points and weights provide the approximate quadrature, 1 N k / w(z)f(z)dz~Y if( i) Ji=0 w (-) z 3 t 1 8 where N is the number of points. Since any piecewise continuous function, / € Z^O, 1] can be expanded in a Chebyshev polynomial series that is convergent in the mean of the norm, we have that JV f(z)*f {z) N = Y *kT {z) J k (3.9) fc=o 73 Chapter 3. Pseudospectral method of solution ... where f 2 1 a = / k (3.10) w(z)f(z)T {z)dz k W i t h Eqs. (3.8) to (3.10) we obtain the interpolation algorithm N /AKZ)-^/^)/^) (3.11) where the interpolating polynomials, Ij(z), are given by 2 N fc=0 where vn = = 1/2 with v = 1 if k ^ 0, N and where Ij(zi) = 5ij is the cardinal condition such that //v(-Zj) = f(zi) [193]. The second derivative of f(z) at the quadrature points is then given approximately by k N / # ^ * ) * £ / f (3.13) If we denote by f, the ./V dimensional vector of the function evaluated at the Chebyshev-Lobatto points, E q . (3.13) can be rewritten as f( )=D( )-f 2 (3.14) 2 where the second derivative matrix is given explicitly by ^ ? = ' j 2 , ( ^ ^ l ~ f (3-15) c This is the basis for the Chebyshev pseudospectral method, whereby E q . (3.5) can be reduced to a set of O D E ' s with the representation of the second derivative operator with T)( \ 2 Thus, with the application of the pseudospectral method based on E q . (3.15) to E q . (3.5) i n 2D we obtain, duij dt = A £fl0 D$u ~~ "™ L~ik=0^ik kj "*J + B Eto *ikD$ 1 " £^k=0 "iK^jk + \f(u ,v i3 J\™t]'"VJ 1 e l: (3 16) where u {t) = u(x ,t). In E q . (3.16), A = and B = and they appear as a consequence of the linear transformations [XL, XR] and i3 uyj 74 Chapter 3. Pseudospectral method of solution [VLJVR] to [—1,1] and include the respective diffusion coefficients. No-flux boundary conditions are implemented by solving the system EL><^ = 0 ( 3 - 1 7 ) for UQJ and upfj, j = 1,N — 1 and similar conditions for Uio and Uin, i = 1,N1. In order to apply the Chebyshev pseudospectral method, we employ a multidomain approach used previously [226]. It consists of dividing the intervals [XL,XR\ and \yi,yp\ into N overlapping subintervals, 1^ = [XQ,X^} and I = [yoiyjw]; respectively, with M = N — 1 and u and v = 1,N , and in each dimension all the subintervals have the same length. For each subinterval, we apply the procedure described in Eqs. (3.9) to (3.13) with the resulting system of coupled O D E s given by E q . (3.16) with A = £ ^ s v crl S Dl X and B = y^ ^ l y a n M 0 x d the indexes in Eqs. (3.16) and (3.17) going from 0 to (Nch - 2 ) x J V , + l . The first and second derivative matrices D ^ , D^ ) in Eqs. (3.17) and (3.16), respectively, for the Chebyshev multidomain ( C M D ) method, become a block diagonal matrix as shown on Section 2.3. 2 The application of Chebyshev multidomain in the solution of E q . (3.5), requires a choice of two parameters, the number of subdomains A^ and the number of Chebyshev points per subdomain N h- In order to attain convergence, the number of collocation points has to be increased. For Chebyshev multidomain it is possible to increase both N and N or fixing one while increasing the other. In this work, we focus on increasing N and fixing to a relatively small value. The reason for this is discussed in Section 3.2.1. Increasing the value of N h, requires a decrease in the time step size to preserve stability, thus increasing the computation time. s c cri s s c 3.1.2 The Fourier method. The Fourier method is a standard technique used in fluid dynamics [195], heat conduction [208] and more recently used in reaction diffusion systems [188, 211, 212]. The Fourier method has been used by K a r m a [212] to solve the F H N equations and to study the transition to meandering in a spiral, whereas Jones and O'Brien [211] used a Fourier method to solve the Gray Scott equations. 75 Chapter 3. Pseudospectral method of solution To apply the Fourier method to solve Eq. (3.5) in two dimensions we apply the formalism followed in Canuto et al. [195] (Page 78). We require that Eq. (3.5) is satisfied at the collocation points Xj = y = j, r = 0, TV — 1, r dU du 2 W at dv dt = (3.18) g{u,v)\ where D = D\ (^) , appears from the linear transformation of [—L, L] to [0, 2TT] in each dimension. Eq. (3.18) can be rewritten as ,JV\ dt dv" dt (3.19) The differentiation matrix D^) is given by U (xj,y ,t). where Uj (t) A~ CA. The matrix A is the transformation from physical to frequency domain, C is the diagonal matrix with elements —k , representing the second derivative in the frequency domain, and A~ is the operator matrix that returns the solution to the physical domain. The transformation from physical to frequency domain was computed using the FFT. For a complete discussion of the matrices A and C, refer to page 126 in Canuto et al [195]. 2 N r r l 2 l 3.1.3 Barkley's m e t h o d . The method proposed by Barkley [190, 198] to solve Eq. (3.5), is a modified version of the usual finite difference (FD) scheme and was designed for the family of kinetic models with a form analogous to Eq. (3.3). This method, which has been widely used by several researchers [185, 187, 209, 220, 222, 223, 242], takes advantage of the special form of the nonlinearity given by / in Eq. (3.3). Figure 3.IB shows that u(t) exhibits fast transitions along the segments a-b and c-d, followed by the slow time variation along d-e for which u(t) « 0. It is thus possible to define a quantity <5 << 1 such that for times t in d-e if u ( x , y , t k ) < 5 then it is redefined as u ( x , y , t k ) = 0. In [190], Barkley considered 5 = 0.0001 as a reliable value for his computations. With this technique, Barkley avoids unnecessary computations of the values u [ x , y , t k ) , during the time segment d-e. This technique is not directly applicable to kinetic model II, since u(t) is not small in the segment d-e as shown in Fig. ID. k 76 Chapter 3. Pseudospectral method of solution 3 . 1 . 4 I n i t i a l a n d b o u n d a r y c o n d i t i o n s . The numerical schemes presented in the previous section are employed to study propagating wave solutions of E q . (3.5), specifically the propagation of one dimensional (ID) plane waves and two dimensional (2D) spiral waves. The solutions are studied over the interval 0, = [-L, L] i n I D and over the square domain Q, = [—L,L] x [—L,L] in 2D. The initial conditions considered have the general form u {x,y) = vo(x,y) = c 0 ( l + e x p ( 4 ( | a | - c ) ) ) - - ( l + exp(4(|2 |-c )))2 ; 3 1 ; 2 2 Vy Vx,\/y (3.20) The functions uo(x,y) and vo(x,y), which are also considered to define initial conditions in I D , are redefined over some regions i n order to obtain a propagating pulse in I D and a spiral wave in 2D. The methods considered in this section for solving E q . (3.5), are C M D , F D the Fourier and Barkley methods. For the first three methods, a second order Runge K u t t a algorithm is used to advance the time as it provides much more accurate solutions than would a simple Euler integration. In the present chapter, no flux boundary conditions are considered for F D , the Barkley method and C M D as in [190, 228, 237, 241]. Comparisons between the Fourier method and F D are considered with different boundary conditions. The Fourier method is a special case due to its applicability to periodic problems and requires periodic boundary conditions as reported by K a r m a [212]. In order to extend the dynamics to the non periodic case, it is necessary to remove periodicity i n the solution obtained with the Fourier method by implementing a different type of boundary condition. The boundary conditions implemented are a variation of the Dirichlet type and has been considered previously for kinetic model I [234]. It consists of fixing the value of the numerical solutions of E q . (3.5), u(x, y, t) and v(x, y, t) at the boundary, such that u(x,y,t) is in its refractory state as shown by the position 'd' i n Figs. 3.1A and 3.1C. This constraint will eliminate the passage of a propagating front through the boundary. The effects of the choice of boundary conditions i n two dimensions for the Fourier method, is shown in F i g . 3.2, where the dynamics take place on the squared domain [—30,30] x [—30,30]. Figure 3.2A is the result of the integration of E q . (3.19) which includes diffusion and Dirichlet boundary conditions and shows the contours of the solution u(x,y,t) for some time t. 77 Chapter 3. Pseudospectral method of solution The spiral shown in F i g . 3.2A will continue to rotate and propagate radially outwards and eventually annihilate when it hits the boundaries of the domain, simulating the effects of considering no flux boundary conditions. B y contrast, the behavior in F i g . 3.2B is obtained for periodic boundary conditions. Therefore, comparisons using Dirichlet boundary conditions for the Fourier method are presented for F D only. It is important to notice that with the procedure explained above, it is possible to use the Fourier method to solve any F H N type equation for the non periodic case. 3.2 3.2.1 N u m e r i c a l studies with kinetic model I. Plane I D waves. In order to test the accuracy of the numerical methods, we first solved E q . (3.5) with (3.3) for the I D case. The domain considered is x £ [—30,30]. No flux boundary conditions are applied in conjunction with F D , C M D and the Barkley method, whereas Dirichlet boundary conditions are used for the Fourier method and F D . Figure 3.3 shows successive profiles of the pulse u(x,t), which depends only on x and t, at times t* = 0,1.66,3.33 and 5, for different numerical methods. The parameters considered are e = 0.005, a = 0.3, b = 0.01, D\ = 1 and total integration time t = 5. The profile at t = 0 is the initial condition for u and is given by » M = { O ° i l l W x < -25 x > -25 «-•<»-{ (3.21) with c i = 25 and c = 21 in E q . (3.20). 2 Figure 3.3A shows a comparison of the numerical solutions obtained with the Fourier method (•) and F D (x) with N = 512 points and Dirichlet boundary conditions. Figure 3.3B, compares C M D (•) with N = 40 subdomains and N h = 13 points per subdomain with a total of N = 442 points, and F D ( x ) with N = 512. We show in F i g . 3.3C, a comparison between C M D (•) with N = 442 and the Barkley method (0) with N = 512. s c In Fig. 3.3, it is qualitatively observed that the speed of the pulse obtained with the Fourier method is larger than the one obtained with F D (Fig. 3.3A), whereas the speed of the pulse obtained with C M D is larger than the speed obtained with either F D or the Barkley method (Figs. 3.3B-C). 78 Chapter 3. Pseudospectral method of solution A measure of the global error is given by the rate of convergence of the wave speed. The speed of the pulse is calculated by finding x* such that u{x*,t) = 0.5 at the wavefront. A quantitative analysis of the convergence of the wave speed is presented in Tables 3.1-3.5 for the different numerical schemes. We compare the rate of convergence of all methods and use c = 9.059835 as the converged value of the wave speed as given by C M D and the Fourier methods. In Tables 3.1 to 3.3, the accuracy of the speed c is independent of the value of At, which varies from I O to 1 0 . Large values of N implies a reduction of A t for stability purposes. Table 3.1 shows that with the Fourier method, N = 512, four figure accuracy is obtained. The converged wave speed is obtained with N = 1024. Table 3.2, shows the convergence of the wave speed for C M D with N = 40 and different values of N h- Similarly, Table 3.3 shows the convergence of the wave speed for C M D with N h = 13 and different values of N . The total number of grid points is N = (N h — 2)N +2. Four significant figures are obtained with N = 772 i n both Tables 3.2 and 3.3. Convergence to c is obtained with N = 1922 and N — 1982 in Tables 3.2 and 3.3, respectively. The convergence in Tables 3.2 and 3.3 versus the total number of points, N, appears to be very similar. Reduction of N h increases the sparsity of the second derivative operator D^ ), reducing the number of operations involved i n evaluating the second derivative. Therefore, our computations will consider relatively small values of - 4 - 6 s c c s c s c 2 For F D , the convergence to the speed c, which shows a dependence on A t , is much slower than either C M D or the Fourier method as shown i n Table 3.4. Four figure accuracy is barely achieved with N = 8192 and A t = I O . Similar results were obtained with Dirichlet boundary conditions and not reported. The poor convergence of the wave speed with the Barkley method is shown i n Table 3.5. The value of the wave speed closest to c is 9.0386 with N = 4096, A t = I O " and <5 = I O " . From Table 3.5, it is shown that for a fixed number of points, the accuracy of the speed deteriorates as A t is decreased. When the parameter 5 is reduced from S = 1 0 to 5 — 10~ the accuracy of the wave speed improves. - 7 4 5 - 4 5 A convergence analysis of the numerical solution is presented in F i g . 3.4 for t* = 5, where the local error given by E ( ) = log (|ur (^.**) - i ^ t e . O I ) act N Xi 10 (3-22) 79 Chapter 3. Pseudospectral method of solution is calculated with the different methods. In E q . (3.22), u^ (t*) is the converged numerical solution of E q . (3.5) with (3.3) given with the Fourier method with TV = 8192 for Dirichlet boundary conditions and considered as the "exact" solution. For no flux boundary conditions, this exact solution is obtained with C M D , with TV = 9902 points. A t t* = 5, the pulse is located in the interval x £ [14,30] as shown by the last profile on the extreme right of Fig. 3.3. In F i g . 3.4, Ejy(x) is shown for x £ [14,30] as the major source of error comes from this region. xact In F i g . 3.4A the Fourier method gives for TV = 512, E^(x) « 10~ (dashed curve), and for TV = 4096 (solid bold curve), E (x) w 1 0 " . In F i g . 3.4B, the C M D calculations give EN(X) ta 10~ with TV = 442 (dashed curve) and EN{X) « 1 0 - with TV = 1982 (solid curve). Figure 3.4C shows the convergence of the F D scheme with no-flux boundary conditions. To obtain EN{X) « 1 0 , TV = 8192 points are required (dashed bold curve). Finally, the convergence analysis for the Barkley method is shown in F i g . 3.4D with no-flux boundary conditions. We find that E^(x) « 10°, for values of 7V up to 8192 (solid bold curve). Also, the solution loses accuracy when TV is increased from 2048 (solid thin curve) to 4096 (dashed bold curve). 2 15 N 2 -7 - 2 We thus conclude that the pseudospectral methods C M D and Fourier, outperform those based on finite differences, F D and the Barkley method. In order to corroborate that C M D gives the correct solution, we find that the asymptotic speed of a planar wave is given by [224] (3.23) In F i g . 3.5, we show that the converged value of c with C M D does indeed tend to c as e —> 0. a 3.2.2 Spiral waves. We solve the 2D Fitzhugh Nagumo equations, E q . (3.5), so as to give spiral wave solutions. Thus, E q . (3.5) with (3.3) was solved with e = 0.005, a = 0.3, b = 0.01, Di = 1 as done in [190], over the square domain fl = 80 Chapter 3. Pseudospectral method of solution [—30,30] x [—30,30]. The initial condition considered is u(x,y,0) = | (3.24) v{x,y,0) = | with c i = 5 and-C2 = 1 in E q . (3.20), where (J and P|, refer to the union and intersection of sets, respectively. W i t h this initial condition, spiral waves are generated. Although a spiral wave develops at t = 1 for the initial conditions given by E q . (3.24), we study the spiral at t* = 8 which develops over the whole domain and has completed two rotations. Comparisons of the three methods, C M D , F D and the Barkley method, with no flux boundary conditions are presented in this section, and the results are summarized in Figs. 3.6 to 3.8. The Fourier method gives quite similar results as C M D and will not be considered for the convergence analysis. In Fig. 3.6, we show some contour plots of the numerical solution u(x, y, t*). In F i g . 3.6A, a comparison between F D with Dirichlet boundary conditions and TV = 512 points ( x ) , and the Fourier method with TV = 1024 (•) is shown. The solutions obtained with the Fourier method with TV = 512 and TV = 1024 points, are undistinguishable from each other at this level of comparison. Therefore, for TV = 512, the Fourier method gives a converged solution whereas F D has not converged yet. Figure 3.6B shows that for C M D , the solution with TV = 442 (•) is indistinguishable from the one with TV = 1432 (o). In Figs. 3.6C and 3.6D, the solution obtained with C M D with TV = 1432 (o) is compared with F D ( x ) and the Barkley method (0), both with TV = 512. The solutions given by F D and the Barkley method, converge more slowly than the results with CMD. A n alternate comparison of the different numerical methods is given by the I D profiles obtained from a slice of a spiral at y = —15 (horizontal dashed line in F i g . 3.6B) for a particular t. The slices shown in Figs. 3.7A, 3.7C and 3.7E, correspond to the spiral solutions shown in Figs. 3.6B, 3.6C and 3.6D, respectively. Figure 3.7A confirms that the solution with C M D converges rapidly as the profiles with TV = 442 (•) coincide with the profiles with TV = 1432 (o). The results in Figs. 3.7C and 3.7E show that the numerical solutions with F D and the Barkley method have not converged with 81 Chapter 3. Pseudospectral method of solution N = 512. The departures from the converged solution are significant. The errors E^{xi) given by E q . (3.22), for the I D solutions shown in Figs. 3.7A, 3.7B and 3.7E with different numbers of points and the different numerical schemes, are presented in Figs. 3.7B, 3.7D and 3.7F. In this case, u | ( t * ) in E q . (3.22) is the numerical solution obtained with C M D and N = 1432. xact In F i g . 3.7B, E^r(xi) is shown for C M D . The errors were calculated for N = 332 (dashed), 442 (solid), and 1102 (dashed bold). For F D and the Barkley method in Figs. 3.7D and 3.7F, respectively, the number of points considered were AT = 512 (dashed), 1024 (solid) and 1500 (dashed bold). From F i g . 3.7B, the error obtained with C M D with N = 1102 (dashed bold) is EN(X) « 10~ . The same accuracy cannot be obtained with F D nor the Barkley method with AT up to 1500 (Figs. 3.7D and 3.7F), where EN(X) > 1 0 . The errors obtained with the Barkley method and F D with N = 1024 (solid lines in Figs. 3.7D and 3.7F), are comparable to the error obtained with C M D N = 442 (solid line in F i g 3.7B). In order to illustrate the evolution with time of the accuracy of the solution, we show in F i g . 3.8 spirals for an integration time up to t* = 100. It is clear from F i g . 3.8 that the solutions for C M D with AT = 442 (•) and N = 882 (o) coincide. Thus, the solutions for C M D remain converged for this solution at long times. B y contrast, the solutions for F D ( x ) and the Barkley method (0), both with N = 512 points, differ significantly and also in comparison with the solution obtained with C M D with N = 882. 2 _ 1 A n important feature of a spiral wave is given by the trajectory of its tip as discussed in [191, 240]. From the results shown in F i g . 3.8, it is observed that the tip of the spirals obtained with F D and the Barkley method, do not coincide with the spiral tip given by the C M D method. This implies that the trajectories of the tip might not agree with each other. Based on the numerical comparisons presented, we conclude that the results obtained with the C M D method are the most reliable. The computational time ( C P U ) for the different numerical schemes was determined versus N and At. The C P U was found to vary in accordance with C P U = + B. For kinetic model I, the values of the coefficients in this relationship are A = 4.69, and B = -1.97 for C M D , and A = 3.48, and B = 1.61 for F D , respectively. Thus, C M D with N = 442 is approximately 5 times faster than F D and twice as fast as the Barkley method 82 Chapter 3. Pseudospectral method of solution with TV = 1024, with all three schemes providing the same accuracy. The variation of the C P U time versus ^ for the Barkley method lies below the line giving C P U versus ^ for the F D method. The C P U time decreases for the Barkley method with increasing S as expected. Angular speed. For particular choices of the parameters a, b and e, for kinetic model I, we obtain different trajectories of the spiral of the tip. The position of the tip is defined as i n [189] by Barkley, that is, the intersection of the contours w = 5 and f(u = 2, v) = 0. The trajectory of the tip depends on the local kinetics, as shown by Barkley [189, 191]. It goes from simple patterns, such as a circle, to more complicated flower like patterns [189, 191]. The trajectory of the spiral tip is referred to as spiral meandering [191] and is considered as an important aspect of cardiac arrhythmias. A n accurate description of this phenomenon is thus important. To calculate the angular speed, we have to restrict the parameters a, b and e, such that the tip of the spiral is a circle for which Barkley [191] used a = 0.52, b = 0.05 and e = 0.02. However, in this chapter we considered a square domain, different from the circular domain in [191], and our simulations employed e = 0.0208, so as to give a good approximation to a circular trajectory. Equation (3.5) with (3.3) was solved over the rectangular domain [—7.5,7.5] x [—7.5, 7.5] up to a time t* = 40. The initial condition is given by otherwise (3.25) otherwise with c i = 3 and c = 1 in E q . (3.20). 2 We measure the position of the tip for the interval t G [20,40], with one unit of time between measurements. The data obtained was fitted to a circle. The results obtained for the rotational speed versus N and At for the different methods appear to converge faster than the wave speed in the I D situation, although to much fewer significant figures. In each case, we obtained for each method very similar results for pairs of values of iV 83 Chapter 3. Pseudospectral method of solution and A t . For example, for C M D at t* = 30 we get for N and A t equal to (482,2.5 x 10~ ), (642, IO" ) and (884, I O " ) , rotational speeds of 0.7608, 0.7578 and 0.7614, respectively. For t* = 40, the rotational speeds are 0.7460, 0.7474 and 0.7467. For F D and the Barkley method, are used N and A t equal to (512, I O ) and (1024,5 x I O ) . The rotational speeds for F D were 0.7527 and 0.7582 at t* = 30 and 0.7538 and 0.7449 at t* = 40, respectively. For the Barkley method, the rotational speeds were 0.7818 and 0.7840 for t* = 30 and 0.7625 and 0.7677 for t* = 40 respectively. Although this comparison is useful, our results do not permit a definitive statement regarding the convergence of the different numerical methods. This could be due to the fact that the trajectory of the tip is not exactly circular. 5 5 - 4 3.3 5 - 5 N u m e r i c a l studies for kinetic model II. In this Section, we present a numerical study of the classical cubic F H N model, E q . (3.4), with diffusion, which has been used by several authors [206, 212, 214, 227]. The parameters considered are e = 0.001, a = 0.1, b = 0.5 and A = 1. The numerical methods considered in the present section are F D and C M D . For kinetic model II, the Barkley method [190], cannot be applied for this model as stated in [198] because for kinetic model II the refractory period of u(t) (segment d-e in Fig. 3.IB) has an overshoot and setting u(t) = 0 if u(t) < S with 5 « 1 does not apply. 3.3.1 Plane waves. Analogous to the study of I D waves for kinetic model I, we analyze the propagation of a I D pulse obtained with the initial condition, x>0 x <-17 otherwise otherwise (3.26) with c i = 17 and c = 13 in E q . (3.20). In F i g . 3.9, we show the profiles from left to right obtained for t* = 0,0.6,1.2 and 1.8 respectively. 2 In Figure 3.9A, the agreement of the results with C M D for N = 552(») and N = 8802(o), is a demonstration of the convergence of the solution. We take C M D with N = 8802 as the "exact" converged solution. In F i g . 3.9B, the results with F D with = 512 (x) are compared with those attained with 84 Chapter 3. Pseudospectral method of solution C M D with N = 8802(o). It is clear that F D gives a very poor solution, as the pulse travels at much lower speed than the one obtained with the C M D method. W i t h N increased to 1024, the accuracy of the solution obtained with F D ((x) in F i g . 3.9C), improves noticeably, but the solution differs from the solution with C M D and TV = 552. A quantitative analysis of the speed of the pulse is shown in Tables 3.6 and 3.7, for C M D and F D , respectively. For both methods, the convergence of the speed c is independent of the value of A i , provided a stable solution is obtained. The results in Table 3.6, confirm the rapid convergence of the results with C M D . We get the wavefront speed c to five significant figures with C M D for N = 1322, and to eight significant figures for N = 2752. B y contrast, in Table 3.7 we show only four significant figures for the speed with F D for N = 8192. A study of the error given by E q . (3.22) is shown in F i g . 3.10. The last profile in F i g . 3.9, obtained with C M D and N = 8802 is considered as u | ( i * ) (t* = 1.8). The variation of the error, EN(X), is shown versus TV for C M D and F D in Figs. 3.10A and 3.10B, respectively. The convergence with C M D with N = 552 (dashed), TV = 1322 (solid) and N = 3852 (solid bold) is clearly much faster than the convergence with F D with N = 512 (dashed), N = 1024 (solid) and N = 8192 (solid bold). W i t h these choices for the number of grid points the error with C M D can be made as small as 10~ over the whole interval whereas with F D it is only 10~ within the interval of interest. TOCt 8 2 To validate the numerical values of c in Table 3.6, we show the variation of c / c versus e in F i g . 3.11 where the asymptotic speed c is determined as done in [224] and given by a a (3.27) 3.3.2 Spiral waves. We also consider the simulation of spiral waves in 2D with the different numerical methods for kinetic model II. A spiral wave was initiated by con- 85 Chapter 3. Pseudospectral method of solution sidering the initial conditions u(x,y,0) = v(x,y,0) = 0 u 0.15 0 0 {x < 0} \J{y > 5} otherwise {a;<l}n{y<-10} otherwise (3.28) with c i = 5 and c = 1 in E q . (3.20). The parameters considered are e = 0.001, a = 0.1, b = 0.5 and A = l over the domain f i = [-20, 20] x [-20, 20]. 2 In F i g . 3.12, we present an analysis of the numerical accuracy of a spiral solution obtained with C M D and F D at time t* = 5, analogous to Figs 3.6 and 3.7 for kinetic model I. A comparison of contour plots (u(x, y, t*) = 0.5) with the two different numerical methods is shown in Fig. 3.12A. The contour solutions in Fig. 3.12A, were obtained with C M D with N — 1102 (o), N = 662 (•) and F D with N = 1024 (dashed). It follows that the C M D solution is converged with N = 662, whereas the solution with F D and N = 1024 is not converged. For the same solutions, a cut of the spiral at y = —10 and x £ [—20,20] (dashed straight line in F i g . 3.12A) is shown in Fig. 3.12B. In Fig. 3.12B, it is confirmed that the solutions obtained with C M D agree, whereas the solution obtained with F D with N — 1024 does not coincide with the C M D solution with N = 1102. In Figs. 3.12 C and D, the error E^{x) versus x, for C M D and F D , respectively, at t* = 5 and y = —15 is shown for different values of N. In this case, u1 (t*) is the numerical solution given by C M D with N = 1102. In F i g . 3.12C, the error with C M D and N = 662 (solid) is of the same order as the error given by F D with N — 1500 (dashed bold) shown in Fig. 3.12D. xact In the same way as done for kinetic model I, the C P U time satisfies the relationship, C P U = + P>, with A = 7.49, and B = -6.01 for C M D and A = 5.18 and B = 1.154 for F D , respectively. In this case, C M D with Af = 662 is two times faster than F D with A^ = 1500. A more accurate solution has been achieved with C M D than F D with half the number of points and faster by a factor of 2. Angular speed for kinetic model II. We conclude the studies of convergence, by measuring the angular speed of the wave given by E q . (3.5) with kinetic model II (Eq. 3.4) as done in 86 Chapter 3. Pseudospectral method of solution Section 3.2.2. We considered parameters in E q . (3.4), such that the spiral tip trajectory is as close as a circle as possible. In this case the spiral tip was defined by taking the intersection of the two level curves, u = 0.5 and v = v /2, where v is the local maximum of f2(u,v) = 0 (Eq. 3.4), i.e. v = Au(l — u)(u — a). The initial condition is max max otherwise (3.29) otherwise with c i = 5 and c = 1 in E q . (3.20). The parameters considered in E q . (3.4) are e = 0.097, a = 0.098,6 = 0.55, A = 4 and the simulations were carried over the square domain = [—5,5] x [—5,5]. The simulations were carried out for a time t = 30. 2 In this case, the results obtained with C M D give a faster rate of convergence than F D . However, similar to the studies for kinetic model I, the converged values have fewer significant figures than the values obtained for the speed of the I D pulse. As an example, for C M D at t* = 24 we get for the pairs (TV, A t ) equal to (162,10~ ), (322, 10~ ), (482,10" ) and (642, 5 x 1 0 " ) , rotational speeds equal to 1.4824,1.4813,1.4822 and 1.4826, respectively. For t* = 30, the rotational speeds are 1.4910,1.4917,1.4923 and 1.4912. For F D at t* = 24 and pairs of (TV, At) equal to (256,10~ ), (512,5 x 10~ ) and (1024,10 ) we obtain rotational speeds equal to 1.4555,1.4780 and 1.4826, respectively. For t* = 30, the rotational speeds are, 1.4503,1.4666 and 1.4920. 4 5 5 6 4 5 -5 3.4 Operator splitting method for C M D . In this section, we apply the operator splitting method (OS) [232] to solve E q . (3.5) for kinetic model I and II. Particularly, we restrict the use of the OS technique to C M D . The OS method has been proposed as an efficient algorithm to integrate reaction diffusion equations [233], and particularly excitable systems [232, 236]. It consist of solving the diffusion and the reaction part in E q . (3.16) as two uncoupled processes as described i n [232], where a general description of this method is given. The values of A t = 10~ and 10~ for kinetic models I and II, used in the C M D method, were chosen for stability purposes. These values of A t are 4 5 87 Chapter 3. Pseudospectral method of solution chosen by considering the range of the eigenvalue spectrum in the C M D second derivative operator matrix D^ ). W i t h N = 442, the eigenvalues vary from 1 to 10 and this range increases as we increase the number of Chebyshev points per subinterval Nch with N„ fixed. The eigenvalue spectrum of D^) appears to be independent of N . Therefore, the choice of A t for solving E q . (3.5) numerically is restricted by the fast modes in D^ ). 2 3 2 s 2 In order to solve E q . (3.5) for a time step A t with the OS method, we apply the Euler method to the diffusion term for k time steps with each step equal to whereas the reaction part was solved in time with forward Euler and size step A t . The boundary conditions were applied after advancing one time step ^ in the diffusion part. The spiral solution for kinetic model I at t* = 8 is shown in F i g . 3.13A, where the accuracy and efficiency of the OS method is evaluated. In F i g . 3.13A, we show the solutions for C M D without OS (N = 442; solid) and with OS (N = 442 (•) and N = 386 (x)). If we choose N = 442, A t = 0.002 and k = 6, the solution with OS is ten times faster than without OS, but the accuracy is compromised. Comparable accuracy is obtained for OS with N = 442, and OS with N = 386, A t = 0.004 and k = 8. From a comparison of Fig. 3.13A and Fig. 3.6D, the solution obtained with the OS method with N = 386 gives a solution with the same accuracy as the Barkley method and N = 512, but the OS method is twice as fast. To get an idea of the computational time for kinetic model I, we consider the OS method with N — 386. The computational time to simulate one complete rotation of the spiral was 45 s, corresponding to 2.4 time units in the model. The results for kinetic model II are shown in Fig. 3.13B, where we compare the numerical solutions for C M D without (N = 662; solid), and with OS (AT = 662; (•) and N = 602 (x)). The numerical solution with OS and N = 662, A t = 2 x I O and k = 4, gives a very good approximation to the converged C M D solution without OS with only a quarter of the computational time. A solution with the same accuracy as the one obtained with OS and N = 662 is given by OS with N = 602, A t = 0.0008 and k = 2. From Figs. 3.13B and 3.12A, the solution obtained with the OS method is somewhat less accurate than F D with N = 1024 (dashed in Fig. 3.12A). However, the computational time obtained with the OS method with N = 602 is up to 30 times faster than F D .with N = 1024. - 4 88 Chapter 3. Pseudospectral method of solution Similarly, an estimate of the computational time for kinetic model II is considered for the OS method with N = 602. The computational time to simulate a complete rotation of the spiral was 130 s, corresponding to an approximate of 0.55 units of time. 3.5 Summary. The objective of the present chapter was to develop accurate and fast algorithms to solve partial differential equations of the F H N type. Such equations present propagation of shock like waves as in [226]. The Chebyshev multidomain algorithm proposed was based on a pseudospectral approach with Chebyshev-Gauss-Lobatto points considered previously by Olmos and Shizgal [226] and Don et al [197]. We corroborate that the Chebyshev multidomain approach is an efficient algorithm for solving equations involving shock like waves as done previously with Fisher's equation [226]. It appears that there has not been a comparable study of the convergence of different numerical methods for reaction diffusion equations of the F H N type. We compared Chebyshev multidomain with finite difference and the method proposed by Barkley [190]. Simulations were carried out for two F H N equations with different local dynamics. The first equation, named kinetic model I, was proposed by Barkley [190], whereas the second equation, kinetic model II, is the classic cubic F H N model discussed by Keener [215]. We have demonstrated the superiority of the Chebyshev multidomain approach with regard to accuracy and computational time in comparison with finite difference and the Barkley method as discussed in Sections 3.2 and 3.3. The method developed by Barkley applies to a particular type of local dynamics, E q . (3.3), and can be very fast but the results are compromised due to the decreased accuracy. The Chebyshev multidomain approach can be used to solve equations of the F H N type with more general local dynamics as shown in Section 3.3. In Section 3.4, we implemented an operator splitting method for the Chebyshev multidomain approach for kinetic model I and II. W i t h the operator splitting scheme, the computational time for the solution of the F H N equations is considerably shorter than both the finite difference and the Barkley 89 Chapter 3. Pseudospectral method of solution method. For kinetic model I, for solutions with similar accuracy, C M D becomes two times faster than the Barkley method and 36 times faster than F D . Similarly for kinetic model II, C M D is 30 times faster than F D . Work is in progress on F H N equations with more complex ion kinetics as defined by the Beeler-Reuter equations [192] are considered in Chapter 4. 90 Chapter 3. Pseudospectral method of solution Table 3.1: Convergence of the wavefront speed versus N, for kinetic model I with the Fourier method and Dirichlet boundary conditions, e = 0.005, 30, x = 30. ,b = 0.01, X R L 128 8.9 N c 256 9.0 512 9.0598 1024 9.059835 2048 9.059835 4096 9.059835 Table 3.2: Convergence of the wavefront speed versus N, for kinetic model I with C M D and no flux boundary conditions, e = 0.005, a = 0.3, b = 0.01, XL — —30, x = 30. N = (N h — 2)N + 2 is the total number of grid points, with N = 40 subdomains and N h is the number of collocation points per subdomain. R c s S c Nch N c 13 442 9.05 20 722 9.059 30 1122 9.05983 50 1922 9.059835 Table 3.3: Convergence of the wavefront speed versus N, for kinetic model I with C M D and no flux boundary conditions, e = 0.005, a = 0.3, b = 0.01, XL = - 3 0 , x = 30. N = (N - 2)N + 2 is the total number of grid points, with N subdomains and N h = 13 collocation points per subdomain. R ch s N s N c S c 40 442 9.05 70 772 9.0598 130 1432 9.05983 180 1982 9.059835 250 2752 9.0598358 350 3852 9.0598358 91 Chapter 3. Pseudospectral method of solution Table 3.4: Convergence of the wavefront speed for kinetic model I with F D and no flux boundary conditions, e = 0.005, a = 0.3, b = 0.01, XL = —30, x = 30. R At \ IO" IO" IO" IO" N 4 5 6 7 128 8 8 - 256 8 8 - 512 S\9 8.9 8.9 8.9 1024 9T67) 9.02 9.02 9.02 2048 9.033 9.049 9.051 9.051 4096 9.04012 9.0559 9.0575 9.0577 8192 9.05754 9.05913 9.05929 Table 3.5: Convergence of the wavefront speed for kinetic model I with the Barkley method and no flux boundary conditions, e = 0.005, a = 0.3, b = 0.01, x L At \ 10"-3 io--4 10"•5 10"-6 = -^30,0^ = 30. 128 N S= 6= 5= 5= S= S= 5= 5= 10" IO" 10~ 10~ 10" 10~ 10" io- 4 5 4 5 4 5 4 5 7 7 7 7 7.489 - 256 8 8.430 8 8 6 8.3 6 512 8.74 8.7 8.797 8.900 7.67 8.81 7.67 1024 8.84 8.85 8.944 9.003 8.28 8.96 8.29 2048 8.9983 9.0310 8.628 9.0138 4.85 8.629 4096 - 9.0213 9.0386 8.8252 9.0369 6.718 8.826 8192 8.93663 9.04782 7.7692 8.93819 92 Chapter 3. Pseudospectral method of solution Table 3.6: Convergence of the wavefront speed versus TY, for kinetic model II with C M D and no flux boundary conditions, e = 0.001, a = 0.1, b = 0.5, A = 1, x = - 2 0 , x = 20. TV" = (N ~ 2)TV + 2 is the total number of grid points, with N subdomains and Nch = 13 collocation points per subdomain. L R ch S s 40 442 17.5 TV TV c S 50 552 17.59 120 1322 17.625 150 1652 17.6251 250 2752 17.625152 350 3852 17.625152 Table 3.7: Convergence of the wavefront speed versus TV, for kinetic model II with F D and no flux boundary conditions, e = 0.001, a = 0.1, A — 1, b = 0.5, x = - 2 0 , x = 20. L R TV c 256 15 512 16.8 1024 17.42 2048 17.574 4096 17.612 8192 17.621 93 Chapter 3. Pseudospectral method of solution ! c b - ^- "(0 ! - \ t a 0.16 1 1 1 1 i d 0.12 e i i -b c v(0 i i i i i ' i (B) 1.2 r 1 e=0.0001 / d i (A) \ ' \ \ 0.8 i c 0.08 uit) 0.4 0.04 h v(0 a • a -0.04 -0.4 _i b i 0 i 0.4 u (C) i i 0.8 i _ - 1.2 -0.4 • e d . • 0 0.4 i i 0.8 i 1.2 (D) Figure 3.1: Local dynamics of the spatially homogeneous case of E q . (3.2) for kinetic models I and II. In A and C the phase portraits for kinetic model I (a = 0.3,6 = 0.01, e = 0.0001) and II (a = 0.1,6 = 0.5, e = 0.0001), respectively are shown. The origin is stable but excitable. For a proper initial condition (position a) the solution undergoes a long excursion before returning to the origin. In B and D , u{t) and v(t) for kinetic model I and II, respectively, are plotted versus time. 94 t Chapter 3. Pseudospectral method of solution Figure 3.2: Numerical solution for kinetic model I, at a specific time, in a square domain of length 60 x 60 with (A) Dirichlet boundary conditions; (B) Periodic boundary conditions. The same initial condition, E q . (3.20), was considered for both cases. 95 Chapter 3. Pseudospectral method of solution Figure 3.3: Time dependent profiles u(x, t*) versus x for t* = 0,1.66, 3.33, 5, for the propagation of a I D pulse with different numerical schemes for kinetic model I (Eq. (3.5) with (3.3)) and a = 0.3, b = 0.001, e = 0.005. (A) Fourier ( • ) , F D ( x ) N = 512; (B) C M D with N = 442 points (•) (N = 40 and N = 13) and F D with N = 512 points ( x ) ; (C) C M D with N = 442 (•) and the Barkley method with N = 512 (0). For (A), Dirichlet boundary conditions u(—30, t) = u(30,t) = 0 and v(—30,t) = v(30,t) = 0.7 were applied. For (B,C), no flux boundary conditions. s ch 96 Chapter 3. Pseudospectral method of solution Fourier x z -12 -16 -20 1 1 0 i i i i ' Barkley -2 if -4 v7/ \ -8 -10 -12 14 i i 18 / 22 x (D) 26 Figure 3.4: Error, EN{X), (Eq. 3.22) versus x for the I D solution of E q . (3.5) for kinetic model I and a = 0.3, b = 0.001, e = 0.005 at t* = 5. (A) Fourier method with N equal to 512 (dashed), 1024 (solid), 2048 (dashed bold), 4096 (solid bold); (B) C M D with N = 13 and N and N equal to (40, 442; dashed), (180,1982; solid) and (5.00,5502; solid bold); (C) F D with N equal to 512 (dashed), 2048 (solid) 8192 (dashed bold) and 16384 (solid bold); (D) Barkley method with N equal to 512 (dashed), 2048 (solid thin), 4096 (dashed bold) and 8192 (solid bold). For the Fourier method, . exact it*), is the last profile i n F i g . 3.3A obtained with Fourier and Af = 8192 points. For the rest of the cases, uf (t*), is taken to be the last profile i n F i g . 3.3C given by C M D with N = 9902 (A^ = 900, N = 13). CH 3 xact CH 97 30 Chapter 3. Pseudospectral method of solution 0.00001 0.0001 0.001 0.01 8 Figure 3.5: Convergence of the numerical wavefront speed obtained with C M D to the asymptotic speed c as e —> 0, with a = 0.3, b = 0.01, D\ = 1. The asymptotic speed c is given by E q . (3.23). a a 98 Chapter 3. Pseudospectral method of solution (C) (D) Figure 3.6: Contour plots of the solution u(x,y,t*) of E q . (3.5) for kinetic model I and a = 0.3, b = 0.001, e = 0.005, u = 0.5 at t* = 8. (A) F D , N = 512 (x) and Fourier TV = 1024 ( • ) ; (B) C M D with N = 442 (N„.= 40, N = 13) (•) and C M D with N = 1432 (JV = 130, N = 13) (o); (C) F D with TY = 512 ( x ) and C M D with N = 1432 (o); (D) Barkley method with iV = 512 (0) and C M D with N = 1432 (o). (A) Dirichlet boundary conditions; (B-D) No flux boundary conditions. ch S ch 99 Chapter 3. Pseudospectral method of solution CMD-CMD > • - - - -30 10 -10 (A) x FD-CMD Figure 3.7: Graphs on the left: I D profiles from cuts of the spirals at y = —15 and x 6 [—30,30], as given by the different numerical schemes at time t* = 8 for kinetic model I. For C M D N = 13. (A) C M D with N and N equal to (442,40; (•)) and (1432,130; (o)); (C) F D with N = 512 ( x ) and C M D with N = 1432 (o); (E) Barkley method with N = 512 (0) and C M D with N = 1432 (o). Graphs on the right: Analysis of EN(X) error versus x for the profiles shown on the left. (B) C M D with N and 100 N equal to (332, 30; dashed), (442, 40; solid) and (1102,100; dashed bold); (D) F D with N equal to (512, dashed), (1024, solid) and (1500, dashed bold); (F) Barkley method with N equal to (512, dashed), (1024, solid) and (1500, dashed bold). The errors were calculated with respect to the numerical solution given by C M D with N = 1432. CH S S Chapter 3. Pseudospectral method of solution -30 -10 10 x 30 Figure 3.8: Comparison of contour plots (u(x,y,t*) = 0.5) of the numerical solution of kinetic model I (Eq. (3.5) with (3.3)) for different numerical schemes at t* = 100. C M D with TV = 442 (N = 40, N = 13) (•) and TV = 882 (N = 80, N = 13) (o), F D (x) and Barkley method (0) with N = 512. s s ch ch 101 Chapter 3. Pseudospectral method of solution 1.2 FD-CMD n 0.8 * 0.4 0.4 -20 -10 10 20 -0.4 -20 10 -10 20 (B) (A) 1.2 FD-CMD ft f] n 0.8 0.4 "ILL -0.4 -20 10 -10 20 (C) Figure 3.9: Time dependent profiles versus x at t* = 0,0.6,1.2,1.8, for the propagation of a I D pulse with different numerical schemes for kinetic model II and a = 0.1, b = 0.5, e = 0.001, A = 1. (A) C M D with N = 13 and AT and N equal to (552,50; (•)) and (8802,800; (o)) (B) C M D with AT = 8802 (o) and F D N = 512 ( x ) ; (C) C M D with N = 8802 (o) and F D , AT = 1024 ( x ) . ch s ; 102 Chapter 3. Pseudospectral method of solution X -4 .-j 2 1 10 1 1 12 1 1 1 1 14 1 1 16 1 18 20 x (B) Figure 3.10: Error, E^(x), (Eq. 3.22) versus x for the I D solution of E q . (3.5) for kinetic model II and a = 0.1, b = 0.5, e = 0.001, A = 1 at t* - 1.8. (A) C M D with N = 13 and values of N and N equal to (552,50; dashed), (iV = 1322,120; solid), (N = 3852,350; solid bold); (B) F D with N equal to (512, dashed), (1024, solid), (4096, dashed bold) and (8192, solid bold). ch s 103 Chapter 3. Pseudospectral method of solution 1.016 1.012 « 1.008 o 1.004 1.000 j 1—i—i—i—> 0.00001 i 111 0.0001 1—i—i—i—i i 111 0.001 8 Figure 3.11: Convergence of the numerical wavefront speed obtained with C M D to the asymptotic speed c as e —* 0 with a = 0.1, b = 0.5, D\ = 1 and A = I. The asymptotic speed c is given by E q . (3.27). a a 104 Chapter 3. Pseudospectral method of solution (C) x x (D) Figure 3.12: Convergence analysis for the 2D solution u(x, y, t*) of E q . (3.5) with (3.4) (kinetic model II) with e = 0.001, a = 0.1, A = 1 and b = 0.5 and t* = 5. (A) Comparison of contour plots (u(x,y,t*) = 0.5) between F D N = 1024 (dashed) and C M D N = 662 (N = 60, N = 13) (•) and the converged solution obtained with C M D , N = 1102 (N = 100, N = 13) (o); (B) I D profiles of u(x,-10,t*) versus x for F D N = 1024 (dashed), C M D TV = 662 (•) and the converged solution obtained with C M D , N = 1102 (o); (C) .Ejv(x) versus x for the profile shown in (B) for C M D with N h = 13 and values of N and N equal to (552,50 ; dashed), (662,60; solid), and (N = 882,80; dashed bold); (D) E (x) versus x for F D N equal to (512, dashed), (1024, solid) and (1500, dashed bold). s ch a ch c s N Chapter 3. Pseudospectral method of solution (B) X Figure 3.13: Comparison of contour plots (u = 0.5) of the numerical solution of (A) kinetic model I for different numerical schemes at t* = 8. C M D with N = 442 (N = 40, N = 13) (solid), OS with N = 442 (•) and OS with N = 386 (N = 128, N = 5) ( x ) ; (B) kinetic model II for t* = 5. C M D with N = 662 (N = 60, N = 13) (solid), OS with N = 662 (•) and OS with N = 602 (N — 200, N h = 5) ( x ) . 106 s ch s ch s a ch c Bibliography [184] F . Amdjadi, J . Gomatam, Spiral waves on static and moving spherical domains, J . Comp. Appl. Math. 182 (2005) 472-486. [185] I. Aranson, H . Levine, L.Tsimiring, Controlling spatiotemporal chaos, Phys. Rev. Lett. 72 (1994) 2561-2566. [186] M . Argentina, P. Coullet, V . Krinsky, Head-on collisions of waves in an excitable Fitzhugh-Nagumo system: a transition from wave anihilation to classical wave behavior, J . Theor. Biol. 205 (2000) 47-52. [187] M . Bar, N . Gottschalk, M . Eiswirth and G . E r t l , Spiral waves in a surface reaction: model calculations, J . Chem. Phys. 100 (1994) 12021214. [188] E . Barillot, J . Boissonade, Asymptotic pseudospectral method for reaction-diffusion systems, J . Phys. Chem. 97 (1993) 1566-1570. [189] D . Barkley, M . Kness and L . Tuckerman, Spiral-wave dynamics in a simple model of excitable media: The transition from simple to compound rotation, Phys. Rev. A 42 (1990) 2489-2493. [190] D . Barkley, A model for fast computer simulation of excitable media, Physica D , 49 (1991) 61-70. [191] D . Barkley, Spiral meandering, in Chemical Waves and Patterns pp 163-189 by R. Kapral and K . Showalter, (Kluwer Academic 1995). [192] G . W . Beeler and H . Reuter, Reconstruction of the action potential of ventricular myocardial fibres, J . Physiol. 268 (1977) 177-210. [193] J . P. Boyd, Chebyshev and Fourier Spectral Methods (Dover, New York, 2000). [194] F . Brauer and C . Castillo-Chavez, Mathematical Models i n Population Biology and Epidemiology, Texts in Applied Mathematics (Springer, New York, 2000). 107 Bibliography [195] C . Canuto, M . Y . Hussaini, A . Quarteroni, T . A . Zang, Spectral Methods in Fluid Dynamics, (Springer, New York, 1988). [196] E . M . Cherry, H . S. Greenside, C . S. Henriquez, A space-time adaptive method for simulating complex cardiac dynamics, Phys. Rev. Lett. 84 (2000) 1343-1346. [197] W . S. Don, D . Gottlieb, J . - H Jung, A multidomain spectral method for supersonic reactive flows, J . Comput. Phys. 192 (2003) 325-354. [198] M . Dowle, R . Mantel, D . Barkley, Fast simulations of waves in threedimensional excitable media, Int. J . Bifurcation and Chaos 7 (1997) 2529-2545. [199] E . Cytrynbaum, J.P. Keener, Stability conditions for the traveling pulse: Modifying the restitution hypothesis, Chaos 12 (2002) 788-799. [200] C . Diks, B . Hoekstra, J . Degoede, Spiral wave dynamics, Chaos Sol. Frac. 5 (1995) 645-660. [201] J . C. Eilbeck, A collocation approach to the numerical calculation of simple gradients in reaction-diffusion systems solutions, J . Math. Biol. 16 (1983) 233-249. [202] I. R . Epstein, J . A . Pojman, A n Introduction to Nonlinear Chemical Dynamics (Oxford University Press, Oxford, 1998). [203] F . Fenton, E . Cherry, H . M . Hastings, S. J . Evans, Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity, Chaos 12 (2002) 852-892. [204] B . Fornberg, A Practical Guide to Pseudospectral Methods (Cambridge University Press, Cambridge, 1996). [205] R . Fitzhugh, Impulses and physiological states in theoretical models of nerve membrane, Biophy. J . 1 (1961) 445-466. [206] L . Glass, M . Josephson, Resetting and annihilation of reentrant abnormally rapid heartbeat, Phys. Rev. Lett. 75 (1995) 2059-2062. [207] D . Gottlieb, J.S. Hestheaven, Spectral methods for hyperbolic problems, J . Comput. Appl. Math. 128 (2001) 83-131. [208] D . Gottlieb and S. A . Orszag, Numerical Analysis of Spectral Methods ( S I A M , Philadelphia, 1977). 108 Bibliography [209] G . Gottwald, A . Pumir, V . Krinsky, Spiral wave drift induced by stimulating wave trains, Chaos 11 (2001) 487-494. [210] R . Gray, Termination of spiral wave breakup in a Fitzhugh-Nagumo model v i a short and long duration stimuli, Chaos 12 (2002) 941-951. [211] W . B . Jones, J . J. O'Brien, Pseudo-spectral methods and linear instabilities i n reaction-diffusion fronts, Chaos 6 (1996) 219-228. [212] A . Karma, Meandering transition in two-dimensional excitable media, Phys. Rev. Lett. 65 (1990) 2824-2828. [213] A . Karma, Scaling regime of spiral wave propagation in single-diffusive media, Phys. Rev. Lett. 68 (1992) 397-400. [214] J . P. Keener, A geometrical theory for spiral waves in excitable media, S I A M J . Appl. Math. 46 (1986) 1039-1056. [215] J . P. Keener, Mathemathical Physiology, Interdisciplinary Applied Mathematics (Springer, New York, 1998) [216] D . A . Kessler, H . Levine and W . N . Reynolds, Theory of the spiral core i n excitable media, Physica D 70 (1994) 115-139. [217] V . Krinsky, A . Pumir, Models of defibrillation of cardiac tissue, Chaos 8 (1998) 188-203. [218] Y . Kuramoto, Chemical Oscillations, Waves and Turbulence (Dover, New York, 1984). [219] K . Lee, Wave pattern selection in an excitable system, Phys. Rev. Lett. 79 (1997) 2907-2910. [220] A . Z. L e i , et al, Spiral wave maintained by noise, Fluctuation and Noise Letters, 4 (2004) 447-452. [221] C . - H . Luo, Y . Rudy, A dynamic model of the cardiac ventricular action potential, Circ. Res. 74 (1994) 1071-1096. [222] D . Margerit, D . Barkley, Selection of twisted scroll waves in threedimensional excitable media, Phys. Rev. Lett. 86 (2001) 175-178. [223] I. Mitkov, I. Aranson, D . Kessler, Meandering instability of a spiral interface i n the free boundary limit, Phys. Rev. E 54 (1996) 6065-6069. 109 Bibliography [224] J . D . Murray, Mathematical Biology I and II, Interdisciplinary Applied Mathematics, (Springer, New York, 2002). [225] J . Nagumo, S. Arimoto, S. Yoshizawa, A n active pulse transmision line simulating nerve axon, Proc. I R E 50 (1962) 2061-2070. [226] D . Olmos, B . Shizgal, A Spectral Method of Solution of Fisher's Equation, J . Comp. Appl. Math. 193 (2006) 219-242. [227] N . Otani, A primary mechanism for spiral wave meandering, Chaos 12 (2002) 829-842. [228] A . Panfilov, P. Hogeweg, Spiral breakup in a modified FitzhughNagumo model, Phys. Lett. A 176 (1993) 295-299. [229] A . Panfilov, J . P. Keener, Effects of high frequency stimulation on cardiac tissue with an inexcitable obstacle, J . Theor. Biol. 163 (1993) 439-448. [230] V . Petrov, S.K. Scott, K . Showalter, Excitability, wave reflection and wave splitting i n a cubic autocatalysis reaction diffusion system, P h i l . Trans. R. Soc. Lond. A 347 (1994) 631-642. [231] R . Peyret, Spectral Methods for Incompressible (Springer, New York, 2002). Viscous Flow [232] Z. Qu, A . Garfinkel, A n advanced algorithm for solving partial differential equation i n cardiac conduction, I E E E Trans. Biomed. E n g . 46 (1999) 1166-1168. [233] D . L . Ropp, J . N . Shadid, Stability of operator splitting methods for systems with indefinite operators: reaction-diffusion systems, J . Comp. Phys. 203 (2005) 449-466. [234] T . Shardlow, Nucleation of waves in excitable media by noise, M u l t i scale Model. Simul. 3 (2004) 151-167. [235] J . M . Starobin, C . F . Starmer, Common mechanism links spiral wave meandering and wave-front-obstacle separation, Phys. Rev. E , 55 (1997) 1193-1196. [236] J . A . Trangenstein, C . K i m , Operator splitting and adaptive mesh refinement for the Luo-Rudy I model, J . Comp. Phys. 196 (2004) 645-679. 110 Bibliography [237] K . H . W . J . ten Tusscher, A . V.' Panfilov, Wave propagation in excitable media with randomly distributed obstacles, M u l t . M o d . Simul. 3 (2005) 265-282. [238] J . Tyson, J . Keener, Singular perturbation theory of travelling waves in excitable media (A review), Physica D 32 (1988) 327-361. [239] J . Tyson, What everyone should know about the BelousovZhabotinsky reaction, Lecture notes in biomathematics, Springer, V o l . 100 (1994) 569-587. [240] A . T . Winfree, Varieties of spiral wave behavior: A n experimentalist approach to the theory of excitable media, Chaos 1 (1991) 303-334. [241] F . X i e , Z. Q u , J . Weiss, A . Garfmkel, Coexistence of multiple spiral waves with independent frequencies in a heterogeneous excitable medium, Phys. Rev. E 63 (2001) 031905 1-4. [242] Z. Zhu, H . X i n , Dynamics of spiral waves under the modulation of noise pulses, Phys. Rev. E , 64 (2001) 056124(1-5) 111 Chapter 4 Annihilation and reflection of spiral waves at a boundary for the Beeler-Reuter model. A n understanding of the propagation of waves in excitable media is a very important research area [252, 257, 270, 274] . Particular attention has been given to the study of spiral waves [251, 252, 272, 273, 274]. A spiral wave is a self sustained wave that rotates freely or around some obstacle. In cardiology, it is thought that spiral waves, which have been observed experimentally in isolated cardiac tissue [249, 267], play an important role in the development of certain arrhythmias as well as in the onset of fibrillation [252, 263]. 4 A very important feature of spiral waves is the behaviour of its tip, as discussed by Fenton et al [252] and Comtois et al [248]. For a particular set of parameters characterizing the medium, the spiral tip executes different trajectories which can be circular or more complicated patterns [243, 261, 273]. This phenomena is referred to as spiral meandering and has been studied with different kinetic models [244, 251, 252, 273]. A second feature of the spiral tip, called spiral wave drift, is the response of the spiral wave to an external perturbation. Some examples of external perturbations, include interaction between two spirals [275], interaction of spirals with a boundary [256, 276] and the drift induced by light as in the Belousov-Zhabotinsky reaction [258]. In the present chapter, we study numerically a spiral wave in the meandering regime, and the drift effects of the boundary of the domain on the tip trajectory of the spiral. In the simulations, we observe two different phenomena A version of this chapter has been submitted for publication. D a n i e l O l m o s and Bernie Shizgal, A n n i h i l a t i o n and reflection of spiral waves at a boundary for the Beeler-Reuter model, P h y s i c a l Review E. 4 112 Chapter 4. Annihilation and reflection of spiral waves due to the interaction of the spiral wave with a boundary, which are reflection, and annihilation of the spiral wave. Experiments in isolated cardiac tissue [249, 259, 260, 267], have been carried out to study the interaction of spiral waves with obstacles and the boundaries of the tissue. Davidenko et. al. [249], Pertsov et. al. [267], and Ikeda et. al. [259] showed annihilation of spiral waves at the boundary. Also, Ikeda et. al. [260], observed attachment of meandering spirals to an obstacle of some minimum size. Therefore, the phenomena of annihilation and reflection observed in the computer simulations in this chapter, suggest that the interaction of spiral-boundary and spiral-obstacle not necessarily end in annihilation at the boundary or attachment to the obstacle as observed experimentally. This is a very important issue in cardiology, and a proper understanding of the conditions that lead to the annihilation of the spiral at some inexcitable obstacle, can help in the development of techniques on how to eliminate arrhythmias generated by the spiral behaviour [251, 259]. Spiral drift due to boundary effects has been considered previously. It has been found experimentally [256] and studied numerically with no-flux boundary conditions [276]. Gomez-Gesteira et al. [256] considered the Belousov-Zhabotinsky reaction and found that the boundary affected the trajectory of the spiral tip. The trajectory moved along the boundary, whereas in other cases the spiral was annihilated at the boundary [256]. Yermakova and Pertsov [276], analyzed the effects of the boundary on the trajectory of the spiral tip that followed a circular path. B y considering no flux boundary conditions, they showed that the period of the spiral increases when the core of the spiral is close to the boundary. Also, they showed that the center of the circular trajectory drifted at constant speed along the boundary giving as a result a trajectory resembling the shape of a trochoid. However, the case when the trajectory of the spiral tip meanders and traces a more complex pattern other than a circle has not been considered. The phenomena of meandering, also referred as compound rotation, was first noted by Winfree [271]. Figure 4.1 illustrates the two types of trajectories of the tip of a spiral wave constructed by considering the motion of the tip of the arrow attached to the small circle of radius h that rotates either on the inside or outside of the large circle of radius R. When meandering occurs, the trajectory of the tip executes a flower like pattern [243, 273], where the petals lie on a circle of radius R. When the petals lie outside the circle, the trajectory resembles a curve called a hypotrochoid (Fig. 4.1A), 113 Chapter 4. Annihilation and reflection of spiral waves whereas in the case the petals lie inside the circle, the trajectory resembles an epitrochoid (Fig. 4.IB). B y considering different parameter values in a particular model with excitable kinetics [243, 251], it is possible to take the limit R —> oo. In this case, which we refer to as the limiting Roc case with R » 1, the flower has almost an infinite radius and the petals lie essentially on a straight line. In the present chapter, we consider the case when a spiral meanders and interaction with the boundary takes place. The meandering effects observed are different and more complex compared to the case when the spiral tip follows a circular motion. When the spiral meanders and is close to the boundary (with no flux boundary conditions), it is observed that the trajectory is annihilated or reflected by the boundary. In the case where the trajectory is reflected, the angle of reflection is not necessarily equal to the angle of incidence. Also, the reflection angle is very sensitive to the position along the spiral tip trajectory, at which the trajectory hits the boundary. Therefore, the main question we address here is to find the conditions for which the spiral is annihilated at the boundary. In order to analyze the effects of the boundary on a meandering spiral, we focus attention on the degenerate i?oo limiting case. The infinite radius regime is just a transition from the outward petal to the inward petal flower tip trajectory and therefore is not a generic behaviour [243]. The analysis of the ROQ case, is considered due to its simplicity compared to the case of finite R. Near the boundary, the behaviour of the tip of a spiral can be approximated by the R^ case, and therefore, the results obtained for this limit may also provide an understanding for the case when R is finite. In Section 4.1, we present the model equations used i n the numerical experiments and we also provide a description of the numerical method employed in their solution. In Section 4.2, we present the results of numerical simulations for different values of R, including the Roo case. Some of the boundary effects observed in compound rotation are shown for the i?oo m Section 4.3. In Section 4.4, we present the phenomenon of annihilation at a boundary as a function of the incident angle of the trajectory obtained with the Roo case. In Section 4.5, an argument based on the reactivation variable j for the sodium channels is discussed to explain the phenomenon of annihilation and reflection of a spiral at a boundary. Due to the complexity of the problem, we present a qualitative rationalization to explain why the probability of annihilation varies with respect to the incident angle of the c a s e 114 Chapter 4. Annihilation and reflection of spiral waves trajectory for the i?oo case (Section 4.6). Finally, we present some conclusions i n Section 4.7. 4.1 Membrane models and numerical methods. Spiral wave dynamics in excitable media and particularly trajectories of the spiral tip, have been extensively studied with reaction diffusion partial differential equations [251, 252, 273] as well as reduced O D E models [243, 246]. The simplest model used is based on the Fitzhugh-Nagumo equations which have been discussed by'different authors [250, 254, 262, 273]. In this chapter we considered the kinetic model given by Beeler and Reuter [245], which has more realistic dynamics than the Fitzhugh Nagumo equations. The Beeler-Reuter (BR) model, which is based on the HodgkinHuxley formalism, consists of a set of eight coupled nonlinear ordinary differential equations described in [245]. The transmembrane potential V satisfies the equation dV 1 = ~~~Q~ (lion ~ lapp) (4-1) where C is the membrane capacitance, I is a stimulus current, and Ii corresponds to the sum of four ionic currents I , I , INO. and Ic where their form is given by Eqs. (A.1-A.4) in Appendix A . J/v is the fast inward current and is carried by sodium ions, whereas Ic , is the slow inward current and is carried by calcium ions. I^ and Ic are voltage and time dependent currents. I , I are the time independent and time dependent outward currents, respectively, carried mostly by potassium ions. I , I^ and Ic are controlled by gate variables. Each of the six gate variables satisfies the relationship given by E q . (A.6) in Appendix A . The B R model also contains the equation that controls the intracellular calcium concentration of the cell and is given by E q . (A.5). The B R model is extended to spatiotemporal dynamics by adding a diffusion term in the membrane potential equation (Eq. 4.1) giving m app on kl X1 a a a a kl a X1 X1 ^ = D V V - ^ - (I 2 im - I p) ap a a (4.2) where D = a/S C is the diffusion coefficient for the isotropic case where a refers to the conductivity of the medium and S is the ratio of cell surface area per unit of volume [253, 262]. For this Chapter the values of D and C are 0.1 m m m s and 1/nFcm respectively [251, 252]. v m v 2 _ 1 -2 m 115 Chapter 4. Annihilation and reflection of spiral waves The numerical method used in this chapter is a variation of the method considered by Olmos and Shizgal [266]. The present calculations were carried out over a square domain Q, = [—L,L] x [—L,L]. Due to the multiple spatial and temporal scales in the B R equations (Eq. 4.2), we considered the operator splitting method [268, 269] to separate the diffusion and reaction processes. For spatial discretization, which is considered in the diffusion process, we used a non uniform grid as defined with a multidomain method based on Chebyshev collocation points [265, 266]. In order to solve the B R equations, we considered a procedure similar to the operator splitting algorithm presented by Qu and Garfinkel [268] with some minor differences. The calculations for the reaction part were speeded up by dividing the domain Q into smaller squares II of size 30 points in each dimension. In each square II, we calculated the values \d V/dx \ and \d V/dy \ and as soon as either of these values were larger than 5 at some point (x*,y*) € II, we defined the whole square as a place where the front might be located. In the same way, all the eight squared neighbors of II were considered in the same way as the region II to assure correct propagation of the wavefront. We finally obtain two main regions, one where the fast changes in time and space occur and one with no fast changes. In the region where V changes very slowly (second derivatives less than 5), we solve the dynamics as in [268]. Because V changes slowly, we consider V as constant for a time step and solve E q . (A.6) for V respect to time analytically, by taking Yi and T j constants for that interval of time. In order to calculate the a's and /3's in Eqs. (A.7) and (A.8), we proceed as in [251] where tabulated values given by step function approximations of ct(V) and P(V) as a function of V are considered. In LT, the region where fast changes of V occur, we solve E q . (A.6) numerically with an explicit second order Runge-Kutta scheme with a time step A t = 0.01, four times smaller than in the region with no fast changes. For the multidomain approach used in this chapter [265, 266], we choose Ni = 150 and Ni = 180 subintervals with N h = 5 points per subinterval, giving a total of N — 452 and N = 542 points in each dimension, respectively. The size of the domain is 2 L x 2L with L = 40 mm. Convergence to the solution for the one dimensional problem was obtained with N = 752 points. However, the results obtained with TV = 452 and N = 542 points provided a very good approximation to the converged solution. In this chapter, we present studies of the evolution of the spiral wave tip in time. There are different ways to define the tip of the spiral and a good summary can be found in [252]. In this chapter, the tip of the spiral was calculated by taking a level curve of the V variable (V = —25 m V ) , and then finding the point with largest curvature along this curve. The evolution in 2 2 2 2 c 116 Chapter 4. Annihilation and refection of spiral waves ... time of the position of the tip defines a curve called the spiral tip trajectory. The computational time for the simulations considered in this Chapter is calculated for the time that a spiral traces a complete rotation. When N = 452, a complete rotation takes about two hours to compute and represents about 110 ms of simulation time. For ./V = 542, the same calculation takes almost four hours to compute. 4.2 Numerical Results. Figure 4.2 shows the different spiral tip behaviours by taking the calcium conductance gs = 0.03 m S c m as done in [251], and varying Ej^ and g^ over a range of values. A circular trajectory, also known as simple rotation, is shown in Fig. 4.2A. The reminder of the figures portray the phenomenon referred as compound rotation. The trajectories shown in Figs. 4.2(G-J) have a shape very similar to a family of curves called hypotrochoids and are referred to have outward petal flower patterns as in F i g . 4.1 A [243]. O n the other hand, the trajectory shown in Fig. 4.2E has an inward flower like motion as in F i g . 4.IB, very similar to an epitrochoid. Figures 4.2B and 4.2F represent the transition from an inward to an outward flower like tip trajectory and they correspond to the case when the radius R in either case the hypotrochoid or the epitrochoid, goes to infinity. Figures 4.2(G-I) show that the trajectories undergo transitions from a 5 to 4 to 3 outward petal flower. The transition from Fig. 4.2G to 4.2H and the one from 4.2H to 4.21 occurs gradually and intermediate trajectories are similar to those shown in 4.2C and 4.2D. - 2 a a I focus attention to the parameter values gs = 0.03 m S c m ~ , g^ = 2.38 m S c m and N = 452 collocation points in each dimension, which give the Rev case. W i t h these parameters, the tip of the spiral wave meanders as shown in F i g . 4.2B, that is, it follows a straight line. Even though there appears to be pattern following a straight line, it is not possible to be certain that the R^ has been achieved. For practical purposes when R » 1 will be considered as the R^ case. In F i g . 4.3A, five trajectories are shown for different initial conditions for an integration time of t* « 1.75 s. The initial conditions were constructed following the procedure given in Appendix B with t = 50 ms, c\ = 30, C2 = 25, a = 0 and varying the parameter yo from —10 m m to —14 mm with increments of 1 mm. As a result, the dif2 a - 2 g 117 Chapter 4. Annihilation and reflection of spiral waves ferent initial conditions are a copy of each other translated horizontally by a distance of approximately 0.8 mm. In F i g . 4.3A, all the trajectories hit the boundary and get reflected with the same angle for the first reflection but not the second. This response is due to the fact that each trajectory is an horizontal translation of the other. Therefore, the distance from the initial position of the tip trajectory to the boundary is exactly the same for all five trajectories. The time required for all the trajectories to hit the boundary for the first time is t — 0.596 s. However, when the trajectories hit the boundary a second time (Fig. 4.3A), the angles of reflection of each trajectory is different. Each trajectory has to travel a different distance from each other before hitting the boundary the second time. The corresponding times when the tip trajectories hit the boundary for the second time are t* = 1.0326,1.0376,1.0426,1.0476 and 1.0526 s, respectively. The trajectories shown in F i g . 4.3A are shown in Figs 4.3(B-F) for longer integration times up to 5 seconds. The initial position of the trajectory is marked with an asterisk. It is clear that the solutions follow completely different paths after the second reflection from the boundary. Moreover, in Figs. 4.3(B-D) it is shown that the spiral tip instead of being reflected at the boundary, annihilates as shown by the arrows at times t « 3.62,1.82 and 3.22 s, respectively, ending the spiral motion. B y contrast, the trajectories shown in Figs. 4.3E and 4.3F, which are shown for times up to t ss 4.3 and 3.62 s, respectively, remain inside the domain even for times as large as 5 s. The spiral tip in Figs. 4.3E and 4.3F was followed until the points indicated by solid circles. Notice also that in Figs. 4.3(B-D), the spiral tip hits the boundary 5, 3 and 5 times, respectively, before they annihilate at the boundary. The spiral tip in Figs. 4.3E and 4.3F gets reflected at the boundary 6 and 5 times respectively and they still remain in the domain. It is important to mention that spiral tips never actually touch the boundary when the spiral tip is reflected. The tip of a spiral touches the boundary only when the spiral gets annihilated. However, we refer to these two events in general as hitting the boundary. In the next Section, we discuss this concept in a more precise way. The phenomenon of reflection and annihilation at the boundary also occurs for finite ii!, as shown in F i g . 4.4. In Figs. 4.4A and 4.4B, an outward petal flower pattern with R « 30 mm is shown, whereas in Figs. 4.4C and 4.4D an inward petal flower with R ss 16 mm occurs. The radius R in these trajectories was estimated from the semicircles in bold in Figs. 4.4A and 118 Chapter 4. Annihilation and reflection of spiral waves 4.4C, for R « 30 m m and R « 16 mm, respectively. In order to obtain the outward petal flower, the parameters in E q . (4.2) are the same as i n the i?oo case but with gjvo = 2.41 m S c m . We used analogous to the R^ case, two different initial conditions with yo = 10 mm and t equal 50.5 and 51.5 ms, such that each trajectory is a short vertical translation from the other. The corresponding trajectories are plotted in Figs. 4.4A and 4.4B for an integration time t* = 5 s. The initial positions of the spiral tip trajectories are marked with an asterisk, and the final position either by a solid circle or an arrow as in F i g . 4.3. From Figs. 4.4A and 4.4B, it is clear that after the first interaction with the right boundary, the response is exactly the same. However, at the second reflection, the angles of reflection at the bottom boundary are different. When the third reflection takes place at the left boundary, the trajectories in Figs. 4.4A and 4.4B get reflected a third time. Annihilation of the spiral wave occurs in F i g . 4.4B when the spiral tip hits the boundary for the fourth time (arrow), whereas for 4.4A, the trajectory still remains inside the domain (solid circle). - 2 g In Figs. 4.4C and 4.4D the same phenomenon is shown for an inward petal flower pattern. To obtain this flower pattern, we considered the same parameters as in the i?oo d set g^ = 2.34 m S c m . In Figs. 4.4C and 4.4D, we show two spiral tip trajectories, where the initial conditions were calculated with the same procedure as in Figs. 4.4A and 4.4B and the integration time is t* = 5 s. The parameters were taken as yo = — 10 m m and t equal to 85 and 83 ms, respectively. The trajectories for these two initial conditions are a small vertical translation from each other, until the tips of the spiral hit the boundary for the second time, where the response is completely different as shown in Figs. 4.4C and 4.4D (bottom boundary). In Fig. 4.4C the spiral tip is reflected the second time from the lower boundary and then follows the bold semicircle and is reflected a third time also from the lower boundary. The trajectory in Fig. 4.4C involves six reflections with the spiral tip disappearing at the boundary ending with the rotating motion. B y contrast, the spiral tip trajectory in F i g . 4.4D, after the first reflection experiences reflections from the bottom and left boundaries and then subsequently appears to move laterally the right and upper boundaries without being annihilated (solid circle). The radius R for the movement of the spiral tip laterally along the boundary is difficult to calculate. c a s e a n - 2 a g 119 Chapter 4. Annihilation and reflection of spiral waves . . . 4.3 Boundary effects on the rotation period in compound rotation. The study of compound rotation is much more complex than in the case of simple rotation. In simple rotation, when the system has no external perturbations, a trajectory far from the boundary remains circular. However, when such a trajectory is close enough to the boundary, drift occurs [276] and the trajectory changes from a circle to a curve called a trochoid. Another phenomenon observed in [276] is that the rotation period decreases when the distance of the center of the circular trajectory to the boundary is decreased. Annihilation of the spiral wave when its tip trajectory is circular occurs when the tip trajectory is too close to the boundary as described in [264]. In compound rotation we focus on the R^ case to study the annihilationreflection phenomenon. B y considering the Roc case we remove the variable R as a parameter, avoiding the possible complicated behaviours shown in Fig. 4.4D. Also, the trajectories obtained with the Roc case can be seen as a local linear approximation of a trajectory with finite R when the spiral tip is about to interact with the boundary. Therefore, studies with the Roo case will provide information about the annihilation-reflection properties for R finite. Then, for the following analysis we solve the standard B R equations [251] with gN = 2.37 m S c m and gs — 0.03 m S c m with a fine grid in each spatial dimension and iV = 542 collocation points. W i t h these parameters we obtain a trajectory of the tip with R » 1 which for practical effects is considered to be in the Roo case. - 2 - 2 a Figure 4.5 shows two different reflections of a linear spiral tip trajectory from a boundary of the domain. In F i g . 4.5A we show a reflection of such a spiral tip and define the angle of incidence 6i, and the angle of reflection 0 with the dashed lines which are tangent to the petals of the trajectory. For the Roo case, we consider that the trajectory tip has reflected from the boundary when there is a change in the direction of the line on which the petals lie and this direction is maintained for two or more spiral rotations. The tip is considered to hit the boundary at the first point along the trajectory that reaches a minimum distance with the boundary (filled circle in Fig. 4.5A). A different type of reflection is shown in F i g . 4.5B, where there are two positions along the tip trajectory at which the minimum is reached (asterisk and filled circle). In this case, we consider the second point (filled r 120 Chapter 4. Annihilation and reflection of spiral waves circle) as the place where the tip hits the boundary. Figure 4.6A, drawn for one rotation of the spiral wave, shows the trajectory starting from the solid circle to the arrow. This curve defines a unit of trajectory and consists of a petal and an arc. The period of rotation for the chosen parameters is To = 120 ms and the arc length of the unit of trajectory is Lo = 24.8 mm. The trajectories are thus taken to be made up of consecutive unit trajectories. However, close to the boundary, the shape, length and period of rotation of a unit of trajectory are not the same as for those units far from the boundary as shown in Figs. 4.6(B-D). In Fig. 4.6B the unit of trajectory is conserved for the first two rotations of the spiral (starting from the asterisk) and for the third unit of trajectory (bold) the shape is changed, owing to the reflection from the boundary. The unit of trajectory after reflection (dashed) appears to be once again coincident with the first rotation before interaction with the boundary. A similar behaviour is observed in F i g . 4.6C. In F i g . 4.6D the same behaviour as in F i g . 4.6(B-C) is observed, but in this case the trajectory disappears at the boundary. In Table 4.1 we compare the period of rotation TO and the arc length Lo for the first two units of trajectory (starting from the asterisk) which are far from the boundary, the unit of trajectory in bold, and the unit of trajectory in dashed, given in Figs. 4.6(B-D). In Table 4.1, the period of rotation and arc length the two trajectories (bold and dashed) in F i g . 4.6B have a reduced period of rotation and arc length respect to the first two units, which are far from the boundary. The opposite effect is observed in F i g . 4.6C, i.e. the period of rotation and the arc length of the unit trajectories in bold and dashed lines have increased relative to the first two units far from the boundary. The last case, given in Fig. 4.6D, shows that the unit trajectory in bold has shorter period and length than a normal unit (Table 4.1). In this case, the tip trajectory disappears at the boundary and no dashed unit exists. 121 Chapter 4. Annihilation and reflection of spiral waves . . . 4.4 Annihilation and reflection as a function of the incident angle. The results in Fig. 4.3 show that annihilation and reflection of a spiral wave can occur for the same #j. In order to understand the factors that determine absorption and reflection of spiral waves, we consider a fixed angle Oi with different initial conditions. A spiral tip trajectory with a specific incident angle 8i can be constructed as described in Appendix B . For the results that follow, we take g^ = 2.37 m S c m and gs = 6.03 m S c m with N = 542 points, which we take as the Roo case. In order to obtain a trajectory with incident angle 9i = 120°, we choose initial conditions such that yo = 0, ci = 75 and c = 72 and a = 1.53. W i t h this procedure, the spiral tip trajectories shown in F i g . 4.7A, with t = 105 ms and t = 150 ms are prepared with the same 9i. In both cases the spiral tip disappears at the boundary. - 2 - 2 a 2 g g The trajectory at the right in F i g . 4.7A is an horizontal translation of the one at the left but with an extra petal. The extra petal is due to the different values of t for each of the two trajectories. W i t h t = 150 ms, the initial position (asterisk) of the tip trajectory is farther from the lower boundary by a distance A„ = 8.78 mm, than the initial position (asterisk) of the tip trajectory, with t = 105 ms. The quantity A„ gives the distance travelled by the tip in the vertical direction during a complete rotation of the spiral wave. As a consequence an extra petal is observed for the trajectory on the right. Since A„ = A sin the wavelength A is the distance between two petals that are far from the boundary. If a third tip trajectory is generated with the same 9i but with the initial position of the tip 2A„ m m farther from the lower boundary than the initial position of the tip trajectory with t = 105 ms, this new trajectory will have two extra petals compared to the trajectory obtained with t = 105 ms. This third trajectory is found with t = 195 ms. g g g g g g Two trajectories with initial positions separated by a vertical distance X such as in F i g . 4.7A, exhibit the same behaviour when the corresponding spirals hit the boundary, analogous to the situation in F i g . 4.3A for the first reflection. In F i g . 4.3A, the initial position of the tip trajectory is at the same distance from the boundary for all the trajectories and the same behaviour is observed. Therefore, the outcome of a spiral wave hitting the boundary is periodic as a function of the distance of the initial position of v 122 Chapter 4. Annihilation and reflection of spiral waves its tip trajectory with respect to the boundary, and the period is In order to obtain trajectories such that the initial position of their tip trajectory lie within the period X , we take values of t € [105,150]. Different values of t inside this interval, will give trajectories whose initial distance from the boundary is different from each other. The distance of the initial positions of the trajectories with t = 105 and t = 150 from the boundary, will differ exactly by A„ = 8.8 mm. Therefore, all the possible cases by which a spiral wave can hit the boundary for 9i = 120°, can be obtained by taking t £ [105,150]. Taking a smaller interval for t , will leave out cases by which the spiral wave can interact with the boundary, and a larger interval will repeat cases by which the spiral wave can interact with the boundary. v g g g g g g The trajectories with the spiral tips calculated with t = 122,129 and 136 ms are shown in Figs. 4.7(B-D), respectively. In Figs. 4.7(B-D), the initial positions (asterisks) of the trajectories with t = 122,129 and 136 ms, are farther from the boundary by 3.3, 4.6 and 6 mm, respectively, than is the initial position (asterisk) of the tip trajectory with t = 105 ms. g g g W i t h t = 105 ms, the tip trajectory changes its direction abruptly and then disappears at the boundary (Fig. 4.7A). When t is increased to 122 ms (Fig. 4.7B), the trajectory does not disappear at the boundary immediately but a new petal, shown in bold, is formed first. W i t h a further increase in t , this new petal gets closer to where the trajectory hits the boundary. From Fig. 4.7C, the new petal moves until it appears before the tip hits the boundary. In frame D , a new spiral unit (bold) has been formed. The new spiral unit is completely formed when t — 150 ms (Fig. 4.7A), as expected, by the definition of A^. The phenomenon of the creation of a new petal repeats as we increase further the value of t above 150 ms. A new petal will be completely formed when t = 195 ms, which gives a trajectory that starts A„ = 8.8 mm farther from the boundary than the trajectory obtained with t = 150 ms. g g g g g g g In Fig. 4.8, the trajectories were generated with the procedure in Appendix B , with y = - 2 0 mm, c = 35, c = 33, a = 0.123 such that 6 = 70°. Following the same procedure as for di = 120°, we consider two tip trajectories given by the values t = 44 and t = 71 (Fig. 4.8A). Analogous to the case with 9i = 120°, the trajectory at the right in Fig. 4.8A is an horizontal translation of the one at the left but with an extra petal. In this case, the initial position (asterisk) of the tip trajectory with t = 71 ms is X — 8.9 mm farther from the bottom boundary than is the initial position (asterisk) 0 x g 2 t g g v 123 Chapter 4. Annihilation and reflection of spiral waves of the trajectory with t = 44 ms. Following the same reasoning as in the previous example (0; = 120°), we take representative values of t G [44,71]. In Figs. 4.8B to 4.8D, we show three trajectories with intermediate values of t equal to 52,60 and 67 ms, respectively. These three trajectories have the initial position (asterisk) of their trajectories, 2.7, 5.4 and 7.8 m m farther from the boundary compared to the initial position of the trajectory given by t = 44. g g g g In F i g . 4.8A, where t = 44 ms, the tip hits the boundary and the tip trajectory disappears. However, when t is increased to 52 ms (Fig. 4.8B), the spiral wave gets reflected at the boundary. Therefore, annihilation and reflection of the spiral wave is observed with the same angle Bi. In Figs. 4.8(B-D), we show the evolution of the change of the shape of the trajectory unit that interacts with the boundary. For t = 71 ms (Fig. 4.8A), the distance from the initial position of the trajectory to the boundary, is farther A„ = 8.9 mm than the same distance for the trajectory with t = 44 ms. Therefore, the response of the trajectories with t = 44 and t = 71 ms is the same. Increasing the value of t beyond 71 ms will repeat the observed behaviour in Figs. 4.8(B-D), which illustrates the periodicity with period X . g g g g g g g v Annihilation of the spiral wave with Bi = 120° and annihilation and reflection of the spiral with Bi = 70° shown in Figs. 4.7 and 4.8, respectively, for different values of t , suggests to question for which angles Bi annihilation is observed. Moreover, are there particular values of 0i that make more favorable the phenomenon of annihilation than others? Are there values Bi where annihilation does not occur? These questions are addressed in F i g . 4.9, where we show the fraction F^Bi) of the spirals that are annihilated at the boundary as a function of the incident angle In this graph, we considered the range of angles from 20° to 160° with steps of 10°. For each angle 9i, we found a pair of parallel tip trajectories with a pair of values tg and t , such that the initial position of the tip trajectory with t is a distance A„ = Asin#i mm, farther from the boundary than the initial position of the tip trajectory with t . As discussed for 9i = 120°, all the possible cases by which a spiral wave hit the boundary can be obtained by taking t £ [tj,* ,] only. For this interval, the spiral has completed exactly one rotation. Therefore, we sample 20 to 30 tip trajectories for each angle Bi with equidistant values of t between t € [ig,^ ,]. The equidistant values of t E [*g,*g], correspond to different trajectories which their initial position is farther from the boundary than the initial position of the trajectory obtained with t and the distance by which they differ are multiples of \ / g g g g 2 g 2 g g g g v 124 Chapter 4. Annihilation and reflection of spiral waves (number of sample trajectories). Notice that if we take a smaller interval than [tg,tp] for our sampling, we will leave out many cases by which the spiral tip can interact with the boundary as for this case, the spiral will have not yet completed a rotation. If we sample with equidistant values of t over an interval larger than [tj,^], then because of the periodicity given by A„, will imply that our sampling is not uniform anymore over the interval [£i,£ ], biasing the results in Fig. 4.9. g 2 In F i g . 4.9, we show that for 9 6 [20,50] and 8 G [150,160], every spiral wave is reflected at the boundary. There is thus a range of values 6i G (50,150), for which at least one spiral wave is absorbed at the boundary. For 9i = 120, we find that all the trajectories considered are absorbed by the boundary. t t In order to understand the mechanisms that influence the annihilation and reflection of the spirals, it is necessary to understand the physical mechanisms occurring near the tip of the spiral during the trace of a unit of trajectory. In the next section, we present a study of the propagation of a wavebreak to understand the shape of a spiral tip trajectory; the same study will be helpful to understand what why we obtain annihilation and reflection of a spiral at a boundary. A rationalization of the results presented in Fig. 4.9, is considered in Section 4.6 after discussing the physical mechanisms involved in annihilation and reflections of spiral waves at a boundary. 4.5 The role of excitability in controlling reflection and annihilation at the boundary. The reflection and annihilation of spiral waves at a boundary can be explained by considering the gate variable j G [0,1] i n the B R equations (Appendix A ) . The variable j in the B R model controls the reactivation of the sodium channels, responsible for the initiation of an action potential (AP) [245]. The term m hj controls the opening and closing of the sodium gate, where at rest m « 0, h « 1 and j « 1. A t the beginning of the A P , m approaches 1 on a very fast time scale while h and j remain near 1, such that m / i j « 1. A t this stage the g = gN m hj conductance is maximized, giving a large influx of sodium ions, responsible for the depolarization of the cell. However, after a few milliseconds, the inactivation variables h and j approach zero, so that w?hj « 0, and then the sodium 3 3 z N a a 125 Chapter 4. Annihilation and reflection of spiral waves conductance ~g « 0, terminating the depolarizing current. Another excitation is prevented until the parameters h and j get close enough to 1. The reactivation parameter j sets the time when the medium is ready to accept another A P . This suggests that when the tip of a spiral hits the boundary, the spiral wave will be annihilated if the gate variable has not sufficiently recovered at regions near to the tip of the spiral wave, so as to accept an A P . Na In order to clarify this assertion, we present an analysis of the propagation of an A P in one dimension and the response to an external stimulus. In F i g . 4.10A, we show the propagation of an A P , from left to right, where V (solid line) and j (dashed line) are plotted versus x. In this case, the propagating pulse is generated by applying a stimulus current I , as given in E q . ( B . l ) , with c\ = 100, C2 = 95 and a = 0, for the first millisecond. app Following the study of Glass and Josephson [255] for the propagation of a pulse on a ring, we apply a second stimulus at x = —47.5 (filled circle) after the first propagating pulse has passed the location x = —47.5. The results of applying the stimulus at t* = 280 ms and t* = 240 ms are summarized in Figs. 4.10B,C and 4.10D,E, respectively. The stimulus has the form of E q . ( B . l ) with c i = 50, c = 45 and a = 0. In Figs. 4.10B and 4.10D, we show the plot of V and j after the stimulus has been applied, at times t* = 280 ms and t = 240 ms, respectively. In both cases (Figs. 4.10B and 4.10D), the variables V and j at the location where the stimulus was applied show a peak, product of the stimulation, where V raises up to 0 m V and j approaches to zero. For the case t* = 280 (Fig. 4.10B), the stimulation occurred when j is more than 98% recovered as shown by the value of j inside the small box (Fig. 4.10B) next to the stimulus. Due to the high recovery level of the j gate next to the point where the stimulus is applied, two new A P s are generated at the stimulus location (Fig. 4.10C), one that propagates in the direction of the original A P and a second one going in the opposite direction, as indicated by the bold arrows i n F i g . 4.10C. A different result is obtained when the stimulus is applied at t = 240 ms (Figs. 4.10D,E). From F i g . 4.10D, the stimulation occurred when the j variable is recovered within the range [75%, 95%], as shown by the values of j inside the small box (Fig. 4.10D) next to the stimulus location. In this case, one A P in the opposite direction to the original pulse is generated, as indicated by the bold arrow in Fig. 4.10E. However, an A P could not be generated in the same direction as the original pulse (Fig. 4.10E). The reason is that at the region on the left of the location stimulus, j is recovered above 95% and an A P propagates, whereas on the right of the location stimulus, the gate j 2 126 Chapter 4. Annihilation and reflection of spiral waves is less than 75% recovered, blocking the propagation of an A P . The principle that an A P cannot be generated if j is not completely recovered is applied to study spiral waves and their interaction with a boundary. Prior to the analysis of the effects of the tip trajectory near the boundary, a discussion of the evolution of a wavebreak, as a function of the excitability of the medium, similar as the one in [251], is presented. A n analogy of the propagation of a wavebreak and the front near the tip is considered afterwards. In F i g . 4.11, a wavebreak [251], which consists of a plane wave with a free end, is shown. The initial profile of the wavebreak is shown at the left side on Fig. 4.11 (Panel A ) . The direction of propagation of the front is from left to right as shown by the arrows. The shaded area behind the propagating plane front indicates the region that is in excited state. The point 'c' indicates the position where the wavebreak is located. Each point in the propagating front far from point 'c', influences the excitation of places that are completely recovered, located to the right of the front, within a distance of p/2 as shown in F i g . 4.11 (Panel A ) , by the semicircle with center at the point 'a'. In the same way, a point 'b' to be excited by the front, is affected by the points located at the propagating front that are a distance less than p/2 of point 'b'. Then, each point at the front far from the broken end (point c) propagates at the same speed, giving a uniform advance of the front. A t the broken end (Point 'c' i n F i g . 4.11, Panel A ) , the points near 'c' (within a distance of p/2 at the front) have to excite a larger area (shaded region), than points at the front far from 'c'. Each of the panels B , C and D in Fig. 4.11, shows the evolution in time of the propagating wavebreak in Panel A at two different times, t\ and t , for values of g^a equal to 1.9 m S c m , 2.0 m S c m and 2.1 m S c m , respectively. As discussed in [251], the excitability of the medium, which is the inverse of the threshold of excitation, increases as the value of g^ increases. The case in panel B , where g^ = 1-9 m S c m , represents the low excitability regime. In this case, the depolarizing current generated at a neighborhood of point 'c' is not large enough to excite the p o i n t ' d ' , and by logic, any other part that lies at the bottom part of the three quarters of circle. Therefore the propagating front shrinks (Panel B) until it disappear at the boundary. In the critical case, g^ = 2.0 m S c m , the depolarizing current at a neighborhood of 'c' in panel A , can excite point 'd', but not any other point below 2 - 2 - 2 - 2 a - 2 a - 2 a 127 Chapter 4. Annihilation and reflection of spiral waves 'd', giving as a result the propagation of a wavebreak that neither shrinks nor grows as shown in panel C, with the two consecutive fronts. Finally, in Fig. 4.11 (panel D ) , the case of high excitability (g^ = 2.1 m S c m ) is considered. In this case, the depolarizing current at a neighborhood of 'c' is able to excite p o i n t ' d ' plus, part of the shaded area below p o i n t ' d ' , giving as a result a developing of a spiral wave (panel D ) . - 2 a The analysis of Fig. 4.11 for the propagation of a wave in a medium, where the conductance provides a measure of excitability, can be thought of as the propagation of a wave in regions where the activation variable j in the B R model is in different stages of recovery. In regions where the sodium channels are in the process of reactivation (0 < j < 0.7), can be seen as regions of low excitability, whereas in regions where the sodium channels are activated (j > 0.96) can be seen as regions of high excitability. In Figs. 4.12 and 4.13, the contour plots of the recovery variable j are shown. The area in white is the part where the j gate is at least 85% recovered and can accept propagation of an action potential. The colored areas denote different stages of j in the process of recovery, where the widest colored area denotes a recovery of at most 20%, which means that the gate j is closed. The narrower bands are intermediate stages of between 20% and 85% and their range is 20% to 40%, 40% to 70% and 70% to 85%. The black dot filled in white denotes the location of the tip of the spiral at a given time. In Figs. 4.12 and 4.13, 9i = 80 and the value of t is equal to 68 ms and 56 ms, respectively. g In F i g . 4.12A, the tip of the spiral is located at a petal. In this case, the front propagates over a region that is almost completely recovered facilitating the propagation and therefore, the tip follows a pattern with large curvature. This is analogous to the case of high excitability in the medium for the wavebreak propagation (Fig. 4.11). The tip of the spiral can be seen as the point 'c' in F i g . 4.11 (Panel A ) , propagating over a highly excitable region. In Fig. 4.12B, the tip now is located at an arc. In this case, notice that the front propagates through a medium that is not completely recovered. As a consequence, the large curvature of the tip trajectory cannot be executed by the spiral. This is analogous to the case where the wavebreak propagates on a medium low or moderate excitable. 128 Chapter 4. Annihilation and reflection of spiral waves In F i g . 4.12C, the tip of the spiral gets really close to the boundary of the medium and the tip trajectory gets deformed from its original periodic shape. Again, we consider the case where the wavefront near the tip behaves like the wavebreak near the point 'c' in F i g . 4.11 (Panel A ) . Very close to the boundary, the region given by the shaded area in F i g . 4.11 (Panel A ) gets shortened as part of the three quarters of the circle lies outside the computational domain. Therefore, for this smaller region, the amount of current used for activation is the same as if the three quartered circled region were still complete. Therefore the chance of activating a region near the boundary increases. This explains the deformation of the spiral tip trajectory in Fig. 4.12C. Moreover, the front next to the tip of the spiral continues to propagate as the region ahead of it is now recovered and then reflection at the boundary is observed as shown in Fig. 4.12D. In F i g . 4.12C, when the region near the tip, that in the absence of the boundary cannot be excited, gets excited, the unusual front propagating in this region encounters the waveback of the previous excitation of the spiral. Then, the propagation of the new front near the tip takes place on a region that has low excitability, just like being in an axe of the unit of a trajectory (Fig. 4.12D). Under this conditions, the angle of the line at which the petals lie, is changed (Fig. 4.12D). The phenomenon of annihilation is summarized in Fig. 4.13. In F i g . 4.13A, we show again that when there is a region almost or completely recovered ahead of the front, the tip of the spiral wave traces a long curvature curve generating a petal. In F i g . 4.13B, the tip of the spiral has almost arrived at the boundary of the domain and is propagating over a region that is not recovered. Although there is more current available to stimulate the region near the boundary (part of the shaded three quarters of circle in F i g . 4.11 falls outside the domain), it is not sufficient to excite what has not been able to recover. Therefore the tip of the spiral wave leaves the domain as shown in F i g . 4.13C. Finally, in F i g . 4.13D the spiral is not sustained anymore as the tip has disappeared at the boundary giving an end to the the spiral wave motion. This analysis is consistent with the result presented by Yermakova and Pertsov [276] for the interaction of a circular spiral tip trajectory with a boundary. They observed that the curvature of the trajectory (a circle) increased near the boundary. This is consistent with the idea whereby the spiral tip trajectory near the boundary behaves as a wavebreak as shown in 129 Chapter 4. Annihilation and reflection of spiral waves Fig. 4.11 (Panel A ) , where part of the shaded region lies now outside the integration domain. The immediate consequence as discussed throughout this section is an increase in the curvature in the tip trajectory. 4.6 A rationalization of the fraction of trajectories annihilated, F A ( ^ ) , Fig. 4.9. Figure 4.9, shows the variation of the fraction of trajectories annihilated versus the incident angle di. It was shown in the previous section that annihilation was due to the incomplete recovery of the j variable in regions where the front has to propagate in order to stay inside the domain as shown in Fig. 4.13. From this observation, we suggest that when annihilation occurs, the intersection of the tip trajectory with the boundary takes place at an arc of a unit of a spiral. This is anticipated, as the propagating region near the tip is in its less recovered regime when the tip is located at an arc, making more difficult for the front to propagate into such a region. Figure 4.14A, shows a series of trajectories for different angles for which annihilation is observed. Notice that in all cases, the trajectory hits the boundary with an arc. From F i g . 4.14A, the trajectory with 9i = 140° has a supplementary incident angle of 8 = 40°. When 0i is larger than 150°, the trajectory, hits the boundary with the arc portion of the spiral unit and the supplementary incident angle between is small (6 < 30°). From the discussion at the end of the previous section, we know that the tip trajectory gains curvature near the boundary. Therefore, in this case the boundary effects does not allow the trajectory to leave the domain. S S When 6i = 140°, reflection and annihilation are observed (Figs. 4.14A and 4.14B). From F i g . 4.14A, the tip hits the boundary when it is located at a very early stage of an arc (bold curve), giving annihilation. For the same angle, the tip hits the boundary when it is located at a later stage i n the arc (bold), giving reflection of the tip (Fig. 4.14B). From the previous section, we know that at an arc the excitability of the medium is low whereas at a petal it is high. However, Efimov et. al. [251] show that the threshold of excitation takes its largest value at the beginning of an arc and gradually decreases until it reaches a minimum value at a petal, increasing the chance of excitation. From the same paper [251], we know that an increase of the threshold of excitation is equivalent to a decrease of excitability. Therefore, 130 Chapter 4. Annihilation and refection of spiral waves boundary effects at later stages of the axe in a unit of trajectory are more effective bending the tip trajectory than at the beginning of the arc, where excitability has its lowest value. Figure 4.14B, where different examples of trajectories when reflection occurs, shows that the intersection of the tip trajectory and the boundary, occurs at a later stage of the arc or at a petal. The same phenomenon occurs for € [90°, 150°]. In the case 9i = 120° (Fig. 4.7), annihilation of the spiral wave occurs regardless of the stage of the arc at the moment of hitting the boundary. In Fig. 4.7A, annihilation of the spiral occurs because the two trajectories hit the boundary at early stages of an arc, so that bending due to boundary effects is not strong enough to change the direction of the trajectory. In Fig. 4.7B, for t = 122 ms, the trajectory turns to the left almost parallel to the boundary due to boundary effects. A t this stage, the medium where the wave propagates is now ready to be activated and a petal is formed. A t the end of the petal (end of bold line in Fig. 4.7B), the region where the front needs to propagate to remain inside the domain, is now at its maximum level of refractoriness (lowest excitability) as a new petal has just formed, leading to the annihilation of the spiral. The same phenomenon is observed in Figs. 4.7C and 4.7D, for t = 129 and 136 ms, respectively. g g When 8i < 50°, no annihilation of the spiral tip is observed. A representative example is given by the trajectory with 6\ = 40° in F i g . 4.14B. When 9 < 50°, the spiral tip generally touches the boundary when it is located at a petal portion or at a very late stage of the arc portion of the spiral unit of trajectory. When the tip moves along an arc portion, because neither the incidence angle nor the length of the arc are too large, the trajectory does not have the time to reach the boundary (in the arc portion), when a petal is suddenly formed (0; = 40° in Fig. 4.14B). When the tip, located at the petal portion of a spiral unit, hits the boundary, reflection of the tip is observed as discussed in Section 4.5 for Fig. 4.12. Finally for 9i e [50°, 90°], the probability of the spiral being annihilated at the boundary increases as Ot increases (Fig. 4.9). When 0; is close to 50° and annihilation is observed, the tip hits the boundary with the later stage of an arc due to the presence of the petals which in this case face the boundary (Fig. 4.14A). Near 9{ = 90°, the tip can hit the boundary at an early and a later arc portion of a spiral unit. Then, there is a higher probability for the tip to hit the boundary, when the tip is located at an arc, for 9{ close to 90° compared to 9i close to 50°. Because the probability of annihilation in131 Chapter 4. Annihilation and reflection of spiral waves . . . creases when the tip is located at an arc, it gives the increasing dependence of F (0i) on Oi for 6>; e [50°, 90°] in Fig. 4.9. A Therefore, we can summarize in a general way when annihilation is observed. For incident angles smaller than 80° a trajectory has larger probability of hitting the boundary with a petal or with the later stage of the arc, than for angles larger than 100°. W i t h only this argument, the probability of annihilation is lower for Bi < 80° than for Bi > 100°. However, this is not the case for angles larger than Bi = 150° (Fig. 4.9). In this case, the trajectories always hit the boundary with an arc, but as discussed in previous paragraphs, boundary effects are able to keep the tip inside the domain. 4.7 Summary. In this chapter, the phenomenon of annihilation and reflection of a spiral wave at a boundary for the Beeler-Reuter model, was considered. The Roo case, in which the petals of the tip trajectory lie on a straight line, was studied numerically in detail. The case where a spiral tip traces a trochoid of finite R (Fig. 4.1), is very complex to study compared to the R^ case, as can be clearly seen in Figs. 4.3 and 4.4. However, the tip trajectory near the boundary for the case with finite R, where reflection and annihilation take place, can be approximated by the tip trajectory given by the R^ case. Therefore, the results obtained with the Roo case can be applied to understand annihilation and reflection when R is finite. In Section 4.4, by considering trajectories in the Roo case, it was shown that annihilation of the spiral depends strongly on the angle of incidence Bi with respect to the boundary. For angles above Bi = 150° and below 50°, no annihilation was observed. For Bi £ [50°, 150°], both, annihilation and reflection of the tip of the spiral was observed (Fig. 4.9). In order to understand the results in F i g . 4.9, an analysis based on the recovery gate j in the B R model was considered. This variable is responsible for the local reactivation of the Na channels, and a lack of recovery will forbid the activation of an A P at such location. Therefore, in Section 4.5 we present an analysis to understand the role of the j variable in accepting an A P . In the same section, it was considered the effects of the impermeable boundary on the propagation of A P in two dimensions. This analysis was 132 Chapter 4. Annihilation and reflection of spiral waves based in the understanding of propagation of a wavebreak in a medium with different degrees of excitability. This section (4.5), ends with a qualitative explanation for annihilation and reflection of a spiral at the impermeable boundary. Based on the arguments presented in Section 4.5, a qualitative explanation of Fig. 4.9 was shown in Section 4.6. In this section, we discussed how the different angles of incidence of the spiral tip trajectories affect the way the tip interacts with a boundary. A quantitative understanding of this phenomenon is a very difficult task. This is mainly due to the effects of the impermeable boundary on the spiral tip trajectory. The analysis presented in this chapter bridges a couple of observations about simple rotation. In the case of simple rotation, when the tip hits the boundary the spiral wave annihilates as shown by Krinsky et al [264]. However, the studies presented by Yermakova and Pertsov [276], show that if the circular trajectory is close enough to the boundary an increase in the curvature is observed. From the studies presented in this chapter, the results in [264] and [276] correspond to annihilation and reflection of the wave, respectively, with the corresponding arguments as discussed in Section 4.5. The studies presented in this chapter play a very important role when a spiral wave interacts with an anatomic obstacle [260, 267]. It has been observed experimentally that a spiral wave, in the presence of an obstacle of some minimum size, anchors to the obstacle [260, 267]. In the same way, it has been observed that annihilation of a spiral wave occurs at a boundary of isolated tissue [267]. The results presented in this chapter, extend these ideas to the possibility of the avoidance of attachment of the spiral to the obstacle and the annihilation of the spiral at the boundary. 133 Chapter 4. Annihilation and reflection of spiral waves . . . Table 4.1: Measure of the total length Lo in millimeters, and period of rotation TO in milliseconds of a unit of trajectory for three cases shown in Figs. 4.6(B-D) Notice that for case D corresponding to F i g . 4.6D, there is no information available for the dashed unit of trajectory, as the trajectory was annihilated at the boundary. B First two Bold Dashed C D Lo TO Lo TO Lo TO 24.8 16.46 22.38 120 78 110.4 24.8 25.2 26.58 120 121.7 133.14 24.8 22.5 120 108.7 - - 134 Chapter 4. Annihilation and reflection of spiral waves ... • B I A " rx ( ' ' \ 1 i j 'T\ / X \ - \ \ ^ JC^ i \ ' * if ^ ' / i 1 1 / / / ' \ \ I V 1 / X ^ \ "3y. •^Jr" I —"-w 1 1 1 ^ \ 1 ) / ^ ^ s / i i- I Figure 4.1: Compound rotation or flower like trajectories, resemble two different curves (A) Hypotrochoid for outward petal flowers (B) Epitrochoid for inward petal flowers. A hypotrochoid is the trajectory traced by the arrow at the end of a small line segment of length r attached to the center of a small circle of radius h, that rotates on the inner side of a circle of radius R. A n epitrochoid is obtained in the same way, but the rotation of the circle of radius h takes place on the outer side of the circle of radius R. 135 Chapter 4. Annihilation and reflection of spiral waves 2.18 2.38 2.6 2.8 9Na Figure 4.2: Spiral tip trajectories obtained with the B R model with respect to two parameters, g^ and E^a- gs was taken as 0.03 m S c m in all cases. Trajectory E is an inward petal flower, C , D , G , H , I are outward petal flowers whereas B and F represent the i?oo case. The units for the conductance g and E^ are in m S c m and in m V , respectively. The calculations considered N = 332 points in each dimension. - 2 a - 2 s a 136 Chapter 4. Annihilation and reflection of spiral waves (A) (B) (C) (D)" (E)" (F) Figure 4.3: Spiral tip trajectories with the B R model with g^ = 2.38 m S c m and gs = 0.03 m S c m . This is the case. (A) Five tip trajectories. Each trajectory is a a small horizontal translation of the others; (B-F) The same trajectories as in (A) but for a longer integration time. The initial point of each trajectory is marked with (*) and the end of the trajectory is shown with (•). The filled arrow indicates the place at which the trajectory leaves the domain; x and y are in mm; N = 452 points in each dimension. a - 2 - 2 137 Chapter 4. Annihilation and reflection of spiral waves (Cf (D) Figure 4.4: Spiral tip trajectories with the B R model with gs = 0.03 m S c m and ( A , B ) 9Na = 2.41 m S c m ; (C,D) g = 2.34 m S c m . For the pairs ( A , B ) and (C,D), each of the trajectories is a small vertical translation of the other. The initial point of the trajectory is marked with (*), the end of the trajectory is shown with (•). The filled arrow indicates the place at which the trajectory leaves the domain. For ( A , B ) yo = 10 mm and t = 50.5 and 51.5 ms and for (C,D) yo = —10 m m and t = 85 and 83 ms, respectively; x and y are in mm; N = 452 points in each dimension. - 2 - 2 - 2 Na g g 138 Chapter 4. Annihilation and reflection of spiral waves . . . -10 -20 i A i •s. i \ /—H—. S( J i i 1 ^\ -30 -40 -10 •r; i 1 V I < I -20 / / f I \ \ } * , e" r I * - -30 -40 -40 i \ * ^ • -20 1 0 x 1 1 1 20 1 40 Figure 4.5: Quantities associated with reflections at the boundary. Dashed lines indicate the direction of the trajectory of the tip (A) &i and 9 are the incident and reflected angle is 95°. The filled circle is the place where the tip, hits the boundary; (B) A special case of reflection at the boundary. There are two points where a minimum with the boundary is reached (asterisk and filled circle). The filled circle point is the one considered to be the point at which the tip hits the boundary; g^ = 2.37 m S c m , gs = 0.03 m S c m and the rest of the parameters are as in the standard B R model. r - 2 - 2 a 139 Chapter 4. Annihilation and reflection of spiral waves . . . ^-20 Figure 4.6: (A) A unit of a spiral tip trajectory far from the boundary for the Roo case. It consist of a petal and an arc. The total arc length is Lo = 24.8 mm; (B-D) Effects of the boundary on the spiral tip units. The first two units are unaffected by the boundary, whereas the bold and dashed units are the ones affected by the boundary. The filled circle indicates the place where the tip hits the boundary; Parameters as in Fig. 4.5. 140 Chapter 4. Annihilation and reflection of spiral waves -20 20 -10 -10 x (mm) x (mm) Figure 4.7: Spiral interactions with the boundary for an incident angle of Bi = 120°, and different values of t . (A) Trajectories with t = 105 (left) and t = 150 (right); (B-D) T i p trajectories with t = 105 ms (left) and t = 122,129 and 136 ms, respectively. The filled circle in each plot, indicates the position along the trajectory at which the tip hits the boundary. Parameters as i n Fig. 4.5. g g g g (J 141 Chapter 4. Annihilation and reflection of spiral waves ... Or - ~1\ -10 - 1 1 \=&.9 mm - -20 :v=44 A : - c - -30 - _ - -40 Or • # I I -10 - ?B -20 - >.. - -30 -40 - -20 0 x (mm) 40 Figure 4.8: T i p trajectories with = 70°, and different values of t . (A) T i p trajectories with t = 44 (left) and t = 71 ms (right), respectively; (B-D) T i p trajectories with t = 44 ms (left) and t — 52,60 and 67 ms, respectively. The filled circle in each plot, indicates the position along the trajectory at which the tip hits the boundary. Parameters as in F i g . 4.5. g g g g g 142 Chapter 4. Annihilation and reflection of spiral waves Chapter 4. Annihilation and reflection of spiral waves 20h \ \ • 0.5 =5 -20 i. > -60 -100 : & Direction of * propagation & = C O 0 fif -0.5 • j recovered above 98% next to the stimulus location -100 -100 100 -50 0 50 100 50 100 x (C) 100 -50 0 (E) Figure 4.10: (A) Propagating pulse for the B R model in I D . g = 0.03, g = 2.37 m S c m . V (solid) and j (dashed) are plotted versus x at a fixed time. A stimulus is applied at x = —47.5 mm (•) at times (B) t* = 280 and (D) t* = 240 ms; (B) The gate variable j is above 98% recovered at the position at which the stimulus is applied; 144 (C) Response of the stimulus applied at t* = 280. Two new single pulses are generated (bold arrows). (D) The gate variable j has recovered 75 to 85%. (E) Only one A P was generated with the stimulus (bold arrow). The propagation of the new A P in the direction of the original pulse was blocked as j has not recovered completely. s - 2 Na Chapter 4. Annihilation and reflection of spiral waves ... / / / / / / / / / / / / t, / /. / / t tj (B) g =1.9 low excitability Na (A) / / / / / / / / A / / // / / // / / / A / / / t / / / 2 t, 2 ^ ^ 4) / (C) g a=2-0 N moderate excitability (D)g 2.1 high excitability Figure 4.11: (A) A plane front at a fixed time, with a free end (point c) that propagates from left to right. The point 'a' is located at the wavefront and 'b' a n d ' d ' ahead of the front. The semicircle of radius p/2 with center in 'a' is the region of influence of the point 'a'. The shaded region is the region of influence of the point 'c'; (B) Two snapshots of the propagating front when the medium has low excitability g^a = 1-9 m S c m . The front's length shortens; (C) g^ — 2.0 m S c m . Moderate excitability. The length of the front neither shrinks nor grows; (D) The medium is highly excitable gNa = 2.1 m S c m , and then a spiral wave is formed from the free end. Values of g^ taken from [251]. - 2 - 2 a - 2 a 145 Chapter 4. Annihilation and reflection of spiral waves Figure 4.12: Reflection of a spiral. 0< = 80, t = 68. The contour plots shown for different integration times, represent different stages of recovery of the j variable. The black dot filled in white is the location of the tip of the spiral for that specific time. g 146 Chapter 4. Annihilation -40 -40 |\I -20 20 40 and reflection of spiral waves . . . 0 -20 0 A (A) x x 0 20 0 40 Figure 4.13: Annihilation of a spiral. 6i = 80, t = 56. The contour plots shown for different integration times, represent different stages of recovery of the j variable. The black dot filled in white is the location of the tip of the spiral for that specific time. In this case, the tip of the spiral has left the domain on C . g 147 Chapter 4. Annihilation and reflection of spiral waves . . . A =6O 0i 0 e 80 e-ioo i= 0 B -10 -20 " * a=40 / 110 AO D -20 1 4 0 ^ a = 1 2 0 Q\a=70 - e.=i40cl sr\ -30 -40 -40 i = 0 20 40 Figure 4.14: Particular examples of trajectories of spirals that (A) annihilate at the boundary. 8i = 60° at the extreme left and for the others with increments in 6i by 10°; (B) Are reflected at the boundary. 148 Bibliography [243] D . Barkley, Spiral meandering, Chemical and wave patterns, (Kluwer academic publishers, Netherelands, 1995). [244] J . Beaumont, N . Davidenko, J . M . Davidenko, J . Jalife, Spiral waves in two-dimensional models of ventricular muscle: formation of a stationary core, Biophys. J . 75 (1998) 1-14. [245] G . W . Beeler, H . Reuter, Reconstruction of the action potential of ventriculal myocardial fibres, J . Physiol. 268 (1977) 177-210. [246] V . N . Biktashev, Control of reentrant vortices by electrical stimulation, Computational Biology of the heart, (John Wiley & Sons, Chichester, 1997). [247] T . R . Chay, Proarrhythmic and antiarrhythmic actions of ion channel blockers on arrhythmias in the heart: model study, A m . J . Physiol. 271 (1996) 329-356. [248] P. Comtois, J . Kneller, S. Nattel, Of circles and spirals: Bridging the gap between the leading circle and spiral wave concepts of cardiac reentry, Europace 7 (2005) S10-S20. [249] J . Davidenko, A . V . Pertsov, R . Salomonsz, W . Baxter, J . Jalife, Stationary and drifting spiral waves of excitation in isolated cardiac muscle, Nature, 355 (1992) 349-351. [250] L . Edelstein-Keshet, Mathematical models in biology, Classics in A p p . Math., ( S I A M , New York, 1988) [251] I. R . Efimov, V . I. Krinsky, J . Jalife, Dynamics of rotating vortices in the Beeler-Reuter model of cardiac tissue, Chaos, Sol. & Fract. 5 (1995) 513-526. [252] F . Fenton, E . Cherry, H . M . Hastings, S. J . Evans, Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity, Chaos 12 (2002) 852-892. 149 Bibliography [253] F . Fenton, E . M . Cherry, A . Karma, W . J . Rappel, Modeling wave propagation in realistic heart geometries using the phase-field method, Chaos 15 (2005) 013502(1-11). [254] R . Fitzhugh, Impulses and physiological states i n theoretical models of nerve membrane, Biophy. J . 1 (1961) 445-466. [255] L . Glass, M . E . Josephson, Resetting and annihilation of reentrant abnormally rapid heartbeat, Phys. Rev. Lett. 75 (1995) 2059-2062. [256] M . Gomez-Gesteira, A . P . Munuzuri, V . Perez-Munuzuri, V . PerezVillar, Boundary imposed spiral drift, Phys. Rev. E , 53 (1996) 54805483. [257] R . A . Gray, J . Jalife, Ventricular fibrillation and atrial fibrillation are two different beasts, Chaos 8 (1998) 65-77. [258] S. G r i l l , V . S. Zykov, S. C . Miiller, Spiral wave dynamics under pulsatory modulation of excitability, J . Phys. Chem. 100 (1996) 1908219088. [259] T . Ikeda, T . Uchida, D . Hough, J . J . Lee, M . C . Fishbein, W . J . Mandel, P. S. Chen, S. Karagueuzian, Mechanisms of spontaneous termination of functional reentry in isolated canine atrium, Circulation 94 (1996) 1962-1973. [260] T . Ikeda, M . Yashima, T . Uchida, D . Hough, M . C . Fishbein, W . J . Mandel, P. S. Chen, H . Karagueuzian, Attachment of meandering reentrant wavefronts to anatomic obstacles in the atrium, Circ. Res. 81 (1997) 753-764. [261] A . Karma, Meandering transition i n two dimensional excitable media, Phys. Rev. Lett. 65 (1990) 2824-2827. [262] J . P. Keener, J . Sneyd, Mathemathical Physiology, Interdisciplinary Applied Mathematics (Springer, New York, 1998). [263] A . G . Kleber, Y . Rudy, Basic mechanisms of cardiac impulse propagation and associated arrhythmias, Physiol Rev, 84 (2004) 431-488. [264] V . I. Krinsky, V . N . Viktashev, A . M . Pertsov, Autowave approaches to cessation of reentrant arrhythmias, A n n . New York Academy Sci. 591 (1990) 232-246. 150 Bibliography [265] D . Olmos, B . Shizgal, A Spectral Method of Solution of Fisher's Equation, J . Comp. Appl. Math. 193 (2006) 219-242. [266] D . Olmos, B . Shizgal, Pseudospectral method of solution of the Fitzhugh-Nagumo equations, (submitted) [267] A . M . Pertsov, J . M . Davidenko, R. Salomonsz, W . T. Baxter, J . Jalife, Spiral waves of excitation underlie reentrant activity in isolated cardiac muscle, Circ. Res. 72 (1993) 631-650. [268] Z. Qu, A . Garfinkel, A n advanced algorithm for solving partial differential Equation i n cardiac conduction, I E E E Trans. Biomed. Eng. 46 (1999) 1166-1168. [269] J . A . Trangenstein, C . K i m , Operator splitting and adaptive mesh refinement for the Luo-Rudy I model, J . Comp. Phys. 196 (2004) 645-679. [270] J . Tyson, What everyone should know about the BelousovZhabotinsky reaction, Lecture notes in biomathematics, Springer, Vol. 100 (1994) 569-587. [271] A . T . Winfree, When times breaks down, (Princeton University Press, Princeton, 1987), Pages 181-183. [272] A . T . Winfree, Electrical instability in cardiac muscle: phase singularities and rotors 138 (1989) 353-405. [273] A . T . Winfree, Varieties of spiral wave behavior: A n experimentalist approach to the theory of excitable media, Chaos 1 (1991) 303-334. [274] A . T . Winfree, Evolving perspectives during 12 years of electrical turbulence (Focus issue), Chaos 8 (1998) 1-19. [275] F . X i e , Z. Qu, J . N . Weiss, A . Garfinkel, Interactions between spiral waves with different frequencies in cardiac tissue, Phys. Rev. E . 59 (1999) 2203-2205. [276] Ye. A . Yermakova, A . M . Pertsov, Interaction of rotating spiral waves with a boundary, Biophys., 31 (1986) 932-940. 151 Chapter 5 Conclusions One of the objectives of the present thesis was to develop accurate and efficient pseudospectral solutions of reaction diffusion equations with excitable dynamics (Eq. 1.13). The problem with this type of equation is that the solution develops shock-like wave behaviour giving rise to a multiple spatial and temporal scales problem [283]. The numerical schemes in this thesis were based on a collocation method with Chebyshev-Gauss-Lobatto quadrature points. In order to develop accurate numerical methods for these type of equations, we considered in Chapter 2 and 3 the Fisher's and the FitzhughNagumo ( F H N ) equations, respectively. In Chapter 2, we studied the Fisher's equation, a prototypical reaction diffusion equation. The collocation method used the Chebyshev-Gauss-Lobatto quadrature points. The solutions of F E are characterized by propagating fronts that can be steep depending on the value of the reaction rate coefficient, p (Fig. 2.6). In order to obtain an accurate solution, the main domain was subdivided into smaller subintervals as proposed by Shizgal and coworkers [293, 302, 303]. The multidomain method provided stable and accurate solutions of F E for values of p as large as 10 (Fig. 2.6). We compared the present numerical treatment with the D S C approach of Zhao and Wei [306] who employed an interpolation based on the sine function. For the value p = 10 , the accuracy obtained with the multidomain method was better than the one reported by Zhao and Wei (Chapter 2). 6 6 Due to the excellent performance of the Chebyshev multidomain method in the solution of the Fisher's equation, we considered this scheme for the solution of differential equations of the F H N type (Chapter 3). The solutions of the F H N equations provide the simplest P D E to model excitable media (Eq. 3.2), and has been used extensively for studies of spiral waves [278, 288, 289, 290, 294, 300]. The solutions of the F H N equations present propagation of shock like waves as in Fisher's equation. In Chapter 3, we presented a convergence study for reaction diffusion equa152 Chapter 5. Conclusions. tions for the F H N model, where we compared a Chebyshev multidomain method with finite difference based methods, especially the method proposed by Barkley [278]. Simulations of plane waves in I D and spiral waves in 2D were carried out for two F H N equations with different local dynamics. The first equation, named kinetic model I, was proposed by Barkley [278], whereas the second equation, kinetic model II, is the classic cubic F H N model discussed by Keener [289]. In the same chapter, we demonstrated the superiority of the Chebyshev multidomain approach with regard to accuracy and computational time in comparison with finite difference and the Barkley method. The method developed by Barkley applies to a particular type of local dynamics, E q . (3.3), and can be very fast but the results are compromised due to the decreased accuracy. The Chebyshev multidomain approach can be used to solve equations of the F H N type with more general local dynamics than the Barkley method. For the solution of the F H N equations, we implemented an operator splitting method for the Chebyshev multidomain approach for kinetic model I and II. W i t h the operator splitting scheme, the computational time for the solution of the F H N equations is considerably shorter than both the finite difference and the Barkley method. Therefore, in Chapters 2 and 3, we developed numerical schemes based on pseudospectral methods for reaction-diffusion equations with excitable local dynamics. The methods developed with pseudospectral methods were superior in accuracy with respect the number of points used, and the computational time over conventional finite difference numerical methods. The solution of reaction-diffusion equations with spectral methods has been considered by a number of people in the field [281, 288, 305]. Bueno et. al. [281] considered a spectral method for irregular domains to solve reactiondiffusion problems, K a r m a [288] considered the classical Fourier method for the solution of the F H N equations and Zhang and N g [305] considered the Chebyshev pseudospectral method to analyze convergence of the solution of the Luo-Rudy equations for a very short integration time. The Chebyshev multidomain method presented in this thesis represents another forward step in the use of spectral methods for solving equations of reaction-diffusion type with excitable dynamics. A n important concern when solving shock-like waves in equations like the 153 Chapter 5. Conclusions. Fisher's and the F H N equations is the phenomenon referred to as Gibbs oscillations [280, 282, 302]. This type of oscillation arises in spectral approximations of functions with jump discontinuities at the location of the jump [280, 282]. In the numerical solutions for the Fisher's and F H N equations at the locations where the propagating fronts are located, Gibbs oscillations did not play a significant role in the numerical solutions. As discussed in [302], the use of the Chebyshev multidomain approach for functions with rapid variations in the spatial variable helps to avoid changes in two scales within each sub-domain, providing better convergence results. The numerical results obtained for the Fisher's and F H N equations extends the use of spectral methods to the case of non smooth problems of the reaction diffusion type, adding to the use of spectral methods in other areas where propagation of shock waves is present, such as in fluid dynamics [282, 284, 302]. The excellent results obtained for the solution of the Fisher's and F H N equations provided confidence for solving the Beeler-Reuter equations, a more realistic model of reaction-diffusion with excitable dynamics. The accuracy of the numerical method considered in Chapter 4 was tested in a I D travelling pulse. In this case, the solutions with N = 542 collocation points provided a very good approximation to the converged solution. The second objective of this thesis was the study of the phenomenon of annihilation and reflection of a spiral wave in the meandering regime at a boundary, for the Beeler-Reuter model. The most practical way to study this phenomenon was to consider the limiting i ? o o case, in which the petals of the tip trajectory lie on a straight line. The case where a spiral tip traces a trochoid of finite R (Fig. 4.1) is very complex to study compared to the limiting Roo case (See Figs. 4.3 and 4.4). However, the tip trajectory near the boundary for the case with finite R, where reflection and annihilation take place, can be approximated by the tip trajectory given by the limiting Roo case. Therefore, the results obtained with the limiting Roo case can be applied to understand annihilation and reflection when R is finite. In Chapter 4, we have shown that annihilation of a spiral wave was due to the lack of recovery of the medium to accept propagation of action potentials in regions where the front of the spiral wave near its tip needs to propagate in order to remain in the domain. In the Beeler-Reuter model, the gate variable j , which is called the sodium reactivation variable, controls the time at which the sodium channels are ready to accept another excitation. Therefore, the analysis for the annihilation and reflection of spiral waves at 154 Chapter 5. Conclusions. a boundary was explained in terms of this gate variable. B y considering trajectories in the limiting case, we showed that annihilation of the spiral depends strongly on the angle of incidence di with respect to the boundary. For angles above 6i = 150° and below 50°, no annihilation was observed. For di G [50°, 150°], both, annihilation and reflection of the tip of the spiral with exception of &i = 120°, was observed (Fig. 4.9). A n analysis of annihilation and reflection and its dependence of the incident angle 9i, as a function of the reactivation sodium channel gate j, is provided. The study presented in Chapter 4 is an extension of the work by Yermakova and Pertsov [304], where they analyzed the effects of the boundary on the trajectory of a tip in the simple rotation regime. Moreover, the analysis presented in this chapter bridges a couple of observations about simple rotation. In the case of simple rotation, when the tip hits the boundary the spiral wave annihilates as shown by Krinsky et al [292]. However, the studies presented by Yermakova and Pertsov [304], show that if the circular trajectory is close enough to the boundary an increase in the curvature is observed. From the studies presented in this chapter, the results in [292] and [304] correspond to annihilation and reflection of the wave, respectively, with the corresponding arguments as discussed in Chapter 4. The situations studied in this chapter play very important roles when a spiral wave interacts with an anatomic obstacle [287, 296]. It has been observed experimentally that a spiral wave, in the presence of an obstacle of some minimum size, anchors to the obstacle [287, 296]. In the same way, it has been observed that annihilation of a spiral wave occurs at a boundary of isolated tissue [296]. The results presented i n this chapter extend these ideas to the possibility of the avoidance of attachment of the spiral to the obstacle and the annihilation of the spiral at the boundary. 5.1 5.1.1 Future directions. The role of obstacles and application to cardiac arrhythmias. In excitable systems such as cardiac tissue, propagation does not occur in an homogeneous media. The presence of small vessels inside the myocardium provide regions in the cardiac tissue with different degrees of excitability 155 Chapter 5. Conclusions. [295]. Also, the presence of natural obstacles such as the inferior vena cava or pulmonary veins in the atria, may interfere in the normal propagation of A P s [277, 285] or in the meandering of spiral waves [287, 296]. These heterogeneities can generate spiral waves [295], can destabilize them [299], or if the heterogeneity is given by an obstacle of some minimum size, a meandering spiral might attach to the obstacle giving a more controlled and predictable behaviour in the spiral dynamics. [287]. A n extension of the studies we presented in this thesis about reflection and annihilation of spirals at a boundary, lead to different questions such as: When a spiral interacts with an obstacle of determined size, what is the role played by the shape and size of the obstacle in order to observe attachment? What is the role played by the radius R of the meandering trajectory compared to the size of the obstacle, to assure attachment of the spiral? 5.1.2 Quantitative factors explaining the annihilation-reflection phenomenon. A different issue arising from Chapter 4 consists of providing a more quantitative analysis of the factors that give annihilation and reflection of a spiral wave at the boundary. In Chapter 4, we could only give a qualitative explanation of annihilation-reflection, based on two parameters, the angle of incidence 9i and the level of recovery of the region where the tip trajectory needs to propagate. However, measured quantities that allow us to know with certainty when annihilation or reflection occurs have not been possible to obtain at this stage. Further studies in this area will be considered. It is important to mention that no previous studies in this area have been reported in the literature. A third extension of this thesis arises from the results shown in F i g . 4.3. Due to the sensitivity of initial conditions we ask: It is possible to determine whether or not the trajectory of the spiral tip follows a chaotic pattern? Is there some measure of the time that can tell us with some certainty that the spiral will be annihilated? Under which conditions? These questions are very difficult to answer now, but might be more doable when quantitative parameters are found to explain in a more precise way the phenomenon of annihilation-reflection of a spiral wave. 156 Chapter 5. Conclusions. 5.1.3 The Bidomain model. As discussed in Chapter 1, one of the simplest ways to study spatial propagation of waves is by considering a R D equation type as i n E q . (1.13). However, wave propagation in cardiac tissue has some factors to be considered. This is the case of anisotropic conductivity i n the intracellular and extracellular spaces of myocardial cells. The difference between the conductance properties in each of the two spaces leads to consider a generalization of E q . (1.13) studied in this thesis, that keeps track of the currents in and and between the extra and intracellular spaces. W i t h these new considerations, equation (1.13) is generalized as V • (cxiVVi + <r VV ) e e = 0 (5.1) (5.2) where o~i and a are conductivity tensors, x is the membrane surface to volume ratio, Vi and V are the intracellular and extracellular potential respectively, and lion depends on the ionic model considered. This set of equations is referred to as the bidomain model [286]. e e The bidomain model has gained acceptance among elecrocardiologists and bioengineers modelling cardiac electrical activity [291] and has been considered by different authors [291, 297, 298, 301]. In [297], Roth describes the new phenomena found by Winfree, for cardiac wave propagation, by considering the bidomain equations, whereas in [291, 301], numerical schemes to solve such equations were proposed. Therefore, i n the ongoing effort of developing spectral methods for equations with excitable local dynamics, I will adapt the Chebyshev multidomain method developed in this thesis to solve Eqs. 5.1 and 5.2 in two dimensions. One of the main difficulties, other than the multiple scales problem, arises with the conductivity tensors and a , which might be space dependent. e 157 Bibliography [277] E . M . Azene, N . A . Trayanova, E . Warman, Wavefront-Obstacle interactions in cardiac tissue: A computational study, A n n . Biomed. Eng., 29 (2001) 35-46. [278] D . Barkley, A model for fast computer simulation of excitable media, Physica D , 49 (1991) 61-70. [279] G . W . Beeler and H . Reuter, Reconstruction of the action potential of ventricular myocardial fibres, J . Physiol. 268 (1977) 177-210. [280] J . P. Boyd, Chebyshev and Fourier Spectral Methods (Dover, New York, 2000). [281] A . Bueno Orovio, V . Perez Garcia, F . Fenton, Spectral methods for partial differential equations in irregular domains: The spectral smoothed boundary method, S I A M J . Sci. Comp. 28 (2006) 886-900. [282] C . Canuto, M . Y . Hussaini, A . Quarteroni and T . A . Zang,'Spectral methods in fluid dynamics, Springer Series in Computational Physics, (Springer, New York, 1988). [283] E . M . Cherry, H . S. Greenside, C . S. Henriquez, A space-time adaptive method for simulating complex cardiac dynamics, Phys. Rev. Lett. 84 (2000) 1343-1346. [284] W.-S. Don, D . Gottlieb, J . - H . Jung, A multidomain spectral method for supersonic reactive flows, J . Comp. Phys. 192 (2003) 325-354. [285] R . A . Gray, J . Jalife, Ventricular fibrillation and atrial fibrillation are two different beasts, Chaos 8 (1998) 65-77. [286] C . S. Henriquez, Simulating the elctrical behavior of cardiac tissue using the bidomain model, Crit. Rev. Biomed. Eng. 21 (1993) 77. [287] T . Ikeda, M . Yashima, T . Uchida, D . Hough, M . C . Fishbein, W . J . Mandel, P. S. Chen, H . Karagueuzian, Attachment of meandering 158 Bibliography reentrant wavefronts to anatomic obstacles in the atrium, Circ. Res. 81 (1997) 753-764. [288] A . Karma, Meandering transition in two-dimensional excitable media, Phys. Rev. Lett., 65 (1990) 2824-2828. [289] J . P. Keener, Mathemathical Physiology, Interdisciplinary Applied Mathematics (Springer, New York, 1998). [290] J . P. Keener, A geometrical theory for spiral waves in excitable media, S I A M J . Appl. Math. 46 (1986) 1039-1056. [291] J . P. Keener, K . Bogar, A numerical method for the solution of the bidomain equations in cardiac tissue, Chaos 8 (1998) 234-241. [292] V . I. Krinsky, V . N . Viktashev, A . M . Pertsov, Autowave approaches to cessation of reentrant arrhythmias, A n n . New York Academy Sci. 591 (1990) 232-246. [293] G . Mansell, W . Merryfield, B . Shizgal, U . Weinert, A comparison of differential quadrature methods for the solution of partial differential equations, Comp. Meth. Appl. Mech. Engrg. 104 (1993) 295-316. [294] N . Otani, A primary mechanism for spiral wave meandering, Chaos 12 (2002) 829-842. [295] A . V . Panfilov, J . P. Keener, Effects of high frequency stimulation on cardiac tissue with an inexcitable obstacle, J . theor. Biol, 163 (1993) 439-448. [296] A . M . Pertsov, J . M . Davidenko, R . Salomonsz, W . T . Baxter, J . Jalife, Spiral waves of excitation underlie reentrant activity in isolated cardiac muscle, Circ. Res. 72 (1993) 631-650. [297] B . J . Roth, A r t Winfree and the bidomain model of cardiac tissue, J . Theo. B i o l , 230 (2004) 445-449. [298] S. Takagi, et al, A physical approach to remove anatomical reentries: a bidomain study, J . Theo. Biol., 230 (2004) 489-497. [299] K . H . W . J . ten Tusscher, A . V . Panfilov, Influence of nonexitable cells on spiral breakup in two-dimensional and three-dimensional excitable media, Phys. Rev. E , 68 (2003) 062902(1-4) 159 Bibliography [300] A . T . Winfree, Varieties of spiral wave behavior: A n experimentalist approach to the theory of excitable media, Chaos, 1 (1991) 303-334. [301] E . J . Vigmond et. al., Computational Techniques for solving the bidomain equations i n three dimensions, I E E E Trans. Biomed. E n g . 49 (2002) 1260-1269. [302] H . Yang, B . Shizgal. Chebyshev pseudospectral multidomain technique for viscous flow calculation, Comput. Methods A p p l . Mech. Eng. 118 (1994) 47-61. [303] H . Yang, B . Seymour, B . Shizgal, A Chebyshev pseudospectral multidomain method for steady flow past a cylinder, up to Re=150, Comput. fluids 23 (1994) 829-851. [304] Ye. A . Yermakova, A . M . Pertsov, Interaction of rotating spiral waves with a boundary, Biophys., 31 (1986) 932-940. [305] Z. Zhang, K . T. Ng, Two dimensional Chebyshev pseudospectral modelling of cardiac propagation, Med. and Biol. Eng. and Comp. 38 (2000) 311-318. [306] S. Zhao and G . W . Wei, Comparison of the discrete singular convolution and three other numerical schemes for solving Fisher's equation, S I A M J . Sci. Comput. 25 (2003) 127-147. 160 Appendix A Beeler-Reuter equations, The I currents in E q . (4.2) satisfy, ion I kl = 1-4 exp(0.04(V + 85)) - 1 exp(0.08(V + 53)) + exp(0.04(V + 53)) I (A.2) exp(0.04(V + 35)) INU = (9N m hj + g ){V 3 a Ica = g fd(V 1 - e x p ( - 0 . 0 4 ( V + 23)) (A.l) exp(0.04(V + 77)) - 1 = 0.8xi X1 V + 23 +0.07 - E ) NaC (A.3) Na + 82.3 + 13.0287 ]n[Ca ]) (A.4) ++ s where g^ = 4 m S c m , gNaC = 0.003 m S c m , E^ = 50 m V and g = 0.09 m S c m . In this case, the ionic calcium concentration [ C a ] , in the cytosol (Eq. A.4) satisfies - 2 - 2 a a s - 2 + + d[Ca++] = - I O " I a + 0.07(10- 7 [Ca ]) 7 ++ C dt (A.5) The variables x\,m,n,j,d and / , are called gating variables and they are voltage and time dependent. They take values between zero and one. Each gating variable satisfies an O D E of the form dX Xoo — X (A.6) ~~dt and TX are given by where X = x\, m, n,j, d and / . X^ -^oo — a ax + Px x where T X C e°x( + x) v ax,(3x and the constants C , l x x C ^ = + C {V X 1 az + Px (A.7) X + C) x (A.: i = 1,7, can be found in [251]. 161 Appendix B Initial conditions used to generate spiral waves. To generate an initial condition, the variables in the B R model took initially the values given by the steady state in the standard B R equations with the original parameters [251], i.e. x\ = 0.0056, m = 0.011, h = 0.99, j = 0.97, d = 0.003, / = 1 and [Ca] = 1 x 10~ . A propagating front i n the positive x direction that evolves into a spiral wave is generated by applying for the first millisecond a stimulus current of the form 7 7app(x, y) — i l (l+exp{2.5(\x+ay\-c W (l+exp(2.5(|x+aj/|-c ))) 1 x <0 i! 2 0 x >0 (B.l) with cy = 30 and c = 25. W i t h a = 0, there is no y dependence in E q . ( B . l ) and there are two roots to I = 0.5, at x\ = —22 and x = —30.5 as shown by the horizontal parallel lines i n Fig. B . 1 A . After t = 100 ms and for a time At = 20 ms, the conductances gfja and gs are taken to be zero over the region ax > y — yo, where a = 0 and yo = —20 mm, giving the region y < —20 which is delimited by the line y = —20. The pulse propagates in the positive x direction, but only for y < —20 (shaded region). After t — 100 + Atg = 120 ms, g^a and gs are taken as their original values, resulting i n the propagation of a front with a free end. The propagating front with the free end evolves into a spiral wave. The trajectory of the tip of the spiral is shown in Fig. B . 1 A . 2 app 2 g g g Spiral tip trajectories with incident angle 6^. For the Roo case, I generate a spiral solution with a tip trajectory that has an incident angle 9i, as discussed in the previous paragraph with g^ = 2.37 m S c m , gs = 0.003 m S c m and N — 542 points i n each dimension. W i t h these parameters, the tip hits the boundary with a trajectory having an incident angle of &i = 63° as shown i n Fig. B . 1 A . a - 2 - 2 162 Appendix B. Initial conditions used to generate spiral waves. In order to construct spirals such that their trajectories hit the boundary with a particular incident angle Oi / 63°, it is necessary to rotate I and y > yo in the (x,y) plane, and modify c\,c ,yo and t , depending on 9i. app 2 g Rotation of I for a particular is obtained by taking a — tan (^I^QO"^ in E q . ( B . l ) . Thus for di = 120°, a = 1.53 and I = 0.5 gives the linear relationship between x and y shown by the parallel lines at the left lower corner of F i g . B . 1 B . The values of c\, c , yo and t are chosen by trial and error, in order to get a trajectory to hit the lower boundary. Once the values of c\, c , yo and t are found for each angle, I move only one of the two parameters, t or yo, in order to generate all the trajectories for that specific angle. For 0j = 120°, c\ — 75, c = 72, j/o = 0 mm and t — 150 ms. The conductances are taken as zero at the region x > A front with a free end (shaded region), evolves into a spiral, that hits the lower boundary with a tip trajectory having Oi = 120°. app app 2 2 g g g 2 g 163 Appendix B. Initial conditions used to generate spiral waves. 164
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Pseudospectral solutions of reaction-diffusion equations...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Pseudospectral solutions of reaction-diffusion equations that model excitable media : convergence of… Liceaga, Daniel Olmos 2007
pdf
Page Metadata
Item Metadata
Title | Pseudospectral solutions of reaction-diffusion equations that model excitable media : convergence of solutions and applications |
Creator |
Liceaga, Daniel Olmos |
Publisher | University of British Columbia |
Date Issued | 2007 |
Description | In this thesis, I develop accurate and efficient pseudospectral methods to solve Fisher's, the Fitzhugh-Nagumo and the Beeler-Reuter equations. Based on these methods, I present a study of spiral waves and their interaction with a boundary. The solutions of Fisher's equation are characterized by propagating fronts with a, shock-like wave behavior when large values of the reaction rate coefficient is taken. The pseudospectral method employed for its solution is based on Chebyshev-Gauss-Lobatto quadrature points. I compare results for a single domain as well as for a subdivision of the main domain into subintervals. Instabilities that occur in the numerical solution for a single domain, analogous to those found by others, are attributed to round-off errors arising from numerical features of the discrete second derivative matrix operator. However, accurate stable solutions of Fisher's equation are obtained with a multidomain pseudospectral method. A detailed comparison of the present approach with the use of the sinc interpolation is also carried out. Also, I present a study of the convergence of different numerical schemes in the solution of the Fitzhugh-Nagumo equations. These equations, have spatial and temporal dynamics in two different scales and the solutions exhibit shock-like waves. The numerical schemes employed are Chebyshev multidomain, Fourier pseudospectral, finite difference methods and in particular a method developed by Barkley. I consider two different models of the local dynamics. I present results for plane wave propagation in one dimension and spiral waves for two dimensions. I use an operator splitting method with the Chebyshev multidomain approach in order to reduce the computational time. I conclude this thesis by presenting a study of the interaction of a meandering spiral wave with a boundary, where the Beeler-Reuter model is considered. The phenomenon of annihilation or reflection of a spiral at the boundaries of the domain is studied, when the trajectory of the tip of a spiral wave is essentially linear. This phenomenon is analyzed in terms of the variable j, which controls the reactivation of the sodium channel in the Beeler-Reuter model. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-02-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080413 |
URI | http://hdl.handle.net/2429/31453 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2007-318997.pdf [ 7.83MB ]
- Metadata
- JSON: 831-1.0080413.json
- JSON-LD: 831-1.0080413-ld.json
- RDF/XML (Pretty): 831-1.0080413-rdf.xml
- RDF/JSON: 831-1.0080413-rdf.json
- Turtle: 831-1.0080413-turtle.txt
- N-Triples: 831-1.0080413-rdf-ntriples.txt
- Original Record: 831-1.0080413-source.json
- Full Text
- 831-1.0080413-fulltext.txt
- Citation
- 831-1.0080413.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080413/manifest