O S C I L L A T O R Y W A V E S IN A N I N H O M O G E N E O U S EXCITABLE MEDIUM by A L A I N PRAT B.Sc, The University of British Columbia, 2002 A THESIS SUBMITTED IN PARTIAL F U L F I L M E N T OF THE REQUIREMENTS FOR THE D E G R E E OF M A S T E R OF SCIENCE in THE F A C U L T Y OF G R A D U A T E STUDIES ( Mathematics) THE UNIVERSITY OF BRITISH C O L U M B I A August 2005 © Alain Prat, 2005 Abstract A n excitable medium is a nonlinear dissipative dynamical system able to sustain undamped wave propagation. Examples of excitable media include cardiac tissue and neural media. This thesis investigates a class of novel dynamical phenomena in excitable media, known collectively as oscillatory waves. The examples of oscillatory waves considered here are Tango waves, breathers and pulse generators. Tango waves, also known as oscillatory fronts, are characterized by a back-and-forth motion of a wavefront. Breathers, also known as oscillatory pulses, are characterized by an expanding and contracting motion of a pulse-shaped pattern. Pulse generators are breathers that periodically emit traveling pulses. These oscillatory waves are caused by the presence of either a timeindependent forcing or spatial variation in a model parameter. The main point of this thesis is that oscillatory fronts or pulses can be understood as emerging from the Hopf bifurcation of a stationary front or pulse, respectively. The secondary point of this thesis is to investigate secondary bifurcations connected to the Hopf bifurcation. Possible applications of oscillatory waves will also be discussed. Analytical and numerical methods are used to compute the equilibrium solutions and the spectrum of their linearizations for piecewise linear and cubic Fitzhugh-Nagumo models of excitable media. Combining this with numerical simulations allows us to link oscillatory waves with a Hopf bifurcation. A center manifold reduction allows us to also connect the oscillatory waves to a generalized Hopf bifurcation. ii Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures vi Acknowledgements Chapter 1. vii Introduction 1 1 Excitable Media 2 Two Examples of Excitable Media 2.1 Cardiac Tissue 2.2 Action Potential Propagation 3 Oscillatory Waves 4 Related Works 4.1 Calcium Waves 4.2 Neural Media 4.3 Belousov-Zhabotinsky Reaction 4.4 The Complex Ginzburg-Landau Equation Chapter 2. Equilibria 11 1 Piecewise Linear Model 1.1 Stationary Fronts 1.2 Stationary Pulses 2 Cubic Model 2.1 Stationary Fronts 2.2 Stationary Pulses Chapter 3. 11 12 15 19 19 21 Linear Analysis 1 Piecewise Linear Model 1.1 Eigenvalues for Stationary 1.2 Eigenvalues for Stationary 2 Cubic Model 2.1 Eigenvalues for Stationary 2.2 Eigenvalues for Stationary 1 2 3 4 5 6 6 7 9 10 26 Fronts Pulses Fronts Pulses iii 29 29 32 35 36 36 Table of Contents Chapter 4. Numerical Simulations 39 1 Oscillatory Fronts 2 Oscillatory Pulses 2.1 Piecewise Linear Model 2.2 Cubic Model 39 40 40 42 Chapter 5. Bifurcation Analysis 48 Chapter 6. Summary and Discussion 55 Bibliography 57 Appendix A . M A P L E Programs 1 Stationary fronts 2 Eigenvalues for stationary fronts 3 Stationary pulses 4 Eigenvalues for stationary pulses 5 First Liapunov coefficient 6 Domain sizes used 59 59 60 60 60 61 62 Appendix B . 63 Numerical Methods iv List of Tables 3.1 Possible local bifurcation of equilibria 29 A . l Domain sizes used in M A P L E calculations 62 B. l Domain size and mesh points for numerical simulations 63 v List of Figures 1.1 Cubic nonlinearity and piecewise linearity for F - N model 1.2 Traveling front and back 1.3 Traveling pulse 2.1 2.2 2.3 2.4 2.5 2.6 Stationary front in the piecewise linear model Cubic and quadratic bifurcation scenarios in the piecewise linear model Stationary front in the cubic model Cusp bifurcation in the cubic model Quadratic and cubic bifurcation scenarios in the cubic model Example of a stationary pulse in the cubic model 2 3 4 16 18 20 22 24 25 3.1 Stability diagram for fronts in e-a parameter space for the piecewise linear and cubic models 32 3.2 Stability diagrams for pulses in e-h parameter space (piecewise linear model) . . . 34 3.3 Bifurcation diagrams for stationary pusles in the piecewise linear model 35 3.4 Stability diagrams for stationary pulses in e-h parameter space (cubic model) .. 37 3.5 Bifurcation diagrams for quadratic and cubic scenarios in the cubic model 38 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 Supercritical Hopf bifurcation of a front (piecewise linear model) Super and subcritical Hopf bifurcation of a pulse (piecewise linear model) Snapshots of a pulse generator Bifurcation diagrams for stationary and oscillatory pulses (cubic model) Compex oscillations for the pulse generator Diverging period of oscillations of pulse generator Oscillatory front and pulse Pulse generators 5.1 Example of a generalized Hopf bifurcation of pulses 5.2 Example of a generalized Hopf bifurcation of pulses 5.3 First Liapunov coefficient for Hopf bifurcation of a front vi 41 42 43 44 45 46 46 47 52 53 54 Acknowledgements I would like to thank my thesis supervisor Dr. Wayne Nagata for being an excellent teacher, and for his patience and support. I would also like to thank Dr. Yue-Xian L i for introducing me to his "Tango wave" problem, which led to the development of the topics in this thesis. Thanks also to Dr. Paul Bressloff of the University of Utah for a useful collaboration on this project. Chapters 2, 3 and 4 in this thesis result from joint work with Dr. Yue-Xian L i . This work has been supported by an N S E R C P G S A scholarship. vii Chapter 1 Introduction 1 Excitable M e d i a One of the most commonly observed dynamical phenomena in nature is that of wave propagation. Examples of propagating waves are sound waves in air, water waves on the surface of the ocean, and waves of electrical excitation in the heart. In these and many other examples, the waves propagate through a medium, and the properties of the medium determine how the wave evolves in time. In the presence of dissipation, the medium sustaining wave propagation can be said to be either passive or excitable. In passive media, the amplitude of the propagating wave is gradually damped in time. An example of passive wave propagation is sound waves passing through air. In excitable media, on the other hand, the amplitude of the wave remains constant, even though dissipation is present. An example of an excitable medium is the heart, where waves of electrical activity propagate from the top of the heart and spread downwards, causing the heart to contract. One of the simplest models of excitable media is the FitzHugh-Nagumo model: dv = f(y) -w + z(x) + W v dt dw = bv — w ~dt el 2 1 (1.1) (1.2) Chapter 1. Introduction f(v) f(v) (a) (b) Figure 1.1: a) Cubic nonlinearity f(v) = v(l — v)(v — a) b) Piecewise linearity f(v) = —v + H(v — a), used as an approximation of the cubic nonlinearity. Throughout this thesis, we take f(v) as in a) or b) above. where v(x, t) and w(x, t) are defined for (x, t) 6 W x R , with n = 1, 2 or 3. e > 0 and 1 + b > 0 are constants, f(v) is a cubic or piecewise linear function as shown in fig. 1.1 and z(x) represents a stationary inhomogeneous input. Examples of traveling wave solutions of the above are traveling fronts, backs and pulses. These are shown in figures 1.2a, 1.2b and 1.3 respectively. In two or three dimensions, one can also have spiral waves, target waves and scroll waves. 2 Two Examples of Excitable M e d i a We briefly discuss two examples of wave propagation in excitable media: waves of excitation in cardiac tissue and action potential propagation in axons (i.e. nerve fibers). 2 Chapter 1. Introduction v(x,t) v (a) (x,t) (b) Figure 1.2: a) Traveling front b) Traveling back 2.1 Cardiac Tissue During each heart beat, a wave of excitation passes through the heart, compressing first the atria which pushes blood into the ventricles, and then compressing the ventricles which push blood into the body. This wave of excitation is a traveling pulse with a velocity of about 0.5 m/s and a characteristic size of about 10 cm. It is conducted by electrical coupling between neighbouring cells of the heart muscle. Every cell in the cardiac muscle can become electrically excited when a perturbation exceeds a certain threshold. Such a perturbation alters the membrane potential to about 70 mV, after which the potential returns to its equilibrium of about -90 mV. The coupling between cells causes the excitation to rapidly spread throughout the cardiac tissue. 3 Chapter 1. Introduction v(x,t) X Figure 1.3: Traveling pulse 2.2 Action Potential Propagation Intricate communication between nerve cells takes place by transmitting signals, known as action potentials, along a long slender projection of a nerve cell, called an axon. Axons are, in effect, the primary transmission lines of the nervous system. They are about one micrometer across and can sometimes be up to a meter or more in length. Action potentials are abrupt, short-lived reversals in the electrical charge of an axon due to the opening and closing of voltage-gated ion channels. They propagate as traveling pulses along the axon, with speeds of up to 120 meters per second. In this context, the axon membrane is the excitable medium. The energy sustaining the wave propagation is stored in the ion channels. After the wave passes, there is a refractory phase, where the potential energy is restored using ion pumps. 4 Chapter 1. Introduction 3 Oscillatory Waves This thesis investigates a recently discovered class of novel dynamical phenomena in excitable media, called oscillatory waves. The examples of oscillatory waves discussed here are Tango waves, breathing patterns and pulse generators. The Tango wave, also known as an oscillatory front, is a wave which moves in a back-and-forth pattern, alternating between a front and a back (see fig. 4.7a). The breathing pattern, also refered to as an oscillatory pulse, consists of a unimodal pulse shaped pattern whose amplitude oscillates in time (see fig. 4.7b). The pulse generator is a breathing pattern which periodically emits traveling pulses (see figs. 4.3 and 4.8). The first mentions of Tango waves appear in [15] and [12]. The first mentions of breathing patterns and pulse generators in excitable media appear in [4] and [16]. Oscillatory waves are caused by the inhomogeneous term z(x) in (1.1)-(1.2). Tango waves are caused by a monotonic z(x), such as a linear ramp. Breathing patterns and pulse generators are caused by a unimodal z(x), such as a Gaussian. Note that throughout this thesis, we are considering (1.1)-(1.2) in one space dimension. The main point of this thesis is that oscillatory waves can be traced back to a Hopf bifurcation of an equilibrium. A linear z(x) with a large slope stabilizes a front, and a reduction of the slope releases the front through a Hopf bifurcation. Similarly, a Gaussian z(x) with a large amplitude stabilizes a pulse, and a reduction in amplitude induces a Hopf bifurcation. A secondary point of this thesis is to investigate bifurcations related to the Hopf bifurcation, including a generalized Hopf bifurcation. The remaining chapters of this thesis are organized as follows. In chapter 2, we find equilibrium solutions. These equilibria are then used in chapter 3 as the starting point for a linearization of the equations. Computing the spectrum of the eigenmodes of the 5 Chapter 1. Introduction linearized equations reveals that for certain parameter values, we have a pair of purely imaginary eigenvalues, signaling the prescence of a Hopf bifurcation. In chapter 4, we link the Hopf bifurcation found in the previous chapter with oscillatory waves found in numerical simulations. A detailed look at numerical simulations suggest that the oscillatory waves are also linked to a generalized Hopf bifurcation. This point is studied in the next chapter by performing a bifurcation analysis using center manifold reduction. Chapter 6 is a summary and discussion of the thesis. 4 Related Works 4.1 Calcium Waves The work in this thesis was initially motivated by a recent study of calcium waves observed in mature egg cells when activated by a sperm at fertilization, called fertilization calcium waves. Numerical simulations of a realistic model of fertilization calcium waves revealed the occurrence of an interesting wave phenomenon, originally called "Tango waves", when the medium is intrinsically inhomogeneous [12]. Tango waves are characterized by a back-and-forth pattern of front propagation. They occur in a bidomain model of traveling calcium fronts observed in mature egg cells triggered by the injection of a large amount of calcium. The equations of this model, in a scaled form, are given by I• dh — C = g[l-(l + c)h], . . (1.5) where c is the scaled calcium concentration in the cytosol and z is a weighted sum of the calcium concentrations in the cytosol and in the calcium stores, called the total calcium. h is a variable describing how the opening of the channels through which the stores release Chapter 1. Introduction calcium is regulated by c. See [12] for more detailed description and study of the model. Because calcium diffuses slowly in the intracellular medium due to the existence of calcium buffers, 8 is typically very small. As a result, the injection of a large amount of calcium at a localized site disperses slowly. This creates an inhomogeneous distribution of the total calcium z in the cellular medium that has a Gaussian shaped spatial distribution. This distribution changes very slowly because 5 << 1 which implies that z is a slow variable. Results in [12] showed that the back-and-forth movement of Tango waves occurs in spatial domains in which the slope of the inhomogeneity is within an appropriate range. Since z varies slowly in time, we can replace z(x, t) by a time-independent term z(x), resulting in the following reaction-diffusion type of equations containing a spatially distributed parameter. ^ (1.7) = Tango waves also occur in this simpler model. Both Gaussian and linear profiles of z(x) have been shown capable of generating Tango waves in this model. It was demonstrated that Tango waves also occur in the F-N model in the presence of an inhomogeneity z(x) that is either Gaussian or linear [12]. These studies bring up an interesting theoretical question: Why are Tango waves related to space inhomogeneities? In this thesis, we answer this question by showing that introducing inhomogeneity induces a Hopf bifurcation of an equilibrium, thus leading to oscillatory Tango waves. 4.2 Neural Media Although excitable media are often modelled by P D E s such as (1.1)-(1.2), in some cases it is more appropriate to use integro-differential equations (IDEs). One example is the 7 Chapter 1. Introduction following model of excitatory neural networks [14]: V t = -v-w w = e(bv — w), t + I(x) + k{x)*S{y), (1.8) ( —oo < x < oo, t>0), (1.9) where oo / k(x-y)S(v(y,t))dy (1.10) •oo is the convolution between the kernel k(x) and an output or gain function S(v). v(x,t) is a neural field that represents the local activity of a population of excitatory neurons at position x € M, the kernel k(x — y) > 0 (k(s) = k(—s)) describes the strength of the connection from neurons located at y to neurons at x and I(x) is an external input current. The neural field w(x,t) represents some form of negative feedback mechanism such as spike frequency adaptation or synaptic depression and the gain function S(v) describes how the output from each neuron depends on the state variable v for this neuron. This type of model of neural media was extensively studied by Amari [1] and has gained revived interest in recent years [14]. It models a neural medium with two interacting neuronal layers in the presence of an inhomogeneous and time-independent external forcing, I(x). The nonlinear function S(v) is usually taken to be a sigmoid function S(v) — 1/(1 + e -7(u-«;)^ w n e r K e j the threshold. In the limit 7 s —> 00, S(v) approaches a Heaviside function. In the simplest case, the kernel k(x) is taken to be the exponential function k(x) = ±e~W . d The I D E system (1.8)-(1.9) and the P D E system (1.1)-(1.2) have many common features. For example, in the homogeneous case (i.e. in the case that I(x) = 0, z(x) = 0) they are both two variable models of excitable media capable of sustaining traveling fronts and pulses [14]. Futhermore, for many of the results that we will present in further 8 Chapter 1. Introduction chapters about oscillatory waves, there are analogous results which have been discovered for (1.8)-(1.9) [4], [3], [7]. This similarity brings up an interesting question: why are the I D E system (1.8)-(1.9) and the P D E system (1.1)-(1.2) so similar in their dynamical behavior? This question has not yet been fully answered. However, it is clear that results about one of these systems can be used as guidance for work on the other. 4.3 Belousov-Zhabotinsky Reaction The Belousov-Zhabotinsky (BZ) reaction gained prominence in the 1960's as the first example of an oscillating chemical reaction. Apart from this important historical role, the B Z reaction can also serve as a prototypical excitable medium. A uniform layer of the BZ reagents sustains propagating waves such as target waves and spiral waves. In the presence of illumination, we have the photosensitive B Z reaction, which can be described by a modified Oregonator model [10]: du ~dt dv ~dl dw ~dt ^ [u - v? - w(u - q)] + D V u (1.11) u-v (1.12) 2 u ^[fv-w(u + q) + (j)(x)} + D V w 2 w (1.13) The inhomogeneous illumination <p(x) in the above is analogous to the inhomogeneity z(x) in (1.1)-(1.2). Therefore, it is natural to ask: can a stationary inhomogeneous illumination <b(x) result in oscillatory waves? The first hint at an answer has been provided in [13], where it was discovered experimentally that illuminating a small circular subregion of a two-dimensional layer of B Z reagents results in target waves periodically eminating from the illuminated area. This raises a number of other questions: how is the oscillatory target wave in the BZ reaction related to the pulse generator discovered here in numerical simulations of (1.1)-(1.2)? Can other oscillatory waves, such as Tango 9 Chapter 1. Introduction waves, be observed in the photosensitive B Z reaction? If so, are these oscillatory waves related to a Hopf bifurcation? 4.4 The Complex Ginzburg-Landau Equation The complex Ginzburg-Landau equation (CGLE) is one of the most important equations in dynamical systems. It describes extended media in which a homogeneous state is oscillatory and near a Hopf bifurcation [6]. The homogeneous C G L E reads: dA — = LLA - (1 + ia)\A\ A 2 + (1 + i(3)V A (1.14) 2 In [9], it was discovered through numerical simulations that if one modifies the above by introducing the inhomogeneity /i(r) = 1 + (u + i(5)exp [— ( r / r ) ] , then \A\ is either 2 0 a stationary pulse shaped pattern or an oscillatory breathing pattern. This leads us to ask: Why does introducing an inhomogeneity result in oscillatory waves? How are the oscillatory waves in the inhomogeneous C G L E related to oscillatory waves in inhomogeneous excitable media? Is the transition from stationary to oscillatory pattern in the inhomogeneous C G L E related to a Hopf bifurcation? 10 Chapter 2 Equilibria We begin by solving for equibria of (1.1)-(1.2). These will then be used as the starting point for the linear analysis in the next chapter. Since equilibria are defined as timeindependent states of a dynamical system, we find them by setting (v(x,t),w(x,t) (v (x),w (x)) s a = in (1.1)-(1.2). Doing this gives a single equation for v (x): s v" + f(v ) - bv + z(x) = 0, - o o < x < co s (2.1) s with w (x) = bv (x). Below, we consider two types of equilibria: stationary fronts and s s stationary pulses. Each of these is caused by a different type of inhomogeneity: fronts are caused by a monotonic z(x) and pulses are caused by a unimodal z(x). 1 Piecewise L i n e a r M o d e l For the piecewise linear model, we can construct explicit analytical forms of the equilibrium solutions. These are found by solving (2.1) with f(v) = — v + H(v — a) (see fig. 1.1b): (2.2) Notice that due to the Heaviside function, (2.1) is not defined when v = a. Therefore, s if a function v (x) statisfies v (x ) = a for some x G R, we cannot claim that it is a s s 0 0 11 Chapter 2. Equilibria solution for all x G E . Because of this, we define an equilibrium solution as a function x i—• v (x) : E —» E satisfying the following: s (El) The set D = {x G E | v (x) = a} is countable with isolated elements s (E2) v (x) G C for x G R\D 2 s (E3) v (x) satisfies (2.1) for all x G E \ D s (E4) V (XQ) S = V (XQ) = S a for all x 0 £ D (E5) ^ ( z j , ) = v'(x$) for all x G D 0 Equilibrium solutions are found by solving (2.1) on the intervals in E \ D and then matching at points x G D using (E4) and (E5). 0 1.1 Stationary Fronts We define a stationary front as an equilibrium solution satisfying the following conditions: (Fl) v (x) - z(x)/(b + 1) - » 0 as x - * oo s (F2) v {x) - (z(x) + l)/(6 + 1) - » 0 as x -> - o o s (F3) There exists a unique xo G E such that v (x) < a for x > XQ and u (x) > a for s x < s XQ. A n example of a stationary front is shown in fig. 2.1. The first two conditions above are equivalent to requiring that v"(x) — > 0 as |x| —> oo (see (2.3)). This is analogous to boundary conditions and specifies the limiting behavior of v (x), whilst accounting for s the fact that v (x) may be unbounded as x —> oo or x —> —oo. The third condition above s states that the front solution only crosses the threshold a at one location. This condition is satisfied for monotonic fronts such as the one infig.2.1. To find the stationary front solution, we solve for x < XQ and x > x and then match at x — x . Under the assumption 0 0 12 Chapter 2. Equilibria (F3), (2.2) becomes: { — z(x) if x > x 0 { (2.3) 1 —1 — z(x) if x < XQ Conditions (E4) and (E5) will give us the necessary matching conditions. For simplicity, we also impose the following restrictions on z(x): (ZO) z(x) e C 1 (Zl) l i m _ 2 ' ( x ) exists x 0 0 (Z2) l i m _ _ o o z'(x) exists x (Z3) z'[x) < 0 and l i m ^ o o z(x) < a(b + 1) - 1/2 < l i m ^ _ x 0 O z{x) Note that we allow the cases l i m ^ o o z(x) = —oo and l i m ^ - o o z(x) = oo. The conditions (Z1)-(Z3) above enable us to ensure that conditions (F1)-(F3) will be satisfied. For x < x , we solve (2.3) to get: 0 v {x) = v {x) + j^s p + Ae px b + Be'?*, x < x (2.4) Q where ft = y/b + 1, A, B are constants and V p { x ) = 2r3 J J*"M*)ds (2.5) Taking B = 0, we have that (2.4) satisfies (F2). To see this, integrate by parts: ^ - FTi WiT) (" "* L " + e e As)ds "F " +e ) e As)ds ,£ e ' (2 6) then substitute the above in (2.4) and rearrange terms: Vs{x) _ ^2+1 , _ L _ (_ -* £ , =Ae + e e A s ) d s + e -, V(s)ds ) (2.7) 13 Chapter 2. Equilibria Now let x —> —oo: lim v (x) z(x) + 1 6+1 - s I—» —CO (2.8) /"OO X — lim X—• 1 lim A e - ^ - h , 2(6 +1) V = = n 1 L / V z'(x) lim -V V ^^-^ 2(6 + 1) 2(6+1) 1 N , + 7771—TT - 0 0s —e 2(6+1) — oo 2 e z'{s)ds + e / e? J 1» , l i m + lim 1 . e -^ V x ' (2.9) / 2.10 r°e-^V(s)ds — (2.11) n x x^-oo ^(s) — lim — h p Jx x^-oo e-^ l i n w ^ *'(*)//? if P e-^V(s)ds JX •OO lim ~ ° ° + / 0x 0 if |(_~ e ^ V ( ) d S S < oo = oo J (2.12) 0 if 2 ( f e + 1 ) ] — oo - l i m ^ o o z ' ^ ) / / ? (2.13) if < oo (2.14) = 0 In going from (2.10) to (2.11) and (2.11) to (2.12), we have used l'Hospital's rule. To show that (2.14) follows from (2.13), note that: > / Mds J —oo (2.15) where M = sup (2.16) e~ z'{x) 0x z€(-oo,0] Now suppose l i m _ _ x f^° e~ z'(s)ds l3s oo 0 0 z'(x) 7^ 0. Then since z'(x) < 0 we have M < 0, which gives — oo. This establishes the implication (2.13)=>(2.14). Note also that in going from (2.9) to (2.10) and (2.12) to (2.13) we are using the fact that all the limits involved exists, which follows from (Z2). For x > x , we solve (2.3) to get: 0 v (x) = v (x) + Ce a p 14 l 3 x , x > x 0 (2.17) Chapter 2. Equilibria where C is a constant. Using the same reasoning as outlined in (2.8)-(2.14), one can show that the above satisfies ( F l ) . Putting (2.4) and (2.17) together gives: v (x) = v (x)+{ s p (2.18) 1 + 6 y Ce~ , 0x x>x 0 Using the matching conditions (E4) and (E5) gives A and C: p-0XQ A "2(6+1) «' ' =2(6TT} <' ' = C 2 19 2 20 and the front solution is: 1 ^ ) ^ = )+ f 2 - 0(*-*«>) e 276TTT s< 2{o+l) ^ -/?O-* ) < x x . e 0 ( 2 5 - 2 1 ) >xo x The value of x is determined by setting v(xo) = a: 0 -^—^+v (x ) p 0 (2.22) =a It follows from (Z3) that the above has a unique solution for x and that furthermore, (F3) 0 is satisfied. Since z(x) € C , it follows that v (x) G C so that (E2) is satisfied. Therefore, 1 2 p (2.21) is an equilibrium solution satisfying (F1)-(F3). For the case that z(x) = —ax, the stationary front is plotted infig.2.1. 1.2 Stationary Pulses We define an even stationary pulse as an equilibrium solution v (x) satisfying the following s conditions: (PI) v (x) is even s (P2) limia-i^oo v (x) s =0 (P3) Either (i) v (x) < a for all x G R or (ii) there exists x G R such that v (x) > a for s 0 15 s Chapter 2. Equilibria Figure 2.1: Stationary front in the piecewise linear model. Parameters used: a = 0.25, 6 = 1 and a = 0.1 x G (—xo, XQ) and v (x) < a for x G (—oo, — XQ) U ( X , OO) s 0 A n example of a stationary pulse is shown in fig. 2.2D. Condition (PI) is for simplicity. (P2) are the boundary conditions and (P3) ensures that v (x) — a does not have more s than one solution for x > 0. As before, we have the matching conditions (E4) and (E5). We also impose the following restrictions on z(x): (Z0) z(x) G C 1 (Zl) z(x) is even (Z2) lim^i^oo z(x) = 0 (Z3) z'(x) < 0 for x > 0 The above are sufficient to ensure that there exists a stationary pulse statisfying ( P l ) (P3). Pulses satisfying case (i) of (P3) are called subthreshold and those satisfying case (ii) are called superthreshold pulses. In order to solve for pulse solutions of any amplitude, we must treat these cases separately. 16 Chapter 2. Equilibria Subthreshold Pulses Since v (x) < a for all (2.1) becomes the following linear equation: a 0 = - ( 1 + b)v - +z(x) + v' s The solution of the above is v (x) s t —oo < X < (2.23) oo — v (x), where v (x) is as defined in the previous p p section. Notice that conditions (PI) and (P2) are satisfied and follow from (Zl) and (Z2). Superthreshold Pulses Due to (PI), we can find superthreshold pulses by solving (2.1) on the half space [0, oo), provided that we require solutions to statisfy ^(0) = 0. Solving over the intervals [0,x ) 0 and (x , oo) and matching at x using (E4)-(E5) gives the solution: 0 0 0< x < x 0 (2.24) where /3 = yb + T. XQ is determined by requiring that i> (xo) = a: s 1 _ v (x ) = a p e -2/3xo ^r— 0 = R(x ) (2.25) 0 Notice that for fixed parameter values and a given inhomogeneity z(x), the above equation may have several solutions for x , implying that multiple superthreshold pulses can 0 coexist. If (2.25) has no solutions, then there are no superthreshold pulses. Gaussian Inhomogeneity As an example, consider the Gaussian inhomogeneity z(x) = he~ / , x2 2a2 where h > 0 is the amplitude and cr > 0 is the width. In this case, the subthreshold pulses are given by: v/2^rW^ / 2 v (x) p ' 4^ 2 [ ' 17 Chapter 2. Equilibria 0 1 x 2 3 -4 -3 -2 -1 0 x Q 1 2 3 4 Figure 2.2: (A) v (x ) (solid) and R(x ) (dashed) which are repesctively the Ihs and rhs of (2.25). v (xo) is plotted for three values of h: 0.5 (lower solid), 0.9 (middle solid), and 1.3 (upper solid). R(x ) is plotted for two values of 6: 0.6 (lower, thin dashed) and 3 (upper, thick dashed). Each symbol on a solid curve represents a pulse solution for that value of h and for 6 = 3. (B) The quadratic bifurcation scenario for b = 0.6. Each solution is represented by its peak value, v (0). (C) The cubic bifurcation scenario for 6 = 3. Each symbol corresponds to the crossing point denoted by the same symbol in (A). (D) Profiles of the three coexisting pulse solutions for h = 0.9 and 6 = 3. In (A) and (D), the dash-dotted line represents the threshold a. Parameter values: a = 0.25, a = l. p 0 Q p 0 s The superthreshold pulses are found by solving (2.25) for x , with v (x) as given above, 0 p and then using these values of XQ in (2.24). The existence of stationary pulses is summarized by plotting the relation between v (0), the maximum of the pulse, and h, the s strength of the inhomogeneity. This gives diagrams as shown in figs. 2.2B and 2.2C. In fig. 2.2C, parameter values satisfy a — l/(2/3 ) > 0 and we have what we call the 2 cubic scenario. In fig. 2.2B, a — l/(2/? ) < 0 and we have the quadratic scenario. The 2 18 Chapter 2. Equilibria shape of the bifurcation diagram in the cubic and quadratic scenarios can be deduced by looking at fig. 2.2A. Explanations are given in the figure captions. Notice that in the cubic scenario we have a saddle-node bifurcation at h = hi. 2 Cubic Model In the case of the cubic model, equilibria are given by: v'' + f(v,)-bv + z(x) = 0, x e l a (2.27) where f(v) = v(l — v)(v — a). Since analytical solutions of the above are very difficult to obtain, we resort to numerical methods. Consequently, we must use the finite do- main / = [—a,a], with the hope that a large enough domain can approximate the case of / = (—00,00). In both the case of a stationary front and pulse, we solve for equi- librium solutions using M A P L E . These results are then confirmed independently using asymptotic methods (for some limiting cases). 2.1 Stationary Fronts For stationary fronts, we consider the linear inhomogeneity z(x) = z — ax. c When implementing numerical and asymptotic methods for this inhomogeneity, it is useful to rescale variables in order to simplify (2.27). Let: v {x)=v a + y/p-bv{x) c (2.28) where v = (1 + a)/3 and p = 3v% — 3v + 1. In the new variable u, (2.27) now reads: c c ( 19 2 - 2 9 ) Chapter 2. Equilibria Figure 2.3: The stationary front computed using M A P L E is the thick line and the asymptotic front is the thin line. Parameters used: o — 0.02, = 0.1 and o = 0.2 where we have also taken z = bv — f(v ). Now let v(x) — V(y), where y = xa/(p — b) ^ . 3 c c c 2 dV 2 Then a 2 Our boundary conditions for the above are V(y) —> V±(y) as x —> Too, where V+(y) and V-(y) are, respectively, the upper and lower branch solutions of V — V + y = 0. Notice 3 that these are analogous to the boundary conditions (F1)-(F2) used for the piecewise linear model. Equation (2.30) is solved numerically using M A P L E (details are given in the appendix). 20 Chapter 2. Equilibria Asymptotic Solution In order to provide evidence in support of our numerical result, we can solve (2.30) in the case that a <C 1. Solving (2.30) in the outer regions and matching at the transition layer gives the following asymptotic solution: (2.31) A n example of an asymptotic solution alongside a M A P L E solution is shown in fig. 2.3. 2.2 Stationary Pulses For stationary pulse solutions, we consider the Gaussian inhomogeneity z(x) = he~ / . x2 2 Stationary pulse solutions are found numerically by using M A P L E to solve (2.27) with Dirichlet boundary conditions (details are given in the appendix). As with the piecewise linear model, we can illustrate the existence of stationary pulses by plotting v (0) as 3 a function of h. This gives the bifurcation diagrams shown in figs 2.5, which can be again separated into the cubic and quadratic scenarios. These diagrams are qualitatively similar to the piecewise linear case, the main difference being that in the cubic scenario, there is a cusp bifurcation (see fig. 2.4). As with stationary fronts, we can support our M A P L E solutions using asympototic methods. Asymptotic Method Unlike the case of the front, we do not have an explicit asymptotic formula for the pulse profile. However, the shape of the bifurcation diagrams (i.e. v (0) as a function of h) can s be determined asymptotically for a <C 1. To do this, multiply (2.1) by v' and integrate s 21 Chapter 2. Equilibria h (a) (b) Figure 2.4: a) The "S-shaped" curve has saddle-node bifurcation points at both knees. As the parameter b is increased, these knees coalesce until we are left with a monotonic curve. Parameters: a = 0.35, b = 0.146 (lower curve), b = 0.126 (middle curve) and b = 0.106 (upper curve), b) Each point represents a "knee" of the curves such as in a). The points join together in a singularity in the h — b parameter space. This is indicative of a cusp bifurcation. Parameter: a = 0.35. from x — —oo to 0 to get: / (f(v) - bv)dv = -h JO e- v' {x)dx (2.32) x2/2a2 s J-oo where v$ = v (0) and we have used the fact that v' (—oo) = ^(0) = 0. For a <C 1, s s the integral on the right can be evaluated asymptotically (assuming that v (x) remains s sufficiently regular as a —> 0): f e- v' {x)dx~ x2/2cT2 s 22 -a v' \0) 2 s (2.33) Chapter 2. Equilibria So for cr <C 1, we have rvo / (f(v) - bv)dv ~ /ICT <(0) (2.34) 2 But we know that v"(x) + f(v ) — bv + z(x) = 0, so we have v"(0) + /(fo) — bvo + h = 0. s s Using this in the above gives: / {f(v) - bv)dv - -h(f(v ) Jo 0 - bv + h)o 2 0 (2.35) Treating the above as a equation, we can compare this asymptotic formula with our M A P L E results. This is shown in figure 2.5. 23 Chapter 2. Equilibria Figure 2.5: Relation between the amplitude of the pulse, v (0), and the strength of the inhomogeneity, h, for pulses in the cubic model, a) is the quadratic bifurcation scenario and b) is the cubic bifurcation scenario. Parameter values: a = 0.35 and b = 0.036 (for a)) or b — 0.086 (for b)). The lines are the asymptotic results and the crosses are points obtained using M A P L E . a 24 Chapter 2. Equilibria X X (a) (b) Figure 2.6: a) Stationary pulse found using M A P L E b) Stationary pulse found in numerical simulation. When placed on the same figure, the curves in a) and b) are virtually indistinguishable. Parameter values: a — 0.35, b — 0.3 and h — 0.384. e = 0.115 for the simulation in b). 25 Chapter 3 Linear Analysis Most nonlinear dynamical systems are described by equations which are, for the most part, analytically intractable. However, it is possible in many cases to describe the dynamics analytically in the vicinity of an equilibrium solution. The starting point for such an analysis is a linearization of the equations about the equilibrium. This then gives a linear equation which, in most cases, can be readily analyzed by separating solutions into a superposition of eigenmodes. The eigenvalues are the growth rates of these eigenmodes, and these can be used to determine the stability and possible bifurcations of the equilibrium. In finite dimensional dynamical systems, the Hartman-Grobman theorem (along with a theorem stating that two linear systems are topologically equivalent only if the real parts of their eigenvalues are of the same sign) ensures that bifurcations of the equilibrium are only possible when one or more critical modes are present. Furthermore, if the critical modes (i.e. corresponding to eigenvalues on the imaginary axis) can be perturbed under a parameter variation, then one has a bifurcation. The details of the bifurcation can then be analyzed using the center manifold theorem, and the topological dynamics are determined using the reduction principle. Combining these previous tools allows one to determine (modulo topological equivalence) all local dynamics in a neighborhood of the equilibrium. In many cases, bifurcation branches of equilibria can be extended beyond the bifurcation point and basins of attraction of stable equilibria can 26 Chapter 3. Linear Analysis be extended beyond the neighborhood of the equilibrium. In such cases, local dynamics (such as local bifurcations and stability) can be of help in understanding the global dynamics. For infinite dimensional systems, it is usually much more difficult to rigorously prove stability and bifurcation of equilibria. Many technical requirements must be met in order to apply the principle of linearized stability and the center manifold theorem. Since a consideration of these technical requirements is outside the scope of this thesis, our approach is use the modes of the linearization to look for possible bifurcation points, and then compare with numerical simulations. The motivation for this is based on the idea that either our infinite dimensional system can be approximated by a finite dimensional system, or that the technical requirements which have been ignored could be shown to be satisfied (with sufficient work). The result of our linear analysis is the discovery of a pair of critical modes with complex conjugate eigenvalues, suggesting a Hopf bifurcation, and thus the emergence of oscillatory solutions. This Hopf bifurcation occurs with a change in a model parameter, most notably the reduction of the strength of the inhomogeneity. Let (v (x),w (x)) s s be an equilibrium solution of the Fitzhugh-Nagumo equations. Let us introduce the new variables v(x,t) = v(x,t) — v (x) and w(x,t) = w(x,t) — w (x). In s s these new variables, the F - N equations read: dv 2 f(v + v)- s dw = bv ~dt — f(v ) a e + dx w 2 (3.1) (3.2) Linearizing the above equations about (v, w) dv -w (0, 0) gives: dv dx 2 vf'(v ) s dw = bv — w ~dt 27 w + 2 (3.3) (3.4) Chapter 3. Linear Analysis In order to find the eigenmodes of the above equations, we separate variables by writing (v,w) = e $(x), where Xt = (cb(x),ib(x)). This gives the following eigenvalue/eigenfunction problem for A and $: £ $ = A$ (3.5) where f'(v )+d /dx 2 £ =5 -1 2 s b/6 (3.6) -1/5 with 8 = 1/e. Simplifying (3.5) gives: 4>(x) ${x) = (3.7) 6/(1 + A) with cj)(x) given by the following Schrodinger eigenvalue problem: <t>" + f'(v )<j> = <i>, xei a (3.8) x For piecewise linear f(v), we have / = [—a, a] with the boundary conditions 0(±a) = 0. For cubic f(v), we have / = (-co, oo) with the requirement that 0(±oo) be bounded. A is given as the two (possibly complex) roots of: eA + (3.9) X A+l Let A and A be the solutions of the above. Then provided that A ^ A , every eigen+ - + - value/eigenfunction pair of the Schrodinger problem yields two eigenvalue/eigenfunction pairs of the original eigenvalue problem (3.5). A and A are given by: + _ ix ~ 0 l± A* = - ± V(X - e) - 4e(6 2e 2 X ) (3.10) Using the above with the spectrum of the Schrodinger eigenvalue problem will allow us to identify critical modes. For the problems considered in this thesis, the eigenvalues 28 Chapter 3. Linear Analysis of the Schrodinger problem will be real and nondegenerate. This limits the eigenvalue spectrum so that critical modes can only give rise to three possible local bifurcations, as shown in table 3.1. Of the three bifurcations listed, we will be mainly interested in the Hopf bifurcation. The table can be summarized by saying that a Hopf bifurcation occurs if e < 6 and x = - If e Xmax is the largest eigenvalue in the spectrum, then an equilibrium destabilizes through a Hopf bifurcation at Xmax = e < 6. Table 3.1: Possible local bifurcations of equilibria. The condition for the occurrence refers specifically to the value of the ratio 6/e. Values of A Type of bifurcation Condition for the occurrence Value of x Saddle-node (SN) 6/e = any value X = b A" = 0, A+ = 6/e - 1 Hopf (HP) 6/e > 1 X=« A± = iy/b/e - 1 Double-zero (DZ) 6/e = 1 X= b=e A± = 0 1 1.1 ± Piecewise L i n e a r M o d e l Eigenvalues for Stationary Fronts If /(y) = —v + H(v - a), then f'(v ) s = - 1 + 8(v„ - a), where 8(-) is the Dirac delta function. Using a property of the delta function and the fact that v (x) = a iff x — x , s 0 we have 8(v -a)= s 8(x - s ) IW o ) | 1 0 . (3.H) so that the eigenvalue problem (3.8) becomes: 4' + . * 8(x-x )(f)= (x + l)0, 0 29 -oo<x<oo (3.12) Chapter 3. Linear Analysis This is a Schrodinger problem with a 5 function "potential well" located at x = XQ. If an eigenvalue of (3.12) is isolated, it is said to be part of the discrete spectrum. Otherwise, the eigenvalue is part of the continuous spectrum. Thus the eigenmodes can be classified as being part of either the continuous spectrum or discrete spectrum. It is a property of Schrodinger problems that eigenvalues in the continuous spectrum satisfy x < ~ 1> - d an those in the discrete spectrum satisfy x > ~ 1- Since eigenmodes with x < ~ 1 cannot be critical (see table 3.1), we can focus on modes with x > — 1 (there is in fact only one such mode). This mode must satisfy the following conditions: (i) lim.ia.i_oo (f)(x) = 0 (ii) <j>(x) is continuous (iii) 4>(x) is smooth everywhere except at x = XQ Condition (i) is a property of eigenfunctions corresponding to discrete eigenvalues for Schrodinger problems. Condition (iii) reflects the fact that because of the delta function, we expect a jump discontinuity in <j)'{x) at x = XQ. In order to find the magnitude of this jump, we integrate (3.12) from XQ — e to XQ + e, where e > 0. This gives I <t>"dx+ 8(x - x )(pdx = ( x + 1) / <t>dx 0 (3.13) so that (j)'(x + e)-(j)'(xo-e)+ 0(s ) = ( x + 1) / 0 <pdx o IW o ) | (3.14) Jxo-e Letting e —> 0 gives the jump condition: -* -iSi W) - w= (3 i5) The eigenfunction is now found by solving (3.12) on x < x and x > XQ and then matching 0 at x = XQ using the above and condition (ii). Solving (3.12) for x < x and x > x gives: 0 Ae rx <Kx)={ + Be~ , x < x 0 Ce rx 0 rx + De~ ,x rx 30 > x 0 (3-16) Chapter 3. Linear Analysis where r = \Jx + 1 and A, B, C and D are constants. Condition (i) above implies that B = C = 0. Using the continuity condition 4>(XQ) = <J>(XQ) then gives Ae ° rx = De~ , rXo so that: f *- , x < x <p(x) = K{ [ e- - ,x > x r ( s o ) e 0 r(x (3.17) Xo) 0 where K ^ 0 is an arbitrary constant. Finally, using the jump condition (3.15) gives the eigenvalue x : X = - l + ,) M ( 3 2 1 8 ) Since the above is the only eigenvalue satisfying x > — 1, it is the largest eigenvalue. Hopf Bifurcation of a Front Suppose we are considering stationary fronts for the inhomogeneity z(x) = —ax. Then (2.21) gives |fs(z )| = 0 (/2/2 + o~)/(b + 1) and the eigenvalue is: =- {jTra) X <' > 1+ 3 19 To emphasize that x depends on a, we write x = xi )- We have a Hopf bifurcation at a a critical slope a = a > 0 if: c (1) e < 6 (see table 3.1) (2) x(cr ) = e (see table 3.1) c (3) ^|rj=a 7^ 0 (eigenvalue crossing condition) c If we set: then for e < b conditions (l)-(3) above are satisfied, with a stable front for a > a and c an unstable front for a < a . The stability and bifurcation of a stationary front are c summarized in fig. 3.1a, where we have plotted the Hopf bifurcation line in the e — a parameter space. 31 Chapter 3. Linear Analysis Figure 3.1: Hopf bifurcation line in the e — a parameter space for a) the piecewise linear model and b) the cubic model. Parameter values: a) b = 1. b) a = 0.2, b = 0.1. 1.2 Eigenvalues for Stationary Pulses Subthreshold Stationary Pulses For subthreshold stationary pulses, we have v (x) < a and the eigenvalue problem (3.8) s becomes: 4>" = (X + 1)0 (3-21) with the condition that 0(±oo) be bounded. Solving the above gives a continuous spectrum x < ~ 1- Since none of these eigenvalues correspond to critical modes, there are no local bifurcations of subthreshold stationary pulses. Furthermore, subthreshold pulses are always stable. 32 Chapter 3. Linear Analysis Superthreshold Stationary Pulses Each superthreshold pulse solution v (x) crosses the threshold a twice at x = ±XQ. For s these pulses, the S function in the Schrodinger problem (3.8) can be expressed in the following form: 5{v {x) - a) = s rr + — (3.22) |u'(-a;o)| |u'(a;o)| This corresponds to a "double delta function" potential for the Schrodinger eigenvalue problem (3.8). For this eigenvalue problem, we have a continuous spectrum x < ~ 1 a n d a set of discrete eigenvalues satisfying x > ~ 1 (this can be deduced either by solving the eigenvalue problem explicitly or as a general property of Schrodinger eigenvalue problems with appropriate potentials). Since there are no critical modes in the continuous spectrum, we focus on the discrete spectrum. The eigenvalues of the discrete spectrum are found by solving for the eigenfunctions as in section 1.1 above (using conditions (i)—(iii) as before). Doing this gives the following equations: 1 ± -2Vx+ix e 0 2K(zo)| Each equation above has a single solution. By picking the '+' ('-') sign in (3.23), we obtain the eigenvalue x+ (X-) t n a u corresponds to the even (odd) eigenfunction. By a property of Schrodinger eigenvalue problems, we know that the leading eigenvalue must correspond to the even eigenfunction. Hopf Bifurcation of Stationary Pulses In the previous chapter, we constructed explicit stationary pulse solutions for the Gaussian inhomogeneity z(x) = he~ . x2/2 The existence of such solutions was summarized using diagrams such as figs. 2.2. In this section, we show that in the cubic scenario, the upper branch of solutions on these diagrams can undergo a Hopf bifurcation as the 33 Chapter 3. Linear Analysis Stable SN _J 0 2 1 3 I 4 5 e B 6 HP 4 Stable 2 0 SN I 0 L_ 1 2 3 4 5 6 7 e Figure 3.2: Stability diagram for stationary pulses in the two-parameter e — h space. In both panels, H P denotes the curve of Hopf bifurcation points and SN represents the line of saddle-node points. The two intersect at a codimension two double-zero bifurcation marked by a filled circle. In (A), the H P happens at a finite value of h, irrespective of the value of e. Parameter values: a = 0.35, b — 4. Thus, x+(oo) « -0.6 < 0. In (B), the H P cannot happen at finite values of h for small values of e. Parameter values: a = 0.17, 6 = 6. Thus, x+(oo) « 0.24 > 0. parameter h is varied. This gives the bifurcation diagram as shown in fig. 3.3B. In order to show that a Hopf bifurcation occurs, we must solve (2.25) for x > then use this in 0 (3.23) to solve for x- To illustrate the dependence of x+ X+{h)- If X+(h ) = e < b and x+(M c o n the parameter h, we write 0, then we have a Hopf bifurcation at h = h c (see table 3.1). For the cubic scenario (fig. 2.2C), we can calculate: (3.24) X±(oo) = If e < b and Aa (b+ 1) 2 (3.25) • — 1 < e, then by continuity of x+W we must have a Hopf bifurcation on the upper branch of the bifurcation diagram. Stability and bifurcation diagrams in 34 Chapter 3. Linear Analysis 0.6 0.4 o >" 0.2 0 0.5 0.8 _ o 0.6 >" 0.4 0.2 1 0.4 2.5 1.5 3.5 j- C 0.3 S > 0.2 0.1 0 Figure 3.3: Bifurcation of the stationary pulse solutions shown by v (0) as a function of h. The filled circles stand for stable solutions and empty circles for unstable ones. In all cases, the middle branch is unstable and the lower branch is stable. In (A), the upper branch is also stable. Thus, bistability occurs for all h in (hi, h ). This happens for a = 0.2, 6 = 2, e > b = 2. In (B), a Hopf bifurcation point divides the upper branch into an unstable and a stable segment. The parameter values are: a = 0.35, b = 4, e = 0.2. In (C), the upper branch is unstable for all h values. The parameter values are: a = 0.17, 6 = 6, e — 0.2. s r the e-cr parameter space are shown in fig. 3.2. Note that the exact location of the Hopf points was determined by solving (3.23) numerically. In later sections, the location of the Hopf points will be compared with the results from numerical simulations. 2 Cubic Model In the cubic model, we do not have either explicit analytical expressions for stationary solutions, nor do we have analytical methods to compute the eigenvalues of (3.8). Thus, 35 Chapter 3. Linear Analysis we use numerically determined stationary solutions and solve the eigenvalue problem (3.8) numerically using M A P L E (details are given in the appendix). Note that because we will be working on the finite domain x G [—a, a] with Dirichlet boundary conditions, the eigenvalue problem (3.8) is a Sturm-Liouville problem. 2.1 Eigenvalues for Stationary Fronts For the linear inhomogeneity z(x) = —ax, we computed stationary front solutions in the previous chapter. For these stationary solutions, the eigenvalues of the linearization are found by solving the Sturm-Liouville problem (3.8) with v (x) s as shown in fig. 2.3. If the largest eigenvalue is critical, then we have a bifurcation where the stationary front changes from stable to unstable. We thus focus on the largest eigenvalue. It is known that the eigenfunction with the least number of zeros is the one corresponding to the largest eigenvalue (this follows from Sturm's comparison theorem). In our case, this leading eigenfunction has no zeros (as determined numerically). This fact is used to compute the leading value (see appendix for more details). Computing the resulting eigenvalues for several parameter values and joining these points together gives the Hopf bifurcation line in e — h parameter space, as shown in figure 3.1b. This diagram is analogous to the one for the Hopf bifurcation of fronts in the piecewise-linear model. The numerically determined Hopf points will be compared with numerical simulations in later sections. 2.2 Eigenvalues for Stationary Pulses In order to find the eigenvalues of the linearization about stationary pulses (for the Gaussian inhomogeneity z(x) = he~ / ), x2 2(j2 with v (x) s we must solve the Sturm-Liouville problem as shown in figure 2.6. Following the same strategy as for the stationary front, we look for an even eigenfunction without zeros and compute the corresponding eigenvalue. We then use these eigenvalues, along with the corresponding values of h, to 36 Chapter 3. Linear Analysis form a set of Hopf points in the e — h parameter space. Connecting these points in a curve gives the diagrams shown in fig. 3.4. Bifurcation diagrams are shown infig.3.5. Notice that in contrast to the piecewise linear model, in the cubic scenario we have Hopf bifurcations at two values of h and in the quadratic scenario, we have a Hopf bifurcation on the lower branch of solutions. In later sections, we will supplement these diagrams with numerically determined periodic solutions, giving a fuller picture of the bifurcations. 0.4 b=0.166 Stable 0.3 B 0.4 |t>=0.126 Stable 0.3 h Unstable " ^& 0.2 . O - C H X ^ " HP >:ttr 0.1 0.06 0.12 0.2 0.18 U nstable^ 0.1 HP ;-^v ^x,;^ 0.05 0.1 0.15 e e 0.28 b=0.086 0.34 b=0.106 r 4 0.26 0.18 0.18 0.15 0.08 Figure 3.4: Bifurcation diagrams in e — h space for four different values of b. The location of the Hopf points is denoted by HP. SN denotes the line of Saddle-node bifurcation points. The saddle-node points occur at the knees of bifurcation diagrams such as the one in 2.5b. Other parameters: a = 0.35. 37 Chapter 3. Linear Analysis Figure 3.5: Qualitative diagrams showing stability/bifurcation changes for the quadratic and cubic bifurcation scenarios, a) In the quadratic scenario, the upper branch is unstable and the lower branch is divided into stable and unstable by a Hopf bifurcation point denoted HP. b) In the cubic scenario, we have Hopf bifurcations on both the upper and lower branches. 38 Chapter 4 Numerical Simulations We now come to the main theme of this thesis, which is to relate oscillatory waves discovered in numerical simulations with the Hopf bifurcation discovered in the previous chapter. The two oscillatory waves considered are the oscillatory front and the oscillatory pulse, each of which can be traced back to the Hopf bifurcation of an equilibrium. The oscillatory front or "Tango wave", which in the simplest case is caused by a linear inhomogeneity z(x) = —ax, can be understood as emerging from a supercritical Hopf bifurcation of a stationary front, with the bifurcation parameter being the slope of the inhomogeneity a. For large values of the slope, the stationary front is stable and "pinned" at the origin. As the slope is decreased, the front destabilizes and is released into a front moving in a back and forth manner, called an oscillatory front or Tango wave. The oscillatory pulse (or breathing pattern), which in the simplest case is caused by a Gaussian inhomogeneity, also emerges from a Hopf bifurcation as the strength of the inhomogeneity is decreased. 1 Oscillatory Fronts A n oscillatory front can be described as a wave which is a traveling front (moving to the right, see fig. 1.2a) for half of its period and a traveling back (moving to the left, see fig. 1.2b) for the other half. During the transition from a front to a back, the wave slows 39 Chapter 4. Numerical Simulations down, stops and then reverses direction. Figure 4.7a shows an example of v(x,t) for an oscillatory front caused by the inhomogeneity z(x) = —ox. It is observed numerically that the amplitude of oscillations decreases as the strength o of the inhomogeneity is increased, and that for sufficiently large inhomogeneities, initial conditions approach a stationary front. This suggests that a Hopf bifurcation can be found at a critical value of the slope o. This is confirmed by plotting the amplitude of oscillatory fronts for several values of the parameter o, as in fig. 4.1, and comparing this to the Hopf point determined in the previous chapter. For all cases tested numerically, we found that the Hopf bifurcation was supercritical. 2 Oscillatory Pulses For oscillatory pulses, the bifurcation behavior is slightly different in the piecewise linear and cubic models, and thus we treat them separately. 2.1 Piecewise Linear Model An oscillatory pulse can be described as a bump shaped pattern whose amplitude oscillates in time. A n example of an oscillatory pulse is shown in figure 4.7b. In this and other examples in this section, we are using the inhomogeneity z(x) = he~ ^ . The x2 2 amplitude of the oscillations can be quantified by finding the maximum and mininum v(0,t). Doing this for several values of the parameter h gives the bifurcation diagrams as shown in figure 4.2. In both of these diagrams, as we increase the parameter h the oscillatory pulses emerge through a bifurcation and then disappear through another. In fig. 4.2a, the gradual decline in amplitude of oscillations, along with the termination point coinciding with the Hopf point, points to a supercritical Hopf bifurcation. Infig.4.2b, on the other hand, the oscillatory solutions continues noticably beyond the Hopf point, 40 Chapter 4. Numerical Simulations 3 2 1 X o -1 -2 a -3 0.008 0.01 0.012 0.014 Figure 4.1: Amplitude of front oscillations (open circles) found in numerical simulations compared with the Hopf bifurcation found in the previous chapter (cubic model). The thin line represents the unstable stationary solution and the thick line the stable one. The boundary between the thin and thick horizontal line segments is the Hopf point. The location x(t) of the moving front is found by setting v(x(t),t) = v , assuming this equation has only one solution. In the above, a is the bifurcation parameter. Parameter values: a = 0.2, b = 0.1 and e = 0.02. c suggesting a subcritical Hopf bifurcation. The bifurcations occuring at the lower values of h do not appear to be related to any local bifurcations. The most likely explanation is that these are global bifurcations, possibly homoclinic, and related to the codimension two double-zero bifurcation (see fig. 3.2). In our simulations, we discovered that for e <C 1, the dynamical behavior is no longer simply a pulse oscillation, but instead a pulse shaped pattern periodically emitting traveling pulses (see fig. 4.3). We call this spatio-temporal structure a pulse generator. 41 Chapter 4. Numerical Simulations v 1 • • * 0.8 1.2 1.6 rB 1.5 2 h 2.5 Figure 4.2: Bifurcation diagrams for pulses in the piecewise linear model. For each oscillatory pulse, both the maximum and the minimum values of v(0, t) are plotted (filled dots). Unstable and stable stationary pulses are in dotted and solid curves, respectively. Only the peak value v (0) is plotted for each stationary pulse. The H P appears to be super-critical in (A) at hn = 1.89096 and sub-critical in (B) at hn — 2.34.. Parameter values: a = 0.17, b = 6 and e = 1.75 in (A); a = 0.35, b = 4 and e = 0.5 in (B). s 2.2 Cubic Model Recall from the previous chapter that in the cubic model, there were (in some cases) two Hopf bifurcations as the parameter h was varied. This gave bifurcation diagrams such as fig. 3.5. We have numerically computed the amplitude of pulse oscillations for various values of h and have plotted these on the bifurcation diagram, with the hope of connnecting the oscillatory behavior with the two Hopf bifurcations. In figures 4.4A and 4.4B this connection is immediately apparent. In addition to these supercritical Hopf points, in figures 4.4B and 4.4C, we have what appears to be subcritical Hopf bifurcations, although this can only be inferred indirectly. Provided that subcritical Hopf bifurcations are indeed present, then combining figures 4.4A and 4.4B would allow us to conclude that 42 Chapter 4. Numerical Simulations l=D.O/ i 1 1 1 1 -50 -25 0 25 50 Figure 4.3: A n oscillatory pulse which emits traveling pulses, called a pulse generator. This generator is pinned at the origin by the inhomogeneity. The spatio-temporal profile of this solution is shown in snap-shots at six time moments during one oscillation cycle. The inhomogeneity z(x) is plotted as a dotted curve in each panel. Parameter values: a = 0.35, b = 4, e = 0.01, h = 2.5, a = 1. a transition from a supercritical to subcritical Hopf bifurcation is present, thus signalling the presence of a generalized Hopf bifurcation. The connection between the oscillatory waves and the generalized Hopf bifurcation will be explored in more detail in the next chapter. As in the piecewise linear model, if 6 < 1 we find that the oscillatory pulse emits a pair of traveling pulses periodically. A n interesting aspect of the pulse generator is that the frequency with which it emits pulses varies with the strength of the inhomogeneity. 43 Chap 0 0.2 h 0.4 0.6 0 0.2 h 0.4 0.6 Figure 4.4: Bifurcation diagrams for pulses in the cubic model. Periodic pulses obtained in numerical simulations are marked by filled dots. (A) The periodic branch emerges and disappears through a supercritical Hopf bifurcation at both ends of the domain for e = 0.115, called a sup-sup combination. (B) A sup-sub combination in which the periodic branch starts at a supercritical Hopf and ends at a subcritical Hopf for e = 0.02. In (A), (B), and (D) a = 0.25, b = 0.3. (C) A sub-sub combination. Parameter values: a = 0.2, b = 0.17 and e = 0.08. (D) Bifurcation diagram for a pulse generator for e = 0.01. This is illustrated in fig. 4.5. Fig. 4.6 also shows how the period of oscillations appears to diverge as the strength of the inhomogeneity approaches a critical value. The exact nature of this bifurcation is unknown. The mode-locking behavior seen in fig. 4.5 is indicative of a torus bifurcation (this was suggested already in [7]). On the other hand, the diverging period is characteristic of a homoclinic bifurcation. Another interesting phenomenon is the spatio-temporal pattern shown in fig. 4.8b. In this case, pulses are emitted every half-period, with emissions in one direction at one half of the cycle and then in the other direction in the other half. This type of pattern could also be qualified as a pulse generator, although the pattern emitting traveling pulses 44 Chapter 4. Numerical Simulations i i i i i i i i i i 0 10 20 Time, t 30 40 0 10 20 Time, t 30 40 Figure 4.5: Complex oscillations of v(0, t) for a pulse generator. As we can see, the period of oscillations increases with h. With each large excursion from the stationary solution near v = 1, the oscillatory solution emits a pair of traveling pulses. Parameter values: a = 0.25, 6 = 0.3, e = 0.01, a = 1. is no longer a breathing pattern. Instead, it resembles more what is sometimes called a "wiggling" pattern, characterized by side-to-side movement of a pulse. This type of pattern could be related to a Hopf bifurcation of an unstable pulse. 45 Chapter 4. Numerical Simulations 50 r 40 •o30| o — i Q -20|o 101 3.4 <*» I I I l 0.42 0.44 Parameter, h 0.46 0.48 Figure 4.6: Period of oscillations as a function of inhomogeneity strength h. The period diverges as h approaches 0.47. Parameters used: a = 0.25, b = 0.3, e = 0.01. Figure 4.7: a) Fully developed oscillatory front (Tango wave). Parameter values: a = 0.2, b = 0.1, e = 0.02 and a = 0.007. b) Breathing pattern. Parameters: a = 0.25, b - 0.3, e = 0.02 and h = 0.3. 46 Chapter 4. Numerical Simulations Figure 4.8: a) Pulse generator, b) Pulse generator where pulse emissions are out of phase with one another. Parameter values for both a) and b): a = 0.25, b = 0.3, e = 0.01 and h = 0.28. The pattern in a) is obtained using the initial condition v(x,0) = w(x,0) = 0 and the pattern in b) is obtained using the initial condition v(x, 0) = H(—x), w(x, 0) = 0. 47 Chapter 5 Bifurcation Analysis In the previous chapter, figures 4.4A and 4.4B suggest that oscillatory pulses can be linked not only to a codimension one supercritical Hopf bifurcation but also to a codimension two generalized Hopf bifurcation. On the other hand, oscillatory fronts seem to be only linked to a supercritical Hopf bifurcation. To reinforce these claims, in this chapter we determine the reduced dynamics on the centre manifold at the Hopf points, which are characterized by a single number, called the first Liapunov coefficient, denoted by h(0). If Zi(0) > 0 (Zi(0) < 0), then the Hopf bifurcation is subcritical (supercritical). If / (0) = 0, then we are at a generalized Hopf point. In what follows, we compute ^(0) for x the Hopf bifurcation of fronts and pulses in the cubic model. It should be noted that, as in previous chapters, the justification for our approach depends upon comparisons with numerical simulations. Our approach is to approximate the first Liapunov coefficient and ignore most (if not all) technical requirements normally needed to apply center manifold reduction. First, it is useful to introduce new variables v(x,i) = v(x,t) — v (x) a and w(x,t) = w(x,t) — bv (x). In these new variables, (1.1)-(1.2) becomes s = CU + F{U), xe[-a,o] (5.1) at where U = (v,w) T and we have Dirichlet boundary conditions at both ends. The linear 48 Chapter 5. Bifurcation Analysis operator £ is: 2(a + l)v (x) - 3v (x) - a + D dx 2 a £ s 2 (5.2) b/8 -1/8 and the nonlinear terms are: (a + 1 - 3v (x))v 2 F(U) = 8 s - v 3 (5.3) 0 In the above, 8 = 1/e. Recall from chapter 3 that the spectrum of the linearization is given by the eigenvalue problem £ $ = A<J>, which can reduced to the following SturmLiouville eigenvalue problem: 4>"(x) + {2{a + l)v (x) - 3v (x) 2 8 s - a) cj>{x) = <l>(x), *(x) = C b/(X + 1) <j>(x) x (5.4) with x — ^ + b/(X + 1) and C is a (possible complex) normalization constant. If e < b, e then at x = e we have a pair of pure imaginary eigenvalues X — iui, X = — iu, with u — \Jbjt — 1, and all other eigenvalues with 3ft(A) < 0. Throughout what follows, we assume that parameters are fixed so that X = 6 a n d we are at the Hopf bifurcation point. Let {$ >^o} be the eigenfunctions with pure imaginary eigenvalues (that is, the 0 eigenfunctions of the leading eigenvalues). Also, let E° = span{$ , ^o} and let E be the s 0 space spanned by the eigenfunctions with 5ft(A) < 0. B y properties of Sturm-Liouville eigenvalue problems, our eigenfunctions { $ „ , $ „ } form a complete basis for some Hilbert space H (say, L ), which we take as our solution space. Thus we have TC = E © E and 2 c s we can write: (5.5) 49 Chapter 5. Bifurcation Analysis where V(x, t) € E ,z = ($ , U) and z = (l> , [/). Here $ (a;) is the adjoint eigenfunction, s f 0 0 which is the solution to the adjoint eigenvalue problem C & = \& and is given by: T &{x) = 1 D 0(z) (5.6) -5/CX + l) where D is a complex normalization constant. The eigenfunction $ 0 and the adjoint eigenfunction $ J are normalized so that they satisfy the relative normalization ( $ o i $0) 1. The inner product here for any two vector valued complex functions v — (v\, v ) T 2 u = is (ui,U2) T = (v,u) J^ (viUi + V2U2) a = and dx. The normalization of $ 0 is arbitrary, as long as the relative normalization is kept the same. We choose C = 1, which gives 5 = (l-i/w)/(2(0 0». l z(t), z(t) and ip(x,t) can be thought of as new coordinates in our space H. To get the time evolution of z(t), z(t) and tp(x,t), we project onto (5.1) to get: z = iuz + ( $ , F ( z $ + z$o + *)) 0 (5.7) 0 ^ = £47 + F(z% + z$ Q + *) - ($ , F(z$ 0 + z$o + * ) ) $ o - ($o> ( ^ o + z$o + * ) ) * o F 0 (5.8) The center manifold theorem states that in the phase space H there is a local (i.e. near U = 0) manifold, Wf , which is tangent to E . Furthermore, this manifold is invariant c oc and locally exponentially attracting under the flow of (5.1) [5]. The tangency condition allows us to write the center manifold V(z, z) : E° —> E as: s V(z,z) = \w z 2 Q2 + W zz+ -W z ^0{\z\ ) with W20 = W 2- Since W £ is invariant, Q c l n 2 z 2Q (5.9) = V(z,z) must be a solution to (5.7)-(5.8). Substituting ^ = V(z, z) into (5.8) and collecting terms at 50 0(|2|2) gives us the following Chapter 5. Bifurcation Analysis boundary value problems for W 2, Wu and W : 0 V + l + z6a; \ . , 3 5(1 + 4a; ) ) 2 W ° 2 ^ = 20 , , r L w ° 2 ^ + bw (x) = Lw (x) n u .. n 2 g ( >' x W ° 2 w 2{x) = 0 (5.10) 6/(1 + 2iu) 1 . + 2g(x), W = (5.11) w (x) u n a - a + Df^] and g(x) = (a + 1 - 3v )4> - (0, (a + 2 where L = [2(a + l)v (x) - 3v (x) 2 a 2 s s 1 - 3v )4> )(f)l'(0,0). Finally, substituting ^ = V(.z, £) into (5.7) gives: 2 s i 1 = zu;z + -g z + g zz + 920 = 902 = 9u = 2D<5(0, (a + 1 - 3v )q> ) (5.13) gi = 2D<5(0,(a + l - 3 y ) ( « ; o + 2tyi )0 - 30 ) (5.14) 2 1 2 20 u -^902? 1 + -g iz z + ... 2 2 2 s 3 s 2 1 (5.12) The above is called the restricted equation. It describes the dynamics on the center manifold. Since the center manifold is locally exponentially attracting, any nontrivial local dynamics must lie on the center manifold and thus be captured by (5.12). Using a series of near identity coordinate transformations, (5.12) may be transformed into 1 z = iuz + h(0)z z + ... , Zi(0) = 2~j23&(*02o0ii +ug i) 2 (5.15) Truncating the above gives: z = iuz + h{0)z z 2 (5.16) Provided that Zi(0) ^ 0, the above is topologically equivalent to (5.12) and thus captures the local nontrivial dynamics (algebraic decay or growth) at the Hopf bifurcation point. If /i(0) < 0 then we have algebraic decay at the Hopf point and the bifurcation is supercritical. If li(0) > 0, then we have algebraic growth at the Hopf point and the bifurcation is subcritical. If ^(0) = 0, then we are at a generalized Hopf point. To determine /i(0), equations (5.4), (5.10) and (5.11) are solved numerically using M A P L E . Details are given in the Appendix. 51 Chapter 5. Bifurcation Analysis 0.8 H 0.6 o o o o 0.4 0.28 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 (a) 0.2 0.28 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 (b) 0.8 0.6 0.4 H O O 0.2 0.28 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 h (c) (d) Figure 5.1: (a)-(c) show the amplitude of oscillatory pulses found in numerical simulations for the parameter values a) e = 0.115 b) e = 0.09 c) e = 0.02. (d) shows the first Liapunov coefficient as a function of e, as computed using centre manifold reduction. The results in d) predict that the generalized Hopf point should occur at about e = 0.09, in agreement with the results of (a)-(c). Other parameter values: a — 0.35, b = 0.3. 52 Chapter 5. Bifurcation Analysis L2n 1.2 'OoOCBg 0 . 8 / 0.8 V0.6 V0.6 0.4 0.4 0.2 H 0.2 0.15 0.2 0.25 0.3 0.35 0.4 0.45 {/" 0.15 0.2 h 0.25 0.3 0.35 0.4 0.45 h (a) (b) 1.2 0.41 0.8 V0.6 0.4 0.2 0.085 0.15 0.2 0.25 0.3 0.35 0.4 0.45 -0.1 J h (c) (d) Figure 5.2: A second example illustrating the connection between the generalized Hopf point of (d) and the periodic solutions found in (a)-(c). Parameter values: a = 0.2, b = 0.17 and (a) e = 0.08 (b) e = 0.069 and (c) e = 0.065. 53 Chapter 5. Bifurcation Analysis 0.02 0.04 0.06 0.08 Figure 5.3: First Liapunov coefficient for Hopf bifurcation of a stationary front. For all values of e shown above, we have Z (0) < 0, meaning that the Hopf bifurcation is supercritical. Parameter values: a = 0.2, b = 0.1. We also computed l\(0) for other parameter value sets, for example a = 0.1 and b = 0.25 and found that once again Zi(0) < 0 always. x 54 Chapter 6 Summary and Discussion The phenomena of oscillatory fronts in inhomogeneous media was initially discovered in [18] and [8] by performing a perturbation analysis around a codimension two bifurcation point, refered to as an NIB bifurcation in [8]. In a homogeneous medium, the NIB bifurcation consists of a stationary front bifurcating into two traveling fronts through a pitchfork bifurcation. When a weak inhomogeneity is added another degree of freedom is introduced in parameter space. This leads to a rich variety of dynamical behaviors such as front reflections and oscillations. Here we have emphasized the connection between the Hopf bifurcation and oscillatory waves in inhomogeneous media. In the case of both the front and pulse oscillations, the Hopf bifurcation occurs with a change in the strength of the inhomogeneity. Thus, the oscillatory waves are seen as a direct consequence of inhomogeneity. A connection between oscillatory pulses in inhomogeneous media and a Hopf bifurcation had previously been discovered in [17] for a model in optics. In the area of oscillatory waves in inhomogeneous media, directions for future work include a number of theoretical and applied questions. On the mathematical side, there are a number of questions about the bifurcation theoretic nature of the waves which have not been fully answered. For example, are oscillatory waves always related to a Hopf 55 Chapter 6. Summary and Discussion bifurcation? Are oscillatory fronts always related to a supercritical hopf bifurcation? If so, can it be shown analytically? Are oscillatory fronts always related to an NIB point? What is the bifurcation of the pulse generator shown in fig. 4.4D? How is it related to the Hopf bifurcation of oscillatory pulses? Is the bifurcation of pulse generators related to a torus bifurcation, or a homoclinic bifurcation? If oscillatory pulses are related to a double-zero point (see fig. 3.4), where is the associated homoclinic bifurcation of pulses in the cubic model? Continuing with mathematical questions, the are many generalizations one could take with waves in inhomogeneous media. For example, what is the influence of an inhomogeneity moving with constant velocity (eg. z{x + ct) in (1.1)-(1.2))? What would happen if we had a two bump shaped inhomogeneity for z(x)7 Would two pulse generators interact? Do oscillatory waves exist in 2-d and 3-d models of excitable media, or on circular and spherical domains? Do patterns analogous to pulse generators exist in 2-d and 3-d excitable media? If so, would they only generate target waves? Generally speaking, given a simple unimodal inhomogeneity, how many different types of pulse generators could there be? As we have shown in fig. 4.8, there are at least two distinct pulse generator patterns in one dimensional excitable media. In terms of applications, the main question seems to be: can oscillatory waves (including pulse generators) be observed in experiments, or in nature? If so, what kind of inhomogeneous forcing would give rise to them? What role could they play in biological systems? Another interesting question would be: could oscillatory waves be used to control patterns? For example, to suppress undesirable patterns? Could introducing a localized inhomogeneity result in target waves which could dominate over other waves normally observed, such as spiral waves? 56 Bibliography [1] S. Amari. Dynamics of pattern formation in lateral-inhibition type neural fields. Biological Cybernetics, 27:77-87, 1977. [2] M . Bode. Front-bifurcations in reaction-diffusion systems with inhomogeneous parameter distributions. Physica D, 106:270-286, 1997. [3] P. C. Bressloff and S. E . Folias. Front bifurcations in an excitatory neural network. SIAM J. Appl. Math, 65:131-151, 2004. [4] P.C. Bressloff, S. E . Folias, A . Prat, and Y - X L i . Oscillatory waves in inhomogeneous neural media. Phys. Rev. Lett, 91:178101, 2003. [5] J. Carr. Applications of Centre Manifold Theory. Springer-Verlag, Berlin/New York, 1981. [6] M . C. Cross and P. C. Hohenberg. Pattern formation outside of equilibrium. Rev. Mod. Phys., 65:851-1112, 1993. [7] S. E. Folias and P. C. Bressloff. Breathing pulses in an excitatory neural network. SIAM J. Appl. Dyn. Syst, 3:378-407, 2004. [8] A . Hagberg, E . Meron, I. Rubinstein, and B . Zaltzman. Controlling domain patterns far from equilibrium. Physical Review Letters, 76:427-430, 1996. [9] M . Hendrey, K . Nam, P. Guzdar, and E. Ott. Target waves in the complex ginzburglandau equation. Physical Review E, 62:7627-7631, 2000. 57 Bibliography [10] H - J Krug, L . Pohlmann, and L . Kuhnert. Analysis of the modified complete oregonator accounting for oxygen sensitivity and photosensitivity of belousov-zhabotinsky systems. J. Phys. Chem., 94:4862-4866, 1990. [11] Y . Kuznetsov. Elements of Applied Bifurcation Theory. Springer, New York, 1995. [12] Y . - X . L i . Tango waves in a bidomain model of fertilization calcium waves. Physica D, 186:27-49, 2003. [13] V . Petrov, Q. Ouyang, G . L i , and H . Swinney. Light-induced frequency shift in chemical spirals. J. Phys. Chem., 100:18992-18996, 1996. [14] D. Pinto and B . Ermentrout. Spatially structured activity in synaptically coupled neuronal networks: I. traveling fronts and pulses. SIAM J. Appl Math., 62:206-225, 2001. [15] A . Prat and Y . - X . L i . Stability of front solutions in inhomogeneous media. Physica D, 186:50-68, 2003. [16] A . Prat, Y . - X . L i , and P.C. Bressloff. Inhomogeneity-induced bifurcation of stationary and oscillatory pulses. Physica D, 202:177-199, 2005. [17] J. Rubin and C . K . R . T . Jones. Bifurcation and edge oscillations in the semiconductor fabry-perot interferometer. Optics Communications, 140:93-98, 1997. [18] P. Schutz, M . Bode, and H.-G. Purwins. Bifurcation of front dynamics in a reactiondiffusion system with spatial inhomogeneities. Physica D, 82:382-397, 1995. 58 Appendix A M A P L E Programs A l l programs described below were executed in M A P L E 7. Approximate solutions to the boundary value problems below were found numerically using M A P L E ' s "dsolve" algorithm, which consists of a trapezoid finite difference scheme with a Richardson extrapolation method enhancement scheme. 1 S t a t i o n a r y fronts For stationary fronts, we solve the equation v" + v — v — ax = 0. This equation is 3 equivalent to (2.27) after an appropriate change of variables. The solution to v" + v — s v — ax — 0 is odd and so we can solve on the half space [0, a] and take v(0) = 0 as a 3 boundary condition. The other boundary condition, v(a) = vr, is requiring that v(a) be the numerically obtained solution of v — v — a a = 0. 3 vr:=fsolve(v-v" 3-sigma*a=0,v=-100..-le-10); eq:=((2-c)" 2)*diff(v(x),x,x)+v(x)-v(x)' 3-sigma*x/((2-c)2)=0; ans:=dsolve(eq,v(0)=0,v(a)=vr, numeric,abserr=0.1e-5, maxmesh=512,continuations): 59 Appendix A. MAPLE Programs 2 Eigenvalues for stationary fronts eq2 below is equivalent to the eigenvalue problem (3.8) after a change of variables. eq2 is solved simultaneously with eq as defined above. The key to finding the leading eigenvalue is to find an even nodeless eigenfunction. Due to the evenness of the eigenfunction, we solve on the half space [0, a] and take a zero-flux boundary condition at x — 0. To ensure that the eigenfunction is nodeless, the easiest way is to plot it. eq2:=((2-c)"2)*difT(y(x),x,x) + (l-3*v(x)"2)*y(x) = chi*y(x); ans:=dsolve(eq,eq2,v(0)=0,v(a)=vr,D(y)(0)=0,y(0)=l,y(a)=0, numeric, abserr=0.1 e- 5, maxmesh=512, continuations): odeplot(ans,[x,y(x)],0..a); 3 Stationary pulses Since pulses are even solutions, they are solved on the half space [0, a] with zero flux boundary conditions at x = 0. Note how v (0), denoted vO below, is treated as an s input parameter and that h is solved for by the program. This enables one to compute bifurcation diagrams such as 2.5, where several pulses exist for a given value of h. f:=v- > v*(l-v)*(v-alpha); eq3:=((2-c)' 2)*difF(v(x),x,x)+f(v(x))-b*v(x)+ h*exp(-((x/(sigma*(2-c))) ~ 2)/2)=0; ans:=dsolve(eq3,v(0)=vO,v(a)=0,D(v)(0)=0, numeric,abserr=1.0e-5, maxmesh=512,continuation=c): 4 Eigenvalues for stationary pulses eq3 and eq4 are solved simultaneously, as with the case of fronts. 60 Appendix A. MAPLE Programs eq4:=((2-c) ~ 2)*diff(y(x),x,x) + (2*(alpha+l)*v(x)-3*v(x)" 2-alpha)*y(x) =chi*y(x); ans:=dsolve(eq3,eq4,v(0)=vO v(a)=0,D(v)(0)=0,D(y)(0)=0,y(0)=l y(a)=0, ) I numeric,abserr=0.1e-5,maxmesh=512,output=listprocedure,continuation=c); 5 First Liapunov coefficient We provide the program for finding the Liapunov coefficient in the case of the Hopf bifurcation of pulses. The program for fronts is almost identical. Below we assume that the stationary pulse v(x) and its leading eigenfunction y(x) are contained in "ans", as in the section above. Note that we solve for the center manifold coefficient 1002(2) in (5.11) by splitting it up into real and imaginary parts. V:=subs(ans,v(x));Y:=subs(ans,y(x));setattribute(V,_syslib); setattribute(Y,_syslib); ch:=subs(ans,chi) (sigma); omega:=sqrt(b/ch-1); A:=(omega" 2+l)/(l+4*omega" 2) ;B:=(6*omega" 3)/(l+4*omega ~ 2); kd:=evalf(Int(Y(x)~2,x=0..a)); Y N : = x - > Y(x)/((2*kd)" (1/2)); k:=evalf(Int((alpha+l-3*V(x))*YN(x)~ 3,x=0..a)); eq5:=ch*(A*vr(x)-B*vi(x)) =diff(vr(x),x,x) +(2*(alpha+l)*V(x)-3*V(x) '2-alpha)*vr(x) +2*((alpha+l-3*V(x))*YN(x) ~ 2-2*k*YN(x)); eq6:=ch*(B*vr(x)+A*vi(x)) = diff(vi(x),x,x) +(2*(alpha+l)*V(x)-3*V(x) "2-alpha)*vi(x); eq7:=difT(g(x),x,x)+(2*(alpha-)-l)*V(x)-3*V(x)"2-alpha)*g(x)-b*g(x) +2*((alpha+l-3*V(x))*YN(x)" 2-2*k*YN(x))=0; ans2:=dsolve(eq7,D(g)(0)=0,g(a)=0,numeric,output=listprocedure); G:=subs(ans2,g(x)); ans3:=dsolve(eq5,eq6,D(vr)(0)=0,D(vi)(0)=0,vr(a)=0,vi(a)=0, 61 Appendix A. MAPLE Programs numeric, output=listprocedure); VR:=subs(ans3,vr(x)); VI:=subs(ans3,vi(x)); g20:=(l-I/omega)*(l/ch)*(2*k); g21:=(l-I/omega)*(l/ch)*(2*evalf(Int((alpha+l-3*V(x))*(VR(x)-I*VI(x) +2*G(x))*YN(x)' 2-3*YN(x)' 4,x=0..a))); ll:=(l/(2*omega"2))*Re(I*(g20)'2+omega*g21); 6 Domain sizes used Table A . l : Domain sizes used in M A P L E calculations. The numbers below refer to the size of the half space [0,a]. Figure Domain size 2.4, 2.5, 3.4, 5.1d 15 3.1 10 5.2d 20 62 Appendix B Numerical Methods A l l numerical simulations were performed using a finite difference method in space, an Euler method in time and zero-flux boundary conditions at both ends. Domains are taken to be [—a, a], with 2a being the domain size. The integration time step is given by At = 0.2e(2a/n) , where n is the number of mesh points. The domain sizes for simulations 2 of oscillatory pulses were chosen large enough so that v(±a, t) ~ 0 and w(±a, t) ~ 0. The domain sizes for simulations of oscillatory fronts were chosen large enough so that v(±a,t) fa constant and w(±a,t) fa constant. Table B . l : Domain size and mesh points for numerical simulations. Figure Domain size Mesh points 4.2 5 >/I6 2000 4.4A,B,C 7.5 VlO 2000 4.1 20 1000 5.1b 30 1000 5.1a,c 7.5 v/10 2000 4.4D 15 v/10 2000 5.2a,b,c 40 1000 63
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Oscillatory waves in an inhomogeneous excitable medium
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Oscillatory waves in an inhomogeneous excitable medium Prat, Alain 2005
pdf
Page Metadata
Item Metadata
Title | Oscillatory waves in an inhomogeneous excitable medium |
Creator |
Prat, Alain |
Date Issued | 2005 |
Description | An excitable medium is a nonlinear dissipative dynamical system able to sustain undamped wave propagation. Examples of excitable media include cardiac tissue and neural media. This thesis investigates a class of novel dynamical phenomena in excitable media, known collectively as oscillatory waves. The examples of oscillatory waves considered here are Tango waves, breathers and pulse generators. Tango waves, also known as oscillatory fronts, are characterized by a back-and-forth motion of a wavefront. Breathers, also known as oscillatory pulses, are characterized by an expanding and contracting motion of a pulse-shaped pattern. Pulse generators are breathers that periodically emit traveling pulses. These oscillatory waves are caused by the presence of either a timeindependent forcing or spatial variation in a model parameter. The main point of this thesis is that oscillatory fronts or pulses can be understood as emerging from the Hopf bifurcation of a stationary front or pulse, respectively. The secondary point of this thesis is to investigate secondary bifurcations connected to the Hopf bifurcation. Possible applications of oscillatory waves will also be discussed. Analytical and numerical methods are used to compute the equilibrium solutions and the spectrum of their linearizations for piecewise linear and cubic Fitzhugh-Nagumo models of excitable media. Combining this with numerical simulations allows us to link oscillatory waves with a Hopf bifurcation. A center manifold reduction allows us to also connect the oscillatory waves to a generalized Hopf bifurcation. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2009-12-15 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0092145 |
URI | http://hdl.handle.net/2429/16670 |
Degree |
Master of Applied Science - MASc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2005-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2005-0601.pdf [ 2.56MB ]
- Metadata
- JSON: 831-1.0092145.json
- JSON-LD: 831-1.0092145-ld.json
- RDF/XML (Pretty): 831-1.0092145-rdf.xml
- RDF/JSON: 831-1.0092145-rdf.json
- Turtle: 831-1.0092145-turtle.txt
- N-Triples: 831-1.0092145-rdf-ntriples.txt
- Original Record: 831-1.0092145-source.json
- Full Text
- 831-1.0092145-fulltext.txt
- Citation
- 831-1.0092145.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0092145/manifest