UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Analysis of models of bursting electrical activity in pancreatic beta cells De Vries, Gerda 1995

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1995-060969.pdf [ 3.9MB ]
Metadata
JSON: 831-1.0079910.json
JSON-LD: 831-1.0079910-ld.json
RDF/XML (Pretty): 831-1.0079910-rdf.xml
RDF/JSON: 831-1.0079910-rdf.json
Turtle: 831-1.0079910-turtle.txt
N-Triples: 831-1.0079910-rdf-ntriples.txt
Original Record: 831-1.0079910-source.json
Full Text
831-1.0079910-fulltext.txt
Citation
831-1.0079910.ris

Full Text

ANALYSIS OF MODELS OF BURSTING ELECTRICAL ACTIVITY IN PANCREATIC BETA CELLS By Gerda de Vries B. Math. (Applied Math. and Comp. Science), University of Waterloo, 1989 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF MATHEMATICS INSTITUTE OF APPLIED MATHEMATICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA August 1995 © Gerda de Vries, 1995 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of M o1Lry’cJs The University of British Columbia Vancouver, Canada Date 31, k19 DE-6 (2/88) Abstract Pancreatic 18-cells are responsible for producing and secreting insulin, a hormone which is essential in the regulation of blood glucose level. In the presence of glucose, 3-cells exhibit periodic bursting electrical activity (BEA) consisting of active and silent phases. During the active phase, the membrane potential of these cells oscillates rapidly, whereas the membrane potential changes slowly during the silent phase. The plateau fraction, which is the ratio of the active phase duration to the period of the bursting phenomenon, has been correlated to the insulin-glucose response of these cells. There are a number of mathematical models describing BEA in pancreatic /3-cells. In the first part of this thesis, the class of first-generation models which consist of three coupled first-order ordinary differential equations is analyzed. Numerical and analytical techniques are presented for determining the approximate plateau fractions from the model equations as a function of the glucose-dependent parameter. A consistent nondimensionalization of the model equations is proposed which permits all of the models to be written in a standard form. The equations then are separated into a second-order fast subsystem and a first-order slow subsystem. A bifurcation analysis of the fast subsystem in which the slow variable is treated as the bifurcation parameter reveals that the transition from the active phase to the silent phase occurs near a homo clinic bifurcation. The use of several numerical and analytical methods is demonstrated for the determination of the approximate location of this bifurcation, which is needed in the computation of the approximate silent and active phase durations and, subsequently, the approximate plateau fractions. 11 Leading-order problems for the silent and active phases are obtained using a combina tion of singular perturbation and multiple scales analyses and averaging techniques. The corresponding leading-order silent and active phase durations and the resulting leading- order plateau fractions are reduced to quadrature. Finally, the approximations are com pared with the respective exact (numerically computed) silent and active phase durations and plateau fractions over a range of the glucose-dependent parameter for which the mod els exhibit BEA. In the remainder of this thesis, a detailed study of a polynomial analog model of BEA in pancreatic /3-cells is carried out. Depending on the values of the model parameters, the model exhibits a wide variety of oscillatory solution behaviour, including types of bursting observed in excitable cells other than pancreatic /3-cells. A bifurcation map in the fast subsystem parameter space is produced by computing curves which represent codimension-2 bifurcations. These curves bound regions on the map within which bifurcation diagrams of the fast subsystem are qualitatively the same. The importance of the bifurcation map is that it shows the relationships between the various types of oscillatory behaviour and, hence, provides a basis for an extension of the classification of bursting oscillations. Since the analog model consists of polynomial functions, the curves on the bifurcation map can be derived analytically. 111 Table of Contents Abstract ii Table of Contents iv List of Tables viii List of Figures ix Acknowledgements xiv Dedication xv 1 Introduction 1 1.1 Thesis Outline 4 1.2 The Mathematics of Cell Electrophysiology 8 2 Biophysical Models and Preliminary Analysis 14 2.1 The First Mathematical Model for Bursting Electrical Activity in the Pan creatic Beta Cell 15 2.1.1 The Biophysical Model 15 2.1.2 The Chay-Keizer Model 16 2.1.3 Rinzel’s Fast/Slow Subsystem Analysis 21 2.2 Further Developments 26 2.2.1 The Chay (1986) Model 28 2.2.2 The Himmel-Chay Model 31 iv 2.2.3 The Chay (1987) and Chay-Kang Models 32 2.2.4 The Sherman-Rinzel-Kejzer Model 33 2.2.5 The Chay-Cook Model 34 2.3 Nondjmensionalizatjon and Problem Reformulation 35 3 Escape from the Active Phase 41 3.1 Approximation of z8 with the AUTO Software Package 43 3.2 Approximation of Zesc with Bisection Techniques 44 3.3 Approximation of Zesc with Melnikov’s Method 48 3.3.1 Reformulation of the Fast Subsystem as a Near-Integrable System 49 3.3.2 Derivation and Evaluation of the Melnikov Integral 56 3.4 Approximation of z with the Fredhoim Alternative 64 3.4.1 Derivation of the Equations 64 3.4.2 Numerical Solution 68 3.4.3 Why The First-Order Correction is Zero for Models with p = 1 . . 71 3.5 Comparison of the Different Methods 74 4 Determination of the Plateau Fraction 77 4.1 Range of Values of 9 78 4.2 Definition of the Exact Silent and Active Phase Durations and the Plateau Fraction 80 4.3 Leading-Order Approximation of the Silent Phase Duration 82 4.4 Leading-Order Approximation of the Active Phase Duration 87 4.5 Leading-Order Approximation of the Plateau Fraction 96 4.6 Examination of the Transition Phases 100 5 Multiple Bifurcations in a Polynomial Analog Model 103 V 5.1 Pernarowski’s Analog Model and Various Solution Behaviours 105 5.2 Symmetry Considerations 112 5.3 Preliminary Dynamical Systems Analysis 113 5.4 Constraints on the Hopf Bifurcations 114 5.4.1 Number of Hopf Bifurcations 116 5.4.2 Location of the Hopf Points 118 5.4.3 Criticality of the Hopf Bifurcations 121 5.5 Constraints on the Homoclinic Bifurcations 124 5.5.1 The Melnikov Distance Functions 127 5.5.2 Homoclinic Orbits at the Local Extrema of the G = 0 Nulicline 131 5.5.3 Double Homoclinic Orbits 135 5.5.4 Neutral Homoclinic Orbits 137 5.5.5 Coalescing Homoclinic Orbits 140 5.6 Constraints on the Saddle-Node of Periodics Bifurcations 148 5.7 The Complete Bifurcation Map 149 5.8 Two-Parameter Bifurcation Diagrams 154 5.8.1 Numerical Two-Parameter Bifurcation Diagrams 155 5.8.2 Analytical Two-Parameter Bifurcation Diagrams 159 5.9 On the Classification of Bursting Oscillations 161 5.10 Effect of e on the Solution Behaviour of the Full System of Equations . . 163 6 Conclusions and Discussion 167 Bibliography 176 Appendices 184 vi A Model Equations A.1 The Reduced Chay-Keizer Model A.2 The Chay (1986) Model A.3 The Himmel-Chay Model A.4 The Chay-Kang Model A.5 The Sherman-Rinzel-Keizer Model A.6 The Chay-Cook Model B Bifurcation Diagrams for the Polynomial Analog Model 205 184 184 188 192 196 199 202 vii List of Tables 2.1 Comparison of the first-generation models consisting of three ordinary dif ferential equations 29 3.1 Approximate values of ZHC and, subsequently, Zegc obtained with AUTO, bisection, Melnikov’s method, and the Fredhoim alternative method . 3.2 Comparison of and 6, showing that 6 >> e > 0 5.1 Summary of the qualitative features of the Hopf, homoclinic, and SNP bifurcations exhibited by bifurcation diagrams of the fast subsystem in the labelled regions on the complete bifurcation map 5.2 Summary of the qualitative features of the homoclinic and SNP bifurca tions and the nature of the detached branches of periodic orbits exhibited by bifurcation diagrams in regions F1-F4 on the bifurcation map 187 191 195 198 201 204 44 56 152 153 A.1 Variables and parameters of the reduced Chay-Keizer model A.2 Variables and parameters of the Chay (1986) model . . . A.3 Variables and parameters of the Himmel-Chay model . . A.4 Variables and parameters of the Chay-Kang model A.5 Variables and parameters of the Sherman-Rinzel-Keizer model A.6 Variables and parameters of the Chay-Cook model viii List of Figures 1.1 Regular pattern of bursting electrical activity induced by 10 mM glucose in a mouse pancreatic /3-cell 2 1.2 Correlation of the insulin release rate with the plateau fraction 5 1.3 Three-dimensional view of a small section of a cell membrane 9 1.4 Schematic drawing of a channel protein or ion channel in its open and closed states 10 2.1 Graphs of the functions m(V), h(V), n(V), and r(V), describing the channel kinetics of the Chay-Keizer model 21 2.2 Numerical solution of the reduced Chay-Keizer model 22 2.3 133(V; Cat) and the curve of steady states of the fast subsystem for the reduced Chay-Keizer model 24 2.4 Bifurcation diagram, generated with AUTO, of the fast subsystem of the reduced Chay-Keizer model 25 2.5 Bifurcation diagram with the Ca, nullcline and the projection of the nu merical solution of Figure 2.2 superimposed 27 2.6 Illustration of the two types of bifurcation diagrams exhibited by the fast subsystem of the nondimensional /3-cell models 40 3.1 Typical bifurcation diagram of the fast subsystem with the numerical so lution of the full system of equations superimposed (I) 42 3.2 Solution trajectories of the fast subsystem 46 ix 3.3 Phase portraits of the fast subsystem showing the stable and unstable orbits of the saddle point 47 3.4 Cumulative effect of the damping term in the active phase, for models withp=1 50 3.5 Comparison of the solutions of the full system of equations and the same system of equations where E(u, z, z) has been replaced by the first two terms of its Taylor series about ü = 0 52 3.6 Cumulative effect of the damping term in the active phase, for models withp1 (I) 53 3.7 Cumulative effect of the damping term in the active phase, for models withp1(II) 55 3.8 Phase portrait of the unperturbed fast subsystem 57 3.9 Geometrical representation of the Melnikov distance 59 3.10 Graphs of the Melnikov distance as a function of z 63 4.1 Dependence of the -nullcline on /3 79 4.2 Typical bifurcation diagram of the fast subsystem with the numerical so lution of the full system of equations superimposed (II) 81 4.3 Dependence of the number of spikes per burst on /3 83 4.4 Comparison of the solution of the full system of equations and the solution of the leading-order silent phase problem 85 4.5 Comparison of the exact and leading-order silent phase durations . . . 88 4.6 Period To(z) and averaged calcium current 1h(z) 93 4.7 Comparison of the solution of the full system of equations and the solution of the leading-order active phase problem 95 4.8 Comparison of the exact and leading-order active phase durations . . . 97 x 4.9 4.10 4.11 The G = 0 nulicline of the fast subsystem Square-wave bursting Tapered bursting Bursting with an undershoot The bifurcation map: TB and CHB curves Parabolic-amplitude bursting The bifurcation map: HBL curves The bifurcation map: DHB curve 5.9 Phase portraits of the fast subsystem showing the stable and unstable or bits of the saddle points and illustrating the existence of a large homoclinic orbit 5.10 Illustration of a series of SNIC bifurcations at the local minimum of the G 0 nullcline, resulting in the deformation of a right homoclinic into a large homoclinic orbit 5.11 The bifurcation map: SNL curves 5.12 The bifurcation map: DSL curve 5.13 Bifurcation diagrams illustrating that the periodic orbits near a homoclinic orbit can be either stable or unstable 5.14 The bifurcation map: NSL curves 98 99 102 Comparison of the exact and leading-order plateau fractions Transition phase duration as a percentage of the burst period Projection of the numerical solution of the full system of equations onto the (u, t) phase plane with the leading-order approximations of the transition phases superimposed 106 108 109 110 118 119 121 124 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 126 133 134 136 137 139 xi 5.15 Bifurcation diagrams illustrating the existence of coalescing homoclinic orbits 141 5.16 The bifurcation map: CSL curves 142 5.17 Phase portrait of the fast subsystem showing that a large periodic orbit can exist in addition to right and left periodic orbits for certain values of the model parameters 144 5.18 Bifurcation diagrams illustrating the existence of a branch of periodic or bits which is detached from the main branch for certain values of the parameters 145 5.19 Bifurcation diagrams illustrating that large left and large right homoclinic orbits do not coexist 147 5.20 The bifurcation map: CSNP curve 149 5.21 The complete bifurcation map 150 5.22 Enlarged view of a section of the bifurcatioll map 151 5.23 Two-parameter bifurcation diagram generated with AUTO for ii = 1.5. 155 5.24 Two-parameter bifurcation diagrams generated with AUTO for ?2 = 0.8 and ü=0 158 5.25 Analytical two-parameter bifurcation diagram for i% = 1.5 and comparison with AUTO results 160 5.26 Illustration of an apparent paradox: the bifurcation diagram of the fast subsystem is that of a tapered burster, yet the numerical solution of the full system of equations exhibits square-wave bursting 165 5.27 Resolution of the apparent paradox 166 B.1 Bifurcation diagram for region A 205 B.2 Bifurcation diagram for region Bi 206 xii B.3 Bifurcation B.4 Bifurcation B.5 Bifurcation B.6 Bifurcation B.7 Bifurcation B.8 Bifurcation B.9 Bifurcation B.10 Bifurcation B.11 Bifurcation B.12 Bifurcation B.13 Bifurcation B.14 Bifurcation B.15 Bifurcation B.16 Bifurcation B.17 Bifurcation B.18 Bifurcation B.19 Bifurcation B.20 Bifurcation B.21 Bifurcation B.22 Bifurcation B.23 Bifurcation B.24 Bifurcation B.25 Bifurcation B.26 Bifurcation diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for diagram for region B2 206 region B3 207 region B4 207 region B5 208 region Cl 209 region C2 210 region C3 210 region C4 211 region C5 211 region C6 212 region Dl 213 region D2 213 region D3 214 region D4 214 region D5 215 region El 216 region E2 217 region E3 217 region E4 218 region E5 218 region Fl 219 region F2 219 region F3 220 region F4 220 xlii Acknowledgements First and foremost I thank my supervisor, Robert Miura, for his continued guidance and encouragement during the course of this research, both in times of frustration and success. Appreciation also goes to Leah Edelstein-Keshet for her immeasurable moral support throughout the years. Thanks also to Mark Pernarowski for so gracefully allowing me to use his work as a springboard. xiv foar myn leave en moediche heit en mem xv Chapter 1 Introduction Cells of the islets of Langerhans in the pancreas secrete several important hormones in order to maintain the glucose level in the blood within a narrow operating range. The two most important hormones are insulin, which is produced and secreted by /3-cells, and glucagon, which is produced and secreted by cr-cells. Insulin helps to decrease the blood glucose level, whereas glucagon helps to increase the blood glucose level. The disease called diabetes mellitus is a result of the inability of the body to properly utilize and/or produce insulin. If left untreated, diabetes can lead to serious medical complications. In a short preliminary paper in 1968 [24], Dean and Matthews provided evidence that the /3-cell is electrically excitable. In their more detailed follow-up paper in 1970 [25], they reported that /3-cells exhibit a characteristic periodic pattern of electrical activity in the presence of glucose, namely bursting electrical activity (BEA), an example of which is shown in Figure 1.1. Bursting electrical activity is a periodic phenomenon, during which the membrane potential of the cell undergoes silent and active phases. The silent phases are characterized by a slowly changing membrane potential and the active phases are characterized by a rapidly oscillating membrane potential. The term ‘burst’, referring to the repetitive spiking activity during the active phase, was coined by Dean and Matthews. During subsequent studies, it was found that the bursting pattern depends highly on the concentration of external glucose. When the glucose concentration is low (<5 mM), the membrane potential remains at a resting value of approximately -68 mV [3]. When the glucose concentration is increased, the /3-cell membrane exhibits BEA. As 1 Chapter 1. Introduction 2 mV -20 silent phase Figure 1.1: Regular pattern of bursting electrical activity induced by 10 mM glucose in a mouse pancreatic j3-cell. (Adapted from Ozawa and Sand [66].) the concentration is increased further, the duration of the active phase relative to the burst period increases also. That is, the higher the glucose concentration, the longer the active phases become and the shorter the silent phases [6, 8, 15, 23, 62]. When the glucose concentration is raised beyond 16.6 mM, the bursts disappear and the membrane potential exhibits continuous spiking. Insulin release induced by the presence of glucose occurs over the same concentration range [53]. Moreover, the rate of insulin release and the duration of the active phase are strongly correlated experimentally [4, 61, 63]. Therefore, in order to understand the regulation of insulin release, it is necessary to understand the electrophysiological properties of the /3-cell. Since the landmark work by Dean and Matthews, this has been a subject of intense investigation. The first biophysical model for BEA was proposed by At.water et al. [6]. They believed that three types of ion channels in the membrane of the /3-cell played an important role in producing the electrical activity and that the intracellular calcium concentration was the agent controlling the switch from the active phase to the silent phase and vice versa. Based on the ideas of Atwater et al., Chay and Keizer [18] developed the first mathe matical model of BEA in 1983. Their model was a modification of the classical Hodgkin active phase Chapter 1. Introduction 3 and Huxley model [47] for electrical activity in the squid giant axon. Although little electrophysiological data, such as channel densities, activation curves, or time constants, were available at that time, their model was very successful in simulating BEA. As new channels are discovered, or new data becomes available, the biophysical the ories are modified, resulting in a wealth of mathematical models of BEA. In this thesis, we restrict our attention to deterministic models. These models are assumed to represent the average behaviour of a single ,8-cell clustered together in an islet or in a petri dish. For single isolated /3-cells, the stochastic opening or closing of single channels has an important influence on the membrane potential. As a result, single isolated cells do not produce bursts [14]. However, /3-cells are known to be in close electrical contact with their neighbours through gap junctions [53]. For tightly coupled clusters of cells, the electrical effects of single-channel events are shared by all the cells. When the cluster is large enough, the effects of single-channel openings are less of a factor, stochastic noise disappears, and the islet behaves deterministically. For models addressing the issues related to stochasticity and the effects of coupling through gap junctions, see [86, 87, 89]. On the basis of the underlying physiological assumptions, the deterministic models that have been reported in the literature can be divided into two major groups. The class of models which we call first-generation models [13, 14, 16, 17, 44, 54, 87], developed from 1983 to 1989, are based on the assumption that the slow variable in the system is the intracellular calcium concentration. The essential differences between these models are the specific ionic currents that are thought to be important for producing BEA. Since 1989, there has been a development of second-generation models [15, 55, 88]. These models are characterized by the fact that they no longer depend on the assump tion that a slowly varying intracellular calcium concentration is responsible for burst termination. In the first two models, advanced by Chay [15] and Keizer and Smolen [55], the slow process is replaced by the slow voltage-dependent inactivation of the calcium Chapter 1. Introduction 4 current, based on a hypothesis put forward by Satin and Cook [83, 84]. In these models, the intracellular calcium concentration no longer oscillates slowly, but varies rapidly and in synchrony with the rapidly oscillating membrane potential. Recent experiments by Valdeolmillos et al. [91] and Santos et al. [81], in which the intracellular calcium concen tration is measured using fluorescence techniques, indeed show that intracellular calcium does not accumulate slowly during the active phase, but rather rises rapidly at the onset of the active phase, and reaches a steady state level after the first few spikes. In [88], Smolen and Keizer note that the slow voltage-dependent inactivation hypothesis is not compatible with the available data. Instead, they propose that the slow process depends on ATP (adenosine triphosphate) or ADP (adenosine diphosphate) concentration. The identity of the slow variable remains unknown. Similarly, the ionic currents playing an important role in producing BEA and related phenomena remain under inves tigation. For example, the implications of the hypothesis that a calcium release activated current (CRAC) is present in /3-cells is examined in a very recent paper by Bertram et al. [11]. Rather than focusing on the controversies surrounding these models, the goal in this thesis is to present techniques to study problems associated with BEA, and to add to the mathematical understanding of BEA in pancreatic /3-cells in particular, and the bursting phenomenon in general. In Section 1.1, we outline the organization of this thesis and introduce the major issues discussed in the following chapters. The necessary electrophysiological background for this thesis is provided in Section 1.2. 1.1 Thesis Outline This thesis can be divided into two parts. In the first part, consisting of Chapters 2, 3, and 4, we concentrate on biophysical models of bursting electrical activity in pancreatic Chapter 1. Introduction 5 0 C z C 0 a Figure 1.2: Correlation of the rate of insulin release from isolated islets with the plateau fraction as functions of glucose concentration. (From Ozawa and Sand [66].) /3-cells and develop analytical techniques to study a physiologically important problem, namely the determination of the plateau fraction, which is the ratio of the active phase duration to the period of the bursting phenomenon. In the second part, consisting of Chapter 5, we study the bursting phenomenon in a much broader context, concentrating on the bifurcation behaviour of a polynomial analog model of bursting electrical activity. The plateau fraction is physiologically relevant because it appears to be strongly cor related experimentally [4, 61, 63] to the rate of insulin release, as shown in Figure 1.2. This correlation suggests that the relative lengths of the active and silent phases are important in the regulation of insulin release. From a physiological point of view, viable models of pancreatic /3-cells are those that not only reproduce the observed electrical activity, but also predict a plateau fraction as a function of the glucose-dependent pa rameter(s) which is in agreement with experimental observations. Himmel and Chay [44] obtained plots of the plateau fraction as a function of the glucose-dependent parameter by numerically integrating the model equations for various values of the parameter and £ C C C C C C 5 10 15 20 25 Glucose concentration. mM Chapter 1. Introduction 6 extracting the plateau fraction from the solution. Pernarowski [67] and Pernarowski et al. [70] developed an analytical, and much more efficient, method for determining the plateau fraction in the context of the Sherman-Rinzel-Keizer model [87], based on formal perturbation techniques. In [67], Pernarowski notes that the method can be applied to all models which can be written in the same (perturbed Liénard) form as the Sherman Rinzel-Keizer model. The question that was not answered is whether other 3-cell models also can be written in this form. One of the aims of the research described in this thesis is to extend Pernarowski’s analysis to the majority of other models, i.e., to find a unified approach to calculating the plateau fraction. We restrict our attention to those first-generation /3-cell models that consist of three ordinary differential equations. Although these models have been superseded, their analysis remains fruitful for two reasons. First, all models, including the most recent second-generation models, are based on the assumption that there is a slow variable in the system that controls the termination of the burst during BEA. The identification of the slow variable is extremely important from a biophysical point of view. However, we believe that from a mathematical point of view, whether the slow process is governed by the intracellular calcium concentration, the slow voltage- dependent inactivation of the calcium current, the ATP or ADP concentration, or some other variable is irrelevant. We believe that the analysis presented in this thesis can be generalized to other models of BEA. Second, the bursting phenomenon is observed not only in pancreatic /3-cells, but also in the electrical behaviour of many nerve cells. Examples include hippocampal pyramidal neurons [50, 93], thalamic neurons [27], and AB neurons [41]. Furthermore, bursting has been observed in biochemical systems of enzyme reactions [26, 48], as well as in chemical systems, such as the Belousov-Zhabotinskii reaction [49]. Any new insights that can be gained from the analysis of the aforementioned models adds to the understanding of the bursting phenomenon in general. Chapter 1. Introduction 7 In Chapter 2, we introduce all first-generation models of BEA consisting of three ordinary differential equations. Most of the discussion is centered around the development of the first mathematical model, due to Chay and Keizer [18], which consists of five first- order differential equations, and the subsequent reduction to a third-order system [19]. The reasons as to why the Chay-Keizer model gives rise to the bursting phenomenon were provided by Rinzel [76]. We summarize his analysis in this chapter since it is fundamental to the work presented in the remaining chapters. In addition to the Chay Keizer model and the Sherman-Rinzel-Keizer model, the models being considered are due to Chay [13, 14], Chay and Cook [16], Chay and Kang [17], and Himmel and Chay [44]. Finally, we present a consistent nondimensionalization of the model equations so that all models can be written in the same standard form. Before being able to determine the plateau fraction from the model equations, we need to know how to predict the beginnings and ends of the silent and active phases. Based on the bifurcation analysis of the model equations, the end of the silent phase and the beginning of the active phase are easily determined a priori. However, predicting when the active phase ends and the silent phase begins is more difficult. In Chapter 3, we examine various numerical and analytical approaches to determine when the escape from the active phase occurs. The results of Chapter 3 are used in Chapter 4, where we derive leading-order expres sions for the silent and active phase durations and the plateau fraction, using perturbation analysis, including a multiple scales analysis of the active phase and averaging. The silent and active phase durations are shown to be in excellent agreement with the corresponding exact durations. The leading-order plateau fraction is shown to consistently overestimate the exact plateau fraction, due to neglecting the time spent in the transitions from the silent to the active phase and vice versa. We end the chapter and the part of the the sis concerned with the determination of the plateau fraction with a discussion of these Chapter 1. Introduction 8 transitions. In Chapter 5, we study bursting oscillations in a much broader context. To avoid the complex nonlinearities in the biophysical models, we use the polynomial model due to Pernarowski [68]. Depending on the parameter values, this model exhibits a wealth of different types of solution behaviour, including square-wave bursting which is the type of bursting observed in pancreatic ,8-cells. An examination of the bifurcation diagram cor responding to different types of solution behaviour reveals that the bifurcation diagrams themselves are undergoing bifurcations as parameters are varied. We create a bifurcation map in the fast subsystem parameter space, consisting of regions within which the bi furcation diagrams are qualitatively the same. The bifurcation map is significant in that it shows the relationships between the various types of solution behaviour and, hence, provides a basis for a comprehensive classification scheme of bursting oscillations. Finally, in Chapter 6, we make some general remarks about the work presented in this thesis and discuss possibilities for future work. 1.2 The Mathematics of Cell Electrophysiology All living cells are enclosed by a thin membrane, which is made up of phospholipid and protein molecules, as shown in Figure 1.3. The phospholipids provide the basic structure of the membrane and the proteins carry out specific functions. Each phospholipid consists of a hydrophyllic (water-loving, polar) head and two hy drophobic (water-hating, nonpolar) fatty acid tails. To form the basic structure of the membrane, phospholipids line up spontaneously into two layers (bilayer) with the tails sandwiched between the heads (cf. Figure 1.3). Because of the special architecture of the bilayer, namely that the inside is strongly hydrophobic, the bilayer acts as an imperme able barrier to water-soluble molecules. Chapter 1. Introduction 9 iIipid Figure 1.3: Three-dimensional view of a small section of a cell membrane. (From Alberts et al. [1], page 275.) Embedded in the phospholipid bilayer are many different types of protein molecules (cf. Figure 1.3). These molecules carry out the specific functions of the membrane. Some proteins act as receptors for chemical signals, others are enzymes that catalyze reactions which are membrane-related, and still others act to transport specific molecules or ions into and out of the cell. The ability of proteins to control the flow of ions (and, therefore, electrical current) is the underlying reason for the electrical activity exhibited by many types of cells. The proteins which are concerned with the movement of ions have three very important characteristics. First, they form tiny pores and, hence, these proteins are usually called channel proteins or ion channels. Second, these channel proteins exhibit ion selectivity. That is, each type of ion channel permits a specific ion to pass through it, for example, potassium (Kj or sodium (Na+). Third, the channels are not always open. They have switches, or gates, which open only when stimulated. The stimulation can be a change in membrane potential, in which case the channel is said to be voltage-gated. Alternatively, the stimulation can be the binding of a ‘messenger’ molecule to a membrane receptor or /‘L) lipid molecule Chapter 1. Introduction 10 CLOSED OPEN EXTRACELLULAR SPACE I•, ‘‘Ilipid bilayer •0 oil / CYTOPLASM / aqueous pore ion-selective filter Figure 1.4: Schematic drawing of a channel protein or ion channel in its open and closed states. The outside of the channel is hydrophobic and interacts with the phospholipid bilayer. The inside is hydrophyllic and facilitates the passage of water-soluble ions. When the channel gate opens, the protein changes conformation to form a tiny aqueous pore. The pore is of atomic dimensions and it is here that the ion selectivity is determined. The location of the gate and the ion selectivity filter are shown here to be on the outside of the membrane. For most channels, this detailed information is still unknown. (From Alberts et al. [1], page 313.) the binding of an ion (corresponding to the terms ligand-gated and ion-gated channel, respectively). Figure 1.4 shows a schematic drawing of a gated ion channel. The driving force for ionic movement through open channels stems from the differences in the electric potential and the ionic concentrations across the cell membrane. The membrane potential is the voltage difference that exists across all biological membranes. It arises from an imbalance of electrical charges across the membrane, as well as ionic current flow through the membrane. In a resting cell, this imbalance is maintained by passive ionic diffusion through open channels and active, energy-consuming mechanisms such as ionic pumps. Large concentration differences across the membrane constitute the second factor influencing the movement of ions. Typically, potassium ions are at a much higher con centration inside the cell relative to the outside, and sodium, chloride, and calcium ions Chapter 1. Introduction 11 are at a much higher concentration outside the cell. Thus, there are two forces acting on any ion, namely an electrostatic force due to the membrane potential and a diffusional force due to the concentration differences. The equivalent electric potential difference (the Nernst potential, V,) required to maintain the concentrations for a specific ion on the two sides of the membrane is given by the Nernst equation RT C0 — —ln-—, (1.1) where R = universal gas constant (8.31 joules/mole/K), F = Faraday’s constant (96.487 coulombs/mmole), = absolute temperature (K), Z = valence (charge) of the ion, C0,C = outside and inside concentrations of the ion, respectively. That is, if the membrane potential is V, then the electrical gradient just balances the concentration difference for that ion and there is no net flow of the ion across the membrane. When the membrane is made more permeable to a given ion, the membrane potential tends to be driven towards the Nernst potential for that ion. As a result, practically any change in the membrane permeability, due to the opening of a channel in response to some stimulus, causes a change in the membrane potential. Resting membranes are typically permeable primarily to K ions. Therefore, the resting potential equilibrates to a value near the Nernst potential for K+, generally about -75 mV. Precisely how the membrane potential varies with the flow of ionic currents can be simulated through the use of mathematical models. The importance of the classical work Chapter 1. Introduction 12 of Hodgkin and Huxley [47] in this area cannot be overestimated. The underlying theory used to develop their model still is used extensively today. In particular, the membrane can be thought of as a leaky capacitor. The membrane current m thus is the sum of the capacitive current and all ionic currents, viz., dV = Cm + Iion, (1.2) where Cm is the more-or-less constant specific capacitance of the membrane, V is the membrane potential, T is the time, and I is the sum of all ionic currents, depending on the types of channels present in the membrane. In the present context, we ignore the contributions due to electrogenic pumps. Normally, m is zero, due to Kirchhoff’s Law, which essentially restates the law of conservation for charge. That is, the capacitive and ionic currents are in balance, CmZ ‘ion. (1.3) The above expression forms the starting point for all mathematical models of mem brane electrical activity following the Hodgkin-Huxley modelling theory. The models differ in the ionic currents included in Ij. An important assumption is that all ionic currents flow independently of each other. That is, each type of ion channel contributes a term to I. In practice, only the channels which are thought to be important in producing the desired electrical activity are included. The net driving force (due to the electrical gradient and concentration difference) for each ionic type is proportional to the difference between the membrane potential and the Nernst potential of that ion, V — V. The actual current carried by the ion not only depends on this driving force, but also on how easy or difficult it is for the ion to pass through the membrane. In particular, by Ohm’s Law, I = (1.4) ion Chapter 1. Introduction 13 where R0 is the resistance of the membrane to the ion. However, instead of resistance, physiologists prefer to use conductance, which is the reciprocal of resistance. Defining the conductance of the ion, gj0,,, to be gion = __, (1.5) 20fl (1.4) becomes I =g0,(V — (1.6) The conductance gjo in effect measures the ease with which the ion passes through the membrane or the permeability of the membrane to the ion. From experiments, it is known that g often depends on the state of the membrane (for example, the value of the membrane potential or whether ‘messenger’ molecules are bound to the membrane, etc.). For instance, for a voltage-gated channel, the conductance may increase (i.e., channel gates are switched to the open state) as the potential difference across the membrane decreases. Therefore, mathematical models also must include descriptions of the dependence of 9ion on the state of the membrane. There are two major difficulties with the Hodgkin-Huxley modelling theory that should be pointed out. First, it may not always be clear which ion channels are the important players in producing the desired electrical activity. In fact, it may well be that a particular ion channel hasn’t even been identified, let alone characterized as an important player. Second, the dependence of the conductance gjo on the state of the membrane may not be well understood. In the case of bursting electrical activity in pancreatic /9-cells, the controversies surrounding these difficulties have led, and continue to do so, to many different biophysical theories and corresponding mathematical models, as will become evident in the next chapter. Chapter 2 Description of the First-Generation Biophysical Models and Preliminary Analysis In this chapter, we introduce the first-generation models of bursting electrical activity in pancreatic /3-cells consisting of three ordinary differential equations. The biophysical theories upon which the models are based essentially differ in the types of channels that are thought to be important for producing the bursting activity and how these channels are thought to function. Therefore, the mathematical models mainly differ in those terms which represent the currents passing through the channels. However, all models have in common the assumption that there is a net efflux of Ca2+ from the cell during the silent phase and a gradual buildup of Ca2 during the active phase. In Section 2.1, we discuss the first model in detail. Although both the biophysical theory and corresponding mathematical model have been superseded, the discussion is useful, as it serves as a basis for understanding the bursting phenomenon in general. The remaining models of bursting electrical activity in pancreatic /3-cells which are considered in this thesis are introduced in Section 2.2. In Section 2.3, we nondimensionalize the models so that they can be written in a standard form and reformulate the equations for the purpose of the upcoming analysis in the following chapters. 14 Chapter 2. Biophysical Models and Preliminary Analysis 15 2.1 The First Mathematical Model for Bursting Electrical Activity in the Pancreatic Beta Cell We begin by discussing the first biophysical model of bursting electrical activity in pan creatic /3-cells, due to Atwater et al. [6] in Section 2.1.1. In Section 2.1.2, we discuss the corresponding first mathematical model, due to Chay and Keizer [18]. Since we exten sively use Rinzel’s fast/slow subsystem analysis [76] throughout this thesis, a summary of this analysis is provided in Section 2.1.3. 2.1.1 The Biophysical Model By the late 1970’s, there was a strong belief that three types of channels play an important role in the control of the electrical activity. These channels were a voltage-gated Ca2+ channel [60], a voltage-gated K channel [7, 42, 43], and aCa2-activ ted K channel [5, 75]. Based on these beliefs and the hypothesis that the calcium ion is a control agent for the membrane potential [5], Atwater et al. [6] proposed the first biophysical model explaining the bursting electrical activity. In particular, Atwater et al. proposed that changes in the membrane potential which are induced by the presence of external glucose can be explained in terms of changes in the intracellular calcium concentration, [Ca2+]. When no glucose is present, [Ca2’ji is maintained at a relatively high level. Thus, the Ca2-activ ted K channel is open, resulting in a membrane potential near the Nernst potential for K+. When glucose is present, the cell is induced to lower [Ca2+]j through the activation of energy-consuming processes. As [Ca2j, is lowered, theCa+activ ted K+ channel is inhibited, causing the membrane potential to increase. When [Ca2jis sufficiently reduced, theCa2+activ ted K+ channel has a low open probability and a rapid depolarization occurs. The voltage gated K+ and Ca2+ channels are then activated, resulting in a burst of spikes or action Chapter 2. Biophysical Models and Preliminary Analysis 16 potentials. During the spiking activity, Ca2 flows into the cell. As [Ca2j increases, the Ca2+activ ted K+ channel becomes active again, repolarizing the membrane and inhibiting a further influx of calcium. Once the cell has lowered [Ca2], the cycle starts again. 2.1.2 The Chay-Keizer Model The model of Atwater et al. [6] describes how the interaction of membrane processes (ionic currents) and the changing intracellular calcium concentration can produce burst ing electrical activity in a qualitative sense. To study this interaction quantitatively, Chay and Keizer [18] developed the first mathematical model in 1983. Chay and Keizer included the four basic features of the biophysical model, namely the three ionic channels described in the previous section and the dynamic calcium concen tration. Since much of the required electrophysiological data for18-cells was not available at that time, they used the classical Hodgkin and Huxley [47] model for electrical activity in the squid giant axon as a basis for their model. The most important modifications to the Hodgkin and Huxley model are the addition of the Ca2-activ ted K channel, the replacement of the inward sodium current by an inward calcium current, and the addition of a dynamical equation for the intracellular calcium concentration. Following the scheme introduced by Hodgkin and Huxley, the conductance of the voltage-gated K+ channel is written as = gKn4, (2.1) where K is the maximum conductance of the channel and is the fraction of K+ activation. Similarly, the conductance of the voltage-gated Ca2+ channel is written as gca = gcam3h, (2.2) Chapter 2. Biophysical Models and Preliminary Analysis 17 where ca is the maximum conductance, m is the activation parameter of the channel, and h is the inactivation parameter. These equations may be interpreted physically as follows. The K+ channel is assumed to be open if four similar, but independent, channel subunits are in the active state. The probability that one subunit is active is n; the probability that four are active is n”. Similarly, the Ca2+ channel is assumed to be open if three similar channel subunits are active (with probability m) and a different subunit is not active (with probability h). Since m, h, and n represent probabilities, their values range between 0 and 1. To model the Ca2+activ ted K+ channel, a scheme proposed by Plant [72] was adopted. This scheme assumes that Ca2+ binds to the channel according to first-order Michaelis-Menten kinetics (i.e., one Ca2 ion bound to the channel will activate it). Thus, the conductance for this channel is written as - Ca:/Kd gK,ca = gK,ca 1 + Ca/Kd’ (2.3) where gK,Ca is the maximum conductance of the channel, Ca, is the intracellular calcium concentration, and Kd is the dissociation constant for Ca2+ bound to the channel gate. The flow of other ions is reflected in a leakage current of non-specific nature. The leakage current is assumed to have a constant, relatively small, conductance, g. Based on (1.6) and (2.1)-(2.3), the ionic currents thus are Ia(V) = .oam3h(VVca), (2.4) IK(V) = KTt4(V — VK), (2.5) - a IK,Ga(V, Ca:) = gK,ca rz (V — VK), (2.6) lid + IL(V) = gL(V — VL), (2.7) where ‘Ca, ‘K, ‘K,Ca, and IL are the currents flowing through the voltage-gated Ca2+ and K+ channels, theCa2+activ ted K+ channel, and the leakage channel, respectively. Chapter 2. Biophysical Models and Preliminary Analysis 18 VCa, VK, and VL are the Nernst potentials for calcium, potassium, and the leakage, respectively. From equation (1.3), the rate of change for the membrane potential V then satisfies 4Tr2Cm = [Ia(V) + IK(V) + IK,Ca( V, Ca) + IL(V)], (2.8) where r is the radius of a single /3-cell. The factor 4rr2, representing the cell surface area, needs to be included because conductances typically are given in Siemens for the entire /3-cell rather than Siemens/unit area. The activation and inactivation variables m, h, and n, satisfy the relaxation equations introduced by Hodgkin and Huxley, namely, dm — m00(V)—m 29dT Tm(V) dli — h(V)—h 210 - Th(V) ‘ (. dii — n(V)—n 211dT - r(V) (. The voltage dependencies of the steady-state functions m, h, and nm, and the func tions Tm, Th, and T? corresponding to the time constants have exactly the same form as in the Hodgkin and Huxley model, but they are shifted along the voltage axis. So far, the Chay-Keizer model consists of four equations (2.8)-(2.11). However, be cause of the proposed importance of the changes in the intracellular calcium concen tration, a fifth equation is required. The rate of change of the intracellular calcium is essentially the difference between influx and efflux of calcium. The influx of calcium is due to the voltage-gated Ca2 channel (‘Ca). Since little information about the metabolism of calcium was known at the time, the efflux of calcium was modelled as —kcaCaj, where ka represents the glucose-dependent rate constant for the removal of Ca:. The higher Chapter 2. Biophysical Models and Preliminary Analysis 19 the concentration of glucose, the higher the value of ka. The equation for Ca thus is a — 1 dT = 47rr3FaO”) — kcaCai] , (2.12) where F is Faraday’s constant and f is the ratio of free to bound intracellular calcium ions. The inclusion of f reflects the fact that the majority of calcium ions entering the cell are immediately absorbed by high affinity binding sites and that the bound calcium ions do not take part in the processes leading to the electrical bursting activity. It is worth noting at this point that, strictly speaking, Ca, therefore represents the intracellular concentration of free calcium ions, rather than the total intracellular calcium concentration. Finally, the factor of accounts for the cell volume and converts ionic4Kr3 charge to concentration. Equations (2.8)-(2.12) define the five-dimensional model for the bursting electrical activity as presented in [18]. Chay and Keizer present a reduced model in [19], where the activation and inactivation variables of the voltage-gated Ca2+ channel, m and h in (2.2), are replaced by their respective steady-state values (m(V) and h(V)) and the relaxation equations for these variables, (2.9) and (2.10), are omitted. With some minor changes to the parameter values, the reduced model produces essentially the same results. Therefore, for the remainder of this thesis, we use the reduced Chay-Keizer model, instead of the original five-dimensional one. The definitions of the steady-state functions m(V), hco(V), and n(V), and the function r(V) completes the description of the reduced Chay-Keizer model. They are 2 13(V)+(V)’ for y = m, h, and n, where m m cm(V) = B (V V) , (2.14) exp [(Vm — V)] — 1 Chapter 2. Biophysical Models and Preliminary Analysis 20 1 /3m(V) = 4.Oexp [-(Va — v)], (2.15) ri crh(V) 0.O7exp — V) , (2.16) h(V) 1__Sh , (2.17) exp —(Vb—V)]+1 a(V) = - V) , (2.18) exp [_(v — V)] —1 /3(V) = 0.l25exp [_(v — V)] , (2.19) and r(V) = a(V) + (V) (2.20) The graphs of the functions m, h, n, and r,, are shown in Figure 2.1. We observe that the steady-state functions for the activation variables m and n increase with V, whereas the steady-state function for the inactivation variable h decreases with V. The values of all parameters present in the model can be found in Appendix A. The reduced Chay-Keizer model can be integrated using a standard computer code for stiff differential equations (e.g., LSODE from ODEPACK [45]). Figure 2.2 illustrates the numerical solution computed using the parameter values listed in Table A.1 in Appendix A. It is clear that the solution of the mathematical model mimics the bursting electrical activity observed experimentally (cf. Figure 1.1). The model also is compatible with the biophysical explanation of bursting proposed by Atwater et al. [61. In particular, when küa is low (representing a low concentration of external glucose), there is little removal of Ca. Thus, Ca remains relatively high and the membrane potential remains near VK, the Nernst potential for K+, due to the dominant ‘K,Ca For higher values of ka, Ca is gradually reduced, ‘K,Ca becomes less dominant and the membrane potential increases. When ‘K,Ca is negligible, ‘K and ‘Ca are activated, Chapter 2. Biophysical Models and Preliminary Analysis 21 (a) 0.7 0.6 _—. 0.5 0.4 8 0.3 0.2 0.1 0 -70 (c) i i i i (b) —L--” 06 - 8 0.4 - 0. -60 -50 -40 -30 -20 -70 -60 -50 -40 -30 -20 V (mV) V (mV) 0.5 0.4 - 17 0.3 (d) 18 S— 16 0.2 - 150.1 - 0 14 -70 -60 -50 -40 -30 -20 -70 -60 -50 -40 -30 -20 V (mV) V (mV) 8 Figure 2.1: Graphs of the functions m(V), h(V), n(V), and r(V), respectively, describing the channel kinetics of the Chay-Keizer model. resulting in the spiking activity. During the spiking, ‘Ca causes Ca to increase. Thus, ‘K,Ca is slowly activated again, eventually reducing the membrane potential back to VK. When this occurs, -1Ca and ‘K have been inactivated. That is, Ca no longer increases, and the cycle starts again. 2.1.3 Rinzel’s Fast/Slow Subsystem Analysis Even though the Chay-Keizer models put the biophysical model by Atwater et al. [6] on a precise mathematical footing, insights as to why the system of equations should give rise to the bursting phenomenon were still lacking. Rinzel [76] provided these insights with his lucid analysis in 1985. Chapter 2. Biophysical Models and Preliminary Analysis 22 (b) 1.2 1.1 1 0.9 0.8 0.7 ) 0.6 0.5 0.4 0.3 0 10000 20000 30000 40000 0 10000 20000 30000 40000 T (msec) Figure 2.2: Numerical solution of the reduced Chay-Keizer model. (a) Membrane poten tial V. (b) Intracellular calcium concentration Ca:. (a) -20 -30 -40 -50 -60 T (msec) Chapter 2. Biophysical Models and Preliminary Analysis 23 The key idea in Rinzel’s treatment of the system is to identify a fast and a slow subsystem. Due to the smallness of the parameter f in (2.12), the time course of Ca2 is much slower than that of the remaining variables, V and n. In fact, the dynamics of Ca2 determine the time scale of the slow overall periodic behaviour (on the order of seconds). For this reason, the third equation of the reduced Chay-Keizer model, (2.12), is defined to be the slow subsystem. The remaining two equations, (2.8) and (2.11), define the fast subsystem. This subsystem is responsible for generating the spikes during the active phase, which occur on a much faster time scale. The slow behaviour of Ca2 is exploited by performing a bifurcation analysis of the fast subsystem, treating Ca as the bifurcation parameter. Once the bifurcation structure of the fast subsystem is understood, the dynamics of Ca2 are included in order to account for the overall bursting behaviour. The first step in the bifurcation analysis is to determine the steady states. Setting ii = n(V) in (2.8) and using (2.4)-(2.7), we thus require 133(V; Ca2) = cam(V)h(V)(V Vga) +Kfl0(V)(V — VK) + gK,Ca Ca (V — VK) + gL(V — VL) = 0. (2.21)Kd+Ca Note that we denote Ca2 as a parameter in I by use of the semicolon. Figure 2.3a shows the behaviour of I as a function of V for several biophysical values of Ca. For low values of Ca, there is one steady state (I = 0) at a high value of V. Similarly, for high values of Ca1, there is one steady state at a low value of V. For intermediate values of Ca1, there are three steady states. Figure 2.3b shows the location of the steady states plotted in (V, Ca2) space. The stability of these steady states can be determined numerically. Rinzel used AUTO [281 to obtain the bifurcation diagram and the results are reproduced in Figure 2.4 for illustrative purposes. The low V steady states, or the steady states on the left branch of the I = 0 curve, are stable nodes (indicated by a solid curve). The steady states on Chapter 2. Biophysical Models and Preliminary Analysis 24 (a) 4 — 3 - - IL -3 - - 4 I I I I -70 -60 -50 -40 -30 -20 V (mV) (b) 3 I = 0 2 - — 1 - c) 0 - -1 I I I -70 -60 -50 -40 -30 -20 V (mV) Figure 2.3: 183(V; Car) and the curve of steady states of the fast subsystem for the reduced Chay-Keizer model. (a) 133(V;Ca1)/4rr2 as a function of V for Ca2 = 0, 1, 2, and 3 M. The arrow indicates the direction of increasing Ca2. (b) Curve of steady states of the fast subsystem. Chapter 2. Biophysical Models and Preliminary Analysis 25 0 •1 -70 Figure 2.4: Bifurcation diagram, generated with AUTO, of the fast subsystem of the reduced Chay-Keizer model. Stable steady states are indicated by a solid curve; unstable steady states by dashes. Saddle-node bifurcations (SN) occur at the local extrema of the = 0 curve. Immediately above the Hopf bifurcation (RB) point, the unstable steady states are surrounded by stable periodic orbits. The maximum and minimum values of the periodic orbits are indicated by the solid curve surrounding the right branch of the I = 0 curve. The branch of periodic orbits terminates at a homoclinic bifurcation (HC). the middle branch are saddle points and, hence, unstable (dashed curve). The steady states on the right branch are spirals, except near the local maximum of the I = 0 curve, where the steady states are unstable nodes. Saddle-node bifurcations, at which nodes become saddles and vice versa, occur at the local extrema of the I = 0 curve and are labelled by SN. For high values of Ca, the spirals on the left branch are unstable (dashed curve); they become stable at a Hopf bifurcation, labelled by RB. The Hopf bifurcation is supercritical and, hence, stable periodic orbits surrounding the unstable spirals emanate from the RB point. The maxima and minima of the orbits are plotted as a solid curve surrounding the right branch. This branch of periodic orbits terminates at 3 2 1 -60 -50 -40 -30 -20 V(mV) Chapter 2. Biophysical Models and Preliminary Analysis 26 a homoclinic bifurcation, labelled by HC. At this point, the periodic orbit has become a homoclinic orbit, with infinite period, connecting the saddle point on the middle branch of the I = 0 curve to itself. Thus, for values of Ca between the local minimum of the I = 0 curve and the homoclinic bifurcation, the fast subsystem exhibits bistability. The low V steady state is stable, whereas the high V steady state is unstable, but surrounded by a stable periodic orbit. This bistability is the key to explaining the bursting behaviour. Recall that during the bursting activity, the membrane potential switches between a low voltage state (the silent phase) and a high voltage oscillation (the active phase). The mechanism for switching between the silent and active phases becomes clear when the dynamics of Ca1 are included in the analysis. Figure 2.5 shows the bifurcation diagram of Figure 2.4 together with the nullcline (Ca1 = 0) for Ca and the projection of the numerical solution from Figure 2.2 superimposed. When the system is in the low V steady state to the left of the Ga nullcline, where Ga < 0, Ca decreases and V increases. When the local minimum of the I = 0 curve is reached, the solution trajectory switches to the oscillation about the high V steady state to the right of the Ca nullcline. Here Ca3 is positive, so that Ca3 increases until the homoclinic bifurcation is reached. At this point, the trajectory switches back to the low V steady state. It is worth pointing out that for the full system of equations, there is one unstable steady state, occurring at the intersection of I = 0 and the Oa3 nullcline, and a stable limit cycle, namely the periodic bursting solution. 2.2 Further Developments In this section, we introduce the remaining first-generation models for bursting electrical activity which consist of three ordinary differential equations. These are the models due Chapter 2. Biophysical Models and Preliminary Analysis 0 V (mV) 27 Figure 2.5: As Figure 2.4, with the solution of Figure 2.2 superimposed. the right. Ca nulicline and the projection of the numerical Ca < 0 to the left of the Ca nullcline; Ca, > 0 to to Chay [13, 14], Chay and Cook [16], Chay and Kang [17], Himmel and Chay [44], and Sherman, Rinzel, and Keizer [87]. Table 2.1 outlines the different currents included in each particular model and com pares some salient features of the models, namely the form of the various currents in cluded, whether a leakage current is present or not and whether the Nernst potential for the Ca2 ion depends explicitly on its concentration difference or not (cf. equation (1.1)). As will become apparent in Section 2.3 and Chapter 3, the most important distinction between the models from an analytical point of view is the value of the integer exponent of the activation variable n in the voltage-gated K channel (cf. Table 2.1). We note that the first three models, namely the reduced Chay-Keizer, Chay (1986), and Himmel-Chay models, have values of the exponent greater than one, whereas the remaining models, 3 2 1 1 -70 -60 -50 -40 -30 -20 Chapter 2. Biophysical Models and Preliminary Analysis 28 namely the Chay-Kang, Sherman-Rinzel-Keizer, and Chay-Cook models, have an expo nent equal to one. Table 2.1 and the discussions below reflect the controversies surrounding the currents thought to play an important role in the bursting electrical activity at the time when these models were developed. Some controversies have already been resolved, while others remain under investigation. Rather than focusing on these controversies, the purpose of this section is to briefly introduce the different forms of the equation for the membrane potential V and the hypotheses underlying the different forms. The complete definition of each model can be found in Appendix A. 2.2.1 The Chay (1986) Model During the early 1980’s, researchers began applying the relatively new patch-clamping technique to study the function of single ion channels in pancreatic /3-cells. Two sepa rate K+ channels were identified, namely a low-conductance and a high-conductance K+ channel. Experiments by Ashcroft et al. [3], Cook and Hales [21], and Findlay et al. [32] showed that the low-conductance K channel is directly inhibited by a metabolite of glucose, namely ATP. Ashcroft et al. [3] showed that the current due to this channel dominates in the absence of bursting and is responsible for the resting potential of the membrane in the absence of glucose. Experiments by Cook et al. [22] and Findlay et al. [33] showed that the high-conductance K+ channel can be activated by both membrane depolarization and intracellular Ca2+ ions. Based on these results, it was suggested that the two separate K+ channels as used in the Chay-Keizer models, namely the voltage-gated K+ channel and theCa2+activ ted K+ channel, are in fact the same channel. Chay [13] subsequently modified the reduced Chay-Keizer model and demonstrated Chapter 2. Biophysical Models and Preliminary Analysis 29 Ca2+ channels K+ channels 1 Leakage Model voltage- voltage-gated, voltage- Ca2+ current Vca(Caj) gated Ca2-inhibited gated activated included - Ca1 gcamhoo gKn4 yes no 2 9K(V,Ca)fl2 yes yes — Ca13 camhoo KTi3 yes yes 4 gca(V,ca)moohc. YK yes no - Ca15 gcamc,Dhc,D gK no no 6 gcamco gCa(V,Ca)8 yes no Table 2.1: Comparison of the first-generation models consisting of three ordinary dif ferential equations. The table indicates the forms of the conductances of the Ca2+ and K currents included in the equation for the membrane potential V (cf. (2.8)), whether a leakage current is included or not, and whether the Nernst potential for calcium is allowed to vary with Ca1 or not. The models are 1: Reduced Chay-Keizer [19]; 2: Chay (1986) [13]; 3: Himmel-Chay [44]; 4: Chay-Kang [17]; 5: Sherman-Rinzel-Keizer [87]; 6: Chay-Cook [16]. Chapter 2. Biophysical Models and Preliminary Analysis 30 that a single K+ channel, namely a K+ channel activated by both membrane depolar ization and intracellular Ca2+ ions, along with a voltage-gated Ca2+ channel also can account for the bursting activity observed in /3-cells. Chay omitted the modelling of the low-conductance ATP-inhibited K+ channel since the current due to this channel is negligible during bursting. The modified equation for the membrane potential V (cf. (2.8)) is 4r2Cm = — [Ia(V, Ca) + IK(V,Ca)( V, Ca) + IL(V)], (2.22) where ‘Ca represents the voltage-gated Ca2+ current, IK(V,Ca) represents the current flowing through the voltage-gated, Ca2+sensitive K+ channel, and ‘L represents the leakage current as in the Chay-Keizer models. The ‘Ca and1rK(V,Ca) currents are modelled by ICa(V, Ca) = gCam(V)(V — Vga), (2.23) ‘K(V,Ca)(V, Ca1) = .K(V,Ca)T12O” — VK), (2.24) and the IL(V) current is as defined by (2.7). The function m(V) represents the open probability of the Ca2 channel and n represents the probability of opening the K channel. The reason for raising m to the power of 2, rather than 3 as in the Chay Keizer models, is that it is the lowest power giving rise to spiking. The reason for raising n to the power of 2, instead of 4 as in the Chay-Keizer models, is based on the observations in [22, 33] that the slopes of the probability curves for opening are not very steep. That is, a power of 2 fits the curves better than a power of 4. The inactivation variable h (or heo) for the Ca2 current (cf. (2.2)) has been omitted in light of experimental results by Satin and Cook [82] which showed that inactivation is never complete and is not reproduced by Hodgkin-Huxley type kinetics. In addition, the Nernst potential for Ca2, VCa, is allowed to vary with the intracellular calcium concentration, rather than being Chapter 2. Biophysical Models and Preliminary Analysis 31 kept at a constant value as in the Chay-Keizer models. Finally, the calcium dependence of ‘K(V,Ca) comes in through the steady-state value of ii, n(V, Ca1), and the relaxation function r(V, Cat). The complete model can be found in Appendix A. 2.2.2 The Himmel-Chay Model In 1987, Himmel and Chay [44] developed a model to test the hypothesis that the cur rent through the low-conductance ATP-inhibited K channel [3, 21, 32] dominates when the glucose concentration is low, while the current through the high-conductance K+ channel(s) dominates the silent phase during bursting. The model is very similar to the reduced Chay-Keizer model, with the following excep tions. First, Himmel and Chay add a term which represents the current flowing through the low-conductance ATP-inhibited K channel to the equation for the membrane po tential V (cf. (2.8)). Since it is assumed that the conductance of this channel is negligible when the glucose concentration is high enough so that the membrane potential exhibits the bursting phenomenon, we omit this extra current. Second, the activation variable n is raised to the power of 3 instead of 4 and the value of the maximum conductance g< is decreased. Himmel and Chay conclude that since there are no appreciable differences in the results, this aspect of the model must not be very important. Finally, the Nernst potential Vca is allowed to vary with the intracellular calcium concentration, as for the Chay (1986) model. The result of this change is that the amplitude of Ca1 is smaller and in better agreement with experiment. For a precise definition of the equations and the values of all parameters, the reader is referred to Appendix A. Chapter 2. Biophysical Models and Preliminary Analysis 32 2.2.3 The Chay (1987) and Chay-Kang Models Up until now, the mathematical models and, in particular, the functions related to the activation and inactivation variables, were all based on the Hodgkin and Huxley model for the squid giant axon [47]. Since little was known about the properties of the ion channels in ,8-cells, the functions describing the channel kinetics (m, h, etc.) were taken directly from the Hodgkin and Huxley model and shifted along the voltage axis until the new model exhibited the bursting phenomenon. In 1986, Rorsman and Trube [80] published voltage-clamp data that allowed mathe matical modellers to represent the properties of the voltage-gated K+ and Ca2+ channels more realistically. Rorsman and Trube found that the outward current appears to be carried mainly through a K+ channel that is insensitive to Ca2+ ions and that the inward current is carried mainly through a voltage-gated Ca2+ channel. This current seems to be inactivated slowly by an inhibitory effect of intracellular Ca2 itself (the slow inactivation also was observed by Satin and Cook [82]). The first mathematical models incorporating the data from Rorsman and Trube were those by Chay [14] and Chay and Kang [17]. The two models are practically identical and we consider only the Chay-Kang model. The equation for the membrane potential Vis 4r2Cm = — {ICa(V,Ca)(V, Ca1) + IK(V) + IL(V)j, (2.25) where ICa(V,Ca) represents the voltage-gated Ca2+ current that is slowly inactivated by Ca2+, ‘K represents the voltage-gated K+ current, and ‘L represents the usual leakage cur rent, defined by (2.7). Chay and Kang omitted the term representing theCa2-activ ted K+ current (cf. (2.6)) based on the observation thatCa2+activ ted K+ channels are open very rarely under physiological voltages and intracellular calcium concentrations [22, 33]. The calcium feedback mechanism is replaced by the Ca2+inactiv tion of the calcium Chapter 2. Biophysical Models and Preliminary Analysis 33 channels. The ICa(V,c’a) and ‘K cnrrents are modelled by ICa(V,Ca)(V, Ca) = gGa(vCa)m(V)h(V, CctJ(V — Vca), (2.26) IK(V) = Kri(V — VK), (2.27) where m, h, and n have the usual meaning. Note that all activation and inactivation variables, in particular, the activation variable ii, appear linearly. Furthermore, the voltage-clamp data of Rorsman and Trube are represented by the much simpler definitions of m, h, n, and r, (see Appendix A). 2.2.4 The Sherman-Rinzel-Keizer Model In 1988, Sherman, Rinzel, and Keizer [87] modified the channel kinetics of the voltage gated K+ and Ca2+ channels in the reduced Chay-Keizer model to incorporate the voltage-clamp data of Rorsman and Trube [80]. Therefore, the equation for the mem brane potential V almost is identical to that of the reduced Chay-Keizer model, except that the leakage current has been omitted. The model still is based on the premise that the calcium feedback mechanism is provided by the Ca2+activation of a K+ channel, rather than by theCa2+inactiv tion of the Ca2+ channels as suggested by Rorsman and Trube. We mention that the model was not developed to incorporate all of the observations of Rorsman and Trube. Instead, it was developed to provide a more realistic basis than the reduced Chay-Keizer model for a subsequent extension which includes stochastic Ca2+ activated K+ channels to study the effect of channel sharing and explain the irregular electrical behaviour of single cells and small clusters. We do not consider the stochastic version of the model in this thesis. Chapter 2. Biophysical Models and Preliminary Analysis 34 The equation for the membrane potential V in the deterministic version of the Sherman-Rinzel-Keizer model is 4r2Cm = — [10a(V) + IK(V) + IK,Ca( V, Cad)], (2.28) where ‘Ga and ‘K are defined by Iüa(V) = gcam(V)h(V)(V — Va), (2.29) IK(V) = Kfl(V — VK), (2.30) and ‘K,Ca by (2.6). In comparison with the reduced Chay-Keizer model (cf. (2.8)), we note that the leakage current has been omitted. Furthermore, the activation function m and the activation variable n both appear linearly (cf. (2.4) and (2.5)). Finally, as for the Chay-Kang model, the definitions of m, h, n, and r are much simpler (see Appendix A). 2.2.5 The Chay-Cook Model In 1988, Chay and Cook [16] developed a unified model that can give rise to two types of bursting, namely the so-called square-wave bursting (as is observed in pancreatic 3- cells) and parabolic bursting (as is observed in the Aplysia abdominal ganglion R15 cell [73]). Parabolic bursting differs from square-wave bursting in that the spike frequency increases towards the middle of a burst and then decreases again. The name parabolic refers to the plot of instantaneous spike period versus spike number, which appears as a concave-upward parabola. For the purpose of this thesis, the set of parameters that produces square-wave bursting is used. The model is very similar to the Chay-Kang model, except that the voltage-gated, Ca2+inhibited Ca2+ current has been split into two separate currents. That is, the Chapter 2. Biophysical Models and Preliminary Analysis 35 equation for the membrane potential V is given by 4Kr2Cm = — [Ia(V) + ICa(V,Ca)( V, Cad) + IK(V) + IL(V)j, (2.31) where ‘Ca represents the voltage-gated component of the Ca2+ current, ICa(V,Ca) repre sents theCa2+inhibited component of the Ca2+ current, ‘K represents the voltage-gated K+ channel, and ‘L represents the leakage current as usual. As for the Chay-Kang model, the calcium feedback mechanism is provided by theCa2-inhibited component of the Ca2+ current, ICa(V,Ca). The ‘Ca and ‘Ca(V,Ca) currents are modelled by Ia(V) = .camcy(V)(VVca), (2.32) ICa(V,Ca)(V, Ca1) = gCa(V,Ca)soo(V, Caj(V — Vga). (2.33) For precise functional definitions, see Appendix A. 2.3 Nondimensionalization and Problem Reformulation We first nondimensionalize the equations and scale the variables so that they are 0(1). In order to write all models in a standard form, a consistent nondimensionalization of the model equations is required. The following is the basis for the nondimensionalization of all models: —v Ca1 Tv = ----, w = (7K)n, z = ln—j--, t = VK ‘- ______ —V —Si = 2’’ = , = -n-—, = —BVK,4irr Lm VK VK - 3VKCmfkja Tn, 9 = 11 i — Tm rrlccarna where p is the integer exponent of the activation variable n for the voltage-gated K+ channel, X is a parameter with the dimension of M (there is at least one such param eter in each model, usually the dissociation constant Kd), and is approximately the Chapter 2. Biophysical Models and Preliminary Analysis 36 maximum value of r(V) over the physiological range of V. The subscripts i take on the appropriate symbols, such as K, Ca, L, etc., depending on the particular model being nondimensionalized. Letting v = —V/VK ensures that the nondimensional membrane potential has a value near —1 in the silent phase, since at that time V is nearly the Nernst potential for potassium, VK. Note that v takes on negative values to maintain the same orientation of the nondimensional (mathematical) and dimensional (physical) membrane potential. The nondimensionalization of i’z is such that in attains a maximum of approximately 1 during the active phase. The nondimensionalization of time T ensures that the maximum relaxation time associated with the in variable is approximately equal to 1. For the nondimensionalization of Ca1, it is not sufficient to simply take the ratio of Ca1 with any parameter X having the dimension of pM. For some models, notably the Chay (1986) model, the only choice for X is a parameter several orders of magnitude larger than Ca1 (for example, the extracellular calcium concentration, about 1000 times larger than the intracellular concentration). The form of the differential equation for Ca1 Calleads to the choice z = ln This general nondimensionalization is complete for all models, except the two models which allow the calcium Nernst potential to vary with the intracellular Ca2+ concentra tion, namely the Chay (1986) and the Himmel-Chay models. Two or three additional nondimensional parameters are required for these cases. The reader is referred to Ap pendix A for details. With this general nondimensionalizatiori, all first-generation 3-cell models which are summarized in Table 2.1 can be written in the standard form ica(v,z)—t9(v+1) —g(z)(v+1) + 1(v), (2.34) din — th0(v,z)—w 235dt — f(v,z) = e [fle_zica(v, z) — 11, (2.36) Chapter 2. Biophysical Models and Preliminary Analysis 37 where a(V, z) represents the Ca2+ current(s), w”(v + 1) represents the voltage-gated K current (also activated by Ca2 in the case of the Chay (1986) model), g(z)(v + 1) represents theCa2+activ ted K+ current, and 1(v) represents the leakage current. Note from Table 2.1 that g(z) and l(v) can be zero, depending on the particular model. The functions ü3(v, z) and z) are the nondimensionalized n(V, Cat) and r(V, Ca), respectively. The parameter E is small, approximately iO to 1O, depending on the specific model. The parameter /3 is (inversely) related to the blood glucose level. Finally, p is an integer, ranging from 1 to 4, again, depending on the particular model. Although p looks innocuous, it plays an important role in the analysis, as we will see shortly. Next, we simplify the equation for v by using the transformation u = ln(v + 1). Since v is always > —1 with the nondimensionalization described above, this transformation is well-defined. The general nondimensional model then can be rewritten in the generic form = A(u, z) — wi’, (2.37) = w(n,z) —w (2.38) r(u, z) E [/3e_zh(u, z) — 1], (2.39) where A(u, z) = e_u[ica(v(u), z) + l(v(u))] — g(z), (2.40) w(u,z) = th(v(u),z), (2.41) r(u,z) = (v(u),z), (2.42) h(u,z) = ica(v(u),z). (2.43) Finally, we combine the first two differential equations, (2.37) and (2.38), into a single second-order equation by eliminating vi. This is achieved by differentiating (2.37) with Chapter 2. Biophysical Models and Preliminary Analysis 38 respect to t and using equations (2.37)-(2.39). The resulting system is ii + E(u, it, z) + F(u, z)it + G(u, z) — H(u, z), (2.44) = EK(u,z), (2.45) where E(u, it, z) = pw(u, z) [A(u, z) — itJ, (2.46) r(u, z) F(u,z) = [ p — 8A(u,z)] , (2.47)r(u,z) G(u,z) = (2.48) H(u,z) = öA(uz)K(UZ) (2.49) K(u,z) = /3e_zh(u,z) —1. (2.50) We note that for those models for which p = 1, the term involving it in (2.46) reduces to unity, so that E(u, it, z) no longer depends on it. For these models then, the only damping term in (2.44) is the F(u, z)it term, which makes the analysis in Chapter 3 slightly easier than for the remaining models. More generally, since is a small parameter, z changes on a much slower time scale than u and it. Hence, the slow subsystem consists of (2.45). The fast subsystem is given by (2.44), with e = 0 and z treated as a parameter, i.e., it + E(u, ii; z) + F(u; z)it + G(u; z) = 0. (2.51) From (2.51), we deduce that the fast subsystem nullcline is given by (u; z) = 0, (2.52) where we define ã(u; z) = G(zt; z) + E(u, 0; z). (2.53) Chapter 2. Biophysical Models and Preliminary Analysis 39 Using the parameter values listed in Appendix A, the bifurcation diagrams of (2.51) have the same qualitative features for all of the models under consideration, as shown in Figure 2.6. The fast subsystem nullcline G = 0 has a cubic-like shape in all cases. Steady states on the left branch are always stable nodes; steady states on the middle branch are saddle points; steady states on the right branch are either spiral sources or spiral sinks, except near the local maximum of the G = 0 nullcline, where the steady states are nodes. In Figure 2.6a, the stability of the spirals changes only once, at a supercritical Ropf bifurcation lying below the local minimum of the G 0 nullcline. The branch of periodic orbits emanating from the Hopf point terminates at a homoclinic bifurcation on the middle branch. In Figure 2.6b, there are two supercritical Hopf bifurcations. The lower Hopf bifurcation (RB 1) lies below the local minimum of the G = 0 nullcline, whereas the upper Hopf bifurcation (HB2) lies in the region of bistability, near the local maximum of the G = 0 nullcline. The branches of periodic orbits emanating from both Hopf points terminate at homoclinic bifurcations (Rd and HC2, respectively) on the middle branch. For models of this type, the upper Hopf bifurcation and corresponding homoclinic bifurcation do not play a role in the bursting behaviour, since z never attains a high enough value to be affected by these bifurcations. The appearance of a second Hopf bifurcation on the right branch is the result of a codimension-2 Takens-Bogdanov bifurcation, which we investigate in more detail in Chapter 5. In Chapter 3, we investigate methods to determine the location of the homoclinic bifurcations for all models. In Chapter 4, we use the approximations from Chapter 3 as a starting point for computing approximate plateau fractions. Chapter 2. Biophysical Models and Preliminary Analysis 40 (a) SN N SN HB*\ U (b) \ G=0 HC2..B2 N HB1\ U Figure 2.6: Illustration of the two types of bifurcation diagrams exhibited by the fast subsystem, (2.51), of the nondimensional /3-cell models. (a) Bifurcation diagram for the reduced Chay-Keizer and Himmel-Chay models. (b) Bifurcation diagram for the Chay (1986), Chay-Kang, Sherman-Rinzel-Keizer, and Chay- Cook models. Chapter 3 Escape from the Active Phase In this chapter, we investigate various methods for approximating the value of z = Zesc at which the solution trajectory of the full system of equations, (2.44)-(2.45), escapes from the active phase and returns to the silent phase. The reason for doing so is as follows. In Chapter 4, we derive leading-order approximations of the silent and active phase durations and, subsequently, a leading-order approximation of the plateau fraction. We thus need to know when the silent phase starts and ends and, similarly, when the active phase starts and ends. Figure 3.1 shows the projection of the numerical solution to the full system of equa tions onto the bifurcation diagram of the fast subsystem, (2.51). We introduce (tim, Zm) and (tiM, zM) to denote the local minimum and maximum of the G = 0 nullcline, re spectively. We observe that the solution trajectory undergoes the transition from the silent to the active phase near (tim, Zm). This point is easy to calculate. Therefore, we know approximately when the silent phase ends and the active phase begins. However, determining a priori when the solution undergoes the transition from the active phase back to the silent phase, i.e., when the active phase ends and the silent phase begins, is more difficult and the subject of this chapter. We define Zesc to be the value of z at which the solution trajectory projected onto (ti, z) space intersects the middle branch of the G = 0 nullcline (cf. Figure 3.1). Evaluations of Zesc indicate that it depends on the value of the glucose-dependent parameter 3. Since we are interested in the plateau fraction as a function of j3 in Chapter 4, ideally, we also 41 Chapter 3. Escape from the Active Phase 42 N Figure 3.1: Typical bifurcation diagram of the fast subsystem, (2.51), with the numerical solution of the full system of equations, (2.44)-(2.45), superimposed. The projection of the numerical solution intersects the middle branch of the C = 0 nullcline at (uesc, Zesc). The branch of periodic orbits terminates at z ZHC and the corresponding homoclinic orbit connects the saddle point at (UHC, zHC) back to itself. The local minimum and maximum of the C = 0 nullcline occur at (urn, Zrn) and (UM, zM), respectively. wish to determine the dependence of Zesc on To obtain an approximation of Zese, we consider the fast subsystem. From Figure 3.1, we observe that the solution trajectory of the full system of equations undergoes the transition from the active phase back to the silent phase near the homoclinic bifurcation of the fast subsystem. Therefore, letting ZHC be the value of z at which the homoclinic bifurcation of interest occurs, namely HC in Figure 2.6a or HC1 in Figure 2.6b, we may use Z€3 ZHC. Since the glucose-dependent parameter /3 only appears in the slow subsystem, none of the methods described below can be used to determine the dependence of Zesc on /3. However, the results in Chapter 4 are not affected significantly by using the constant approximation of Zesc obtained in this chapter. it Chapter 3. Escape from the Active Phase 43 In Section 3.1, we describe how the bifurcation analysis of (2.51) with the AUTO software package can be used to approximate ZHC. In Section 3.2, we use bisection approaches which are based on the bifurcation of the solution behaviour of (2.51) in the (u, ii) phase plane as z is varied. In addition to these numerical methods, we discuss two analytical methods. In Section 3.3, we compute the leading-order approximation of zHc using the Melnikov [64] integral. However, this method yields sufficiently accurate results for only two of the six models under consideration in this thesis. To obtain higher- order corrections to the Melnikov results, we use a version of the Fredholm alternative in Section 3.4. Finally, a comparison of the various approaches and a discussion of their advantages and disadvantages is provided in Section 3.5. 3.1 Approximation of z3 with the AUTO Software Package In this section, we use the AUTO [28] software package to approximate ZHC and, subse quently, Zesc. When AUTO computes branches of periodic orbits emanating from a Hopf point (cf. Figure 2.6), it also determines the period of the respective periodic orbits. The period increases as the bifurcation parameter z increases, until it is infinite at the homoclinic bifurcation. Although numerical difficulties at values of z near the homoclinic bifurcation prevent the exact determination of ZHC, we can obtain a very accurate ap proximation of ZHC, ZAUT0, by continuing the computation of a branch of periodic orbits until the period is sufficiently large. For the first-generation models under consideration, approximations were sought by requiring that the period at z ZAUTO was at least 1000. The results of this approach are summarized in Table 3.1. Chapter 3. Escape from the Active Phase 44 Model Zm Zesc ZAUTO ZBIS ZMEL Z0 Z0 + 6Z1 1 -0.85967 -0.030 1.00247 1.00180 0.85174 0.85174 0.99065 2 -7.94215 -7.9324 JF 1.03691 1.03691 0.43092 0.43194 1.15586 3 -1.87072 -1.368 1.04949 1.04951 0.92087 0.92087 1.08760 4 -5.48746 -5.062 1.02649 1.02621 0.99297 0.99293 0.99293 5 -5.22656 -5.098 1.04846 1.04753 1.02683 1.02676 1.02676 6 -1.25558 -0.596 1.00027 1.00020 0.90639 0.90638 0.90638 Table 3.1: Approximate values of ZHC and, subsequently, Zest obtained with AUTO (ZAUTO), bisection (zBIs), Melnikov’s method (zMEL), and the Fredhoim alternative method (zo and z0 + Szi), for all first-generation models under consideration in this thesis. To facilitate comparison of the methods, the values have been normalized with Zm corresponding to 0 and Zesc corresponding to 1. As in Figure 3.1, Zm is the value of z at the local minimum of the G = 0 nullcline, and e3c is the approximate median of exact values of Zesc obtained for values of /3 near 13HC (cf. Figure 4.1 in Chapter 4). The models are 1: Reduced Chay-Keizer; 2: Chay (1986); 3: Himmel-Chay; 4: Chay-Kang; 5: Sherman-Rinzel-Keizer; 6: Chay- Cook. 3.2 Approximation of Zesc with Bisection Techniques In this section, we focus on the dependence of the solution trajectories of (2.51) on the value of the parameter z in the (u, ii) phase plane. We discuss two bisection techniques based on the bifurcation of the solution behaviour at z ZHC. Even though both methods are less efficient and more cumbersome than the use of AUTO, we include a discussion of them, as this provides valuable insights for the techniques described later. We note that the (u, ‘it) phase plane is perpendicular to the (u, z) plane shown in Figure 3.1. For each value of z between Zm and ZM, there are three fixed points, corre sponding to the three branches of the C = 0 nullcline. The left fixed point is a stable node, the middle fixed point is a saddle, and the right point is, in general, a spiral. We are interested in the behaviour of solution trajectories of (2.51) in relation to these fixed points for different values of z. Chapter 3. Escape from the Active Phase 45 Figure 3.2 shows trajectories obtained numerically for a value of z z1 < ZHC and a value of z = 22 > ZHC. It is evident that trajectories starting near the right fixed point approach a stable periodic orbit when z < ZHC, whereas trajectories with the same initial conditions eventually approach the left fixed point when z> ZHC. We now use a standard bisection procedure, starting with 21 and z2, to zero in on the exact value of zHC. The resulting approximations of ZHC and, subsequently, z are referred to as ZBIS, and are summarized in Table 3.1. As expected, ZBIS is in close agreement with ZAUTO. Alternatively, we can examine the relative position of the stable and unstable orbits of the saddle point, or the middle fixed point, in the (u, it) phase plane of (2.51). Stable and unstable orbits can be computed with the software package DSTOOL [39]. The results corresponding to Figures 3.2a and b are shown in Figure 3.3. It is clear that for the case z = z1 < zHC, the right stable orbit lies on the outside of the right unstable orbit, whereas for the case z = z2 > zHC, the situation is reversed. When z = ZHC, the two orbits join, creating the homoclinic orbit which connects the saddle point to itself. Starting with z1 and z2, the same standard bisection approach as above then can be used in conjunction with DSTOOL to approximate ZHC and, subsequently, 2esc• The applicability of this method was verified for the Sherman-Rinzel-Keizer model and the resulting approximation was within 1O of zBIS. Since the resulting approximations for the remaining models are expected to be essentially identical to 2BIS as well, this method has not been implemented for the remaining models. Chapter 3. Escape from the Active Phase 46 (a) (b) Figure 3.2: Solution trajectories of the fast subsystem, (2.51), for (a) z = z1 < zHC and (b) z = z2 > zHC. The bullets indicate the location of the three fixed points of (2.51), corresponding to the three branches of the G = 0 nullcline. For both cases, the left fixed point is a stable node, the middle fixed point is a saddle, and the right fixed point is a spiral source. Chapter 3. Escape from the Active Phase 47 (a) (b) Figure 3.3: Phase portraits of the fast subsystem, (2.51), calculated with DSTOOL, showing the stable (solid curves) and unstable orbits (dashed curves) of the saddle point. (a) The right stable orbit lies outside the right unstable orbit for z = z1 < ZHC. (b) The right stable orbit lies inside the right unstable orbit for z = z2 > ZHC. U ‘U Chapter 3. Escape from the Active Phase 48 3.3 Approximation of z with Melnikov’s Method In this section, we use an analytical method due to Melnikov [64] to determine the leading- order approximation of ZHC. The method is based on the perturbation analysis of a near integrable system. The unperturbed system is assumed to be integrable and to possess a separatrix solution, smoothly joining the stable and unstable orbits for a hyperbolic saddle point. When the system is perturbed, the separatrix breaks and, in general, the stable and unstable orbits no longer smoothly join. For autonomous perturbations, the stable orbit either lies completely inside the unstable orbit or vice versa. In [64], Melnikov derived an integral to calculate the distance between the stable and the unstable orbits. The sign of this integral determines the relative orientation of the orbits. In our case, the relative orientation depends on the bifurcation parameter z (cf. Figure 3.3). The relative orientation changes, to leading order, at z ZMEL, corresponding to the homoclinic orbit of the fast subsystem at z = zHC. In order to apply Melnikov’s method, it is necessary to rewrite the fast subsystem equation (2.51) as a perturbed integrable system. From (2.46), we see that the factor in volving it in the E(u, it; z) term is unity for models with p = 1. That is, only the F(u; z)ii term is responsible for any damping in the system. For the Sherman-Rinzel-Keizer model, Pernarowski [67] and Pernarowski et al. [69, 70] showed that this damping term is nu merically small relative to the remaining terms, so that (2.51) is in perturbed Lienard form. Writing (2.51) in precisely this form allowed Pernarowski [67] and Pernarowski et al. [70] to apply Melnikov’s method and obtain a good approximation of the location of the homoclinic bifurcation. In general, it is not necessary that (2.51) is written in perturbed Liénard form. It is sufficient to write (2.51) as a perturbed Hamiltonian system, which we do for all models in Section 3.3.1. Having obtained a near-integrable system for all models, we then apply Chapter 3. Escape from the Active Phase 49 Melnikov’s method in Section 3.3.2. 3.3.1 Reformulation of the Fast Subsystem as a Near-Integrable System In this section, we show that the fast subsystem, (2.51), of all models can be rewritten as a perturbed Hamiltonian system. We initially follow the argument used by Pernarowski [67] and Pernarowski et al. [70] for the Sherman-Rinzel-Keizer model to show that the damp ing term is numerically small relative to the remaining terms. Attempting to extend the same kind of reasoning to the models with p 1 reveals a minor flaw in the argument presented in [67]. The same flaw is not presented in [70]. However, the discussion there may be misleading. In this section then, we also present a revised argument. We begin with the models for which p = 1, namely the Chay-Kang, Sherman-Riuzel Keizer, and Chay-Cook models. As mentioned earlier, the E(u, IL; z) term does not depend on IL for these models, so the F(u; z)IL term is the only dissipative term. For this special case, we write (2.51) as ii + F(u; z)IL + G(u; z) = 0, (3.1) where (u;z) = 1 [w(u;z) — A(u;z)] (3.2) r(u, z) from (2.46) and (2.48) and F(u; z) is defined by (2.47). Figure 3.4 shows the graphs of j c(u(s), z(s))ds versus IL from the solution of the full system of equations, (2.44)-(2.45), for several oscillations in the active phase. We observe that the curves deviate little from a straight line with slope —1. In fact, this behaviour persists over the entire active phase, with the small drift of the ioops being due to the slowly varying oscillations as z changes. We note that for the computations shown in Figure 3.4, t = 0 occurs in the middle of a silent phase and, thus, reflect the accumulation due to the (u, z) term during the remainder of the silent phase as well as Chapter 3. Escape from the Active Phase -0.05 0 0.05 0.1 n -0.05 0 0.05 0.1 n Figure 3.4: Comparison of j (u, z)dt versus ii over several oscillations in the ac tive phase, for (a) Chay-Kang model; (b) Sherman-Rinzel-Keizer model; (c) Chay-Cook model. In each case, the curve deviates little from the straight line with slope —1 through out the entire active phase. 50 (a) -0.3 -0.36 - I -0.42 - -0.48 -0.1 I I I I -0.2 (b) 0.4 0.3 - o 0Io.: (c) -0.24 -0.28 - N -0.32 - -0.36 - -0.4 - -0.44 -0.1 0.2 I I I I I Chapter 3. Escape from the Active Phase 51 during the first portion of the active phase. That is, we may vertically shift all curves to the origin, and the active phase then can be approximated by z + jt O(u(s); z(s))ds 0, (3.3) implying that the cumulative effect of the damping term is small in the active phase. In [67], Pernarowski differentiates (3.3) with respect to t, giving ii + G(u; z) 0, and concludes that the local effect of the F(u; z)iz term in (3.1) also must be small. Below we show that this argument is false. The subtleties of the argument are circumvented in a subsequent paper by Pernarowski et al. [70], where they mention that the observation of the active phase approximation in (3.3) motivated them to compare the (local) magnitude of the damping term with the remaining terms for the Sherman-Rinzel-Keizer model. This comparison showed that the damping term is indeed small. Motivated similarly, we find that F(u; z)ñ is numerically small relative to the remaining terms in (3.1) for all models with p = 1, 50 that (3.1) is in perturbed Liénard form for these models as well. We now examine (2.51) for the models for which p 1, namely the reduced Chay Keizer, Chay (1986), and Himmel-Chay models. For these models, both the E(u, ü; z) and F(u; z)ñ terms contribute to the dissipation in the system. Based on the definition of O(u; z) in (2.53), we hypothesize that the E(u, t; z) term must be split into separate components. Writing the Taylor series expansion of E(u, ii; z) about ii 0, we obtain E(u, ü; z) = E(u, 0; z) + (u, 0; z)ñ + .... (3.4) Numerical experiments reveal that the first two terms in (3.4) are the most important ones. Replacing E(u, it; z) in (2.44)-(2.45) by the first term of the Taylor series does not give a bursting solution. Bursting is recovered when the second term of (3.4) is included as well, and the solution already is in good agreement with the solution of the full system of equations, as shown in Figure 3.5. Adding subsequent terms of (3.4) merely improves Chapter 3. Escape from the Active Phase 52 -0.2 -0.4 -0.6 -0.8 -1 0 2000 t Figure 3.5: Comparison of the solutions u(t) of the full system of equations (solid curve), (2.44)-(2.45), and the same system of equations where E(u, ii, z) has been replaced by the first two terms of its Taylor series about zi 0 (dashed curve). the agreement. We group E(u, 0; z) with G(u; z) and (u, 0; z)z with F(u; z)ü in (2.51). The question remains where the other terms in (3.4) belong. We investigate two possibilities. The first is to include the remaining terms in (3.4) with G(u; z), i.e., to write (2.51) as ii + [F(u; z) + (u, 0; z)] z + O(u, t; z) = 0, (3.5) where O(u, ; z) = G(u; z) + E(u, z; z) — (u, 0; z)z (3.6) includes all terms of the Taylor series except the main damping term (u, 0; z)it. Fig ure 3.6 shows the plots of j O(u(s), a(s), z(s))ds versus . As for the models with p = 1, the relationship is approximately linear with slope —1, so that we can infer that the cumulative effect of the middle term in (3.5) is small in the active phase. However, 500 1000 1500 Chapter 3. Escape from the Active Phase N N (a) -0.5 53 -0.59 -0.68 -0.77 -0.86 -0.95 -0.2 I I I I I I (b) -0.08 -0.1 0 0.1 0.2 n -0.16 - -0.24 - -0.32 -0.12 I I I I I I -0.06 0 0.06 0.12 t (c) -0.24 -0.29 - - -0.34 - - L -0.39 - - -0.44 I I I -0.1 -0.05 Figure 3.6: Comparison of j O(u, z)dt versus over several oscillations in the active phase, for the (a) reduced Chay-Keizer model; (b) Chay (1986) model; (c) Himmel-Chay model. As in Figure 3.4, the curve lies along the straight line with slope —1 throughout the entire active phase. 0 U 0.05 0.1 Chapter 3. Escape from the Active Phase 54 the examination of the local magnitude of this term relative to ü and G reveals that, in fact, it is 0(1). That is, (3.5) is not in perturbed form and, therefore, of no use for the application of Melnikov’s method. The reason why this 0(1) term essentially has no cu mulative effect in the active phase is that it is roughly symmetric about zero during each oscillation. That is, any positive cumulative effect during the first half of the oscillation is cancelled by the negative cumulative effect during the second half of the oscillation. The second possibility is to include the remainillg terms in (3.4) with F(u; z)it, i.e., to write (2.51) as ü + [F(u; z)i.i + E(u, t; z) — E(u, 0; z)] + G(u; z) = 0, (3.7) where is defined by (2.53). Figure 3.7 shows the plots of j (u(s), z(s))ds versus We observe that during each oscillation, the relationship again is approximately linear with slope —1, but that there is a small drift. The small drift remains when z is held constant. That is, there is a cumulative effect due to the damping term in (3.7). However, an examination of the local magnitude of this term reveals that it is small relative to ii and (u; z) in the active phase and that (3.7), in fact, is written in perturbed form. The cumulative effect of the damping term is due to the fact that, even though it is small, it is not symmetric about zero and there is no cancellation as in the previous case. We note that (3.7) includes all cases, p = 1,2,3,4, and as such, we use this form of the fast subsystem equation for the derivation of the Melnikov integral. In light of the smallness of the damping term, we infer the existence of an artificial small parameter, 6, so that the damping term is 0(6) and rewrite (3.7) as ii + 6E(u, ‘a; z) + (u; z) = 0, (3.8) where fr(u, ‘a; z) = [F(u; z)’a + E(u, ‘a; z) — E(u, 0; z)] (3.9) z)d t O 1 • C. 31 Co c O (u, z)d t I I I I I I p jO (u ,z )d t p I. c d • CD CD CD CD I- CD — .C D ‘ Co ‘ fl CD ( Co C p CD < -- ‘ C o n C C D C D c, -• N —.CD . _ _ Co C o O C D Co > S : - : - ’ Co o Co — . - ‘C D CD CD ‘ C C CD Co Co Co C — C D ; C D ÷ C D 0 ’ 0i U 0 — . — , P C D I CD 1• © p I. p I. I I I I I I p I I I (.5 1 (.5 1 Chapter 3. Escape from the Active Phase 56 r Model 6 Reduced Chay-Keizer 0.00126 0.060 Chay (1986) 0.000024 0.025 Himmel-Chay 0.00108 0.030 Chay-Kang 0.00054 0.012 Sherman-Rinzel-Keizer 0.00066 0.080 Chay-Cook 0.00072 0.015 Table 3.2: Comparison of e and 6, showing that 6>> 6> 0. and is defined by (2.53). With this formulation, the fast subsystem nulicline still is given by O 0. Furthermore, (3.8) can be referred to as a strongly nonlinear oscillator, i.e., as 6 —* 0, (3.8) remains nonlinear. In fact, when 6 = 0, (3.8) describes a Hamiltonian system (cf. Section 3.3.2). However, we keep 6 fixed in this thesis. We will see below that we do not need to know the value of the artificial small parameter 6 explicitly. However, to obtain an estimate of the size of 6, we define [70] I du6 = sup max F(u; z)—j— , (3.10) Zm<Z<ZHC at j where the maximum is evaluated over at least one oscillation of the periodic orbit. Nu merical evaluations of the value of 6 are summarized in Table 3.2 and compared to the size of e. We observe that 6>> 6 > 0 for all models, verifying that we may indeed keep the damping term &(u; z)’à, while ignoring the EH(u, z) term in (2.44). 3.3.2 Derivation and Evaluation of the Melnikov Integral We first investigate the nature of the solution of the unperturbed version of (3.8), namely the Hamiltonian equation ii+G(u;z) = 0. (3.11) Chapter 3. Escape from the Active Phase 57 . 0 Cs Figure 3.8: Phase portrait of the unperturbed fast subsystem, (3.11), produced by DSTOOL, for a value of z between zm and zM. The separatrix (dashed curve) is de noted by (u0,o) and crosses the u-axis at u = a3, u = b3, and u c5. Like the corresponding fast subsystem, the fixed points of (3.11) are given by (n; z) = 0. For values of z between Zm and ZM, (3.11) has three fixed points and its corresponding phase portrait is shown in Figure 3.8. We introduce the following notation. The hyper bolic fixed point, or the saddle point, occurs at (u,à) = (a3z),0), where a3(z) is defined by the middle root of O(a3;z) 0, (3.12) and the separatrix, (uo(t; z), üo(t; z)), crosses the u-axis at u = b3(z) to the right of a3 and at u = c3(z) to the left of a3. When the perturbation due to the 6F(u, a; z) term is added back into the system, the separatrix itself is perturbed and breaks into a stable and unstable orbit. For values of z below ZHC, the unstable and unstable orbits separate so that the stable orbit lies outside the unstable orbit, whereas for values of z above zHC, the relative orientation of these (u0,tt) a3 u Chapter 3. Escape from the Active Phase 58 orbits is reversed (cf. Figure 3.3). The separatrix is perturbed but persists at z = ZHC. We now adapt the derivation of the Melnikov integral for a two-dimensional au tonomous system of equations perturbed by a time-periodic function (cf. Lichtenberg and Lieberman [58]). When the perturbation is non-autonomous, the Melnikov integral is used to predict the onset of chaotic motion, at some time t0. In (3.8), the perturbation is autonomous, and as such, the stable and unstable orbits either do not intersect or are entirely coincident [37]. Rather than obtaining a Melnikov distance function which depends on a time to then, we obtain a number. However, the number depends on the parameter z. Letting (x, y) = (u, ), (3.8) can be rewritten as th = y, (3.13) = —G(x; z) — F(x, y; z). (3.14) Perturbing the stable and unstable orbits, (xs(t; z), y3(t; z)) and (xu(t; z), yu(t; z)), re spectively, off the separatrix, (xo(t; z), yo(t; z)), we write xsu(t; z) = xo(t; z) + Sxr(t; z) + 0(62), (3.15) ysU(t; z) = yo(t; z) + 6y(t; z) + 0(62). (3.16) Substituting (3.15)-(3.16) into (3.13)-(3.14), expanding i’(x, y; z) and G(x; z) about (x, y) = (xo(t; z), yo(t; z)), and collecting like powers of 6 yields the integrable leading order equations tho = yo, (3.17) go = —O(xo; z), (3.18) Chapter 3. Escape from the Active Phase 59 N Figure 3.9: Geometrical representation of the Melnikov distance. The dashed curve is the separatrix of the unperturbed fast subsystem, as in Figure 3.8. There is a saddle point at (x,y) = (a,0). The solid curves are the stable and unstable orbits of the fast subsystem. The vectors d and N represent the difference between the stable and unstable orbits and the normal to the separatrix, respectively. (Adapted from [70].) defining the separatrix shown in Figure 3.8, where (x0, yo) replaces (u0,t0). The first- order equations are given by = Yi, (3.19) dt = —G(xo; z)x1 — F(xo, yn; z). (3.20) The stable and unstable solutions of (3.19)-(3.20) differ by d(t; z) = = ) + O(6), (3.21)Il —Y YiYi and the Melnikov distance D(t; z) then is defined to be the projection of dalong a normal N to the separatrix (x0, yo). The geometrical representation of the Melnikov distance is shown in Figure 3.9. From (3.17)-(3.18), we take Chapter 3. Escape from the Active Phase 60 -. I xo;z) ‘ N= I (3.22) \ Yo ) so that D(t;z) — N•d = S[yo(y -yr) +O(xo;z)(x _x)J +0(62) = 6 [(yoy + G(xo; z)x) - (yoy + (xo; z)x)] + 0(62). (323) To find an explicit expression for D, we write D = SD1 +0(62) 6[D —Dr] +0(62), (3.24) where D,U = SU + G(xo; z)xr. (3.25) Differentiating (3.25) with respect to time, we obtain ]8,U = + (xo; z)u + oYU + G(xo; z)hoxr. (3.26) Using (3.17)-(3.18) and (3.19)-(3.20) in (3.26) gives Ds,U = —F(xo,yo;z)yo. (3.27) Integrating (3.27) then yields D(; z) — D(to; z) = — j P(xo, yo; z)yodt, (3.28) D(to; z) — D(—; z) = — j F(xo, yo; z)yodt, (3.29) for the stable and unstable orbits, respectively. We point out that since we are interested in the value of ZHC at which there is a homoclinic orbit surrounding the right fixed point, we restrict xo(t) a3 in (3.28)-(3.29). Chapter 3. Escape from the Active Phase 61 In light of (3.12) and (3.17)-(3.18), we have 1o(t;z) = a3(z), (3.30) urn yo(t;z) = 0, (3.31) t—±oo Using these limits in (3.25), we thus have D(oo; z) = D(—; z) = 0, (3.32) so that from (3.24) and (3.28)-(3.29), we obtain the Melnikov integral up to 0(62), D(z) = 6J(xo,yo;z)yodt, = 6LF(u,n;z)dt, (3.33) where we have returned to u and ü variables and omitted the subscripts. However, we remember that the integral is to be evaluated around the right loop of the separatrix solution of (3.11) (cf. Figure 3.8). We now show that it is not necessary to determine the dependence of u and zt on t explicitly by converting the improper integral in t to a line integral in u. Integrating (3.11) from u = a8 to u and using the fact that ii = 0 when u a8, we have 1 U—2 + j G(u; z)dü = 0. (3.34) Hence, it = ±—2V(u; a3), (3.35) where V(u; a3) = j G(u; z)dã. (3.36) We note that (3.35) implies that the separatrix is symmetrical about it = 0 (cf. Fig ure 3.8). Since it = 0 when u = b8, b3 thus is given by V(b3;a3) = 0. (3.37) Chapter 3. Escape from the Active Phase 62 Using (3.30) and (3.35) in (3.33), we have D(z) = S F(u, it; z)du = [fa F(n, +—2V(u; a3); z)du + j fr(u, ——2V(u; a5); z)du]. (3.38) Finally, using the symmetry of the separatrix and the definition of (3.9) in (3.38), we obtain the S-independent Melnikov distance, D(z) 2fF(u;z)—2V(u;a3)du + x:8 [E(u, +—2V(u; a3); z) — E(u, ——2V(u; a3); z)] du. (3.39) Since E(u, 0; z) is even in it it does not contribute to the Melnikov integral. Furthermore, D(z) reduces to 2f F(u;z)—2V(u;a3)dufor models with p = 1, since E(u,it;z) does not depend on it for these models. Even though (3.39) is an analytical expression for the Melnikov distance, the nonlin earities in F and V force us to evaluate D(z) numerically. The algorithm to evaluate D(z) for each value of the parameter z is as follows: 1. Solve (3.12) to determinea3(z). We use a standard root solver (such as the RTSAFE routine from [74]) consisting of a combination of Newton’s method and bisection to guarantee a solution. 2. Solve (3.37) to determine b3(z), using the same root solver as above. 3. Compute D(z) from (3.39), using a standard numerical quadrature routine from QUADPACK [71]. Note that each evaluation of the integrand of D(z) requires the computation of V(u;a3(z)), which is another quadrature. Figure 3.10 shows two typical graphs of the Melnikov distance as a function of z, computed according to the algorithm described above. In Figure 3.lOa, the Melnikov Chapter 3. Escape from the Active Phase 63 (a) (b) ZM z Figure 3.10: Graphs of the Melnikov distance as a function of z, for two typical situations, cf. Figures 2.6a and b, respectively. distance changes sign once, at ZMEL given by D(ZMEL) = 0, corresponding to models for which the fast subsystem undergoes exactly one homoclinic bifurcation (cf. Figure 2.6a). In Figure 3.lOb, the Melnikov distance changes sign twice, corresponding to a bifurcation diagram where there are two homoclinic bifurcations (cf. Figure 2.6b). In this case, it is the first change of sign that is of interest to us. We note that for z < ZMEL, the Melnikov distance is positive, corresponding to the case where the stable orbit lies on the outside relative to the unstable orbit, so that stable periodic orbits surrounding the right fixed point are possible. The relative orientation of the orbits reverses near z = ZMEL. That is, zMEL is the leading-order approximation of zHC and, subsequently, Zesc. The values of ZMEL are determined by a bisection procedure on D(z) for each of the biophysical models, using an accuracy of at least 10—6 in each step of the algorithm described above and for the bisection. The results are presented in Table 3.1 and are in good agreement with typical z values for the Chay-Kang and Sherman-Rinzel-Keizer models. However, the results for the remaining models, especially the Chay (1986) model, can benefit from a higher-order correction, which is the subject of the section below. 0 Zm ZMEL 0 z Zm ZMEL ZM Chapter 3. Escape from the Active Phase 64 3.4 Approximation of Zesc with the Fredhoim Alternative In this section, we determine the first-order correction to the leading-order value of ZHC obtained with Melnikov’s method above. We note that Melnikov’s method can be used to derive the integrals required in each term of the (5-power series expansion of the distance between the stable and unstable orbits and, thus, it also can be used to derive the first- order correction to ZMEL. Instead, we demonstrate the use of a method based on a version of the Fredholm alternative due to Chow, Hale, and Mallet-Paret [20]. In the case of perturbed Hamiltonian systems, the Fredhoim alternative method is well-known to be equivalent to Melnikov’s method [51], in that it yields the same integrals. A clear introduction to the Fredhoim alternative is given in [52]. In Section 3.4.1, we derive leading-order and first-order problems. From the leading- order problem, we recover the Melnikov integral of the previous section, as expected. The first-order problem determines the required correction. We solve the equations nu merically in Section 3.4.2. The results indicate that the first-order correction is equal to zero for models with p = 1. The reason for this is exposed in Section 3.4.3. 3.4.1 Derivation of the Equations We rewrite the fast subsystem, (3.8), as ii + O(u; z) = —SE(u, it; z) (3.40) and seek an asymptotic solution of (3.40). Letting u(t) = uo(t) + Sui(t) +82u(t) + O(8), (3.41) = zo+6zi+(52+O(63), (3.42) substituting (3.41)-(3.42) into (3.40), expanding F and G about (u, it; z) = (u0, it0; zo), and collecting terms of equal magnitude gives an infinite sequence of equations. Chapter 3. Escape from the Active Phase 65 The leading-order equation is 9 + (zto; zo) = 0. (3.43) We thus have recovered the integrable Hamiltonian system defined by (3.11). Similarly, the first-order equation is üi + G(uo; zo)ui = —F(uo, z0; zo) — G(uo; zo)zi. (3.44) Although we determine the location of the homoclinic orbit up to first order only, we also require the second-order equation, which is U2 + (uo; zo)u2 = —O(uo; zo)z2 — [uu(uo; zo)u + 2(uo; zo)uizi + (uo;z0)z] — [fr(u0,it0;z0)u1 + (uo, it0; zo)iii + (u0,u0;z0)z1]. (3.45) From (3.44) and (3.45), we observe that u, for n 1, satisfies the general equation of the form ü, + G(uo, zo)u = T(u0,u1, ..., it0, itt, ..., it,; zo, z1, ..., zn_i), (3.46) where T is determined by F and and their derivatives. We define the differential operator L by Lu = ii + (uo; zo)u, — <t < , (3.47) so that (3.46) can be rewritten as Lu,., = T. (3.48) According to the version of the Fredhoim alternative for solutions bounded on R due to Chow et al. [20], solutions to (3.48) exist if and only if (T,v) = 0 (3.49) Chapter 3. Escape from the Active Phase 66 for all v for which L*v = 0, where L* is the adjoint of L and (u, v) is the inner product defined by (u, v) = t: u(t)v(t)dt. (3.50) In other words, solutions exist if and only if T is orthogonal to any v in the nulispace of L*. Since L is self adjoint, we require a v so that Lv = i + G(uo; zo)v = 0. (3.51) By differentiating (3.43) with respect to t, it is clear that Üo is always a solution of the homogeneous equation (3.51). From (3.49) and (3.50), solutions to (3.48) thus exist if and only if L T(t)it0(t)d = 0. (3.52) We now proceed by investigating how the Melnikov integral condition can be recov ered from the above systems of equations. As mentioned in the previous section, the Hamiltonian system, (3.43), has a separatrix for every value of z0 between Zm and zM. It therefore is not sufficient to simply consider this system of equations by itself. To obtain a restriction on the values of z0, it is necessary also to consider the solvability condition for the u1 equation. Applying (3.52) to (3.44), the required solvability condition is f [(u0,ü0; zo) + (uo;z0)z1]t0dt = 0, (3.53) which is an integral equation with both z0 and z1 unknown. However, since G(uo;z0)z1 does not depend explicitly on Ito, the contribution of this term to the integral is zero, and (3.53) reduces to i_: fr(0,Ito; zo)Itodt = 0. (3.54) Thus, we obtain exactly the same integral condition as was obtained using the Melnikov distance approach by setting D(z) = 0 in (3.33), which has a solution for precisely one Chapter 3. Escape from the Active Phase 67 value of z0. In fact, z0 zMEL. However, we retain the z0 notation to distinguish between the different methods used to determine these two values of z. Similarly, to determine the value of z1, it is necessary to consider the solvability condition for the u2 equation. Applying (3.52) to (3.45), and remembering that there is no net contribution to the integral by terms which depend on o and z0 only, we require = 0, (3.55) where M(uo, u1, t0,z1; zo, zi) = [Ouu(uo;z0)u + 2(uo; zo)uizi + (uo;z0)z] + [P(u0,z0; zo)u1 + P(u0,zo; zo) + P(uo, o; zo)z1]. (3.56) We have just discussed the conditions that need to be satisfied to uniquely determine ZO and z1. Similarly, of course, we need conditions to uniquely determine u0 and Ui. Recall that we are interested in a solution of (3.40) which is homoclinic to the saddle point on the O = 0 nullcline at a8(z), for some value of z. Therefore, we require lim u(t) =a3(z). (3.57) t—±cc Using (3.12), and (3.41)-(3.42) in (3.57), the boundary conditions at infinity are limuo(t) = a8(zo) = a, (3.58) da3(zo) G(a;zo)lim Ui(t) = z1 = —z1 , (3.59) t—*±co dz G(a;zo) for the leading-order and first-order problems, respectively. Since the equation for each u is a second-order equation, these conditions seem sufficient on first glance. However, as discussed above, ü0 always is a solution of the homogeneous equation (3.51). As such, Chapter 3. Escape from the Active Phase 68 any multiple of It can be added to the solution u, of (3.46), for n 1. To fix the multiplicative constant, an additional condition is required. To obtain this condition, we use the fact that the homoclinic orbit must pass through It = 0 for some t between —00 and co. Without loss of generality, we take It(0) = 0, (3.60) yielding It(0) = 0, (3.61) for all values of n. We now return to the integral defined by (3.55). This integral cannot be converted to an integral in the phase plane (cf. Section 3.3) since it is not possible to obtain a solution for u1 and It1 as functions of u0 and It0. Because Ito is a solution of the homogeneous equation of (3.44), in theory, it is possible to write down an analytical solution for u1 and It1 as functions of t, using reduction of order and the method of variation of coefficients [12]. However, numerical methods then must be used to evaluate u1(t) and It1(t) in (3.55). We prefer to obtain a numerical solution directly, discussed in the next section. 3.4.2 Numerical Solution We pose the problem as a boundary-value problem (BVP) on the domain tL t tR, where tL and tR are as large as allowed by available computer resources, and use the COLSYS code written by Ascher et al. [2]. We start by setting up the BVP for determining the value of z0. Although z0 is of course more easily determined by converting (3.54) into a line integral in u, as was done in Section 3.3, solving this problem numerically is an important step towards determining z1 below. Naturally, we require the leading-order equation, (3.43). The corresponding Chapter 3. Escape from the Active Phase 69 boundary conditions for u0 from (3.58) become uo(tL) = uo(tR) = a, where a is an unknown constant. To let COLSYS determine the value of a3, we introduce the differen tial equation ã = 0, with the corresponding condition G(a, z0) = 0 at either the left or right boundary. To uniquely determine z0, we include the solvability condition (3.54) as follows. We define h(t) = J fr(uo(r), o(r); zo(r))o(r)dr, (3.62) and use the Fundamental Theorem of Calculus to introduce an artificial differential equa tion for h(t), with corresponding boundary conditions, h(tL) = h(tR) = 0. To summarize, we thus arrive at the following system of equations o + (uo; z0) = 0, (3.63) = 0, (3.64) = 0, (3.65) = P(uo, zo;z0ü, (3.66) for the five unknowns u0, it0, zo, a, and h. The required five boundary conditions are / “ / “ 0 zLojtL) = UOI\tR) = a3, O(a;zo) = 0, att=tL, (3.68) k(tL) = h(tR) = 0. (3.69) This system of equations uniquely determines the value of z0, with essentially identical results obtained from COLSYS (within the requested accuracy of 10—6) and from the Melnikov distance approach in Section 3.3 (see Table 3.1). To determine z1, we include the first-order equation, (3.44), and the corresponding boundary conditions for u1, (3.59). Since we now explicitly include the equation for u1, we omit (3.66) and (3.69) and replace them by the solvability condition for z1, (3.55). Chapter 3. Escape from the Active Phase 70 We define k(t) = f M(uo(r), ui(r), o(r), i(r); zo(r), zi(r))o(r)dr (3.70) and, as before, use the Fundamental Theorem of Calculus to introduce a differential equation for k(t), with the two boundary conditions k(tL) = k(iR) = 0. The net result is an addition of three equations, but only two boundary conditions. We now recall that any multiple of ho can be added to the solution for u1, as it is a solution of the homogeneous equation for u1. To uniquely determine the multiplicative constant, we use the additional condition (3.61). To summarize, we thus have the following system of eight equations, iio + uo; ZO) = 0, (3.71) = 0, (3.72) = 0, (3.73) h + (uo; zo)ui = —(zi0,h0; zo) — G(uo; zo)z1, (3.74) = 0, (3.75) k = M(uo,ui,ño,h1;zo,zi)ü , (3.76) for the eight unknowns u0, h, zo, a, u1, iti, z1, and k. The corresponding eight boundary conditions are uo(tL) = uo(tR) = a, (3.77) O(a;zo) = 0, att=tL, (3.78) = 0, (3.79) u1(tL) = u1(tR) = _1z(a;zo), (3.80) G(a8;zo) k(tL) = k(tR) = 0. (3.81) Chapter 3. Escape from the Active Phase 71 This system of equations uniquely determines the value of z1. The results obtained are presented in Table 3.1. We see that for the models with p 1, namely the reduced Chay-Keizer, Chay (1986), and Himmel-Chay models, the first-order correction does indeed improve the leading-order results, as expected, so that zo + 6z1 (see Table 3.2 for values of 6) agrees much better with typical values of z€8 than z0 = ZMEL alone. For the models with p 1, namely the Chay-Kang, Sherman-Rinzel Keizer, and Chay-Cook models, z1 = 0 within the requested accuracy. We shall see in the next section that this is not coincidental and that, in fact, z1 is identically equal to zero. The leading-order results already were in excellent agreement with typical values of Zesc for the Chay-Kang and Sherman-Rinzel-Keizer models and these are not changed with the addition of 6z1 = 0. Unfortunately, the results also remain unchanged for the Chay-Cook model. To explain the remaining discrepancy for this model, we return to the investigations in Section 3.3.1. We hypothesize that the damping term 6(u, ü; z) in (3.9) is not small enough for this model. To test this hypothesis, we replace 6fr(u, t; z) by 6P(u, i; z). The corresponding modified exact value of z for a typical value of j9 is approximately -0.694, yielding a normalized ZMEL = z0 value of 1.06419, which is better than the 0.90639 of Table 3.1, as expected. 3.4.3 Why The First-Order Correction is Zero for Models with p = 1 In this section, we consider the symmetry of the leading-order and first-order solutions for u, u0, and u1, respectively, and examine the effect on the solvability condition for u2 to show that z1 0 for models with p = 1. For these models, we can rewrite (3.40) as ii + (u; z) = —S1’(u; z), (3.82) Chapter 3. Escape from the Active Phase 72 where F(u; z) z) (3.83) from (3.9), since E(u, ; z) — E(u, 0; z) 0 from (2.46). The corresponding leading-order equation for u, (3.43), does not depend on P, so that the separatrix solution is not affected. Changing t to —tin the leading-order problem, (3.43), it is clear that u0 is even in t and, hence, ii is odd in t. Due to (3.82), the first-order equation for u, (3.44), can be rewritten as üi + u(uo; zo)ui = —E’(uo; zo)üo — (uo; zo)zi. (3.84) Letting i(t) be the bounded solution of i + O(uo; zo)1’ = —E(uo; zo)io (3.85) and ‘(t) be the bounded solution of + u(uo;z0) = —(uo; z0), (3.86) we have = + z1. (3.87) Since u0 is even and it0 is odd in t, it is clear from (3.85) that Y is odd in t and, hence, is even. Similarly, 1’ is even in t from (3.86) and, hence, is odd. Furthermore, Y = To see this, we differentiate (3.43) with respect to z0 to obtain + (uo;0)- = —(uo; z0). (3.88) That is, is a particular solution of (3.86), which is bounded because a separatrix solution exists for all z0 between Zm and zM. The complete solution ‘ = then follows c9Zo by requiring that ‘(t) is bounded [20] and Y(O) = 0 from (3.61). Chapter 3. Escape from the Active Phase 73 Due to (3.82), the second-order equation for u, (3.45), can be rewritten as U2 + (uo; zo)u2 = —(uo;z0)z2 — [uu(no;z0)u + 2(uo; zo)uizi + (uo;z0)z] [u(uo; zo)toui + P(uo;z0)z1 + E(uo; zo)zozi], (3.89) yielding the modified solvability condition L Af(uo, Ui, 1o, it1; zo, zi)itodt = 0, (3.90) for u2, where M(uo, u1, 1t, itj; z0, zi) = [uu(uo; zo)u + 2(uo; zo)uizi + (uo; zo)zj + [Eu(uo; zo)itoui + P(uo; zo)it1 + P(uo; zo)itozi], (3.91) similar to (3.55)-(3.56). Using (3.87) in (3.91) gives A! = MaZ + Mbzl + (3.92) where = (uo;z0)/2 + (uo; zo) + O(uo; zo), (3.93) Mb = + E(uo; zo)ito + E(uo;z0)+ E(uo;z0)it, (3.94) i2i = O(uo;z0)i2+ E(uo; zo)itoi’ + E(uo; zo)/. (3.95) Apparently, (3.92) is quadratic in z1. However, Ma is orthogonal to it0, i.e., 11ait0dt = 0, (3.96) Chapter 3. Escape from the Active Phase 74 so that from (3.90) and (3.92), z1 is uniquely defined by f Il[ztodt = — . (3.97) jMbüodt We note that since IIa does not depend on P, the same reduction holds for models with p 1 as well, so that z1 also is uniquely defined for those models. To prove (3.96), we differentiate (3.88) with respect to z0, giving - 92u0 2 +G(uo,2o) = — [uu(uo; zo) (a)2+ 2O(uo;z0)+ (uo;z0)] = — {ãuu(uo;z0)2 + 2G(uo; zo) + G(uo; z0)] = 2Ma. (3.98) Since is bounded, the right-hand side of (3.98) is orthogonal to o and, hence, (3.96) follows. Using the symmetry of u0, Y, and Y in (3.95), we find M is even, so that the numerator in (3.97) vanishes. Similarly, it is clear that every term in I[b is odd and contributes to the integral in the denominator, so that it is non-zero in general. Therefore, z1 0 for all models with p = 1. 3.5 Comparison of the Different Methods In this chapter, we have investigated several methods for approximating zHC, which is the value of z at which the fast subsystem, (2.51), exhibits the homoclinic bifurcation of interest. The approximate value of ZHG in turn approximates Zesc, which is the value of z at which the numerical solution of the full system of equations undergoes the transition from the active phase back to the silent phase. Chapter 3. Escape from the Active Phase 75 We demonstrated the use of two numerical methods to determine approximations to ZHC, as well as Melnikov’s method and the Fredhoim alternative method. The latter two methods are analytical, but need to be combined with some numerical calculations in order to yield explicit approximations. We mention again that the dependence of Zesc on the glucose-dependent parameter /3 cannot be determined with any of these methods since /3 does not appear in the fast subsystem. That is, all methods have resulted in a constant, approximate value for Zesc. The results from the two numerical methods, namely the use of AUTO to approximate the homoclinic orbit by a periodic orbit with very large period and bisection based on the bifurcation of the solution behaviour of the fast subsystem in the (n, it) phase plane as z is varied, are in excellent agreement with typical exact values of Zesc for all models. Since consistently accurate results can be obtained very quickly, we conclude that these methods are the methods of choice when one is interested in the location of the homoclinic bifurcation for a given set of parameters in the fast subsystem, as is the case in this thesis. In Chapter 4, we thus use Zesc ZAUTO ZBJS. In general, however, the numerical methods are inefficient for determining the de pendence of the location of the homoclinic orbit on the values of the fast subsystem parameters. In that case, approximating zHc with the leading-order Melnikov method or, equivalently, the leading-order Fredholm alternative method, is very efficient, since the corresponding leading-order improper integral in t can be converted into a line integral in the (u, ii) phase plane. Since for the purposes of the next chapter we are interested in the dependence of ZHC on the parameter /3 only, which does not appear in the fast subsystem and, therefore, does not affect the analytical approximations of ZHC, the advantage of using this method is not fully demonstrated here. However, we demonstrate its efficient use in Chapter 5. We note that the drawback of this method is twofold. First, the fast subsystem needs to be written in perturbed form. We were able to do so for all models Chapter 3. Escape from the Active Phase 76 under consideration in this thesis, but there is no guarantee that this also is the case for other models. Second, the methods are asymptotic methods which do not always provide accurate results, and higher-order corrections may be required. Higher-order corrections of the leading-order approximations of zHC can be obtained using either Melnikov’s method or the Fredholm alternative method. However, the cor responding higher-order improper integrals no longer can be converted to simple line integrals in the (u, ‘a) phase plane. We posed the problem of obtaining the first-order correction as a boundary-value problem and solved it with COLSYS. The disadvantage of this approach is that it is computationally intensive and, in fact, less practical than using the AUTO and DSTOOL software packages. Chapter 4 Determination of the Plateau Fraction In this chapter, we investigate the dependence of the plateau fraction, which was defined as the ratio of the active phase duration to the burst period, on the glucose-dependent parameter • One can obtain an exact graph of the plateau fraction as a function of /3 by numerically integrating the full system of equations, (2.44)-(2.45), over a range of 3 values and extracting the plateau fraction from the solution for each value of /3. The problem with this approach is that it takes a relatively large amount of computer time to obtain an accurate graph. Moreover, this approach does not provide any physiological insight into the dependence of the plateau fraction on 3, let alone any other model parameters. Instead, we derive an expression for the leading-order approximation of the plateau fraction using asymptotic methods, assuming that the durations that the solution spends in the transitions between the active and silent phases are very short. We extend the techniques of Pernarowski [67] and Pernarowski et al. [70], which were applied to the Sherman-Rinzel-Keizer model, to all models under consideration. In Section 4.1, we decide on a range of values for 3 over which to compute the plateau fractions. The exact silent and active phase durations and the plateau fraction are defined in Section 4.2. In Section 4.3, we derive a leading-order approximation for the silent phase duration, T3°, and compare the results to the silent phase duration obtained numerically. A leading-order approximation of the active phase duration, is derived in Section 4.4, and also compared to the numerical results. In Section 4.5, we discuss the resulting leading-order approximation of the plateau fraction, p. Finally, 77 Chapter 4. Determination of the Plateau Fraction 78 the assumption that the transition phases are very short is re-examined in Section 4.6. 4.1 Range of Values of 3 To decide on a range of values for 3 over which to compute the plateau fractions, we examine the effects of varying 3 on the solution of the full system of equations, (2.44)- (2.45). From (2.50), we note that K(u, z) depends on 3. Varying the value of /3 affects the location of the -nullcline, given by K(u, z) = 0 from (2.45), and, subsequently, the location and stability of the fixed point of the full system of equations. For high values of 3, corresponding to low values of the dimensional ka parameter, the -nul1cline intersects the O = 0 nulicline on the left branch. For these cases, the fixed point of (2.44)-(2.45) is stable. As a consequence, the solution of (2.44)-(2.45) remains in the silent phase and exhibits no bursting behaviour. Decreasing /3 shifts the -nullcline and, therefore, the fixed point of the full system of equations to the right. Bursting behaviour becomes possible when this fixed point occurs on the middle branch of the fast subsystem nullcline and the fixed point of the full system of equations becomes unstable. We use (2.50) to compute the critical value of /3, = e_zmh(um,zm)’ (4.1) where (Urn, Zm) represents the local minimum of the fast subsystem nullcline. For /3 = the -nullcline intersects the G = 0 nullcline at (Urn, Zrn), as shown in Figure 4.1. We note that if /3 is just slightly smaller than j3rn, the -nullcline lies entirely to the left of the branches of periodic orbits emanating from the Hopf bifurcation. That is, z strictly increases during the active phase. Decreasing /3 further eventually causes part of the -nullcline to lie inside the branches of periodic orbits, so that z no longer strictly increases during the active phase. When Chapter 4. Determination of the Plateau Fraction 79 Figure 4.1: Dependence of the -nullcline on 3. For 3 = /3m, the -nu11cline intersects the C = 0 nulicline at its local minimum, (am, Zm). For /3HC, the -nullc1ine intersects the C = 0 nulicline at approximately (ZLHC, zHc). Note that I3HC <urn, i.e., the -nul1c1ine moves to the left as /9 increases. /3 is small enough, the decrease in z balances the increase in z during one oscillation, resulting in continuous spiking and no escape from the active phase back to the silent phase. We restrict our analysis to the values of 3 for which an escape from the active phase is assured. We therefore choose the minimum value of 3 to be /9HC = () , (4.2)e Hch(uHczHC) where z0 is the best approximation of ZHC obtained in Chapter 3, namely ZAUTO or ZBIS, and satisfies G(u, z) = 0 on the middle branch. For /3 = !3HC, the - nulicline thus intersects the = 0 nulicline at approximately (UHC, zHc), as shown in Figure 4.1. U Chapter 4. Determination of the Plateau Fraction 80 4.2 Definition of the Exact Silent and Active Phase Durations and the Plateau Fraction We now define the exact silent and active phase durations, as well as the plateau fraction. Figure 4.2 shows the bifurcation diagram of (2.51) for a typical 3-cell model, with the projection of the numerical solution of (2.44)-(2.45) onto the (u, z) plane superimposed. Based on this figure, we identify four phases, namely the silent and active phases and two transition phases. Capture refers to the transition from the silent to the active phase and escape refers to the transition from the active phase back to the silent phase. During the silent phase, the solution trajectory closely follows the left branch of the G = 0 nullcline. We define the exact silent phase duration to be T3 = tcap — tsji, (4.3) where tcap > t8j, and tcap and t81 are defined as follows. t is the time at which the solution trajectory crosses the left branch of the G = 0 nulicline after escaping from the active phase. More precisely, t91 is the time at which u uj and z z81, where and are related through (u8j, = 0 on the left branch. Similarly, t, is the time at which the solution trajectory passes below the local minimum of the fast subsystem nullcline. In particular, at t = tcap, U = Ucap = Urn, and z = Zcap. During the active phase, the solution trajectory exhibits oscillations surrounding the right branch. Let t. be the time at which u reaches a local maximum during the jth spike. We then can define the exact active phase duration to be Ta = tN — ti, (4.4) where tN > t1, N is the number of spikes per burst and t1 is the time at which u reaches the first local maximum. To compute Ta accurately, the solution of (2.44)-(2.45) must be reported at small step sizes. We have used step sizes of 0.001 for all of the models. Chapter 4. Determination of the Plateau Fraction 81 -2.2 -2.4 -1.15 -0.15 Figure 4.2: Typical bifurcation diagram of the fast subsystem, (2.51), with the projection of the numerical solution of the full system of equations, (2.44)-(2.45), superimposed. Results are shown for the Himmel-Chay model. Corresponding graphs for the remaining models under consideration in this thesis are qualitatively identical. The labels (Urn, Zrn) and (UM, ZM) indicate the local minimum and maximum of the G nullcline, respectively. The silent phase is defined to begin when the solution trajectory crosses the left branch of the G = 0 nullcline at (u3ü,z3) and defined to terminate when the solution trajectory passes below the local minimum of the G = 0 nullcline at (ucap, zcap). The active phase is defined to begin when u reaches its first local maximum after the silent phase, at label 1, and ends when u reaches its last local maximum before returning to the silent phase, at label N. The transition phases between the active and silent phases are referred to as capture and escape, as labelled. Also shown is (Uec, Zesc), at which the active phase undergoes the transition from the active phase back to the silent phase. -1 -1.2 -1.4 -1.6 -1.8 -2 -0.95 -0.75 -0.55 -0.35 U Chapter 4. Determination of the Plateau Fraction 82 The number of spikes per burst, N, is dependent on /3, as shown in Figure 4.3. We can explain the decrease of N with /3 as follows. Since > 0 in the active phase, (2.45) and (2.50) imply that an increase in /3 corresponds to an increase in during the active phase. That is, z = Zesc is approached more rapidly, resulting in a reduction in the number of spikes per burst, as observed. We note that a theoretical treatment of the transition from n to n + 1 spikes per burst as parameters are varied is provided in [90]. The silent and active phases are separated by two transition phases, during which the solution trajectory switches from one to the other. The transitions occur very fast, and we ignore the time that the solution spends in the transition phases in our approximation of the plateau fraction. Letting T be the burst period, we thus can define the exact plateau fraction, p, to be Ta Ta 6Ta+Ts (4.5) Graphs of the exact silent and active phase durations, burst period, and plateau fraction as a function of /3 are presented in the following sections, where comparisons between the exact and corresponding leading-order quantities are made for all models under consideration. 4.3 Leading-Order Approximation of the Silent Phase Duration To derive the leading-order approximation for the silent phase duration, T9°, we examine the full system of equations, (2.44)-(2.45). From (2.45), the change in z is 0(e). From Figure 4.2, we note that du/dz = 0(1) during the silent phase. Thus the change in u also must be 0(e) during the silent phase. For this reason, we introduce a new slow time, , and new dependent variables, U and Z, as follows: = et, (4.6) U(Z;e) = u(t;e), (4.7) Chapter 4. Determination of the Plateau Fraction 83 (a) 28 , i (b) 14 8l1L 2k111 15 20 25 30 35 40 0.01 0.0175 0.025 0.0325 (c) 25 ________________ (d) 35 iii 4 6 8 10 12 14 0.07 0.11 0.15 0.19 0.23 0.02 0.025 0.03 0.035 0.04 6.6 7 7.4 7.8 8.2 8.6 1 Figure 4.3: Dependence of the number of spikes per burst, N, on 3, for /3HC /9 /3m. All other parameter values are as listed in Appendix A. Results are shown for the (a) reduced Chay-Keizer model; (b) Chay (1986) model; (c) Himmel-Chay model; (d) Chay-Kang model; (e) Sherman-Rinzel-Keizer model; (f) Chay-Cook model. Chapter 4. Determination of the Plateau Fraction 84 Z(1;E) = z(t;e), (4.8) so that both U and Z both undergo 0(1) changes on the time scale. Substituting (4.6)-(4.8) into (2.44)-(2.45), we obtain e2U” + E(U,U’, Z) + eF(U, Z)U’ + G(U, Z) = eH(U, Z), (4.9) Z’ = K(U,Z), (4.10) where’ denotes differentiation with respect to L System (4.9)-(4.10) is a singular perturbation problem. We seek an outer solution by substituting the expansions U() Z(1;e) Zo(t)+eU1(t)+..., into (4.9)-(4.10), expanding F, F, G, H, and K using the Taylor series about the leading- order solution (U, U’, Z) = (U0,U, Z0), and expanding E(U0,eU, Z0) using the binomial series. Collecting the leading-order terms yields G(U0,Z) = 0, (4.11) Z = K(U0,Z), (4.12) where is as defined by (2.53). The solution of the leading-order silent phase equations, (4.11)-(4.12), thus follows the G = 0 nullcline, as expected from Figure 4.2, with the change in Z0 governed by (4.12). The accuracy of the time dependence of the leading order problem is confirmed in Figure 4.4. The solution of the leading-order silent phase problem terminates when (u, z) = (tt,, Zm). We see that the duration of the transition from the silent phase to the active phase, starting when u tim and terminating when u reaches its first local maximum, appears to be significant. We examine the transition durations in more detail in Section 4.5. Chapter 4. Determination of the Plateau Fraction 85 (a) -0.2 -0.4 - -0.6 - -0.8 - - -1 - - - - -- - -- - - - — - - -1.2 I I I I I 0 200 400 600 800 1000 1200 t (b) -5.05 -5.15 - - -5.25 - - -5.35 - - -5.45 - - 555 I I I I I 0 200 400 600 800 1000 1200 I Figure 4.4: Comparison of the solution of the full system of equations (solid curves), (2.44)-(2.45), and the solution of the leading-order silent phase problem (dashed curves), given by (4.11) and (4.13). The initial conditions are (u, z) = z3). Results are shown for the Chay-Kang model for /3 = 8.5301e2 All other parameter values are as listed in Appendix A. Corresponding graphs for the remaining models under consideration in this thesis are qualitatively the same. Chapter 4. Determination of the Plateau Fraction 86 To obtain the leading-order silent phase duration, T8°), we differentiate (4.11) and use (4.12) to substitute for Z. The resulting differential equation for Uo is u’ — K(U0,Z) 413° —(Uo,Z)/O(U, Finally, we separate variables in (4.13), integrate, and return to the original variables t and u via (4.6)-(4.7) to give T° — 1 [Urn du (4 14 — 6Ju° N(u)’ where N(u) = K(u,q(u)) (4.15) —G(u, g())/G(, (u))’ and c(u) satisfies O(u, g(u)) = 0. The limits of integration in (4.14) are chosen as follows. According to the definition of the exact silent phase duration in (4.3), the limits of integration for u should be from u31 to Ucap. The upper limit of integration presents no problems, as it is a straight forward exercise to compute the local extrema of the fast subsystem nullcline and obtain Ucap = Urn (cf. Figure 4.2). However, u3 cannot be determined a priori since it depends on where the solution escapes the active phase. We define the leading-order approximation of u3 to be u, satisfying G(u), z) = 0 on the left branch, where z is the best approximation of Zese, the value of z at which the solution of the full system of equations undergoes the transition from the active phase back to the silent phase, determined in Chapter 3. Due to the complicated functional form of N(u), T° must be evaluated numerically. We use a standard numerical quadrature routine from QUADPACK [71]. In general, for each evaluation of the integrand in (4.14), a standard root finding procedure must be used to determine z = g(u). Rather than integrating (4.13), it seems just as (dis)advantageous to integrate (4.12) and use the root finding procedure to determine u as a function of Chapter 4. Determination of the Plateau Fraction 87 z instead. However, for four of the models under consideration, namely the reduced Chay-Keizer, Chay-Kang, Sherman-Rinzel-Keizer, and Chay-Cook models, z = q(u) can be found explicitly. For these models, N(’u) can be reduced to N(u) = K(u;c(u)) (4.16) greatly simplifying the computations. None of the models permit a solution for u as a function of z. In Figure 4.5, the exact and approximate silent phase durations as a function of /3 are compared for all six models under consideration. The jaggedness of the exact silent phase duration curves is due to the discontinuities in the N curves illustrated in Figure 4.3. The agreement is very good for all models (although one may argue that the results for the Himmel-Chay model are not quite as good as for the remaining models). We explain the general increase in the silent phase duration with /3 as follows. Using (2.50) in (4.15), we differentiate (4.14) with respect to 3 to obtain dT8°) — 1 [Urn e)h(u,c(u))Ou(u,c(u))/z(u,c(u))d > 0 (417)d/3 — e JU? L8e_c(U)h(u, (u)) — 112 since O(u, c(u))/ã(u, g(u)) = —dz/du > 0 on the left branch of the = 0 nullcline, and h(u, c(u)) > 0. Alternatively, we recall that < 0 during the silent phase. From (2.39), an increase in /3 thus implies slower dynamics for z and, hence, an increase in the silent phase duration. 4.4 Leading-Order Approximation of the Active Phase Duration To derive a leading-order approximation of the active phase duration, it is necessary to first derive appropriately scaled equations. As the solution enters the active phase, the change in u is no longer O(e), but 0(1). Due to the existence of oscillations on the 0(1) time scale, the ii term is also important. Therefore, (2.44)-(2.45) are appropriately Chapter 4. Determination of the Plateau Fraction 88 (a) 1000 i (b) 1400 900 950 1200 - E 850E 800 750 - - 700 400 15 20 25 30 35 40 0.01 0.0175 0.025 0.0325 1 13 (c) 800 i i i (d) 1600 1500 -750 1400 -700 1300 -650 500 900 1200 - 600 1100 - 550 1000 - 4 6 8 10 12 14 0.07 0.11 0.15 0.19 0.23 1 3 (e) 900 (f) 2000 700 - 1800 - 600 - 1700 - _ >1800_,) 1900-500 - 1600 -400 - 1500 -300 1400 0.02 0.025 0.03 0.035 0.04 6.6 7 7.4 7.8 8.2 8.6 Figure 4.5: Comparison of the exact silent phase duration (jagged curve), T5, and the leading-order silent phase duration (smooth curve), as a function of 3, for 13Hc /3m All other parameter values are as listed in Appendix A. Results are shown for the (a) reduced Chay-Keizer model; (b) Chay (1986) model; (c) Himmel-Chay model; (d) Chay-Kang model; (e) Sherman-Rinzel-Keizer model; (f) Chay-Cook model. Chapter 4. Determination of the Plateau Fraction 89 scaled for the active phase. However, from the analysis in the previous chapter, we know that the active phase problem can be viewed as a slowly varying oscillator problem. The oscillations in u occur on a fast, 0(1), time scale, but the slow evolution of z and, hence, the slow variation of the oscillations occurs on a much slower, 0(e), time scale. Therefore, it is natural to use a multiple scales analysis [9, 56] of the active phase, taking both time scales into account simultaneously. Kuzmak [57] and Luke [59] developed a multiple scales perturbation method for strongly nonlinear oscillators to determine solutions which are valid for times up to 0(1/c). In light of the reformulated fast subsystem equation, (3.5), and its unperturbed version, (3.11), the active phase indeed can be viewed as a strongly nonlinear oscillator so that we can formally apply the Kuzmak-Luke method. Let u(t;e) = U(r,E;e) u0(r,E) + eu1(r,) +..., (4.18) z(t;e) = Z(r,t;e) zo(r,t) +ezi(r,i) +..., (4.19) where the slow time ‘ and the fast strained time r are defined by =et (4.20) and =, (4.21) respectively. The frequency function w() is to be determined by the requirement that ZL0 be periodic in r with a period which does not depend on 1. Furthermore, u(r, ‘) and z(i-, ) are assumed to be periodic in T, for i = 0, 1,2 Even though the exact solutions u and z are functions of t alone, we seek solutions which are functions of both and r treated as independent variables. Using the chain rule for partial differentiation, Chapter 4. Determination of the Plateau Fraction 90 the time derivatives of u are du OU OU — — Or 9t 0u I 0u1 0u Or + + 0(e2) (4.22)= and d2u 02U dc1..’OU 20U — + 2ew + e —- — + e —-OrOt dt Or Ot 20u / 20u1 02u &,Ouo) = w_—-+ew-ã---+ -+ +0(e2), (4.23)OTOt dt Or with similar expressions for the time derivatives of z. Expanding E(u, it, z) and F(u, z) Ouo about (uo, w-—, zo) and (u0, zo), respectively, gives E(u, it, z) = E° + e (Eo)U1 + E° r 0u1 + E0)z1) + 0(e2), (4.24)Iw—+L Or 2’F(u, z) = F° + e (Fo)u1 + F0)z1)+ 0(e ), (4.25) Ouo where the superscript (0) means evaluated at (uo, w--—, z0) or (u0, zo), as appropriate. Substituting (4.22)-(4.23), (4.24)-(4.25), and corresponding expansions for G, H, and K into (2.44)-(2.45), yields the leading-order equations 20u OU0 Ouj OT2 uo,zo)_ã__+G(uo,zo) = 0, (4.26)Or Ozo — = 0, (4.27)Or and the first-order equations O2u1 W Or2 + ‘‘ [0) + F(°)] + [Eo + wF° + G°1 U1 =Or uj H° — 2w Ouo dw OU0 OrO — [E°) + F(°)] — [Eo’ OUO (0) + G0)] 1, (4.28)Or Ozi 0z + -- = K(uo, zo). (4.29) Chapter 4. Determination of the Plateau Fraction 91 From (4.27), we note that z0 is a function of only. With the restriction imposed on the glucose-dependent parameter in the introduction to this chapter, z0 increases monotonically during the active phase. We thus can write as a function of z0 and infer the existence of some function w0 so that w(t) = WO(ZO). (4.30) The change of variables r = s, (4.31) uo(r,t) = Uo(s,t), (4.32) transforms (4.26) into + E(Uo, zo) + F(U0,z0)9- + G(U0,z0) 0. (4.33) For a fixed value of z0 between Zm and ZHC, (4.33) is equivalent to (2.51) and known to have a periodic solution U0 in s from the results in Chapter 3. Letting the period be To(zo) and wo(zo) = 1 , (4.34)T0(ZO) (4.31) implies r = s/To(zo). That is, (4.26) has a periodic solution u0 in T, with constant period 1, as required. Next, we integrate (4.29) with respect to r, from r = 0 to r = 1. Since z1 is assumed to be periodic in r as well, the contribution due to the first term is zero. We thus obtain the equation governing the slow evolution of zo, = k(zo()) = e0(zo()) -1, (4.35) where we have used (2.50) to rewrite K and = f h(uo(r, ), ZO(T, ))dr (4.36) Chapter 4. Determination of the Plateau Fraction 92 is the average value of Ii over one oscillation of the periodic solution. Returning to the original variables u, z, and t, the leading-order active phase problem thus is given by (2.51) and e [eh(z) - ij, (4.37) where — 1 To(z) h(z) = To(z) J h(u(t),z)dt. (4.38) For any value of z between Zm and z, the period To(z) can be obtained easily from the periodic solution of (2.51). To determine i(z), we define I(t, z) = h(u(), z)dt (4.39) and use the Fundamental Theorem of Calculus to introduce the differential equation = h(u(t),z). (4.40) We solve (2.51) and (4.40) simultaneously, and take (z) = urn &(t + To(z), z) - (t, z) (4.41) t—oo To(z) The period To(z) and averaged calcium current iz(z) so obtained are shown in Fig ures 4.6a and b, respectively. We note that as z approaches z, To(z) increases signif icantly, as expected, since the larger z becomes, the closer the trajectory passes to the saddle points on the middle branch of the G = 0 nulicline. Because h(u, z) is smallest for values of u near the middle branch, i(z) decreases significantly as z approaches Zj. In order to obtain a numerical solution of the leading-order active phase problem, ideally, i(z) should be known for every value of z between Zm and zj. However, as we have just seen, obtaining an accurate value for i(z) is computationally expensive, and as such, we look for a more efficient, yet accurate, approximation scheme. Whereas Chapter 4. Determination of the Plateau Fraction 93 (a) 24 z (°)m ZHC (b) 0.52 0.34 ‘ (J)Zm ZHC Figure 4.6: (a) Period To(z) of the periodic solution of (2.51). (b) Averaged calcium current h(z) over one oscillation of the periodic solution. Results are shown for the Sherman-Rinzel-Keizer model, for Zm z 4° and /3 = 2.7179e All other parameter values are as listed in Appendix A. Corresponding graphs for the remaining models under consideration in this thesis are qualitatively identical. Chapter 4. Determination of the Plateau Fraction 94 Pernarowski [67] and Pernarowski et al. [70] use a quadratic approximation of the form h2z + h1z + h0, (4.42) based on a least squares fit over 15 evenly spaced values of z, we choose to compute i(z) for 50 evenly spaced values of z between Zm and z, and use a simple linear interpolation to approximate i(z) for any other value of z. A linear approximation scheme is deemed sufficient based on the observation that the typical graph of (z) in Figure 4.6b, created with the same 50 data points, is smooth. In Figure 4.7 then, the numerical solution of the leading-order active phase problem, given by (2.51) and (4.37) with h(z) approximated as just described, is compared to the numerical solution of the full system of equations, (2.44)-(2.45). From Figure 4.7a, we observe that the two solutions are in excellent agreement over the first half of the active phase. As t increases, the oscillations become increasingly out of phase. We note that the oscillations are never out of phase by more than half a period for any of the models under consideration in this thesis. To obtain the leading-order active phase duration, we separate variables in (4.37) and integrate, yielding (0) T° = j e(z) • (4.43) The limits of integration are chosen as follows. According to the exact definition of the active phase duration in (4.4), the lower limit of integration should correspond as closely as possible to the value of z at the time when u reaches a local maximum during the first spike in the active phase. Since the solution trajectory leaves the silent phase near the local minimum of the G = 0 nulicline and the value of z does not change significantly during the transition from the silent to the active phase, we take z = Zm. Similarly, the upper limit of integration should correspond to the value of z at the time when u reaches a local maximum during the last spike in the active phase. Since z does not change Chapter 4. Determination of the Plateau Fraction 95 (a) o 100 200 300 400 500 600 t (b) -5 -5.2 - - - -5.3 - - -5.4 - - 55 I I I I I 0 100 200 300 400 500 600 t Figure 4.7: Comparison of the solution of the full system of equations (solid curves), (2.44)-(2.45), and the solution of the leading-order active phase problem (dashed curves), given by (4.33) and (4.35). Initial conditions are on the limit cycle solution of (2.51) at z = zm. Results are shown for the Chay-Kang model for /3 = 8.5301e2 All other parameter values are as listed in Appendix A. Corresponding graphs for the remaining models under consideration in this thesis are qualitatively the same. Chapter 4. Determination of the Plateau Fraction 96 significantly during the last spike, we take z = z. For an efficient computation of i(z) again should be approximated as described above. The exact and approximate active phase durations as a function of /3 are compared in Figure 4.8. As for the silent phase durations, the agreement is very good for all models. The decreases in the active phase duration as a function of /3 can be explained as follows. Differentiating (4.43) with respect to /9, we obtain dT° 1 zg eh(za = - 2dz < 0, (4.44)d/3 6 Zm [/9e_zh(z) — 1] since i(z) > 0. Alternatively, we recall that > 0 during the active phase. Therefore, from (2.39), an increase in 3 implies faster dynamics for z and, hence, a decrease in the active phase duration, as observed. 4.5 Leading-Order Approximation of the Plateau Fraction From (4.5), we define the leading-order plateau fraction, to be = (4.45) The exact and leading-order plateau fractions as a function of 3 are compared in Fig ure 4.9. For all models, the leading-order results overestimate the exact plateau frac tions. Since T(°) and T° are in excellent agreement with T3 and Ta, respectively, the discrepancy must be due to the omission of the duration of the transition phases in the denominator of (4.45). As observed in Figure 4.4, the duration of the transition from the silent to the active phase may be significant. The same is true for the transition from the active phase back to the silent phase (cf. Figure 4.7). The time spent in the transition phases, Tb — (T3 + Ta), as a percentage of the burst period, Tb, is shown in Figure 4.10. We observe that, in general, the transition phase Chapter 4. Determination of the Plateau Fraction 97 0 0.01 (d) 700 600 500 a 400 F—i 300 200 100 0 0.07 (f) 650 600 550 a 500 F—i 450 400 350 300 6.6 Figure 4.8: Comparison of the exact active phase duration (jagged curve), Ta, and the leading-order active phase duration (smooth curve), as a function of /3, for /3HC 3 /3m. All other parameter values are as listed in Appendix A. Results are shown for the (a) reduced Chay-Keizer model; (b) Chay (1986) model; (c) Himmel-Chay model; (d) Chay-Kang model; (e) Sherman-Rinzel-Keizer model; (f) Chay-Cook model. (b) 400 300 F— 200 100 (a) 240 210 180 E— 150 120 90 60 (c) 300 250 200 E— 150 100 50 0 15 20 25 30 35 40 /3 Si 0.0175 0.025 0.0325 /3 I I I I I 0.11 0.15 0.19 0.234 6 8 10 12 14 /3 (e) 300 250 200 a 150 100 50 0.02 0.025 0.03 0.035 0.04 /3 7 7.4 7.8 8.2 8.6 /3 Chapter 4. Determination of the Plateau Fraction 98 (a) 0.12 (b) :b 15 20 25 30 35 40 0.01 0.0175 0.025 0.0325 13 /3 (c) 0.35 i i i i (d) 0.4 _________ 4 6 8 10 12 14 0.07 0.11 0.15 0.19 0.23 1 1 (e) 0.45 ____ _ ____ _ ____ _ ° : 0.02 0.025 0.03 0.035 0.04 6.6 7 7.4 7.8 8.2 8.6 1 Figure 4.9: Comparison of the exact plateau fraction (jagged curve), p, and the lead ing-order plateau fraction (smooth curve), 1o), as a function of /3, for /3H0 j3 /3m. All other parameter values are as listed in Appendix A. Results are shown for the (a) reduced Chay-Keizer model; (b) Chay (1986) model; (c) Himmel-Chay model; (d) Chay-Kang model; (e) Sherman-Rinzel-Keizer model; (f) Chay-Cook model. Chapter 4. Determination of the Plateau Fraction 99 (a) 100 i i i i (b) 100 75 75 15 20 25 30 35 40 0.01 0.0175 0.025 0.0325 /3 (c) 100 i (d) 100 75 75 50 50 25 25 0 I I 0 I I 4 6 8 10 12 14 0.07 0.11 0.15 0.19 0.23 /3 (e) 100 (f) 100 75 - - 75 50- - 50 25 - - 25 - - 0 I I 0 I I 0.02 0.025 0.03 0.035 0.04 6.6 7 7.4 7.8 8.2 8.6 /3 3 Figure 4.10: Transition phase duration, Tb — (T8 + Ta), as a percentage of the burst period, Tb, for 13Hc 3 /3m• All other parameter values are as listed in Appendix A. Results are shown for the (a) reduced Chay-Keizer model; (b) Chay (1986) model; (c) Himmel-Chay model; (d) Chay-Kang model; (e) Sherman-Rinzel-Keizer model; (f) Chay-Cook model. Chapter 4. Determination of the Plateau Fraction 100 duration is small, but not insignificant. We explain the increase in the percentage with /3 as follows. For 3 near /3m there are significantly fewer spikes per burst (cf. Figure 4.3) and the active phase is much shorter (cf. Figure 4.8). More importantly though, the fixed point for the full system of equations, (2.44)-(2.45), lies close to the local minimum, (Urn, Zm), of the 0 nullcline. Therefore, the solution trajectory travels near this fixed point in its transition from the silent phase to the active phase, resulting in a significant increase in the duration of this transition phase. This is reflected in Figure 4.10. We note that the time spent in the two transition phases in roughly equal for all values of /3, except for values of /3 near /3m• Finally, we observe that the percentage of time spent in the transition phases is much higher for the Chay (1986) model than for the remaining models. We attribute this to the fact that the number of spikes per burst is much less than for the other models (cf. Figure 4.3), even for the smaller values of /3 near /3HC. We return to Figure 4.9 and observe that the plateau fraction decreases with /3. As for the silent and active phase durations, we can explain this trend by differentiating (4.45) with respect to /3 and using the signs derived in (4.17) and (4.44) to obtain d (°)Pf =_ <0, (4.46) as observed. 4.6 Examination of the Transition Phases We saw in the previous section that the transition phase durations are not insignificant, especially when 3 is near /3m, and it thus is desirable to obtain leading-order approxima tions of these durations. In this section, we follow the approach of Pernarowski [67] and Pernarowski et al. [70] to derive the equations of two surfaces which correctly describe each of the transition phases in (a, ii) space. However, as we will see below, the equations Chapter 4. Determination of the Plateau Fraction 101 cannot be used to derive a leading-order approximation of the transition phase durations. We return to the general nondimensional model, (2.37)-(2.39), and briefly summarize the argument used by Pernarowski et at. in their derivation of the surfaces just mentioned for the Sherman-Rinzel-Keizer model. During both transition phases, z changes little and can be regarded as constant, namely z Zm during the transition from the silent to the active phase (capture) and z during the transition from the active phase back to the silent phase (escape). Furthermore, they argue that w is near w(u, z) and that T(u, z) is 0(1) throughout both transition phases. Using the same argument for all models under consideration, the leading-order prob lem for the capture phase thus is = A(u; Zm) — w(u; Zm), (4.47) th = 0, (4.48) = 0. (4.49) The equations describing the escape phase are exactly the same, except Zm is replaced by z. Projecting the solution of the full system of equations, (2.44)-(2.45), onto the (u, ü) plane, and superimposing the curves defined by (4.47), as shown in Figure 4.11, confirms the accuracy of (4.47) in (u, it) space. However, (4.47) cannot be used to determine u as a function of t, nor can it be used to reduce the transition phase duration to quadrature, similar to what was done for the silent and active phase durations. The reason is that we wish to start and end the integrations on or very near the 0 = 0 nulicline. However, 0 = 0 is equivalent to A(u, z)—w(u, z) = 0 via (2.46)-(2.48) and (2.53). That is, (4.47) induces non-integrable singularities at the endpoints. In order to resolve this issue and successfully obtain the leading-order transition phase durations, we expect that it is necessary to perform a higher-order analysis, requiring Chapter 4. Determination of the Plateau Fraction 102 • 0 Figure 4.11: Projection of the numerical solution of the full system of equations, (2.44)-(2.45), onto the (u, i) plane (solid curve) with the leading-order approximations of the transition phases, (4.47), superimposed. For the upper dashed curve, corresponding to the escape, Z = Zm and for the lower dashed curve, corresponding to the capture, z = z. Results are shown for the reduced Chay-Keizer model. Corresponding graphs for the remaining models under consideration in this thesis are qualitatively identical. the method of matched asymptotic expansions [9] (where asymptotic solutions in the transition phases are matched to the asymptotic solutions in the silent and active phases to obtain a global, uniformly valid solution for the entire bursting phenomenon). Fur thermore, the slow passage through a limit point (in our case the local minimum of the G = 0 nulicline) has been analyzed for several problems, including bursting phe nomena [31, 40, 48]. We expect that a similar analysis would be useful in eliciting the dependence of the duration of the transition from the silent to the active phase on the parameter, and help explain the significant increase in this duration as 3 approaches /3m (cf. Section 4.5). However, we do not pursue these avenues of research further. U Chapter 5 Multiple Bifurcations in a Polynomial Analog Model In this chapter, we study bursting oscillations in a much broader context. To avoid the complex nonlinearities in the biophysical models, we concentrate on a simplified model exhibiting the bursting phenomenon, namely the polynomial model introduced by Pernarowski [68]. Pernarowski’s model is analogous to the biophysical models of bursting electrical activity in pancreatic beta cells, in the same way that the Fitzhugh-Nagumo model [34] of nerve membrane is analogous to the Hodgkin-Huxley model [47] describing the electrical impulses in the squid giant axon. Depending on the parameter values, Pernarowski’s model exhibits a wealth of oscilla tory behaviour, including square-wave bursting which is the type of bursting observed in pancreatic /3-cells. In [68], Pernarowski used a local stability analysis of the fixed points of the fast subsystem, a Hopf bifurcation analysis, and Melnikov’s method in connec tion with homoclinic bifurcations to determine the region in parameter space where the model equations exhibit square-wave bursting oscillations. We ask the question: Given any set of parameter values, can we predict the bifurcation structure of the fast subsystem and, hence, anticipate what kind of oscillatory behaviour the full system of equations will produce? In answering this question, we create a bifurcation map in parameter space, consisting of regions within which the bifurcation diagrams of the fast subsystem are qualitatively the same. The bifurcation map is significant in determining the relationship between the various types of oscillatory behaviour and, hence, provides a basis for an extension of the 103 Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 104 classification schemes of bursting oscillations by Rinzel [77] and Bertram et al. [10]. In Section 5.1, we introduce Pernarowski’s model and various types of solution be haviour. Some preliminary analysis is contained in Section 5.2, where we use the sym metry of the model equations to simplify the bifurcation map and in Section 5.3, where we perform a standard fixed point analysis of the fast subsystem. The details of the development of the bifurcation map are presented in Sections 5.4— 5.7. We develop analytical constraints on the various bifurcations exhibited by the fast subsystem of the polynomial model. The constraints correspond to curves on the bifur cation map bounding the different regions. Hopf bifurcation constraints are developed in Section 5.4 and homoclinic bifurcation constraints are developed in Section 5.5. Con straints associated with another type of bifurcation exhibited by the polynomial analog model, namely the saddle-node of periodics bifurcation, introduced below, are investi gated in Section 5.6. In that section, we need to resort to numerical software to complete the bifurcation map. The complete bifurcation map is presented in Section 5.7. We note that each type of bifurcation diagram predicted by the analysis is shown in the figures in Appendix B. In the development of the bifurcation map, which is very qualitative in nature, we lose some pertinent quantitative information. To recover some of this information, we examine two-parameter bifurcation diagrams in Section 5.8. A second reason for presenting this type of diagram is because of the surprising discovery that it is essentially identical to the two-parameter bifurcation diagram generated for the biophysically based Chay Cook model in [10]. That is, the polynomial analog model not only superficially mimics various types of solution behaviour of models of excitable cells, but also possesses the rich bifurcation structure of such models. We discuss the implications of this discovery for the classification of bursting oscillations in Section 5.9. Finally, knowledge of the bifurcation structure of the fast subsystem, in general, Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 105 is sufficient for predicting the type of oscillatory behaviour of the solution of the full system of equations. However, this is not true near the boundaries of the regions on the bifurcation map. The dependence of the type of oscillatory behaviour exhibited by the full system of equations on the size of the asymptotic parameter e governing the dynamics of the slow subsystem is examined briefly in Section 5.10. 5.1 Pernarowski’s Analog Model and Various Solution Behaviours Pernarowski’s analog model [68] is given by + F(u)ü + G(u, c) = —eH(u, c), (5.1) é = eH(u,c), (5.2) where F(u) a [(u — — 2] , (5.3) G(u,c) = c+u3—3(u+1), (5.4) H(u, c) = 3(’u — n) — c. (5.5) We note the similarity of (5.1)-(5.2) to the transformed biophysical models of bursting electrical activity in pancreatic /3-cells in Chapter 2 (cf. (2.44)-(2.45)). Based on this similarity, the variable u may be interpreted as the membrane potential of an excitable cell. Since the parameter e is small (0 < e << 1), (5.2) models the dynamics of a slowly changing quantity. In the case of pancreatic /3-cells, c may be interpreted as the intracellular calcium concentration or some other slow variable. As in Chapter 3, we divide (5.1)-(5.2) into slow and fast subsystems. The slow subsystem is given by (5.2) and the fast subsystem by ii + F(u)ii + G(u; c) = 0, (5.6) Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 106 6 5 4 c) 3 2 1 0 U Figure 5.1: The G 0 nulicline of the fast subsystem, (5.6). The local minimum is at (Urn, Cm) = (—1, 1) and the local maximum is at (UM, CM) = (1,5). Fixed points lying between u = —2 and u = 2 are in the region of bistability. The C = 0 nullcline is a reflection of itself through the point (0, 3) (cf. Section 5.2). where the use of the semicolon in G indicates that c is treated as the bifurcation pa rameter. The fast subsystem nullcline is given by G = 0, which traces out a cubic in (u, c) space, as shown in Figure 5.1. Since G does not depend on any model parameters, the nullcline and its local extrema remain fixed. We introduce (Urn, Cm) = (1, 1) and (UM, CM) (1, 5) to denote the local minimum and maximum of the nulicline, respec tively. We refer to the region lying between Cm and CM as the region of bistability. For most choices of the parameter values, the fast subsystem indeed exhibits bistability in this region, although parameter values can be found for which this is not the case. The value of U defining the region of bistability on the right branch of the G = 0 nullcline is given by the largest root of G(U, Cm) = 0, which is U = 2. Similarly, the value of u defining the region of bistability on the left branch is given by the smallest root of G(U, CM) = 0, which is U = —2, as indicated in Figure 5.1. The stability of the fixed points of the fast subsystem and, hence, the behaviour of the -3 -2 -1 0 1 2 3 Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 107 solution of the full system of equations, (5.1)-(5.2), is determined by the damping term, F(u)zt. Hence, depending on the values of the parameters in F, the model exhibits a wide variety of oscillatory behaviour. Numerical solutions of (5.1)-(5.2) and the bifurcation diagrams of (5.6) for three sets of parameter values are shown in Figures 5.2-5.4. Figure 5.2 shows square-wave bursting similar to that observed in the electrical activity of pancreatic ,6-cells. The bifurcation diagram of this burster exhibits a supercritical Hopf bifurcation (RB) on the right branch of the C = 0 nulicline. Periodic orbits emanating from the Hopf point are stable and increase in size as c increases, until they disappear at a homoclinic bifurcation (RC), as in Chapter 2. Tapered bursting is shown in Figure 5.3 and is similar to the electrical activity observed in the pyramidal cells of the cat hippocampus [50]. Compared to the bifurcation diagram of the square-wave burster, a second supercritical Ropf bifurcation has appeared on the right branch of the G = 0 nullcline. As c increases, the periodic orbits emanating from the lower Hopf point (not shown) first increase in size and then decrease. The periodic orbits disappear through the upper Hopf bifurcation, rather than through a homoclinic bifurcation. A third type of bursting is shown in Figure 5.4. In this case, the individual spikes do not ride on a plateau, as for the previous bursters, but undershoot the value of u in the silent phase. This type of bursting was named nearly-parabolic bursting by Pernarowski [68]. However, the name is slightly misleading, as parabolic bursting is a different type of bursting, requiring at least two slow variables [73, 78]. Instead, we refer to this type of bursting as bursting with an undershoot. As for the square-wave burster in Figure 5.2, the bifurcation diagram for this burster also exhibits one Hopf bifurcation on the right branch of the G = 0 nullcline. But in this case, the Ropf bifurcation is subcritical, rather than supercritical. Periodic orbits emanating from the Hopf point are unstable at first, but become stable at a saddle-node of periodics bifurcation (SNP). In 108Chapter 5. Multiple Bifurcations in a Polynomial Analog Model (a) III’ (b) 0 50 100 150 200 250 t I I I I SN ci 6— 3- 0- -3 I I I -3 -2 -1 0 1 2 3 4 U Figure 5.2: Square-wave bursting, as observed in pancreatic /3-cells ( = 0.7, = 1.6, a = 0.25, 3 = 4, ü = —0.954, and E = 0.0025). (a) Numerical solution of the full system of equations, (5.1)-(5.2). (b) Bifurcation diagram of the fast subsystem, (5.6), with the projections of the solution from (a) superimposed. As in Chapter 2, solid lines indicate stability, dashed lines indicate instability, SN denotes a saddle-node bifurcation, RB denotes a Hopf bifurcation, and RC denotes a homoclinic bifurcation. Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 109 (a) 2 1 0 -1 -2 0 100 200 300 400 500 600 t (b) 6 5 4 c) 3 2 1 0 -3 -2 -1 0 1 4 U Figure 5.3: Tapered bursting (q = 0.7, i = 2.1, and remaining parameter values as in Figure 5.2). (a) Numerical solution of the full system of equations, (5.1)-(5.2). (b) Bifurcation diagram of the fast subsystem, (5.6), with the projection of the solution from (a) superimposed. Not shown is a supercritical Hopf bifurcation on the right branch of the C = 0 nullcline at (u, c) (2.8, —10.55). Ii,,.. sj% I I SN I I / I’ 2 3 Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 110 (a) i 0 50 100 150 200 t (b) 2 HC HC1- SN c) 0- -1 - L. SNP SNP -2 I I I I I -3 -2 -1 0 1 2 3 4 U Figure 5.4: Bursting with an undershoot ( = 1.2, i 1.0, and remaining parameter values as in Figure 5.2). (a) Numerical solution of the full system of equations, (5.1)-(5.2). (b) Bifurcation diagram of the fast subsystem, (5.6), with the projection of the solution from (a) superimposed. Unstable periodic orbits become stable at the saddle-node of periodics (SNP). Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 111 addition, the periodic orbits grow very quickly in size, and surround all three fixed points when c> Cm. The periodic orbits disappear through a homoclinic bifurcation. We study these large periodic orbits and corresponding homoclinic bifurcations in more detail in the following sections. At this point, it suffices to mention that the undershoot observed in the numerical solution of (5.1)-(5.2) is simply the result of the fact that in the region of bistability, the fast subsystem exhibits periodic orbits with minimum u values smaller than the u values on the left branch of the G = 0 nullcline. We note that the values of all parameters except and i% were kept constant to produce the three different types of oscillatory behaviour in Figures 5.2-5.4. Depending on the values of these parameters, the bifurcation diagram of the fast subsystem may or may not exhibit Hopf, homoclinic, and saddle-node of periodics bifurcations. Furthermore, Hopf bifurcations may be supercritical or subcritical, and homoclinic bifurcations may correspond to a homoclinic orbit surrounding the fixed point on the right branch of the G = 0 nullcline only, or to a homoclinic orbit surrounding all three fixed points, and so on. That is, as the parameters are varied, the bifurcation diagram itself is undergoing bifurcations. It is these bifurcations that are of interest in this chapter. The flopf, homoclinic, saddle-node, and saddle-node of periodics bifurcations exhib ited in Figures 5.2-5.4 are classified as codimension-1 bifurcations. The codimension of a bifurcation is the smallest dimension of a parameter space which contains the bifurcation in a consistent way [37]. More precisely then, the codimension-l bifurcations themselves are undergoing bi furcations as and % are varied. For example, supercritical Hopf bifurcations become subcritical at a codimension-2 degenerate Hopf bifurcation (DHB). Similarly, stable ho moclinic orbits become unstable at a codimension-2 neutral saddle loop bifurcation (NSL). The codimension-2 bifurcations associated with the Hopf, homoclinic, and saddle-node of periodics bifurcations are studied in Sections 5.4-5.6. Each codimension-2 bifurcation Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 112 gives rise to a curve in (, i%) space. These curves in turn bound regions within which the bifurcation diagrams of the fast subsystem, (5.6), are qualitatively the same. Before studying the codimension-2 bifurcations, we perform preliminary symmetry and fixed point analyses in the next two sections. 5.2 Symmetry Considerations In this section, we briefly examine the symmetries contained in the polynomial model. We use the symmetry in F to reduce the bifurcation map to the first quadrant of the (, ) plane. The symmetries will be explored further in Sections 5.4-5.6. From (5.4), we deduce that G(u, c) = —G(—u, 6 — c), (5.7) so that the G = 0 nullcline is a reflection of itself through the point (u, c) = (0, 3), as indicated in Figure 5.1. From (5.3), we deduce that F(u;,q) = F(u;,—i), (5.8) F(u; ii, ) = F(—u; —ii, ). (5.9) Since the parameter i appears in F only, it is clear from (5.8) that changing to — has no effect on (5.1)-(5.2). Therefore, we eliminate the left half-plane of the bifurcation map and consider values of 0 only. Even though changing u to —u and ft to —ft does not affect F (cf. (5.9)), this transformation does affect the bifurcation diagram, in that the bifurcations occurring on the right branch of the G = 0 nullcline now occur on the left branch and vice versa. Curves of codimension-2 bifurcations thus originate in both the upper and lower half-planes. In fact, each curve in the first quadrant of the bifurcation map is reflected across ft = 0. As such, we only show the first quadrant of the bifurcation Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 113 map. However, since curves originating in the fourth quadrant may extend into the first quadrant, we consider both positive and negative values of % during the derivation of the equations. 5.3 Preliminary Dynamical Systems Analysis We apply standard dynamical systems methods (see for example [37, 46]) to the fast subsystem, (5.6), to determine the stability of the fixed points on the C = 0 nullcline. Even though the analysis is elementary, we include it here as a prerequisite for the following sections and derive some quantities that are referred to later. We note again that the location of the fixed points is independent of the parameter values in this simple polynomial model. However, the stability of the fixed points on the C 0 nullcline is affected by the model parameters in the damping term, determined by F. Writing (5.6) as a system of two first-order differential equations by letting = (x1,2)= (u,ü), we obtain (5.10) where f1(;c)f(;c) = I I = I I . (5.11)\ f(; c) ) \ —G(xi; c) — F(xi)x2 ) Fixed points of (5.10) are of the form i = (1,0), where is a root of G(i;c) = 0. The Jacobian matrix of (5.11) at the fixed point = i is / 0 1 DJ(;c) = 1 V2 = I, (5.12) I \ —G(i;c) —F(1) ) with determinant and trace given by DetDf(; c) = G(i; c) = 3( — 1), (5.13) TrDJ(;c) = —F(1), (5.14) Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 114 respectively. The stability of the fixed points is deduced from the eigenvalues of Df(; c), given by TrDf ± Tr2Dj_ 4DetD (5.15) If the determinant is negative, the eigenvalues i1,2 are real and opposite in sign, and hence the fixed point is a saddle. Along the middle branch of the G = 0 nulicline, is between m = —1 and tiM = 1 (cf. Figure 5.1), so that the determinant is negative from (5.13). Thus, all fixed points along the middle branch are saddles. Similarly, the determinant is positive along the left and right branches of the G = 0 nulicline. In this case, the fixed points are either nodes or spirals, depending on the relative magnitude of the trace and the determinant. If the discriminant in (5.15) is positive, the fixed points are nodes. Otherwise the fixed points are spirals. Their stability in turn is determined by the sign of the trace. Stability is obtained if the trace is negative, and instability if the trace is positive. 5.4 Constraints on the Hopf Bifurcations In this section, we examine the codimension-2 bifurcations associated with Hopf bifur cations. Hopf bifurcations occur when stable spirals become unstable spirals and vice versa. Based on the analysis in the previous section, Hopf bifurcations thus occur when the determinant in (5.13) is positive and the trace in (5.14) is zero. Written in term of ti and c, we require G(u; c) = 0, (5.16) G(u;c)=3(u2—1) > 0, (5.17) F(u) = 0. (5.18) Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 115 The first condition simply states that there needs to be a fixed point. The remaining conditions together ensure that the fixed point is indeed a Hopf point. In particular, (5.17) implies that the Hopf bifurcation(s) may occur anywhere on the left or right branch of the G = 0 nulicline, since u < ‘Urn —1 on the left branch and u > UM 1 on the right branch. Finally, since F is quadratic in u (cf. equation (5.3)), (5.18) implies at most two Hopf bifurcations. The symmetry about i = 0 introduced in Section 5.2 is explored further here. Using (5.7)-(5.9) in (5.16) and (5.18) implies G(—u; 6 — c) = 0, (5.19) F(—u; —ii, ) 0, (5.20) so that if there is a Hopf bifurcation at (u, c) = (UHB, CHB), then there is also a Hopf bifurcation at (u, c) = (—uHB, 6 — CHB) when the sign of is reversed. As a result, we expect codimension-2 curves associated with Hopf bifurcations on the right branch of the G = 0 nulicline to have a reflection, representing the codimension-2 curves associated with the corresponding Hopf bifurcation on the left branch, about = 0. From now on, we refer to Hopf bifurcations on the right branch of the 0 = 0 nulicline as right Hopf bifurcations. Similarly, Hopf bifurcations on the left branch of the G 0 nulicline are referred to as left Hopf bifurcations. We examine the codimension-2 bifurcations affecting the number of Hopf bifurcations in more detail in Section 5.4.1 and investigate the codimension-2 bifurcations affecting the location of the corresponding Hopf points relative to the local extrema of the G = 0 nulicline in Section 5.4.2. In Section 5.4.3, we are concerned with the stability of the periodic orbits emanating from the Hopf points. We note that some of this work has been done previously by Pernarowski in [68]. We extend his work here by allowing left Hopf bifurcations and by allowing Hopf points to coalesce. Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 116 5.4.1 Number of Hopf Bifurcations In Figures 5.2-5.4 we observed that the fast subsystem, (5.6), can exhibit one or two right Hopf bifurcations, depending on the values of i and i. For certain values of the parameters, it is also possible that there are no Hopf bifurcations. To determine the number of Hopf bifurcations, we examine F(u) in more detail. From (5.3), we note that the graph of F is a parabola (opening upward when a is positive) and always has two roots, except when = 0, in which case there is a double root at u = %. We discuss this special case a little later. For positive values of , the left and right roots are given by UL = — , (5.21) UR = ‘ + , (5.22) respectively. We are interested in the location of the two roots relative to the u values at the local extrema of the G = 0 nullcline, Urn = —1 and UM = 1 (cf. Figure 5.1), and wish to determine the effect of changing and %. i ii and ‘2 are small, both roots lie between Urn and UM. That is, F is zero along the middle branch of the G = 0 nulicline, where all fixed points are saddles. In this case, there can be no Hopf bifurcations. Varying the value of causes a horizontal shift in the graph and roots of F(U). Increasing t causes the roots to be shifted to the right. When UR is to the right of UM, but UL is still between Urn and UM, F has one zero on the right branch and, thus, there is one corresponding right Hopf bifurcation. When it is increased further, so that UL is to the right of UM as well, F has two roots on the right branch, corresponding to two right Hopf bifurcations. Similarly, left Hopf bifurcations emerge when the roots are shifted to the left of Urn by a decrease in it. Varying the value of causes a vertical shift in the graph of F(U). Increasing i causes a vertical downward shift, resulting in an increase in the distance between the left and Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 117 right roots of F. For q large enough, it is thus possible that UL lies to the left of Urn and UR to the right of uM, corresponding to one right and one left Hopf bifurcation. The four curves representing the emergence of Hopf bifurcations at the local extrema (saddle-nodes) of the G = 0 nulicline, corresponding to UR UM, UL = UM, UR = Urn, and UL = Urn, respectively, are ft = 1+i, (5.23) ft = 1—a, (5.24) ft = —1+, (5.25) ft —1 — . (5.26) These codimension-2 bifurcations are known as Takens-Bogdanov bifurcations (TB) (see [35, 36, 38J and references therein). We note that the curves defined by (5.23) and (5.26) are reflections of each other across ft = 0, as are the curves defined by (5.24)-(5.25), as expected from the symmetry discussion earlier. The four curves are referred to as the minor TB (right), major TB (right), major TB (left), and minor TB (right) curves, respectively. The right and left descriptions refer to whether the emergence of the Hopf bifurcation occurs on the right or left branch of the G = 0 nulicline, and the major and minor descriptions refer to the first (major) or second (minor) Takens-Bogdanov bifurcation on the respective branch. The minor TB (right) curve and the portions of the major TB (right) and major TB (left) curves lying in the first quadrant of the bifurcation map are shown in Figure 5.5, and the number of Hopf bifurcations on the right and left branches of the G = 0 nullcline is indicated in each region defined by the curves. In addition to the Takens-Bogdanov bifurcation, there is a second codimension-2 bifurcation through which Hopf bifurcations can (dis)appear, due to the double root of F, which corresponds to the coalescence of two Hopf points. From (5.3), F(U) has a Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 118 Figure 5.5: Curves of Takens-Bogdanov bifurcations, minor TB (right), major TB (right), and major TB (left), representing the emergence of Hopf bifurcations at the local extrema of the G = 0 nulicline, and the CHB (right) curve, representing the coalescence of two Hopf points on the right branch of G = 0. double root u = f when = 0. From (5.17), we require ii < —1 or i > 1 for Hopf bifurcations to exist. The curve given by q = 0 and ii < —1 is referred to as the CHB (right) curve, and is shown in Figure 5.5. The curve given by i = 0 and i > 1 is the reflection of the CHB (right) curve across = 0, and is referred to as the CHB (left) curve (not shown). 5.4.2 Location of the Hopf Points In this sectioll, we are concerned with the location of the Hopf points relative to the c values at the local extrema of the G(u, c) = 0 nulicline, Cm 1 and CM 5 (cf. Figure 5.1). The vertical location of the Hopf points can dramatically affect the solution behaviour, as shown in Figure 5.6. Figure 5.6b shows a bifurcation diagram similar to 3 2 1 0 0 0.5 1 1.5 2 2.5 77 Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 119 0 100 200 300 400 500 -3 -2 -1 0 1 2 3 U Figure 5.6: Parabolic-amplitude bursting (q = 0.47, G = 1.3, and remaining parameter values as in Figure 5.2). (a) Numerical solution of the full system of equations, (5.1)-(5.2). (b) Bifurcation diagram of the fast subsystem, (5.6). (a) 3.5 2.5 1.5 0.5 -0.5 -1.5 -2.5 (b) 6 5 4 03 2 1 0 : — — — Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 120 the one shown in Figure 5.2b, but in which the Hopf point has moved up along the right branch of the G = 0 nulicline so that it lies in the region of bistability, between Cm and CM (cf. Figure 5.1). The corresponding oscillatory behaviour produced by (5.1)-(5.2) is shown in Figure 5.6a. Rather than a square-wave burster as in Figure 5.2, we now have a burster where the amplitude of the spikes during the active phase changes dramatically. Pernarowski [68] refers to this type of bursting as parabolic amplitude bursting. We distinguish between bifurcation diagrams such as shown in Figures 5.2b and 5.6b by imposing the constraint that a Hopf point lies either inside or outside the region of bistability. From Figure 5.1, there is a bifurcation in the location of the Hopf point when the corresponding root of F occurs at u = —2 or u = 2. The curves corresponding to UR = 2, uL = 2, UR = —2, and UL = —2 are given by the following equations, = 2+, (5.27) 2—n, (5.28) ‘ii = —2 + , (5.29) ft = —2 — ,, (5.30) and are referred to as the minor HBL (right), major HBL (right), major HBL (left), and minor HBL (left) curves, respectively. Moving across the major HBL (right) curve on the bifurcation map corresponds to moving the major Hopf point on the right branch of the G = 0 nulicline into or out of the region of bistability in the bifurcation diagram of the fast subsystem, and so on. The curves are superimposed onto Figure 5.5 in Figure 5.7, and the locations of the Hopf points (if there are any) are indicated in each region. As in the previous section, the symmetry about ft = 0 is evident in the pair of curves defined by (5.27) and (5.30), and the pair defined by (5.28) and (5.29). Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 121 3 I nor - - HBL (right) - - - - - 2 ,mjor right RB below Cm, - - pTihor right RB between Cm and CM, - - left RB between Cm and CM -, major - - - - 1 ( HBL (right) righffrnajo HB’s between - below Cm, HBL 0 above CM (left) Figure 5.7: As Figure 5.5, with the minor HBL (right), major HBL (right), and major HBL (left) curves, representing bifurcations in the location of the respective Hopf point (relative to the c values at the local extrema of the G = 0 nullcline), superimposed. 5.4.3 Criticality of the Hopf Bifurcations In Figures 5.2-5.4, we observed that Hopf bifurcations can be either supercritical or subcritical. That is, the periodic orbits emanating from a Hopf point can be either sta ble or unstable, respectively. In this section, we compute the curves which represent codimension-2 degenerate Hopf bifurcations (DHB), distinguishing regions on the bifur cation map in which the bifurcation diagrams exhibit either supercritical or subcritical Hopf bifurcations. We first apply the center manifold theorem (see for example [37]) to (5.10). We assume that (5.10) has a Hopf point at = (XHB, 0) and c = CHB. Letting = (yi, y2) = (x1 — XHB, x2) and writing the Taylor expansions for G(xi; c) and F(xi) about u = = XHB Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 122 and c = CHB, we obtain dil (5.31) where 0 1 B = , (5.32) Gu(XHB;CHB) 0 = (iW) (5.33) 92(Y)) gi() = 0, (5.34) g2(y) = Guuyi2— Guuuyi3— Fuy1y2 Fuuyi2y, (5.35) and the bar over the G and F denotes evaluation at u XHB, c cHB. The eigenvalues of B are +i, where w = iJ. The eigenvector corresponding to the wi eigenvalue is (1\ fo I I+I I. (5.36) 0) \w) We use the two components of this eigenvector as two new basis vectors, and make the transformation yi(l,O) +y2(O,l) = zi(0,w) +z2(1,0). (5.37) Letting = (zi,z2) we write (5.31) in the standard form (2 0 —wI 1= -‘-I I (5.38) \\Z2) w 0 where h1() --g2W) = — — puz1z2 — (5.39) = giQ7) = 0. (5.40) Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 123 Using results in [37], the stability of the periodic orbits emanating from the Hopf point is determined by the sign of aHB, where 1 03h1 93h1 1 82h1 102h 82h1\ aHB = jj €9z1 + 9z1ô2 + ii-:; ãz12 ‘ -a;-} (5.41) and all derivatives of h1 are evaluated at (0, 0). Using (5.3)-(5.4), XHB = 41 ± 27 from (5.21)-(5.22), and (5.39)-(5.40) in (5.41), we obtain aHB = 1 2 [uuu — UUU]16w = 6a[x—2xH41+1] = 6a [272 j2 + 1]. (5.42) If <0, the Hopf bifurcation is supercritical and if aHB > 0, the Hopf bifurcation is subcritical. We are interested in the degenerate case, aHB = 0. Setting aHB 0 in (5.42), we obtain the two curves representing degenerate Hopf bifurcations in (ii, 41) space, 41 + 772, (5.43) again reflecting the symmetry about 41 = 0. The curves are referred to as the DHB (right) and DHB (left) curves, respectively. The DHB (right) curve is plotted in Figure 5.8. The DHB (left) curve is the reflection of the DHB( right) curve across 41 = 0, and lies entirely in the fourth quadrant. In the region lying between the two DHB curves, all Hopf bifurcations are subcritical. In the regions lying outside these curves, all Hopf bifurcations are supercritical. It is interesting to note that the minor TB curves (cf. Figure 5.5), representing the emergence of a second right or left Hopf bifurcation, lie in the regions of supercriticality. That is, for bifurcation diagrams exhibiting two Hopf bifurcations on one branch, the periodic orbits emanating from the corresponding Hopf points are always stable. Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 124 Figure 5.8: As Figure 5.7, with the DHB (right) curve, representing a degenerate right Hopf bifurcation, superimposed. 5.5 Constraints on the Homoclinic Bifurcations In this section, we examine the codimension-2 bifurcations associated with homoclinic bifurcatjons. We note that homoclinic bifurcations are sometimes referred to as saddle loops (SL) [36, 38], and some of the terminology used later will reflect this. As in Chapter 3, the analysis of homoclinic bifurcations is based on Melnikov’s method [64]. In Chapter 3, we were only interested in a homoclinic orbit surround ing fixed points on the right branch of the G = 0 nullcline. In this chapter, we also are interested in homoclinic orbits surrounding fixed points on the left branch of the G = 0 nulicline, as well as large homoclinic orbits surrounding three fixed points. From now on, we refer to the homoclinic orbits surrounding fixed points on the right branch as right homoclinic orbits. Similarly, homoclinic orbits surrounding fixed points on the left branch are referred to as left homoclinic orbits and homoclinic orbits surrounding 3 2 1 0 0 0.5 1 1.5 2 2.5 Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 125 three fixed points are referred to as large homoclinic orbits. Large homoclinic orbits can be associated with a Hopf bifurcation on the right or left branch of the C = 0 nulicline. When we need to distinguish between these, we refer to large right or large left homoclinic orbits, respectively. In the same vein, we refer to right, left, large (right), and large (left) periodic orbits. The existence of a homoclinic orbit is inferred from a change in the relative position of the stable and unstable orbits of the saddle point in the (u, zt) phase plane as the bifurcation parameter c is varied. The situation for a right homoclinic orbit is as shown in Figure 3.3. For a low value of c the right stable orbit lies outside the right unstable orbit and there is a right periodic orbit. For a high value of c the stable orbit lies inside the unstable orbit and there are no periodic orbits. There is a right homoclinic orbit at the value of c when the relative position of the orbits reverses and the stable and unstable orbits coincide. The situation for a left homoclinic orbit is similar, but in this case the relative position of the left stable and unstable orbits reverses. The situation for a large homoclinic orbit is shown in Figure 5.9. In this case, the orbits of interest are the large ones, connecting to the saddle point in the upper half plane. For a low value of c the large stable orbit lies inside the large unstable orbit and there is a large periodic orbit. For a high value of c the stable orbit lies outside the unstable orbit and there are no periodic orbits. The large periodic orbits cease to exist when the large stable and unstable orbits coincide, creating a large homoclinic orbit connecting the saddle point to itself. The Melnikov distance functions, representing the distance between the stable and unstable orbits at a given value of c, are derived in Section 5.5.1. Since the definitions of F(u) and G(u, c) are simple, we can determine the Melnikov distance functions analyti cally. We use the Melnikov distance functions to derive the equations of the codimension-2 curves in the remaining sections. In Section 5.5.2, we distinguish between regions on the bifurcation map where bifurcation diagrams exhibit right, left, or large periodic orbits. Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 126 (a) 6 4 2 0 -2 -4 -6 -8 -4 0 4 (b) 6 4 2 0 . -2 -4 -6 -8 -4 4 U Figure 5.9: Phase portraits of the fast subsystem, (5.6), calculated with DSTOOL, show ing the stable (solid curves) and unstable orbits (dashed curves) of the saddle point and illustrating the existence of a large homoclinic orbit. (a) For c 2.8, the large stable orbit lies inside the large unstable orbit and there is a large stable periodic orbit sur rounding all three fixed points. (b) For c = 3.0, the large stable orbit lies outside the large unstable orbit and there are no periodic orbits. The remaining parameter values are = 2.0, i = 1.1, and a 0.25. I I -2 2 U -2 0 2 Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 127 In Section 5.5.3, we determine the conditions required to obtain a double homoclinic orbit, that is, a right and left homoclinic orbit at the same time. In Section 5.5.4, we are concerned with the stability of the periodic orbits near homoclinic bifurcations. Fi nally, we examine the conditions required for two homoclinic bifurcations to coalesce and disappear in Section 5.5.5. 5.5.1 The Melnikov Distance Functions In this section, we derive the Melnikov distance functions measuring the separation be tween the stable and unstable orbits of the saddle point in the (u, it) phase plane of the fast subsystem, (5.6), at a given value of c between Cm and CM. The derivation is essen tially identical to the derivation outlined in Chapter 3, and we omit most of the details here. In [68], Pernarowski showed that the value of F(u)it is numerically small relative to ii and G(u, c) in the active phase near the homoclinic bifurcation. We thus can treat the F(u)ii term as a perturbation to the Hamiltonian system ü + G(u; c) = 0, (5.44) so that Melnikov’s method can be applied to the fast subsystem, (5.6). As in the case of the biophysical models for bursting electrical activity in pancreatic /3-cells, the un perturbed system has three fixed points for every value of C between Cm and CM, and its phase portrait at a typical value of c is as in Figure 3.8. The saddle point occurs at (u, it) = (a8, 0), and the separatrix crosses the u-axis at (b8, 0) and (C3, 0), where the values of a3, b8, and c. depend on C. From Chapter 3, the Melnikov distance function is given by D(c) = J F(u(t))it(t)2d , (5.45) Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 128 where (u, it) needs to evaluated along the separatrix solution of (5.44). To convert the integral in (5.45) to a line integral in u, we first integrate (5.44) along the separatrix from a3 to u, and use the fact that it = 0 when u = a3, to give it = +/—2V(u;a3(c)), (5.46) where V(u;a3(c)) = J G(u; c)dü, (5.47) as(c) as before. However, whereas in Chapter 3 we had to resort to a numerical evaluation of V, the simple definition of G allows us to determine V exactly. Using (5.4) in (5.47), we obtain V(u;a3(c)) = (u — a3) [4c + (u2 + a)(u + a3) — 6(u + a3) — 12]. (5.48) Even though a3 can be expressed as a function of c explicitly, for convenience, we choose to express c as a function of a3 instead. Since a3 is given by G(a3,c) = 0, we have c 3(a + 1) — a from (5.4). Then V(u; a3) = (u — a3) [3 + a3u2 + (a — 6)u — 3a — 6a3] = (u — a3)2 {u2 + 2a3u + 3a — 6]. (5.49) Since it = 0 when u = b3 and u = C3 as well, it follows from (5.46) that b3 and c3 are given by V(b3;a3) = 0 and V(c3;a3) = 0, respectively. It is clear from (5.49) that V(u; a3) has a double root at a3. The remaining two roots give b3 and c3. Requiring C3 < a3 < b3, we obtain = —a3 + ./6 — 2a, (5.50) c3 = —a3 — iJ6 — 2a. (5.51) Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 129 We now have all the information to convert the integral in (5.45) to a line integral in u. We obtain three different Melnikov distance functions, DR, DL, and Ds, depending on whether we are interested in right, left, or large homoclinic orbits, respectively. When we wish to examine a right homoclinic orbit, we need to evaluate the integral in (5.45) along the right loop of the separatrix solution of (5.44) (cf. Figure 3.8). Using (5.3), (5.46), and (5.49) in (5.45) gives the right Melnikov distance function, bs DR(a3) 2f F(u)—2V(u;a8)du b I 1 = 2 j F(u)V_(u — as)2 (u2 + 2a3u+ 3a — 6)du = F(u)lu —a3I— (u2 + 2a8u + 3a — 6)du (5.52) = ajb3 [(u — )2 — 2] (u —a3)— (u2 + 2a3u + 3a — 6)du. (5.53) Using standard tables of integrals and (5.50) in (5.53), we find DR(a3)= a [eR2(as)(i2 — 2) + eR1(a5)i+ eRo(as)j, (5.54) where eRo(as) = —v’(a — 2a — 4)1 — a + 6V’a3(a — 3)(a8), (5.55) eRl(a3) = 6/a3(3 — a)1 — a + 3(a — 3)(a + 1)(a3), (5.56) eR2(as) = 4v”/1 — a + 2s../a3(a — 3)(a), (5.57) = ( 2a3 (5.58) — 2a) Similarly, when we wish to study a left homoclinic orbit, we need to evaluate the integral in (5.45) along the left loop of the separatrix solution of (5.44) (cf. Figure 3.8). Using (5.3), (5.46), and (5.49) in (5.45) then gives the left Melnikov distance function, DL(a3)= —a x:s [(u — )2 — 2j (u —a3)— (u2 + 2a8u + 3a — 6)du. (5.59) Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 130 The minus sign is due to the fact that Iu — a3 = —(u — a3) along the left ioop of the separatrix (cf. equation (5.52)). Using (5.51) in (5.59), we find DL(a3)= a [eL2(ag)(112 — 2) + eLl(a3)u + eLo(as)], (5.60) where eLo(a8) = —/(a — 2a — 4)1 — a — 6v’a3(a — 3)(a), (5.61) eLl(a8) = 6’a3(3 — a)i/l — a — 3V’(a — 3)(a + 1)(a3), (5.62) eL2(a8) = 4//1 — a — 2/a3(a — 3)(a), (5.63) (a8) = cos’ (_2as). (5.64) Finally, when we are interested in a large homoclinic orbit, we need to evaluate the integral in (5.45) along both loops of the separatrix solution of (5.44) (cf. Figure 3.8). The corresponding large Melnikov distance function, D, is simply the sum of DR and DL. That is, Ds(a3)= a [es2(as)(2 — 2) + esi(a3)+ eso(a3)J, (5.65) where eso(a3) = eRo(a8)+ eLo(a8), (5.66) esi(as) = eRl(a8)+ eLl(as), (5.67) es2(a5)= eR2(a8)+ eL2(a8). (5.68) We now return to the discussion of the symmetry induced by changing % to —ü. From (5.55)-(5.58) and (5.61)-(5.64), we find eRo(a8) = eLo(—aS), eRl(as) = —eLl(—a3), eR2(a8) = eL2(—as). Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 131 Then (5.54) and (5.60) imply DR(a; 1%, ) = DL(—a8;—, q). (5.69) From (5.7), we have a3(c) = —a3(6 — c), (5.70) which together with (5.69) implies that if there is a right homoclinic orbit at (uHU, CHC), then there must be a left homoclinic orbit at (—Un-C, 6 — cHC). Therefore, we expect the codimension-2 curves associated with homoclinic bifurcations discussed below to exhibit the same symmetry about it = 0 as the codimension-2 curves associated with Hopf bifurcations. 5.5.2 Homoclinic Orbits at the Local Extrema of the G = 0 Nulicline In this section, we distinguish between regions on the bifurcation map where the bifurca tion diagrams of the fast subsystem exhibit right, left, or large periodic orbits. We defer the investigation into the existence of corresponding homoclinic orbits to Section 5.5.5. Right periodic orbits undergo the transition to large right periodic orbits and vice versa when the corresponding homoclinic orbits contact the local minimum of the C = 0 nullcline (cf. Figures 5.2 and 5.4). A right homoclinic orbit occurs at the local minimum of the G = 0 nullclirie when DR(Um) = 0. Recalling that Urn = —1 and using (5.54), we thus require eR2(—1)(u — 2) + eRl(—1)u + eRo(—1) = 0. (5.71) Using (5.55)-(5.58), we find eRO(—1) = 12y’r, (5.72) eRl(—1) = —12/r, (5.73) eR2(—1) = 4V’r, (5.74) Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 132 so that (5.71) reduces to 2 / 3N2 3Ti — — = . (5.75) This curve defines a hyperbola in (Ti i%) space, representing a codimension-2 saddle-node loop bifurcation (SNL) [36, 38, 85] at the local minimum of the G = 0 nullcline. Similarly, left periodic orbits become large left periodic orbits and vice versa when the corresponding homoclinic orbits connect to the local maximum of the G = 0 nullcline. A left homoclinjc orbit occurs at the local maximum of the G = 0 nullcline when DL(uM) 0. Recalling that UM = 1 and using (5.60), we thus require eL2(1)(u — Ti2) + eLl(1)u + eLo(1) = 0. (5.76) Using (5.61)-(5.64), we find eLo(1) 12Vr, (5.77) eLl(1) = 12v’r, (5.78) eL2(1) = 4v”ir, (5.79) so that (5.76) reduces to — ( + 3)2 = . (5.80) This curve is the reflection of the hyperbola defined by (5.75) across = 0 and represents the saddle-node ioop bifurcation at the local maximum of the G = 0 nulicline. It is well-known that the transition from a right or left homoclinic orbit to a large homoclinic orbit is not immediate [10, 30]. Corresponding to each of the saddle-node loop bifurcations just discussed, there is a second saddle-node loop bifurcation at which the large homoclinic orbits contact the local extrema of the G = 0 nullcline, so that there is a pair of saddle-node loop bifurcations near each of the local extrema. For parameter values between those at which two corresponding saddle-node loop bifurcations occur, Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 133 Figure 5.10: Illustration of a series of SNIC bifurcations at the local minimum of the G = 0 nullcline, c = 1, resulting in the deformation of a right homoclinic orbit into a large homoclinic orbit. At c = 1, the fast subsystem has two fixed points. The fixed point at U = Urn = —1 is the saddle-node and the fixed point at u = 2 is a spiral source. The direction of the large arrow indicates increasing values of , which are are 0.88, 1.0, 1.1, and 1.2, respectively. The parameters i% and a were held constant at 1.6 and 0.25, respectively. the homoclinic orbit deforms continuously, as shown in Figure 5.10. The deformation is the result of a series of codimension-1 saddle-node on an invariant circle bifurcations (SNIC) [30]. How fast the deformation occurs depends on the value of the parameter a. The smaller a, the faster the deformation occurs and the smaller the distance in parameter space between the two saddle-node loop bifurcations. In the degenerate case, a = 0, the two saddle-node loop bifurcations merge. We now attempt to derive the equation of the curves in parameter space which rep resent the saddle-node loop bifurcations corresponding to the connection of the large homoclinic orbits to the local extrema of the C = 0 nulicline. Large homoclinic orbits Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 134 — - ight and left - - periodic _/‘ - large -riht - orbits - - - .peiodic orbits, - - - / - - right and left - . - - - - - periodic orbits - - Figure 5.11: As Figure 5.8, with the SNL (right) and SNL (left) curves, representing homoclinic bifurcations at the local extrema of the G = 0 nullcline, superimposed. occur at the local minimum of the G = 0 nullcline when D(m) = Ds(1) = 0. We recall that Ds = DL + DR. From (5.61)-(5.64), we find eLo(—1) = eLl(—1) = eL2(—1) = 0, (5.81) so that DL(—1) 0 for all values of the parameters from (5.60). The condition Ds(—1) o thus reduces to DR(—1) = 0, giving the same result as in (5.75). Similarly, large homoclinic orbits occur at the local maximum of the G = 0 nullcline when Ds(uM) = Ds(1) 0, which reduces to DL(1) = 0, giving the same results as in (5.80). We therefore conclude that the Melnikov theory predicts the location of the saddle-node loop bifurcations in the degenerate case a = 0. The two hyperbolae given in (5.75) and (5.80) are referred to as the SNL (right) and SNL (left) curves, respectively, and are plotted in Figure 5.11. To the left of the SNL (right) and SNL (left) curves, bifurcation diagrams of the fast subsystem exhibit 3 2 1 0 0 Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 135 right and left periodic orbits, respectively. To the right of these curves, the bifurcation diagrams exhibit large right and large left periodic orbits, respectively, as indicated in Figure 5.11. Near the two hyperbolae, we keep in mind that there is a small region where homoclinic orbits contact the saddle-nodes at the local extrema of the G = 0 nullcline. We can think of this region as a boundary layer, where the thickness of the layer depends on the value of a, as well as and i. The boundary layer is non-existent in the case a=0. 5.5.3 Double Homoclinic Orbits In the region on the bifurcation map where bifurcation diagrams exhibit one right and one left Hopf bifurcation, there is the possibility of having a codimension-2 double homoclinic bifurcation, also known as a double saddle loop bifurcation (DSL) [36]. At this bifurcation, right and left periodic orbits are homoclinic to the same saddle point on the middle branch of the G = 0 nullcline, resulting in a figure eight. The necessary conditions for this situation are DR(a3,1) = 0, (5.82) DL(a3,2) = 0, (5.83) a8,1 = a8,2. (5.84) From (5.54) and (5.82), we obtain 2 2 eRlu+eRou — =— . (5.85) eR2 Substituting (5.85) in (5.60) and (5.83), we obtain an equation for i in terms of a3. Solving for i and subsequently for from (5.85), we obtain the following parametric equations (a3 being the parameter), eRoeL2 — eLOeR2 u = , (5.86) eLleR2 — eRleL2 Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 136 3 .--i -- __-1 2 — — — — 1 - - - - - - - right homOlinic - right homoclinic ‘S or1iJ-hès below left odit -lies above left ‘::- - - - ilomoclinic orbji Figure 5.12: As Figure 5.11, with the DSL curve, representing simultaneous right and left homoclinic bifurcations, superimposed. The curve is parameterized by a8, for Urn a3 UM, as indicated. L. eRlq = /u2 + —u + —, (5.87) v eR2 eR2 where we remember that eRo, etc., are functions of a3. The curve defined by (5.86)-(5.87) and Urn < a3 < UM is referred to as the DSL curve, and is shown in Figure 5.12. We note that the DSL curve precisely divides the region on the bifurcation map in question, namely the small triangular region bounded by the a-axis, the major TB (left) curve (cf. Figure 5.5), and the SNL (right) curve (cf. Figure 5.11). The DSL curve is symmetric about = 0, crossing the q-axis at a3 0, as expected from the symmetry of the G = 0 nullcline about (u, c) = (0, 3). Staying in the triangular region, bifurcation diagrams of the fast subsystem exhibit the left homoclinic orbit below the right homoclinic orbit when parameter values are chosen to the left of the DSL curve. The relative position of the two homoclinic orbits reverses when parameter Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 137 Figure 5.13: Bifurcation diagrams illustrating that the periodic orbits near a homoclinic orbit can be either stable, as in (a) and (c), or unstable, as in (b) and (d). (a) = 0.93, = 0.97, a = 0.25. (b) i = 0.93, It 0.57, a = 0.25. (c) = 2.0, ft = 2.0, a = 0.02. (d) = 2.03, 1% 1.1, a = 0.02. values are chosen to the right of the DSL curve. 5.5.4 Neutral Homoclinic Orbits In this section, we distinguish between bifurcation diagrams in which the periodic orbits near a homoclinic bifurcation are stable and those in which the periodic orbits are unsta ble, as shown in Figure 5.13. The stability of periodic orbits near homoclinic bifurcations is determined by the trace of the Jacobian matrix at the saddle point. If the trace is positive, the periodic orbits are stable, whereas if the trace is negative, the periodic or bits are unstable [36, 38]. We compute curves in the fast parameter space along which (a) 6 i I I I (b) 6 5 .........,. - 5 4 \ .\ . 4 ci 3 •1 ci 3 2 \...-- 2 1 \. 10 I I I 0 -3 -2 -1 0 1 2 3 U (c) 6 I I (d) 6 5- - 5 4- - 4 ci \ - ci 3 2- 2 1- / \./ 10 I 0 -4 -2 0 2 4 j- //\• :(/\ 3 U -4 -2 0 2 4 U Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 138 the trace of the Jacobian matrix changes sign with a simple zero, representing neutral homoclinic orbits, also known as neutral saddle loops (NSL). In light of (5.14), the required conditions for a neutral right homoclinic orbit are F(a3) = 0, (5.88) DR(a3) 0. (5.89) Using (5.3), (5.88) can be written as — = 2a3’Ci — a. (5.90) Substituting (5.90) into (5.54) and (5.89), we obtain an equation for i in terms of a3. Solving for and subsequently for from (5.90), we find the following parametric equa tions (again with a3 being the parameter), eR2a — eRo = , (5.91) 2a3eR2 + = Ia3 — (5.92) where we remember that eRo, eRl, and eR2 are functions of a3. The conditions for a neutral left homoclinic orbit are as (5.88)-(5.89), but rather than DR(a8) = 0, we require DL(a3) = 0. The resulting parametric equations are exactly the same as (5.91) and (5.92), but is replaced by e, eRl by eLi, and cR2 by eL2. Similarly, the parametric equations for the case of a neutral large homoclinic orbit are also the same as (5.91) and (5.92), but eFo is replaced by eso, cR1 by esi, and eR2 by es2. The curve defined by (5.91)-(5.92) and Urn <a3 <UM is referred to as the NSL (right) curve, and its analogs for the left and large homoclinic orbits are referred to as the NSL (left) and NSL (large) curves, respectively. The curves are added to the bifurcation map in Figure 5.14. We make the following observations about the new curves. The NSL (right) lies entirely in the first quadrant of the bifurcation map. The NSL (left) curve Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 139 2.5 Figure 5.14: As Figure 5.12, with the NSL (right), NSL (left), and NSL (large) curves, representing neutral homoclinic bifurcations, superimposed. The curves are parameter ized by a3, as indicated. is the reflection of the NSL (right) curve across = 0 and lies entirely in the fourth quadrant. The NSL (large) curve meets the NSL (right) curve on the SNL (right) curve (cf. Figure 5.11) at (, i) = (, ) when a3 = Urn. This point has special significance, as will be discussed in Section 5.8 below. As a3 increases, the curve remains in the first quadrant of the map, until there is a discontinuity at a3 = 0, at which point it reappears in the fourth quadrant. As a3 increases further, the curve traces the reflection of the corresponding portion in the first quadrant and meets the NSL (left) curve on the SNL (left) curve at (, It) = (, —) when a3 = UM (not shown). For parameter values between the NSL curves, the bifurcation diagrams of the fast subsystem exhibit unstable homoclinic orbits, whereas the homoclinic orbits are stable for parameter values lying outside the NSL curves. 3 2 1 0 0 0.5 1 1.5 2 1 Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 140 5.5.5 Coalescing Homoclinic Orbits From Figures 5.2 and 5.3, we observed that the bifurcation diagram of the fast subsystem may or may not exhibit right homoclinic orbits when there are right periodic orbits. Similarly, the existence of left and large homoclinic orbits is not necessarily implied by the existence of corresponding periodic orbits. In this section, we distinguish between bifurcation diagrams exhibiting homoclinic orbits and those that do not. Figure 5.15 shows bifurcation diagrams of the fast subsystem for two sets of parame ters. In both cases, there are two right Hopf bifurcations and corresponding right periodic orbits. However, in Figure 5.15a there are two right homoclinic orbits, whereas in Fig ure 5.15b the homoclinic orbits have disappeared. Based on these observations, we look for a curve on the bifurcation map along which two right homoclinic orbits coalesce. The right Melnikov distance functions, DR(a3), for the two cases just described are plotted in Figure 5.15c. We recall that DR(a3) = 0 gives the leading-order value of a3 at which a right homoclinic orbit occurs. Thus, it is clear that right coalescing homoclinic orbits, or right coalescing saddle loops (CSL), must occur when the conditions DR(a3) = 0, (5.93) dDR(a3) = 0, (5.94) hold. Technically, we also required2DR(a3) 0, but we do not insist on this condition in the derivation below, since it is generically satisfied. We solve (5.93) and (5.94) simultaneously. From (5.54), we obtain dDR(a8) da3 =42(a3)(i% — ,2) + e’Rl(as)u +40(a3), (5.95) where0(a3), e’Rl(a$), and e2(a3) are the derivatives of eRo(a3), eRl(a3), and eR2(aS), Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 6 6 5 5 4 4 c) 3 c) 3 2 2 1 1 0 0 Figure 5.15: (a) Bifurcation diagram exhibiting two right Hopf bifurcations and two corresponding right homoclinic orbits ( = 1.25, ft = 2.54, and a 0.25). (b) Bifurcation diagram also exhibiting two right Hopf bifurcations but no homoclinic orbits ( = 1.25, ft = 2.6, and a = 0.25). (c) Right Melnikov distance function for the two above cases. The solid curve has two zeroes for 1 < c < 5 and corresponds to (a). The dashed curve has no zeroes for 1 <c < 5 and corresponds to (b). respectively, easily obtained from (5.55)-(5.58). Using (5.54) and (5.93), we obtain — — eR1u + eRo (5.96) eR2 Substituting (5.96) into (5.94) and (5.95), we again obtain an equation for ii. Solving for ft and subsequently for from (5.96), we obtain the following parametric equations for ft and i, = eoeR2 — eRoe2 (a) 141 (b)I I I I I \ \: -3 -2 -1 0 1 2 3 4 (c) -3 -2 -1 0 1 U c) 23 4 0.6 0.4 0.2 0 -0.2 eR1e?2 —e1R2 (5.97) Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 142 3 2 as = UM 1 a5 — UJJ 0 Figure 5.16: As Figure 5.14, with the curves of coalescing homoclinic orbits, CSL (right), CSL (left), and CSL (large), superimposed. The curves are parameterized by a3, as indicated. = + + . (5.98) cR2 eR2 The conditions for left and large coalescing homoclinic orbits are as in (5.93)-(5.94), with DR replaced by DL and D, respectively. The resulting parametric equations for these two cases are as in (5.97)-(5.98), with eRo, eRl, and eR2 replaced accordingly. The curve defined by (5.97)-(5.98) and tim <a8 <UM is referred to as the CSL (right) curve, and its analogs for the left and large homoclinic orbits are referred to as the CSL (left) and CSL (large) curves, respectively. The curves are added to the bifurcation map in Figure 5.16. The portion of the CSL (right) curve in the first quadrant of the bifurcation map lies entirely in the region where bifurcation diagrams exhibit two right Hopf bifurcations and the corresponding periodic orbits are small instead of large. That is, this portion of a5 urn ‘7 Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 143 the curve represents the situation inferred in Figure 5.15, as expected. The remaining portion of the curve lies in the fourth quadrant, opposite the CSL (left,large) curve across = 0. We discuss the meaning of this portion a little later. The CSL (large) curve is symmetric about a3 0 on ii = 0, and follows the SNL (right) and SNL (left) curves very closely (cf. Figure 5.11), except for values of a3 near zero. For a3 near zero, the curve deviates significantly to the left of the SNL curves. We now have an apparent paradox, in that we deduced earlier that large periodic orbits exist only to the right of the SNL curves, up to some thin boundary layer. Yet the curve of coalescing large homoclinic orbits lies well to the left of these curves. We use DSTOOL and AUTO to resolve the paradox numerically. We focus our attention on the same triangular region as discussed in connection with the DSL curve (cf. Section 5.5.3), namely the region bounded by the axis, the major TB (left) curve, and the SNL (right) curve. In this region, each bifurcation diagram is guaranteed to exhibit both right and left periodic orbits, with corresponding right and left homoclinic orbits. In order for there to be a curve of coalescing large homoclinic orbits in this region, large periodic orbits also must exist. DSTOOL indeed confirms this for , = 1.64, % = 0, and a = 0.25, as shown in Figure 5.17. A regular one-parameter continuation with AUTO, using c as the bifurcation parame ter and keeping and ft fixed at 1.64 and 0, respectively, fails to reveal these large periodic orbits. We hypothesize that any branch of periodic orbits must be detached from the main bifurcation diagram, which AUTO, therefore, cannot find directly. We make use of a two-parameter continuation, discussed in more detail in Section 5.8 below, to find an initial condition on the detached branch from which start the required one-parameter continuation. Bifurcation diagrams so obtained are shown in Figure 5.18. Figure 5.18a shows an open, detached branch of large periodic orbits terminating at large homoclinic bifurcations Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 144 6 4 2 . 0 -2 -4 -6 U Figure 5.17: Phase portrait of the fast subsystem, (5.6), at c = 3, showing that a large periodic orbit can exist in addition to right and left periodic orbits for certain values of the model parameters. The right and left periodic orbits are both unstable and the large periodic orbit is stable. The remaining parameter values are i 1.64, i% = 0, and a = 0.25. at both endpoints, for parameter values just to the left of the SNL (right) curve. When is decreased, corresponding to moving leftward on the bifurcation map, the detached branch shrinks in size, and the large homoclinic orbits approach each other, as shown in Figure 5.18b. Figure 5.18c shows that when q is decreased further, the large homoclinic orbits coalesce as predicted by the CSL (large) curve, resulting in a closed, but still detached branch of periodic orbits. For even smaller values of , the closed branch further shrinks in size and disappears when the two saddle-node of periodics bifurcations at the top and bottom of this branch coalesce, leaving only the main branch of the bifurcation diagram. The coalescence of the SNP bifurcations is significant and discussed in more detail in Section 5.6 below. We note that the relative position of the right and left homoclinic orbits is reversed in Figures 5.18a and 5.18b. This reversal is the result of -4 -2 0 2 4 Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 145 (a) 6 i i i i i (b) 6 i i i U U (c) 6 i i i i 4 : _____ Figure 5.18: Bifurcation diagrams illustrating the existence of a branch of large periodic orbits which is detached from the main branch for certain values of the parameters. (a) = 1.64. (b) = 1.54919. (c) = 1.5255. The other parameters are = 0 and a 0.02. Areas enclosed by the boxes are enlarged in the figures to the right of the bifurcation diagrams in (b) and (c). Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 146 crossing the DSL curve (cf. Section 5.5.3) as is decreased. We now briefly investigate the dependence of the CSL curves on the parameter a3, as indicated in Figure 5.16. The CSL (right) curve starts at (ii, i) (0, 1) when a3 = As a3 decreases, both and i approach infinity, and the curve approaches the upper half of the SNL (right) curve. However, as a3 approaches Urn, there is a discontinuity and the curve unexpectedly reappears in the fourth quadrant, virtually indistinguishable from the lower half of the SNL (right) curve there. When a3 = Urn, the curve ends at the intersection of the SNL (right) and SNL (left) curves, at (, 2) = 0). Similarly, the reflection of the CSL (right) curve, namely the CSL (left) curve, lies mostly in the fourth quadrant, but appears in the first quadrant for values of a3 very close to UM and is virtually indistinguishable from the SNL (left) curve here. The CSL (large) curve is virtually indistinguishable from the SNL curves for all values of a3, except values of a3 close to zero, as noted above. For values of a3 between 0 and Urn, the curve lies mainly in the first quadrant near the upper half of the SNL (right) curve, whereas for values of a3 between 0 and UM, the curve lies mainly in the fourth quadrant near the lower half of the SNL (left) curve. When a3 approaches Urn or UM, there is a discontinuity at both ends of the curve. They reappear in the opposite quadrants, essentially at the same location as the ‘tails’ of the CSL (right) and CSL (left) curves. We interpret the bifurcations represented by these unexpected portions of the CSL curves with the aid of Figure 5.19. Figure 5.19a shows the bifurcation diagram of the fast subsystem for parameter values just above the CSL (left,large) curve which, for our pur poses, is equivalent to the SNL (left) curve in the bottom right corner of the bifurcation map. We note that are left and right Hopf bifurcations, with corresponding left and large right homoclinic orbits, respectively. According to the analysis in Section 5.5.2, the left homoclinic orbit undergoes the transition to a large left homoclinic orbit when moving inside the SNL (left) curve. Considering this codimension-2 bifurcation alone, we thus Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 147 (a) 8 i , i i i —10 I I I I I I I -5-4-3-2-10 1 2 34 5 6 U (b) 10 i i i i i i i i i i —15 I I I I I I -5-4-3-2-1 0 1 2 3 4 5 6 U Figure 5.19: Bifurcation diagrams illustrating that large left and large right homoclinic orbits do not coexist. (a) Bifurcation diagram for = 2.2, i% = 0.46. There are right and left Hopf bifurcations, with corresponding large right and (small) left homoclinic orbits. (b) Bifurcation diagram for i = 2.0, €i = 0.4. There no longer are any homoclinic orbits. For both diagrams, a = 0.25. should expect bifurcation diagrams to exhibit both large left and large right homoclinic orbits. However, this is not the case, as shown in Figure 5.19b. The bifurcation diagram in this figure still exhibits left and right Hopf bifurcations, but there no longer are any ho moclinic orbits. Moving across these curves has caused the left and large right homoclinic orbits to contact the local maximum of the G = 0 nuficline at u = a8 = UM essentially simultaneously, corresponding to coalescing left homoclinic orbits and coalescing large homoclinic orbits, as predicted. Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 148 Based on the interpretation of all portions of the CSL curves, we now can infer the number of homoclinic orbits in each region of the bifurcation map, as indicated in Figure 5.16. 5.6 Constraints on the Saddle-Node of Periodics Bifurcations In this section, we examine the codimension-2 bifurcations associated with the saddle- node of periodics bifurcation. We can deduce the emergence of some SNP bifurcations by examining the criticality of a Hopf bifurcation and the stability of the corresponding homoclinic bifurcation (if it exists), and use AUTO to continue the study of the coalescing SNP bifurcations mentioned in the previous section. From Figures 5.8 and 5.14, we observe that the DHB (right) curve, representing de generate Hopf bifurcations, lies above the NSL curves, representing neutral homoclinic orbits. Furthermore, the regions in which Hopf bifurcations are subcritical and homo clinic orbits are stable overlap. The existence of subcritical Hopf bifurcations and stable homoclinic orbits at the same time necessarily implies the existence of at least one saddle- node of periodics (cf. Figure 5.4). That is, the DHB (right) and NSL curves also represent the emergence of an SNP bifurcation. The curves are redrawn in Figure 5.20 and the number of saddle-node of periodics is indicated in the regions bounded by these curves. Below the NSL curves, there are either no or two saddle-node of periodics bifurca tions. To divide this region, we return to the observation in the previous section that the detached branch of periodic orbits, exhibited by bifurcation diagrams in the region just below where the NSL curves meet, disappears through the coalescence of two SNP bifurcations (cf. Figure 5.18). Based on this observation, we infer the existence of the codimension-2 CSNP curve, for coalescing SNP bifurcations, which precisely divides this region in two. Points on this curve were computed with AUTO for several values of . Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 149 3 I !::‘ - - - - - - - - DHB no adcl1e-nodes - - - -‘ (right) -ofperiodics 2 - SL --4- - -. (large) - - - onèsàddle-node - - - -1 . ofpe?kdics ---- - XTSL tw-add1e-nodes - • - of periodics - -(right) -- - no sadd1e-nde - - - - - - - of periodics -.. - - csIp: - - -‘ - 0 I -..L- -- .1- 0 0.5 1 1.5 2 2.5 •17 Figure 5.20: As Figure 5.16, with the CSNP curve, representing coalescing SNP bifurca tions on the detached branch of periodic solutions (cf. Figure 5.18), superimposed. Also highlighted are the DHB (right), NSL (right), and NSL (large) curves, which represent degenerate Hopf bifurcations on the right branch of the G = 0 nullcline, and neutral right and large homoclinic orbits, respectively, as well as the emergence of a saddle-node of periodics. We note that this curve is the only curve for which we do not have an analytic equation. Finally, we need to slightly adjust the qualitative features of the periodic orbits in dicated in Figure 5.11. In the region between the CSNP and SNL (right) curves, large periodic orbits exist in addition to the right and left periodic orbits. We note that the large periodic orbits are not associated with a right or left Hopf bifurcation. 5.7 The Complete Bifurcation Map At this point, all constraints on the codimension-1 bifurcations have been derived. The complete bifurcation map is shown in Figures 5.21 and 5.22. Based on the superposition of the qualitative features of the Hopf, homoclinic, and SNP bifurcations indicated on Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 150 Figure 5.21: The complete bifurcation map. Corresponding bifurcation diagrams are shown in the figures in Appendix B. The qualitative features of the Hopf, homoclinic, and SNP bifurcations exhibited by the bifurcation diagrams in each labelled region are summarized in Table 5.1. The regions not labelled in the box surrounding the CSNP curve are shown in Figure 5.22. The horizontal lines at % = 0.8 and 1.5 are discussed in Section 5.8. LD 00 fl Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 151 0.5 0.25 0 Figure 5.22: Enlarged view of a section of the bifurcation map (cf. Figure 5.20), iden tifying regions F1-F4. Corresponding bifurcation diagrams are shown in the figures in Appendix B. The qualitative features of the Hopf, homoclinic, and SNP biftircations and the nature of the detached branch of periodic orbits exhibited by the bifurcation diagrams in these regions are summarized in Table 5.2 and its caption. the intermediate versions of the map, we can deduce the qualitative structure of the bifurcation diagram in each region. We thus have answered the question posed in the introduction to this chapter, namely, can we deduce the bifurcation structure of the fast subsystem given the values of the two parameters 2 and ? The qualitative features of the Hopf, homoclinic, and SNP bifurcations exhibited by the bifurcation diagrams in the regions identified in Figure 5.21 are summarized in Table 5.1. In regions F1-F4, identified in Figure 5.22, all bifurcation diagrams exhibit right and left subcritical Hopf bifurcations in the region of bistability, with corresponding unstable right and left homoclinic orbits, respectively. The relative position of the right and left homoclinic orbits is summarized in Table 5.2, as well as the nature of the detached branch of large periodic orbits. We mention again that an example of each type of bifurcation diagram is included in Appendix B. The regions are roughly grouped according to the number of Hopf bifurcations exhib ited by the bifurcation diagrams and whether corresponding homoclinic orbits are small 1 1.25 1.5 1.75 2 Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 152 Type Sub- or Does HB lie Nature of Number Diagram Region of super- in region of corresponding of shown RB critical bistability? HC SNP’s in Figure A - - - - - - B1 major right sub yes small unstable 0 5.13b B2 major right sub yes small stable 1 5.13a B3 major right sub no small stable 1 B4 major right super no small stable 0 5.2 B5 major right super yes small stable 0 5.6 Cl major right super no large stable 0 minor right super yes small stable C2 major right super no small stable 0 5.15a minor right super yes small stable C3 major right super yes small stable 0 minor right super yes small stable C4 major right super yes — 0 minor right super yes — C5 major right super no — 0 5.3, 5.15b minor right super yes — C6 major right super no — 0 minor right super no — Dl major right sub no — 2 major left sub no — D2 major right sub no — 2 major left sub yes — D3 major right sub yes — 2 major left sub yes — D4 major right sub yes large unstable 2 major left sub yes small unstable D5 major right sub no large unstable 2 major left sub yes small unstable El major right sub no large unstable 2 5.13d E2 major right sub yes large unstable 2 E3 major right sub yes large stable 1 E4 major right sub no large stable 1 5.4, 5.13c E5 major right super no large stable 0 Table 5.1: Summary of the qualitative features of the Hopf, homoclinic, and SNP bifur cations exhibited by bifurcation diagrams of the fast subsystem in the labelled regions on the complete bifurcation map in Figure 5.21. Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 153 Region Relative location Nature of the Number Diagram of the left and right detached branch of of shown homoclinic orbits periodic orbits SNP’s in Figure Fl left below right — 0 closed branch, F2 left below right no large 2 5.18c homoclinic orbits open branch, with F3 left below right two large unstable 2 homoclinic_orbits open branch, with F4 right below left two large unstable 2 5.18a homoclinic_orbits Table 5.2: Summary of the qualitative features of the homoclinic and SNP bifurcations and the nature of the detached branches of periodic orbits exhibited by bifurcation di agrams in regions F1-F4 on the bifurcation map. In each region, there is a major right and left subcritical Hopf bifurcation (in the region of bistability) with corresponding right and left homoclinic bifurcations. or large. For example, bifurcation diagrams in regions B1-B5 all exhibit a major right Hopf bifurcation, and corresponding branches of periodic orbits terminate at a small right homoclinic orbit. Across the SNL (right) curve, bifurcation diagrams in regions E1-E5 also exhibit a major right Hopf bifurcation, but branches of periodic orbits ter minate at large right homoclinic orbits. Within each group, moving from one region to a neighbouring region across a codimension-2 curve generally corresponds to changing just one qualitative feature of the bifurcation diagram. For example, moving from C5 to C6 across the minor HBL (right) curve corresponds to changing the location of the minor right Hopf bifurcation from being inside the region of stability to being outside, and so on. In light of the results in Section 5.6, moving across the DHB or NSL curves corresponds to changing two qualitative features of the bifurcation diagram, namely the Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 154 criticality of a Hopf bifurcation or the stability of a homoclinic orbit, respectively, and the number of SNP bifurcations. 5.8 Two-Parameter Bifurcation Diagrams Up until now, we have concerned ourselves with the qualitative features of bifurcation diagrams only. However, quantitative information pertaining to the exact locations of the four types of codimension-1 bifurcations exhibited by the fast subsystem, namely the saddle-node, Hopf, homoclinic, and saddle-node of periodics bifurcations, cannot be obtained from the bifurcation map. In this section, we generate two-parameter bifurcation diagrams to assist us in obtain ing more quantitative information, without resorting to the one-parameter bifurcation diagrams. We can think of the two-parameter bifurcation diagrams as being a step in be tween the very quantitative one-parameter bifurcation diagrams and the very qualitative bifurcation map, as follows. One-parameter bifurcation diagrams indicate the location of the fixed points and their stability, and show the locations of the codimension-1 bifurca tions. Two-parameter bifurcation diagrams show how the locations of the codimension-1 bifurcations depend on a second parameter. Curves on the two-parameter bifurcation diagrams are thus curves of codimension-1 bifurcations. We can infer codimension-2 bi furcations (the locations of which we indicate with a bullet in the figures below) from the two-parameter bifurcation diagram. The bifurcation map in turn consists of curves of codimension-2 bifurcations. One of the parameters in the two-parameter bifurcation diagrams is c, our primary bifurcation parameter. The other parameter can be either i or ü, and we choose the former, as the results then can be most easily compared to the results obtained by Bertram et al. [10] in Section 5.9 below. In Section 5.8.1, we examine two-parameter Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 155 4 3 - 2 1 0 Figure 5.23: Two-parameter bifurcation diagram generated with AUTO for i = 1.5 and a = 0.25. bifurcation diagrams computed with AUTO and relate them to the bifurcation map. However, equations for all curves in these diagrams, except those representing the SNP bifurcations, can be derived analytically based on the results in Sections 5.4 and 5.5. We do this in Section 5.8.2, and also compare the analytical and numerical two-parameter bifurcation diagrams. 5.8.1 Numerical Two-Parameter Bifurcation Diagrams Figure 5.23 shows the two-parameter bifurcation diagram computed with AUTO for = 1.5. The curves of left SN’s and right SN’s represent the location of the left and right saddle-nodes at the local extrema of the G = 0 nullcline, respectively. Since the location 3 C Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 156 of the saddle-nodes does not depend on the model parameters, these curves are vertical lines at C = Cm = 1 and c = CM = 5. The curves of major right RB’s, minor right RB’s, and left HB’s represent the locations of the respective Hopf bifurcations. For example, for 0.25, there are two right Hopf bifurcations, namely a major one at c 2.9 and a minor one at c 4.8. Similarly, the curves of right IIC’s, left HC’s, and large HC’s represent the locations of the respective homoclinic orbits and the curve of SNP’s represents the location of the saddle-node of periodics. We note that homoclinic bifurcations only occur at values of C in the region of bistability and, hence, the corresponding curves in the two-parameter bifurcation diagram remain between c = 1 and c = 5. In contrast, the curves of RB’s and SNP’s continue outside this region, reflecting the fact that Hopf and SNP bifurcation can occur at any value of c. We relate the two-parameter bifurcation diagram in Figure 5.23 to a horizontal cut through the bifurcation map at z2 = 1.5 (cf. Figure 5.21) as follows. Each intersection of the horizontal line with a curve on the bifurcation map corresponds to a codimension-2 bifurcation, reflected by a bullet in the two-parameter bifurcation diagram. For example, the intersection with the CSL (right) curve corresponds to the emergence or disappear ance of two right homoclinic orbits. In the two-parameter bifurcation diagram, this event is reflected by the bullet at the local minimum of the curve of right HC’s. As mentioned in Section 5.5.2, the transitions from small to large homoclinic orbits at the local extrema of the G = 0 nullcline or saddle-nodes are not immediate, but occur over a range of parameter values in conjunction with SNIC bifurcations. This is clearly evident in Figure 5.23, where the curves of right and left HC’s and the curve of large HC’s join the curves of right and left SN’s at different values of j. The analysis in Section 5.5.2 failed to yield curves on the bifurcation map corresponding to the codimension-2 SNL bifurcations at which large homoclinic orbits contact saddle-nodes. This failure is reflected in the two parameter bifurcation diagram by the omission of bullets at the points where the curve Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 157 of large HC’s joins the curves of right and left SN’s. Finally, we observe that the curve of SNP’s emanates from the bullet corresponding to the codimension-2 NSL (large) bifurcation, as expected from Figure 5.20 and the discussion in Section 5.6. The bullet corresponding to the DHB (right) bifurcation is not shown but occurs at c —7.21. Another curve of SNP’s emerges from this point, also as expected. We now examine the dependence of the two-parameter bifurcation diagram on Returning to the bifurcation map, we observe two points at which curves meet in a cusp- like fashion, namely (, i) = (0, 1) and (, ü) = ). We discuss the changes in the two-parameter bifurcation diagram as they relate to these points. As t decreases from ‘11 = 1.5, the minor TB (right) bullet moves down along the curve of right SN’s towards q = 0, and is exchanged for a major TB (right) bullet when = 1, as expected from the TB (right) curves on the bifurcation map. Similarly, the DHB (right) bullet (not shown in Figure 5.23) moves down along the curve of right HB’s, and is exchanged for an NSL (right) bullet when i = 1. For < 1, the curve of SNP’s originally emanating from the DHB (right) bullet emanates from the NSL (right) bullet instead. In addition, both CSL (right) and CHB (right) bullets move towards the minor TB (right) bullet as i decreases and vanish when i = 1. Figure 5.24a shows the resulting two-parameter bifurcation diagram for ?% = 0.8. As i’t decreases further, the NSL (right) and NSL (large) bullets approach each other, until they meet on the curve of left SN’s when (, ?%) = (, ). At this point, the two curves of SNP’s join, and immediately detach from the curves of HC’s for values of ü < This detachment is clearly shown in the two-parameter bifurcation diagram for i = 0 in Figure 5.24b. We can identify a local minimum on the curve of SNP’s, corresponding to the coalescence of two saddle nodes of periodics, reflected by the CSNP bullet. Sim ilarly, the curve of large HC’s now also clearly has a local minimum, corresponding to Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 158 (a) (b) Figure 5.24: Two-parameter bifurcation diagrams generated with AUTO for (a) i2 = 0.8 and (b) = 0. In both cases, a = 0.25. 3 2 1 0 0 1 2 3 4 5 6 C 2 1 0 1 2 3 4 5 6 C Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 159 the coalescence of two large homoclinic orbits, and reflected by the CSL (large) bullet. Finally, double homoclinic orbits are indicated by the DSL bullet at the intersection of the curves of right and left HC’s. We note that Figure 5.24b exhibits symmetry about c = 3 in light of the discussion in Section 5.2. Furthermore, small imperfections in the curves of HC’s are the result of AUTO having some difficulties with the computations. 5.8.2 Analytical Two-Parameter Bifurcation Diagrams In this section, we derive the equations of the curves of RB’s and HC’s based on the results of Sections 5.4 and 5.5, and compare the analytical two-parameter bifurcation diagram so obtained with the numerical version. For the curves of RB’s, we use (5.3)-(5.4) to rewrite (5.16) and (5.18) as c = 3(u + 1) — t, (5.99) = i.—uI, (5.100) so that c and are parameterized by u. Values of u UM generate the curves of right RB’s, whereas values of u Urn generate the curve of left HB’s. For the curve of right HC’s, we require DR(a3) = 0. Using (5.54) to isolate and U = a3 in (5.99), we again obtain parametric equations for c and , = 3(a+1)—a, (5.101) = j2+ eR1(aS)+ eRo(a3) (5.102) N eR2(a3) eR2(a8) where now a3 is the parameter, for Urn a3 ‘aM. The curves of left and large HC’s are given by (5.101)-(5.102) with the eR’s replaced by 6L5 and es’s, respectively. The curves defined by (5.99)-(5.100) and (5.1O1)-(5.102) and the trivial curves of right and left SN’s are plotted in Figure 5.25a for it = 1.5. Figure 5.25b shows the comparison of the analytical two-parameter bifurcation diagram with the numerical versions generated Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 1 2 3 Figure 5.25: (a) Analytical two-parameter bifurcation diagram for ii = 1.5. (b) Compar ison with AUTO results (solid curves) for a = 0.25, a = 0.15, a = 0.05, and a 0.01. The arrow indicates the direction of increasing a. 160 major HBL! , (left leftHB :::NL (left) major TB (left) :-, : SNL (right):: CSL major HBTJ NB’s (right) mnor TB (right) - - - - (right) I I 0 (a) 4 3 :- 2 1 0 (b) 4 3 s::- 2 1 0 C CHB ‘ (right) 5 6 C Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 161 with AUTO for several values of the parameter a. We observe that the agreement between the respective curves of right and left HB’s is exact for all values of a. The agreement between the respective curves of right and left HC’s also is very good. However, the agreement between the respective curves of large HC’s is highly dependent on the value of a. The smaller the value of a, the better the agreement. From (5.3), (5.6), and (5.44) we see that the value of a affects the size of the pertur bation to the Hamiltonian system. Recalling that the curves of HC’s are derived from Melnikov’s method, which is a leading-order method based on the perturbation of the Hamiltonian system, the dependence of the agreement between the analytical and nu merical results on a is not surprising. What is interesting, however, is the observation that the value of a hardly affects the curves of right and left HC’s, but significantly affects the curve of large HC’s. The reason why this is so lies in the effect of the size of the perturbation on the deviation of the relevant stable and unstable orbits of the saddle point from the separatrix of the Hamiltonian system (cf. Figure 3.8). For the curve of right (left) HC’s, the relevant orbits are the right (left) stable and unstable orbits which are hardly affected by the size of the perturbation (cf. Figure 3.3). In contrast, for the curve of large HC’s, the relevant orbits are the large ones and these are significantly affected by the size of the perturbation term (cf. Figure 5.9). Finally, we note that decreasing the value of z% also decreases the size of the pertur bation term. That is, the analytical and numerical two-parameter bifurcation diagrams for it = 0.8 and i% = 0 correspond very closely even for relatively large values of a. 5.9 On the Classification of Bursting Oscillations The classification of bursting oscillations was begun by Rinzel [77], who identified three types of bursting oscillations. Rinzel’s scheme is based on an examination of the spike Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 162 frequency profile during the active phase and a characterization of the initiation and termination of the active phase in terms of a one-parameter bifurcation diagram. In [10], Bertram et al. elaborate on and extend Rinzel’s classification scheme. Their scheme is put into the context of a sequence of horizontal traversals through a two- parameter bifurcation diagram. For example, in Figure 5.23, a horizontal bar at = 0.75 starting at the curve of left SN’s and terminating at the curve of right TIC’s corresponds to a bursting solution in which the active phase starts at a saddle-node and terminates due to the passage through a right homoclinic bifurcation. The two-parameter diagram shown in [10] is generated with the biophysically based Chay-Cook model. Surprisingly, this two-parameter bifurcation diagram is qualitatively identical to the one shown in Figure 5.23, except that the secondary parameters on the vertical axes, i and i, respectively, are inversely related. In the Chay-Cook model, 1tt affects the time constant of the relaxation equation for the activation variable n (cf. (2.11)). Drawing an analogy between the two parameters, we may interpret in the same way. From the similarity of the two-parameter bifurcation diagrams generated with the Chay-Cook and polynomial models, we conclude that the polynomial analog model not only superficially mimics various types of solution behaviour of models of excitable cells, but also possesses the rich bifurcation structure of such models. Therefore, we suggest that the classification scheme of Bertram et al. can be taken a step further by considering bifurcation maps. In particular, such maps provide information about the relationships among the various types of oscillatory behaviour as evident from the bifurcation diagrams in the different regions. In terms of the bifurcation map for the polynomial model, the types of bursting behaviour classified by Bertram et al. only form a subset of the set of possible types of oscillatory behaviour. In particular, none of the behaviour implied by bifurcation Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 163 diagrams in the lower portion of the bifurcation map have been classified. The question of precisely how the classification of Bertram et al. can be extended further based on the different qualitative features of the bifurcation diagrams in each of the regions on a bifurcation map is not pursued in this thesis, but is left as an opportunity for future work. 5.10 Effect of e on the Solution Behaviour of the Full System of Equations In this final section, we briefly examine the effect of e on the solution behaviour of the full system of equations, (5.1)-(5.2). The reason for initiating this study is based on an observation made during the analysis of one of the biophysical models of bursting electrical activity in pancreatic /3-cells studied in the previous chapters, namely the Chay (1986) model. When the values of the parameters given in [13] are used, square-wave bursting is obtained, as for the other biophysical models under consideration in this thesis. However, Melnikov’s method fails to determine the location of any homoclinic bifurcation, since the Melnikov distance for this model is positive for all values of the bifurcation parameter z between Zm and zM. Upon further investigation with AUTO, we find that the fast subsystem of the Chay (1986) model has the bifurcation structure of a tapered burster (cf. Figure 5.3), rather than the expected bifurcation structure of a square-wave burster (cf. Figure 2.6). We point out that the positive Melnikov distance for all values of z between Zm and zM is in accordance with the bifurcation structure of a tapered burster. Thus, the failure to determine the location of the homoclinic bifurcation is not the result of incorrectly applying Melnikov’s method, but rather, it is the result of inappropriately applying the method to a fast subsystem which exhibits no homoclinic orbits. We note that for the purposes of Chapters 2-4, we have adjusted the values of the parameters slightly so that Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 164 the bifurcation diagram of the fast subsystem is as shown in Figure 2.6b, as indicated in Appendix A.2. We now have an apparent paradox, in that the fast subsystem of the Chay (1986) model is that of a tapered burster, yet the numerical solution of the full system of equations exhibits square-wave bursting. The same phenomenon is easily reproduced with the polynomial analog model, as shown in Figure 5.26. However, when the value of the parameter E is decreased, the numerical solution indeed exhibits tapered bursting, as shown in Figure 5.27. From (5.2), we observe that e governs the rate at which the slow variable c changes. For the smaller value of e, the dynamics of c are sufficiently slow that the solution trajectory remains close to the envelope of periodic orbits and, thus, can ‘sneak through’ the small opening between the middle branch of the G = 0 nulicline and the left branch of periodic orbits and reach the upper Hopf point, as expected. We believe that the nongeneric behaviour observed above is typical near boundaries of regions on the bifurcation map. That is, although the bifurcation diagrams of the fast subsystem do change as predicted by the codimension-2 curves on the map, the transitions in the corresponding numerical solutions of the full system of equations are e-dependent. Wiggins [92] includes a word of caution on the interpretation and application of bifurcation diagrams relating to the same type of phenomenon. Dynamical systems having parameters that change in time, no matter how slowly, and that pass through bifurcation values may exhibit behaviour that is different from the situation where the parameters are constant. For references to detailed analyses of such problems, see [92]. Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 165 6 5 4 c) 3 2 1 0 U Figure 5.26: Illustration of an apparent paradox: the bifurcation diagram of the fast subsystem is that of a tapered burster, yet the numerical solution of the full system of equations exhibits square-wave bursting. The bifurcation diagram shown is for = 0.7, = 1.95, and a = 0.25, lying in region C5 of the bifurcation map, just above the CSL (right) curve. The projection of the numerical solution of the full system of equations, (5.1)-(5.2), is superimposed for = 0.0025. The values of the remaining parameters are as in Figure 5.2. -3 -2 -1 0 1 2 3 4 Chapter 5. Multiple Bifurcations in a Polynomial Analog Model 166 6 5 4 c) 3 2 1 -3 U Figure 5.27: Resolution of the apparent paradox. When the value of e is decreased to = 0.00 1, the numerical solution of the full system of equations exhibits tapered bursting, as expected from the bifurcation diagram of the fast subsystem. All other parameters are as in Figure 5.26. 0 -2 -1 0 1 2 3 4 Chapter 6 Conclusions and Discussion In this thesis, we have presented a variety of mathematical methods and techniques to study models of bursting oscillations. In particular, we have focused on models of bursting electrical activity (BEA) in the membrane potential of pancreatic /3-cells and a polynomial analog model of the bursting phenomenon. In the first part of this thesis, consisting of Chapters 2, 3, and 4, we investigated the determination of the plateau fraction for several models of BEA in pancreatic /3-cells. The plateau fraction is physiologically important because it appears to be strongly correlated to the rate of insulin release. We have extended the techniques of Pernarowski [67] and Pernarowski et al. [69, 70], developed in the context of the Sherman-Rinzel-Keizer model, to a class of first-generation models which describe BEA in pancreatic /3-cells and which consist of three ordinary differential equations. These models have in common the assumption that intracellular Ca2+ concentration is the agent controlling the initiation and termination of the active phase during the bursting activity. Although these models have been superseded by more recent second-generation mod els, the methods and techniques developed remain useful, since all these subsequent models of BEA also are based on the assumption that there is a slow variable in the system that controls the initiation and termination of the active phase. From a math ematical point of view, the identity of the slow variable is irrelevant and, therefore, we believe that the analyses presented in this thesis can be generalized to other models of BEA (see below). 167 Chapter 6. Conclusions and Discussion 168 In Chapter 2, we proposed a consistent nondimensionalization of the model equations which carried all of the models into a standard form. These equations then were refor mulated to facilitate the analysis in subsequent chapters. In particular, the equations were rewritten as a second-order fast subsystem coupled to a first-order slow subsystem. Based on the form of the fast subsystem equation, we discovered that from an analytical point of view, the most important distinction between the models is the value of the integer exponent p of the activation variable n in the voltage-gated K+ current. The value of p affects the form of the damping terms (see below), however, it does not affect the qualitative bifurcation structure of the fast subsystem. The bifurcation diagrams of the fast subsystem, using the slow variable z as the bifurcation parameter, were shown to have at least one Hopf bifurcation and one corresponding homoclinic bifurcation. To determine the approximations of the plateau fractions from the model equations, we needed to be able to predict when the silent and active phases begin and end. Project ing the numerical solution of the full system of equations onto the bifurcation diagram of the fast subsystem revealed that the transition from the silent to the active phase oc curs near the local minimum of the fast subsystem nullcline, which is easy to determine. However, the transition from the active to the silent phase occurs near a homoclinic bifurcation and its location is much more difficult to determine. In Chapter 3 then, we approximated Z3c, which is the value of z at which the nu merical solution of the full system uudergoes the transition from the active phase back to the silent phase, by ZHC, which is the value of z at which the fast subsystem exhibits the homoclinic bifurcation of interest. The value of ZHG in turn was approximated by using two numerical methods and two analytically based methods which required some numerical computations. The two numerical methods were the use of AUTO to approximate the homoclinic orbit by a periodic orbit with very large period and bisection based on the bifurcation of Chapter 6. Conclusions and Discussion 169 the behaviour of the solution trajectories in the phase plane of the fast subsystem. Both methods consistently provided approximations of ZHC which were very close to typical exact Zesc values. The two analytically based methods we used were Melnikov’s method and the Fred hoim alternative method. We reformulated the fast subsystem equation as a perturbed Hamiltonian system. For perturbed Hamiltonian systems, the two methods yield the same results, since the methods produce the same sequence of improper integrals. The methods were extremely efficient for determining the leading-order approximation of ZHC since the corresponding leading-order improper integral could be converted into a line integral in the phase plane of the fast subsystem. However, they provided accurate leading-order approximations of zHC for only two of the six models under considera tion, namely the Chay-Kang and Sherman-Rinzel-Keizer models. For the remaining four models, a higher-order correction was required. The problem of obtaining a first-order correction to the leading-order approxima tions of ZHC was posed as a boundary-value problem which was solved numerically with COLSYS. The resulting first-order correction significantly improved the leading-order ap proximations of ZHC for the three models with p 1, namely the reduced Chay-Keizer, Chay (1986), and Himmel-Chay models. The first-order correction for models with p = 1, namely the Chay-Kang, Sherman-Rinzel-Keizer, and Chay-Cook models, was shown to be identically equal to zero. We expect that a few higher-order corrections to further improve the approximations of zHC for all models can be obtained in the same way. How ever, obtaining higher-order corrections is computationally very expensive and, therefore, not very practical in light of the availability of software packages such as AUTO. We concluded that for the approximation of ZHC and, subsequently, Zesc, which were required for a fixed set of fast subsystem parameters, the numerical methods are the methods of choice, especially the use of AUTO, since they efficiently provided the most Chapter 6. Conclusions and Discussion 170 accurate approximations of Zesc. Although the two analytically based methods for approximating ZHC were not as successful as the numerical methods in terms of efficiency and accuracy, they did give rise to an interesting observation. In order to apply Melnikov’s method and the Fredholm alternative method, we rewrote the fast subsystem equation of Chapter 2, (2.51), as a perturbed Hamiltonian system in Chapter 3, (3.8). The relative ease with which this reformulation could be done depended on the form of the damping terms E(u, it; z) and F(u; z) in (2.51) and especially on the value of the parameter p. The modified damping term SF(u, it; z) in (3.8) was shown to be numerically small relative to the remaining terms for all models under consideration. This result was essential for the extension of the analysis of the Sherman-Rinzel-Keizer model by Pernarowski [67] and Pernarowski et al. [70] to these models. The question remains why this modified damping term is small for all models and whether this observation is generic. The finding that this damping term is not sufficiently small for the Chay-Cook model to yield an accurate leading- order approximation of ZHC by Melnikov’s method and the Fredholm alternative method provides a starting point for future research in this direction. All four methods for approximating ZHC were based on an analysis of the fast sub system. Since the glucose-dependent parameter /9 appears only in the slow subsystem, these methods have provided approximations for zHC and, subsequently, z3 which are independent of /9. In order to elicit the dependence of z68 on 3, the effects of the slowly varying bifurcation parameter z must be included. In [79], Robinson shows how Mel nikov’s method can be extended so that it applies to perturbed systems of equations in which the parameters vary slowly with time. Robinson considers systems for which the perturbation is of the same order as the time scale on which the parameters vary. For the models of bursting electrical activity under consideration in this thesis, the parameter or slow variable z changes on a time scale of 0(E), whereas the perturbation in (3.8) is 0(6). Chapter 6. Conclusions and Discussion 171 In light of Table 3.2, 5 >> > 0, and so Robinson’s method is not directly applicable. The adaptation of Robinson’s method to the type of models studied provides another opportunity for future work. In Chapter 4, we investigated the dependence of the approximate plateau fraction on the glucose-dependent parameter 3. To determine the leading-order plateau fraction, we required the leading-order silent and active phase durations. To obtain the leading-order silent phase duration, we first rescaled the equations to the slow time scale and obtained a singular perturbation problem. The leading-order solution of this problem yielded a leading-order silent phase duration in excellent agreement with the exact silent phase duration for each biophysical model considered here. To obtain the leading-order active phase duration, we used a multiple scales analysis and averaging of the active phase problem. The resulting approximations were shown to be in excellent agreement with the exact numerical results as well. The subsequent leading-order approximations of the plateau fraction were shown to consistently overestimate the exact plateau fractions. This discrepancy is due to the omission of the duration of the transition phases in the leading-order burst period, thereby reducing the denominator in the plateau fraction. Finally, we derived the equations of two surfaces that correctly describe each of the transition phases in (u, t) space. However, the equations could not be used to derive a leading-order approximation of the transition phase durations due to the occurrence of singularities in the corresponding quadratures. We expect that to successfully obtain the leading-order transition phase durations, it is necessary to perform a higher-order analysis, requiring the method of matched asymptotic expansions. In addition, we believe that an analysis of the slow passage through a limit point (the local minimum of the fast subsystem nullcline), similar to that in [31, 40, 48], will be useful in eliciting the dependence of the duration of the transition from the silent to the active phase on the parameter and help explain the significant Chapter 6. Conclusions and Discussion 172 increase in this duration as 3 is decreased. Determining the dependence of the leading-order plateau fractions on the glucose- dependent parameter is a first step in predicting theoretical rates of insulin release from models of bursting electrical activity in pancreatic ,8-cells. In [65], Miura and Pernarowski demonstrate how leading-order plateau fractions for the Sherman-Rinzel-Keizer model can be matched with experimental plateau fractions to obtain the relationship between 9 and the external glucose concentration. Knowledge of this relationship and the rela tionship between the experimental plateau fractions and measurements of rates of insulin release from islets in turn permits the determination of theoretical rates of insulin release from the model. This may help to determine the range of glucose concentrations over which the models are valid. As a final note, we comment on the general applicability of the analytical methods presented in the first part of this thesis. The study of the Sherman-Rinzel-Keizer model by Pernarowski [67] and Pernarowski et al. [67, 70] identified the damping term of the second-order fast subsystem equation as being numerically small. It was the success of that study that led to a search for similar small terms in the corresponding second- order equations for the other models. The question remains whether the reduction that led to this equation is in fact essential for the application of Melnikov’s method or the Fredhoim alternative method to approximate zHC in Chapter 3. The answer to this question is important for the extension of these analytical methods to the study of second generation models of BEA in pancreatic 9-cells, since preliminary investigations of these models indicate that they do not immediately lend themselves to the same reduction. However, the techniques presented in Chapter 4 do not depend on the fast subsys tem being written as a perturbed Hamiltonian system. Therefore, we expect that these analytical techniques can be applied to the second-generation models. Since values of ZHC for these more recent models can be approximated numerically with AUTO, the Chapter 6. Conclusions and Discussion 173 dependence of the plateau fraction on the glucose-dependent parameter can be elicited in an efficient manner for these models as well. In the second part of this thesis, consisting of Chapter 5, we studied the bursting phenomenon in a broader context using Pernarowski’s polynomial analog model [68] of bursting electrical activity. This model exhibits a wide variety of oscillatory solution behaviour depending on the values of the model parameters, including several types of bursting behaviour similar to those observed in excitable cells other than pancreatic /3- cells. Many of the different types of solution behaviour also can be obtained by varying parameter values of biophysical models of bursting, in particular, the Chay-Cook [16] model for bursting electrical activity in pancreatic /3-cells [10]. However, the advantage of using Pernarowski’s polynomial analog [68] model rather than a biophysical model to study the different bifurcation structures of the fast subsystem underlying the various types of solution behaviour is that results can be obtained analytically. We created a bifurcation map in (ii, ‘) space, the fast subsystem parameter space, by considering the codimension-2 bifurcations associated with Hopf, homoclinic, and saddle- node of periodics bifurcations. Each codimension-2 bifurcation corresponds to a curve on the bifurcation map. These curves bound regions within which the qualitative structure of the bifurcation diagrams of the fast subsystem is the same. Moving across one of these curves on the bifurcation map corresponds, in general, to making one qualitative change in the bifurcation diagram. The existence of most curves on the bifurcation map can be inferred from a care ful investigation of bifurcation diagrams and phase portraits of the fast subsystem at various values of and ii. The numerical software packages AUTO and DSTOOL, cre ated precisely for this purpose, are invaluable in this task. For example, the existence of the NSL and CSL (right) curves, representing codimension-2 neutral homoclinic bi furcations and coalescing right homoclinic bifurcations, respectively, was inferred in this Chapter 6. Conclusions and Discussion 174 way. However, we stress the importance of analytical techniques, in particular, the use of Melnikov’s method, to infer the existence of the remaining curves. For example, the CSL (large) curve, representing coalescing large homoclinic bifurcations, was found in this way, leading to the investigation of detached branches of large periodic orbits. One observation made during the development of the bifurcation map is that the NSL (large) and CSL curves exhibit discontinuities. At these discontinuities, the curves disappear from the first quadrant of the bifurcation map and reappear in the fourth quadrant or vice versa. We note that all discontinuities occur as both q and ü tend to infinity. An asymptotic analysis should be used to further understand the meaning of the discontinuities in terms of the corresponding bifurcation diagrams. This could lead to interesting future work. Whereas (one-parameter) bifurcation diagrams are very quantitative in nature, the bifurcation map is very qualitative. In order to recover some quantitative information about the location of the codimension-1 bifurcations, we also examined two-parameter diagrams, consisting of curves of codimension-.1 bifurcations, generated with AUTO, and related them to horizontal cuts through the bifurcation map. We found that the poly nomial model exhibits three qualitatively different two-parameter bifurcation diagrams, depending on the location of the horizontal cuts relative to the points on the bifurcation map where curves meet in a cusp-like fashion, namely (ii, ii) = (0, 1) and (ii, i’t) = (, ). A careful analysis of the various bifurcations near these loci should be used to obtain a more detailed understanding of the transitions between the bifurcation diagrams in neighbouring regions. Surprisingly, we discovered that the two-parameter bifurcation diagram generated with the polynomial analog model for % = 1.5 is essentially identical to the diagram generated with the biophysically based Chay-Cook model in [10]. That is, the polynomial analog model not only mimics the wealth of oscillatory behaviour of models of excitable Chapter 6. Conclusions and Discussion 175 cells, but it also possesses the same rich bifurcation structure of such models. This raises the question whether this diagram is generic for models of bursting oscillations. A related question is whether the bifurcation map also is generic. In [10], Bertram et al. suggest that the two-parameter bifurcation diagram occurs often in models of neuronal spiking and that the underlying cause is the existence of a codimension-3 degenerate Takens-Bogdanov bifurcation of focus type [29]. Since these models typically have a large number of parameters, they should have sufficient flexibility to exhibit this type of bifurcation and, hence, at least support the types of bursting implied by these diagrams. In light of the qualitatively different two-parameter bifurcation diagrams obtained with the polynomial analog model for 2 = 0.8 and i = 0, we suggest that the classifi cation scheme of bursting oscillations of Bertram et al. [10] can be taken a step further by considering bifurcation maps. In particular, these maps provide information about the relationships among the various types of oscillatory behaviour as evident from the bifurcation diagrams in the different regions. Finally, we showed that the size of the asymptotic parameter can alter the solution behaviour of the full system of equations near boundaries of regions on the bifurcation map. Specifically, we showed that a tapered burster can be made to appear to be a square wave burster by increasing sufficiently. We believe that this non-generic behaviour is typical near boundaries of regions. In conclusion, the research presented in this thesis extends the mathematical frame work with which to study a large class of bursting oscillators. The analysis adds to the understanding of specific models of bursting electrical activity in pancreatic ,8-cells and bursting phenomena in general. Bibliography [1] B. ALBERTS, D. BRAY, J. LEwIs, M. RAFF, K. ROBERTS, AND J. D. WATSON, Molecular Biology of the Cell, Garland Pubi., New York, 2nd ed., 1989. [2] U. ASCHER, J. CHRISTIANSEN, AND R. D. RUSSELL, Collocation software for boundary-value ODEs, ACM Trans. Math Software, 7 (1981), pp. 209—222. [3] F. M. ASHCROFT, D. E. HARRISON, AND S. J. H. ASHCROFT, Glucose induces closure of single potassium channels in isolated rat pancreatic /3-cells, Nature, 312 (1984), pp. 446—448. [4] S. J. H. ASHCROFT, J. BASSETT, AND P. J. RANDLE, Insulin secretion mech anisms and glucose metabolism in isolated islets, Diabetes, 21 (Suppi. 2) (1972), pp. 538—545. [5] I. ATWATER, C. M. DAWSON, B. RIBALET, AND E. ROJAS, Potassium perme ability activated by intracellular calcium ion concentration in the pancreatic /3-cell, J. Physiol. (London), 288 (1979), pp. 575—588. [6] I. ATWATER, C. M. DAWSON, A. SCOTT, G. EDDLESTONE, AND E. ROJAS, The nature of oscillatory behaviour in electrical activity from pancreatic /3-cell, in Biochemistry and Biophysics of the Pancreatic /3-cell, W. J. Malaisse and I. B. Täljedal, eds., Georg Thieme Verlag, Stuttgart, 1980, pp. 100—107. Hormone and Metabolic Research Supplement Series No. 10. [7] I. ATWATER, B. RIBALET, AND E. ROJAS, Mouse pancreatic cells: tetraethyl blockage of the potassium permeability increase induced by depolarization, J. Physiol. (London), 288 (1979), pp. 561—574. [8] P. M. BEIGELMAN, B. RIBALET, AND I. ATWATER, Electrical activity of mouse pancreatic beta-cells, J. Physiol. (Paris), 73 (1977), pp. 201—217. [9] C. M. BENDER AND S. A. ORszAG, Advanced Mathematical Methods for Scientists and Engineers, McGraw-Hill, New York, 1978. [10] R. BERTRAM, M. J. BUTTE, T. KIEMEL, AND A. SHERMAN, Topological and phenomenological classification of bursting oscillations, Bull. Math. Biol., 57 (1995), pp. 413—439. 176 Bibliography 177 [11] R. BERTRAM, P. SMOLEN, A. SHERMAN, D. MEARS, I. ATWATER, F. MARTIN, AND B. S0RIA, A role for calcium release activated current (crac) in cholinergic modulation of electrical activity in pancreatic /3-cells, Biophys. J., (1995). In press. [12] W. E. BoYcE AND R. C. DIPRIMA, Elementary Differential Equations, Wiley, New York, 3rd ed., 1977. [13] T. R. CHAY, On the effect of the intracellular calcium-sensitive K channel in the bursting pancreatic /3-cell, Biophys. J., 50 (1986), pp. 765—777. [14] , The effect of inactivation of calcium channels by intracellular Ca2+ ions in the bursting pancreatic /3-cells, Cell Biophysics, 11(1987), pp. 77—90. [15] , Effect of compartmentalized Ca2 ions on electrical bursting activity of pan creatic /3-cells, Am. J. Physiol., 258 (1990), pp. C955—C965. [16] T. R. CRAY AND D. L. COOK, Endegenous bursting patterns in excitable cells, Mathematical Biosciences, 90 (1988), pp. 139—153. [17] T. R. CHAY AND H. S. KANG, Multiple oscillatory states and chaos in the en dogenous activity of excitable cells: pancreatic /3-cell as an example, in Chaos in Biological Systems, H. Degn, A. V. Holden, and L. F. Olsen, eds., Plenum Press, 1987, pp. 173—181. [18] T. R. CRAY AND J. KEIZER, Minimal model for membrane oscillations in the pancreatic /3—cell, Biophys. J., 42 (1983), pp. 181—190. [19] , Theory of the effect of extracellular potassium on oscillations in the pancreatic /3-cell, Biophys. J., 48 (1985), pp. 815—827. [20] S.-N. CHOW, J. K. HALE, AND J. MALLET-PARET, An example of bifurcation to homoclinic orbits, J. Differential Equations, 37 (1980), pp. 351—373. [21] D. L. COOK AND C. N. HALEs, Intracellular ATP directly blocks K channels in pancreatic B-cells, Nature, 311 (1984), pp. 271—273. [22] D. L. CooK, M. IKEUCHI, AND W. Y. FUJIMOTO, Lowering of pH inhibits Ca2+ activated K channels in pancreatic B-cells, Nature, 311 (1984), pp. 269—271. [23] D. L. COOK, L. S. SATIN, AND W. F. HOPKINS, Pancreatic B-cells are bursting, but how?, TINS, 14 (1991), pp. 411—414. [24] P. M. DEAN AND E. K. MATTHEWS, Electrical activity in pancreatic islet cells, Nature, 219 (1968), pp. 389—390. Bibliography 178 [25] , Glucose-induced electrical activity in pancreatic islet cells, J. Physiol. (Lon don), 210 (1970), pp. 255—264. [26] 0. DECROLY AND A. GOLDBETER, From simple to complex oscillatory behavior: analysis of bursting in a multiply regulated biochemical system, J. Theor. Biol., 124 (1987), pp. 219—250. [27] M. DESCHENES, J. P. Roy, AND M. STERIADE, Thalamic bursting mechanism: an inward slow current revealed by membrane hyperpolarization, Brain Res., 239 (1982), pp. 289—293. [28] E. DOEDEL, AUTO: A program for the automatic bifurcation analysis of au tonomous systems, Cong. Num., 30 (1981), pp. 265—284. [29] F. DUMORTIER, R. ROUSSARIE, AND J. SOTOMAYOR, Generic 3-parameter fami lies of planar vector fields, unfoldings of saddle, focus and elliptic singularities with nilpotent linear parts, in Bifurcations of Planar Vector Fields: Nilpotent Singularities and Abelian Integrals, F. Dumortier, R. Roussarie, J. Sotomayor, and H. Zoladek, eds., Springer-Verlag, Berlin, 1991, pp. 1—164. Lecture Notes in Mathematics, Vol. 1480. [30] G. B. ERMENTROUT AND N. K0PELL, Parabolic bursting in an excitable system coupled with a slow oscillation, SIAM J. Appi. Math., 46 (1986), pp. 233—253. [31] T. ERNEUX AND J. P. LAPLANTE, Jump transion due to a time-dependent bifur cation parameter in the bistable iodate-arsenous acid reaction, J. Chem. Phys., 90 (1989), pp. 6129—6134. [32] I. PINDLAY, M. J. DUNNE, AND 0. H. PETERSEN, ATP-sensitive inward rectifier and voltage- and calcium-activated K+ channels in cultured pancreatic islet cells, J. Membrane Biol., 88 (1985), pp. 165—172. [33] , High-conductance K+ channels in pancreatic islet cells can be activated and inactivated by internal calcium, J. Membrane Biol., 83 (1985), pp. 169—175. [34] R. FITZHUGH, Impulses and physiological states in theoretical models of nerve mem brane, Biophys. J., 1 (1961), pp. 445—466. [35] J. GUCKENHEIMER, Multiple bifurcation problems of codimension two, SIAM J. Math. Anal., 15 (1984), pp. 1—49. [36] J. GUCKENHEIMER, Multiple bifurcation problems for chemical reactors, Physica D, 20 (1986), pp. 1—20. Bibliography 179 [37] J. GUCKENHEIMER AND P. HoLMEs, Nonlinear Oscillations, Dynamical Sys tems, and Bifurcations of Vector Fields, vol. 42 of Applied Mathematical Sciences, Springer-Verlag, New York, 1983. [38] J. GUCKENHEIMER AND I. S. LABOURIAU, Bifurcations of the Hodgkin and Huxley equations: a new twist, Bull. Math. Biol., 55 (1993), pp. 937—952. [39] J. GUCKENHEIMER, M. R. MYERS, F. J. WIcKLIN, AND P. A. WORFOLK, DSTOOL: A Dynamical System Toolkit with an Interactive Graphical Interface, Cen ter for Applied Mathematics, Cornell University, Ithica, New York, 1994. [40] R. HABERMAN, Slowly-varying jump and transition phenomena associated with al gebraic bifurcation problems, SIAM J. Appl. Math., 37 (1979), pp. 69—106. [41] R. M. HARRIS-WARRICK AND R. E. FLAMM, Multiple mechanism of bursting in a conditional bursting neuron, J. Neurosci., 7 (1987), pp. 2113—2128. [42] J. C. HENQuIN, Tetraethylammonium poteritiation of insulin and inhibition of ru bidium efflux in pancreatic islets, Biochem. Biophys. Res. Commun., 77 (1977), pp. 551—556. [43] , The potassium permeability of pancreatic islet cells: mechanism of control and influence on insulin release, in Biochemistry and Biophysics of the Pancreatic 3- cell, W. J. Malaisse and I. B. Täljedal, eds., Georg Thieme Verlag, New York, 1980, pp. 66—73. Hormone and Metabolic Research Supplement Series No. 10. [44] D. M. HIMMEL AND T. R. CHAY, Theoretical studies on the electrical activity of pancreatic /3-cells as a function of glucose, Biophys. J., 51(1987), pp. 89—107. [45] A. C. HINDMARSH, ODEFACK, a systemized collection of ODE solvers, in Scientific Computing, R. S. Stepleman et al., eds., North-Holland, Amsterdam, 1983, pp. 55— 64. [46] M. W. HIRSCH AND S. SMALE, Differential Equations, Dynamical Systems, and Linear Algebra, Academic Press, New York, 1974. [47] A. L. HODGKIN AND A. F. HUXLEY, A quantitative description of membrane cur rent and its application to conduction and excitation in nerve, J. Physiol. (London), 117 (1952), pp. 500—544. [48] L. H0LDEN AND T. ERNEUX, Understanding bursting oscillations as periodic slow passages through bifurcation and limit points, J. Math. Biology, 31(1993), pp. 351— 365. Bibliography 180 [49] J. L. HuDSON, M. HART, AND D. MARINKO, An experimental study of multiple peak periodic and nonperiodic oscillations in the Belousov-Zhabotinskii reaction, J. Chem. Phys., 71(1979), PP. 1601—1606. [50] E. KANDEL AND W. SPENCER, Electrophysiology of hippocampal neurons II. after potentials and repetitive firing, J. Neurophysiol., 24 (1961), pp. 243—259. [51] T. KAPER, 1994. Personal communication. [52] J. P. KEENER, Principles of Applied Mathematics: Transformation and Approxi mation, Addison-Wesley, New York, 1988. [53] J. KEIZER, Electrical activity and insulin release in pancreatic beta cells, Math. Biosciences, 90 (1988), Pp. 127—138. [54] J. KEIZER AND G. MAGNUS, ATP-sensitive potassium channel and bursting in the pancreatic beta cell, Biophys. J., 56 (1989), pp. 229—242. [55] J. KEIZER AND P. SMOLEN, Bursting electrical activity in pancreatic /3-cells caused by Ca2- and voltage-inactivated Ca2 channels, Proc. Nati. Acad. Sci., 88 (1991), Pp. 3897—3901. [56] J. KEVORKIAN AND J. D. COLE, Perturbation Methods in Applied Mathematics, vol. 34 of Appi. Math. Sci., Springer, New York, 1981. [57] G. E. KUZMAK, Asymptotic solutions of nonlinear second order differential equa tions with variable coefficients, P.M.M., 23 (1959), PP. 515—526. [58] A. J. LICHTENBERG AND M. A. LIEBERMAN, Regular and Stochastic Motion, vol. 38 of Applied Mathematical Sciences, Springer-Verlag, New York, 1983. [59] J. C. LUKE, A perturbation method for nonlinear dispersive wave problems, Proc. Royal Soc. London, Ser. A, 292 (1966), pp. 403—412. [60] E. K. MATTHEWS AND Y. SAKAMOTO, Electrical characteristics of pancreatic islet cells, J. Physiol. (London), 246 (1975), pp. 421—437. [61] H. P. MEISSNER, Electrical characteristics of the beta cells in pancreatic islets, J. Physiol. (Paris), 72 (1976), pp. 757—767. [62] H. P. MEISSNER AND M. PREISSLER, Ionic mechanisms of the glucose-induces membrane potential changes in B-cells, in Biochemistry and Biophysics of the Pancreatic 3-cell, W. J. Malaisse and I. B. Täljedal, eds., Georg Thieme Verlag, Stuttgart, 1980. Hormone and Metabolic Research Supplement Series No. 10. Bibliography 181 [63] H. P. MEIssNER AND H. SCHMELZ, Membrane potential of beta-cells in pancreatic islets, Pflugers Arch., 351 (1974), pp. 195—206. [64] V. MEL’NIIcov, On the stability of the center for time-periodic perturbations, Trans. Moscow Math. Soc., 12 (1963), pp. 1—57. [65] R. M. MIURA AND M. PERNAROWSKI, Correlations of rates of insulin release from islets and plateau fractions for /3-cells, J. Math. Biology, 57 (1995), pp. 229—246. [66] S. OZAWA AND 0. SAND, Electrophysiology of endocrine cells, Physiol. Rev., 66 (1986), pp. 887—952. [67] M. PERNAROWSKI, The mathematical analysis of bursting electrical activity in pan creatic beta cells, PhD thesis, University of Washington, Seattle, 1990. [68] , Fast subsystem bifurcations in a slowly varying Lie’nard system exhibiting burst ing, SIAM J. Appi. Math., 54 (1994), pp. 814—832. [69] M. PERNAROWSKI, R. M. MIURA, AND J. KEVORKIAN, The Sherman-Rinzel Keizer model for bursting electrical activity in the pancreatic /3-cell, in Differen tial Equations Models in Biology, Epidemiology and Ecology, S. Busenberg and M. Martelli, eds., Springer-Verlag, New York, 1991, pp. 34—53. Lecture Notes in Biomathematics, Vol. 92. [70] , Perturbation techniques for models of bursting electrical activity in pancreatic beta-cells, SIAM J. Appi. Math., 52 (1992), pp. 1627—1650. [71] R. PIESSENS, E. DEDONCKER KAPENGA, C. UBERHUBER, AND D. KAHANER, Quadpack: A Subroutine Packag.e for Automatic Integration, vol. 1 of Series in Com putational Mathematics, Springer-Verlag, 1983. [72] R. E. PLANT, The effects of calcium on bursting neurons, Biophys. J., 21(1978), pp. 217—237. [73] R. E. PLANT AND M. KIM, Mathematical description of a bursting pacemaker neuron by a modification of the Hodgkin-Huxley equations, Biophys. J., 16 (1976), pp. 227—244. [74] W. H. PRESS, B. P. FLANNERY, S. A. TEUKOLSKY, AND W. T. VETTERLING, Numerical Recipes (Fortran version), Cambridge University Press, Cambridge, 1989. [75] B. RIBALET AND P. M. BEIGELMAN, Cyclic variation of K+ conductance in pan creatic /3-cells: Ca2 and voltage dependence, Am. J. Physiol., 237 (1979), pp. C137— C 146. Bibliography 182 [76] J. RINzEL, Bursting oscillations in an excitable membrane model, in Ordinary and Partial Differential Equations, B. D. Sleeman and R. J. Jarvis, eds., Springer, New York, 1985, PP. 304—316. Lecture Notes in Mathematics, Vol. 1151. [77] J. RINZEL, A formal classification of bursting mechanisms in excitable systems, in Mathematical Topics in Population Biology, Morphogenesis, and Neurosciences, E. Teramoto and M. Yamaguti, eds., Springer-Verlag, New York, 1987, pp. 267—281. Lecture Notes in Biomathematics, Vol. 71. [78] J. RINzEL AND Y. S. LEE, On different mechanisms for membrane potential burst ing, in Nonlinear Oscillations in Biology and Chemistry, H. G. Othmer, ed., Springer- Verlag, New York, 1986, pp. 19—33. Lecture Notes in Biomathematics, Vol. 66. [79] C. ROBINSON, Sustained resonance for a nonlinear system with slowly varying co efficients, SIAM J. Math. Anal., 14 (1983), pp. 847—860. [80] P. RORSMAN AND G. TRUBE, Calcium and delayed potassium currents in mouse pancreatic /3-cells under voltage-clamp conditions, J. Physiol. (London), 375 (1986), pp. 531—550. [81] R. SANTOS, L. RosARlo, A. NADAL, J. GARCIA-SANCHO, B. S0RIA, AND M. VALDEOLMILLOS, Widespread synchronous Ca2+ oscillations due to bursting electrical activity in single pancreatic islets, Pflügers Arch., 418 (1991), pp. 417—422. [82] L. S. SATIN AND D. L. CooK, Voltage-gated Ca2+ current in pancreatic B-cells, Pflügers Arch., 404 (1985), pp. 385—387. [83] , Evidence for two calcium currents in insulin-secreting cells, Pflügers Arch., 411 (1988), pp. 401—409. [84] , Calcium current inactivation in insulin-secreting cells is mediated by calcium influx and membrane depolarization, Pffügers Arch., 414 (1989), pp. 1—10. [85] S. SCHECTER, The saddle node separatrix loop, SIAM J. Math. Anal., 18 (1987), Pp. 1142—1156. [86] A. SHERMAN, Anti-phase, asymmetric and aperiodic oscillations in excitable cells - I. Coupled bursters, Bull. Math. Biol., 56 (1994), pp. 811—835. [87] A. SHERMAN, J. RINZEL, AND J. KEIZER, Emergence of organized bursting in clusters of pancreatic /3-cells by channel sharing, Biophys. J., 54 (1988), pp. 411— 425. Bibliography 183 [88] P. SM0LEN AND J. KEIZER, Slow voltage inactivation of Ca2+ currents and bursting mechanisms for the mouse pancreatic /3-cell, J. Membrane Biol., 127 (1992), pp. 9— 19. [89] P. SMOLEN, J. RINzEL, AND A. SHERMAN, T’Vhy pancreatic islets burst but single /3-cells do not: the heterogeneity hypothesis, Biophys. J., 64 (1993), pp. 1668—1680. [90] D. TERMAN, Chaotic spikes arising from a model of bursting in excitable mem branes, SIAM J. Appi. Math., 51(1991), pp. 1418—1450. [91] M. VALDEOLMILLOS, R. M. SANTOS, D. CONTRERAS, B. S0RIA, AND L. M. RosARTo, Glucose-induced oscillations of intracellular Ca2+ concentration resem bling bursting electrical activity in single mouse islets of Langerhans, FEBS Letters, 259 (1989), pp. 19—23. [92] S. WIGGINS, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer-Verlag, New York, 1990. [93] R. S. K. W0NG AND D. A. PRINCE, Afterpotential generation in hippocampal pyramidal cells, J. Neurophysiol., 45 (1981), pp. 86—97. Appendix A Equations and Parameter Values of the First Generation Pancreatic Beta Cell Models The following sections contain the equations and parameter values of the first-generation models which are analyzed in this thesis, as well as the details of the nondimensionaliza tion used for each model. A.1 The Reduced Chay-Keizer Model The dimensional reduced Chay-Keizer model [19] is as follows: 4r2Cm = — [10a(V) + IK(V) + IK,Ca( V, Ca) + IL(V)], c/n — n(V)—n dT — r(V) dCa —3 dT 4irrF’07 — kcaCat where Ia(V) = gcam(V)hoo(V)(V — Vga), IK(V) = .qKn(V—VK), JK(v,ca)(V,Gai) = KCa (V — VK), IL(V) = .qL(V—VL), and Bm(Vm - V) m(V) = —1 184 Appendix A. Model Equations 185 /3m(V) = 4.Oexp [(Va - v)J, 0m(V) m00 — fT7 1? (T7\’V ) i Pm V) ah(V) = O.O7exp [-—(Vh — h(V) = ____________ exP[-(Vb_V)] +1 h v ah(V) - B(V-V) = 1 firexp SVn_ — /3(V) O.l2Sexp — crn(V) ?2 — r(V) o(V) + ,B(V) The dimensions of all variables and parameters, as well as the values of the parameters are given in Table A.1. We note that instead of the A used in [19], we use r = A’. In [19], the value of A’ is given to be 3. In this case, the bifurcation diagram of the fast subsystem is that of a parabolic-amplitude burster (cf. Chapter 5). For consistency among the first-generation models considered in this thesis, we modify A’ to be 3.25, so that the bifurcation diagram is that of a square-wave burster. The qualitative behaviour of the numerical solution of the full system of equations is not affected significantly. The nondimensionalization of the model equations is as described in Section 2.3, with X being equal to the dissociation constant Kd. The value of is taken to be 18 msec. The resulting nondimensional model, written in standard form, then is dv . 4 = ia(V)w (v+1)—g(z)(v+1)+l(v), dw ü3(v)—w cit — (v) Appendix A. Model Equations 186 dz = [/3e_zica(v) — 11, where a(V) 7Camo(v)hoo(V)(VCaV), ez g(z) = 7KCai+z 1(v) = 7L(vL—v), and äm(V) = bm(Vmv) exp{—(vm—v)] 1’ ri m(V) 4.Oexp [_(va — ãm(V) = ri ah(v) = 0.07 exp L—(vh — v)],Sh 1 /3h (v) exp{(vb_v)] +1’ L(v) ah(v) Lh(V) +/3h(v) ä(v) b(v — v) exp[(vn_v)] —1’ = O.l25exp [-_(vc — Sc — _________ — A = ä(v) + (v) The valnes of the parameters are given in Table A.1. Appendix A. Model Equations 187 Table A.1: Variables and Chay-Keizer model. parameters of the dimensional and nondimensional reduced Dimensional Model Nondimensional Model Variable Dimension Variable V mV v Ti - W Ca z T msec t Parameter Dimension Value Parameter Value gca gK gK(v,ca) gL VCa VK VL Vcz Vb vc Vh Vm vTh Sa Sb Sc Sm Sn Bm Bn Cm ka Kd f r T pS PS pS PS mV mV mV mV mV mV mV mV mV mV mV mV mV mV mV 1/mV 1/mV 107F/cm2 msec1 msec 7Ca 7K 7K(V,Ca) 7L UCa VL Va Vb vc Vh Vm vn Sc Sn bm b /3 8000 7500 50 30 100 -75 -40 -50 -20 -30 -50 -25 -20 18 10 80 20 10 10 0.1 0.01 10 0.01 1 0.007 6 3.25 31.8310 29.8416 0. 1989 0.1194 1.3333 -0.5333 -0.6667 -0.2667 -0.4 -0.6667 -0.3333 -0.2667 0.24 0. 1333 1.0667 0.2667 0. 1333 0. 1333 7.5 0.75 21.5919 1 .26e3 0. 1806 Appendix A. Model Equations 188 A.2 The Chay (1986) Model The dimensional Chay (1986) model [13] is as follows: 4r2Cm = — [1c0(V, Ca) + IK(V,Ca)( V, Ca) + IL(V)j, c/n — n(V,Ca)—n dT — r(V,Ca1 dCa1 = [4rF0’ Cad) — kcaCai], where ICa(V, Ca) = gcam(V) (V — Vca(Caj)), IK(V,Ca)(V, Ca) = gK(V,Ca)fl(V — VK), IL(V) = gL(V—VL), and R’T Ca0Vca(Cai) = Vc(Ca) = Aln(i), Bm(Vm - V)0m(V) exp {(Vm — v)] — 1 13m(V) = 4.Oexp [(Va - am(V) m(V) am(V) + m(V)’ c(V, Ca:) B(V + Vc(Ca) — V) exp [(v + Vc(Ca) — V)] —1 /3,(V,Ca) = O.l25exp [-(Vb + Vc(Ca1)— c(V, Ca:) n(V, Ca) = a(V, Ca) + Ca)’ r(V, Ca) a(V, Ca) + Cad) Appendix A. Model Equations 189 The dimensions of all variables and parameters, as well as the values of the parameters are given in Table A.2. We note that instead of the A used in [13], we use r = A’. In [13], the value of A’ is given to be 1.35. Although this value does give rise to square-wave bursting, the bifurcation diagram of the fast subsystem of this model is that of a tapered burster (cf. Chapter 5). For consistency among the first-generation models considered in this thesis, we modify A’ to be 1.4, so that the bifurcation diagram of the fast subsystem is that of a square-wave burster. The qualitative behaviour of the numerical solution of the full system of equations is not affected. The nondimensionalization of the model equations is as described in Section 2.3, with X being equal to Ca0, 7K replaced by 7K(V,C), and the definition of three additional nondimensional parameters as follows: -A = -A L/ — VK R’T = 2FVK The value of is taken to be 8 msec. The resulting nondimensional model, written in standard form, then is dv ica(v,z)—w2(v+1)+l(v), dw — _______ dl — (v,z) = E [/3eica(v, z) — 1], where = 7Camo(v)(11ZV), 1(v) = 7L(VLV), Appendix A. Model Equations 190 and = ñ(v) (v,z) = i3(v,z) = ü(v,z) (v,z) = The values of the parameters are bm(Vm — v) exp[i(vm_v)] 1’ 4.Oexp — v)] m(V) m() + m(V)’ b(v + [L + ziz — v) exp[-(vn+i+vz — v)] —1’ ri0.l25exp + ,u + vz — v) Sb 7K(v,Ca)U)(V z) (v,z)+ ,(v,z)’ A ä(v,z)+ $(v,z) given in Table A.2. Appendix A. Model Equations 191 Dimensional Model Nondimensiorial Model Variable ] Dimension Variable V mV v n - w Ca z T msec t Parameter 1_Dimension Value Parameter Value .Ca pS 250 7Ga 0.4421 gK(V,Ca) pS 1500 7K(V,Ca) 2.6526 9L S 10 7L 1.7683e Vi- mV -87 VL mV -45 V -0.5172 Va mV 55 Va -0.6322 Vb mV -35 V -0.4023 Vm mV 30 Vm 0.3448 V, mV -25 v -0.2874 Sa mV 18 a 0.2068 Sb mV 80 0.9195 Sm mV 10 ‘m 0.1149 S, mV 10 s 0.1149 Bm 1/mV 0.1 bm 8.7 B 1/mV 0.01 b 0.87 Cm 107F/cm2 10 ka msec1 0.03 7.5140e3 f - 0.0001 2.4e5 r itm 6 r msec 1.40 0.175 Ca0 2500 T K 310 A mV -50 -4.4968 ii -0.5747 __________ ___________ _______ 77 -0.1535 Table A.2: Variables and parameters of the dimensional and the nondimerisional Chay (1986) model. Appendix A. Model Equations 192 A.3 The Himmel-Chay Model The dimensional Himmel-Chay model [44] is as follows: 4r2Cm = — [Ia(V, Cad) + IK(V) + IK,Ca( V, Cad) + IL(V)], dn — n(V)—n dT — r(V) dCa [4rF0 — kcaCai], where 1Ca(V, Ca:) gcam(V)h(V) (V — Vca(Caj)), IK(V) = gKfl3(VVK), IK,Ca(V,Cai) — VK), IL(V) = gL(V — VL), and RT Ca0 Vca(Caj) = Bm(Vm - V)0m(V) exp [(Vm — v)] — 1 /3m(V) = 4.Oexp[-(Va_V)], am(V) m(V) = h(V) = O.O7exp [i(Vh — h(V) exp {-(Vb — V)] +1 h 1v’ ah(V)oo — h(V)+/3h(V) Appendix A. Model Equations 193 = B(V-V) ‘Iv Vexp s—’ — i3(V) = O.l25exp — V)], c(V) I — r(V) = a(V) + (V) The dimensions of all variables and parameters, as well as the values of the parameters are given in Table A.3. We note that in [44], the value of VL is misprinted as 40 mV. We use VL = —40 mV here. Furthermore, ka is given to be 0.045 msec’ in the appendix of [44]. However, this value gives rise to continuous bursting. Different values of ka, varying between 0.01 and 0.045, are used in the text of [44] and we use k0a = 0.03 msec’ here. Finally, instead of the A used in [44], we use r = A. In [44], the value of A is misprinted as 1/0.003. We use the value of A given in [19], namely A = 1/3, giving r 3. The nondimensionalization of the model equations is as described in Section 2.3, with X being equal to the dissociation constant Kd, and the definition of two additional nondimensional parameters as follows: R’T Kd It = 2FVKCcLQ’ RT = 2FVK The value of r, is taken to be 18 msec. The resulting nondimensional model, written in standard form, then is dv . 3 = zc(v,z)—w(v+1)—g(z)(v+1)+l(v), dw — ______ cit — (v) dz = e[/3e_zzca(v,z)_1], Appendix A. Model Equations 194 where a(V,Z) = 7cafizo(v)ihoo(v)(1t+vzv), g(z) = 7KCal+z and äm(V) bm(VmV) exp[—(vm—v)j 1’ ri = 4.Oexp — v)], àm (v) — ri ah(v) = 0.07 exp — 1 13h (v) exp [*(vb _v)j + 1’ h(V) hc(v) = ah(v) + 13h(v)’ — v) exp[i-(vn—v)] —1’ = 0.l25exp — Sc 7w(v) = = (v) +$(v) The values of the parameters are given in Table A.3. Appendix A. Model Equations 195 Dimensional Model Nondimensional Model Variable Dimension Variable V mV v n - w Ca z T msec t Parameter Dimension Value Parameter Value gCa PS 7000 7Ca 27.8521 gK p5 2000 7K 7.9577 gK,Ca pS 60 7K,Ga 0.2387 gL pS 28 7L 0.1114 V- mV -87 VL mV -40 VL 0.4598 Va mV 50 Va -0.5747 Vb mV -20 0b -0.2299 V mV -30 Vc 0.3448 Vh mV -50 Vh -0.5747 Vm mV -25 Vm 0.2874 V mV -20 vn 0.2299 Sa mV 18 a 0.2069 Sb mV 10 Sb 0.1149 Sc mV 80 Sc 0.9195 Sh mV 20 0.2299 Sm mV 10 m 0.1149 S mV 10 Sn 0.1149 Bm 1/mV 0.1 bm 8.7 Bn 1/mV 0.01 bn 0.87 Cm 107F/cm2 10 ka msec1 0.03 /9 4.1744 Kd 2 f - 0.002 1.08e3 r 6 r msec 3 0.1667 Ca0 2000 K 310 1.0605 __________ ___________ ______ ii -0.1535 Table A.3: Variables and parameters of the dimensional and the nondimensional Him mel-Chay model. Appendix A. Model Equations 196 A.4 The Chay-Kang Model The dimensional Chay-Kang model [17] is as follows: 4r2Cm = — [ICa(V,Ca)(V, Ca) + IK(V) + IL(V)], dn — n(V)—n dT r(V) dCa: _____ dT 4rrF a)O” Car) — kcaCaj where ‘Ca(V,Ca)(V, Cad) = gCa(V,Ca)moo(V)hoo(V, Ca)(V — Vga), IK(V) = gKfl(VVK), IL(V) = gL(V—VL), and m(V) 1 1 +exp [—(Vm — V)] h(V,Ca) 1 r 1 Ca i—v+ —--exp -- n(V) 1 1+ exp — V)] r(V) T 1+ exp — V)] The dimensions of all variables and parameters, as well as the values of the parameters are given in Table A.4. We note that instead of the A used in [17], we use r = In [17], the value of A is misprinted as 0.017 msec’. The value used in the companion paper [14] is 0.117 msec’, giving r 8.547 msec. The nondimensionalization of the model equations is as described in Section 2.3, with X being equal to Kh. The value of is taken to be 9 msec. The resulting nondimensional Appendix A. Model Equations 197 model, written in standard form, then is dv = ica(v,z)w(v+1)+l(v), dw — zZ(v)—w dt — dz —- = e[/3e_zzca(v, z) — 1], where ia(v, z) = 7Ca(V,Ca)moo(V)Iloo(v, Z)(Va — v), 1(v) = 7L(VL—v), and 1 1 + exp [—(vm — — 1h(v,z) = — 1 + ezexp [] - 7K WV) 1+exp [—(v—v)] 1+exp [zi(v_v)] The values of the parameters are given in Table A.4. Appendix A. Model Equations 198 Dimensional Model Nondimensional Model Variable Dimension Variable V mV v 72 - W Ca: Z T msec t [ Parameter Dimension_[_Value Parameter Value gCa(V,Ca) PS 250 7Ca(V,Ca) 0.4974 gK PS 1300 7K 2.5863 gL PS 10 7L 1.9894e VCa mV 100 VGa 1.25 VK mV -80 VL mV -40 VL -0.5 Vm mV -22 Vm 0.275 Vn mV -9 v -0.1125 Sh mV 10 0.125 Sm mV 7.5 Sm 9.375e2 S mV 10 Sn 0.125 Cm 107F/cm2 10 ka msec1 0.06 8.5301e2 Kh 90 f - 0.001 5.4e4 r 6 r msec 8.547 A 0.9497 Table A.4: Variables and parameters of the dimensional and the nondimensional Chay-Kang model. Appendix A. Model Equations 199 A.5 The Sherman-Rinzel-Keizer Model The dimensional Sherman-Rinzel-Keizer model [87] is as follows: 4r2Cm = — [Ia(V) + IK(V) + ‘K,Ca(V, Ca)], dn — n(V)—n ciT — r(V) dCa —3 dT = 4KrF — kcaCaj where 10a(V) = gcamoo(V)hc.o(V)(V—Vca), IK(V) = gKn(V—VK), IK,Ca(V,Gai) = K,caKa(V_VK) and 1 m(V) 1+exp[_-(Vm _v)] h(V) 1 1 +exp — V)] nco(V) 1 + exp - V)] r(V) exp - V)] + exp [-(Vb - V)] The dimensions of all variables and parameters, as well as the values of the parameters are given in Table A.5. We note that the parameter r used here is equivalent to the c/A used in [87]. In [87], c is given to be 60 msec and A is either 1.6 or 1.7. We use A = 1.6, • • 60giving r = = 37.5 msec. Appendix A. Model Equations . 200 The nondimensionalization of the model equations is as described in Section 2.3, with X being equal to the dissociation constant Kd. The value of is taken to be 22 msec. The resulting nondimensional model, written in standard form, then is dv = zca(v)—w(v+l)—g(z)(v+1), dw — _____ dt — dz —— = e[/3e_zzca(v) — 11, where = ez g(z) = 7KCai+Z and 1 m(v) = 1+exp[—(vm—v)] — 1 h(v) 1 +exp [-(vh—v)] - 7K w(v) = 1 +exp [—(v _v)j A r(v) = exp [(vb — v)] + exp [*vb — v)] The values of the parameters are given in Table A.5. Appendix A. Model Equations 201 Dimensional Model Nondimensional Model Variable Dimension Variable V mV v n - w Ca: 1IM Z T msec t Parameter Dimension Value Parameter Value gCa PS 1400 7Ca 5.8012 gK p5 2500 “fK 10.3592 gK,Ca P5 30000 7K,Ca 124.3104 VCa mV 110 VCa 1.4667 VK mV -75 Vb mV -75 -1 Vh mV -10 Vh -0.1333 Vm mV 4 V, 5.3333e2 V mV -15 v -0.2 Sa mV 65 0.8667 Sb mV 20 0.2667 Sh mV 10 0.1333 Sm mV 14 m 0.1867 S mV 5.6 7.4667e2 Cm 107F/cm2 10 ka msec1 0.06 2.7179e Kd 100 f - 0.0005 6.6e4 r 6.5 r msec 37.5 1.7045 Table A.5: Variables and parameters of the dimensional and the nondimerisional Sher man-Rinzel-Keizer model. Appendix A. Model Equations 202 A.6 The Chay-Cook Model The dimensional Chay-Cook model [16] is as follows: where 2 dV4irr Cm dn dT dCa dT = — [Ia(V) + ‘Ca( V,Ca)( V, Ca) + IK(V) + IL(V)], — n(V)—n — r(V) = f + ICa(V,Ga)( V, Ca)) — kcaCai], 1a(V) 1Ca(V,Oa)O”, Cat) IK(V) IL(V) = .crnoo(V)(V — Vga), = gCa(V,Ca)scc(V, Ca)(V — Vga), = gKn(V—VK), = gL(V—VL), 1 1 +exp — V)] 1 + exp (-(v8 - V)) and m00(V) s(V, Caj n(V) l+exp r(V) = T 1+ exp — V)] The dimensions of all variables and parameters, as well as the values of the parameters are given in Table A.6. The nondimensionalization of the model equations is as described in Section 2.3, with X being equal to K. The value of r, is taken to be 9 msec. The resulting nondimensional Appendix A. Model Equations 204 Dimensional Model Nondimensional Model] Variable Dimension Variable V mV v 72 - W Ca z T msec t Parameter Dimension Value Parameter Value gCa pS 250 7Ca 0.4974 gCa(V,Ca) PS 10 7Ca(V,Ca) 1.9894e2 gK p5 1300 7K 2.5863 9L p5 50 yi 9.9472e VCa mV 100 VGa 1.25 VK mV -80 VL mV -60 vL -0.75 Vm mV -22 Vm 0.275 V mV -9 -0.1125 V mV -22 v8 -0.275 Sm mV 7.5 3m 9.375e2 Sn mV 10 0.125 S3 mV 10 .s3 0.125 Cm 107F/cm2 10 ka msec1 0.08 /3 5.7578 K, 1 f - 0.001 7.2e4 r im 6 r msec 9.09 ). 1.01 Table A.6: Variables and parameters of the dimensional and nondimensional Chay-Cook model. Appendix B Bifurcation Diagrams for the Polynomial Analog Model 6 5- 4- c) 3 - 2- 1- 0 I I -3 -2 -1 0 1 2 3 U Figure B.1: Bifurcation diagram for region A. Parameter values used to generate this diagram are i = 0.2, i = 0.5, a = 0.25. 205 Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 206 6 ...... Figure B.2: Bifurcation diagram for region Bi. Parameter values used to generate this diagram are = 0.93, i = 0.57, a = 0.25. 6 5- 4- •1•’ c 3 - I 2- 1 1- 0 I I I I I -3 -2 -1 0 1 2 3 4 U Figure B.3: Bifurcation diagram for region B2. Parameter values used to generate this diagram are = 0.93, ‘Ii = 0.97, a = 0.25. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 207 4 Figure B.4: Bifurcation diagram for region B3. Parameter values used to generate this diagram are = 0.9, 2 = 1.15, a = 0.25. 6 4 Figure B.5: Bifurcation diagram for region B4. Parameter values used to generate this diagram are = 0.7, €t = 1.6, a = 0.25. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 208 6 5- 3 Figure B.6: Bifurcation diagram for region B5. Parameter values used to generate this diagram are = 0.47, Li = 1.3, a = 0.25. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 209 (a) 15 6 40 2 3 Figure B.7: (a) Bifurcation diagram for region Cl. Parameter values used to generate this diagram are = 1.45, % = 2.56, a = 0.02. (b) Enlarged view of the region of bistability, showing the upper flopf bifurcation and corresponding branch of periodic orbits, and that the large periodic orbits are stable near the corresponding homoclinic bifurcation. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 210 6 Figure B.8: Bifurcation diagram for region C2. Parameter values used to generate this diagram are = 0.6, z% = 1.77, a = 0.25. 6 3 Figure B.9: Bifurcation diagram for region C3. Parameter values used to generate this diagram are = 0.43, i = 1.56, a = 0.25. Appendix B, Bifurcation Diagrams for the Polynomial Analog Model 211 6 5- 4- cj 3 2- 1- 0 I I I I -3 -2 -1 0 1 2 3 U Figure B.10: Bifurcation diagram for region C4. Parameter values used to generate this diagram are i = 0.2, €t = 1.5, a = 0.25. Figure B.11: Bifurcation diagram for region C5. Parameter values used to generate this diagram are i = 0.7, i2 = 2.1, a = 0.25. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 212 2: . Figure B.12: Bifurcation diagram for region C6. Parameter values used to generate this diagram are = 0.3, % = 2.6, a = 0.25. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 212 Figure B.12: Bifurcation diagram for region C6. Parameter values used to generate this diagram are = 0.3, ft = 2.6, a = 0.25. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 213 15 i i i i i i i i i 45 6 Figure B.13: Bifurcation diagram for region Dl. Parameter values used to generate this diagram are = 2.3, ft = 0.25, a = 0.25. 45 1 i i Figure B.14: Bifurcation diagram for region D2. Parameter values used to generate this diagram are = 2.2, ft = 0.46, a = 0.25. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 214 Figure B.15: Bifurcation diagram for region D3. Parameter values used to generate this diagram are = 1.87, % = 0.11, a = 0.25. Figure B.16: Bifurcation diagram for region D4. Parameter values used to generate this diagram are = 1.7, ?2 = 0.1, a = 0.25. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 215 40 56 Figure B.17: Bifurcation diagram for region D5. Parameter values used to generate this diagram are = 2.0, = 0.4, a = 0.25. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 216 (a) 1J i i V 25 6 .. Figure B.18: (a) Bifurcation diagram for region El. Parameter values used to generate this diagram are = 2.0, ii = 1.1, a = 0.02. (b) Enlarged view of the region of bistability, showing that the large periodic orbits are unstable near the homoclinic bifurcation. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 217 /\\ I I ..........,.. Figure B.19: Bifurcation diagram for region E2. Parameter values used to generate this diagram are i = 1.45, i = 0.5, a 0.02. 6 \I 5- 4- c) 3- 4 Figure B.20: Bifurcation diagram for region E3. Parameter values used to generate this diagram are = 1.1, ‘z2 = 0.85, a = 0.02. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 218 Figure B.21: Bifurcation diagram for region E4. Parameter values used to generate this diagram are 1.2, i% = 1.0, a = 0.25. ............................. \. 4 Figure B.22: Bifurcation diagram for region E5. Parameter values used to generate this diagram are = 1.0, = 1.45, a = 0.02. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 219 6 5- 4- 1’ 3 ./••• 2- 1- 0 I I I -3 -2 -1 0 1 2 3 U Figure B.23: Bifurcation diagram for region Fl. Parameter values used to generate this diagram are = 1.4, i = 0.2, a = 0.25. 6 5- . 4 ....... 2- 1 0 I I I -3 -2 -1 0 1 2 3 U Figure B.24: Bifurcation diagram for region F2. Parameter values used to generate this diagram are = 1.51, ii = 0.0, a = 0.02. Appendix B. Bifurcation Diagrams for the Polynomial Analog Model 220 6 5 4 / L) 3 •1 0 I I I I I I -4 —3 -2 -1 0 1 2 3 4 U Figure B.25: Bifurcation diagram for region F3. Parameter values used to generate this diagram are = 1.447, i = 0.3, a = 0.02. 0j Figure B.26: Bifurcation diagram for region F4. Parameter values used to generate this diagram are = 1.64, i% = 0.0, a = 0.02.

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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-0079910/manifest

Comment

Related Items