Georgia State University Digital Archive @ GSU Physics and Astronomy Dissertations Department of Physics and Astronomy 5-7-2011 Mechanisms of Multistability in Neuronal Models Tatiana Malashchenko Georgia State University, tmalashehenko1@student.gsu.edu Recommended Citation Malashchenko, Tatiana, "Mechanisms of Multistability in Neuronal Models" (2011). Physics and Astronomy Dissertations. Paper 46. http://digitalarchive.gsu.edu/phy_astr_diss/46 This Dissertation is brought to you for free and open access by the Department of Physics and Astronomy at Digital Archive @ GSU. It has been accepted for inclusion in Physics and Astronomy Dissertations by an authorized administrator of Digital Archive @ GSU. For more information, please contact digitalarchive@gsu.edu.

MECHANISMS OF MULTISTABILITY IN NEURONAL MODELS by TATIANA MALASHCHENKO Under the direction of Gennady S. Cymbalyuk ABSTRACT Multistability is a fundamental attribute of the dynamics of neuronal systems under normal and pathological conditions. The mechanism of bistability of bursting and silence is not well understood and to our knowledge has not been experimentally recorded in single neurons. We considered four models. Two of them described the dynamics of a leech heart interneuron: the canonical model and a low-dimensional model. The other two models described mammalian pacemakers from the respiratory center. We investigated the low-dimensional model and identified six different types of multistability of dynamical regimes. We described six generic mechanisms underlying the co- existence of oscillatory and silent regimes. The mechanisms are based either on a saddle equilibrium or a saddle periodic orbit. The stable manifold of the saddle equilibrium or the saddle orbit sets the threshold between the regimes. In the two models of the leech interneuron

the range of the controlling parameters supporting the co-existence of bursting and silence is limited by the Andronov-Hopf and homoclinic bifurcations (Malashchenko, Master Thesis 2007). The bistability was found in a narrow range of the leak currents' parameters. Here, we introduced a propensity index to bistability as the width of the range on a bifurcation diagram; we investigated how the propensity index was affected by modifications of the ionic currents, and found that conductances of only two currents substantially affected the index. The increase of the conductance of the hyperpolarization-activated current, Ih, and the reduction of the fast Ca2+ current, ICaF, notably increased the propensity index. These findings define modulatory conditions under which we suggest the bistability of bursting and silence could be experimentally revealed in leech heart interneurons. We hypothesize that this mechanism could be commonly found in a large variety of neuronal models. We applied our techniques to models of vertebrate neurons controlling respiratory rhythm, which represent two types of inspiratory pacemakers of the Pre-Bӧtzinger Complex. We showed that both types of neurons could exhibit bistability of bursting and silence in accordance with the mechanism which we described. INDEX WORDS: Multistability, Bistability, Bursting, Silence, Co-existence, Dynamical systems, Computational neuroscience, Doctor of Philosophy, Georgia State University

MECHANISMS OF MULTISTABILITY IN NEURONAL MODELS by TATIANA MALASHCHENKO Dissertation Submitted in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy in the College of Arts and Sciences Georgia State University 2011

Copyright by Tatiana Igorevna Malashchenko 2011

MECHANISMS OF MULTISTABILITY IN NEURONAL MODELS by TATIANA MALASHCHENKO Committee chair: Gennady Cymbalyuk Committee: Andrey Shilnikov Donald Edwards Vadym Apalkov Mukesh Dhamala Electronic Version Approved: Office of Graduate Studies College of Arts and Sciences Georgia State University May 2011

iv ACKNOWLEDGEMENTS It is my pleasure to thank those who made this dissertation possible. First, I would like to thank my adviser Dr. Gennady S. Cymbalyuk who have guided and supported me from the initial to the final level of my graduate study. I would like to thank members of my scientific committee: Dr. Andrey Shilnikov, Dr. Donald Edwards, Dr. Vadym Apalkov and Dr. Mukesh Dhamala and the chairmen of the Physics and Astronomy department Richard Miller. I owe my deepest gratitude to Dr. Andrey Shilnikov who helped me to learn the bifurcation analysis and for his support and advisement during my work on the master thesis. I would also like to thank Dr. Donald Edwards and Dr. Paul Katz for inspiration in neurobiology classes. I thank William Barnett, Ariel Hart and Lev Dashevskiy for a careful reading of the dissertation and numerous comments. I offer my regards and blessings to all my family members who supported me in any respect during my study and the completion of the dissertation. This work was supported by the Brain and Behavior program of Georgia State University, Atlanta, GA. and NSF grant PHY-0750456.

v TABLE OF CONTENTS ACKNOWLEDGEMENTS iv LIST OF TABLES vii LIST OF FIGURES viii 1 INTRODUCTION 1 1.1 Bistability of Bursting and Silence in the 14D Model 6 1.2.1 The Unstable Orbit Disappearing at a Homoclinic Bifurcation 8 1.2 Multistability Induced by Currents Modulation 15 2 RESULTS 17 2.1 Chapter I: Coexistence of Oscillatory and Silent Regimes in a Simplified Model of a Neuron 17 2.1.1 Model {ICaS INa} 18 2.1.2 Canonical Parameters of {ICaS INa} 11 2.1.3 Equilibria and Oscillatory Regimes: The (gleak Eleak)-Parameter Bifurcation Diagram 23 2.1.4 Switching Between Bursting and Silent Regimes by a Pulse of Current 37 2.1.5 Conclusion I 42 2.2 Chapter II: Similarity Between the Half-Center Oscillator and the 14D Models Dynamics 43 2.2.1 14D Model of Leech Heart Interneuron 44 2.2.2 Switching Between Bursting and Silent Regimes 46 2.2.3 Transient Bursting Activity 51 2.2.4 Conclusion II 54

vi 2.3 Chapter III: Propensity to Bistability of Bursting and Silence of the Leech Heart Interneuron 54 2.3.1 Current Conductances Differently Affect Bistability 55 2.3.2 Changes in 2D Bifurcation Diagram of Stationary and Oscillatory States When Ih and ICaF are Modified 64 2.3.5 Conclusion III 68 2.4 Chapter IV: Modeling of Neurons of Respiratory Center and Revealing Bistability by Current Modulation 68 2.4.1 Model of the Cd2+ - Insensitive Pacemaker 69 2.4.2 Bistability of Bursting and Silence of the Cd2+ - Insensitive Pacemaker with Ih Modulation 73 2.4.3 Model of the Cd2+ - Sensitive Pacemaker 78 2.4.4 Bistability of Bursting and Silence of the Cd2+ - Sensitive Pacemaker 83 Induced by Blockade of Ip 2.3.4 Conclusion IV 86 3 DISCUSSION 86 4 REFERENCES 97 5 APPENDICES 107 5.1 Model of Half Center Oscillator 91107

vii LIST OF TABLES Table 1: Different types of multistability in {ICaS INa} model; separating regimes and mechanisms supporting multistability 42 Table 2: Minimum and maximum values of the propensity index for the different voltage-gated conductances 64 Table 3: Initial conditions of HCO model 109

viii LIST OF FIGURES Figure 1: Switching the regime of activity from silence into bursting by a square pulse of current 8 Figure 2: Bifurcation diagram of the stationary states and hyperpolarized subthreshold oscillations of the model 11 Figure 3 Co-existence of bursting and silence 13 Figure 4 Bifurcation diagram of the oscillatory and stationary regimes 14 Figure 1.1: Evolution of temporal characteristics of the bursting activity with half- inactivation voltage of the fast sodium current, Bh, varied 22 Figure 1.2: Transition from tonic spiking into bursting and evolution of the bursting waveforms in response to the increase of the half-inactivation voltage of the slow calcium current, BhCaS 24 Figure 1.3: Dependence of the burst duration, interburst interval and spike frequency of the bursting activity on BhCaS 25 Figure 1.4: Bursting activity of the canonical 4D model 26 Figure 1.5: The (gleak Eleak)-parameter bifurcation diagram for the oscillatory and stationary regimes of the model 27 Figure 1.6: Inset from the (gleak Eleak)-parameter bifurcation diagram for the oscillatory and stationary regimes of the model 29 Figure 1.7: The (gleak V)-bifurcation diagrams for two values of the leak reversal potential 30

ix Figure 1.8: Bifurcation diagram showing the evolution of equilibria and oscillatory regimes for two values of the leak reversal potential 32 Figure 1.9: Voltage traces of the co-existence of bursting and subthreshold oscillations and the co-existence of bursting and tonic spiking 36 Figure 1.10: Stable subthreshold oscillations coexist with bursting and silent attractors 28 Figure 1.11: Examples of switching between bursting and silence produced by square pulses of injected current 40 Figure 1.12: Revealing of the unstable oscillations 41 Figure 2.1: Unstable sub-threshold oscillations recorded after the perturbation by a pulse of current 49 Figure 2.2: The strength-duration relationship of the minimum depolarizing and hyperpolarizing square pulses of current switching the silent regime into bursting regime 50 Figure 2.3: The switch of bursting into silence in the half-center oscillator (HCO) by 51 a pulse of current delivered to one neuron from the pair Figure 2.4: Long-lasting transient transitions from bursting to silence 53 Figure 3.1: Effect of the variation of the maximum conductances of IP, INa, ICaF, and ICaS currents on the propensity index to the bi-stability of bursting activity and silence 59 Figure 3.2: Effect of the variation of the maximum conductances of IK1, IK2, IKA, and Ih currents on propensity index for bi-stability of bursting and silence 60 Figure 3.3: Dependence of the propensity index on gP, gNa, gCaF and gCaS 61 Figure 3.4: Dependence of the propensity index on gK1, gK2, gKA and gh 63

x Figure 3.5: Bifurcation diagrams of the oscillatory and rest regimes in parameter space (gleak, Eleak) for different gh 66 Figure 3.6: Bifurcation diagrams of the oscillatory and rest regimes in parameter space (gleak, Eleak) as gh and gCaF are varied 67 Figure 4.1: The steady-state curves of the activation gating variable and steady-state conductance of Ih of the Cd2+ Insensitive pacemaker 72 Figure 4.2: Bursting activity exhibited by the model of Cd2+ -insensitive pacemaker 72 Figure 4.3: The gleak-parameter bifurcation diagram of the oscillatory and stationary states of the Cd2+ - insensitive pacemaker 75 Figure 4.4: The gleak-parameter bifurcation diagram of the oscillatory and stationary states of the Cd2+ - insensitive pacemaker when the activation of Ih is changed 76 Figure 4.5: Switches between bursting activity and silence in the model of the Cd2+ - insensitive pacemaker 77 Figure 4.6: Bursting activity of the model of the Cd2+ - sensitive pacemaker 80 Figure 4.7: The gleak-parameter bifurcation diagram of oscillatory and stationary states of the Cd2+ - sensitive pacemaker 82 Figure 4.8: The gleak-parameter bifurcation diagram of oscillatory and stationary states of the Cd2+ - sensitive pacemaker, with Ip blocked 84 Figure 4.9: Bistability of bursting activity and silence in the model of the Cd2+ - sensitive pacemaker with Ip blocked 85

1 1. INTRODUCTION Multistability of oscillatory and silent regimes is a ubiquitous phenomenon exhibited by excitable systems such as neurons and cardiac cells (Fuentealba et al., 2005; Egorov et al., 2002; Hounsgaard et al., 1989; Loewenstein et al., 2005; Jalife and Antzelevitch, 1979). As a feature it appears particularly advantageous for neurons that are part of multifunctional central pattern generators (Marder, 1994; Getting, 1989; Briggman and Kristan, 2008; Venugopal et al., 2007; Shilnikov et al., 2008). In motor control, coexistence of silent regime and tonic spiking activity is found in spinal motoneurons of the cat, rat, turtle and frog (Hounsgaard et al., 1984; Hounsgaard and Kiehn, 1989; Conway et al., 1988; Eken and Kiehn, 1989; Kiehn and Eken, 1998; Perrier and Hounsgaard, 2000). This multistability is involved in posture maintenance (Eken and Kiehn, 1989; Kiehn and Eken, 1998). Some mechanisms underlying bistability- for example, the coexistence of tonic spiking and silence and the coexistence of tonic spiking and bursting - have been intensively studied and are well understood (Rinzel, 1978; Guttman et al., 1980; Gutkin et al., 2009; Fröhlich et al., 2006; Shilnikov et al., 2005; Channell et al., 2007; Shilnikov and Cymbalyuk, 2004). Surprisingly, there is a gap in our knowledge of dynamical mechanisms supporting the bistability of bursting and hyperpolarized silence, bursting and subthreshold oscillations, and multistability of bursting, subthreshold oscillations and silence. The classification of mechanisms supporting the multistability of oscillatory and silent regimes is yet incomplete, and remains a fundamental problem for both neuroscience and the theory of dynamical systems. Bursting activity is the result of an interplay of ionic currents which are voltage-gated on various time scales. We envisage a neuron as a slow-fast dynamical system. The qualitative theory of dynamical systems provides a rigorous description of scenarios producing the co-

2 existence of different attracting regimes, i. e., multistability, in the system’s dynamics. Exemplary studies by Rinzel (1978) and by Guttman, Lewis and Rinzel (1980) formulated and answered a set of questions which describe a basic scenario of bistability of tonic spiking and silence. It is based on the presence of a saddle periodic orbit separating the basin of attraction of the tonic spiking periodic orbit from the equilibrium representing the silent regime. The scenario also describes the modulation of the neuron’s dynamics in response to variations of a bifurcation parameter. According to this scenario, the saddle periodic orbit emerges through a subcritical Andronov-Hopf bifurcation and disappears through a saddle-node bifurcation for periodic orbits; these two bifurcations define the boundaries for the bistability region. This scenario describes hysteresis and catastrophe-like, fast and nonreversible transitions between silence and tonic spiking as the bifurcation parameter is varied. Guttman, Lewis and Rinzel (1980) showed that a switch between these regimes can be executed by a pulse of current in experiments on the squid giant axon in saline with low calcium concentration. This study is precipitated by our keen interest in the dynamics of the leech heart interneuron, which constitutes the core of the leech heartbeat timing network. They are found as pairs of mutually inhibitory neurons located in ganglia 3 and 4 (Calabrese et al., 1995). This preparation provides a unique opportunity for studying cellular and network mechanisms of bursting. A leech heart interneuron can be decoupled from the network with bicuculline (Schmidt and Calabrese, 1992); its endogenous and network activities are well described by the canonical model (Hill et al., 2000; Cymbalyuk et al., 2002). This model has been instrumental in predicting that these neurons have endogenous dynamics supporting bursting activity in a single cell, and has showed the high sensitivity of the bursting regime to variations of the leak current. This sensitivity explained why these interneurons show bursting activity while recorded

3 extracellularly, and tonic spiking while recorded intracellularly. Biophysically, the leak current is mediated by multiple channels. Channels with mixed ionic permeability to Na+, Cl-, K+ and Ca2+ ions, channels pervious to K+ such as TASK channels, or channels with primarily Na+ permeability such as NALCN are examples of leak channels that control neuronal resting membrane potential (Bayliss et al., 2003; Lu et al., 2007; Koizumi et al., 2010). Depending on the membrane composition of these channels, neuromodulations targeting this type of channel can reduce or augment leak conductance. For example, TASK leak channels are pH-sensitive and an increase of the extracellular pH up-regulates leak conductance (Larkman and Perkins, 2005). Neurotransmitters including serotonin and noradrenalin close TASK channels leading to the decrease of total leak conductance (Sirois et al., 2002; Perrier et al., 2003). The model of the leech heart interneuron is a system of 14 ordinary differential equations with variables operating on different time scales. In a previous study by Cymbalyuk et al. (2000), a two-parameter bifurcation diagram has been numerically obtained, mapping oscillatory and stationary regimes on the plane of the leak current’s parameters, which are the conductance and the reversal potential of the leak current (gleak, Eleak). It shows the borders describing transitions from silence to bursting and from bursting to silence, bounding the zone where bursting and silence coexist (Cymbalyuk et al., 2002). A rigorous analysis of the canonical model of a leech heart interneuron is a difficult subject. Therefore we developed a low dimensional model and carried out the analysis of this model (Malashchenko, Master Thesis, 2007; Chapter I). In the master thesis we have described two types of multistability; here, a detailed bifurcation analysis allows us to describe four more types of multistability and dynamical mechanisms underlying them (Chapter I).

4 Our approach toward this aim is to simplify the canonical model so that the core mechanisms could be thoroughly studied while similar changes of the dynamics still could be implemented experimentally. Thus, in addition to a commonly used slow-fast dissection, we designed the model so that it would represent a pharmacologically achievable reduction through a removal of certain ionic currents. This approach maintains possibilities for experimental validation of the obtained results. It proved fruitful in previous theoretical studies by Cymbalyuk et al. of a minimal model of the leech heart interneuron under a blockade of Ca2+ currents along with a partial block of outward currents (Shilnikov et al., 2005; Cymbalyuk et al., 2001; Cymbalyuk and Calabrese, 2001; Shilnikov and Cymbalyuk, 2008). This simplified model, describing the dynamics of the transient sodium current, INa, and non-inactivating slow potassium current, IK2, is given by a system of three differential equations with slow and fast time scales. The inactivation of INa and membrane potential constitute the fast subsystem; and the activation of IK2 is the slow variable. The reduced model introduced here is based on INa; the slow, low-threshold calcium current, ICaS; and the leak current, Ileak. We have the following rationale for this choice of currents. ICaS provides the slowest variable of the 14D canonical model of the leech heart interneuron. It has been suggested for the 14D model that ICaS underlies the bursting activity (Cymbalyuk et al., 2002). Also, it was demonstrated that ICaS introduced through dynamic clamp can reinstate the bursting activity of a tonically spiking leech heart interneuron having all Ca2+- currents blocked and being isolated from other neurons by the application of saline containing Mn2+ (Olypher et al., 2006). The significance of this achievement is apparent in light of past experiments showing that these neurons are sensitive to the leak current’s parameters, so that bursting activity of a neuron pharmacologically singled out from the network could not be

