A study of wave propagation in the FitzHugh Nagumo system by Kelly Marie Paton B.Eng., The University of Saskatchewan, 2008 B.Sc., The University of Saskatchewan, 2008 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Mathematics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2011 c© Kelly Marie Paton 2011 Abstract An excitable medium has two key properties: a sufficiently large stimulus provokes an even bigger response (excitability), and immediately following a stimulus the medium cannot be excited (refractoriness). A large class of biological systems from cardiac tissue to slime mold are examples of excitable media. FitzHugh Nagumo (FHN) is the canonical model of excitable media. Its two variables are the state of excitation and refractoriness of the one- or two-dimensional medium. Although one of the simplest models, FHN ex- hibits complex dynamics that have not been fully explored. For example, it supports a stable traveling pulse solution. However, this pulse can be desta- bilized by large perturbations. In Chapter 2 I explore a one-dimensional example where the perturbation is an increasing refractory profile. This perturbation can lead to collapse of the pulse depending on the steepness of the profile, as conjectured by Keener [Keener, J. (2004) J Theo. Bio. 230(4):459-73]. In Chapter 3 I consider a perturbation in two dimensions which can cause the stable traveling pulse to wrap around the perturbation and generate self-sustaining spiral activity. The one-dimensional example for exponential refractory profiles is ex- plored numerically for a piecewise linear FHN system. Steep profiles lead to collapse while milder profiles allow propagation. The exponential profiles are used as bounds for more general profiles to predict where collapse and propagation will occur. I also make use of a singular FHN system in the limit ǫ → 0 to provide insight into the behaviours of the full FHN system for small ǫ and small diffusion. I conclude this chapter by showing analyti- cally that, in contrast to the full system, a wave in the singular system will propagate for any exponential refractory profile. ii Abstract The two-dimensional case is explored numerically in a FHN system. The use of a temporarily refractory region as a perturbation is a novel mecha- nism for generating spiral activity. Moreover, it is shown to be robust for refractory regions of a large area. This situation models the appearance of abnormal electrical activity in the heart. In particular, it models the appearance of abnormal electricity activity in undamaged cardiac tissue. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Abbreviations and Symbols . . . . . . . . . . . . . . . . . viii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . x Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Modeling excitable media . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Background . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.2 FitzHugh Nagumo ODE dynamics . . . . . . . . . . . 2 1.1.3 FitzHugh Nagumo PDE system . . . . . . . . . . . . 3 1.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2.1 Chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2.2 Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . 8 2 Front propagation in one dimension . . . . . . . . . . . . . . 9 2.1 Introduction: A study of wave propagation in FHN . . . . . 9 2.2 Singular perturbation analysis of FHN . . . . . . . . . . . . 10 2.2.1 Traveling wave solutions . . . . . . . . . . . . . . . . 10 2.2.2 The singular limit: SFHN . . . . . . . . . . . . . . . . 18 2.3 Reaction term selection . . . . . . . . . . . . . . . . . . . . . 21 iv Table of Contents 2.3.1 Piecewise linear FHN system . . . . . . . . . . . . . . 21 2.3.2 Piecewise linear SFHN system . . . . . . . . . . . . . 23 2.4 Convergence to rest or propagation in FHN . . . . . . . . . . 23 2.4.1 Simulations for an initial exponential refractory profile 24 2.4.2 Predictions for a general initial refractory profile . . . 26 2.5 Front propagation analysis for an initial exponential refrac- tory profile in SFHN . . . . . . . . . . . . . . . . . . . . . . . 37 2.6 Conclusions and future work . . . . . . . . . . . . . . . . . . 41 3 Spiral waves in two dimensions . . . . . . . . . . . . . . . . . 46 3.1 Introduction: Modeling cardiac arrhythmias . . . . . . . . . 46 3.1.1 Biological background and medical motivation . . . . 46 3.1.2 Spiral waves . . . . . . . . . . . . . . . . . . . . . . . 48 3.1.3 Mechanism for spontaneous functional reentry . . . . 49 3.2 Two-dimensional simulations . . . . . . . . . . . . . . . . . . 51 3.2.1 Simulation results . . . . . . . . . . . . . . . . . . . . 53 3.2.2 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 61 3.3 Conclusions and future work . . . . . . . . . . . . . . . . . . 66 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 Appendices A Numerical details . . . . . . . . . . . . . . . . . . . . . . . . . . 73 B Determination of collapse or propagation in Section 2.4.1 74 C Singular derivation of dmin in Section 2.4.2 . . . . . . . . . . 76 D Excitable media terminology and properties . . . . . . . . . 79 D.0.1 Terminology . . . . . . . . . . . . . . . . . . . . . . . 79 D.0.2 Properties of stimuli . . . . . . . . . . . . . . . . . . . 80 v List of Figures 1.1 The dynamics of an action potential in the FHN ODEs, (a) over time, and (b) in the phase plane. . . . . . . . . . . . . . 4 1.2 A snapshot of a traveling pulse in the FHN PDE system, (a) over space, and (b) in the spatially-parameterized (u, v) phase plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1 A qualitative (u, v) phase plane for the FHN reaction terms. . 12 2.2 Sample initial conditions in the refractory and excited vari- ables over space for simulations in Section 2.4.1 and analysis inSection 2.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.3 Propagation and collapse results for exponential refractory profiles simulated in the full FHN system, shown in the A−λ plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.4 Sample propagation and collapse behaviour for wavefronts facing exponential refractory profiles in the full FHN system. 28 2.5 Examples of square waves that fit to an exponential upper bound; an illustration of how we determine the minimum propagation distance d̃min for a propagation conjecture for general refractory profiles. . . . . . . . . . . . . . . . . . . . . 32 2.6 A comparison of sample singular and actual minimum prop- agation intervals: singular dmin values vs actual simulated d̃min. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.7 Singularly-derived dmin scores indicating minimum propaga- tion intervals for general refractory profiles with these expo- nential upper bounds. . . . . . . . . . . . . . . . . . . . . . . 35 vi List of Figures 2.8 Comparison of full and linearized leading-order front speed as a function of the refractory variable. . . . . . . . . . . . . . . 38 2.9 Sample propagation behaviour for wavefronts facing exponen- tial refractory profiles in the singular and full systems. . . . . 42 2.10 A sample exponential refractory profile that leads to collapse in the full system but propagation in the singular system. . . 43 3.1 Diagram of the electrical conduction system in the heart. . . 47 3.2 Sample initial conditions in (a) the excitable variable u, and (b) the refractory variable v, for 2D simulations. . . . . . . . 52 3.3 Reentrant activity dependent on refractory rectangle dimen- sions in 2D simulations. . . . . . . . . . . . . . . . . . . . . . 54 3.4 Illustration of the common dynamics in the time interval t=0 to t=t’ for all 2D simulations. . . . . . . . . . . . . . . . . . . 56 3.5 Illustration of the dynamics at time t = t∗ indicating (a) reentry for large refractory rectangles, and (b) no reentry for small refractory rectangles (b). . . . . . . . . . . . . . . . . . 57 3.6 Numerical simulation of a two-dimensional domain showing the development of reentrant activity due to an initially re- fractory region. . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.7 Numerical simulation of reentry in a two-dimensional domain due to an initially refractory rectangular region of small width and intermediate length. . . . . . . . . . . . . . . . . . . . . . 62 vii List of Abbreviations and Symbols 1D one-dimensional, in space 2D two-dimensional, in space ODE ordinary differential equation PDE partial differential equation FHN FitzHugh Nagumo SFHN Singular FitzHugh Nagumo SA sinoatrial, as in sinoatrial node O( ) order of magnitude t fast time scale τ slow time scale x outer spatial scale ξ moving inner spatial scale u,U excitable variable v, V refractory or recovery variable f excitable reaction term g refractory reaction term D constant coefficient of diffusion ǫ small parameter a, b constant parameters c traveling front speed c0 leading order traveling front speed vz a critical refractory magnitude called the zero-speed level A magnitude of an exponential refractory profile viii List of Figures λ gradient of an exponential refractory profile λ∗ critical λ value in SFHN for wave propagation with respect to vz dmin singularly-derived minimum propagation interval d̃min actual minimum propagation interval from FHN simulations γp numerically-determined propagation curve in the A− λ plane γp̄ numerically-determined collapse curve in the A− λ plane γ ideal propagation/collapse curve between γp and γp̄ h height of a step or square wave xstep position of a step or square wave τstep time when a front reaches xstep xf time-dependent position of a wave front xf∞ time-dependent asymptotic expression for wavefront position for large τ vf refractory variable value at a front xz singularly-derived time-dependent position of vz d(λ,A, τ) distance between xf and xz d∞(λ,A) asymptotic distance between xf and xz for large τ cz wavefront speed linearized around vz c∞ asymptotic wavefront speed for large τ β an intermediate variable used in a singular calculation yf an intermediate variable used in a singular calculation t′ time when the wavefront reaches the refractory rectangle v′ refractory magnitude of refractory rectangle at time t′ v∗ refractory magnitude denoting highest relatively refractory value t∗ time at which a refractory rectangle decays to the magnitude v∗ ix Acknowledgements My thanks to . . . Eric Cytrynbaum, for his patience, enthusiasm, knowledge, and always-open door. Additional thanks for the extra things: sharing in outdoor adventures, allowing me the freedom to travel, and being a mentor and a friend (whether consciously or not). Josh Zukewich, for providing encouragement, energy, clarity and support — both mathematical and personal, long-term and last-minute — tough love, and everything else in between. Michael Ward and Eric, for their comments and guidance on this thesis and all the work leading up to it. Jessica Conway, for stepping up to make the first round of edits on the rough drafts, and for providing support and mentorship. The UBC math biology group, for the welcoming and supportive environ- ment. Lee Yupitun, for answering countless questions about deadlines and details. My parents and sister for their constant support and unconditional love. This work was supported financially by NSERC, the PIMS IGTC program in Mathematical Biology, and UBC. x For Mom and Dad. xi Chapter 1 Introduction 1.1 Modeling excitable media 1.1.1 Background The first successful mathematical description of electrophysiology was the foundational model developed by Hodgkin and Huxley in the early 1950’s [17]. Hodgkin and Huxley performed voltage clamp experiments on a squid giant axon and turned their observations of transmembrane potential, cur- rents, and conductance, into a circuit-like model. The result was a system of four ordinary differential equations (ODEs) that accurately described sig- nal propagation along an axon. Although the complex Hodgkin-Huxley equations have proven to be a quintessential and revealing model of signal propagation along a nerve, they are difficult to analyze. FitzHugh [14] and Nagumo [31] addressed this issue a decade later when they reduced the original system of four variables down to a simpler model of only two variables. Their simple model is much more amenable to analysis and it still captures the key phenomena of the dynamics: (1) a sufficiently large stimulus will trigger a large response, and (2) after such a stimulus and response, the medium requires a period of recovery time before it can be stimulated again. These two properties are described as excitable and refractory. Excitation occurs rapidly while recovery is a slow process. A medium that exhibits excitability and refractoriness is classified as an ex- citable medium. The FHN model has been adapted to model a wide range of excitable media from the laboratory-observed Belousov-Zhabotinsky chem- ical reaction [34] to slime mold amoeba [43] to cardiac tissue [1] [15] [27] [32]. 1 1.1. Modeling excitable media The following sections will introduce the FitzHugh Nagumo equations, system properties, and basic dynamics. We will conclude this chapter with an outline of the following two chapters of this thesis, each of which explores a specific dynamic in a certain instance of the FitzHugh Nagumo (FHN) system. 1.1.2 FitzHugh Nagumo ODE dynamics The general form of the non dimensional FitzHugh Nagumo system consists of two ordinary differential equations (ODES): du dt = f(u, v) + I(t), (1.1) dv dt = ǫg(u, v), The variable u represents the state of excitation of the medium and v rep- resents the refractoriness of the medium. The small positive parameter ǫ separates the time scales of the slow refractory variable and the fast ex- citable one. The dynamics of u and v are described by the reaction terms, f and g. The time dependent current term, I(t), models the standard method of applying a stimulus. The excitable reaction term f is a nonlinear cubic-like function in u that is monotone decreasing in v. For fixed v within some bistable interval (vmin, vmax), it has three zeros: the low (or recovered) state uL(v), the threshold um(v), and the high (or excited) state uH(v). The outer two zeros are stable and the middle zero is unstable. For high refractory values outside the bistable interval, f only has the low zero uL(v). This represents the inability of the system to undergo a large excursion when it is in a refractory state. The common model for the recovery variable which I use throughout this thesis is linear dynamics of the form g(u, v) = u− bv (1.2) The constant b > 0 is chosen small enough so that the level curves f(u, v) = 0 2 1.1. Modeling excitable media and g(u, v, ) = 0 (called nullclines) intersect only once in the phase plane, resulting in only one stable fixed point that we refer to as the rest state. The kinetics in g ensure that an excited system tends to become refractory, and a refractory system, once it returns to an unexcited state, will eventually recover toward rest. When the recovery effects of g are combined with the bistable thresh- old effects of f , the FHN ODE system effectively models the characteristic excited-recovery response of an excitable system to an above-threshold stim- ulus. This response is called an action potential. See the description of an action potential and its projection in the phase plane in Figure 1.1. 1.1.3 FitzHugh Nagumo PDE system Spatially-distributed excitable media exhibit the property of active wave propagation, which is accurately modelled by reaction-diffusion equations (see the review paper [37] for an excellent overview of wave propagation in excitable media). A classic example of an actively propagating wave in ex- citable media is a forest fire spreading through a dry forest. A sufficiently large stimulus starts a fire. One burning tree ‘excites’ the next, and the fire propagates through space. Once the fire has passed, a new flame front can- not propagate through the barren ground until the forest has had sufficient time to grow back; this is the refractory period. Thus, spatially-distributed media not only experience action potentials, they also propagate them. A more general term for this stable propagating action potential solution is a traveling pulse. The spatially distributed version of the FHN system assumes that the excitation variable diffuses. This gives a reaction-diffusion equation of the form ∂u ∂t = ∇ · (D · ∇u) + f(u, v), (1.3) ∂v ∂t = ǫg(u, v). This set of coupled partial differential equations (PDEs) is a simple model 3 1.1. Modeling excitable media 0 50 100 150 200 250 −0.5 0 0.5 1 t, time (a) u , v v, refractory u, excited −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −0.1 0 0.1 0.2 u, excitable (b) v, r e fra ct or y Figure 1.1: The dynamics of an action potential in the FHN ODEs (1.1). (a) The changing profiles of u, the state of excitation (blue solid line), and v, refractoriness (green dashed line), over time during an action potential. (b) The (u, v) phase plane depicting the trajectory (solid red line) of the action potential over time. A red dot is plotted every second time step. The system begins at the origin (rest). An applied current of I(t) = 0.25 · H(t − 50) · H(52 − t) stimulates the system through a coun- terclockwise rotation around the loop and back to rest. The spacing of the dots illustrates the fast changes in excitability and slow changes in refrac- toriness. Note the particularly close spacing of the dots on the left-hand side of the loop, indicating the slow return to rest from the recovered state. The cubic function is the level curve of f(u, v) = 0 and the linear function is the level curve of g(u, v) = 0 (dashed black lines). Reaction terms are f(u, v) = u(1− u)(u− a)− v, and g(u, v) = u− bv, with constants a = 0.1, b = 0.5, ǫ = 0.01. Numerically integrated using an explicit Euler method with a time step of dt = 0.4. 4 1.1. Modeling excitable media of generic excitable media that supports active wave propagation. In this thesis, we work with the nondimensional system and assume that D is a small nondimensional constant parameter, which is the case on sufficiently large domains [23]. From here on, we refer to this PDE system simply as the general FHN system. Traveling wave solutions When we consider small ǫ, the traveling wave profile is characterized by an abrupt rise to the excited state and a drop back down to the refractory state. We call these transitions a front and a back, respectively (see Figure 1.2). The ‘spatially-parameterized’ (u, v) phase plane is a useful tool for vi- sualizing a traveling wave profile in the PDE system (Figure 1.2). At any one step in time, the standard phase plane is a representation of the current state of the ODE system. When considering a spatially-distributed system, we can project the state of each point in space onto the phase plane, re- sulting in a spatially-parameterized phase plane. Projecting the points of a full pulse solution at one step in time onto the spatially-parameterized phase plane results in a loop. The top of the loop is the pulse back; the bottom of the loop is the pulse front. As we let time run, we can follow the trajectory of each point. A right-moving traveling pulse is comprised of points that move around the loop clockwise as time runs. See Figure 1.2 for a snapshot of a traveling pulse profile over space and its projection onto the spatially-parameterized phase plane at one step in time. The traveling pulse can exhibit rich behaviours when faced with an ob- stacle. An obstacle can be an inhomogeneous region that is non-excitable or less excitable than the rest of the domain. This results in a heterogeneous domain. A homogeneous domain can experience an obstacle in the form of a perturbation such as a refractory bump or an added transient stimulus. In one spatial dimension, inhomogeneities can lead to propagation failure or more exotic wave-reflection behaviours [28]. In two spatial dimensions, in- homogeneities can lead to rich spatiotemporal patterns. I consider obstacles in the form of refractory perturbations in a homogeneous domain. Chapter 5 1.1. Modeling excitable media −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 0.05 0.1 u, excitable (a) v, r e fra ct or y front back 0 0.5 1 1.5 2 2.5 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x, space (b) u ,v frontback v, refractory u, excited Figure 1.2: A traveling pulse at one point in time during a simula- tion of a one-dimensional FHN system of the form (1.3). (a) The profile of a right-moving traveling pulse in u (blue solid line) and v (green dashed line) over space. (b) The spatially-parameterized (u, v) phase plane depicting a parametric projection (solid red line) of the traveling pulse in (a). Every (u, v) point in space maps to a point in the phase plane. The cu- bic function is the nullcline of f(u, v) and the linear function is the nullcline of g(u, v) (dashed black lines). As the pulse moves from left to right, the points rotate clockwise around the loop in the spatially-parameterized phase plane. Reaction terms are f(u, v) = u(1−u)(u−a)−v, and g(u, v) = u−bv, with constants a = 0.139, b = 2.54, D = 0.001, ǫ = 0.008. Parameter values within the appropriate range for human cardiac tissue were taken from [15]. Numerically integrated using a semi-implicit Crank-Nicholson method [7] with grid spacing dx = 0.0224 and time step dt = 0.1, on a periodic domain. 6 1.2. Outline 2 considers propagation behaviours in the one-dimensional case. Chapter 3 considers the two-dimensional case. 1.2 Outline 1.2.1 Chapter 2 In one spatial dimension (1D), the steady-state solutions of the FHN system on a homogeneous domain are a traveling pulse or rest. If we begin with an initial condition that includes a traveling pulse and a refractory obstacle, the pulse can either experience (a) continued propagation, or (b) collapse to rest [22]. When will an initial condition converge to one or the other? Keener states without proof in [22] that the convergence of a traveling pulse in the FHN system to either propagation or collapse depends on the steepness of the recovery profile near the wave. Since the propagation behaviour of a pulse is largely determined by whether the front can propagate or not, I explore Keener’s conjecture in Chapter 2 by considering the simplified condition of a wavefront approaching an increasing refractory profile in the FHN system. I begin Chapter 2 with an asymptotic discussion of the traveling pulse in the FHN system for small ǫ, concluding by applying the limit of ǫ → 0 to derive a simplified singular FHN system that is conducive to analytical expressions. These are referred to as the FHN or full system (for small ǫ) and the SFHN or singular system (for ǫ → 0). Although my goal is to provide insight into the convergence of a wavefront in the full system, analytical estimations provided by the singular system are invaluable. As a first step, I classify a family of initial conditions (increasing refrac- tory profiles of exponential shape) as converging to either front propagation or front collapse in numerical simulations of the FHN system. This work supports Keener’s conjecture: steeper refractory profiles lead to collapse. The propagation results for the family of exponential refractory profiles are combined with analytical estimations from the SFHN system to provide con- ditions for which general classes of initial conditions with a propagating wave 7 1.2. Outline front can be expected to converge to front propagation or front collapse. During the derivation of the singular system, we point out that it is not always an accurate approximation of the full system. There are cases where the two systems display qualitatively different front propagation behaviours. I provide an analytical demonstration of this disagreement for the exponen- tial recovery profile that was explored numerically in the FHN system. We conclude with suggestions for future work. 1.2.2 Chapter 3 In two spatial dimensions (2D), complicated propagating wave formations such as spiral waves — a curved traveling wavefront of excitation spiralling around a non-excited non-refractory core — have been studied numerically and analytically in the FHN system. Spiral waves are the accepted paradigm used to model abnormal electricity activity in the heart [5]. This application is discussed at the start of Chapter 3. The rest of the chapter focuses on the appearance of spiral activity in the presence of a refractory obstacle in a 2D homogeneous medium. This situation is an accepted model for the appearance of certain types of abnormal heartbeats. The analysis in Chapter 3 is numerical and qualitative. Simulations are performed with a particular instance of the FHN equations that is adapted specifically to model human cardiac tissue [15]. I demonstrate numerically that it is possible to generate spiral waves by imposing a temporarily refrac- tory region in front of a propagating wavefront. This is a novel mechanism for spiral formation. Moreover, it is a reasonable suggestion for a cause of abnormal activity in the heart. We discuss the robustness of the mechanism to our modeling choices and perform a basic order-of-magnitude check of bi- ological plausibility. Chapter 3 is concluded with two suggestions of future work that could provide further insight into this mechanism. 8 Chapter 2 Front propagation in one dimension 2.1 Introduction: A study of wave propagation in FHN In this chapter we consider wave propagation in a one-dimensional FitzHugh Nagumo system of the general form (1.3). In particular, we study the steady state behaviour of a large perturbation of a pulse solution (see the illustration in Figure 1.2) on a homogeneous domain. Under certain conditions large perturbations of the traveling pulse solution can cause the pulse to become unstable and collapse [22]. The system converges to a spatially homogeneous rest steady state in the event of pulse collapse. Thus, given a propagating wave pulse, we expect that the final state of the system will either be (a) continued propagation or (b) collapse to the rest state. Keener conjectures that collapse is prevented as long as the spatial gradi- ent of the refractory profile is not too steep [22]. Thus, the large perturbation we consider is an increasing refractory profile in front of the wave. In the limit of small ǫ, the front of the pulse is sharp. To simplify our analysis, we consider only a traveling wave front rather than the full pulse solution (which could be considered a wave front connected to a wave back). If a wave front stalls, then the corresponding front in a pulse would also stall, allowing the back to catch up. The pulse would collapse to rest. All the FHN propagation analysis is carried out numerically as it is difficult or impossible to get exact analytical results from the FHN system of coupled PDEs. We also consider the FHN system in the singular limit, 9 2.2. Singular perturbation analysis of FHN termed the singular FHN (SFHN) system. The considerably reduced SFHN system of ODEs is more amenable to analysis. Thus, as long as the singular system appropriately approximates the full system, it is a very useful tool in understanding the complex full system. We also discuss the limitations of this singular approximation. 2.2 Singular perturbation analysis of FHN In this chapter we take the diffusion constant D = ǫ2, 0 < ǫ << 1, to consider small diffusion on the length scale of the 1D domain. Time is rescaled from the transient time scale t to the slower time scale τ = ǫt; we show below that this is the time scale of wave propagation, and thus the time scale we are interested in examining. Variables are u = u(x, τ) and v = v(x, τ). The general FHN PDE system (1.3) becomes ǫuτ = ǫ 2 d 2u dx2 + f(u, v), (2.1) vτ = g(u, v). We derive leading-order behaviour and singular properties of wavefront solutions of (2.1) in Section 2.2.1 below. Taking the limit ǫ → 0 in Sec- tion 2.2.2 yields a singular system that is simple, analytically tractable, and desirable as an alternative model to the full system (2.1). We keep the re- action terms nonspecific throughout this section so that our results are as general as possible. We choose particular reaction terms in Section 2.3 to allow numerical simulations in Section 2.4 and explicit analytical expressions in Section 2.5. 2.2.1 Traveling wave solutions The existence and properties of a traveling wave solution to system (2.1) are well established (see [21], [23], [18], [13]). Since Section 2.4 and Section 2.5 rely on the singular properties of system (1.3), key parts of the singular perturbation analysis for small ǫ are reproduced here in Section 2.2.1. All 10 2.2. Singular perturbation analysis of FHN of Section 2.2.1 is established work; some readers may choose to skip the singular perturbation analysis and begin with the definition of the singular system in Section 2.2.2. Most of space is described by evolution on an outer spatial scale. Due to the bistable nature of f , there are two stable outer solution branches. Transition layers between these outer branches occur on an inner spatial scale. We separately consider the outer and inner spatial scales, then per- form asymptotic matching to derive (a) an expression indicating the sign of the leading-order speed of a traveling wave solution, and (b) the refractory variable value that results in a front speed of zero. Consideration is given to both the transient time scale, t, and the slower time scale, τ = ǫt. We frequently refer to the dynamics pictured in the phase plane in Figure 2.1. Notation The following asymptotic analysis uses these conventions: • t: fast/transient time scale, O(1). Variables have hats; i.e., û. • τ : slow time scale, O(1/ǫ). Variables are unadorned. • x: outer spatial scale, O(1). Diffusion is O(ǫ2). Variables are lower- case. • ξ: inner spatial scale, O(1/ǫ). Diffusion is O(1). Variables are upper- case. Outer spatial scale On the outer spatial scale, the dynamics of (2.1) take place on two time scales: the transient t = O(1) time scale of dynamics in the excited variable, u, and the slower τ = O(1/ǫ) time scale of changes mainly in the refractory variable, v. We first consider the transient time scale, t = τ/ǫ. The dynamics on this time scale describe the initial behaviour of the system. Expand u and 11 2.2. Singular perturbation analysis of FHN !"!#$ !"!#%&' !"!#%() !"#$%&'( )"#$%&'( *#+,!- *#%,!- *#.,!- !"#$%&'( !"#$%&)( *"#$%&)( *"#$%&'( *+,- *+,- ).*/ ).*/ ! " Figure 2.1: A qualitative (u, v) phase plane for the FHN (1.3) re- action terms. The cubic-like function is the nullcline of f(u, v) and the linear function is the nullcline of g(u, v) (solid lines). Their intersection is a stable fixed point. Values vmax and vmin (dashed lines) enclose the bistable range of f . The middle branch of the f nullcline, um(v), is unstable and acts as a threshold to separate solutions that go rapidly to the outer branches. Solutions on the outer branches slowly evolve upward (on the excited far- right branch uH(v)) or downward (on the refractory far-left branch uL(v)). Solutions on the excited branch eventually fall off the knee near v = vmax and rapidly travel across to the refractory branch. Solutions on the refrac- tory branch travel slowly down to the fixed point. When higher-order effects (i.e., diffusion) are included in the system, solutions at the fixed point can be pushed across the separatrix; this is the mechanism that enables a trav- eling wave. Below the zero–speed level vz (dashed line), fronts travel with positive speed. Above, fronts move at a negative speed. Slow and fast time scales differ by an order of magnitude ǫ. 12 2.2. Singular perturbation analysis of FHN v asymptotically as û(x, t) = û0(x, t) + ǫû1(x, t) + . . . , (2.2) v̂(x, t) = v̂0(x, t) + ǫv̂1(x, t) + . . . , where we’ve denoted the ‘fast’ variables with a hat ˆ (û(x, t) = u(x, ǫt)). Set ǫ = 0 to get the leading order behaviour of (2.1) on the fast time scale: û0t = f(û0, v̂0), (2.3) v̂0t = 0. Given some spatially distributed initial conditions in u and v, the leading- order refractory variable on this transient time scale, v̂0, does not change: v̂0(x, t) = v̂0(x, 0). (2.4) The leading-order excitable variable û0 goes to one of the two steady states described by f(û0, v̂0) = 0: lim t→∞ û0(x, t) = { uH(v̂0), û0(x, 0) > um(v̂0(x, 0)), uL(v̂0), û0(x, 0) < um(v̂0(x, 0)). (2.5) The middle unstable branch of f = 0, um(v), is the separatrix that describes the basin of attraction of each initial condition within the interval v ∈(vmin, vmax). Initial conditions above (below) this interval in v travel horizontally in the phase plane to uL(v̂0) (uH(v̂0)). (See Figure 2.1.) Dependent on initial conditions, we can have several intervals in space that alternate between the excited and refractory steady states. A transition between the two branches is a wave front or back. We next consider the slow time scale τ = ǫt. Expand the slow variables in (2.1) as u(x, τ) = u0(x, τ) + ǫu1(x, τ) + . . . , (2.6) v(x, τ) = v0(x, τ) + ǫv1(x, τ) + . . . . 13 2.2. Singular perturbation analysis of FHN Again setting ǫ = 0, the leading order equations at this slow time scale become 0 = f(u0, v0), (2.7) v0τ = g(u0, v0). The steady-state solutions from the transient time scale (2.4) (2.5) are the initial conditions for this slow time scale. As v0 evolves according to g(u0, v0), u0 moves along its stable branch of f(u0, v̂0) = 0. This can be described concisely as v0τ = { g(uH(v0), v0), u0(x, 0) = uH(v̂0), g(uL(v0), v0), u0(x, 0) = uL(v̂0). (2.8) See the ‘slow’ arrows in Figure 2.1. Solutions on the rest branch, uL(v), move down to the stable fixed point. Solutions evolving upward along the excited branch, uH(v), will eventually ‘fall off’ the knee of the nullcline. When this happens, the steady-state condition from the transient period is no longer satisfied and the system experiences another fast time transition period (described by the fast dy- namics (2.3)) as it moves horizontally to the refractory branch. Once on the refractory branch, it will move down to the stable fixed point in the slow time τ . See these fast and slow trajectories in Figure 2.1. This outer spatial scale analysis considered the singular limit of ǫ = 0. Higher order diffusive effects for nonzero small ǫ can sufficiently perturb solutions from the stable fixed point and induce fast horizontal transitions to the excited branch (again, see Figure 2.1). Diffusive coupling near the knee of the excited branch can also pull excited solutions over to the refractory branch. Thus, discontinuities propagate due to diffusion. Inner spatial scale Wave fronts and backs in the outer scale appear as discontinuities. To resolve what happens in these transition layers, we will define a new inner 14 2.2. Singular perturbation analysis of FHN spatial scale centred at a discontinuity at x = xf . First, we define xf as the location where u attains the midpoint value between the excited and refractory branches: u(xf , τ) = 1 2 (uH(v) + uL(v)). (2.9) The outer spatial scale analysis indicates that discontinuities will propagate. If we assume xf is a function of both the fast t and slow τ , asymptotic analysis reveals xf is constant in the fast time scale. Thus, xf = xf (τ); fronts move in the slow time scale. Since we are interested in the movement of fronts, we will examine dynamics on this inner spatial scale on the slow time scale τ . Define a moving inner coordinate system as ξ = x− xf (τ) ǫ (2.10) to look for a front solution at x = xf of width ǫ. The inner variables U and V in these inner coordinates are u(xf + ǫξ, τ) = U(ξ, τ), v(xf + ǫξ, τ) = V (ξ, τ). (2.11) The inner equations of system (2.1) are Uξξ + cUξ + f(U, V ) =ǫUτ , (2.12) cVξ =ǫ[Vτ − g(U, V )]. (2.13) We have defined dxf dτ , the speed of the front, as c(V ). The dependence on V becomes clear in the following derivation of an expression for the speed of the front. In order to extract leading-order behaviour in the inner spatial scale, we 15 2.2. Singular perturbation analysis of FHN again perform an asymptotic expansion of the system variables. U(ξ, τ) = U0(ξ, τ) + ǫU1(ξ, τ) + . . . , (2.14) V (ξ, τ) = V0(ξ, τ) + ǫV1(ξ, τ) + . . . , c(V ) = c0(V ) + ǫc1(V ) + · · · = c0(V0) + O(ǫ). Thus, the leading order behaviour (ǫ→ 0) of (2.1) in the inner spatial scale (2.10) is given by: U0ξξ + c0U0ξ + f(U0, V0) =0, (2.15) c0V0ξ =0. (2.16) We look for solutions with respect to ξ; τ is simply a parameter. The general case requires c0 6= 0. From (2.16) we see that V0 is spatially uniform: V0 = V0(τ) in (vmin, vmax). This removes V0 as a variable in (2.15), leaving an ODE for U0 with respect to ξ. Asymptotic matching in the singular limit It has been proven that there exists a unique c0(V0) such that the (U0ξ , U0) phase plane for (2.15) has a heteroclinic orbit connecting the two stable zeros of the f nullcline at V0 (see [13], [23] for both existence and uniqueness of this path). Selection of a simple form for the reaction term f allows the heteroclinic trajectory to be found explicitly [23], while existence for more complicated f choices can be proved by a shooting argument [11] [16]. In either case, the heteroclinic orbit represents a front solution which moves at speed c0(V0). Our goal here is to show that there is a refractory level at which the leading-order front speed is zero.1 Matching conditions in space require continuity between the inner vari- ables U0, V0 and their respective outer u0, v0 solutions. Conditions that 1There are actually two heteroclinic orbits for a given V0 less than the zero-speed refractory level. The second heteroclinic orbit represents a front moving in the opposite direction, with speed −c0(V0). 16 2.2. Singular perturbation analysis of FHN match the outer refractory variable across the inner region are simple: lim ξ→±∞ V0 = lim x→x± f v0 (2.17) The constant inner refractory variable is determined by the value of the continuous outer refractory variable at the front. Assuming a transition from the excited branch (for x < xf ) to the refractory branch (for x > xf ), the matching conditions for the excited variable across the inner region are lim ξ→∞ U0 = uL(v0), (2.18) lim ξ→−∞ U0 = uH(v0), lim ξ→±∞ U0ξ = 0. Multiplying (2.15) by U0ξ , integrating over the inner layer, applying the conditions in (2.18), and simplifying, yields the following expression: c0(V0) ∫ ∞ −∞ U0ξ 2 dξ = ∫ uH(V0) uL(V0) f(U0, V0) dU0. (2.19) Since the integral on the left-hand side will always be positive, the sign of c0 will correspond to that of the integral on the right-hand side. Due to the cubic nature of f as defined in 1.1.2, there is a V0 value at which the numerator of the fraction in (2.19) will evaluate to zero. We define it as vz and refer to it as the zero-speed level:∫ uH (v) uL(v) f(U0, vz) dU0 = 0. (2.20) Why is the zero-speed level so important to the question of propagation? We have already noted that if the front of a pulse stalls, the pulse will collapse to rest. The zero-speed level indicates the critical refractory level at which a front will stall. We note that the vz defined in (2.20) only requires that the leading-order 17 2.2. Singular perturbation analysis of FHN speed c0 is zero. This may be misleading if the higher-order terms in the asymptotic expansion for c are nonzero. This situation is discussed below in Section 2.2.2. 2.2.2 The singular limit: SFHN When we apply the singular limit (ǫ→ 0), the FHN system (2.1) of coupled PDEs is reduced to a vastly simplified ODE system that describes changes on the outer spatial scale, x, over the slow time scale, τ . With ǫ = 0, only the leading-order behaviours remain. Ideally, we would like to use the singular FHN system as a substitute for the full FHN system. Advantages include analytical tractability and simpler numerical solving schemes. We introduce this singular system below. In Section 2.4, we employ the singular FHN system as a very useful tool to predict the behaviour of wavefronts in the full system for small ǫ. However, there are certain cases where the singular reduction does not approximate the full system. These limitations are discussed below and a calculation demonstrating a particular instance of disagreement is shown in Section 2.5. Singular FHN system The leading-order behaviour of the full system (2.1) on the outer spatial scale, x, and slow time scale, τ , can be described by: 1. The evolution of the refractory variable, v = v(x, τ), on the excited and refractory branches, and 2. The movement of fronts and backs (discontinuities in the excited vari- able, u, that separate the branches). The slow outer evolution along the branches is described in (2.7). The movement of fronts and backs can be expressed in terms of the speed c for ǫ = 0. The leading-order speed c0 can be found analytically for a simple f or numerically for a more complicated f [23]. If we identify multiple time- dependent front and back locations as xfi = xfi(τ), indexed by i = 1, 2, . . . , 18 2.2. Singular perturbation analysis of FHN where i even (odd) implies a transition from excited (rest) to rest (excited) for increasing x, the singular FitzHugh Nagumo (SFHN) system is expressed as dv dτ = { g(uL(v)), x ∈ (xfi , xfi+1), i even (rest branch) g(uH(v)), x ∈ (xfi , xfi+1), i odd (excited branch) (2.21) dxfi dτ = (−1)i+1c0(v(xfi , τ)). (2.22) The factor of (−1)i+1 in (2.22) ensures back and front speeds are assigned the appropriate sign. We refer to this as the singular FHN (SFHN) system (in contrast to the full FHN system (2.1)). A cautionary tale Elegant and simple though it may be, there are restrictions to the usefulness of the singular system defined in (2.21) and (2.22) as an approximation to the full system. First, it is only valid at the end of the transient time period during which fronts form from the initial conditions. We always choose initial conditions that consist of a fully-formed wavefront in order to bypass the transient period that is not accurately represented by the singular system. Second, the behaviour of the full system (0 < ǫ << 1) diverges from the singular system (ǫ = 0) near the regime of zero-speed wavefronts (as discussed in the Appendix of [8]). As we approach the zero-speed level vz, the leading-order speed c0 goes to zero and the O(ǫ) correction term c1 becomes important. As soon as c0 becomes O(ǫ), then it can be absorbed into c1, and the leading-order speed is now c1 ∼ O(ǫ). The leading-order equations (2.15) and (2.16) change to U0ξξ + f(U0, V0) =0, (2.23) c1V0ξ =V0τ − g(U0, V0), (2.24) and we no longer have two simple ODEs. The recovery dynamics become 19 2.2. Singular perturbation analysis of FHN much more complicated. In particular, (2.24) indicates V0 is no longer nec- essarily constant over the inner layer. As a result, (2.23) may no longer support a simple traveling wave solution in U0. If a standard traveling wave solution still exists, it should be possible to derive the correction term c1. If both variables experience more complicated dynamics, then c1 may no longer have a clear representation as the O(ǫ) correction term to the wave speed. In either case, the singular equations may need to be modified to consider more complicated dynamics in either or both variables, and they will almost certainly lose their uncoupled simplicity. Exploration of the structure and feasibility of a modified singular system near the zero-speed regime is left as future work. We conclude that we must be cautious with our assumption that c0 always represents the speed of a wave in the FHN system, and instead recognize that this assumption is only valid sufficiently far away from the zero-speed situation. Since the SFHN system adopts the O(1) speed c0, it is not always an appropriate approximation to the FHN system. In addition, equation (2.20) for vz is derived by requiring c0 = 0, which ignores the impact of a non-zero O(ǫ) correction term c1. Thus, the definition of the zero-speed level vz too is suspect, and we can expect the actual ‘zero-speed level’ in the full system to be perturbed from the singularly-derived vz in (2.20). Where can we use the SFHN system to approximate the FHN system, and where do we need to exercise caution? The potential discrepancy arises only within the regime of zero-speed traveling wavefronts. For O(1) wave speeds, the description of the singular recovery dynamics in (2.21) is valid outside the O(ǫ)-width of the wavefront layer, and (2.22) accurately de- scribes the translation of the wavefront layers. Section 2.4 carefully and successfully uses the SFHN properties to predict wavefront propagation for certain refractory profiles in the full system. In contrast, Section 2.5 provides an example of when the SHFN system can fail to capture the qualitative propagation behaviour of the full system. 20 2.3. Reaction term selection 2.3 Reaction term selection The following sections numerically investigate wave propagation in FHN and use the analytical properties of the related singular system to predict prop- agation behaviour. Numerical simulations will require an explicitly defined excitable reaction term, f . Below we motivate our choice of a specific f and give the resulting fully-defined instances of the FHN (Section 2.3.1) and SFHN (Section 2.3.2) systems that will be used in the rest of this chapter. Results from Section 2.4 and Section 2.5 may not extend to different choices of f . Both Section 2.4 and Section 2.5 rely on the analytical tractability of the SFHN system. Motivated by previous work such as [23], [26], and [36], we choose to capture the excitable reaction term’s cubic-like behaviour using a simple piecewise linear function f(u, v) = H(u− a)− u− v (2.25) where H(u − a) is the Heaviside step function. Due to symmetry about a = 1/2, we take a ∈ (0, 1/2). The piecewise linearity of f allows us to carry out explicit calculations that would be difficult (or impossible) with a nonlinear system. For example, the high and low stable states of u can be easily calculated as uH(v) = 1− v, (2.26) and uL(v) = −v, (2.27) respectively. 2.3.1 Piecewise linear FHN system With the reaction terms (1.2) and (2.25), the FHN system (2.1) becomes ǫuτ = ǫ 2uxx +H(u− a)− u− v, (2.28) vτ = u− bv, 21 2.3. Reaction term selection System (2.28) is the full FHN system used for the rest of this chapter. The a value must be chosen small enough for the system to be excitable. We choose the b value small enough with respect to a to ensure the system has only one stable fixed point; this restriction on b with respect to a is easily calculated. We choose to fix a = 0.1, b = 0.1, and ǫ = 0.01 for all calculations with this system. Explicit traveling wave speed For the piecewise linear f (2.25), the techniques in [23] can be expanded upon to explicitly find both the leading-order speed of the traveling wave and the related zero-speed level. The speed will be used below in the explicit definition of the singular system corresponding to the piecewise linear FHN system (2.28). Steps for deriving the speed are briefly summarized below. The details are not hard to reproduce. The front position in this instance is explicitly defined via (2.9) as u(xf , τ) = 1 2 − v(xf , τ). (2.29) Leading-order inner spatial scale front profile solutions can be found by solving the ODE for U0 in (2.15) and applying the matching conditions with the outer spatial scale (2.18). Matching the inner spatial scale wavefront solutions across the position of the front leads to an explicit expression for the leading-order front speed c0(V0), c0(V0) = √ 1− V0 − a a+ V0 − √ a+ V0 1− V0 − a, (2.30) a function only of the V0 value at the front xf . This easily-derived an- alytical expression for the wave speed is a distinct advantage of choosing the piecewise linear f over more complicated functions. As a compara- tive example, consider the common choice of f as the cubic polynomial, f(u, v) = u(1− u)(u− a)− v; finding an analytical expression for the speed in this case requires the lengthy analytical expressions for the roots of a cubic. 22 2.4. Convergence to rest or propagation in FHN The zero-speed level for this system can be found by evaluating (2.20) or setting (2.30) to zero and solving for vz: vz = 1 2 − a. (2.31) This expression is used in the FHN convergence predictions in Section 2.4 and the SFHN propagation analysis in Section 2.5. 2.3.2 Piecewise linear SFHN system The SFHN system that corresponds to the piecewise linear full system (2.28) is found by evaluating (2.21) with the linear g (1.2), the explicit excited and refractory branches (2.26) (2.27), and the explicit leading-order speed (2.30): dv dτ = { −(b+ 1)v, x ∈ (xfi , xfi+1), i even 1− (b+ 1)v, x ∈ (xfi , xfi+1), i odd (2.32) dxfi dτ = (−1)i+1 (√ 1− v − a a+ v − √ a+ v 1− v − a ) . (2.33) Again, simulations are performed with fixed a = 0.1, b = 0.1. System (2.32) is the SFHN system that is used for the remainder of this chapter. 2.4 Convergence to rest or propagation in FHN In this section, I investigate front propagation in the piecewise linear FHN system (2.28). The two qualitatively different steady-state behaviours in the presence of an obstacle are collapse or propagation. The ‘obstacle’ consid- ered here is an increasing refractory profile, motivated by Keener’s conjec- ture that collapse will not occur in a homogeneous domain as long as the spatial refractory gradient is not too large in front of the wave [22]. His claim has been supported by numerical simulations [8]. The goal in this sec- tion is to determine clear conditions under which a refractory profile leads to collapse or continued propagation. As a first step toward this goal, I begin this section by numerically sim- 23 2.4. Convergence to rest or propagation in FHN ulating a wavefront in the FHN system (2.28) with a class of 2-parameter exponential refractory profiles. An exponential profile has the variable re- fractory magnitude and gradient required to test Keener’s conjecture. More- over, it allows analytical calculations when considered in the singular system (which is taken advantage of in the next section). I record whether a given exponential refractory profile leads to rest or propagation. The results sup- port Keener’s collapse conjecture. Those results are then directly applied to predict collapse and propagation conditions for general refractory profiles in Section 2.4.2. 2.4.1 Simulations for an initial exponential refractory profile The initial condition for the excitable variable is: u(x, 0) = { uH(v(x, 0)), x < 0, uL(v(x, 0)), x > 0. (2.34) This is a singular wavefront solution with the front at x = 0. This allows for later comparison with the singular system. Without loss of generality, we always restrict the initial refractory value at the front as v(0, 0) < vz so that the front is initially propagating to the right (positive x). Also, note the dependence on v. For the sake of analytical tractability in later comparisons with the sin- gular system, we examine an exponential recovery profile of the form: v(x, 0) = Aeλx, (2.35) with nontrivial λ > 0 and A < vz. See Figure 2.2 for a visual depiction of the initial conditions (2.34) (2.35). A range of exponential profiles is simulated in the FHN system (2.28). Details of numerical simulations can be found in Appendix A. Each profile is classified as leading to collapse or propagation. For details on how we determine collapse versus propagation, see Appendix B. 24 2.4. Convergence to rest or propagation in FHN −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x uH(v(x),t=0), x<0 uL(v(x),t=0), x>0 v(x,t=0)=A eλ x v z u, excited v, refractory Figure 2.2: Sample initial conditions for simulations in Section 2.4.1 and analysis in Section 2.5. A wavefront in u centred at x = 0 (2.34) and an exponential profile in v (2.35). In the full FHN system, the abrupt wavefront smooths out into a steep curve of width O(ǫ) as the simulation runs. The refractory variable on the lower branch (to the right of the front) decays in time. If the wavefront gets sufficiently close to a point in the refractory profile that is at the zero-speed level vz, the front stalls and leads to collapse. Otherwise, it leads to propagation. [Shown for parameters: A = 0.25, λ = 4.] Propagation and collapse results in the A − λ plane are displayed in Figure 2.3. The actual numerical result is two monotonic numerical curves of interest: the points indicating the most extreme initial conditions that lead to propagation, γp(ǫ) (marking the upper boundary of the propagation region), and those marking the most lenient initial conditions that lead to collapse, γp̄(ǫ) (marking the lower boundary of the collapse region). We assume the existence of a curve, γ(ǫ), between these two that marks the 25 2.4. Convergence to rest or propagation in FHN boundary between propagation and collapse. As expected from Keener’s collapse conjecture and the theory of the zero-speed level, both high λ and high A values lead to propagation failure. Front trajectories from two sample A − λ pair FHN simulations indicating propagation and collapse are plotted in Figure 2.4. The trajectory of the singularly-derived zero-speed level (2.31) is also plotted. As the recovery variable decays in a monotonically increasing profile, vz travels to the right. In the case of propagation (Figure 2.4(a)), the front approaches a fixed speed. It also travels at a fixed distance from the zero-speed level. In the case of collapse (Figure 2.4), we notice the front only has to get slightly closer than this fixed distance to vz in order to stall. Since the refractory profile is monotone increasing, the true critical refractory value that results in a stalled wavefront is smaller than the singularly-derived vz. This true non-singular zero-speed level is indicated by the intersection of γ with the vertical A axis in Figure 2.3. We conclude that a wave in the full system has stricter refractory-variable stall conditions than a wave in the singular system. This observation reinforces our note of caution when applying the singularly-derived results for the zero-speed regime to the full system. 2.4.2 Predictions for a general initial refractory profile We’ve determined the propagation behaviour of a front facing an exponen- tial refractory profile in the piecewise linear FHN system (2.28). The next natural step is to determine whether a front will propagate or collapse given an arbitrary initial refractory profile. To our knowledge, this is an open problem. For a discussion of the limitations of classifying system conver- gence dependent on initial conditions, see [8]. In this section, we contribute to the body of knowledge surrounding this problem by developing predic- tions of propagation and collapse in the piecewise linear FHN system (2.28) for refractory profiles under certain restrictions. We consider a wavefront in u as defined in (2.34), and any refractory profile v(x, 0) that is positive and monotonically increasing on the interval [0, x1]. Call the refractory profile the test profile. Results from Section 2.4.1 26 2.4. Convergence to rest or propagation in FHN 0 1 2 3 4 5 6 7 8 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 λ A Full FHN: Propagation Full FHN: Collapse A=v z γ Figure 2.3: Propagation and collapse results for exponential re- fractory profiles simulated in the full FHN system, shown in the A − λ plane. Red dots indicate the numerically-determined FHN propa- gation/collapse curve, γ. Steep profiles and profiles of high magnitude lead to collapse. A dashed blue line indicates the singularly-derived zero-speed level, vz. The intersection of γ with the A-axis indicates the actual zero- speed level for the full system. It is slightly smaller than the singular vz. 27 2.4. Convergence to rest or propagation in FHN 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 tim e, t position, x front zero−speed (a) Propagation 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 tim e, t position, x front zero−speed (b) Collapse Figure 2.4: Sample propagation and collapse behaviour for wave- fronts facing exponential refractory profiles in the full FHN sys- tem. Plots of front (solid line) and singularly-derived zero-speed vz (dashed line) location versus time for two simulations of the full FHN system (2.28) on either side of the propagation/collapse curve γ (pictured in Figure 2.3). The upper panel in (a) displays propagation. The front approaches a con- stant speed, and travels a fixed distance away from vz. The lower panel in (b) displays collapse. Notice that the front stalls when it is close to vz but not yet at it. The true zero-speed level is smaller than the singularly-derived one. 28 2.4. Convergence to rest or propagation in FHN are used in conjunction with the leading-order behaviour of the full system to make predictions for collapse and propagation on the interval [0, x1]. The convergence predictions require the following two properties that are informed by the singular limit. First, all refractory profiles will decay, to leading order, with time constant 11+b > 0 in the region x > xf + O(ǫ). 2 Second, the leading-order front speed is a positive-valued function for v < vz that monotonically decreases as v increases to vz. 3 Call these properties 1 and 2. Convergence to collapse Here we derive a condition for the piecewise linear FHN system that de- scribes conditions on refractory profiles for which the front will collapse. We define a lower bound profile as a profile that is initially less than or equal to the test profile at every point in space within [0, x1]. We then choose an exponential profile of the form (2.35) that is a lower bound to the test profile. The condition for collapse is: Collapse conjecture: The existence of a collapse-inducing exponential lower bound for the test profile within the interval [0, x1] implies that the test pro- file will lead to collapse. The proof of the collapse conjecture is as follows: Simulate two systems simultaneously, one with the exponential lower bound and one with the test profile. The constant decay rate of the refrac- tory variable (property 1) implies that the lower bound profile will always be lesser than or equal to the test profile. Both property 1 and property 2 2This is proved by solving the singular ODE (2.32) describing the evolution of the leading-order refractory variable on the lower branch. The solution is v(x, τ ) = v(x, 0)e−(1+b)τ , x > xf +O(ǫ). (2.36) The test profile has the leading-order solution form (2.36). 3This is proved directly by examination of the expression for the leading-order front speed, (2.30). 29 2.4. Convergence to rest or propagation in FHN (the speed is a monotonically decreasing function of v) together imply that the speed of a front in the lower bound case will always be faster than the speed of a front in the test case. If the front in the lower bound case stalls at a point x′ < x1, then the front in the test case will necessarily also stall at a point x′′ ≤ x′ < x1. As long as the lower bound falls within the collapse region of the A− λ plane (see Figure 2.3), we predict that the test case will also map to collapse. The converse of this condition is not true. Collapse does not require the existence of such a lower bound. A simple counterexample is given by a step function: vstep = { 0, 0 < x < xstep h, x > xstep (2.37) where the constant h > 0 is the height of the step function that starts at x = xstep. Since any exponential function of the form (2.35) is strictly greater than zero, it is not possible to form the required lower bound. However, given a step function with h sufficiently greater than vz and xstep sufficiently close to x = 0, collapse is certain. Although this condition is derived using the properties of the SFHN system, we expect it to hold for the full FHN system. The wavefront in the singular system is a discontinuity. The u value immediately to the right of the wavefront is at the lower steady state. In contrast, the wavefront has a finite width in the full system. To the right of the front, u decays exponentially to the lower steady state over the half-width of the wavefront layer. Thus, the refractory variable in the right half-width of the wavefront layer in the full system would experience a stronger positive feedback than in the singular case due to the positive u term in g (1.2). This correction would actually lead to a stronger condition for collapse. Convergence to propagation Here we derive a condition for the piecewise linear FHN system that pre- dicts when a front will propagate for an arbitrary refractory profile. We 30 2.4. Convergence to rest or propagation in FHN define an upper bound profile as a profile that is initially greater than or equal to the test profile at every point in space within [0, x1]. We then choose an exponential profile of the form (2.35) that is an upper bound to the test profile. For each such upper bound exponential profile that allows propagation, we can calculate a minimum interval [0, d̃min) over which we can predict propagation for the test profile. The condition for propagation is: Propagation conjecture: The existence of a propagation-allowing exponential upper bound for the test profile implies that a wavefront approaching the test profile is predicted to propagate in the interval [0, d̃min), where d̃min ≤ x1 is a calculable numerical property of each exponential profile. If d̃min > x1, then the front will propagate within the interval [0, x1]. The proof of this propagation conjecture follows: Simultaneously simulate one system with the test profile and another system with the exponential upper bound profile. Property 1 (the constant decay rate of the refractory variable) indicates the upper bound will always be greater than the test profile at each point in space. The combined effects of properties 1 and 2 are more complicated than in the collapse condition, and require careful consideration. It is possible for the test profile to be much less than the upper bound for small x and then approach the upper bound at some larger x. A front in this test case would travel faster than the upper bound case for small x (due to property 2), and could potentially reach a high refractory point before it has had time to decay sufficiently (due to property 1), thus ending in collapse. We can avoid this situation by restricting the interval over which we can guarantee propagation with the upper bound principle. We call this the minimum propagation interval, [0, d̃min). It is a property of each exponential upper bound. The logic behind the calculation of d̃min is as follows. Consider the full set of all test profiles that fit under a certain exponential upper bound but lead to collapse. The ‘worst-case scenario’ in this set — the test profile that will lead to collapse within the shortest interval — is a test profile of step 31 2.4. Convergence to rest or propagation in FHN 0 0.5 1 1.5 2 2.5 3 3.5 4 0 1 2 3 4 5 6 x re fra ct or y va ria bl e, v Figure 2.5: Examples of square waves fit to an exponential upper bound to determine the minimum propagation distance d̃min. The point where a front turns around in a simulation of the full system when faced with the smallest such square wave resulting in collapse indicates the d̃min for this exponential. The step position of the smallest such square wave resulting in collapse indicates the singularly-derived approximate dmin. form (2.37). This is because a step profile permits the fastest front propa- gation to the potential stall point, xstep, (by property 2) which ensures the v value at xstep will have decayed as little as possible by the time the front reaches it (by property 1). Thus, the smallest step profile that (a) fits under a certain exponential upper bound, and (b) leads to collapse, provides a lower bound on the interval over which propagation is guaranteed for a test profile with this upper bound. We find the value of d̃min for each exponential curve by fitting step waves of successively greater magnitude under the exponential curve until we find one that results in the front stalling (see Figure 2.5). We denote the stall location of the front as d̃min. This is accomplished through simulations, and it completes the formulation of the propagation conjecture. 32 2.4. Convergence to rest or propagation in FHN As an addition to the propagation conjecture, we note that simulations in the full system are tedious and time-consuming. This makes it difficult to discern trends in d̃min for various exponential profiles. Insight can be gained into d̃min by considering the singular limit. The singular limit allows us to analytically derive an approximate expression for dmin (distinguished by the lack of a tilde). In the singular limit, the stall will occur exactly at the location of the step wave. Details of this derivation can be found in Appendix C. The result is dmin = ln(1/2−aA ) λ− (b+1) √ a(1−a) 1−2a (2.38) A comparison of the singular dmin from (2.38) and the actual d̃min from simulations for a range of exponential profiles can be seen in Figure 2.6. The two are similar enough to allow the use of the singular values to understand how the actual d̃min will depend on A and λ. See Figure 2.7(a) for singular dmin values calculated for range of A, λ for a = b = 0.1. See Figure 2.7(b) for a sample set of level curves in dmin. As expected, steeper profiles with greater magnitudes have smaller minimum propagation intervals. Notice from Figure 2.6 that the error in the singularly-derived dmin with respect to the true d̃min decreases for increasing λ. This results from the root of the discrepancies between the singular speed and the true full system’s speed. We have predicted that any divergent behaviour will occur for speeds of order ǫ. This corresponds to some small interval |v − vz| < µ(ǫ) in which v is close to vz, where µ is a small parameter that captures the discrepancy between the singularly-derived vz and the actual zero-speed level for the full system. When the slope of the profile is shallower, this 2µ-width interval is spread out over a larger region of space. Since the SFHN-FHN discrepancies occur within this interval, it is reasonable that the absolute error will increase as the space covered by this interval increases. We conclude this discussion with the comment that the propagation con- jecture will be weak for many profiles. Remember that we have focused on a calculation of the minimum interval over which we can guarantee propaga- tion in the case of any test profile. If the test profile is a square wave — the 33 2.4. Convergence to rest or propagation in FHN 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 λ d m in Simulation Analytical Figure 2.6: A comparison of sample singular and actual minimum propagation intervals. Singular dmin calculated from (2.38) (red dots) are plotted with the actual minimum intervals d̃min (blue triangles) found by numerically simulating the full system. Values are shown for a sampling of λ’s at A = 0.15, with a = b = 0.1. 34 2.4. Convergence to rest or propagation in FHN 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 λ A 0 0.5 1 1.5 2 2.5 3 3.5 (a) Singularly-derived dmin for a set of exponential profiles. Figure 2.7: Singularly-derived dmin scores, calculated for a = 0.1, b = 0.1. These are the singularly-derived minimum propagation interval size dmin via (2.38) that approximate the actual minimum interval size d̃min. Panel (a) indicates the singular dmin for each set of exponential parameters (λ,A); (b) shows a sample set of calculated level curves in dmin. 35 2.4. Convergence to rest or propagation in FHN 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 λ A d min=0.3 d min=0.5 d min=2.2 d min=4.9 (b) Sample singularly-derived dmin level curves. Figure 2.7: (cont’d) Singularly-derived dmin scores, calculated for a = 0.1, b = 0.1. These are the singularly-derived minimum propagation interval size dmin via (2.38) that approximate the actual minimum interval size d̃min. Panel (a) indicates the singular dmin for each set of exponential parameters (λ,A); (b) shows a sample set of calculated level curves in dmin. 36 2.5. Front propagation analysis for an initial exponential refractory profile in SFHN worst-case scenario — then the propagation conjecture accurately predicts the actual interval of propagation d̃min. On the other end of the spectrum, a test profile identical to its propagation-allowing exponential upper bound is guaranteed to propagate within the entire test interval, not just within the minimum guaranteed interval defined by d̃min. Most test profiles will fall between these two extremes and lead to propagation within an interval somewhat larger than [0, d̃min). 2.5 Front propagation analysis for an initial exponential refractory profile in SFHN This section mimics the study of propagation in Section 2.4.1 — a trav- eling wavefront of the form (2.34) interacting with an initial exponential refractory profile of the form (2.35) — with two major differences. First, we consider the singular system (2.32) (2.33) instead of the full system (2.28). This enables us to directly compare the results from both systems. Second, the analytical tractability of the SFHN system in the case of an exponential refractory profile allows us to approach the problem analytically instead of numerically. We take parameters a, b as unspecified constants, with the condition b < a/(1 − a) to ensure we have a valid FHN-type system with only one stable fixed point. Surprisingly, this analysis will show that the only possible outcome given these initial conditions in the singular system is continued propagation. This clearly demonstrates that although the sin- gular system is a useful tool, as has already been shown in Section 2.4, it does not accurately capture the complex behaviours of the full system that emerge in zero-speed regimes. A front will stall if it reaches a point in the recovery profile that has the value vz from (2.31). We are interested in speeds near this zero-speed location. Thus, we expand the singular piecewise linear speed given in (2.30) 37 2.5. Front propagation analysis for an initial exponential refractory profile in SFHN 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 −4 −3 −2 −1 0 1 2 refractory variable, v v=v z full speed linearized speed Figure 2.8: Comparison of full and linearized leading-order speed as a function of v. A comparison of the full speed (sold blue curve), given by (2.30), and the speed linearized about the zero-speed level (solid red curve), given by (2.39). As long as we are reasonably close to the zero-speed level (dashed black line), the linearized speed is an acceptable approximation. This is shown for fixed a = 0.1 but the agreement holds for all a ∈ (0, 1/2). around v = vz: cz(v) = −4 ( v + a− 1 2 ) + O((v − vz)3), (2.39) a valid expression for v values sufficiently close to the zero speed value vz. From odd symmetry, the quadratic-order term must be zero, leaving us with cubic-order error only. See Figure 2.8 for a visual comparison of the full singular speed c(v) and the leading-order (linearized) speed expanded near the zero speed recovery level, cz(v). This calculation assumes that the front is close enough to the zero-speed recovery level for the linearized speed to be an acceptable approximation of the front speed. 38 2.5. Front propagation analysis for an initial exponential refractory profile in SFHN Deriving an expression for the location of the front, xf (τ), with respect to the location of the zero-speed recovery level, xz(τ), leads directly to the classification of collapse versus propagation. Define the distance between the location of the front and the zero-speed level as d(λ,A, τ) = xz(τ)− xf (τ). (2.40) A positive (negative) distance indicates the zero-speed level is to the right (left) of the front. Thus, a system with an initially right-propagating front that exhibits propagation has the property d(λ,A, τ) > 0. In the following, I show that d(λ,A, τ) > 0 for any initial exponential refractory profile. Substituting the linearized speed (2.39) into the piecewise linear SFHN equation for the evolution of the front position (2.33) gives dxf dτ = −4vf + (2− 4a), (2.41) which requires an expression for the recovery variable at the front, vf (xf , τ). To the right of the wavefront, the recovery variable evolves according to (2.32): dv dτ = −v(1 + b), forx ≥ xf . (2.42) With the exponential initial condition from (2.35), the solution to this ODE is v(x, τ) = Aeλx−(1+b)τ , forx ≥ xf . (2.43) Recall the requirements A < vz and λ > 0 (Section 2.4.1). Evaluated at the front position, the recovery variable at the front is vf (xf , τ) = Ae λxf−(1+b)τ , (2.44) which can be substituted into equation (2.41): dxf dτ = −4Aeλxf−(1+b)τ + (2− 4a). (2.45) This can be solved by shifting xf = yf + (2 − 4a)τ , solving the separable 39 2.5. Front propagation analysis for an initial exponential refractory profile in SFHN ODE in yf , and then translating the solution back to xf . The resulting expression describes the time evolution of the location of the front: xf (τ) = − 1 λ ln ( 4A β (eλβτ − 1) + 1 ) + (2− 4a)τ, (2.46) where β = 2− 4a− b+1λ and we require 4Aβ (eλβτ − 1) + 1 > 0. Two distinct asymptotic results emerge for different ranges of λ. Define the critical λ as λ∗ = b+1 2−4a . For λ > λ∗, β > 0; λ < λ∗ is equivalent to β < 0. In the limit of large τ for β < 0, the transient exponential term goes to zero. Some rearranging of (2.46) for β > 0 allows similar simplifications in the limit of large τ . Thus the front location for large τ becomes xf∞(τ) ∼ (2− 4a)τ − 1 λ ln ( 1− 4Aβ ) , λ < λ∗( b+1 λ ) τ − 1λ ln ( 4A β ) , λ > λ∗ (2.47) with the restriction 1 − 4Aβ > 0 for λ < λ∗. In both cases the front speed asymptotes to some constant value for large enough τ : c∞ ∼ { 2− 4a, λ < λ∗ b+1 λ , λ > λ∗ (2.48) Next, setting v(x, τ) as given by (2.43) equal to vz from (2.31) and solving for xz gives the location of the zero-speed level: xz(t) = b+ 1 λ τ + 1 λ ln ( 1− 2a 2A ) . (2.49) Since we are interested in the ‘steady-state’ behaviour of the system, we look for the behaviour of d(λ,A, τ) in the limit of large τ : d∞(λ,A) = xz(τ)− xf∞(τ), (2.50) d∞(λ,A) = 1 λ ln [ (1− 2a) ( 1 2A − 2β )] − βτ, λ < λ∗ 1 λ ln ( 1 1−λ∗ λ ) , λ > λ∗. (2.51) 40 2.6. Conclusions and future work Evaluation of the log arguments with the restrictions A < vz, λ > 0, and b < a/(1 − a) reveals that the distances are always positive. The critical value λ = λ∗ designates a qualitative change in the relationship between the singular front and zero-speed locations. The front and zero-speed trajecto- ries are plotted in Figure 2.9 for each λ region. For steep profiles (λ > λ∗), the front approaches a fixed distance from vz, asymptotically approaching a constant speed equivalent to the constant rate of decay of the recovery profile. For mellow profiles (λ < λ∗), the front asymptotically approaches a constant speed that is slower than the constant rate of decay of the re- covery profile. Either way, an exponential refractory profile in the singular system always leads to propagation. See Figure 2.10 for a comparison of propagation in the SFHN system with collapse in the FHN system. Singular front trajectories in Figure 2.9 are plotted from simulations of the SFHN system (2.32), (2.33), but the results match the theoretical expressions (2.46) and (2.49) derived above exactly. We have also plotted the front trajectories of the full system for the corresponding simulations. Both sample simulations fall within the propagation region of the A− λ plane in Figure 2.3. Note the agreement of the speed of the front in each system past the initial transient stages. This suggests that asymptotic speeds predicted in (2.48) for the singular system can be extended to describe the full system in this regime. 2.6 Conclusions and future work Our simulations in Section 2.4.1 of a wavefront and an increasing exponential refractory profile provided support for Keener’s qualitative conjecture that a traveling wave in the FitzHugh Nagumo system will continue propagating unless the refractory gradient is too steep [22]. Section 2.4.2 then applied the propagation results from Section 2.4.1 to derive conditions on propagation or collapse for more general increasing refractory profiles. These conjectures can be succinctly summarized as follows: 1. Increasing refractory profiles that can be bounded below by an expo- 41 2.6. Conclusions and future work 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 tim e, t position, x front, FHN zero−speed, FHN front, SFHN zero−speed, SFHN (a) λ > λ∗ 0 0.5 1 1.5 2 2.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 tim e, t position, x front, FHN zero−speed, FHN front, SFHN zero−speed, SFHN (b) λ < λ∗ Figure 2.9: Sample propagation behaviour for wavefronts facing ex- ponential refractory profiles in the full and singular systems. Plots of front (solid lines) and singularly-derived zero-speed vz (dashed lines) lo- cations versus time for the singular (red) and full (blue) systems. Both plots take place in the propagation region of Figure 2.3. The upper panel in (a) displays propagation for steep profiles, λ > λ∗. The lower panel in (b) displays propagation for mellow profiles, λ < λ∗. 42 2.6. Conclusions and future work 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 tim e, t position, x front, FHN zero−speed, FHN front, SFHN zero−speed, SFHN Figure 2.10: A sample exponential refractory profile that leads to collapse in the full system but propagation in the singular system. This figure shows front (solid lines) and singularly-derived zero-speed vz (dashed lines) locations versus time for the singular (red) and full (blue) systems. The chosen exponential profile falls within the collapse region of Figure 2.3. 43 2.6. Conclusions and future work nential profile that was shown to map to collapse in Section 2.4.1 are predicted to also result in wave collapse. 2. Increasing refractory profiles that can be bounded above by an expo- nential profile that was shown to map to propagation in Section 2.4.1 are predicted to allow propagation within at least a minimum distance interval that can be explicitly determined. However, both of these predictions are weak, and could greatly benefit from further developments in future work. Although Section 2.4.2 successfully used the singular system as a guide to understand and predict the behaviour of the full system, Section 2.5 demonstrated that there are situations where the singular system fails to approximate the full system. The root of the discrepancies between the two systems was discussed in Section 2.2.2. Two directions are suggested for future work. First, a quantified degree of agreement between the two systems would help to understand when we are in a regime where propagation in the SFHN system is a good approximation to that in the FHN system. ‘Agreement’ for a set of initial conditions can be loosely defined as a high degree of similarity between the trajectories of the wavefront in each system. For example, one could simulate both systems with the same initial conditions, track the wavefront trajectories at each time point and sum up the absolute distance between them over time. The result would be zero if the trajectories agree perfectly, and would increase for decreasing agreement. The second suggestion for future work is the development of a modified SFHN system that is accurate near the zero-speed regime. This may involve solving the modified leading order equations (2.23) and (2.24) to find the O(ǫ) correction term to the leading-order wave speed, c1. This would involve a solvability condition with adjoint eigenfunctions. Ideally, the expansion of the wave speed to a two term approximation would be sufficient to in- crease the accuracy of the singular system without sacrificing its simplicity. However, we also note that the inner-layer recovery dynamics become more complicated for small speeds. If the recovery variable is no longer a constant 44 2.6. Conclusions and future work over the inner layer, the singular description may become more complicated than the current two ODEs. In either case, a valid singular approximation near the zero-speed regime would be a useful addition to the FHN toolbox. 45 Chapter 3 Spiral waves in two dimensions 3.1 Introduction: Modeling cardiac arrhythmias If left untreated, certain types of cardiac arrhythmias can be fatal. Mathe- matical modeling provides an avenue to gain insight into how arrhythmias form and how they can be treated [5]. We begin this chapter with an in- troduction to the biology of cardiac arrhythmias in Section 3.1.1, a brief discussion of the history of modeling cardiac arrhythmias as spiral waves in excitable media in Section 3.1.2, and an introduction to the mechanism of spiral formation that will be modelled in this chapter in Section 3.1.3. 3.1.1 Biological background and medical motivation A normal heartbeat is driven by pulses of electricity generated by the sinoa- trial (SA) node at regular intervals. When the SA node fires, an independent pulse travels through prescribed conduction pathways in the heart and then spreads through the walls of the ventricles [10] (see Figure 3.1). Its passage signals sequential mechanical contractions of the atria and ventricles in turn. The result is a periodic cardiac rhythm called a sinus rhythm. Defects in this electrical signalling system are the root cause of reentrant cardiac arrhyth- mias, the most common type of dangerously abnormal cardiac rhythms. A reentrant arrhythmia refers to an abnormal heartbeat caused by cir- cuitous electrical activity in the heart, termed reentry. Reentrant activity acts as an additional and abnormal source of activity in the heart that dis- turbs the normal sinus rhythm. When reentry appears in the ventricles (see 46 3.1. Introduction: Modeling cardiac arrhythmias Sinoatrial (SA) node Conduction pathways Figure 3.1: Diagram of the electrical conduction system in the heart. A normal heartbeat is initiated when the sinoatrial (SA) node fires. It discharges an electrical impulse that travels through the conductive path- ways (shown in yellow), causing the atria and ventricles to contract in turn. The impulse dissipates along the bottom of the ventricles. Reentrant cardiac arrhythmias occur when this progression is disturbed by a functional circuit in the heart. Ventricular reentry (a functional circuit in the ventricles) can be fatal. [Image adapted from [10].] 47 3.1. Introduction: Modeling cardiac arrhythmias Figure 3.1), the consequences can be fatal if it is left untreated. Sudden car- diac death due to ventricular reentrant arrhythmias causes one in six deaths in industrialized countries [41], making it one of the leading causes of death. Reentry can occur due to an anatomical or a functional circuit. An anatomical circuit describes a clear physical loop that allows contin- uously circulating activity. Anatomical circuits are created when abnormal conductive paths are formed by physical irregularities. The resultant ar- rhythmias are chronic since the anatomical circuit is an ever-present feature in the heart. In contrast to an anatomical circuit, functional reentry occurs when a wave of electricity interacts with an obstacle and is able to continue circu- lating around the resultant functional circuit. There are two distinct types of obstacles that lead to functional circuits: fixed obstacles such as regions of dead or damaged tissue, and transient obstacles such as additional waves of activity. Fixed obstacles such as regions of dead or damaged tissue are the result of previous heart disease. They can be located by medical imag- ing techniques, making possible both treatment and prevention of future arrhythmias. Transient obstacles may explain spontaneously-arising ven- tricular arrhythmias that appear in patients who have no prior history of heart disease and hence hearts with no discernible fixed obstacles. Reentrant activity arising from transient causes has proven difficult to understand and diagnose [38]. In this chapter, I explore one type of transient obstacle that leads to the appearance of functional reentry. 3.1.2 Spiral waves A spiral wave is the common mathematical paradigm used to describe reen- trant activity within excitable media, and cardiac tissue in particular [38] [42] [45].4 A literature search for spiral waves in cardiac tissue will reveal a well- established field, both in modeling theory [4][24] [25] [40] and experiments 4Spiral waves refer to a two-dimensional phenomenon (in space). The three-dimensional analog to a spiral wave is a scroll wave. In one spatial dimension, it is simply a pulse on a loop [15] [27]. 48 3.1. Introduction: Modeling cardiac arrhythmias [9] [47]. Biologically relevant questions include how spirals form, move and break up [12] [39], and are annihilated [19]; these provide insight into how arrhythmias appear, change, and can be stopped. Here we focus on the first question — the formation of spiral waves in excitable media — to make predictions in how arrhythmias appear in cardiac tissue. Spiral waves arise through a symmetry-breaking mechanism [44]. The formation of spirals has been considered in the presence of a spatial inho- mogeneity (which is how we model a fixed obstacle such as a patch of dead or damaged tissue). For instance, Panfilov and Keener [33] have shown that spirals can arise when a wavefront curls around regions of dead tissue, simi- lar to a water wave swirling around a rock. Lewis has proposed an alternate mechanism based on reflections generated as the wave passes through regions of damaged tissue [28]. However, as discussed above, spontaneous reentrant arrhythmias can ap- pear without warning in hearts that were not previously damaged and thus have no spatial inhomogeneities. How do arrhythmias arise in seemingly healthy hearts? Winfree has shown that a properly-timed external stimulus — modeling an impulse from an implanted electrode — in the neighbour- hood of a propagating wave can result in spirals [46]. This still leaves the question of how a spiral can arise spontaneously, without an external stim- ulus. 3.1.3 Mechanism for spontaneous functional reentry This chapter introduces a novel method for the spontaneous generation of reentrant activity in a homogeneous medium, modeling the spontaneous appearance of a reentrant arrhythmia in healthy tissue. We impose a tran- siently refractory region in front of a propagating wavefront. This concept is unique in that it is dependent only on the initial conditions of the system; the refractory region creates a functional barrier rather than an anatomical one, and does so without external intervention. Below we demonstrate and discuss the conditions for this mechanism to lead to reentry. It is robust for refractory regions of a variety of shapes and sizes. 49 3.1. Introduction: Modeling cardiac arrhythmias As well as being mathematically interesting, our mechanism is biologi- cally reasonable. It is indeed theoretically possible for this type of functional barrier to appear in healthy heart tissue. When a piece of excitable tissue spontaneously goes through a large excursion from rest it forces the local tissue to become refractory and propagates a wave of excitation — called an action potential — outwards. However, a wave will not propagate if the magnitude of the spontaneous stimulus is below a certain threshold or if the spontaneously-firing piece of tissue does not occupy a critical area size [48]. In these cases the local tissue will respond to the excitation sim- ply by becoming refractory. This can happen in the heart if tiny nodes of spontaneously-firing excitable tissue (such as the cells comprising the SA node — see Figure 3.1) exist in irregular locations within the conduction pathways of the heart [10]. Since such nodes are comprised of excitable car- diac tissue, they can be indiscernible from the surrounding tissue [20]. Our mechanism simulates an action potential — such as, but not necessarily, a regular heartbeat generated from the SA node — propagating towards one of these refractory regions. This is a physiologically plausible scenario. In the following section we set up our simulation in a 2D FHN system and then present numerical results, demonstrating that our scenario leads to reentry for a large region of initial condition space. We conclude with a critical discussion of the modeling assumptions and a comparison of the modelled results to order-of-magnitude spatial and temporal characteristics of an arrhythmia in cardiac tissue. The following relies heavily on certain terminology and properties of excitable media; before continuing, the reader is referred to the definitions and subtle properties within Appendix D. 50 3.2. Two-dimensional simulations 3.2 Two-dimensional simulations Following Glass and Josephson’s simple one-dimensional human heart tissue model in [15], we choose a specific two-dimensional FHN system: ∂u ∂t = D ( ∂2u ∂x2 + ∂2u ∂y2 ) + u(1− u)(u− a)− v, (3.1) ∂v ∂t = ǫ (u− bv) . where the diffusion coefficient D is a constant, 0 < ǫ << 1 is a small parameter, and a and b are constants. Values for the reaction term constants are taken directly from Glass and Josephson’s model as a = 0.139 and b = 2.54. Although Glass and Josephson chose ǫ = 8 ·10−3 and D = 1 ·10−3, we choose to decrease both parameters to ǫ = 4 · 10−3 and D = 2 · 10−4. The decrease emphasizes the asymptotic limit discussed throughout this thesis of ǫ → 0 and small D, resulting in simulations that display clearly defined regions of excitation and recovery that change on distinct time scales. Time and space units are milliseconds [10−3 s] and centimetres [10−2 m] (dimensional constants of unity are not shown). Simulations take place on a rectangular two-dimensional domain. We simulate an infinite domain by setting our finite domain to be considerably larger than the region of critical dynamics (which is approximately 2x2 cm2), with Neumann (absorbing) boundary conditions. Below we discuss the ex- pected effect of a finite domain close to the size of the critical dynamics. Details of the numerics can be found in Appendix A. We take for our initial conditions two-dimensional surfaces in the ex- citable (u) and the refractory (v) variables. See Figure 3.2. Every point in the domain is initially in a rest state (defined as u and v both zero) save for two exceptions: • A rectangular region near the centre of the domain that is assigned a positive refractory value of v = 1 (see Figure 3.2 (b)), and • An action potential (defined by a coupled profile of u and v values over a spatial interval of 10−1 cm) spanning the width of the domain at one 51 3.2. Two-dimensional simulations −0.2 0 0.2 0.4 0.6 0.8 (a) Excitable variable, u. 0 0.2 0.4 0.6 0.8 1 (b) Refractory variable, v. Figure 3.2: Sample initial conditions in u and v for 2D simulations. Variable magnitude is indicated by height. The length and width dimensions of the refractory rectangle in (b) are varied. Together, (a) u and (b) v define the initial state of the domain. 52 3.2. Two-dimensional simulations edge of the domain, propagating toward the refractory rectangle (see the peak and slight depression in u at the edge of Figure 3.2 (a) and the concurrent rise in v at the same edge of Figure 3.2 (b)). In the following discussion, we refer to the rectangular region with a positive initial refractory value simply as the ‘rectangle’. For the sake of comparison, all simulations maintain a constant initial distance d between the wave and the nearest edge of the rectangle. Since the refractory rectangle will eventu- ally decay to the rest state, we choose a distance small enough (with respect to the wave speed and the refractory decay rate) such that the rectangle is still absolutely refractory when the wave reaches it. Recall that a point defined as absolutely refractory is not capable of being excited by a stim- ulus such as an incoming wave. Length and width of the refractory region are varied. We expect reentry in the form of spiral waves to occur if the wavefront can wrap around the absolutely refractory region, re-enter it as it becomes relatively refractory (and capable of being excited), and propagate backwards through it. 3.2.1 Simulation results Numerical results mapping refractory rectangle size (length and width) to reentrant activity (reentrant activity versus no reentrant activity) are given in Figure 3.3. The independent parameters are the length and width of the rectangle. Width (in cm) denotes the side of the refractory rectangle parallel to the incoming planar wavefront. Length (in cm) denotes the side perpendicular to the wavefront. See Figure 3.4 (a) for an illustration of the orientation of the rectangle with respect to the wavefront. Reentrant activity occurs on, to the right of, and above the curve; no reentrant activ- ity occurs below and to the left of the curve. The upper right portion of the length-width plane corresponds to larger rectangles and the lower left corresponds to smaller rectangles. Thus, larger rectangles generally lead to reentry while smaller ones generally do not. This general result is ex- plained below. In addition, there is a minimum width required to support the back-propagation of excitation; observe the vertical section of the curve 53 3.2. Two-dimensional simulations 0 0.2 0.4 0.6 0.8 1 1.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 i ii reentry no reentry width [cm] le ng th [c m] Figure 3.3: Reentry activity dependent on refractory rectangle di- mensions in 2D simulations. Width (in cm) denotes the side of the refractory rectangle parallel to the incoming planar wavefront. Length (in cm) denotes the side perpendicular to the wavefront. Discretization is on the scale of 0.03 cm. Reentry occurs on, above, and right of the curve. No reentry occurs to the left and below the curve. Globally, sufficiently long rectangles of at least a minimum width of 0.42 cm lead to reentry. Suffi- ciently wide rectangles of nonzero length (within numerical accuracy) lead to reentry. Regions labelled i and ii indicate regions of the curve where the dependence of reentry on size does not hold locally. that asymptotes to a minimum width of 0.42 cm (for our choice of parame- ters). This width corresponds to the biological ‘critical electrode size’; any region of excitation with a dimension smaller than this will die out instead of propagating. To within numerical accuracy, beyond a certain width, any nonzero length leads to reentry. There are however two exceptions to the general shape of the curve: region i in Figure 3.3, at small width and intermediate length; and region ii in Figure 3.3, at small length and large width. In both these regions the reentry/no-reentry curve breaks its monotonic shape. In the following, we begin with a discussion of the general mechanism that allows larger 54 3.2. Two-dimensional simulations rectangles to support reentry and then move on to the subtleties of these two regions. This entire discussion is qualitative; any quantitative analysis, while expected to be quite difficult, is an open problem for future work. Mechanisms of reentry Reentry is dependent on the location of the wave when the refractory rect- angle has decayed sufficiently to allow excitation. This, in turn, depends on the relationship between three factors: the wave speed outside the refractory rectangle, the decay rate of the refractory variable, and the size of the re- fractory rectangle. The wave speed in a quiescent domain and the refractory decay rate are properties of the system, and thus these first two factors are the same in all our simulations. In all cases, the wave travels the distance d to the refractory region in some time t′ > 0. We have chosen the initial value in the refractory variable within the rectangle (v0 = v(t = 0)) to be sufficiently greater than v∗ such that the decayed refractory value at time t’ (v′ = v(t = t′)) is still absolutely refractory: v′ > v∗. Thus the wave of excitation cannot propagate through the refractory region when it reaches it at time t′. See the progression from t = 0 to t = t′ in Figure 3.4. When the wave reaches the rectangle, it breaks into two wavefronts that propagate around the rectangle (not shown). The deciding factor is the third factor that defines the location of the wave when the rectangle becomes relatively refractory: the dimensions of the refractory rectangle. Define t∗ as the time when the absolutely refractory region will just become relatively refractory. A sufficiently small rectangle allows the two waves to propagate fully around it, rejoin, and move past it by time t∗ (see Figure 3.5 (b)). The wave cannot back-propagate through the rectangle so reentrant activity does not occur. In contrast, the waves can only propagate partially around a large rectangle by time t∗. The broken wavefront is in a position to re-enter the region at two points (see Figure 3.5 (a)). By the time it has propagated backward through the rectangle, it can wrap around again (not shown). In this geometry, reentrant activity is induced in the form of a pair of spiral waves. See Figure 3.6 for a time 55 3.2. Two-dimensional simulations length width d absolutely refractory: v0>>v* excited wavefront t = 0 (a) Common schematic at time t = 0. (less) absolutely refractory: v0>v'>v* excited wavefront t = t' (b) Common schematic at time t = t′. Figure 3.4: Illustration of the common dynamics in the time interval [0, t′] for all 2D simulations. All simulations start with a propagating wavefront that travels the constant distance d to the refractory rectangle by time t = t′. Colours indicate variables: Red is excited, blue is refractory. The excited wavefront has a refractory tail. The refractory rectangle decays in time. 56 3.2. Two-dimensional simulations relatively refractory: v = v* t = t* excited wavefront excited wavefront (a) Reentry at time t = t∗ due to a large refractory rectangle. relatively refractory: v = v* t = t* excited wavefront (b) No reentry at time t = t∗ due to a small refractory rectangle. Figure 3.5: Illustration of the reentry and no reentry dynamics at time t = t∗ for large and small refractory rectangles. At time t = t∗ the refractory rectangle has decayed to the critical relatively refractory value v∗. Colours indicate variables: Red is excited, blue is refractory. The figure in (a) depicts how reentry occurs for a sufficiently large refractory rectangle. The wavefronts are still at the rectangle when it becomes relatively refrac- tory. The case of a sufficiently small rectangle is depicted in (b); the rejoined wavefront is already past the rectangle when it comes relatively refractory. 57 3.2. Two-dimensional simulations lapse of a simulation with a large enough rectangle to result in reentry. This behaviour holds for a wide range of dimensions, except for the parameter regions i and ii shown in Figure 3.3. Regions i and ii (the regions of hysteresis in Figure 3.3) experience a more subtle phenomenon. We first explore region i by choosing a fixed small width in region i (i.e., 0.4 cm) such that as we increase the length from zero, we experience (a) no reentry, (b) reentry, (c) no reentry. See figure Figure 3.4 for an illustration of the orientation of the rectangle and the incoming wave. Case (a): A region in case (a) behaves as in the general ‘sufficiently small rectangle’ case. No reentry can occur. Case (c): For sufficiently long rectangles at this width, the broken wave- front does not propagate fully around the rectangle by time t∗. The two fronts expand into the relatively refractory region through the sides and the facing components of the fronts annihilate each other as they collide. Any remaining back-propagating components cover less than the critical electrode size required for propagation. Hence, back-propagation does not occur. This describes a ‘minimum width’ phenomenon: rectangles of less than a minimum width do not permit back-propagation. Case (b): This intermediate case is particularly interesting because it is the exception to the minimum-width rule. It allows reentry because it has a slightly shorter rectangle than case (c) but not as short as that in (a). See the time-sequence snapshots of a simulation in Figure 3.7 showing the critical moments that lead to reentrant activity. At time t∗ the halves of the wavefront have reached the far corners of the relatively re- fractory rectangle [roughly the snapshot in Figure 3.7 (b)]. The fronts diffuse through the corners and become locally de-excited [Figure 3.7 (c); contrast this snapshot with the one showing the fronts propagating around (as opposed to through) the corners in Figure 3.6 (d)]. This also slows down the fronts. A sufficient amount of the super-threshold 58 3.2. Two-dimensional simulations (a) t = 10 ms. (b) t = 70 ms. (c) t = 130 ms. (d) t = 190 ms. (e) t = 250 ms. (f) t = 310 ms. Figure 3.6: Numerical simulation of a two-dimensional domain showing the development of reentrant activity due to an initially refractory region. Here we show the excitable variable u. Light regions are excited; dark values are depressed, and thus refractory. The domain is 2.4 cm by 2.4 cm. 59 3.2. Two-dimensional simulations (g) t = 370 ms. (h) t = 430 ms. (i) t = 490 ms. (j) t = 550 ms. (k) t = 610 ms. (l) t = 670 ms. Figure 3.6: (cont’d) Numerical simulation of a two-dimensional do- main showing the development of reentrant activity due to an initially refractory region. Here we show the excitable variable u. Light regions are excited; dark values are depressed, and thus refractory. The domain is 2.4 cm by 2.4 cm. 60 3.2. Two-dimensional simulations pulse remains to propagate around each corner; meanwhile, the lag from diffusing through the corners has allowed time for the rectangle to recover a bit more so the pulse can easily propagate backwards into it [Figure 3.7 (d),(e)]. Since the fronts have components that re-enter the rectangle parallel to each other, the effects of annihilation are less severe. A component of the backward-propagating excitation larger than the critical electrode size survives the re-entry into the rectangle [Figure 3.7 (f)] and reentrant activity is able to occur. Next, consider region ii in Figure 3.3. To explore region ii we choose a fixed large width (i.e., 1.06 cm) in region ii so that as we decrease the length to zero, we sequentially experience (a’) reentry, (b’) no reentry, (c’) reentry. These rectangles appear as wide refractory strips, with a high width to length ratio. Refer again to figure Figure 3.4 to see the orientation of the rectangle and the incoming wave. Case (a’): This case describes all sufficiently long rectangles at this width. It behaves as in the ‘sufficiently large rectangle’ case outline above. Reentry occurs as expected. Case (c’): In case (c’), the rectangle is as thin (length-wise) as possible at this width. The wavefront halves propagate around the thin refractory region just in time to re-enter it at time t∗. Although the stimulus is somewhat locally de-excited by the relatively refractory value, diffu- sion is sufficiently strong for the excitation to diffuse through the thin refractory layer. Back-propagation and reentry occur. Case (b’): The dynamics in the intermediate case are similar to those in (c’), except that the refractory layer is just thick enough (or diffusion is just not quite strong enough) to stop the pulse from diffusing through it. Reentry does not occur. 3.2.2 Discussion Our simulations indicate that large refractory obstacles (for example, all rectangles of at least a certain width) will lead to reentry and small ones 61 3.2. Two-dimensional simulations (a) t = 190 ms. (b) t = 220 ms. (c) t = 250 ms. (d) t = 280 ms. (e) t = 310 ms. (f) t = 340 ms. Figure 3.7: Numerical simulation of reentry in a two-dimensional domain due to an initially refractory rectangular region of small width and intermediate length [region i, case (b)]. Here we show the excitable variable u. Light regions are excited; dark values are depressed, and thus refractory. The portion of the domain shown is 1.8 cm by 1.8 cm. 62 3.2. Two-dimensional simulations will not. There are also some more subtle behaviours in regions i and ii for refractory areas with small widths and lengths. We now discuss the persis- tence of these results with respect to changes in the modeling machinery, the refractory region geometry, and the initial conditions. We conclude with a brief discussion of the biological plausibility of our proposed mechanism for reentry. Numerical simulation assumptions Finite domain: We chose to simulate an infinite domain. What if we used a finite do- main of width L and let the width of the refractory rectangle approach L? A refractory rectangle with a width that approaches the size of the domain would not allow the broken wavefront to propagate past the rectangle. The result would be neither propagation nor reentry, but rather a third type of behaviour that we can term wave block. The bifurcation diagram in Fig- ure 3.3 would have an additional vertical line at some width w∗ < L marking a third region resulting in wave block for widths w > w∗. Refractory region geometry : The geometry and composition of the refractory region — a single pos- itive refractory value over a rectangular region — was chosen for ease of numerical implementation. As a consequence, our refractory region had sharp corners and was sharply defined by a jump in the refractory variable. It is reasonable to ask if our results depend on either of these features. The geometry and refractory gradient in particular are of interest since a more biologically reasonable refractory region would have smooth corners and a continuous refractory gradient. We expect the general results to hold for a region with smooth corners: larger areas will still lead to reentry while smaller areas will not. The system behaviour in region ii (Figure 3.3) is due to the unique relationship between the large width and small length of the rectangle in this regime and there- fore the hysteresis in the bifurcation curve is expected to persist even if we 63 3.2. Two-dimensional simulations smooth out the sharp corners. It is unclear whether region i will persist if we smooth out the sharp corners. However, note that each of regions i and ii occur for some particular relationship between the length and the width of the domain. Thus, we expect that regions i and ii could be eliminated by implementing a circular refractory region instead of a rectangular one. This would reduce the dependence from the relationship between two parameters (length, width) down to a single parameter (either circumference or area). Then the general results would hold globally, with no exceptions: large re- gions would map to reentry and small regions would map to no reentry. The abruptness of the transition should also have no effect on the gen- eral idea of reentry as long as the refractory variable is a monotonically decreasing function of its distance from a central point. For example, con- sider the case of a 2D Gaussian-shaped refractory profile and compare it to the jump-type refractory profile we used in our simulations. The loca- tion of the bifurcation line we found for the jump-type profile may shift when we use a profile with a gradual gradient, but we would still expect large (small) regions to map to reentry (no reentry). A continuous but non-monotonic refractory gradient would likely introduce some additional subtleties, especially if it occurs in an irregularly-shaped refractory region. A non-monotonic refractory gradient in an irregularly-shaped region is a much more complex situation worthy of future study. Initial conditions: For our simulations we assumed a fixed initial distance d between the wavefront and the refractory rectangle. We anticipate that varying the dis- tance d would yield qualitatively similar results. As long as the propagating wavefront reaches the refractory region when it is still absolutely refractory (hence unexcitable), the symmetry of the wavefront is broken and reentry is possible. Similarly, the initial refractory value of the region must be chosen so that it is still absolutely refractory by the time the wavefront reaches it. Both of these conditions depend on the wave speed, which is a property of the material. They are also coupled. For example, increasing the distance requires increasing the initial refractory value. As long as the distance d and 64 3.2. Two-dimensional simulations the initial refractory value of the refractory region are chosen so that the wave reaches the region while it is still absolutely refractory, the dynamics of our mechanism with respect to these initial conditions are robust. Biological feasibility Note that our system parameters were chosen to match Glass and Joseph- son’s one-dimensional cardiac tissue model in [15]. Since Glass and Joseph- son chose parameters within the range of acceptable values for a human heart, we expect our system’s spatial and temporal properties to map closely to healthy cardiac tissue. With this in mind, we examine our results criti- cally below with respect to what we would expect to see in cardiac tissue. Spatial scale: We first consider the size of a human heart. The refractory regions that map to reentry have areas of 0.1 cm2 and up. If the areas were on the or- der of square metres, our mechanism would be suspect; areas on the order of cm2 are not indications of implausibility. Next, we consider the natural size of refractory regions in cardiac tissue. Theoretically, a number of sub threshold excitations in a concentrated area (as described in Section 3.1) of any size can lead to a transiently refractory region. Comparison with experimental data is necessary to make any actual claims about whether or not this is plausible in a human heart. Temporal scale: Next, we look at the spiral waves generated and displayed in Figure 3.6. The spirals have been fully established by t = 430 ms, displayed in (h). Im- age (l) taken 240 ms later reveals a nearly identical situation to (h). Due to the known properties of spiral waves, we expect the spirals to persist with this same frequency. Monitoring the excitation activity at one particular point near the edge of the domain would return a pulse approximately every ∼240 ms. In terms of frequency, this translates to roughly 250 pulses per minute. Normal heartbeats have a frequency of 60 - 100 beats per minute. 65 3.3. Conclusions and future work Tachycardia is defined as a relatively rapid heartbeat. Therefore a frequency of 250 pulses per minute is within a reasonable magnitude for a tachycardia caused by reentrant activity. 3.3 Conclusions and future work We conclude that our general mechanism of reentry in a two-dimensional FHN system is robust. It does not depend on the sharp corners of the refrac- tory region, nor on the abruptness of the transition from rest to refractory, nor on the specific arrangement of initial conditions. Moreover, our simple order-of-magnitude analysis does not render it biologically implausible. However, we must remember that the FitzHugh Nagumo model we’ve implemented is only a simple phenomenological model of excitable media that has been adapted for cardiac tissue. As such, any actual statements about biological applicability require this mechanism to be considered in more accurate and specific models. For example, there exist ionic models which capture detailed properties of cardiac tissue [3] [29] [30] [6], and higher- dimensional models that consider the actual geometry of the heart [35]. The implementation and study of this mechanism in an accurate cardiac model is the first suggested direction for future work. Biological applications aside, this chapter is nonetheless an interest- ing demonstration of a spiral-generating mechanism in the FHN system. The second suggested direction for future work is the analytical study of this mechanism in the FHN model. In particular, the concepts of one- dimensional wave propagation studied in Chapter 2 may be applicable to the study of a wavefront approaching the refractory region along a selected one-dimensional contour of this two-dimensional domain. However, it may be too difficult to perform any analytical work that will explain these results more quantitatively. 66 Bibliography [1] Rubin R. Aliev and Alexander V. Panfilov. A simple two-variable model of cardiac excitation. Chaos, Solitons and Fractals, 7(3):293 – 301, 1996. [2] Uri M. Ascher, Steven J. Ruuth, and Brian T. R. Wetton. Implicit- explicit methods for time-dependent partial differential equations. SIAM Journal on Numerical Analysis, 32(3):pp. 797–823, 1995. [3] G. W. Beeler and H. Reuter. Reconstruction of the action potential of ventricular myocardial fibres. The Journal of Physiology, 268(1):177– 210, 1977. [4] Andrew J. Bernoff. Spiral wave solutions for reaction-diffusion equa- tions in a fast reaction / slow diffusion limit. Physica D, 53:125–150, 1991. [5] David J. Christini and Leon Glass. Introduction: Mapping and control of complex cardiac arrhythmias. Chaos, 12(3):732–739, 2002. [6] R.H. Clayton, O. Bernus, E.M. Cherry, H. Dierckx, F.H. Fenton, L. Mirabella, A.V. Panfilov, F.B. Sachse, G. Seemann, and H. Zhang. Models of cardiac tissue electrophysiology: Progress, challenges and open questions. Progress in Biophysics and Molecular Biology, 104:22– 48, 2011. [7] J. Crank and P. Nicolson. A practical method for numerical evalua- tion of solutions of partial differential equations of the heat-conduction type. Advances in Computational Mathematics, 6:207–226, 1996. 10.1007/BF02127704. 67 Chapter 3. Bibliography [8] Eric N. Cytrynbaum. Using Low Dimensional Models to Understand Cardiac Arrhythmias. PhD thesis, The University of Utah, December 2001. University of Utah. [9] Jorge M. Davidenko, Arcady V. Pertsov, Remy Salomonsz, William Baxter, and Jose Jalife. Stationary and drifting spiral waves of excita- tion in isolated cardiac muscle. Nature, 355:349–351, January 1992. [10] Dale Dubin. Rapid Interpretation of EKG’s. COVER Publishing Com- pany, 6th edition, 2000. [11] Christopher P Fall. Computational Cell Biology, chapter 7, pages 189– 194. Springer, 2002. [12] F. H. Fenton, E. M. Cherry, H. M. Hastings, and S. J. Evans. Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity. Chaos, 12:852–892, September 2002. [13] P. C. Fife and J. B. McLeod. The approach of solutions of nonlinear diffusion equations to travelling front solutions. Arch. Rat. Mech. Anal., 65:335–361, 1977. [14] R. Fitzhugh. Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal, 1:445 – 466, 1961. [15] Leon Glass and Mark E. Josephson. Resetting and annihilation of reentrant abnormally rapid heartbeat. Physical Review Letters, 75(10):2059–2062, 1995. [16] Dan Henry. Geometric Theory of Semilinear Parabolic Equations, chap- ter 5, pages 128–130. Springer-Verlag, 1981. [17] A.L. Hodgkin and A.F. Huxley. A quantitative description of membrane current and its application to conduction and excitation in nerve. Jour- nal of Physiology, 117(4):500–544, August 1952. 68 Chapter 3. Bibliography [18] Christopher K. R. T. Jones. Stability of the travelling wave solution of the fitzhugh-nagumo system. Transactions of the American Mathemat- ical Society, 286(2):431–469, December 1984. [19] Kamyar Kamjoo, Rakumi Uchida, Takanori Ikeda, Michael C. Fishbein, Alan Garfinkel, James N. Weiss, Hrayr S. Karagueuzian, and Peng- Sheng Chen. Importance of location and timing of electrical stimuli in terminating sustained functional reentry in isolated swine ventricular tissues. Circulation, 96:2048–2060, 1997. [20] Hart Katz. Spontaneously-arising arrhythmias. Personal communica- tions, September 2010. [21] J. P. Keener. Waves in excitable media. SIAM Journal of Applied Math, 39(3):528–548, December 1980. [22] James P. Keener. The topology of defibrillation. Journal of Theoretical Biology, 230(4):459 – 473, 2004. Special Issue in honour of Arthur T. Winfree. [23] James P. Keener and James Sneyd. Mathematical Physiology, chapter 9, pages 268–295. Springer-Verlag, 1998. [24] James P. Keener and John J. Tyson. Spiral waves in the belousov- zhabotinskii reaction. Physica D: Nonlinear Phenomena, 21(2-3):307 – 324, 1986. [25] J.P. Keener. A geometrical theory for spiral waves in excitable media. SIAM J. Appl. Math, 46:1039–1059, 1986. [26] J.P. Keener. Homogenization and propagation in the bistable equation. Physica D, 136:1–17, 2000. [27] Trine Krogh-Madsen and David J. Christini. Resetting and termination of reentry in a loop-and-tail cardiac model. Phys. Rev. E, 77:011916, Jan 2008. 69 Chapter 3. Bibliography [28] Tim Lewis. The effects of nonexcitable regions on signal propagation in excitable media: Propagation failure and reflection. PhD thesis, The University of Utah, December 1998. [29] CH Luo and Y Rudy. A model of the ventricular cardiac action poten- tial. depolarization, repolarization, and their interaction. Circulation Research, 68(6):1501–1526, 1991. [30] C. Morris and H. Lecar. Voltage oscillations in the barnacle giant muscle fiber. Biophysical Journal, 35(1):193 – 213, 1981. [31] J. Nagumo, S. Arimoto, and S. Yoshizawa. An active pulse transmission line simulating nerve axon. Proceedings of the IRE, 50(10):2061–2070, October 1962. [32] A. Panfilov and J.P. Keener. Re-entry in an anatomical model of the heart. Chaos, Solitons and Fractals, 5(3-4):681 – 689, 1995. Nonlinear Phenomena in Excitable Physiological Systems. [33] Alexander V. Panfilov and James P. Keener. Effects of high frequency stimulation on cardiac tissue with an inexcitable obstacle. Journal of Theoretical Biology, 163(4):439 – 448, 1993. [34] A.M. Pertsov, E.A. Ermakova, and A.V. Panfilov. Rotating spiral waves in a modified fitz-hugh-nagumo model. Physica D: Nonlinear Phenom- ena, 14(1):117 – 124, 1984. [35] Zhilin Qu, Jong Kil, Fagen Xie, Alan Garfinkel, and James N. Weiss. Scroll wave dynamics in a three-dimensional cardiac tissue model: Roles of restitution, thickness, and fiber rotation. Biophysical Jour- nal, 78(6):2761 – 2775, 2000. [36] John Rinzel and Joseph Keller. Traveling wave solutions of a nerve conduction equation. Biophysical Journal, 13:1313–1337, 1973. [37] B M Sager. Propagation of traveling waves in excitable media. Genes and Development, 10(18):2237–2250, 1996. 70 Chapter 3. Bibliography [38] Faramarz H. Samie and Jose Jalife. Mechanisms underlying ventricular tachycardia and its transition to ventricular fibrillation in the struc- turally normal heart. Cardiovascular Research, 50:242–250, 2001. [39] B. Sandstede and A. Scheel. Superspiral structures of meandering and drifting spiral waves. Phys. Rev. Lett., 86:171–174, Jan 2001. [40] B. Sandstede, A. Scheel, and C. Wulff. Bifurcations and dynam- ics of spiral waves. Journal of Nonlinear Science, 9:439–478, 1999. 10.1007/s003329900076. [41] T. K. Shajahan, Sitabhra Sinha, and Rahul Pandit. Spiral-wave dy- namics depend sensitively on inhomogeneities in mathematical models of ventricular tissue. Phys. Rev. E, 75:011929, Jan 2007. [42] John J. Tyson and James P. Keener. Spiral waves in a model of my- ocardium. Physica D: Nonlinear Phenomena, 29(1-2):215 – 222, 1987. [43] B. N. Vasiev, P. Hogeweg, and A. V. Panfilov. Simulation of Dic- tyostelium Discoideum aggregation via reaction-diffusion model. Phys. Rev. Lett., 73:3173–3176, Dec 1994. [44] Arthur T. Winfree. Spiral waves of chemical activity. Science, 175(4022):634–636, 1972. [45] Arthur T. Winfree. Sudden cardiac death: A problem in topology. Scientific American, 248:114–161, 1983. [46] A.T. Winfree. Electrical instability in cardiac muscle: Phase singu- larities and rotors. Journal of Theoretical Biology, 138(3):353 – 405, 1989. [47] Francis X. Witkowski, L. Joshua Leon, Patricia A. Penkoske, Wayne R. Giles, Mark L. Spano, William L. Ditto, and Arthur T. Winfree. Spa- tiotemporal evolution of ventricular fibrillation. Nature, 392:78–82, March 1998. 71 Chapter 3. Bibliography [48] Vladimir Sergeevich Zykov. Simulation of wave processes in excitable media, chapter 5, page 117. Manchester University Press ND, 1987. 72 Appendix A Numerical details All numerical simulations of the FHN system are performed in the commer- cial computing environment of MATLAB. Numerical integration over time is performed using an implicit-explicit method [2] comprised of a semi-implicit Crank-Nicolson scheme for the diffusion term and an explicit Euler scheme for the reaction terms. The Crank-Nicolson scheme for the diffusion term is an unconditionally stable numerical scheme that is highly accurate for a small enough ratio of Ddt/(dx)2 [7], where dt is the time step and dx (or dy) is the grid spacing. A sufficiently small dt is required to ensure accuracy and stability with the explicit Euler scheme. In Chapter 2 (one-dimensional simulations), nondimensional grid spac- ings and time steps are dx = 0.0025 and dt = 0.0001 respectively. In Chapter 3 (two-dimensional simulations), I use grid spacings of dx = dy = 0.03 cm and a time step of dt = 0.1 ms for all 2D simulations. To check the accuracy of my simulations, I replicated some of my results using grid spacings 1/3 the size of my chosen dx and dy and a time step an order of magnitude smaller than my chosen dt (these simulations are not shown). The results presented in Chapter 3 are qualitatively robust to these tested refinements. 73 Appendix B Determination of collapse or propagation in Section 2.4.1 The condition for determining collapse or propagation relies on (a) the rel- ative position of the wave front at each step, and (b) the speed of the wave. Simulations with an initial exponential refractory profile show that if a front travels backward even one grid step, it leads to full collapse. Thus the problem of categorizing collapse can be simplified to the problem of flagging a single backwards movement. The analysis in Section 2.5 reveals that the speed of a propagating front in the SFHN system asymptotes to a constant value for exponential re- fractory profiles of the form (2.35). Numerical comparisons indicate that a propagating front in the FHN system also exhibits this behaviour (see front trajectories in Figure 2.9). Once a front reaches a constant speed, we conclude it will continue propagating. The domain size and simulation time are chosen sufficiently large in each case to conclusively determine the final behaviour of the wavefront. The exact algorithm for determining collapse or propagation is as follows, beginning after the first iteration of the system in time: 1. Calculate the position of the front as defined by (2.29). Linear inter- polation is used to approximate the front position if it falls between two points on the grid. 2. Calculate the difference between the current position and the previous position. (a) A negative difference indicates the front has moved backwards. 74 Appendix B. Determination of collapse or propagation in Section 2.4.1 The initial condition set is categorized as mapping to ‘collapse’. (b) If the difference is positive or zero, we iterate to the next time step and return to step 1. 3. If the front reaches a constant speed, then the front is categorized as ‘propagating’. This method is used to map each initial condition to one of the two states. 75 Appendix C Singular derivation of dmin in Section 2.4.2 Fitting a square wave of the form vstep = { 0, 0 < x < xstep h, x > xstep (C.1) under an initial exponential refractory profile of the form v(x, 0) = Aeλx (C.2) requires the following relationship between the height of the square wave and its location: h = v(xstep, 0), (C.3) h = Aeλxstep . Three such square waves that fit under a given exponential are pictured in Figure 2.5 in Section 2.4.2. We want to calculate the smallest distance xstep at which a step wave of height given by (C.3) will successfully stall a front. This will be our singularly-derived dmin. We use the fact that stall will occur in the singu- lar case if the decayed step wave is of magnitude vz or greater when the wavefront reaches it. The refractory variable solution (2.36) from property 1 can be used to find the decayed height of a fitted square wave when the front reaches xstep 76 Appendix C. Singular derivation of dmin in Section 2.4.2 at τstep: v(xstep, τstep) = v(xstep, 0)e −(1+b)τstep . (C.4) The time it takes a front to reach xstep is simply the distance traveled divided by the speed: τstep = xstep c0(0) , (C.5) where c0(0) is the explicitly–defined speed given in (2.30) evaluated at v = 0: c0(0) = 1− 2a√ a(1− a) , (C.6) This constant front speed is valid for the entire interval preceding the square wave because the v value will not change ((u, v) = (0, 0) is a steady state). The condition for propagation failure is v(xstep, τstep) = vz, (C.7) where vz has been calculated as vz = 1 2 − a. (C.8) Substituting the expression from (C.4) into the left-hand side of (C.7) gives v(xstep, 0)e −(1+b)τstep = vz, Aeλxstepe−(1+b)τstep = vz. (C.9) Now we use (C.5) in the left-hand-side of (C.9) and the expression for vz (C.8) in the right-hand side to get A exp [ xstep ( λ− (1 + b) √ a(1− a)/(1− 2a) )] = 1 2 − a. (C.10) The xstep that satisfies this equation defines the endpoint of our interval of 77 Appendix C. Singular derivation of dmin in Section 2.4.2 assured propagation, dmin: dmin = ln(1/2−aA ) λ− (b+1) √ a(1−a) 1−2a (C.11) 78 Appendix D Excitable media terminology and properties Some introduction to excitable media is provided in Chapter 1. However, simulation descriptions and discussions in Chapter 3 rely heavily on the ter- minology of excitable media. Below we provide definitions of the additional necessary terms. In addition, explanations of the more subtle reentrant phe- nomena in Chapter 3 require knowledge of several properties of excitable stimuli, also stated below. D.0.1 Terminology The state of a medium is defined by its refractory value, v, and its excitable value, u. We describe the possible states as follows: • A medium that has no significant refractory or excited value is said to be quiescent or at rest. In this state, the medium has the capability to propagate excitation if it receives a stimulus greater than or equal to some critical threshold. • If a medium has a positive refractory value and cannot be excited by a stimulus of any magnitude, it is defined as absolutely refractory. These regions are considered not excitable. A positive refractory value couples with a negative excitable value. • A region that has a positive refractory value but has the capability to propagate excitation if given a stimulus sufficiently larger than some threshold is relatively refractory. A relatively refractory medium has a lower refractory value than an absolutely refractory medium. 79 Appendix D. Excitable media terminology and properties • An excited medium is one that is mid- large excursion from the rest state. It has a positive excitable value. The term excitable encompasses both relatively refractory and quiescent media. An unexcitable medium is called depressed. The term refractory indicates a region that is either absolutely or relatively refractory. The crit- ical refractory value marking the boundary between absolutely and relatively refractory (unexcitable and excitable) is defined as v∗. 5 D.0.2 Properties of stimuli To propagate a wave of excitation, an excited stimulus must be greater than or equal to some critical magnitude. This is the so-called ‘threshold’ behaviour of excitable media. In addition, the super-threshold or above threshold stimulus must occupy a space greater than or equal to some critical area [48]. In cardiac tissue, this is often referred to as a critical electrode size. We also need to consider the coupled behaviour of the excitable and refractory variables. A region with a relatively refractory value will locally decrease the excitation value of any stimulus passing through it. However, excitation diffuses. Even if a pulse is de-excited locally by a high refractory value, the pulse may diffuse through a sufficiently thin relatively refractory region. Furthermore, a pulse travels slower through a relatively refractory region than it would through a quiescent region. When two traveling pulses collide in one spatial dimension, they an- nihilate each other. Each pulse runs into the absolutely refractory tail of the other and cannot propagate. This behaviour extends to higher spatial dimensions. 5The critical refractory value defined as v∗ in Chapter 3 is the lower limit of the zero- speed refractory value vz defined in Chapter 2. 80
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A study of wave propagation in the FitzHugh Nagumo...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A study of wave propagation in the FitzHugh Nagumo system Paton, Kelly Marie 2011
pdf
Page Metadata
Item Metadata
Title | A study of wave propagation in the FitzHugh Nagumo system |
Creator |
Paton, Kelly Marie |
Publisher | University of British Columbia |
Date Issued | 2011 |
Description | An excitable medium has two key properties: a sufficiently large stimulus provokes an even bigger response (excitability), and immediately following a stimulus the medium cannot be excited (refractoriness). A large class of biological systems from cardiac tissue to slime mold are examples of excitable media. FitzHugh Nagumo (FHN) is the canonical model of excitable media. Its two variables are the state of excitation and refractoriness of the one- or two-dimensional medium. Although one of the simplest models, FHN exhibits complex dynamics that have not been fully explored. For example, it supports a stable traveling pulse solution. However, this pulse can be destabilized by large perturbations. In Chapter 2 I explore a one-dimensional example where the perturbation is an increasing refractory profile. This perturbation can lead to collapse of the pulse depending on the steepness of the profile, as conjectured by Keener [Keener, J. (2004) J Theo. Bio. 230(4):459-73]. In Chapter 3 I consider a perturbation in two dimensions which can cause the stable traveling pulse to wrap around the perturbation and generate self-sustaining spiral activity. The one-dimensional example for exponential refractory profiles is explored numerically for a piecewise linear FHN system. Steep profiles lead to collapse while milder profiles allow propagation. The exponential profiles are used as bounds for more general profiles to predict where collapse and propagation will occur. I also make use of a singular FHN system in the limit ε→ 0 to provide insight into the behaviours of the full FHN system for small ε and small diffusion. I conclude this chapter by showing analytically that, in contrast to the full system, a wave in the singular system will propagate for any exponential refractory profile. The two-dimensional case is explored numerically in a FHN system. The use of a temporarily refractory region as a perturbation is a novel mechanism for generating spiral activity. Moreover, it is shown to be robust for refractory regions of a large area. This situation models the appearance of abnormal electrical activity in the heart. In particular, it models the appearance of abnormal electricity activity in undamaged cardiac tissue. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2012-01-09 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0072522 |
URI | http://hdl.handle.net/2429/39955 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2012-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2012_spring_paton_kelly.pdf [ 2.18MB ]
- Metadata
- JSON: 24-1.0072522.json
- JSON-LD: 24-1.0072522-ld.json
- RDF/XML (Pretty): 24-1.0072522-rdf.xml
- RDF/JSON: 24-1.0072522-rdf.json
- Turtle: 24-1.0072522-turtle.txt
- N-Triples: 24-1.0072522-rdf-ntriples.txt
- Original Record: 24-1.0072522-source.json
- Full Text
- 24-1.0072522-fulltext.txt
- Citation
- 24-1.0072522.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:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0072522/manifest