5 intracellularly recorded. Olypher et al. showed that an intrinsic mechanism for regulating burst duration might be based on the kinetics of ICaS inactivation, and we further corroborated that assertion in this study (Chapter I). This model makes a complementary study to previous works of Cymbalyuk and Calabrese (2001) on simplified models, since it employs the same fast subsystem, but differs in terms of its slow variable (Shilnikov et al., 2005; Cymbalyuk and Shilnikov, 2005; Shilnikov and Cymbalyuk, 2005; Shilnikov et al., 2008; Channell et al., 2007). It allows us to focus the investigation on the potential roles of ICaS in the dynamics of a leech heart interneuron. We have previously constructed a simplified model described by the set of ionic currents and their kinetics ICaS, INa (Malashchenko, Master Thesis, 2007). Here, we further tune the model; we sweep the parameters determining the kinetics of the currents and choose one set which produced activity with temporal characteristics close to those recorded experimentally from leech heart interneurons (Chapter I). In Chapter I we describe the mechanisms supporting multistability in the simplified model under variation of parameters of the leak current. We show that the simplified model can exhibit six different types of multistability. We focus on five types, involving the co-existence of 1) tonic spiking and silence, 2) tonic spiking and subthreshold oscillations, 3) bursting and subthreshold oscillations, 4) the co-existence of bursting, subthreshold oscillations and silence and 5) the-coexistence of bursting and silence. Two types of bistability, the bistability of bursting and silence and the bistability of tonic spiking and bursting, have been shown for the low-dimensional model in our previous work (Malashchenko, Master Thesis, 2007). We show that all these five different types of multistability occur due to a presence of a separating regime that is either a saddle periodic orbit or a saddle equilibrium. We found that the parameter range wherein multistability is observed is limited by the parameter values, at which the separating regimes emerge and terminate (Chapter

6 I). The knowledge of the dynamical mechanism allowed us to investigate the ranges of parameters supporting multistability. In the following section 1.2 we applied this approach to the canonical leech heart interneuron (Malashchenko, Master Thesis, 2007). 1.1 Bistability of Bursting and Silence in the 14D Model We studied the dynamics of the canonical leech heart interneuron model. Bistability of bursting and silent regimes had been demonstrated previously (Cymbalyuk et al., 2002), but the mechanism supporting it was not determined. The model contained nine voltage dependent currents and is described by a system of 14 differential equations (Hill et al., 2001). The parameters of the 14D model were tuned to produce activity with a waveform close to the one experimentally observed (Hill et al., 2001; Cymbalyuk et al., 2002). The trajectories of this model were obtained using the Matlab ODE solver, ode15s. Absolute tolerance and relative tolerance were 10−9 and 10−8. The integration and bifurcation analysis were performed using the software package CONTENT which is freely available at http://www.staff.science.uu.nl/~kouzn101/CONTENT/ (Khibnik et al., 1993). The integration of equations was done using the Runge-Kutta method of the 4-th order with the tolerance of integration set as 10−9. Bistability was reported in the model with elevated conductance of the leak current, g leak . We were able to elicit a switch from silence to bursting activity by applying a square pulse of current. The switch could be triggered by either a negative or a positive pulse (Figure 1). The duration of the pulse was 0.03 sec. All parameters including Eleak were fixed to the canonical

7 values, but g leak was set to 10.7 nS; the canonical value is 9.9 nS (Hill et al., 2001; Cymbalyuk et al., 2002). The initial conditions of voltage and gating variables were set so that the model initially exhibited silence. In Figure 1 A a hyperpolarizing pulse of current with the amplitude of -0.05 nA switched the activity from silence to a bursting regime. By trying different amplitudes of the pulse, we found that negative pulses within the range -0.0213 nA <Iinj < 0.0 nA did not switch silence into the bursting regime, while negative pulses with an amplitude larger than 0.213 nA switched the activity into the bursting regime. Perturbations of the initially silent neuron by the depolarizing square pulse of current showed that only a pulse with an amplitude larger than 0.0175 nA could switch the activity from silence to bursting. Figure 1.1B shows the example of perturbation by a positive pulse with the amplitude Iinj = 0.05 nA. The critical values, Iinj = -0.0213 nA and Iinj = 0.0175 nA, set two thresholds, negative and positive, for the amplitude of the pulse of current for the switch from silence to bursting. The previous study of the leech heart interneuron model showed that the bistability of bursting and silence is associated with the Andronov-Hopf bifurcation (Cymbalyuk et al., 2002). These results suggest that the separating barrier between the two attractors, bursting and silent is the stable manifold of a saddle periodic orbit.

8 Figure 1: Switching the regime of activity from silence into bursting by a square pulse of current. A) A negative, hyperpolarizing pulse of current Iinj= -0.05 nA switches activity from silence to bursting. The minimum magnitude of Iinj which can induce the switch is -0.02131 nA. B) A positive, depolarizing pulse of current Iinj= 0.05 nA also initiates the bursting activity. To switch the activity, Iinj has to be larger than the value of 0.0175 nA. The parameters of the leak current are Eleak = -0.0635 V and g leak = 10.7 nS. 1.2.1 The Unstable Orbit Disappearing at a Homoclinic Bifurcation We analyzed the stability of the equilibria of the model for E leak = -0.0635 V. Evolution of the equilibria and their stability were analyzed with gleak being varied as the bifurcation parameter in the CONTENT software. The plot of the membrane potential associated with each equilibrium versus the bifurcation parameter gleak exhibited a Z-shaped curve of the equilibria branches. This curve has two branches representing the equilibria, which are the hyperpolarized and depolarized rest states, separated by the saddle. First, we set a large value of gleak, so that the model stayed at the stable hyperpolarized stationary state. We examined how the stability of the hyperpolarized equilibrium changed in response to a decrease in the bifurcation parameter gleak (Figure 2). The hyperpolarized equilibrium lost its stability through the sub-critical Andronov-

9 AH Hopf bifurcation. For Eleak = -0.0635 V, the bifurcation occurred at g leak = 10.67 nS; and the stable rest state became a saddle-focus. The sub-critical Andronov-Hopf bifurcation (AH) gave rise to an unstable periodic orbit with zero amplitude and period of 3.05 sec. As the bifurcation parameter gleak was increased, the amplitude of the unstable orbit grew proportionally to g leak 10.67 while gleak was close to the bifurcation value g leak . As the parameter g leak was AH increased further, the unstable orbit approached the saddle located on the middle branch of the Z- shaped curve. To depict this evolution of the orbit, we plotted the minimum, maximum and the average values of the membrane potential of the unstable orbit against g leak together with the Z- shaped curve. In the vicinity of the saddle, the period of the orbit grew logarithmically fast. For the range of g leak from 10.85 nS to 10.873 nS (Figure 2 B) the period of the orbit grew rapidly from 7 sec to 100 sec. The graph was fitted by the function -1.198∙ln(|10.873 – g leak |). Near g leak = 10.873 nS the amplitude of the oscillations stayed constant while the period grew. At the value g leak = 10.873 nS the homoclinic bifurcation (Hom) occurred in the system. For g leak > 10.873 nS, the stable rest state was the only attractor. Concerning the unstable subthreshold oscillations this analysis suggests that they could be recorded experimentally. In the model, we could compare the unstable periodic orbit obtained using the bifurcation analysis software CONTENT with that recorded after a square pulse of current is applied with an amplitude very close to the threshold value. Let’s consider again the canonical model with Eleak = -0.0635 V and g leak = 10.7 nS. We set the initial conditions, which corresponded to the bursting regime. With these initial conditions the model would exhibit bursting indefinitely; to confirm that the bursting regime is the attractor, the system was integrated for several thousand seconds. To perform the switch from bursting to silence, we

10 applied a negative pulse of current near the first spike of a burst. Pulses with an amplitude of the current pulse close to -0.010445 nA revealed unstable subthreshold oscillations for several periods (Figure 3 B,D). Now we plotted together the saddle periodic orbit obtained with the parameter continuation software CONTENT for g leak =10.7 nS and the subthreshold oscillations recorded after the pulse (Figure 3 E,F). The unstable periodic orbit is presented in red. The trajectory of the unstable subthreshold oscillations is extracted from the trace in Figure 3 B between two vertical lines and is shown in blue in the inset Figure 3 E,F. Minimum and maximum voltages of the unstable periodic orbit are also plotted in Figure 3 B and Figure 3 D. The two trajectories appeared very close to each other when projected on the plane created by the state variables gating the conductance of the slow calcium current (Figure 3 E) and also plotted as membrane potentials versus time (Figure 3 F). Figure 3 illustrates the position of the unstable periodic orbit relative to the rest state and the bursting trajectory projected onto the plane (mCaS, hCaS). Figure 3 allows one to observe how a pulse of current brings the phase point of the bursting trajectory into the basin of attraction to the rest state, which is separated by the stable manifold of the saddle periodic orbit. To summarize, the unstable periodic orbit appears through the sub-critical Andronov-Hopf bifurcation and disappears at the homoclinic bifurcation; and it corresponds to unstable subthreshold oscillations, which could be revealed by a pulse of current.

11 Figure 2: Bifurcation diagram of the stationary states and hyperpolarized subthreshold oscillations of the model. A) The solid green line represents the stable hyperpolarized equilibrium (EP1). Subcritical Andronov-Hopf bifurcation (AH) occurs at g leak =10.67 nS (marked by the left vertical dashed line, Eleak=- 0.0635 mV) where the equilibrium loses its stability and gives rise to the unstable periodic orbit. For the values of g leak smaller than this critical value the equilibrium is unstable (dashed blue line). The dashed blue line located above it depicts the saddle equilibria (EP2). The curve of the depolarized equilibrium has a membrane potential near 0.01 V and is not shown. The unstable periodic orbit corresponds to the unstable subthreshold oscillations (USTO). It is marked by the two dashed red curves and the dashed brown curve locating the minimum, maximum, and average values of the membrane potential of the oscillations correspondingly. The unstable orbit disappears at the homoclinic bifurcation

12 (Hom), marked by the dashed vertical line at the right, g leak = 10.873 nS. B) Red dots graph the values of the period of USTO as the g leak is varied between the two bifurcation values. As the unstable periodic orbit approaches the homoclinic bifurcation, the period grows as -1.198∙ln(|10.873- g leak )| as shown by the blue curve. C) Inset shows the part of the graph in B) taken in the vicinity of the homoclinic Hom bifurcation, where b=1.198, g leak =10.873 nS. We numerically computed a two-parameter bifurcation diagram of oscillatory and stationary regimes (Figure 4). The range of the parameters of the leak current (gleak Eleak) supporting bursting activity has a banana shape and is surrounded by a large area supporting tonic spiking activity on the left side and a silent regime on the right side. The borders between regimes are associated with bifurcations in the system. Here, we explored the transition from bursting into silence and the mechanism supporting bistability of these two regimes. Based on the knowledge that the range of parameters where unstable subthreshold oscillations exist is bounded between the Andronov-Hopf and homoclinic bifurcations, first, we computed the AH curve in (gleak, Eleak) parameter space. Then, in order to locate the homoclinic bifurcation curve (Hom), we increased g leak while following the saddle periodic orbit. We started at the Andronov-Hopf bifurcation where the orbit appeared and followed the orbit until its period became very large, 100 sec, near the homoclinic bifurcation. CONTENT software allows tracking the periodic orbit with the given value of the period; we used g leak and Eleak as the two parameters to find the curve with a fixed period of 100 sec. This curve gives a good estimation for the bifurcation values determining the homoclinic bifurcation. The estimated bifurcation values are shown by the red circles in Figure 4. One can see that the detected points representing

13 the homoclinic bifurcation of the unstable periodic orbit match well with the border between bursting and silence calculated previously by numerical simulations (Cymbalyuk et al., 2002). Figure 3: Co-existence of bursting and silence. The negative square pulse of Iinj perturbed the bursting activity. Bursting is shown in dark blue. A,C,E) 2D projection reveals unstable oscillation (red circle) along with the two attractors: the bursting and the hyperpolarized equilibrium (green dot). A,B) A pulse with the amplitude of -0.01044 nA moved the phase point towards the unstable periodic orbit (in red), and allowed recording the subthreshold oscillations for eight periods, but did not switch the activity from bursting into silence. C,D) The pulse with the amplitude -0.01045 nA again allowed recording the subthreshold oscillations for a few periods, and moved the phase point into the basin of attraction of the equilibrium, thus producing the switch from bursting into silence. The red dashed lines in plots B) and D)

14 mark the minimum and maximum values of the membrane potential of the unstable orbit obtained with the bifurcation analysis. E,F) This is the part of the trajectory with the unstable oscillations revealed by the perturbation (in blue) plotted together with the unstable periodic orbit (in red). The parameters of the leak current were the same as in Figure 1. Figure 4: Bifurcation diagram of the oscillatory and stationary regimes. The green curve is composed of numerically found points where the transition from bursting activity into silence occurs. The sub-critical Andronov-Hopf bifurcation (AH) of the hyperpolarized rest state (silent regime) is shown by the blue curve and marks the boundary where the silent regime loses stability, giving rise to unstable subthreshold oscillations. The red circles correspond to the homoclinic bifurcation (Hom) of the unstable subthreshold oscillations. The area between these two curves (the blue and the red dots) marks the parameter regime where unstable subthreshold oscillations exist. The grey curve marks the period-doubling bifurcation (PD) of large amplitude periodic spiking and marks the transition from tonic spiking to bursting activity. The area of bursting activity is located between the PD and Hom curves.

15 Hence, the transition from bursting activity into silence was associated with the homoclinic bifurcation, and the range of bistability in leak current parameter space was limited by the sub-critical Andronov-Hopf and homoclinic bifurcations for the saddle periodic orbit. The properties of perturbation switching between bursting and silence in a single neuron we investigate in Chapter II. We also showed that similar switches could be induced in and in a half- center oscillator (Chapter II). 1.2 Multistability Induced by Currents Modulation The ability of networks to function in more than one regime can be induced through neuromodulation by affecting the cellular membrane dynamics. For example, in motoneurons bistability results from cellular membrane dynamics, and synaptic input plays the role of a toggle-switch between two regimes (Kiehn and Eken, 1998). Such cellular bistability is an endogenous property of motoneurons that can be recruited when necessary by neuromodulators targeting the dynamics of intrinsic ionic currents (Hounsgaard and Kiehn, 1985; Conway et al., 1988). Numerous studies concentrate on studying the ionic basis for inward and outward currents supporting multistability (Hounsgaard and Kiehn, 1989; Hsiao et al., 1997; Hughes et al., 1999; Williams at el., 2002). The questions we address in Chapter III are: Can different currents be evaluated with respect to the propensity of a neuron to multistability? Which ionic current affects multistability, and how could it be regulated to induce neuronal multistability? The mechanisms underlying multistability specify the unstable regimes that form barriers between attractors. Defining the unstable regimes and the bifurcations bounding the existence of these regimes in the parameter space are the key parts of the description of the mechanisms supporting multistability in neuronal models. For example, in cardiac pacemakers and the squid

16 giant axon such an unstable regime is a saddle orbit, which appears through Andronov-Hopf bifurcation. A saddle orbit represents unstable subthreshold oscillations (Rinzel, 1978; Guttman et al., 1980; Landau et al., 1990). A saddle orbit has also been detected in the mechanism supporting the bistability of tonic spiking and bursting in a neuronal model (Shilnikov and Cymbalyuk, 2005; Shilnikov et al., 2005). Biological systems which were initially exhibiting a single regime, could turn bistable after the modulation of membrane ionic currents. For example, for the narrow range of the leak current parameters, the model of the leech heart interneuron exhibits bistability of bursting activity and silence (Cymbalyuk et al., 2002). We suggest that the interval of this range can be used as an index evaluating a neuronal propensity for bistability. Since the bistability is observed within a small range of the leak conductance, we look for the conditions that expand this range. Here we use the knowledge about the mechanism supporting the bistability of bursting and silence described in the Master Thesis by Malashchenko (2007). In Chapter III, we investigate how every voltage-gated ionic conductance affects the propensity of a neuron to bistability and how it changes the saddle periodic orbit with leak current parameters varied. We suggest that this is a generic approach for the analysis of the dynamics of bistable neurons that highlights the role of identified currents. We found that only two currents substantially affects the index of the propensity to bistability. These data suggests that under conditions reported here the coexistence of bursting and silence in an isolated leech heart interneuron could be revealed experimentally (Chapter III). We hypothesize that the mechanism of bistability of bursting and silence described in this work would be commonly found in the parameter space near the transition from bursting to silence in a large variety of neuronal models. In Chapter IV, we investigate two models of

17 vertebrate neurons controlling respiratory rhythm, which represent two types of inspiratory pacemakers of the Pre-Bötzinger Complex. They are described as Cd2+ sensitive and Cd2+ insensitive endogenous pacemakers. As a basis for these models we use the model of a respiratory neuron developed by Rybak et al. (2003). We incorporate Ih current into both models, and two types of Ca2+ currents into the model of the Cd2+ Sensitive pacemaker, so that the activity of the models is similar to the experimentally recorded ones (Peña et al., 2004). We apply the techniques developed in the previous chapters and show that both types of neurons could exhibit bistability of bursting and silence in accordance with the mechanism which we described. 2. RESULTS 2.1 Chapter I: Coexistence of Oscillatory and Silent Regimes in a Simplified Model of a Neuron Bursting is a primary oscillatory activity of neurons associated with the control of rhythmic movements. We have developed a low-dimensional model typifying the dynamics of a single leech heart interneuron (Malashchenko, Master Thesis, 2007). The model describes the dynamics of the interneuron under pharmacological conditions, which limit membrane ionic currents to the transient Na+ and the slow inactivating, low threshold Ca2+ currents. We carry out a bifurcation analysis of the model and show that it is susceptible to six different types of multistability of dynamical regimes. Two types of multistability described here have been reported for this model under a different parameter regime (Malashchenko, Master Thesis, 2007). All data presented here are newly obtained for this model, being better tuned to

18 experimentally recorded characteristics of the bursting activity. Also, four more types of multistability are found and described. 2.1.1 Model {ICaS INa} We present a simplified leech neuron model containing only the fast sodium (INa) and slow calcium (ICaS) voltage-dependent currents and the leak current (Ileak). We refer to it as model {ICaS INa}, according to the set of voltage-gated ionic currents which it contains. This model is described by the system of the following four equations: I Na I CaS I leak dV C [ g Na f (150,0.028, V )hNa [V E Na ] g CaS mCaS hCaS [V E CaS ] g leak [V E leak ] I inj ], 3 2 dt dhNa f (500, Bh, V ) hNa , dt 0.0405 dmCaS f (420,0.0472, V ) mCaS , dt m CaS dhCaS f (360, BhCaS , V ) hCaS , dt h CaS The conductance and reversal potential of ICaS were fixed to the values of gCaS=80 nS and ECaS=0.135 V, correspondingly. The conductance and reversal potential of INa were chosen as gNa=250 nS and ENa=0.045 V, correspondingly. Parameters describing the leak current (gleak Eleak) were varied and used as the bifurcation parameters in system analysis. Parameter C is the capacitance of the neuron and its value was set to 0.5 nS. Iinj is the injected current which takes on a certain value when a pulse is applied. The pulse duration is set to 0.03 sec. The pulse amplitude is varied during numerical experiments. The other three equations describe the dynamics of a corresponding gating variable of activation and inactivation, where f (A,B,V ) is a steady state activation (inactivation) function of an ionic current given by f (A,B,V ) = 1/[1 +

19 eA(V +B)]. Here B is the half-activation (half-inactivation) membrane potential at which f =1/2. In the model, the activation of INa is considered to be instantaneous, so that mNa = m∞Na ≡ f (−150, 0.028, V). The voltage-dependent time rates for the activation and inactivation variables of the calcium current are taken from Hill et al. 2001: 0.134 m CaS 0.005 - 400(V+ 0.0487) , 1+ e 5.25 h CaS 0.2 - 250(V+ 0.043) , 1+ e One can see that the inactivation of the calcium current, hCaS, is the slowest variable in the model. Bifurcation analysis was performed using the parameter continuation software CONTENT (Kuznetsov et al., 1996) freely available at http://www.staff.science.uu.nl/~kouzn101/CONTENT/. The solutions of this model were obtained using the Runge-Kutta method of the 4th order and a variable-order method based on numerical differentiation formulas, implemented as the ode15s solver in Matlab (MathWorks, Inc.). Absolute and relative tolerances were set to 10−8 and 10−9, respectively. Because the model is based on a subset of inward currents and hence lacks all outward currents except for the leak current, the balance of inward and outward currents has to be restored. It is achieved in the sense that we tune up the available parameters so that the model produces activity with the temporal characteristics of bursting measured experimentally (Cymbalyuk et al., 2002); the instance of such a model found here we call the canonical model. We take into account basic temporal characteristics of bursting waveforms such as the spike frequency, the burst duration, the interburst interval, the period, the number of spikes and the duty cycle. By spike frequency we

20 mean the average frequency of spikes in a burst. The burst duration is measured as the time between the first spike and the last spike in a burst. The interburst interval is the time between the last spike of a burst and the first spike of the following burst. The period of bursting is measured as the time between the first spike of a burst and the first spike of the consecutive burst, which is the sum of the burst duration and the interburst interval. The duty cycle is the fraction of the period occupied by the burst duration, i.e. the ratio of the burst duration to the period. To explore the waveforms of the bursting we swept values of the half-inactivation voltages Bh and BhCaS of the two voltage-gated currents. Then, the parameters determining the leak current, gleak and Eleak, were used as the bifurcation parameters following Cymbalyuk et al., 2002. 2.1.2 Canonical Parameters of {INa ICaS} A considerable concern about the usability of this model was whether we could tune it up to produce bursting with experimentally observable characteristics. The targeted characteristics were experimentally measured (Cymbalyuk et al., 2002): the burst duration was to be between 3.2 and 6.0 sec, the interburst interval was within 1.5 and 3.0 sec, the duty cycle was between 65.8 and 75.5%. The spike frequency ([6.6 14.7] Hz) was not expected to be achievable since the fast potassium currents affecting the shape of spikes, and thus an appropriate spike frequency, were eliminated from the model. The model parameters satisfying these conditions were easily attained in three steps. First, we varied the maximum conductances, and the inactivation kinetics for INa and ICaS. Relative to the 14D model, we set larger values for gNa = 250 nS, and gCaS = 80 nS to make the corresponding currents more accentuated. Tuning the kinetic parameters Bh and BhCaS was motivated by the notion of the window mode of a voltage-gated ionic current (Hughes et al., 1999; Cymbalyuk and Calabrese, 2001). In this mode an inactivating current

21 exhibits the properties of a persistent, non-inactivating current in the interval of membrane potentials where the steady-state activation and inactivation curves overlap. In the 4D model, the “window” mode of INa could play a role similar to that of the sodium persistent current, which supports burst duration in the canonical 14D model (Hill et al., 2001; Cymbalyuk et al., 2002). A similar role would be played by the “window” mode of ICaS. From the studies of the 14D model we can infer that in the 4D model the activation of ICaS would be responsible for the inception of a burst, while the inactivation of ICaS would control the burst termination. This mechanism of burst termination is similar to the one shown for the dynamics of a half-center oscillator assembled from two heart interneurons (Hill et al., 2002; Olypher et al., 2006). Second, we examined the activity of the model in response to variations of the half- inactivation voltage of INa, Bh, in the function h∞Na≡ f (500,Bh,V) = 1/[1+e500(V +Bh)]. As Bh is increased, the curve h∞Na shifts towards the more hyperpolarized values, closing the window voltage interval. This manipulation of Bh directly affects the steady-state conductance g∞Na, which produces the window current supporting the burst phase of the bursting activity. The overall decrease of g∞Na decreases the period, the burst duration and the duty cycle of the model. The bursting activity occurs within the range of Bh [0.02888 0.03692] V (Figure 1.1). For some range of values of Bh smaller than 0.02888 V, the model shows the tonic spiking activity. In the range of Bh between 0.03692 V and 0.03790 V, the model exhibits stable subthreshold oscillations; and it is silent for Bh > 0.03790 V. The parameters, gleak=15.7 nS, Eleak=-0.0505 V and BhCaS=0.06 V, fixed for this sweep, were chosen so that the model exhibits coexistence of bursting and silence. This sweep of Bh allowed us to adjust the burst duration so that it falls within the experimentally measured range, marked by the grey rectangle in Figure 1.1 A, while the interburst interval and spike frequency were not attained. The discontinuity of the graphs of

22 the interburst interval (Figure 1.1 B) and the spike frequency (Figure 1.1 C) near Bh=0.36 V is due to the integer-number difference in the number of spikes. Figure 1.1: Evolution of temporal characteristics of the bursting activity with half-inactivation voltage of the fast sodium current, Bh, varied. (A) Burst duration, (B) interburst interval, and (C) spike frequency. This sweep of Bh allowed us to adjust the burst duration so that it falls within the experimentally measured range, marked by the grey rectangle in A. The leak current parameters are gleak=15.7 nS and Eleak=-0.0505 V.

23 Third, to adjust the interburst interval, we swept the half-inactivation voltage of ICaS, BhCaS. The increase of BhCaS prolongs the interburst interval of the bursting activity. The interburst interval grows monotonically from 0.53 sec to 3.77 sec as BhCaS is swept from 0.048 V to 0.06 V (Figures 1.2 and 1.3). The graph also shows the effect of variation of BhCaS on the spike frequency within a burst (Figure 1.3). Our parameter sweeps of the model were concluded with the following parameters Bh = 0.031 V, BhCaS = 0.06 V. This adjusted model exhibits bursting activity with a burst duration of 4.5 sec, an interburst interval of 3.8 sec, the period 8.3 sec, and the duty cycle around 54.6%; the number of spikes per burst is 26, and the spike frequency is 5.59 Hz (Figure 1.2D). The leak current parameters are gleak = 15.7 nS and Eleak = -0.0505 V. We further tuned the model by setting gleak = 15.2 nS, so that the temporal characteristics of bursting activity fit well to the experimental data (Figure 1.4). 2.1.3 Equilibria and Oscillatory Regimes: The (gleak Eleak)-Parameter Bifurcation Diagram Parameters of the leak current are remarkable targets for modulation of neuronal excitability. At small values of the leak conductance the neuron is overexcited, and stays silent at the depolarized equilibrium. On the other hand, for large values of the conductance the neuron stays silent at the hyperpolarized equilibrium. For different intermediate values of the conductance, it exhibits a plethora of oscillatory regimes such as various bursting, tonic spiking, and subthreshold oscillations. These regimes represent different levels of excitability. To inventory this rich variety of neuronal dynamics, we created a two-parameter (gleak Eleak) bifurcation diagram for equilibria and oscillatory regimes, which is shown in Figure 1.5.

24 Figure 1.2: Transition from tonic spiking into bursting and evolution of the bursting waveforms in response to the increase of the half-inactivation voltage of the slow calcium current, BhCaS. Such a change of BhCaS shifts h∞CaS towards more hyperpolarized values of VM, changing the activity from tonic spiking (A) to bursting (B-D). (A) For BhCaS=0.047V the model exhibits a periodic tonic spiking activity. (B-D) The increase of BhCaS up to 0.048V shifts h∞CaS towards the hyperpolarized value of VM, thus changing the activity from tonic spiking to bursting. (C) The increase of BhCaS to 0.056 V expands the interburst interval. (D) BhCaS=0.06 V brings the value of the interburst interval close to the targeted value. The leak current parameters are the same as in Fig. 1.1 Bh was 0.031 V.

25 Figure 1.3: Dependence of the burst duration, interburst interval and spike frequency of the bursting activity on BhCaS. Dependences of these characteristics on BhCaS are depicted in panels A, B and C, correspondingly. The shaded areas limit the values of BhCaS corresponding to the temporal characteristics of the activity measured experimentally. There are two ranges of values of BhCaS, where the interburst interval falls within the scopes measured experimentally (the two grey rectangles in (A)). Other parameters are Bh = 0.031 V, gleak = 15.7 nS, Eleak = -0.0505 V.

26 Figure 1.4: Bursting activity of the canonical 4D model. Parameters of the model are BhCaS=0.06 V, Bh=0.031 V, gleak=15.2 nS and Eleak=-0.0505 V. Burst duration is 6.0 sec, interburst interval is 3.0 sec, duty cycle is 66.4 %, number of spikes is 35, frequency is 5.7 Hz. The insert from the 2D diagram allows us to take a close look at the bifurcation that has occurred in the system (Figure 1.6). On the 2D diagram we map the areas of the hyperpolarized and depolarized silence (equilibria), and the areas of tonic spiking (a stable periodic orbit), bursting activity, and various types of multistability. These areas are determined by these four types of codimension-one bifurcations: the Andronov-Hopf, the saddle-node and homoclinic bifurcations of equilibria, and the saddle-node bifurcations of periodic orbits. The stable depolarized equilibrium loses stability through the supercritical Andronov- Hopf bifurcation. At the critical value, it gives rise to the stable periodic oscillations, see Figures. 1.5, 1.7 and 1.8. This periodic orbit represents the tonic spiking activity of the neuron. The bifurcation occurs at the curve labeled AH1 in the (gleak, Eleak) diagram in Figure 1.5. At the bifurcation, the periodic orbit is born with zero magnitude and a nonzero frequency ω,

27 determined by the imaginary part of the characteristic exponents of the equilibrium state. The sign of first Lyapunov coefficient determines the stability of the new-born periodic orbit (Shilnikov et al., 1998, 2001). It is negative for the supercritical Andronov-Hopf bifurcation. The corresponding bifurcation curve corresponds to the transition from depolarized silence into tonic spiking activity. Figure 1.5: The (gleak,Eleak)-parameter bifurcation diagram for the oscillatory and stationary regimes of the model. The dark blue curve AH1 marks the supercritical Androvov-Hopf (A-H) bifurcation of the depolarized equilibrium. To the left of AH1, the equilibrium is stable. To the right, it becomes unstable,

28 giving rise to stable tonic spiking. The orbit of spiking loses stability at the period doubling bifurcation, the green curve PD. Followed further, the tonic spiking periodic orbit disappears through a homoclinic bifurcation of an equilibrium marked by the red curve H1. The A-H bifurcation curve (AH2) for hyperpolarized equilibrium is shown in light blue. On this curve, the two points AHG1 and AHG2 mark the Bautin bifurcations (’⋄’). They bound the section of the curve where the A-H bifurcation is supercritical and gives rise to stable subthreshold oscillations. The outer sections, above AHG1 and below AHG2, mark the subcritical A-H bifurcation, giving rise to the saddle orbit. The range where the saddle orbit exists is bounded by the homoclinic bifurcation of the saddle equilibrium, the light brown curve H2. Passage through the supercritical section leads to the onset of stable subthreshold oscillations. These oscillations vanish through a saddle-node bifurcation of periodic orbits on the dashed black curve SNo2. The area supporting bursting is obtained numerically (mapped in orange, light-blue, and partially in pink) and its border is marked by ’+’s. This border is bounded by the curves PD and H2. In the pink zone there coexist bursting and stable subthreshold oscillations (Figure 1.8). The bright blue patch corresponds to the bistability of bursting and silence; it is bounded by the curves AH2, H2 and H1. The yellow area between the curves AH2 and PD corresponds to the coexistence of the tonic spiking and the hyperpolarized silent regime. The dotted lines indicate the four levels of Eleak used in the diagrams: -0.048 V (1, Figure. 1.7), - 0.04938 V (2, Figure. 1.7), -0.04958 V (3, Figure 1.8), and -0.0505 V (4, Figure 1.8). The dark brown curve, SNe corresponds to the saddle-node bifurcation at which the hyperpolarized equilibria disappear (Figures 1.7 and 1.8). The green ’⋆ ’ locates a point of tri-stability (gleak = 15.4 nS, Eleak= -0.0502 V) illustrated in Figure 1.10. The red ’⋆ ’ locates a point of bistability gleak= 15.70 nS, Eleak= -0.0505 V) shown in Figure 1.11.

29 Figure 1.6: Inset from the (gleak Eleak)-parameter bifurcation diagram for the oscillatory and stationary regimes of the model. The notations are the same as in Figure 1.5.

30 Figure 1.7: The gleak-bifurcation diagrams for two values of the leak reversal potential. Parameter Eleak =−0.048 V in panel (1) and Eleak = −0.04938V in panel (2). The solid green and dashed red curves represent the membrane potential at the equilibria, stable and unstable, correspondingly. Upper and lower branches correspond to the depolarized and hyperpolarized states of the interneuron. The branches are bridged by the middle saddle equilibrium, bounded by two folds corresponding to the saddle-node bifurcations. The tonic spiking emerges through a supercritical Andronov-Hopf bifurcation, AH1, that makes the depolarized branch unstable (dashed). Solid violet branches depict the minimum and maximum

31 membrane potential values of the tonic spiking oscillations. The solid brown curve between them corresponds to the average membrane potential of tonic spiking oscillations. Tonic spiking oscillations double the period after the period doubling bifurcation, PD (vertical green line), and disappear through the homoclinic bifurcation, H1. The hyperpolarized equilibrium loses stability through the Andronov- Hopf bifurcation (AH2): subcritical in (1) and supercritical in (2). (1): The unstable periodic orbit increases and terminates through a homoclinic bifurcation (H2). (2): The magnitude of stable oscillations increases as gleak is decreased up to the saddle-node bifurcation, SNo2. After the fold, the branch of unstable subthreshold oscillations is continued as gleak is increased until a homoclinic bifurcation (H2) where it terminates. The pair of dashed blue curves depicts the minimal and maximal values of the membrane potential of the unstable subthreshold oscillations. In both (1) and (2), the yellow rectangle indicates the range of bistability of tonic spiking and hyperpolarized quiescence. In (2) SNo2 and AH2 bound the range of bistability of tonic spiking and the stable subthreshold oscillations (pink area). Numbers 1 and 2 correspond to the dashed lines 1 and 2 in Figure 1.5.

32 Figure 1.8: Bifurcation diagram showing the evolution of equilibria and oscillatory regimes for two values of the leak reversal potential. Parameter Eleak = −0.04958 V in panel (3) and parameter Eleak = −0.0505 V in panel (4) plotted against the bifurcation parameter gleak. Labeling is the same as in Figure 1.7. Numbers 3 and 4 correspond to the dashed lines 3 and 4 in Figure 1.5. Here, the blue rectangles determine the range of bistability of the bursting and hyperpolarized equilibrium. In (3) the pink rectangle marks the range of the coexistence of bursting and the stable subthreshold oscillations. The critical values for the Andronov-Hopf bifurcation (AH2) and the saddle-node bifurcation SNo1 are very close to each other and marked by the single vertical line AH2|SNo1.

33 A large area in the middle of the bifurcation diagram corresponds to the tonic spiking regime of the model. As the parameter gleak is increased from the bifurcation value, the magnitude of the stable periodic orbit increases. The orbit loses its stability through a period-doubling bifurcation on the curve PD in Figure 1.5. With gleak increased further, the unstable orbit disappears through the homoclinic bifurcation. At this value the orbit becomes a homoclinic loop of the saddle equilibrium. The bifurcation occurs with critical values of gHleak and EHleak which are marked by the curve H1 on the diagram (Figure 1.5). More precisely, to obtain this curve we exploited the fact that the period of the orbit grows as −ln |gleak – gHleak| near the homoclinic bifurcation, thus it can be arbitrarily large in the vicinity of the bifurcation (Shilnikov et al., 1998, 2001). The curve H1 marks the parameters’ values corresponding to tonic spiking with a 25 second period. The area between the curves AH1 and PD supports the periodic tonic spiking in the model. Depending on the level, i.e., the value of the other parameter Eleak, a further increase of gleak beyond the border H1 rightward leads to a transition from tonic spiking into either hyperpolarized silence or bursting (Figure 1.5). At a large value of gleak the neuron stays silent at the hyperpolarized equilibrium, which is a stable focus. Lowering gleak makes it unstable through another Andronov-Hopf bifurcation defining the curve AH2 in the parameter plane in Figure 1.5. For the most part, the bifurcation is a subcritical one which means that it gives rise to an unstable, subthreshold periodic orbit of a saddle type. However, the segment on the curve AH2 bounded by the points labeled AHG1 and AHG2 corresponds to a supercritical Andronov-Hopf bifurcation, giving rise to the stable subthreshold oscillations. At these two points the first Lyapunov coefficient is zero. Each point identifies a Bautin bifurcation, which locates the birth of the saddle-node periodic orbit with a zero amplitude and a non-zero frequency and a codimension 2. The feature of such a bifurcation

34 is that its unfolding includes a bifurcation curve of saddle-node periodic orbits that originates from the codimension 2 point. These two points lay on the intersection of the Andronov-Hopf bifurcation AH2 and the saddle-node bifurcation for periodic orbits SNo1 (Figure 1.5 and Figure 1.6). On the curve SNo1, two subthreshold periodic orbits, stable and saddle, coalesce and vanish. Outside the interval between AHG1 and AHG2, the periodic orbit emerges unstable through the subcritical Andronov-Hopf bifurcation at AH2. With gleak increased, it terminates at the homoclinic orbit of the saddle equilibrium. This event defines the curve H2 (Figure 1.5, Figure 1.6). The bursting activity is observed in the model for quite a large range of the leak current parameter values (orange and blue areas in Figure 1.5). Of special interest here are the boundaries of the area, which are marked with +’s. The onset of bursting appears to be in association with a rapid period doubling cascade leading to chaos in the model (Terman, 1992). The final event of the scenario involves a homoclinic bifurcation of an unstable periodic orbit becoming the homoclinic loop of the saddle equilibrium. The reduced 4D model also demonstrates the coexistence of long-period, irregular bursting with chaotic tonic spiking oscillations. This bistability is observed in the narrow stripe bounded by the curves PD and H1. This type of transition has been described for the Hindmarsh-Rose model (Shilnikov and Kolomiets, 2008). The model shows six types of multistability: (1) tonic spiking and the hyperpolarized silent regime; (2) tonic spiking and subthreshold oscillations; (3) tonic spiking and bursting; (4) bursting and subthreshold oscillations; (5) bursting and the hyperpolarized silent regime; (6) bursting, subthreshold oscillations, and silence.

35 In other words: (1) The basin of tonic spiking and the hyperpolarized silent state are separated by the stable manifold of the saddle equilibrium. In the diagram the zone of this bistability is bounded by the subcritical Andronov-Hopf bifurcation curve AH2, and the period-doubling bifurcation curve PD (Figure 1.5, Figure 1.7 (1) for Eleak = -0.048 V). (2) The tonic spiking and subthreshold oscillations are separated by the stable manifold of the saddle equilibrium. This bistability area is determined by the range of parameter values supporting the stable subthreshold oscillations, i.e. it is limited by the saddle-node bifurcation curve for the periodic orbits, SNo2 and the supercritical Andronov-Hopf bifurcation curve AH2 (Figure 1.6, Figure 1.7 (2) for Eleak = -0.04938 V). (3) The tonic spiking and bursting coexist near the transition border between them (Figure 1.5). We leave this phenomenon outside the scope of this paper for future investigation (Figure 1.9 C- D). (4) The bursting and subthreshold oscillations are separated by the stable manifold of the unstable subthreshold oscillations. Similarly to (2), the area of this bistability is determined by the range of parameter values supporting the stable subthreshold oscillations, and hence is bounded by the saddle-node bifurcations’ curves for the periodic orbits, SNo2 and SNo1 (Figure 1.6, Figure 1.8 (3) for Eleak = -0.04958 V, Figure 1.9 A-B).

36 Figure 1.9: Voltage traces of the coexistence of bursting and subthreshold oscillations and the co- existence of bursting and tonic spiking. A-B) Co-existence of bursting and subthreshold oscillations for gleak=12.96 nS and Eleak=-0.04958 V. C-D) Co-existence of bursting and tonic spiking for gleak=12.13195 nS and Eleak=-0.0505 V. (5) The bursting and the hyperpolarized silent regime are separated by the stable manifold of the unstable subthreshold periodic orbit similar to (4). The corresponding zone is bounded by the range of parameters’ values supporting the unstable subthreshold oscillations, i.e. the area is bounded by the subcritical Andronov-Hopf bifurcation curve AH2, and the homoclinic bifurcation curve H2 (Figure 1.6, Figure 1.8 (4) for Eleak = -0.0505 V, Figure 1.11). (6) The bursting, subthreshold oscillations and silence are separated by the stable manifolds of the two saddle orbits. One saddle orbit appears through the subcritical Andronov-Hopf bifurcation and disappears through a saddle-node bifurcation for orbits SNo1. At SNo1 the stable

37 subthreshold oscillations appear (the stable periodic orbit). In turn they disappear at a saddle node bifurcation SNo2, for a smaller value of the bifurcation parameter, where the second saddle orbit appears, which creates the barrier between the bursting and the stable subthreshold oscillations. The second unstable orbit disappears at the homoclinic bifurcation. Thus the tri- stability is bounded by either the Andronov-Hopf bifurcation or the second saddle-node bifurcation SNo2 on one side and by the first saddle-node bifurcation on the other side SNo1. We illustrated the three regimes by switching the neuron’s activity from the subthreshold oscillations in either bursting or silent regimes (Figure 1.10). 2.1.4 Switching Between Bursting and Silent Regimes by a Pulse of Current In the previous section, we applied the bifurcation analysis to find the range of parameters’ values where bursting and silent regimes co-exist. That analysis predicts that a model with parameters taken from this range starting with different initial conditions would demonstrate one regime or the other. These results alone leave a concern as to whether these regimes would be observable. It might be difficult to demonstrate both regimes if the basin of attraction of one of them is too small. For example, for Eleak=-0.0505 V, the coexistence is determined for gleak between 15.466 nS and 15.776 nS. If gleak is picked in the vicinity of the Andronov-Hopf bifurcation, the basin of attraction to the equilibrium is small, and it is hard to achieve the switch from bursting into a silent regime. If gleak is picked close to the homoclinic bifurcation, then the situation reverses and it is harder to switch from the silent regime into bursting. With the following few numerical experiments we seek to test whether it is feasible to demonstrate both regimes. For this section, we have chosen a somewhat intermediate parameter value for the leak conductance gleak=15.7 nS, with Eleak=-0.0505 V.

38 Figure 1.10: Stable subthreshold oscillations coexist with bursting and silent attractors. By a pulse of current one can switch the activity from the subthreshold oscillations into either of the coexisting regimes. The parameters are gleak= 15.4 nS and Eleak= -0.0502 V (marked by the green ⋆ on Figures 1.5 and 1.6). Initial conditions leading to the stable subthreshold oscillations presented are [V mCaS hCaS hNa]=[- 0.04671933 0.5275212 0.01250879 0.9996319]. The pulses were delivered at 38.9 sec with the amplitudes -0.07 nA (A) and -0.047 nA (B). The question is, how to test for multistability? Although we address this question using a rather abstract model, we have a preference for those perturbations which can be easily implemented in an experimental set up, so that the predictions made could be tested in living neurons. A procedure of this sort has been utilized by Guttman, Lewis and Rinzel (1980). Following their approach, we used perturbations of the model made by ”injection” of a square pulse of current into the neuron. We tested whether such perturbation could be used to switch

39 from one regime to the other. The duration of the pulses was set to 0.03 sec, which is close to the effective width of a synaptic current pulse in the leech heart interneuron. We tried different amplitudes of the pulse. To switch from the silent regime into bursting, one needs to provide a pulse of current with sufficiently large amplitude. The pulse can be either depolarizing or hyperpolarizing. For each polarity the critical value of the pulse which is just sufficient for the switch is different. Thus, there are two thresholds for the switch from silence into bursting. The description is more complicated for the switch from bursting into silence; and one more parameter of stimulation has to be taken into account. It is a phase within the period of bursting. Interestingly enough, for the parameters chosen one could switch from bursting into silence again by either a hyperpolarizing or depolarizing pulse of current (Figure 1.11). An important question for experimental testing is whether these unstable oscillations could be recorded. The model study predicts that it might be possible. A projection of a trajectory onto (mCaS, hCaS)-plane before and after a pulse was injected into the neuron showed that the phase point oscillated in the vicinity of the unstable subthreshold orbit sufficiently long to trace a few cycles of the unstable oscillations. Here, the green dot represents the stable hyperpolarized equilibrium, whose basin of attraction is bounded by the stable manifold of the unstable periodic orbit. The orbit is shown in red. Application of a hyperpolarizing Iinj = -0.029 nA disturbs the neuron but the state is still in the basin of attraction of the equilibrium. However, application of a pulse of a larger magnitude Iinj = -0.03 nA moves the model outside of the basin of the silent attractor away into bursting (Figure 1.12). These two values specify a threshold in terms of the critical amplitude of the hyperpolarizing pulse of current.

40 Figure 1.11: Examples of switching between bursting and silence produced by square pulses of injected current. The switch from silence to bursting is depicted in panel (A,B) and the switch from bursting to silence is shown in panel (C,D). Stimulation of the neuron in the silent regime by a pulse of either depolarizing current, Iinj = 0.61 nA (A) or hyperpolarizing current, Iinj= -0.42 nA (B), switched the regimes. For this parameter regime, application of either a depolarizing (Iinj = 0.05 nA) (C) or hyperpolarizing pulse (Iinj = -0.05 nA) (D) of the current could produce a switch from bursting into silence. The parameters of the model are BhCaS = 0.06 V, Bh = 0.031 V, gleak = 15.7 nS, Eleak = -0.0505 V. The two red dashed lines mark maximal and minimal membrane potentials of the subthreshold oscillations.

41 Figure 1.12: Revealing of the unstable oscillations. A pulse of current with a critical value of amplitude applied to a neuron in the silent regime allows us to record oscillations in close proximity to the saddle periodic orbit. The parameters are the same as in the figure 1.10 except for gleak=15.55 nS and Eleak=- 0.0505 V. The grey dot represents the stable equilibrium point. The grey dotted curve marks the unstable periodic orbit. Shown on the right are the corresponding traces after the pulses of magnitudes Iinj=-0.029 nA and Iinj= -0.03 nA. The vertical dashed lines mark the time interval where the frequency of the unstable subthreshold oscillations and resonant frequency of the rest state are calculated.

42 2.1.5 Conclusion I Summing up, we observed up to six types of multistability in {INa,ICaS} model. The separating regimes are based on a saddle orbit and on a saddle equilibrium. The mechanisms supporting multistability are described among five different types of multistability and are presented in Table 1. (Malashchenko et al., submitted to PloS One). Table 1: Different types of multistability in {ICaS INa} model; separating regimes and mechanisms supporting multistability. Type of multistability Separating regime Regime appears Regime disappears Tonic spiking & Saddle equilibrium Limited by Andronov- Limited by Homoclinic silence Hopf Tonic spiking & Saddle orbit Andronov-Hopf Saddle-node subthreshold oscillations Tonic spiking & NA bursting Bursting & Saddle orbit Andronov-Hopf Saddle-node subthreshold oscillations Bursting & Two saddle orbits Andronov-Hopf or Saddle-node subthreshold Saddle-node oscillations &silence Bursting & silence Saddle orbit Andronov-Hopf Homoclinic

43 Chapter II: Similarity Between the Half-Center Oscillator and the 14D Models Dynamics Multistability is ubiquitously observed in neuronal networks and in single neurons. In such systems several observable regimes coexist without change to any biophysical properties. The activity of a bistable neuron or neuronal network can be switched from one state to another by a transient perturbation such as a pulse of synaptic current or random noise. The mechanisms supporting bistability are important to understand from the perspective of its functional roles and the diseases in which it plays a role. The functional role of bistability is implicated in controlling locomotion (Carlin et al., 2000; Perrier and Hounsgaard, 2000) and short and long-term memories (Tahvildari et al., 2007; Durswewitz and Seamans, 2006; Marder et al., 1996). On the other hand, bistability of neuronal networks can accompany pathological conditions such as epilepsy (Hahn and Durand, 2001; Foss and Milton, 2003; Fröhlich and Bazhenov, 2006; Takeshita et al., 2007). It still remains unknown how bistability occurs in biological systems, and what are the controlling factors, and what are the mechanisms supporting it. Preparations with invertebrates provide the advantage of gaining insights into the dynamics of a bistable system on the cellular level. The existence of bistability was predicted theoretically in various cells and later confirmed there experimentally. Classic examples are the theoretical studies that revealed the bistability of tonic spiking and silence in the squid giant axon (Rinzel 1978; Best, 1979) and bistability of tonic spiking and bursting in Aplisia (Canavier et al., 1994). We have shown the bistability of the bursting and silence for the canonical model of the leech heart interneuron (Malashchenko, Master Thesis, 2007). Here we address two questions which were left unanswered before. First, we investigate the strength-duration properties of the current pulses switching silence into bursting. Second, we explore whether the half-center oscillator configuration would protect normal pattern generation against perturbations.

44 2.2.1 14D Model of Leech Heart Interneuron A leech heart interneuron model is represented by a 14D system of ordinary differential equations. It describes a single isopotential compartment with eight voltage-gated membrane conductances and a leak current. Voltage-gated conductances are represented in terms of the Hodgkin and Huxley formalism (Hodgkin and Huxley, 1952). The sum of currents passing though the membrane is described by the equation dV C ( I Na I P I CaF I CaS I h I K 1 I K 2 I KA I leak I inj (t )) , dt where C is the total membrane capacitance (0.5 nF), Iion is an intrinsic voltage-gated current, Iinj is the injected current and Ileak is the leak current. The voltage-gated currents are given by I Na g Na mNa hNa (V E Na ); 3 I P g P mP (V E Na ); I CaF g CaF mCaF hCaF (V ECa ); I CaS g CaS mCaS hCaS (V ECa ); 2 2 I K1 g K1mK1hK1 (V E K ); I K 2 g K 2 mK 2 (V E K ); I KA g KA mKA hKA (V E K ); 2 2 2 I h g h mh (V Eh ); I leak g leak (V Eleak ); I inj (t ) I inj f (t inj ); 2 where g ion is the maximal conductance, and Eion is the reversal potential of a specific ionic current. These currents include: fast calcium current-ICaF, slow calcium current-ICaS, fast sodium current-INa, delayed rectifier-potassium current-IK1, persistent potassium current-IK2, fast transient potassium current-IKA, hyperpolarization activated cationic current-Ih, persistent sodium current-IP, and leak current- Ileak. The value of function f (t inj ) is 1 when the pulse is injected; otherwise f (t inj ) is 0, and m and h are the activation and inactivation gating variables,

45 respectively. These variables describe the dynamics of ionic channels and are governed by the following equations dm K 2 f ( 83, 0.02, V ) m K 2 , dt 200, 0.035, 0.057, 0.043,V dm P f ( 120, 0.039,V ) m P , dt ( 400, 0.057, 0.01, 0.2,V ) dm Na f ( 150, 0.029, V ) m Na , dt 0.0001 dhNa f (500, 0.030, V ) hNa , dt hNa (V ) dmCaF f ( 600, 0.0467, V ) mCaF , dt mCaF (V ) dhCaF f (350, 0.0555, V ) hCaF , dt ( 270, 0.055, 0.06, 0.31, V ) dmCaS f ( 420, 0.0472, V ) mCaS , dt ( 400, 0.0487, 0.005, 0.134, V ) dhCaS f (360, 0.055, V ) hCaS , dt ( 250, 0.043, 0.2, 5.25, V ) dm K 1 f ( 143, 0.021, V ) m K 1 , dt (150, 0.016, 0.001, 0.011, V ) dhK 1 f (111, 0.028, V ) hK 1 , dt ( 143, 0.013, 0.5, 0.2, V ) dm KA f ( 130, 0.044, V ) m KA , dt ( 200, 0.03, 0.005, 0.011, V ) dhKA f (160, 0.063, V ) hKA , dt ( 300, 0.055, 0.026, 0.0085, V ) dmh f h (V ) mh dt ( 100, 0.073, 0.7, 1.7, V ) where the steady-state activation and inactivation functions f ( A, B, V ) are given by a function 1 f (a, b,V ) except for the function of the steady-state activation of Ih which is 1 e a (V b ) 1 described by f h (V ) 180(V 0.047) . Time constants of the gating variables are 1 2e e 500(V 0.047)

46 d governed by the function, (a, b, c, d ,V ) c . The exception includes the 1 e a (V b ) inactivation time constant of INa and the activation time constant of ICaF, which are described by the following functions. 0.006 0.01 hNa (V ) 0.004 500(V 0.028) 1 e cosh(300(V 0.027)) 0.024 mCaF (V ) 0.011 cosh(330(V 0.0467)) The parameter values of the 14D model are chosen from Cymbalyuk et al. (2002). The reversal potentials are ENa= -0.045 V, ECaS= 0.135 V, EK= -0.07 V, Eh= -0.021 V, and the maximal conductances for a single neuron are g Na = 200 nS, g P 7 nS, g CaF 5 nS, g CaS 3.2 nS, g K1 100 nS, g K 2 80 nS, g KA 80 nS, g h 4 nS. The conductance and reversal potential of leak current, g leak and Eleak, were varied since these parameters were used as the controlling parameters in bifurcation analysis of the model. 2.2.2 Switching Between Bursting and Silent Regimes Previously we investigated a mechanism supporting the bistability of bursting and silence in the model of the leech heart interneuron (Malashchenko, Master Thesis, 2007). By performing a bifurcation analysis of the model, we have shown that for a specific set of parameters of conductance and reversal potential of a leak current (gleak, Eleak) a silent regime coexists with bursting regime, and the result was mapped on a 2D diagram of stationary and oscillatory states in the model (Figure 4). An example of bistability shown in Figure 2.1 illustrates the switch between attractors by an appropriately timed perturbation with a negative square pulse of current.

47 The bistability of bursting and silence is based on a saddle periodic orbit. The stable manifold of this orbit separates two basins of attractions and represents a threshold between them. The saddle orbit is shown in red on a 3D projection of the phase space (mCaS, hCaS, V) in Figure 2.1 A,C. The unstable subthreshold oscillations could be recorded for several seconds if the phase point was placed very close to this saddle orbit by a pulse of current (Figure 2.1 B,D). These oscillations approximate the saddle orbit. The perturbation by a negative pulse of current moved the phase point of the system across the threshold into the basin of attraction of the stationary state, and the bursting activity was annihilated (Figure 2.1 A,B). The switch was done by a pulse of current with amplitude 0.1044 nA. If during the perturbation, the phase point did not cross the threshold, the system exhibited transient subthreshold oscillations and then resumed bursting. To demonstrate this, the amplitude of the negative square pulse was slightly increased to 0.1045 nA (Figure 2.1 C,D). Perturbations of the initially silent neuron by the depolarizing and hyperpolarizing square pulses of current, Iinj, showed that there are two thresholds for the amplitude of the pulse of current for the switch from silence to bursting by the negative and positive Iinj. Repeating this analysis we found the two thresholds for pulses with different durations (Figure 2.2). The longer the duration of the pulses was the smaller the two thresholds were. The two curves obtained were I rheobase fitted to the Lapicque formula strength-duration curve I thresh t inj (Brunel and Rossum, tm 1- e 2007), where Irheobase is the asymptotic current's threshold value for the infinitely long pulse, tinj is the pulse duration, tm is a time constant characterizing the process of charging cell's membrane. The leech heart interneurons are organized in an oscillatory neuronal network, a central pattern generator. Within it, the interneuron pairs, located in ganglia 3 and 4, form half-center

48 oscillators that pace the rhythmic activity of the central-pattern generator. Therefore our finding raises the question of whether bistability can be observed in a half-center oscillator under the condition that single neurons are bistable. We considered the model of a half-center oscillator that consists of two identical neurons coupled through synaptic inhibitory connections. The synaptic currents are introduced by two components: graded and spike-mediated (Cymbalyuk et al., 2002). Each neuron had g leak = 11.3 nS, so that it would exhibit bistability. The equations describing the membrane potential of each neuron and the synaptic currents coupling neurons are provided in the Appendix. Interneurons of the half-center oscillator produce alternating bursting. We expected that if the activity of one bistable neuron is switched into a silent regime during the burst phase, the second bistable neuron being inhibited at that time could turn into silent regime as well. Thus, the regime of the half-center oscillator could be switched from bursting to silence by a perturbation applied to just one neuron. To test this hypothesis, first we chose the initial conditions of the membrane potential and gating variables of both neurons that led to the bursting activity of the half-center oscillator. We integrated the activity for 1000 seconds to confirm that bursting activity was the attracting regime. For g leak = 11.3 sec the bursting activity was the attractor, but small irregularity of temporal characteristics of the activity was noticed. For example, the period of the bursting activity in both neurons was around 7.2 ± 0.24 sec, and the spike frequency was 8.37 Hz ± 0.419 Hz. A pulse of current was applied during the spiking phase of the first neuron (labeled as HN(L) in Figure 2.3). Figure 2.3 shows that the perturbation of a single neuron by a negative pulse with the amplitude 0.03 nA and the duration of 0.03 sec switched the regime of activity from bursting to silence. The switch of activity from bursting to silence depends not only on the characteristics of the pulse and the timing when the pulse is

49 applied, but also on the dynamics of the synaptic currents such as the strength of coupling and its time constants. Figure 2.1: Unstable sub-threshold oscillations recorded after the perturbation by a pulse of current. Bursting is shown in dark blue. The 3D projection shows unstable oscillations (red circle) located between two basins of attraction: bursting mode (dark blue trace) and silent mode (green dot). A)B) The switch from bursting into silence is performed by a negative pulse of current with the amplitude 0.1044 nA C)D) A negative pulse with the amplitude of 0.1045 nA does not produce the switch from bursting to silence. Parameters of the model were kept the same as in Cymbalyuk et al (2002), except for gleak=10.7 nS. The pulse of current was applied near the first spike of a burst.

50 Figure 2.2: The strength-duration relationship of the minimum depolarizing (grey) and hyperpolarizing (black) square pulses of current switching the silent regime into bursting regime. To toggle the switch the pulse amplitude must be larger than the thresholds determined by the curves. The parameters of the leak current are the same as in the figure 2.1 (Malashchenko et al., submitted to PRE).

51 Figure 2.3: The switch of bursting into silence in the half-center oscillator (HCO) by a pulse of current delivered to one neuron from the pair. Stimulation of one neuron by a hyperpolarizing pulse of current seized the activity of both neurons. The amplitude of the negative pulse was 0.03 nA. Parameters of the leak current were gleak=11.3 nS, Eleak=-0.0635 V. Initial conditions of HCO are given in the Appendix. 2.2.3 Transient Bursting Activity At the border of the transition from bursting to silence the canonical model of the leech heart interneuron exhibited long-lasting transient bursting activity. In order to locate the area of bistability in the (gleak, Eleak) bifurcation diagram, we needed to detect the exact border indicating the transition from bursting to silence while changing the gleak and Eleak parameters. For values of gleak larger than 10.84 nS the neuronal model exhibited only hyperpolarized silence, the value 10.84 nS indicated the transition point from bursting to silent regimes. In order to calculate the

52 parameter values corresponding to the transition from bursting to silence, we took into account that the system is bistable. Its activity depends on the choice of initial conditions of all state variables. In the (gleak, Eleak) bifurcation diagram each point indicating the transition from bursting activity to silence was found by numerical integration of the system for successive parameter values. We performed the following steps. The initial conditions and gleak corresponding to bursting activity were picked and, in order to confirm that bursting activity is the attractor, the system was integrated for 2000 seconds. If bursting was the attractor, we picked the initial conditions corresponding to the last point from the newly integrated trajectory. Then, we slightly increased gleak and performed integration for 2000 seconds again. If bursting was still the attracting regime, we increased gleak and the procedure was repeated until the transient switch from bursting to silence occurred. The detected gleak was mapped in the 2D diagrams as a transition point from bursting to silence. We assigned to different values Eleak and, by using the described procedure, calculated the corresponding gleak values for the transition from bursting to silence. The border was located next to a homoclinic bifurcation curve that limited the range of parameter values of the existence of the saddle orbit in (gleak, Eleak) parameter space. It is noteworthy that the model of the interneuron can exhibit bursting activity for several hundreds of seconds and then bursting suddenly falls into silence if the parameters' values are close to the homoclinic bifurcation values. In Figure 2.4 (left panel) the neuronal model exhibits bursting activity for 500 sec and after 500 sec the activity ceases (Malashchenko, Master Thesis, 2007). The transient bursting activity can be induced again by application of a perturbing pulse of current. The similar behavior was observed in the model of HCO (Figure 2.4, right panel). If g leak was set to 11.6 nS, the activity of HCO ceased after around 500 sec but was reiterated by the perturbation (Figure 2.4).

53 Figure 2.4: Long-lasting transient transitions from bursting to silence. The left panel shows an example of the spontaneous transition from the bursting activity to silence occurring in the single interneuron after about 500 sec. The leak current parameters were gleak=10.8337 nS, Eleak=-0.0635 V. At this value of gleak the system is located near the homoclinic bifurcation for the saddle periodic orbit. Homoclinic bifurcation appears at gleak=10.87 nS. The right panel depicts the transient bursting activity of the half-center oscillator. The bursting activity terminates in both neurons, HN(L) and HN(L), after about 500 sec. The leak current parameters are gleak=11.6 nS, Eleak=-0.0635 V.

54 2.2.4 Conclusion II We investigated the strength-duration properties of the current pulses switching silence into bursting in the 14D model of the leech heart interneuron. We showed that strength-duration dependence could be curve fit into Lapicque's formula. We showed that the switch by a pulse of current from bursting to silence could be toggled not only in the isolated leech heart interneuron model but also in the half-center oscillator model. In addition to previously reported long transient activity in the model of the single leech heart interneuron (Malashchenko, Master Thesis, 2007), here, we demonstrated that the model of the half-center oscillator also exhibited transient bursting activity in vicinity of the gleak parameter where the transition from bursting to silence occurred. Chapter III: Propensity to Bistability of Bursting and Silence of the Leech Heart Interneuron Here, we focused on a specific type of bistability: the co-existence of bursting activity and silence in a single neuron. In (Malashchenko, Master Thesis, 2007) we predicted that the leech heart interneuron can exhibit bistability of bursting and silence under conditions with elevated leak conductance, although the range of parameters leading to bistability is very narrow. In this chapter we explored the conditions that could increase the propensity of a neuron to bistability. We developed a general method that allows one to identify how up- and down- regulations of currents affect neuronal bistability. Such analysis predicted that under pharmacological conditions reported in this section the bistability of bursting and silence can be revealed in a single leech heart interneuron.

55 2.3.1 Current Conductances Differently Affect Bistability The bifurcation analysis of the bistable system demonstrated that the range of bistability is limited by values of gleak corresponding to the Andronov-Hopf and homoclinic bifurcations at which this saddle orbit originates and terminates (Malashchenko, Master Thesis, 2007). We used this knowledge to explore the effect of each ionic current on evolution of this orbit. We traced the Andronov-Hopf and homoclinic bifurcations and mapped them on 2D diagrams. gleak and a maximum conductance of an ionic current were used as bifurcation parameters. The purpose was to examine the relative position of the sub-critical Andronov-Hopf (AH) and homoclinic (Hom) bifurcation curves while changing gleak and one maximum ionic conductance at a time. In this way, we evaluated how up- and down-regulation of an ionic conductance affected the range of the leak conductance where the unstable orbit existed. On the diagrams we mapped the curve of the transition from bursting to silence (BSi), calculated according to the procedure described in session 2.2.3. For a given ionic conductance value, the difference between gleak, corresponding to the transition from bursting to silence (BSi), and gleak, corresponding to the sub-critical Andronov-Hopf bifurcation (AH), determined the gleak interval where bistability of bursting and silence was observed. We suggested that the interval of the gleak range between the BSi and AH curves can be used as an index the evaluate the neuronal propensity for bistability. The curve of the transition from bursting to tonic spiking activity (BTs) was also calculated. We calculated the 2D bifurcation diagram using CONTENT obtained for successively changing parameters in a parameter space A) (gleak, gP), B) (gleak, gNa), C) (gleak, gCaF) and D) (gleak, gCaS), where g p , g Na , g CaF , g CaS are maximum conductances of the persistent Na+, the fast Na+, the low-threshold fast Ca2+ and the low-threshold slow Ca2+ currents (Panels A-D in Figure 3.1). For bifurcation analysis we set the upper and lower boundaries of the ionic

56 conductances so that the maximum of an ionic conductance did not exceed double the canonical value, and the minimum was set to zero nS. If the bursting regime disappeared at a value smaller (larger) than the established maximum (minimum), this value was considered as the upper (lower) boundary. Panel A of Figure 3.1 shows sub-critical Andronov-Hopf and homoclinic bifurcation curves in (gleak, gP) parameter space; the saddle orbit exists for the range of parameters determined by these curves. The corresponding gleak width (horizontal slice) for the given g p values was very narrow and stayed relatively constant as g p was up- or down-regulated from its canonical value of 7 nS. It predicts that the index of propensity is not noticeably affected by the change of g p conductance. To locate the range of bistability, we calculated the transition from bursting to silence (BSi) and thus defined the right-side border of the bistable area. The BSi curve passes very close to the homoclinic curve. The left border of the bistable range is defined by the AH curve. The area of bistability in the 2D diagram was approximately located between Andronov-Hopf and homoclinic curves, and more precisely it was determined by Andronov- Hopf and BSi curves. The change in g p parameter slightly affects the bursting area that is observed for the range of parameters between BTs (transition from tonic spiking to bursting) and BSi curves. The model exhibited only tonic spiking activity with gleak values smaller than determined by the BTs curve. Let us now consider the role of the fast Na+ current in the model. The increase of g Na from the canonical value of 200 nS up to 400 nS slightly decreased the gleak interval between the Adronov-Hopf and homoclinic curves (Figure 3.1 B). Similar to the analysis performed for g p in Figure 3.4 A we calculated the BTs, AH, BSi, and Hom curves. According to the diagram in

57 Figure 3.1 B the bursting area and the area of bistability were little altered with the change of g Na . The effect of up- or down-regulation of Ca2+ currents on the AH and Hom bifurcation is shown in Figure 3.1 C,D. The down-regulation of g CaF increased the gleak interval between the two bifurcation curves, whereas the down-regulation of g CaS produced the opposite effect. Considering the gleak interval between the AH and BSi curves, the range of bistability increased as g CaF decreased, and in contrast the range of bistability decreased as g CaS decreased. Furthermore, we investigated the effect of four other currents. Figure 3.2 shows 2D bifurcation diagrams in a different projection of the parameter space A) (gleak, gK1), B) (gleak, gK2), C) (gleak, gKA) and D) (gleak, gh), where g K1 , g K 2 , g KA and g h are maximum conductances of the delayed rectifier K+, fast K+, transient K+, and hyperpolarization activated cationic currents, respectively. The increase of g K1 from 0 to 200 nS had almost no effect on the interval of gleak values between the AH and Hom bifurcation curves (Figure 3.2 A); yet the range of (gK1, gleak) parameters supporting the bursting regime in the model increased drastically as g K1 decreased from the canonical value of 100 nS to 0 nS. This range is located between the BTs and BSi curves. Panels B-C in Figure 3.2 show that the changes of g K 2 and g KA conductances had no significant effect on the gleak interval between AH and Hom bifurcation curves. Similarly, we investigated the effect of g h . The elevation of g h strongly affected the gleak interval between the AH and Hom bifurcation curves. For instance, the increase of g h from the canonical value of 4 nS to 8 nS expanded the gleak interval between bifurcations by a factor of 3.2 (Figure 3.2 D). Moreover, the area between the AH and BSi curves and the area between the BTs and BSi curves increased with the elevation of g h .

58 Let us sum up the data presented in Figure 3.1 and Figure 3.2. The value of the propensity index (Δindex) for the canonical set of parameters corresponds to the width of the range of g leak equal to 0.17 nS. The bifurcation diagrams in Figure 3.1 and 3.2 are difficult to analyze visually in terms of the change of the propensity index. Therefore data from these diagrams were used to construct graphs showing the dependence of the propensity index on the maximum ionic conductances (Figure 3.3 and 3.4). A value of Δindex larger than the canonical one (0.17 nS) indicates that the corresponding voltage-gated conductance increases the range of bistability of bursting and silence. For example, let us consider the dependence of the propensity index Δindex on g P (Figure 3.3A). The value of Δindex stays relatively constant near its canonical value and decreases if g P is up- or down-regulated. A similar analysis established that the increase of g Na from 100 nS to 400 nS decreased the propensity index (Panel B in Figure 3.3). The similar analysis showed that the calcium currents affected the propensity to bistability differently. The propensity index had local maximum near a g CaS value of 5 nS. The index was elevated as well by the complete elimination of g CaF (Panels C and D in Figure 3.3).

59 Figure 3.1: Effect of the variation of the maximum conductances of IP, INa, ICaF, and ICaS currents on the propensity index to the bi-stability of bursting activity and silence. For each value of gi an index of propensity was calculated. The blue curves mark the sub-critical Andronov-Hopf bifurcation (AH), the red curves represent the homoclinic bifurcations (Hom). Between these two bifurcation curves the saddle periodic orbit exit. The dashed yellow curves mark the transition from bursting to silence (BSi). The green curves mark the transition from bursting to tonic spiking (BTs). The horizontal dashed black lines plotted across the diagram correspond to the canonical value of a voltage-gated current conductance. A) B) D) Some effect on the propensity index is detected. C) 2D diagram in parameter space of (gCaF, gleak) shows that the decrease of gCaF increased the index of propensity. Eleak is fixed to -0.0635 V.

60 Figure 3.2: Effect of the variation of the maximum conductances of IK1, IK2, IKA, and Ih currents on the propensity index for bi-stability of bursting and silence. The meanings of colors and types of lines are described in Figure 3.4. A) B) C) No significant change of the propensity index was detected. C) The 2D diagram in parameter space of (gh, gleak) shows that the increase of gh increased the index of propensity.

61 Figure 3.3: Dependence of the propensity index on gP, gNa, gCaF and gCaS. The dashed grey vertical line corresponds to the canonical value of a parameter in the model. A) Manipulation of gP almost did not affect the index of propensity. B) The decrease of gNa to 100 nS slightly increased the index of propensity from 0.17 to 0.2. C) The complete reduction of gCaF increased the index from 0.17 to 0.23. D) The increase of gCaS almost did not affect the propensity index.

62 A similar procedure was used to analyze the effects of different values of maximum conductances of g K1 , g K 2 , g KA and g h on the propensity index. The index of propensity varied with the increase of g K1 from 0 nS to 200 nS, and the dependence was presented by a jagged curve (Panel A in Figure 3.4). The index of propensity was larger near the canonical value of g K1 . The increase of g K 2 slightly decreased the propensity index, whereas the increase of g KA almost did not affect the propensity index (Figure 3.4 B,C). Panel D in Figure 3.4 shows the increase of g h from 2 to 8 nS elevated significantly the index of propensity to bistability. Notice, the increase of g h , compared to other conductances, had the biggest effect on the relative distance between the AH and Hom bifurcation curves. The data presented in Figures 3.3 and 3.4 are summarized in Table 2. It contains the maximum (indexMax) and the minimum (indexMin) values of the propensity index for a given conductance and difference between these values (Δindex). To sum up, we sorted currents into three groups based on the magnitude of the propensity index Δindex: small, intermediate, and significant. Comparing the data in Table 1 we may conclude that the increase of Ip and IKA caused the smallest effect on the index of propensity. Manipulation of the maximum conductances of currents ICaS, ICaF, IK1, IK2, and INa produced an intermediate effect on the propensity index and the up-regulation of Ih led to a significant increase of the propensity index for bistability of bursting and silence. These results leave the question of how ionic conductance modulation would be reflected on 2D diagrams in the leak current parameter space. Here, we suggest that the larger the area within the leak current parameters, the more probable it is to find these neurons bistable.

63 Figure 3.4: Dependence of the propensity index on gK1, gK2, gKA and gh. The dashed grey vertical line corresponds to the canonical value of a parameter in the model. A) The canonical value of gK1 was near the local maximum of the propensity index. B) C) The change of gK2 and gKA did not lead to a significant increase of the index of propensity. D) The increase of gh significantly increased the propensity index. The change of gh from 4 nS to 8 nS increased the propensity index from 0.17 nS to 0.31 nS.

64 Table 2: Minimum and maximum values of the propensity index for the different voltage-gated conductances. Minimum and maximum values are marked in the table by indexmin and indexmax. Corresponding conductances to indexmin and indexmax are marked by g(indexmin) and g(indexmax). Values of indexes and conductances were taken from the diagrams in Figures 3.3 and 3.4. The symbol Δindex corresponds to the difference between minimum and maximum values of indexes Δindex =indexmax- indexmin. Based on the magnitude of Δindex the ionic conductances were classified by effectiveness using the following criteria: small effect if Δindex < 0.05 nS, (grey), intermediate effect if 0.05 nS ≤ Δindex<0.2 nS (bold) and significant effect if Δindex ≥ 0.2 nS (underlined bold). The index of propensity for the canonical parameters corresponds to 0.17 nS. gNA, nS gP , nS gK1, nS gK2, nS gKA, nS gCaF, nS gCaS, nS gh, nS (100-400) (4 -14) (0-200) (30-140) (0-160) (0-10) (2-6.5) (2-8) indexmin 0.1366 0.1366 0.0997 0.1054 0.1571 0.1347 0.1223 0.0751 g(indexmin) 400 14 200 120 160 10 2 2 indexmax 0.1993 0.17 0.1834 0.1832 0.1803 0.2337 0.1929 0.3200 g(indexmax) 100 7 85 40 0 0 5 8 Δindex 0.06 0.03 0.08 0.08 0.02 0.10 0.07 0.24 2.3.2 Changes in 2D Bifurcation Diagram of Stationary and Oscillatory States When Ih and ICaF are Modified The parameters (gleak, Eleak) defining the AH and Hom bifurcations limit the range of parameters (gleak, Eleak) where the bistability of bursting and silence can be observed (Figure 4). The area of bistability defined in the (gleak, Eleak) bifurcation diagram for the canonical set of ionic current parameters was used as the control. We compared the ranges of bistability under different ionic current conditions. We assumed that if up- or down-regulation of ionic currents

65 modifies the area of bistability such that it becomes larger than the control, then it is the more probable to observe bistability in an experimental setting. We computed the bifurcation diagrams in leak current parameter space for the value of g h fixed to 8 nS (Figure 3.5 B) and compared it with the control diagram (Figure 3.5 A). Both diagrams map the areas of parameter values supporting bursting activity and bistability of bursting and silence. These diagrams also depict the sub-critical AH and Hom bifurcation curves. The diagram in Panel B indicates that the elevation of g h significantly increases the area of bistability of bursting and silence in (gleak, Eleak,) parameter space. From the data presented in Table 1 there were five currents characterized as intermediately effective, and the strongest impact among conductances of ICaS, ICaF, IK1, IK2, and INa currents on the propensity index was made by g CaF . Therefore, we investigated whether the superposition of the two modulatory conditions could further increase the range of bistability. We calculated the bifurcation diagram in (gleak, Eleak,) parameter space with the g CaF = 0 nS and g h = 8 nS. As a result, Figure 3.6 shows that the range of bistability is further increased compared to the diagrams in Figure 3.5 A and B.

66 Figure 3.5: Bifurcation diagrams of the oscillatory and rest regimes in the parameter space (gleak, Eleak) for different gh. The red dots are numerically found points where the transition from bursting activity into silence occurred. The sub-critical Andronov-Hopf bifurcation of the hyperpolarized rest state (silent regime) is shown by the light blue curve and marks the boundary where the silent regime loses stability, giving rise to unstable subthreshold oscillations. The green curve corresponds to the homoclinic bifurcation of the unstable subthreshold oscillations. The area between these two curves (the blue and the red curve) marks the parameter regime where unstable subthreshold oscillations exist. The red curve

67 locates a period-doubling bifurcation of the large amplitude periodic spiking and detects the transition from tonic spiking to bursting activity. A) The diagram is computed for gh = 4 nS. B) gh = 8 nS, the increase of conductance of Ih increased the index of propensity and the area of bistability. The width of the bi-stable area was increased by a factor 1.8 for the canonical value of Eleak=-0.0635V. Figure 3.6: Bifurcation diagrams of the oscillatory and rest regimes in parameter space (gleak, Eleak) as gh and gCaF are varied. The boundaries between different oscillatory and stationary states are defined as in Figure 8. The diagram was calculated for parameters gh=8 nS and gCaF=0 nS. The elimination of fast Ca2+ current and the increase of Ih current expanded the width of the bistable area.

68 2.3.3 Conclusion III We introduced a propensity index to bistability as the width of the range of gleak supporting the co-existence of the bursting and silence regimes for a given Eleak We investigated the effect of modulation of conductances of ionic currents on the propensity index. The analysis was performed determining the bifurcations that characterize the appearance and the disappearance of the separating regime. These factors impact the range of parameters where bistability can be observed. Only two currents were highlighted that increased the propensity index bistability of bursting and silence: the hyperpolarization-activated cationic current, Ih, and L-type low threshold fast Ca2+ current, ICaF. 2.4 Chapter IV: Modeling of Respiratory Center Neurons and Revealing Bistability by Current Modulation Previously, we described a mechanism underlying the bistability of bursting and silence (Malashchenko, Master Thesis, 2007). It appears to be generic and we hypothesized that it would be commonly found in many different models of bursting neurons. We applied our techniques to models of vertebrate neurons from the central pattern generator controlling respiratory rhythm in mice and rats (Ramirez et al., 2007; Smith et al., 2000). We analyzed two models, which represent two types of inspiratory neurons of the Pre-Bӧtzinger Complex. These two types are determined through analysis of a neuron's activity in isolation from the network. Neurons of both types are endogenously bursting neurons. The neurons of the first type are shown to be dependent on the persistent Na+ currents, since the blockade of this current by a specific blocker, riluzole, annihilates the neuronal rhythmic activity. The Ca2+ currents are not necessary for

69 producing bursting activity because the Ca2+ currents blockers such as Cd2+ do not abolish bursting activity. The neurons of the second type produce bursting activity that does not require the persistent Na+ current, but does require Ca2+ currents. The bursting activity persists if riluzole is applied, but vanishes in the presence of Cd2+. The first type of neuron is called Cd2+ - insensitive (CI) and the second type of neuron is called Cd2+-sensitive (CS) (Ramirez et al., 2004; Peña et al., 2004). It has been shown that both types of pacemakers generate bursting activity during normoxia, but under the hypoxic conditions only the Cd2+ -insensitive pacemakers are active and continue to drive the respiratory rhythm. The Cd2+ -insensitive pacemakers are critically important for generation of gasping. If these neurons fail to generate rhythm during hypoxia, it might be fatal for the animal. Dysfunction of the CI neurons could cause such pathological conditions as SIDS (Peña et al., 2004; Doi and Ramirez, 2010). In this chapter, we investigate whether the pacemaker neurons of these types could exhibit bistability of bursting and silence. 2.4.1 Model of the Cd2+ - Insensitive Pacemaker The Cd2+ -insensitive pacemaker produces bursting activity as a result of its endogenous dynamics. The primary feature of a Cd2+ -insensitive (CI) neuron is the ability to generate bursting that is dependent on the dynamics of the persistent Na+ current and not Ca2+ currents. It continues bursting activity with Ca2+ currents blocked and ceases with a blockade of the persistent Na+ current (Peña et al., 2004). The period of bursting activity of the CI pacemaker was recorded in the range 0.7 sec and 3 sec (Peña et al., 2004). The hyperpolarization-activated cationic current, Ih, has been found in this type of neuron; and this current plays an important role in determining the period of respiratory neurons (Thoby-Brisson et al., 2000). The blockade of Ih with Cs+ or ZD 7288 does not annihilate bursting activity. The CI pacemaker model was

70 assembled based on the model previously developed by Rybak et al. (2003). We incorporated into it Ih current. The gating variable and time constant of Ih were taken from Amini et al. (2005), and the value of the reversal potential was obtained from Thoby-Brisson et al. (2000). The model’s dynamics are described by four voltage-gated currents: the fast Na+ current, INaf, the persistent Na+ current, INap, the delayed rectifier K+ current IK, Ih, and the leak current, Ileak. The system is described by the following equations: I I Naf Nap K I I h dV C [[ g Naf m Naf h Naf [V E Na ] g Nap m Nap h Nap [V E Na ] [ g K m K [V E K ] [ g h m h [V E h ] 3 4 2 dt I leak g leak [V E leak ] I inj ], dm Naf f (0.0438,0.006, V ) m Naf , dt (0.0009,0.0438,0.014, V ) dh Na f (0.0675,0.0108, V ) h Naf , dt (0.0352,0.0675,0.0128, V ) dm Nap f (0.0471,0.0031, V ) m Nap , dt (0.001,0.0471,0.0062, V ) dh Nap f (0.057,0.006, V ) h Nap , dt (7,0.059,0.006, V ) dm K f (0.0445,0.005, V ) m K , dt (0.004,0.0445,0.01, V ) dm h f h (V ) m h , dt h(V ) where capacitance of the cell membrane, C, is 0.0362 nF, the parameters of conductances and reversal potentials of ionic currents are gNaf= 150 nS, gNap= 4 nS, gK= 50 nS, gh= 2 nS, gleak= 1.2 nS, ENa= 0.05865 V, EK= -0.080 V, Eh= -0.0414 V, and Eleak= -0.0635 V.

71 The steady-state activation and inactivation curves are functions of the membrane potential, V, described by f ( A, B, V ) 1 VA 1 exp( ) B Time constants of gating variables except for the time constant of the Ih current are described by the function ( A, B, C , V ) A , and the time constant of Ih is presented by V B cosh( ) C (V 0.04) 2 h (V ) 0.02 exp( ) 0.09. The steady-state function of m h, 0.0016 1 f h (V ) , where mx is the constant that was used to modify the curve of V 0.075 1 mx exp( ) 0.006 activation. The steady-state activation curves of Ih for two values of mx 1 and 1.5 appear as shifted (Figure 4.1); the curve corresponding to mx= 1.5 is moved towards the hyperpolarized membrane potential. The mechanism underlying bursting in this model is mostly based on the interplay of activation and inactivation of the Ip. This current inactivates during the burst phase and slowly activates during the interburst interval (Butera et al., 1999; Rybak et al., 2003). The steady-state conductance of Ip has its maximum value, 0.32 nS, at the potential of -0.046 V. To tune the activity of the CI model satisfying the characteristics obtained experimentally, parameters of the leak current were tuned, and gh was set to 2 nS. Bursting activity of the CI model of period 2.17 sec is shown in Figure 4.2; parameters were fixed to gleak=1.2 nS, and Eleak= 0.0635 V. The blockade of Ip terminated the bursting activity of the model and led to silence. A blockade of Ih modulated the bursting activity in the following way: if the model exhibited irregular bursting, the blockade of Ih decreased the period of activity and bursting became more regular. This effect is similar to the experimental observation (Thoby-Brisson et al., 2000) and it

72 can be demonstrated with parameters gleak = 1.1 nS and wx=1.5. If the CI model produced regular bursting activity the decrease of Ih increased the period of activity. Figure 4.1: The steady-state curves of the activation gating variable (black) and steady-state conductance (blue) of Ih of the Cd2+ -insensitive pacemaker. The dashed grey curve corresponds to the activation curve mh shifted towards hyperpolarized membrane potential, mx=1.5. Figure 4.2: Bursting activity exhibited by the model of Cd2+ -insensitive pacemaker. The period of the bursting activity is 2.17 s. The leak current parameters were gleak =1.2 nS, Eleak=-0.0635 V.

73 2.4.2. Bistability of Bursting and Silence of the Cd2+ -Insensitive Pacemaker with Ih Modulation Following the techniques described in the previous chapters we computed a one dimensional bifurcation diagram for stationary and oscillatory states using gleak as the bifurcation parameter. First, we set the gleak value so that the model exhibited the hyperpolarized rest state. By successively decreasing this parameter, we found that the stable rest state lost its stability at the super-critical Andronov-Hopf bifurcation (AHB) characterized by the negative sign of the first Lyapunov coefficient. The AHB gave rise to stable sub-threshold oscillations (Figure 4.3). Starting from the gleak value corresponding to AHB, we traced the evolution of the stable periodic orbit as we decreased gleak. At the critical value gleak = 2.0285 nS the stable orbit underwent a saddle-node bifurcation for periodic orbits, LPc (Figure 4.3). At this bifurcation value the stable and unstable orbits coalesce and disappear. After we detected this bifurcation, we were able to trace back the unstable orbit. As we increased gleak, the unstable periodic orbit rapidly grew in size and then terminated at the homoclinic bifurcation gleak = 2.02851 nS, HB, located in close vicinity to LPc (Figure 4.3). We found that the bursting could co-exist with the subthreshold oscillations in the tiny range between the critical values for gleak of LPc and HB bifurcations. The range of gleak supporting stable oscillations is limited by two bifurcations: the saddle-node bifurcation for periodic orbit (LPc) and the super-critical Andronov-Hopf bifurcation (Figure 4.3). The transition from the stable oscillations to silence appeared as a smooth transition as gleak crossed the value of 2.02856 nS. The range of gleak parameters from the interval [1.130 2.0285] nS supports the bursting activity (Figure 4.3). The initial CI pacemaker model did not exhibit bistability of bursting and silence. It left the open question of whether this CI pacemaker model could exhibit bistability of bursting and silence under a different parameter regime.

74 Since Ih was commonly found to be a target of modulation, we explored whether a change of the kinetics of this current could give the model a propensity for multistability of bursting and silence. We took into consideration the knowledge from Chapter III that and Ih current increases the propensity index of a neuron to bistability. The dynamics of the Ih current were shown to be controlled by intracellular cAMP or Ca2+ concentration and neuromodulators such as serotonin. For example, the change of the intracellular cAMP in neurons shifts the activation curve and affects its slope (Robinson and Siegelbaum, 2003; Lüthi and McCormick, 1999). We modulated the steady-state activation curve of mh as it is shown in Figure 4.1 (the dashed grey curve) by setting mx = 1.5. We repeated our analysis (Figure 4.4) and constructed the one-dimensional diagram using gleak as a controlling parameter. The diagram showed that the hyperpolarized rest state lost its stability at the Andronov-Hopf (AHB) bifurcation. This time it was the sub-critical bifurcation and gave rise to the saddle periodic orbit. As gleak was increased to back-trace the orbit, the saddle orbit approached the saddle equilibrium and terminated at the homoclinic bifurcation (HB). On the other hand, the range of gleak supporting the bursting activity was determined between the values 1.099 nS and 2.005 nS (Figure 4.4). This way, we detected the bistability of bursting and silence exhibited by the model within the range of gleak from 2.00503 nS to 2.00511 nS. This range of bistability is bounded by two bifurcations: the sub-critical Andronov-Hopf and homoclinic. Thus, we found the same mechanism of bistability of bursting and silence that was described in the simplified and canonical models of the leech heart interneuron.

75 Figure 4.3: The gleak-parameter bifurcation diagram of the oscillatory and stationary states of the Cd2+ - insensitive pacemaker. The solid green line identifies the hyperpolarized stationary state. This rest state represents the silent regime. The rest state loses stability at the Andronov-Hopf bifurcation, AHB, the value of gleak = 2.02856 nS marked by the dashed vertical line. The AHB gives rise to the stable periodic orbit. The orbit represents the stable subthreshold oscillations. The dashed blue curve depicts the saddle stationary state; the depolarized branch of the saddle stationary states is not shown. The stable periodic orbit terminates at the saddle-node bifurcation for periodic orbits (LPc). The unstable periodic orbit appears at LPc corresponding to the value of gleak 2.02851 nS and disappears at the homoclinic bifurcation, marked by the vertical line. The blue rectangle indicates the partial range of gleak corresponding to the bursting activity. The yellow and white rectangles indicate the range of stable subthreshold oscillations and the silent regime.

76 Figure 4.4: The gleak-parameter bifurcation diagram of the oscillatory and stationary states of the Cd2+ - insensitive pacemaker when the activation of current Ih is changed (mx=1.5). The green curve represents the rest state. The Andronov-Hopf bifurcation (AHB) occurs at gleak= 2.00503 nS (marked by the vertical dashed line), where the rest state loses stability and the saddle periodic orbit appears. The blue dashed curve depicts the saddle equilibria. The saddle orbit that corresponds to the unstable subthreshold oscillations is shown by two red curves; they locate the minimum and maximum values of membrane potential of the saddle orbit. The saddle orbit terminates at the homoclinic bifurcation, HB, gleak= 2.005113 nS. The blue rectangle indicates the part of the range where bursting activity is observed. The white rectangle indicates the range of the silent regime. To confirm with numerical experiments that the model has an observable bistability of bursting and silence, we picked a value for gleak from the range between 2.00503 nS and 2.00511 nS, which located AHB and HB, correspondingly. We were able to elicit a switch from bursting to silence by brief hyperpolarizing pulses of current. An example of perturbation by a negative

77 pulse of current with the amplitude 0.2 pS and the duration 0.01 sec is shown in Figure 4.5A. On the other hand, the switch from silence to bursting activity could be also toggled by a perturbation by a pulse of current. For this switch pulses of either polarity, negative or positive, could be used (Figure 4.5 B-C). Hence, bistability of bursting and silence was demonstrated in the model of the CI pacemaker when the dynamics of the Ih current were modified. Figure 4.5: Switches between bursting activity and silence in the model of the Cd2+ - insensitive pacemaker. A) A hyperpolarizing pulse of current Iinj switches the activity from bursting to silence. A negative pulse with the amplitude of 0.2 pA and the duration of 0.01 s was applied. B) The same pulse as in (A) switches the activity from silence to bursting. C) A positive, depolarizing pulse of current with the same amplitude as in (A) and (B) initiates the bursting activity. The parameters of the leak current were gleak = 2.00511 nS, Eleak =-0.0635 V, and mx = 1.5.

78 2.4.3 Model of the Cd2+ - Sensitive Pacemaker We designed this model to mimic the dynamics of the Cd2+-sensitive pacemaker. We tuned the CS model to satisfy the following constraints. The period of bursting activity was 1.5 sec - 4.5 sec (Peña et al., 2004). The blockade of the persistent Na+ current does not annihilate bursting activity of the CS pacemaker although it slightly increases the period. The blockade of Ca2+ currents makes these neurons silent. The blockade of hyperpolarization-activated non- selective cationic current (Ih) does not annihilate bursting activity (Thoby-Brisson et al., 2000). We built the CS pacemaker model based on the CI pacemaker model by incorporating Ca2+ current into it. We used the kinetics of Ca2+ currents measured by Elsen and Ramirez (1998): high-voltage activated (ICaH) and low-voltage activated (ICaL). The half-maximum activation of IK was shifted from the original value of -0.0445 V to -0.050 V. This model contained six voltage-gated ionic currents: INaf, INap, IK, ICaL, ICaH, and Ih. They are described by a system of 11 differential equations:

79 I I CaL Naf I Nap I K dV C [[ g Naf m Naf hNaf [V E Na ] g Napm NaphNap[V E Na ] [ g K mK [V E K ] g CaLmCaLhCaL[V ECa ] 3 4 2 dt I CaH I h I leak g CaH mCaH hCaH [V ECa ] [ g h mh [V E h ] g leak [V Eleak ] I inj ], 2 2 dm Naf f (0.0438,0.006,V ) m Naf , dt (0.0009,0.0438,0.014,V ) dhNa f (0.0675,0.0108,V ) hNaf , dt (0.0352,0.0675,0.0128,V ) dm Nap f (0.0471,0.0031,V ) m Nap , dt (0.001,0.0471,0.0062,V ) dhNap f (0.059,0.006,V ) hNap , dt (7,0.059,0.006,V ) dmK f (0.05,0.005,V ) mK , dt m Nap (0.004,0.0445,0.01,V ) dmh f h (V ) mh , dt m h (V ) dmCaL f (0.05905,0.00236,V ) mCaL , dt 0.2 dhCaL f (0.08072,0.00531,V ) hCaL , dt 1.0 dmCaL f (0.02782,0.00569,V ) mCaH , dt 0.05 dhCaH f (0.05221,0.00523,V ) hCaH dt 0.5 where C is 0.0362 nF, parameters of the conductances and reversal potentials of ionic currents are gNaf= 150 nS, gNap= 4 nS, gK= 10 nS, gCaL= 15 nS, gCaH= 20 nS, gh= 4 nS, gleak= 2.28 nS, ENa= 0.060 V, EK= -0.080 V, ECa= 0.060 V, Eh= -0.0414 V, Eleak= -0.0634 V, and mx=1. The steady-state activation and inactivation functions and time constants of gating variables were described previously in the CI pacemaker model.

80 To achieve the bursting activity of the CS model with a period close to the experimentally recorded value (Peña et al., 2004), the conductance and reversal potential of the leak current were tuned, gleak = 2.28 nS and Eleak = -0.0634 V; and the model exhibited the bursting activity with a period of 4.17 s (Figure 4.6 ). The blockade of the Ip current did not terminate bursting, yet slightly increased the period of activity, which was in accordance with the experimental constraints. The reduction of gh slightly increased the period of the bursting activity. The complete removal of Ih (gh=0) did not terminate the bursting activity. The bursting activity was driven by Ca2+ currents; if we set gCaLto zero, the model exhibited hyperpolarized silence. Figure 4.6: Bursting activity of the model of the Cd2+ - sensitive pacemaker. The period of the bursting activity is 4.17 s. The leak current parameters were gleak=2.28 nS and Eleak=-0.0643 V. We investigated how variation of gleak affected the oscillatory activity of the CS pacemaker model (Figure 4.6). The model with gleak chosen within the interval [1.95 3.478] nS exhibited bursting activity. At the critical value of gleak = 1.95 nS the model underwent the transition from bursting to tonic spiking, whereas the value of gleak= 3.487 defined the transition from bursting to silence. We showed in Chapters I and II that near transition from bursting to

81 silence the models could exhibit bistability of these two regimes. To investigate whether the CS model could be bistable, we again constructed a one-dimensional bifurcation diagram using gleak as the bifurcation parameter. Following our procedure, we detected that at the critical value gleak =3.4782 nS the sub-critical Andronov-Hopf bifurcation occurred characterized by the positive sign of the first Lyapunov coefficient. At this bifurcation the first saddle orbit emerged. As the parameter gleak was successively increased from the AHB bifurcation, the first saddle orbit was back traced until it underwent the saddle-node bifurcation (LPc) at the parameter gleak = 3.4808 nS. At LPc, two saddle orbits merged and disappeared. The bifurcation LPc happened in close proximity to the Andronov-Hopf bifurcation (Figure 4.7). With gleak decreased further from the LPC bifurcation, the second saddle periodic orbit traced until it disappeared at the homoclinic bifurcation, HB (Figure 4.7). For values of gleak close to the parameter value for AHB the period of bursting activity was very long. For example, at gleak = 3.478 nS the period of bursting activity reached 3500 sec. The bistability of bursting and silence was not detected in the initial CS pacemaker model.

82 Figure 4.7: The gleak-parameter bifurcation diagram of oscillatory and stationary states of the Cd 2+- sensitive pacemaker. The diagram is obtained for Eleak = −0.0634 V. The solid green and dashed blue curves represent the membrane potential of the equilibria, stable and unstable, correspondingly. The saddle periodic orbit emerged at the sub-critical Andronov-Hopf bifurcation, AHB, where the rest state loses its stability (dashed). The dashed red curves depict the minimum and maximum membrane potential values of the saddle periodic orbits, first and second. The first saddle orbit emerged at the Andronov-Hopf bifurcation and disappeared at the saddle-node bifurcation for periodic orbit, marked by the two blue dots LPc. The second saddle orbit appeared at LPc, and terminated at the homoclinic bifurcation, HB. The blue and the white rectangles indicate the ranges of the bursting activity and silent regime, respectively.

83 2.4.4 Bistability of Bursting and Silence of the Cd2+ - Sensitive Pacemaker Induced by Blockade of Ip The analysis of the first CS model’s dynamics did not reveal the bistability of bursting and silence. We tested whether it would produce bistability under conditions that are pharmacologically achievable. We investigated the effect of the blockade of the persistent Na+ current. The persistent Na+ current can be blocked by riluzole and under this condition the bursting activity of the CS neuron would not terminate. To reproduce the blockade in the model, the gp parameter was set to 0 nS. We constructed the gleak-parameter bifurcation diagram (Figure 4.8). The analysis showed that the model exhibited bursting activity for the range of gleak parameters located within the interval [1.882 2.293] nS. The bistability of bursting and silence was detected in the interval of gleak [2.270 2.293] nS. The area of bistability is limited by the Andronov-Hopf and homoclinic bifurcations. The mechanism supporting bistability is based on the saddle periodic orbit, which is the same as the one described in the simplified 4D and 14D models of the leech heart interneuron (Malashchenko, Master Thesis, 2007). We also examined the activity of the model in response to perturbations by either a negative or positive pulse of current. For example, a negative pulse of current with the amplitude 0.5 pA and the duration 2 sec switched the activity from bursting to silence (Figure 4.9 A). The switch of the activity from a silent regime to bursting could be triggered by the negative or positive pulses of currents, as it is shown in Figure 4.9 B C.

84 Figure 4.8: The gleak-parameter bifurcation diagram of oscillatory and stationary states of the Cd2+ - sensitive pacemaker, with Ip blocked. The solid green and dashed blue curves represent the membrane potential of the equilibria, stable and unstable, correspondingly. The lower branch corresponds to the hyperpolarized stationary state. The saddle periodic orbit emerges at the sub-critical Andronov-Hopf bifurcation, AHB, that turns the hyperpolarized branch unstable (dashed). Dashed red curves depict the minimum and maximum membrane potential values of the saddle periodic orbit. The dashed brown curve between them corresponds to the average membrane potential of the saddle orbit. The saddle orbit terminates at the homoclinic bifurcation, HB. The blue and white rectangles indicate the range of the bursting activity and silent regime, correspondingly.

85 Figure 4.9: Bistability of bursting activity and silence in the model of the Cd2+ -sensitive pacemaker with gp = 0 nS. The parameters of the leak current are gleak = 2.289 nS and Eleak= -0.0634 V. The amplitude of the pulse was 0.5 pA for all traces. The duration of the pulse was 2 s. A) A hyperpolarizing pulse of current Iinj was applied to switch the activity from bursting to silence. Hyperpolarizing (B) and depolarizing (C) pulses of the current switched the activity from silence to bursting.

86 2.4.5 Conclusion IV In this chapter we developed models of two types of inspiratory neurons using as a basis the model of a respiratory neuron developed by Rybak et al. (2003), Cd2+ -sensitive and Cd2+ - insensitive. By applying bifurcation analysis we found that both types of models can exhibit the bistability of bursting and silence under certain modulations of ionic currents. A shift of activation of Ih towards hyperpolarized membrane potential induced the bistability of bursting and silence in the Cd2+ -insensitive pacemaker. The blockade of persistent Na+ current caused bistability in the model of the Cd2+ -sensitive pacemaker. The mechanism supporting bistability in these models is in accordance with the mechanism of bistability in the simplified and canonical models of the leech heart interneuron (Malashchenko, Master Thesis, 2007). 3. DISCUSSION In a single neuron, interplay between different ionic currents can lead to the coexistence of different regimes of activity, such as tonic spiking, bursting, silence, and subthreshold oscillations. The documented types of bistability include the coexistence of bursting and tonic spiking (Lechner et al., 1996; Canavier et al., 1994; Cymbalyuk et al., 2002; Fröhlich and Bazhenov, 2006), of depolarized and hyperpolarized silent states (Hounsgaard and Kiehn, 1989; Crunelli et al., 2005), of two tonic spiking regimes (Cymbalyuk and Shilnikov, 2005), of tonic spiking and silence (Guttman et al., 1980), and of different patterns of bursting (Canavier et al., 1994; Cymbalyuk et al., 2002; Butera, 1998). Although bistability as a phenomenon has been shown in many different neuronal systems, the biophysical mechanisms supporting bistability are

87 largely not known. Studying neuronal dynamics on the cellular level allows one to investigate the mechanisms supporting bistability which have been the focus of numerous experimental and theoretical studies. One of the most studied types of bistability is the coexistence of tonic spiking and silence (Rinzel, 1978; Guttman et al., 1980). In contrast, the coexistence of bursting and silence has not previously been the focus of any theoretical or experimental study of the dynamics of a single neuron. Here we filled this gap in part, by describing multiple scenarios of such coexistence in the dynamics of neuronal models. The subject of our study was the dynamics of endogenously bursting neurons. We considered neurons controlling the leech heart beat and neurons controlling respiratory rhythm in mice and rats. First, we developed a low-dimensional model of the leech heart interneuron and we detected that the model is capable of producing six types of multistability. We further elaborated the study of neuronal bistability in a high-dimensional neuronal model of the leech heart interneuron that is based on the dynamics of nine voltage-dependent currents. Then, we investigated the effect of current modulation on a propensity to bistability. Such analysis allowed us to predict conditions causing a neuronal propensity to bistability of bursting and silence in the leech heart interneuron. We also analyzed the dynamics of two models of respiratory pacemakers, the Cd2+ - sensitive and the Cd2+ - insensitive. We showed that both types of neurons could exhibit the bistability of bursting and silence under modulation of certain ionic currents. The co-existence of tonic spiking and silence in particular has been extensively studied (Rinzel, 1978). It was detected in various types of neurons, from spinal motoneurons to neuronal networks of the entorhinal cortex, and the thalamocortical network (Egorov et al., 2002; Fuentealba et al., 2005; Hounsgaard and Kiehn, 1989; Williams et al., 2002). Bistability has been

88 found to occur spontaneously or to be induced by neuromodulators. For example, spontaneous bistability was demonstrated in Purkinje cells in the rat and guinea pig (Loewenstein et al., 2005). Turtle and cat motoneurons become bistable after modulation of membrane properties by serotonin (Hounsgaard and Kiehn, 1989; Perrier and Hounsgaard, 2000). Analysis of neuronal Hodgkin-Huxley type models exhibiting coexistence of tonic spiking and silence suggests two main types of mechanisms. First, the depolarized tonic spiking regime and rest state are separated by the stable manifold of a saddle equilibrium (Izhikevich, 2007). For example, this mechanism supports bistability of tonic spiking and silence observed in the simplified model of cerebellar Purkinjie neurons (Loewenstein et al., 2005; Fernandez et al., 2007). The second type, the trajectory of the tonic spiking of larger amplitude goes around the hyperpolarized rest state so that the states are separated by the stable manifold of a saddle periodic orbit (Guttman et al., 1980; Gutkin et al., 2009; Hahn and Durand, 2001; Rinzel, 1978). The regimes of activity that are separated in phase space by the stable manifold of the saddle type regime, either periodic orbit or equilibrium state, have also been shown in previous studies by Cymbalyuk and Shilnikov (2005) in a 3D reduced model of the leech heart interneuron. The presence of a saddle type separating regime was detected in four types of bistability under different parameter regimes: (1) tonic spiking and bursting; (2) tonic spiking and tonic spiking; (3) tonic spiking and hyperpolarized silence; (4) bursting and depolarized silence (Cymbalyuk and Shilnikov, 2005; Shilnikov et al., 2005; Shilnikov et al., 2008). A saddle periodic orbit is a generic component of a number of mechanisms supporting bistability in the dynamics of single neurons (Guttman et al., 1980; Hahn and Durand, 2001; Rinzel, 1978; Fröhlich and Bazhenov, 2006; Cymbalyuk and Shilnikov, 2005; Cressman et al., 2009; Shilnikov et al., 2005). Its stable manifold acts as a threshold between the two regimes; and different perturbations of the state of the model could

89 induce crossing the threshold and a switch from one regime to the other. Analysis of the process of switching between regimes of activity with noise or pulses of injected current applied to a bistable neuron or network is a valuable tool for the investigation of bistability (Guttman et al., 1980; Gutkin et al., 2009; Hahn and Durand, 2001; Paydafar et al., 2006). After the separating regime is identified, an essential part of the description of a mechanism supporting bistability is the description of the appearance and disappearance of the regime under variation of a controlling parameter. This description allows one to identify the factors which cause the bistability and the range of the parameter values supporting it. The choice of parameter is usually dictated by the topic of study. Bifurcation theory provides appropriate tools for this task. For the Hodgkin-Huxley model, the Rinzel mechanism determines that the unstable orbit appears through the Andronov- Hopf bifurcation and disappears through the saddle-node bifurcation for periodic orbits as the polarizing current is varied (Rinzel, 1978). The Hodgkin-Huxley model was also analyzed under variation of the concentration of K+. This analysis was motivated by the hypothesis that the elevation of external K+ concentration can induce seizures. The analysis shows the same mechanism of bistability (Hahn and Durand, 2001). We investigated the mechanisms supporting multistability of the leech heart interneuron. Due to the complexity of the model's dynamics, we simplified the high-dimensional model. As a result, a low dimensional model was developed, which exhibits temporal characteristics close to those recorded from leech heart interneurons. Although the model is simple, it maintains a number of biophysical correspondences to the living counterpart. With this model we explored and accentuated the role of the low threshold slowly inactivating calcium current, as one sufficient to support a plethora of different mechanisms supporting the coexistence of different

90 regimes of activity. It shows six types of multistability: (1) tonic spiking and the hyperpolarized silent regime; (2) tonic spiking and subthreshold oscillations; (3) tonic spiking and bursting; (4) bursting and subthreshold oscillations; (5) bursting and the hyperpolarized silent regime; and (6) bursting, subthreshold oscillations, and silence. We illustrated some of these mechanisms by a series of numerical experiments to show that they are sufficiently robust to be observed. In this study we used the leak conductance as the controlling parameter. It was singled out as such due to the following rationale: in contrast to the other conductances, the leak conductance does not depend on membrane potential, yet can be modulated. Its variation has a strong effect on the oscillatory properties of leech heart interneurons (Cymbalyuk et al., 2002). By using slow-fast decomposition techniques, the universality of leak current parameters has been exploited in analysis of neuronal dynamics in a variety of models. Typically, low dimensional systems that are based on Hodgkin-Huxley formalism can be dissected in different time scales, slow and fast (Rinzel, 1985; Bertram et al., 1995). Such decomposition allows considering the slow variable as constant in the fast time scale and using it in turn as a controlling parameter. It has been shown in neuronal models containing slow currents that the leak current could be used as a tool to classify the dynamic regimes of a model, since slow currents could be approximately considered as constant in a fast time scale, and the slow currents can be lumped together as a component of the leak current (Guckenheimer et al., 2005). Guckenheimer et al. (2005) applied this idea to describe the neuronal dynamics in several models for which the Na+ current is the dominant fast current. The authors considered models of the Aplisia neuron, the thalamic relay neuron, a simplified model of the leech heart interneuron, and a neuron from the respiratory center. They detected bifurcations of a fast subsystem that describe the transitions

91 between spiking and silent phases. The analysis was presented in forms of 2D bifurcation diagrams depicting the regions of oscillatory and stationary states. The leak current plays an important physiological role in controlling the excitability of the cell. For example, it has been shown that neuronal activity of leech heart interneurons is highly sensitive to the leak current. The notion of such sensitivity is supported by experiments showing that the activity of the leech heart interneuron depends on the method of recording. The pharmacologically isolated heart interneuron produces tonic spiking activity if the recording is made intracellularly, whereas it exhibits bursting activity when the recording is performed extracellularly (Cymbalyuk et al., 2002). It has been suggested that the damaged interneuron introduces an additional inward leak current that provides a shunting component to the leak current. Similarly, sensitivity to the method of recording has also been reported in Xenopus neurons (Aiken et al., 2003). The intracellular recording shows that Xenopus neurons fire a single short burst or single spikes in response to a depolarizing pulse of current, whereas the whole-cell recording shows tonic spiking activity in response to stimulation. In this work, the leak current is described as a non-selective cationic current primarily driven by K+ ions. We used leak current parameters to describe five different mechanisms of multistability detected in a low-dimensional model of the leech heart interneuron. We detected separating regimes supporting different types of multistability as leak conductance varied. Then, we described their evolution in terms of bifurcations that occurred in the system. These bifurcations limit the range of leak current parameters where different types of multistability can be observed. We focused on bistability of bursting and silence since this type had not been detected in isolated neurons experimentally. We thoroughly studied the mechanism in the low- dimensional model as well as in the 14D model of the leech heart interneuron (Malashchenko,

92 Master Thesis, 2007). Similarly to the Rinzel mechanism, the mechanism supporting bistability of bursting and silence requires the existence of an unstable periodic orbit allowing separation of the stable rest state (silence) and bursting activity. By applying the bifurcation analysis to models of the leech heart interneuron, we showed that for a given value of E leak the stable hyperpolarized equilibrium and bursting regime co-existed in a certain range of values of gleak. For a smaller gleak, this range is limited by the sub-critical Andronov-Hopf bifurcation at which the unstable periodic orbit appears. In this regard, the orbit appears in accordance with the Rinzel mechanism of the coexistence of tonic spiking and silence. For a larger gleak, the range of bistability is limited by the homoclinic bifurcation of a periodic orbit. The latter distinguishes the mechanism presented here from the Rinzel mechanism, which has the unstable orbit ceasing at the saddle- node bifurcation for periodic orbits. Outside of the above-mentioned range, bistability of bursting and silence was not observed. We performed a two-parameter bifurcation analysis to trace the Andronov-Hopf and homoclinic bifurcations, using the leak current parameters as the bifurcation parameters. These curves were placed together with the numerically integrated borders of bistability on the bifurcation diagram in the (gleak Eleak) parameter space. On the diagrams for the low dimensional model and for the 14D model, we showed that the transition from bursting to silence coincides with the two-parameter curve of the homoclinic bifurcation (Malashchenko, Master Thesis, 2007). This type of analysis of biophysically realistic models could provide help establishing the origin and roles of bistability in the functioning of the nervous system. Similarly to electronic and other technical devices, in neuronal networks the bistable neurons could play the roles of switch, relay, logic, and memory elements (Marder et al., 1996). In the realm of motor control, bistable neurons appear to be a natural choice as elements of multifunctional central pattern

93 generators so that the same network can generate several behavioral patterns (Briggman and Kristan, 2008; Venugopal et al., 2007; Shilnikov et al., 2008). For example, Brigman and Kristan (2008) suggest that crawling and swimming CPGs of the leech could be an example of a multifunctional network. The authors show that the vast majority of neurons are involved in both behaviors when swimming or crawling is produced. In contrast to multifunctional central pattern generators, for the leech heartbeat central pattern generator bistability in the leech heart interneurons seems to be a pathological, life-threatening phenomenon. This study predicts the existence of bistability of the leech heart interneuron under conditions leading to elevated conductance of the leak current. Here, we used the range of leak current parameters leading to bistability as a marker characterizing a neuronal propensity for bistability. In Figure 4 the range of bistability of bursting and silence was detected as relatively narrow, bringing difficulties to detecting this type of bistability experimentally. Therefore, the next question we addressed here is which biological conditions can affect this range. In this study we introduced a new measure characterizing the bistability of bursting and silence, called the propensity index. It was defined as the range of leak conductance leading to bistability. Since the neuron is isolated, we analyzed how up- or down- regulations of ionic currents affect the propensity index. Regulation of ionic current can be controlled by neuromodulators. Moreover, in some neurons it has been shown that neuromodulators could lead to multistability. For example, serotonin induces bistability of spinal motoneurons in different animals. Hounsgaard and Kiehn (1989) found bistability of tonic spiking and silent regimes in the presence of serotonin in an in vitro preparation of the spinal cord of the turtle. They showed that serotonin supports the depolarized silent regime, plateau potential, and the Na+-dependent tonic spiking activity.

94 Specifically, the plateau potential is maintained by an L-type calcium current and a Ca2+- dependent mechanism. It can be revealed by the blockade of Na+ current with TTX. The perturbation of the neuron by a brief pulse of current switches the activity between the hyperpolarized and depolarized states (Perrier, Majia-Gervacio, Hounsgaard, 2000; Perrier and Hounsgaard, 2000). In the leech heart interneuron two L-type Ca2+ currents are described, slow and fast. Our analysis of the model indicated that these currents had different influences on the propensity to bistability. Data analysis showed that the increase of slow Ca2+ conductance increased the propensity to bistability of bursting activity and silence, whereas the increase of fast Ca2+ conductance reduced it. Neuromodulators play an important role in controlling motor behavior. In the leech, endogenous peptide myomodulin has been shown as a modulator of bursting activity of heart interneurons. The myomodulin decreases the period of interneuron activity through the enhancement of the conductance of h-current (Tobin and Calabrese, 2005). Based on our prediction that h-current leads to bistability, the question of whether myomodulin could induce bistability in an interneuron naturally arises. In other biological systems the role of h-current in bistability is alternative depending on the cell type and on the mechanism supporting bistability. An example is the bistability of tonic spiking and silence of cerebellar Purkinje neurons observed in vitro. Some studies show the blockade of this current induces bistability (Williams et al., 2002) whereas others reported bistability of Purkinjie neurons in the presence of h-current (Loewenstein et al., 2005). The ionic mechanism of Purkinjie cell bistability is suggested to rely on non-inactivating inward currents such as the persistent Na+ current that initiates the depolarized tonic spiking regime (Williams et al., 2002). Our data indicated that persistent Na+ current plays an insignificant role in bistability of bursting and silence of the leech heart

95 interneuron. We also showed that in the model of the Cd2+ - sensitive inspiratory pacemaker, the blockade of the persistent Na+ current induced the bistability of bursting and silence. Whereas in the Cd2+ - insensitive pacemaker model of the respiratory center it was the change of kinetics of h-current that caused the bistability. These findings might shed light on differences in the roles of h-current and persistent Na+ current in bistability of Purkinjie cell and the neurons described here. We investigated the effect of up- or down-regulation of each current on the bounds limiting the separating regime and the propensity index. Only two currents were shown to increase significantly the propensity index: fast Ca2+ current and hyperpolarization-activated cationic current. Based on our modeling study, we predicted that by manipulating these currents experimentally one could reveal the bistability of bursting and silence in a single neuron experimentally. The methods described in this work allow one to investigate the effect of currents in bistability and can be applied to other neuronal models if the mechanism and the separating regime supporting bistability are known.

96 Publications: 1. Malashchenko, T., Shilnikov, A., & Cymbalyuk, G. Coexistence of oscillatory and silent regimes in a model of a neuron. Submitted to PloS One. 2. Malashchenko, T., Shilnikov, A., & Cymbalyuk, G. Bistability of bursting and silence regimes in a model of a leech heart interneuron. Submitted to Phys Rev E. 3. Malashchenko, T., & Cymbalyuk, G. Propensity to bistability of bursting and silence in a model of a leech heart interneuron. In preparation. Conference Presentations: 1. [May 2010] Oral presentation at Brain and Behavior retreat. Malaschenko, T., & Cymbalyuk, G. Good currents playing bad: propensity to bi-stability in neuronal dynamics. Atlanta, GA. 2. [Apr 2010] Poster presentation at Workshop “Dynamics of Bursting Activity of Neurons”. Malaschenko, T., & Cymbalyuk, G.. Good currents playing bad: propensity to bi-stability in neuronal dynamic. Atlanta, GA. 3. [Nov 2009] Oral presentation at SESAPS meeting. Malashchenko, T., & Cymbalyuk, G. The propensity to bi-stability of bursting and silence of the leech heart interneuron. Atlanta, GA. 4. [Oct 2009] Poster presentation at SfN meeting. Malashchenko, T., Barnett, W. H., & Cymbalyuk, G. Co-existence of bursting and silence in the dynamics of a leech heart interneuron. Chicago, IL.

97 5. [Jul 2009] Poster presentation at CNS meeting. Malashchenko, T., Barnett, W. H., Burylko, O., & Cymbalyuk, G. Control of bursting activity by modulation of ionic currents. Berlin, Germany. 6. [Apr 2009] Poster presentation at SENN conference. Malaschenko T., & Cymbalyuk, G. The propensity to bi-stability of the leech heart interneuron model. Jacksonville, FL. 7. [Nov 2008] Poster presentation at SfN meeting. Barnett, W. H., Malaschenko, T., & Cymbalyuk, G. Role of slow currents in regulation of bursting activity. Washington, DC. REFERENCES: 1. Aiken, S. P., Kuenzi, F. M., & Dale, N. (2003) Xenopus embryonic spinal neurons recorded in situ with patch-clamp electrodes--conditional oscillators after all? European Journal of Neuroscience, 18(2), 333-343. 2. Bayliss, D. A., Sirois, J. E., &Talley, E. M. (2003). The TASK family: two-pore domain background K+ channels. Molecular Interventions, 3(4), 205-219. 3. Butera, R. J. (1998). Multirhythmic bursting. Chaos, 8(1), 274-284. 4. Berkowitz, A. (2008). Physiology and morphology of shared and specialized spinal interneurons for locomotion and scratching. Journal of Neurophysiology, 99, 2887-2901. 5. Bertram, R., Butte, M. J., Kiemel, T., & Sherman, A. (1995). Topological and phenomenological classification of bursting oscillations. Bulletin of Mathematical Biology, 57(3), 413–439.

98 6. Briggman, K. L., & Kristan, W. B. (2008). Multifunctional pattern-generating circuits. Annual Review of Neuroscience, 31, 271-294. 7. Brunel, N., & Van Rossum, M. C. (2007). Quantitative investigation of electrical nerve excitation treated as polarization. Biological Cybernetics, 97 (5-6), 341–349. 8. Calabrese, R. L., Nadim, F., & Olsen, O. H. (1995). Heartbeat control in the medicinal leech: a model system for understanding the origin, coordination, and modulation of rhythmic motor patterns. Journal of Neurobiology, 27(3), 390-402. 9. Canavier, C. C., Baxter, D. A., Clark, J.W., & Byrne, J. H. (1994). Multiple modes of activity in a model neuron suggest a novel mechanism for the effects of neuromodulators. Journal of Neurophysiology, 72, 872-882. 10. Channell, P., Cymbalyuk, G., & Shilnikov, A. (2007). Origin of bursting through homoclinic spike adding in a neuron model. Phys Review Letters, 98 (13), 134101. 11. Crunelli, V., Toth, T., Cope D., Blethyn K., & Hughes, S. (2005) The 'window' T-type calcium current in brain dynamics of different behavioural states. Journal of Physiology, 562(1), 121-129. 12. Cymbalyuk, G. S., & Calabrese, R. L. (2001). A model of slow plateau-like oscillations based upon the fast Na+ current in a window mode. Neurocomputing, 38-40, 159-166. 13. Cymbalyuk, G. S., Gaudry, Q., Masino, M. A., & Calabrese, R. L. (2002). Bursting in leech heart interneurons: cell autonomous and network based mechanisms. Journal of Neuroscience, 22, 10580-10592. 14. Cymbaluyk, G., & Shilnikov, A. L. (2005). Coexistence of tonic spiking oscillations in a leech neuron model. Journal of Computational Neuroscience, 18(3), 255-263.

99 15. Doi, A., & Ramirez, J. M. (2008). Neuromodulation and the orchestration of the respiratory rhythm. Respiratory Physiology & Neurobiology, 164, 96-104. 16. Durand, D., Jensen, A., & Bikson, M. (2006). Suppression of neural activity with high frequency stimulation. Conf. Proc. IEEE Engineering in Medicine and Biology Society, 1, 1624-1625. 17. Durstewitz, D., & Seamans, J. K. (2006). Beyond bistability: biophysics and temporal dynamics of working memory. Neuroscience, 139(1), 119-133. 18. Egorov, A., Hamam, B., Fransén, E., Hasselmo, M., & Alonso, A. (2002). Graded persistent activity in entorhinal cortex neurons. Nature, 420(6912), 173-178. 19. Elsen, F. P, & Ramirez, J. M. (1998). Calcium currents of rhythmic neurons recorded in the isolated respiratory network of neonatal mice. Journal of Neuroscience, 18(24), 10652-10662. 20. Fall, C. P., & Rinzel, J. (2006). An intracellular Ca2+ subsystem as a biologically plausible source of intrinsic conditional bistability in a network model of working memory. Journal of Computational Neuroscience, 20(1), 97-107. 21. Fernandez, F. R., Engbers, J. D., & Turner, R. W. (2007). Firing dynamics of cerebellar purkinje cells. Journal of Neurophysiology, 98(1), 278-294. 22. Foss J. and Milton J., In: Milton J., & Jung, J. (eds) Epilepsy as a dynamic disease. Springer, Berlin Heidelberg 16, 283 (2003). 23. Fröhlich, F., & Bazhenov, M. (2006). Coexistence of tonic firing and bursting in cortical neurons. Physical Review E, 74, 031922.

100 24. Fuentealba, P., Timofeev, I., Bazhenov, M., Sejnowski, T., & Steriade, M. (2005). Membrane bistability in thalamic reticular neurons during spindle oscillations. Journal of Neurophysiology, 93(1), 294-304. 25. Getting, P. A. (1989). Emerging principles governing the operation of neural networks. Annual Review of Neuroscience, 12, 185-204. 26. Guckenheimer, J., Tien, J. H., & Willms, A. (2005). Bifurcations in the fast dynamics of neurons: implications for bursting. In S. Coombes & P. C. Bresloff (Eds.), Bursting: the genesis of rhythm in the nervous system. New Jersey: World Scientific Publishing Co. 27. Gutkin, B. S., Jost, J., & Tuckwell, H. C. (2009). Inhibition of rhythmic neural spiking by noise: the occurrence of a minimum in activity with increasing noise. Naturwissenschaften, 96(9), 1091-1097. 28. Guttman, R., Lewis, S., & Rinzel, J. (1980). Control of repetitive firing in squid axon membrane as a model for a neuroneoscillator. Journal of Physiology, 305, 377-395. 29. Hahn P., & Durand, D. (2001). Bistability dynamics in simulations of neural activity in high-extracellular-potassium conditions. Journal of Computational Neuroscience, 11(1), 5-18. 30. Heyward, P., Ennis, M., Keller, A., & Shipley, M. T. (2001). Membrane Bistability in Olfactory Bulb Mitral Cells. Journal of Neuroscience, 21(14), 5311-5320. 31. Hill, A., Lu, J., Masino, M., Olsen, O., & Calabrese, R. L. (2001). A Model of a Segmental Oscillator in the Leech Heartbeat Neuronal Network. Journal of Computational Neuroscience, 10,281-302.

101 32. Hodgkin, A., & Huxley, A. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. Journal of Physiology, 117(4), 500- 44. 33. Hounsgaard, J., & Kiehn, O. (1989). Serotonin-induced bistability of turtle motoneurones caused by a nifedipine-sensitive calcium plateau potential. Journal of Physiology, 414, 265-282. 34. Hughes, S. W., Cope, D. W., T´oth, T. I., Williams, S. R., & Crunelli, V. (1999). All thalamocortical neurones possess a T-type Ca2+ “window” current that enables the expression of bistability-mediated. Journal Physiology, 517 (3), 805-815. 35. Izhikevich, E. (2007). Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. The MIT press. 36. Jalife, J., & Antzelevitch, C. (1979). Phase resetting and annihilation of pacemaker activity in cardiac tissue. Science, 206(4419), 695-697. 37. Khibnik, A. I., Kuznetsov, Yu. A., Levitin, V. V., & Nikolaev, E. V. (1993). Continuation techniques and interactive software for bifurcation analysis of ODEs and iterated maps. Physica D, 62, 360. 38. Koizumi, H., Smerin, S. E., Yamanishi, T., Moorjani, B. R., Zhang, R., & Smith, J. C. (2010). TASK channels contribute to the K+-dominated leak current regulating respiratory rhythm generation in vitro. Journal of Neuroscience, 30(12), 4273-4284. 39. Kuznetsov, Y. A., Levitin, V. V., & Skovoroda, A. R. (1996). Continuation of stationary solutions to evolution problems in CONTENT. Report AMR9611. Amsterdam: Centrum/Voor Wiskunde en Informatica.

102 40. Larkman, P. M. & Perkins, E. M. (2005). A TASK-like pH- and amine-sensitive 'leak' K+ conductance regulates neonatal rat facial motoneuron excitability in vitro. European Journal of Neuroscience, 21(3), 679-691. 41. Lechner, H., Baxter, D., Clark, C., & Byrne, J. (1996). Bistability and its regulation by serotonin in the endogenously bursting neuron R15 in Aplysia. Journal of Neurophysiology, 75(2), 957-962. 42. Loewenstein, Y., Mahon, S., Chadderton, P., Kitamura, K., Sompolinsky, H., Yarom, Y., & Husser, M. (2005). Bistability of cerebellar Purkinje cells modulated by sensory stimulation. Nature of Neuroscience, 8(2), 202-211. 43. Lu, B., Su, Y., Das, S., Liu, J., Xia, J., & Ren, D. (2007). The neuronal channel NALCN contributes resting sodium permeability and is required for normal respiratory rhythm. Cell, 129(2), 371-383. 44. Lüthi, A., & McCormick, D. A. (1999). Modulation of a pacemaker current through Ca(2+)-induced stimulation of cAMP production. Nature Neuroscience. 2(7): 634-641. 45. Marder, E. (1994). Invertebrate neurobiology. Polymorphic neural networks. Current Biology, 4, 752-754. 46. Marder, E., Abbott, L. F., Turrigiano, G. G., Liu, Z., Golowasch, J. (1996). Memory from the dynamics of intrinsic membrane currents. Proceedings of the National Academy of Sciences, 93(24), 13481-13486. 47. Marder, E., & Calabrese, R. L. (1996). Principles of rhythmic motor pattern generation. Physiology Review, 76, 687-717.

103 48. Maistrenko, Y. L., Lysyansky, B., Hauptmann, C., Burylko, O., & Tass, P. A. (2007). Multistability in the Kuramoto model with synaptic plasticity. Physical Review E, 75(6- 2), 066207. 49. Nadim, F., Manor, Y., Kopell, N., & Marder, E. (1999). Synaptic depression creates a switch that controls the frequency of an oscillatory circuit. Proceedings of the National Academy of Sciences, 96(14), 8206-8211. 50. Olypher, A., Cymbalyuk, G., & Calabrese, R. L. (2006). Hybrid systems analysis of the control of burst duration by low-voltage-activated calcium current in leech heart interneurons. Journal of Neurophysiology, 96(6), 2857-2867. 51. Paydafar, D., Forger D., & Clay J. (2006), Noisy inputs and the induction of on-off switching behavior in a neuronal pacemaker. Journal of Neurophysiology, 96(6), 3338- 3348. 52. Perrier, J. F., & Hounsgaard, J. (2000). Development and regulation of response properties in spinal cord motoneurons. Brain Res Bull, 53(5), 529-535. 53. Perrier, J. F., Mejia-Gervacio, S., & Hounsgaard, J. (2000). Facilitation of plateau potentials in turtle motoneurones by a pathway dependent on calcium and calmodulin. Journal of Physiology, 1, 107-113. 54. Peña, F., Parkis, M. A., Tryba, A. K., & Ramirez, J. M. (2004). Differential contribution of pacemaker properties to the generation of respiratory rhythms during normoxia and hypoxia. Neuron, 43(1), 105-117. 55. Ramirez, J. M., Tryba, A. K., & Peña, F. (2004). Pacemaker neurons and neuronal networks: an integrative view. Current Opinion in Neurobiology, 14(6), 665-674.

104 56. Ramirez, J. M., Folkow, L. P., & Blix A. S. (2007). Hypoxia tolerance in mammals and birds: from the wilderness to the clinic. Annual Review of Physiology, 69, 113-143. 57. Robinson, R. B, & Siegelbaum, S. A. (2003). Hyperpolarization-activated cation currents: from molecules to physiological function. Annual Review of Physiology, 65, 453-80. 58. Rinzel, J. (1978). On repetitive activity in nerve. Federation Proceedings, 37, 2793- 2802. 59. Rinzel, J. (1985). Bursting oscillations in an excitable membrane model. Lecture Notes in Mathematics, 1151, 304–316. 60. Rybak, I. A., Shevtsova, N. A., St-John, W. M., Paton, J. F., & Pierrefiche, O. (2003). Endogenous rhythm generation in the pre-Bötzinger complex and ionic currents: modelling and in vitro studies. European Journal of Neuroscience, 18(2), 239-257. 61. Rybak, I. A., Ptak, K., Shevtsova, N. A., & McCrimmon, D. R. (2003). Sodium currents in neurons from the rostroventrolateral medulla of the rat. Journal of Neurophysiology, 90(3), 1635-1642. 62. Sirois, J. E., Lynch, C., & Bayliss, D. A. (2002). Convergent and reciprocal modulation of a leak K+ current and I(h) by an inhalational anaesthetic and neurotransmitters in rat brainstem motoneurones. Journal of Physiology, 541, 717-729. 63. Schmidt, J., & Calabrese, R. L. (1992). Evidence that acetylcholine is an inhibitory transmitter of heart interneurons in the leech. Journal of Experimental Biology, 171, 329- 347. 64. Shilnikov, L. P., Shilnikov, A., Turaev, D. & Chua, L. (1998). Methods of Qualitative Theory in Nonlinear Dynamics. Part I . World Sci.

105 65. Shilnikov, L. P., Shilnikov, A., Turaev, D. & Chua, L. (2001). Methods of Qualitative Theory in Nonlinear Dynamics. Part II. World Sci. 66. Shilnikov, A. L., and Cymbalyuk, G. (2004). Homoclinic saddle-node orbit bifurcations en a route between tonic spiking and bursting in neuron models, Invited review. Regular & Chaotic Dynamics, 3(9), 281-297. 67. Shilnikov, A., & Cymbalyuk, G. (2005). Transition between tonic spiking and bursting in a neuron model via the blue-sky catastrophe. Physical Review Letter, 94(4), 048101. 68. Shilnikov, A., Calabrese, R., & Cymbalyuk, G. (2005). Mechanism of bi-stability: tonic spiking and bursting in a neuron model. Physical Review E, 71, 056214, 1-9. 69. Shilnikov, A. L., & Kolomiets, M. L. (2008). Methods of the qualitative theory for the Hindmarsh-Rose model: a case study. Tutorial. Journal of Bifurcations and Chaos, 18 (8), 1-27. 70. Shilnikov, A., Gordon, R., & Belykh, I. (2008). Polyrhythmic synchronization in bursting networking motifs. Chaos, 18(3), 037120. 71. Smith, J. C., Butera, R. J., Koshiya, N., Del Negro, C., Wilson, C. G, & Johnson, S. M. (2000). Respiratory rhythm generation in neonatal and adult mammals: the hybrid pacemaker-network model. Respiratory Physiology, 122(2-3), 131-147. 72. Smith, J. C., Abdala, A. P., Rybak, I. A., & Paton, J. F. (2009). Structural and functional architecture of respiratory networks in the mammalian brainstem. Philosophical Transactions of the Royal Society B: Biological Sciences, 364(1529), 2577-2587. 73. Suffczynski, P., Kalitzin, S., & Lopez Da Silva, F. (2005). Bistability in epileptic phenomena. Bio-Algorithms and med-systems, 1(1/2), 259-266.

106 74. Takeshita, D., Sato, Y., & Bahar, S. (2007). Transitions between multistable states as a model of epileptic seizure dynamics. Physical Review E. 75(5-1), 051925. 75. Tass, P.A., & Hauptmann, C. (2007). Therapeutic modulation of synaptic connectivity with desynchronizing brain stimulation. International Journal of Psychophysiology, 64(1), 53-61. 76. Terman, D. (1992). The transition from bursting to continuous spiking in an excitable membrane model. Journal of Nonlinear Science, 2, 133-182. 77. Tobin, A. E., & Calabrese, R. L. (2006). Endogenous and half-center bursting in morphologically inspired models of leech heart interneurons. Journal of Neurophysiology, 96(4), 2089-2106. 78. Thoby-Brisson, M., Telgkamp, P., and Ramirez J. M. (2000). The role of the hyperpolarization- activated current in modulating rhythmic activity in the isolated respiratory network of mice. Journal of Neuroscience, 20(8), 2994-3005. 79. Venugopal, S., Travers, J. B., & Terman, D. H. (2007). A computational model for motor pattern switching between taste-induced ingestion and rejection oromotor behaviors. Journal of Computational Neuroscience, 22(2), 223-238. 80. Wang, X. J. (1994). Multiple dynamical modes of thalamic relay neurons: rhythmic bursting and intermittent phase-locking. Neuroscience, 59(1), 21-31. 81. Williams, S. R., Christensen, S. R., Stuart, G. J., & H¨ausser, M. (2002). Membrane potential bistability is controlled by the hyperpolarization-activated current I(H) in rat cerebellar Purkinje neurons in vitro. Journal of Physiology, 539(2), 469-483.

107 5. APPENDICES 5.1 Model of Half Center Oscillator A half-center oscillator was modeled by 38 first-order differential equations system. It represents the dynamics of two identical leech heart interneurons coupled through inhibitory synapses. The changes of the membrane potentials on both neurons are described by the following equations: dV1 ( I Na I P I CaF I CaS I h I K 1 I K 2 I KA I leak I inj (t ) I SynS I SynG ); 1 1 1 1 1 1 1 1 1 1 1 C dt dV C 2 ( I Na I P I CaF I CaS I h I K 1 I K 2 I KA I leak I SynS I SynG ), 2 2 2 2 2 2 2 2 2 2 2 dt where C is membrane capacitance and INa, IP, ICaF, ICaS, Ih, IK1, IK2, IKA are voltage-gated currents. Currents are modeled based on Hodgkin-Huxley formalism as Iion=gion*miion*hjion (V-Eion). Equations describing currents are the same as for the model of the isolated leech heart interneuron that was described in Appendices 5.1. Parameters such as capacitance, maximum ionic conductances and reversal potentials were also presented in Appendices 5.1. Synaptic currents are given by the sum of the spike-triggered and graded components, which are determined by the following equations: I SynS g synS M 1 (V1 E syn ) FSY1 1 I SynS g synS M 2 (V2 E syn ) FSY2 2 g synG (V1 E syn ) P23 I SynG 1 P23 Co g synG (V2 E syn ) P13 I SynG 2 P13 Co ,

108 where gsynS = 150 nS and gsynG = 30 nS are the maximum conductances of the spike-triggered and graded synaptic currents, correspondingly. Esyn is the synaptic reversal potential (Esyn = -0.0625 V), and Co = 10-32 [C3]. Dynamical variables FSX and FSY represent activation and inactivation of the spike-triggered synaptic current and M is the modulation variable of the current. Fast synapse variables are described by the following equations: dM 1 tau (1000,0.04,0.1,0.9, V2 ) dt 0.2 dM 2 tau (1000,0.04,0.1,0.9,V1 ) dt 0.2 dFSY1 FSX 1 FSY1 dt 0.011 dFSY2 FSX 2 FSY2 dt 0.011 dFSX 1 f (1000,0.01, V2 ) FSX 1 dt 0.002 dFSX 2 f (1000,0.01,V1 ) FSX 2 dt 0.002 where f ( A, B,V ) (1 e A(V B ) ) 1 and the time constant tau is described by D tau ( A, B, C , D,V ) C . 1 e A(V B ) Graded synapse IsynG is governed by the dynamical variables P and A. P depends on Ca2+ currents and a dynamical threshold A. P and A variables are given by the following equations:

109 dP1 max( 0,( I CaF I CaS )10 9 A1 ) 10 P1 dt dP2 max( 0,( I CaF I CaS )10 9 A2 ) 10 P2 dt dA1 10 10 f (100,0.02, V2 ) A1 dt 0.2 10 dA2 10 f (100,0.02,V1 ) A2 dt 0.2 Function max( 0,( I CaF I CaS )10 9 A1, 2 ) yields value of ( I CaF I CaS )109 A1, 2 for the condition ( I CaF I CaS )109 A1, 2 ≥0, otherwise the value of the function is zero. Bistability of HCO was demonstrated under the condition with elevated leak conductance. Figure 3.2 demonstrates the switch of the half-center oscillator activity from bursting to silent regime triggered by a negative pulse of current. Pulse was applied at the time 20.3 sec. The duration of the pulse was 0.03 sec and the amplitude was set to 0.3 nA. The initial conditions for the activity presented in Figure 3.2 are given in Table 5.1. Table 5.1: Initial conditions of HCO model. coordinates Neuron 1- HN(L) Neuron 2- HN(R) V -0.043285962477163 -0.057036910025725 mCaF 0.932234202216838 0.001725626832112 hCaF 0.006921199822378 0.638046881080213 mCaS 0.916151724457494 0.013972029923988 hCaF 0.187029801239018 0.700079996523472 mK1 0.075217838683608 0.005500600080519 hK1 0.752177336656006 0.961186533947848 mK2 0.231118880189264 0.044830626059398 mKA 0.581580647828857 0.150536165542992 hKA 0.028618531911801 0.281637973529545 mh 0.278261391037481 0.545500724334109

110 mP 0.402773326824277 0.106159406948084 mNa 0.105008802133364 0.014684695603250 hNa 0.967575377360531 0.999998833416576 FSX 0 0.000000001506794 FSY 0 0.011974737168553 M 0.100003371291043 0.454000702006309 P 0.000000000009250 0 A 0.000000000002402 0.000000000016832

## Reviews