Local dynamics of ion based neuron models for cortical spreading depression, stroke and seizures

Local Dynamics of Ion–Based Neuron

Models for Cortical Spreading

Depression, Stroke and Seizures

vorgelegt von

Diplom–Physiker Niklas H¨

ubel

geb. in Unna

von der Fakult¨

at II – Mathematik und Naturwissenschaften

der Technischen Universit¨

at Berlin

zur Erlangung des akademischen Grades Doktor der Naturwissenschaften, Dr. rer. nat.

genehmigte Dissertation

Promotionsausschuss

Vorsitzender: Prof. Dr. Michael Lehmann Gutachter: Prof. Dr. Eckehard Sch¨oll, PhD Gutachter: Prof. Dr. Dr. h.c. J¨urgen Kurths

Tag der wissenschaftlichen Aussprache: 29. August 2014

In this thesis we investigate the phase space structure of biophysical neuron models with dynamic ion concentrations. We start with the classical Hodgkin–Huxley model of an electrically excitable neural membrane with constant intra– and extracellular ion concentrations. The extension of this model to include changes in ion concentrations that result from transmembrane currents is carefully reviewed. This extension describes a closed system, in which no particle exchange with the surroundings is considered, however the neuron contains ion pumps that dissipate energy to keep it far from the thermodynamic equilibrium. Exploiting all symmetries of the model and applying one commonly used and one novel approximation of the gating dynamics, we obtain a reduced ion–based neuron model with only four dynamical variables.

The dynamics of the closed neuron system is investigated in the next part. A new stable fixed point that coexists with the physiological resting state is found. Numerical simu-lations show that the new fixed point is extremely close to the Donnan equilibrium, i.e., the thermodynamic equilibrium of the system. The neuron cannot fire action potentials in this state, because the electric energy that is usually stored in the ion gradients is almost fully dissipated. We refer to this condition as free energy–starvations (FES). This is the first bistability of neuron states with completely different intra– and extra-cellular ion concentrations ever reported. Perturbations that cause the transition from the physiological resting state to FES are, for example, long stimulations with applied currents or a temporary interruption of the pump activity.

We perform the first bifurcation analysis of a unified model for action potentials and ion dynamics. We vary the pump rate as a bifurcation parameter, and thereby prove the coexistence of a physiological resting state and FES in a large number of reduced ion–based model variants. The result is also replicated for a very biophysically detailed model developed by Kager et al.. Most importantly, the bifurcation analysis shows that a closed neuron system can only recover from FES if the rate of the ion pumps is extremely enhanced. This leads to our first major conclusion: ion homeostasis cannot rely on the pumps alone. Spreading depression (SD) is an important example of neural

ion dynamics, in which transient FES is observed. Recovery from FES during SD was long believed to be due to the pumps—our results disprove this hypothesis.

The next part deals with open neuron systems, in which potassium ions can be ex-changed with an external glial or a vascular reservoir. This changes the dynamics from bistability to what we call ionic excitability, i.e., excitability with large changes in the ion concentrations. For the first time all of the many different time scales of dynamics in open neuron systems are identified. This leads to a new slow–fast description of ion dynamics during SD. The gain or loss of potassium ions is the slowest–changing vari-able, and it causes the transition of the neuron from physiological conditions to FES or vice versa. The importance of this quantity was not realized before and a bifurcation analysis, in which it is varied as a parameter provides a phase space explanation of each process involved in local SD dynamics: the ignition of SD, the subsequent prolonged depolarized phase, the abrupt repolarization, and the final recovery. Furthermore we investigate oscillatory dynamics obtained for coupling the neuron to a bath with an elevated potassium ion concentration. The dynamical repertoire contains seizure–like activity (SLA), tonic firing and periodic SD—each associated with specific bifurcations. This is the first time seizure–like dynamics and SD are mathematically described and categorized as genuinely different dynamics.

Thereafter we direct our attention to osmosis–driven cell swelling caused by extreme ion dynamics. We combine earlier models of neural volume dynamics to derive a descrip-tion based on first physical principles—as opposed to the phenomenological approaches normally used in SD modeling. This enables us to identify critical system parameters, and, in particular, to understand the central role of anion channels in volume dynamics. Moreover, we conclude from numerical simulations that an adiabatic approximation is valid for neural volume dynamics. The general phase space structure of the system is shown not to change when cell swelling in included.

Finally, the necessity to model all time scales of neural dynamics is demonstrated in a direct physiological application. To assess the effect of a sodium channel mutation on the SD susceptibility of tissue, we perform simulations with a pure spiking model and with the reduced ion–based model. In the context of neural mass models the former would lead to partly contradicting conclusions, while only the results for the ion–based model are consistent with experiments and show an increased SD susceptibility.

In dieser Arbeit untersuchen wir die Phasenraumstruktur biophysikalischer Nerven-zellmodelle mit dynamischen Ionenkonzentrationen. Wir beginnen mit dem klassis-chen Hodgkin–Huxley–Modell zur Beschreibung elektrisch anregbarer Zellmembranen mit kontstanten intra– und extrazellul¨aren Ionenkonzentrationen. Die Erweiterung des Modells zur Beschreibung von Ionendynamik, die aus den Membranstr¨omen resultiert, wird ausf¨uhrlich erl¨autert. Diese Erweiterung f¨uhrt zu einem geschlossenen System ohne Teilchenaustausch mit der Umgebung. Energiezufluss zur Versorgung der Ionenpumpen, die das System fern vom thermodynamischen Gleichgewicht halten, ist jedoch n¨otig. Unter Ausnutzung aller Symmetrien und durch die Anwendung einer h¨aufig verwende-ten sowie einer neuartigen N¨aherung der Gatingdynamik erhalten wir ein reduziertes ionenbasiertes Neuronenmodell mit lediglich vier dynamischen Variablen.

Im n¨achsten Teil untersuchen wir die Dynamik dieses geschlossenen Systems. Dabei finden wir einen neuartigen stabilen Fixpunkt, der neben dem physiologischen Ruhezu-stand koexistiert. Numerische Simulationen zeigen, dass dieser Fixpunkt große ¨ Ahnlich-keit mit dem Donnan–Gleichgewicht, d.h. dem thermodynamischen Gleichgewicht der Zelle, aufweist. Ist das Neuron in dem neuen Fixpunkt, so k¨onnen keine Aktionspoten-tiale mehr ausgel¨ost werden, da die Energie, die normalerweise durch das große Gef¨alle zwischen intra– und extrazellul¨aren Ionenkonzentrationen bereitgestellt wird, nahezu vollst¨andig verbraucht ist. Daher bezeichnen wir den neuen stabilen Systemzustand mit ‘free energy–starvation’ (FES). Es handelt sich bei unserer Beobachtung um die erste jemals gefundene Bistabilit¨at von Zust¨anden in Neuronenmodellen mit unterschiedlichen Ionenkonzentrationen. St¨orungen, die zum ¨Ubergang von dem physiologischen Ruhezu-stand zu FES f¨uhren, sind unter anderem lange Stimulationen mit angelegten Str¨omen oder eine kurzzeitige Schw¨achung der Ionenpumpen.

Erstmalig f¨uhren wir die Bifurkationsanalyse eines vereinheitlichenden Modells f¨ur Ak-tionspotentiale und Ionendynamik durch. Als Bifurkationsparameter variieren wir die St¨arke der Ionenpumpen, und es gelingt uns die Koexistenz eines physiologischen Ruhezu-stands und eines FES–artigen ZuRuhezu-stands f¨ur eine Vielzahl von reduzierten ionenbasierten

Modellvarianten zu zeigen. Des weiteren reproduzieren wir dieses Resultat mit einem biophysikalisch ¨außerst detaillierten Modell, das Kager et al. entwickelt haben. Ein weiteres wichtiges Ergebnis der Analyse ist, dass ein abgeschlossenes Neuron sich nur von FES erholen kann, wenn die Pumpst¨arke drastisch erh¨oht wird. Dies f¨uhrt uns zur ersten bemerkenswerten Schlussfolgerung: Ionenhom¨oostase in Neuronen kann nicht ausschließlich auf den Pumpen beruhen. ‘Spreading depression’ (SD) stellt ein wichtiges Extrembeispiel f¨ur neuronale Ionendynamik dar, bei dem vorr¨ubergehend FES eintritt. Es wurde bisher vermutet, dass die Erholung von FES durch die Ionenpumpen erkl¨art werden kann—unsere Ergebnisse widerlegen diese Hypothese.

Im n¨achsten Teil betrachten wir das Neuron als ein offenes System, das Kaliumionen mit externen Reservoirs, z.B. Gliazellen oder dem vaskul¨are System, austauschen kann. Die vorherige Bistabilit¨at wird dadurch zu einer ionischen Anregbarkeit, d.h. einer An-regungsdynamik mit großen Ver¨anderungen in den Ionenkonzentrationen. Erstmalig werden die einzelnen Zeitskalen aller relevanten Prozesse in einem solchen System her-ausgearbeitet. Dies f¨uhrt zu einer Zeitskalentrennung und einer damit einhergehenden neuartigen Beschreibung der lokalen Dynamik von SD, bei der die Zu– bzw. Abnahme des Kaliumionengehalts die langsamste dynamische Variable ist. Sie bewirkt die ¨Uberg¨ange zwischen physiologischen Bedingungen und FES, und hat damit eine zentrale, bisher nicht erkannte Bedeutung. Durch eine Bifurkationanalyse, in der diese Gr¨oße als Pa-rameter variiert wird, lassen sich die einzelnen dynamischen Prozesse, die w¨ahrend SD zu beobachten sind, explizit anhand der Phasenraumstruktur erkl¨aren. Die Ausl¨osung von SD, die lang andauernde Depolarisation, die anschließende abrupte Repolarisation und die langsame abschließende Erholung k¨onnen nun mit Hilfe der Bifurkationsstruk-tur verstanden werden. Außerdem untersuchen wir station¨are Oszillationen, zu denen es kommt, wenn das Neuron an ein Kaliumbad mit erh¨ohter Konzentration gekoppelt wird. Das dynamische Repertoire umfasst dabei epilepsieartige neuronale Aktivit¨at, tonisches Feuern und periodische SD, die jeweils mit bestimmten Bifurkationen zusammenh¨angen. Diese Analyse ist die erste grunds¨atzliche mathematische Unterscheidung von SD und epilepsieartiger Dynamik.

Im Anschluss besch¨aftigen wir uns mit dem Anschwellen von Zellen durch Osmose w¨ahrend extremer Ver¨anderungen in den Ionenkonzentrationen. Wir kombinieren ex-istierende Modelle f¨ur dynamische Volumen¨anderungen und erhalten so ein Modell, das—im Gegensatz zu den den in der SD–Modellierung weit verbreiteten ph¨ anome-nologischen Ans¨atzen—auf rein physikalischen Prinzipien beruht. Diese Beschreibung erm¨oglicht es uns die f¨ur Zellschwellen entscheidenden Parameter herauszuarbeiten. Ins-besondere erkennen wir erstmalig die Ins-besondere Wichtigkeit der Leitf¨ahigkeit der Anio-nenkan¨ale. Des weiteren zeigen unsere Simulationen, dass neuronale Volumendynamik

adiabatisch approximiert werden kann. Wir zeigen anschließend, dass die generelle Phasenraumstruktur durch Volumendynamik nicht ver¨andert wird.

Abschließend demonstrieren wir anhand einer physiologischen Anwendung, dass die Modellierung aller Zeitskalen notwendig ist um neuronale Dynamik theoretisch zu un-tersuchen. Um den Einfluss eines genetischen Defekts des Natriumkanals auf die SD– Suszeptibilit¨at des Neurons zu untersuchen vergleichen wir Simulationen von einem reinen Aktionspotential–Modell mit dem ionenbasierten Ansatz. Ersteres w¨urde in einer Interpretation bez¨uglich neuronaler Felder zu widerspr¨uchlichen Schlussfolgerun-gen f¨uhren. Die Ergebnisse f¨ur das ionenbasierte Modell sind hingegen konsistent mit experimentellen Daten und belegen eine erh¨ohte Suszeptibilit¨at.

Firstly, I am very grateful to have had the opportunity to write this thesis. For this I would like to thank Prof. Sch¨oll, extending my thanks to the members of the AG Sch¨oll for creating a relaxed and cheerful working environment. Also, I have great appreciation for the continuous support Prof. Sch¨oll has shown, especially in the last few months, and for his fairness and good judgement.

Markus Dahlem and I have shared many fuitful discussions, which have in the end shaped this thesis. Moreover Markus has bridged the gap between my work and the neuroscience community, helping me to ask relevant questions and for the results of this thesis to have an impact in a wider context. Thank you Markus, not only as a grateful academic student and colleague but also as a friend.

My thanks go to Prof. Kurths for providing a second assessment of the thesis.

Finally to my parents for the ongoing support and to Thomas Isele—it was a joy to share an office with you.

• N. H¨ubel, E. Sch¨oll, and M. A. Dahlem. Bistable dynamics underlying excitability of ion homeostasis in neuron models. PLoS Comp. Biol., 10:e1003551, 2014. doi: 10.1371/journal.pcbi.1003551.

• N. H¨ubel and M. A. Dahlem. Dynamics from seconds to hours in Hodgkin–Huxley model with time–dependent ion concentrations and buffer reservoirs. arXiv:1404.3031, 2014. submitted to PLoS Comp. Biol.

• M. A. Dahlem, J. Schumacher, and N. H¨ubel. Linking a genetic defect in migraine to spreading depression in a computational model. PeerJ, 2:e379, 2014. doi: 10.7717/peerj.379.

• A. Sinha, S. J. Schiff, and N. H¨ubel. Estimation of internal variables from Hodgkin– Huxley neuron voltage. IEEE EMBS Conference on Neural Engineering, 2013.

Abstract i

Deutsche Zusammenfassung iii

Acknowledgements vi

List of Publications vii

1 Introduction 1

2 Mathematical Neuron Models 8

2.1 The Hodgkin–Huxley Model . . . 9

2.1.1 Dynamics and Bifurcations of the Hodgkin–Huxley Model . . . 12

2.2 A Reduced Ion–Based Model . . . 15

2.3 Reduced Model Variants . . . 23

2.4 A physiologically Detailed Model . . . 25

3 The Phase Space of Closed Neuron Models 29 3.1 Overview of Dynamics in the Reduced Model . . . 30

3.2 Bifurcation Analysis of Reduced Model . . . 35

3.3 Ionic Excitability . . . 37

3.4 Bifurcation Analyses of Reduced Model Variants . . . 39

3.5 Bifurcation Analysis of the Detailed Model by Kager et al. . . 45

3.6 Discussion . . . 46

4 Neuron Models with Buffer Reservoirs 50 4.1 Ion Regulation in Neurons . . . 51

4.1.1 Time Scales in Neuron Models with Buffer Reservoirs . . . 56

4.2 Potassium Gain and Loss as a Bifurcation Parameter . . . 58

4.3 Slow–Fast Decomposition of Ionic Excitability . . . 65

4.4 Ionic Oscillations For Bath Coupling . . . 69

4.5 Discussion . . . 75

5 Dynamical Aspects of Osmotic Volume Changes 79 5.1 Model for Volume Dynamics . . . 80

5.2 Most Relevant Parameters for Osmotic Equilibrium . . . 83 viii

5.3 Transient Volume Dynamics and Adiabatic Approximation . . . 86 5.4 Ion Dynamics with Volume Changes . . . 87 5.5 Discussion . . . 89 6 Application: Linking a Genetic Defect to Spreading Depression 92 6.1 Modified Time Constant . . . 93 6.2 Mutation Effect on Action Potential Dynamics . . . 95 6.3 Mutation Effect on Spreading Depression and Anoxic Depolarization . . . 98 6.4 Discussion . . . 100

7 Outlook 103

Introduction

Electrophysiological brain function is the collective result of a huge number of mecha-nisms ranging from processes at the cellular level like electrical activity and ion move-ments across the neural membrane to interactions—local and nonlocal—of neurons by synaptic connections, gap junctional coupling and diffusion of ions in the interstitial space. The neural membrane contains different types of ion channels and metabolic pumps that maintain the physiological resting state, and the neural tissue as a whole is comprised of neurons, blood vessels and glia cells—just to name a few examples of the many components and mechanisms in the brain that interact to make normal brain function possible.

Computational neuroscience complements experimental and clinical neuroscience in pro-viding the researcher with the opportunity to isolate mechanisms that cannot be isolated in an experimental setup. Simulations and analytical insights help us to interpret data and to understand the nervous system in both health and disease. The function of a single neuron, i.e., its ability to produce voltage spikes depends very sensitively on the ion concentrations in the microenvironment of the cell. This thesis will almost entirely be devoted to the computational and mathematical study of the fundamental require-ments of ionic homeostasis in single neurons. Our approach is to use a biophysical mathematical model that allows us to systematically include or leave out the different cellular and extracellular mechanisms that are believed to be relevant for ion homeosta-sis. We can thereby evaluate the importance of the specific processes and in particular develop hypotheses about when and why ion homeostasis fails, which is a central issue in pathophysiology.

The question of how ion homeostasis works and possibly fails is closely related to a phenomenon called spreading depression or spreading depolarization (SD). It was first

described by Le˜ao in 1944 [1] who experimented with cerebral cortex of rabbits and dis-covered a decrease in electrical activity by electroencephalogram (EEG) measurements in response to intense electrical stimulations. He accordingly named the phenomenon “depression of electric activity”. In his experiments the process was triggered by a local stimulus and propagated as a wave through the tissue with a velocity of a few millime-ters per minute. This is why it was called spreading depression. In subsequent studies over the next few decades a more complete list of the characteristic properties of the cellular dynamics during SD could be established. SD was shown to be associated with a slow voltage variation that reaches its peak a few seconds after the wave ignition, with a huge rise in the extracellular potassium ion concentration, and with strong falls in the extracellular concentrations of sodium and chloride ions [2]. In early experiments with isolated retina as carried out by Martins–Ferreira also the intrinsic optical signal was identified as an indicator of SD. The related changes in the transparency of the retina probes are due to cell swelling. Cell swelling could also be shown directly by measur-ing the shrinkage of the extracellular space (ECS) durmeasur-ing SD with the help of trapped impermeant ions in the ECS.

All these observations show that SD is a massive perturbation of ion homeostasis. The fact that all the changes, i.e., the changes in the slow membrane potential, in the electric activity, in the ion concentrations and in the cell volume, are reversed after a transient phase of some minutes makes SD an even more fascinating phenomenon to study with respect to neuronal homeostasis. The discovery of SD started a long series of experiments showing that SD is a very general process that occurs in many types of nervous tissue. Since neuronal homeostasis is so severely disturbed during SD, its understanding is fundamental far beyond its role in pathophysiology, if we want to explain how the brain works in general.

In parallel to the experimental advances also mathematical modeling approaches for SD were developed since the 1960s. Grafstein first described spatial potassium dynamics during SD in a reaction–diffusion framework with a phenomenological cubic rate function mimicing the local potassium release by the neurons in 1963 [3]. Another early model that focuses exclusively on the propagation properties of SD in a cellular automata network was developed by Reshodko and Bures [4]. The detailed mechanisms involved in the instigation of SD waves are not described in such modeling attempts. A model of SD propagation in cortical neuronal structures that is amenable to a more direct interpretation in terms of biophysical quantities was developed by Tuckwell and Miura in 1978 [5]. It contains both ion movements across the neural membrane and ion diffusion in the ECS. In more recent studies Dahlem et al. suggested certain refinements of the spatial coupling mechanisms, e.g., the inclusion of nonlocal and time–delayed feedback terms to explain very specific patterns of SD propagation in pathophysiological situations

like migraine with aura and stroke [6, 7]. Even more spatially extended SD models exist, but let us rather direct our attention to a very different class of neuron models that was developed in parallel without any relation to the SD phenomenology.

In the 1940s and 1950s Hodgkin and Huxley developed the first mathematical model of electrically excitable living cells. In a groundbreaking series of articles [8–12] that “the field of computational neuroscience started with” (Bard Ermentrout, in his text book on mathematical neuroscience [13]) they describe the generation and the propa-gation of action potentials in the squid giant axon. Hodgkin and Huxley describe the electrical properties of the axon by assigning a capacitance to the lipid bilayer mem-brane, and describe the highly nonlinear conductance dynamics by voltage–gated ion channels. Although originally developed for the squid axon, the Hodgkin–Huxley (HH) model became the prototype conductance–based model for excitable cells in general. The original HH model contains a sodium current with activation and inactivation gat-ing variables, and a potassium current that only has an activation variable. However, depending on the particular cell that is modeled all kinds of ion channels can be incor-porated into such a model scheme, and especially many cardiac models were developed in the HH spirit in the following years and decades [14–16]. The application to corti-cal dynamics was first mostly indirect by using the approximative FitzHugh–Nagumo model [17]. This model is a mathematical reduction of the HH model to two dynam-ical variables that captures the threshold and excitation properties of the latter, but does not allow for a direct interpretation in terms of ion channels. The reduced model ansatz of the FitzHugh–Nagumo model and the also HH–related Morris–Lecar model [18] became extremely popular in spatially extended network models for cortical activity. Neural mass and neural field models [19] which were developed since the 1970s represent another popular class of cortical models that assume even simpler local dynamics. In contrast to the original HH axon, mammalian neurons are morphologically far more complex, and also contain a much larger number of different sodium, potassium, and calcium ion channels [20]. Biophysically more detailed models that include those ion channels received an increasing attention since the late 1990s [20]. They were in partic-ular used in single neuron models for epileptoform activity and other slow modulations of spiking dynamics [21–23] that result from the specific interplay between different channel dynamics at different time scales.

The role of ion dynamics in HH–like models was first studied in the 1980s in cardiac models [15, 24], but in cortical models it lasted until 2000 when Kager, Wadman and Somjen looked at this aspect [23] in a computational model with abundant physiologi-cal detail in terms of morphology and ion channels. In fact, their model was designed for seizure–like activity (SLA) and local SD dynamics, and succeeded spectacularly in

reproducing the experimentally known phenomenology. With SD being the slowest ex-ample of ion dynamics in neurons, Kager et al. have bridged the gap between extremely fast spiking activity and a disturbance of ion homeostasis that lasts for minutes. In a series of four articles between 2000 and 2008 Kager et al. account for fundamental neu-ronal processes like spiking–induced ion accumulation, glial ion regulation, and osmotic cell swelling [23, 25–27]. After this first biophysical model of local SD dynamics other detailed models were developed [28], and those complex models were also employed in spatially extended systems [29].

The effect of cortical ion dynamics was in the following studied by Barreto, Cressman et al. in a model for epileptoform bursting modulation. They used a much simpler HH–like description that neglects morphological details entirely and only contains the classical HH ion channels [30, 31]. At the same time also HH–like models of intermediate complexity were developed to describe potassium dynamics during neural bursting [22, 32, 33], but interestingly in none of these considerably simpler models extreme ion dynamics like in SD was studied. The only exception is probably a study by Zandt et al. who describe in the simple framework of Cressman et al. what they call the “wave of death” that follows the anoxic depolarization after decapitation as measured in experiments with rats [34].

Now in this thesis we aim at investigating the entire repertoire of ion dynamics in the simplest possible HH–like model that exhibits SD dynamics. In this sense the paper by Zandt et al. can be seen as the starting point for the work presented in the following chapters, because it has inspired us to study extreme ion dynamics in a simple model and led us to modify the Cressman model so that SD dynamics is obtained. With the minor changes we made, mostly for physical consistency, this is to our knowledge the first conductance– and ion–based model for SD that employs such simple channel dynamics. The main advantage of this minimal biophysical ansatz is that it provides a good overview of the basic neural mechanisms, which makes it possible to develop a clear intuition for the phase space structure of the system. Starting from such a minimal model scheme additional physiological details can be incorporated subsequently to model specific behavior more accurately. Analyzing the interplay of just a few fundamental mechanisms is also very helpful in detecting physical inconsistencies that are, in fact, often contained in mathematical neuron models. From a practical point of view it is also much easier to carry out a mathematical bifurcation analysis of a low–dimensional model variant first and increase the complexity slowly. We are going to present the first bifurcation analysis of a conductance–based neuron model with spiking and ion dynamics. Understanding the numerical problems involved in the analysis of the simple model allows us to proceed to the detailed model by Kager et al. and apply the same

analysis. We remark that also Cressman et al. have done a bifurcation analysis of their model, however they have averaged over the fast dynamics first. So strictly speaking they do not analyze the bifurcation structure of a model containing the full range of time scales, and—more importantly—they do not investigate the fundamental stability of ion homeostasis, but instead focus on the role of ion dynamics during bursting. The thesis is structured as follows. In Chapt. 2 we introduce the mathematical neuron models we are going to analyze in the next chapters. First, the Hodgkin–Huxley model is reviewed, and the direct extension of this membrane model to include ion dynamics is discussed in detail. We comment on a few necessary modifications that we make in comparison to Cressman et al., which will turn out to be crucial for the phase space structure of the system. In addition, also other ways in which the HH model can be extended to describe ion dynamics are listed. We refer to all these models as reduced ion–based models, since they rely on a rather simple HH membrane. Furthermore, as an example of a biophysically detailed model, the model by Kager et al. is included to this chapter. We will often use this model to replicate the results we obtain for the reduced model variants. This shall support our physiological conclusions. The ion–based models in this chapter contain ion pumps, but no other mechanisms that contribute to ion homeostasis. We hence classify them as pure transmembrane ion models. From the thermodynamic perspective such models are called closed, since no ion exchange with the system surroundings is included. We remark that they are still dissipative systems through which the energy flows that is needed for the ion pumps.

In Chap. 3 we analyze the phase space of the closed neuron models and obtain rather surprising results. We find a new stable fixed point that coexists with the physiological resting state. The system can be driven into this condition by long current stimulations or a temporary interruption of the ion pump activity. The new stable state is very similar to the Donnan equilibrium, i.e., the thermodynamic equilibrium of the cell, and the neuron will remain in this state even for an increased pump activity. This is shown by direct numerical simulations and bifurcation analyses. The result is reproduced for all variants of reduced ion–based models with electrically excitable membranes, and also for the Kager et al. model. Bistabilities of neuron states with different ion configurations have not been reported either in computational or experimental studies, and since such behavior is not experimentally observed, we draw the conclusion that ion homeostasis cannot be maintained by the ion pumps alone. Since in the newly found fixed point no electrical activity is possible and the energy normally provided by the ion gradients is almost fully dissipated, we call this condition free energy–starvation (FES). We show that, if an extracellular ion regulation mechanism is added, the system dynamics changes and the response to the above mentioned stimulation types is SD–like dynamics. The system makes an excursion in the phase space that reaches very close to the peculiar

new fixed point, but returns to the initial resting state afterwards. To distinguish this excitation dynamics from the electric excitability of the membrane that underlies action potential dynamics, we refer to it as ionic excitability.

The results from Chap. 3 are extended to develop a full phase space picture of SD dynamics in Chap. 4. While the conclusion in the previous chapter was that ion pumps alone are not sufficient, we now direct our attention to the role of ion regulation in open neuron systems. We identify the ion loss or gain through the external reservoirs to which the neuron is coupled as the central quantity in ionic excitability and succeed to give a complete phase space interpretation of local SD dynamics. The ignition, the prolonged phase of suppressed electric activity, and the final recovery can all be explained in terms of well–defined thresholds in a slow–fast analysis. This allows us to describe local SD dynamics as a sequence of activation and inhibition processes at different time scales. In particular, ion exchange with reservoirs is identified as the inhibitory mechanism that makes recovery from the vicinity of the thermodynamic equilibrium possible. This reduced ion–based model also exhibits oscillatory dynamics, and we can classify three different types: periodic spreading depression, tonic firing, and seizure–like activity. Each of these dynamics sets in via specific bifurcations. This is the first time seizure– like dynamics and SD are mathematically categorized as genuinely different dynamics. Only SD must be explained by the presence of a metastable state of FES.

While the first chapters were primarily devoted to the mathematical understanding of neuronal ion dynamics we include osmotic cell swelling to the model in Chapt. 5 bearing in mind the experimental methods of measuring SD. Volume changes during extreme ion dynamics have been modeled before, but the systematical overview we are giving in this chapter includes several new results. On the one hand we identify critical model parameters that have a strong impact on volume dynamics, and on the other hand we find out that certain parameters surprisingly do not. In particular, the conductances of anion channels in the membrane play an important role, while the details biophysical mechanism that leads to volume changes can be ignored. This observation allows us to eliminate the rate equation for the cell volume adiabatically and thereby derive a model for neuronal ion dynamics and cell swelling that has only five dynamical variables, but is still amenable to a biophysical interpretation. We confirm that the general phase space structure remains qualitatively unchanged.

In Chapt. 6 the practical relevance of the reduced ion–based neuron models that we have reviewed, that we have partly developed and that we have analyzed throughout Chapts. 2–5 is shown in modeling a genetic defect that leads to a higher probability of migraine with aura, a pathology closely related to SD. This is an example for how a biophysical modification can be incorporated in such models, and how this leads to a clinically reported effect. In Chapt. 7 we conclude the thesis with a brief outlook on future research projects that may follow our here presented results. First steps that have been taken to investigate the observability properties of our model are presented.

Mathematical Neuron Models

In this chapter we are going to review some of the most common mathematical models of neuronal dynamics. For this purpose we will first review the classical conductance– based Hodgkin–Huxley (HH) model [12]. This model describes the electrical spiking activity of neurons under healthy physiological conditions (see Sect. 2.1). Its dynamical variables are the membrane potentialV and three so–called gating variables n, m and h. Ion concentrations within and outside the cell are fixed model parameters. The most important dynamical properties and the bifurcation structure are briefly reviewed. The following Sect. 2.2 is devoted to the extension of the HH model to include ion dy-namics. Our presentation is a review of a paradigmic model used in Refs. [30, 31, 34] with some modifications that we explicitly comment on. The presented ion–based description of neural dynamics accounts for changes in ion concentrations due to fluxes across the neural membrane. Such fluxes are induced by the electrical activity of the membrane. We carefully review that such an extension requires the introduction of ATPase–driven ion pumps that exchange intracellular sodium ion for extracellular potassium ions and thereby stabilize the physiological state. So in this model ansatz we describe the ion dy-namics that follows directly from the HH model. In this sense the model differs strongly from physiologically very detailed model types that contain a considerably larger number of different ion channels and sometimes also consider morphological details like dendrites and dendritic branching. One such model is presented in the last section. We hence call the HH–like model of Sect. 2.2 a reduced or a minimal ion–based neuron model. We remark that whenever we refer to sodium, potassium or chloride we mean the ions Na+, K+ and Cl− also if we omit the word ‘ion’ for better readability.

After setting up the model we discuss some immediately obvious symmetries and point out how the different time scales in the model imply that electroneutrality is a constraint for the changes of ion concentrations inside the neuron. This reduced ion–based model

or variants of it will be the object of later detailed phase space analyses that reveal the mechanisms underlying different pathological ion dynamics.

While the reduced model of Sect. 2.2 is obtained directly as an extension of the HH model, the way this extension can be done is not unique. We hence comment in Sect. 2.3 on different variants of reduced ion–based models.

In the last Sect. 2.4 of this chapter we present the model equations of a widely accepted neuron model that Kager et al. [23] have developed to investigate the dynamics of cortical spreading depression (SD) and seizure–like activity. This physiologically very detailed model contains a large number of dynamical variables and physiological parameters, and we will often use this model to replicate the results we obtain for reduced ion–based models.

2.1

The Hodgkin–Huxley Model

Figure 2.1: Illustration of the HH model taken from Ref. [13]. The cell membrane (left) consists of an insulating lipid bilayer with ion channels. In the equivalent electrical circuit the lipid bilayer is a capacitor and the channel is a conductor in series with a battery. Note that this scheme contains a chloride leak current instead of the unspecified

leak current. For details see Ref. [13].

The Hodgkin–Huxley model is one of the most widely used computational models in neuroscience. It is a conductance–based model [12] that describes single neuron dynamics in terms of an electrically active membrane carrying an electric potentialV , and the three gating variablesn, m and h that render the system excitable. Fig. 2.1 shows a scheme of the model taken from Ref. [13]. In the equivalent electrical circuit description the neuronal membrane acts as a capacitor, the ion channels are conductors and the ion gradients are batteries.

Ion species included are sodium Na+, potassium K+, and an unspecified ion carrying a leak current. Altogether the model contains four differential equations describing the dynamics of the membrane potential V according to Kirchhoff’s current law, and the

voltage–gated dynamics of the three gating variables that determine the conductances of the potassium and sodium channels. The HH model has been developed to describe action potentials, also known as nerve pulses or spikes. The mechanism that underlies the neuron’s capability to elicit voltage spikes in response to a sufficient stimulation is the excitability which results from the interplay between the gating variables and the much faster membrane potential.

According to the Kirchhoff current law for electrical circuits a change of the membrane potential is proportional to the current that is flowing across the membrane with the proportionality constant given by the capacitanceCm of the membrane. Cm is assumed

to be constant which might not be always be the case in real neural systems, but is assumed throughout this thesis. The individual currents are computed as the product of the conductance gion of the respective ion channel and the driving force, which is

given by the difference between the membrane potential and the respective ion’s reversal potentialEion, where ion∈ {K+, Na+}.

In addition to the specific ion currents for sodium and potassium there is an unspecific leak current Ileak with a conductance gleak and a reversal potential Eleak. The

conduc-tance gion for a voltage–gated channel, i.e., for sodium and potassium, is given by the

product of the maximal conductance g_ion and a combination of gating variables as de-scribed below. So the currents that the HH model takes into account are a sodium current INa, a potassium current IK, a leak current Ileak that is carried by unspecified

ions, and an applied currentIapp:

dV dt = − 1 Cm (INa+IK+Ileak− Iapp) , (2.1) INa = gNam3h(V − ENa), (2.2) IK = gKn4(V − EK) , (2.3)

Ileak = gleak(V − Eleak) . (2.4)

Here the potassium current is modeled as a delayed rectifier current with an activation variablen, while the sodium current is described by a transient current with an activation variable m and an inactivation variable h. All gating variables are voltage–dependent

and their dynamics is of the following from: dx dt = x_∞− x τx with (2.5) x_∞ = αx αx+βx and (2.6) τx = 1 αx+βx forx_{∈ {n, m, h}.} (2.7)

The steady–state values x_∞ of the gating variables and the time constantsτx are given

in terms of the following exponential functions that Hodgkin and Huxley estimated:

αm = 0.1(V + 40) 1_{− exp((V + 40)/10)} , (2.8) βm = 4 exp((V + 65)/18) , (2.9) αn = 0.01(V + 55) 1_{− exp((V + 55)/10)} , (2.10) βn = 0.125 exp((V + 65)/80) , (2.11) αh = 0.07 exp((V + 65)/20) , (2.12) βh = 1 1 + exp(−(V + 35)/10) . (2.13)

This model is capable of producing action potentials in response to appropriate exter-nally applied, depolarizing currents Iapp. A remark about the Nernst potentials is due

at this point. The Nernst–Planck equation describes the simultaneous effect of diffusive flux and electrical drift on ions in two compartments that are separated by a membrane. It follows from this description that the equilibrium membrane potential for a given set of ion concentrations is

Eion=−

RT

zF ln([ion]i/[ion]e) (2.14) whereR is the ideal gas constant, T is the absolute temperature, F is Faraday’s constant andz is the valence of the considered ions. This potential is normally called the reversal potential or the Nernst potential of the ion channel. The ion concentrations are denoted by [ion]_i/e with subscripts i and e for the intra– and extracellular compartment. We will now omit the square brackets for better readability. At a temperature of 36◦C the fraction of the constants in front of the logarithm becomes RT /F = 26.64mV, so the Nernst potentials are

Eion=−

26.64

z ln(ioni/ione). (2.15) In the HH model ion concentrations do not change, and consequently also the Nernst potentials are constants, too. Normally the latter are given as model parameters rather

Table 2.1: Parameters for Hodgkin–Huxley model

Name Value & unit Description

Cm 1µF/cm2 membrane capacitance

g_Na 120mS/cm2 max. sodium conductance g_K 36mS/cm2 max. potassium conductance gleak 3mS/cm2 leak conductance

ENa 50mV sodium reversal potential

EK −77mV potassium reversal potential

Eleak 54.4mV leak reversal potential

than the ion concentrations themselves. A set of typical model parameters can be found in Tab. 2.1.

2.1.1 Dynamics and Bifurcations of the Hodgkin–Huxley Model

The dynamics and the phase space structure of the HH model have been studied systema-tically and in detail, and we will now only briefly review the most important dynamical features. The plots in Fig. 2.2 show time series for a single spike (upper row) and for oscillatory spiking (lower row). Fig. 2.2a–c shows data of one simulation that was done for a 3msec long current pulseIapp with amplitude 3µA/cm2 (indicated by the grey bar

ranging from 10 to 13msec). In response to the stimulation the membrane potential slightly increases (see Fig. 2.2a). After the stimulus has ceased the potential increases further and the membrane elicits a so–called voltage spike with a peak value of about 40mV. It then repolarizes and returns to its initial state after a refractory phase that lasts about 10msec.

Fig. 2.2b and 2.2c show the gating dynamics during this single spike excursion. In reduc-tions of the HH model a linear dependence of the sodium inactivationh on the potassium activation n is often assumed, and Fig. 2.2b shows that this is a valid approximation. Another often used reduction of the model is to assume that the gating variablem takes its adiabatic value m_∞, and we see that this is a good approximation, too (see black solid vs blue dashed line in Fig. 2.2c). The plots in Fig. 2.2d–2.2f present simulation data for a constant applied current Iappwith an amplitude of 10µA/cm2 which changes

the dynamics to stationary oscillations. The frequency is of the order 20msec and the typical approximations for the gating variables are as valid as for the single spike sim-ulations. We will use both approximations in the below described minimal ion–based model to reduce the numerical effort as much as possible. However, we will have to use a sigmopidal fit rather than a linear relation betweenh and n (see Sect. 2.2).

The transition between different types of dynamics in the model, i.e., the onset and the cessation of oscillations for certain values ofIapphappens in bifurcation scenarios. With

0 5 10 15 20 25 30 t/ msec. −100−80 −60 −40 −200 20 40 V /mV membrane potential 0 20 40 60 80 100 t/ msec. −100−80 −60 −40 −200 20 40 V /mV membrane potential 0.0 0.2 0.4 0.6 0.8 1.0 n 0.0 0.2 0.4 0.6 0.8 1.0 h sodium inactivation 0.0 0.2 0.4 0.6 0.8 1.0 n 0.0 0.2 0.4 0.6 0.8 1.0 h sodium inactivation 0 5 10 15 20 25 30 t/ msec. 0.0 0.2 0.4 0.6 0.8 1.0 m sodium activation 0 20 40 60 80 100 t/ msec. 0.0 0.2 0.4 0.6 0.8 1.0 m sodium activation dynamic adiabatic a) b) c) d) e) f)

Figure 2.2: In the upper row plots (a), (b) and (c) simulation data of the HH model for stimulation with a 3msec long current pulse with amplitude 3µA/cm2is shown. In (a) and (c) the stimulation time is indicated by the grey bar. In the lower row plots (d), (e) and (f ) a constant current of amplitude 10µA/cm2 is applied. (a) and (d) show the characteristic evolution of the membrane potential during spiking activity. (b) and (e) justify the often used approximation that assumes a linear dependence between the gating variables h and n. In (c) and (f ) we compare the value of the gating variable m (solid black line) to its adiabatic value (blue dashed line) and see

that the reduction of replacingm by m∞is also a valid approximation.

continuation software tools like AUTO [35] one can track the fixed point of a dynamical system and detect at which parameter values bifurcations lead to stability changes. The generic cases are saddle–node bifurcations (also called limit point bifurcations, LP) and Hopf bifurcations (HB). In a LP the stability changes in one direction (zero–eigenvalue bifurcation), in a HB it changes in two directions and a limit cycle is created (complex eigenvalue bifurcation). In bifurcation diagrams a limit cycle is usually represented by the maximal and minimal value of the dynamical variables. Amongst other possibilities limit cycles can change their stability in saddle–node bifurcations of limit cycles (limit point bifurcations of limit cycles LPlc).

In the model stability changes of the fixed point occur only in Hopf bifurcations. Fig. 2.3 shows bifurcation diagrams of the HH model for two different bifurcation parameters, the applied current Iapp and the potassium Nernst potential EK. The fixed point is

presented as a black line, stable sections are solid, unstable sections are dashed. The maximal and minimal membrane potential of the limit cycle is drawn as a green line with the same stability convention of line styles. Bifurcations are marked by circles and have numbered labels indicating the type of bifurcation. The insets show blow–ups of the limit cycle bifurcations that occur in a very narrow parameter range.

Both bifurcation diagrams look very similar. For the Iapp–continuation the initially

0 50 100 150 200 250 Iapp/ (µA/cm2) −100 −80 −60 −40 −20 0 20 40 V /mV HB1 HB2 LP1lc LP2lc LP3lc

HH bifurcation diagram for Iapp

stable FP unstable FP stable LC unstable LC −70 −65 −60 −55 −50 −45 −40 −35 −30 EK/ mV −100 −80 −60 −40 −20 0 20 40 V /mV HB1 HB2 LP1lc LP2lc LP3lc

HH bifurcation diagram for EK

6.5 7.0 7.5 8.0 −80 −20 40 LP1lc LP2lc LP3lc −67.7 −67.5 −67.3 −80 −20 40 LP1lc LP2lc LP3lc a) b)

Figure 2.3: Bifurcation diagrams for the HH model. In (a) Iapp is used as the

bifurcation parameter, in (b) EK is used. The fixed point (FP) is presented as a

black line, stable sections are solid, unstable sections are dashed. The maximal and minimal membrane potential of the limit cycle (LC) is drawn as a green line with the same stability convention of line styles. Bifurcations are marked by circles and have numbered labels indicating the type of bifurcation. The insets show blow ups of the limit cycle bifurcations that occur in a very narrow parameter range. In (a) the parameter values for the bifurcation points are 9.78 (HB1), 154.77 (HB2), 7.85 (LP1lc),

7.93 (LP2lc) and 6.28µA/cm2 (LP2lc). In (b) they are -66.87 (HB1), -49.55 (HB2),

-67.39 (LP1lc), -67.30 (LP2lc) and -67.61mV (LP3lc).

unstable limit cycle that is created in HB1 undergoes three LPlc of which only the last,

i.e., LP3lc, makes it stable. The now stable limit cycle continues until it disappears

in the second Hopf bifurcation of the fixed point, i.e., in HB2. It follows from this bifurcation diagram that—due to the Hopf scenario—oscillatory spiking sets in with finite frequency and also—due to the hysteresis between HB1 and LP3lc—with finite

amplitude. The bifurcation analysis in Fig. 2.3 yields exactly the same sequence of bifurcations whereas now the parameter is EK. We have not contained time series

for oscillations caused by shifted Nernst potentials in Fig. 2.2, but it is an important property of the model that also changes in Nernst potentials can lead to oscillations. In fact, for models with dynamic ion concentrations ion fluxes make the Nernst potentials change slowly which then induces spiking activity, because the membrane dynamics is governed by a bifurcation structure as in Fig. 2.3b. For both bifurcation parameters the depolarized regime after the cessation of spiking activity at HB2 is normally referred to as the depolarization block.

ECS

I

I

I

ICS

Figure 2.4: Model scheme of an ion–based neuron model. The membrane (dashed line) separates the intracellular space (ICS) from the extracellular space (ECS). The different potassium and sodium ion concentrations in the compartments are indicated by the number of red and blue circles. Currents through the ion gates, i.e. IK and

INa reduce the ion gradients, the pumpIp maintains them. Ions cannot escape to the

surroundings.

2.2

A Reduced Ion–Based Model

Simple ion–based neuron models can be obtained as a natural extension of the HH model. While in the original HH model ion concentrations are model parameters, in ion–based modeling intra– and extracellular ion concentrations become dynamical variables, which causes the Nernst potentials to be dynamic (cf. Eq. (2.15)).

Any neuronal membrane model can straightforwardly be extended to make ion con-centrations dynamicm, since currents induce ion fluxes. However, under the equilibrium conditions of the HH model discussed in the previous section neitherIK = 0 norINa = 0.

Hence we need to include ion pumps [36] to make sure that the rates of change of the in-tracellular (i) and exin-tracellular (e) ion concentrations vanish in the physiological resting state ( ˙Na_i/e= ˙K_i/e= ˙Cl_i/e= 0). The model is schematically presented in Fig. 2.4. We will see below that with the inclusion of ion dynamics the stability of the physiological resting state relies on the compensation of the sodium and potassium leak currents by the ion pumps. So it is necessary to now specify the before unspecified leak current so that it is a mix of sodium, potassium and chloride currents. Hence we cannot directly extend the HH model from Sect. 2.1, but instead review how this extension was carried out by Cressman et al. [30] with a few modifications that we comment on.

The first difference between Ref. [30] and the classical Hodgkin–Huxley model from Sect. 2.1 concerns the leak currents. A leak–only chloride currentICl replaces Ileak, and

also leak conductances for sodium and potassium are included. Furthermore also the pump currentIp enters the dynamics of the membrane potential V :

dV dt = − 1 Cm (IN a+IK+ICl+Ip− Iapp), (2.16) IN a = (gN al +gN ag m3h)· (V − EN a) , (2.17) IK = (gKl +gKgn4)· (V − EK), (2.18) ICl = glCl· (V − ECl) , (2.19)

with g_ionl/g denoting the leak and gated conductances of the respective ion species. As long as ion dynamics is not considered the usage of a general unspecific leak currentIleak

is mathematically equivalent to specifying it as being partially sodium, potassium and chloride, while in the context of ion dynamics this simplification is no longer possible. The dynamics of the gating variables in the Cressman model is given by slightly shifted exponential functions as compared to Eqs. (2.8)–(2.13):

αm = 0.1(V + 30) 1− exp(−(V + 30)/10) , (2.20) βm = 4 exp(−(V + 55)/18) , (2.21) αn = 0.01(V + 34) 1− exp(−(V + 34)/10) , (2.22) βn = 0.125 exp(−(V + 44)/80) , (2.23) αh = 0.07 exp(−(V + 44)/20)) , (2.24) βh = 1 1 + exp(_{−0.1(V + 14))} . (2.25) Gating approximations form and h like the ones we discussed in Sect. 2.1.1 (cf. Fig. 2.2b, c, e, f) can also be used for this model. However, the assumed functional relation between h and n is usually realized as a linear fit [13] and cannot be adopted for our purpose. We suggest to replace it with the following sigmoidal fit to make sureh is non–negative:

h = hsig(n) = 1−

1

1 + exp(_{−6.5(n − 0.35))} (2.26) This is necessary, because ion–based models contain stable or metastable states with open potassium gates, i.e., large values of n, and a linear best fit would then lead to a negative value for the h–gate. No approximation for h was used in Ref. [30], but it reduces the numerical effort of simulations and bifurcation analyses significantly, and at the same time does not change the ion dynamics qualitatively. In addition to this new

constraint we adopt the adiabatic reduction

m = m_∞(V ) (2.27)

from Ref. [30], where m_∞ is defined as in Eq. (2.6). So the dynamics of the membrane is defined by rate Eq. (2.16) for the membrane potential and

dn dt =

n_{∞− n} τn

, (2.28)

where again n_∞ is the standard combination of exponential functions given Eq. (2.6), butτn contains an additional time scale parameter ϕ in comparison to Eq. (2.7):

τn=

1 ϕ(αn+βn)

(2.29)

The two rate Eqs. (2.16) and (2.28) and the gating constraint Eqs. (2.26) and (2.27) fully describe the membrane dynamics. So after the gating reduction the only remaining variables for the membrane areV and n.

The pump current Ip that occurs in the rate equation for the membrane potential

Eq. (2.16) is the result of the ATPase–driven exchange of sodium from the intracel-lular space (ICS) with potassium from the extracelintracel-lular space (ECS) at a 3/2–ratio. The necessity of ion pumps can be seen as follows. If changes of ion concentrations were only proportional to the currents listed in Eqs. (2.17)–(2.18) thenIion= 0 for each

ion _{∈ {Na}+, K+_{, Cl}−_{} would be a necessary condition for stable ion concentrations.}

Since the conductances or sums of conductances in Eqs. (2.17)–(2.18) are all strictly positive the only way to satisfy this condition is when all Nernst potentials and the membrane potential are equal, which would in fact be a thermodynamic equilibrium. This is not what happens in neurons under physiological conditions, but instead the ion pumps maintain a resting state which is most importantly characterized by different Nernst potentials for sodium and potassium. The ion pumps prevent the convergence of the initial Nernst potentials to one common intermediate value and thus keep the system far from thermodynamic equilibrium by dissipation of chemical energy, i.e., by ATP hydrolysis. Their activity increases in response to an elevation of the concentration of sodium ions in the ICS and potassium ion in the ECS. Chloride is not pumped. The pump model in Ref. [30] is given by

Ip(Nai, Ke) =ρ  1 + exp 25− Nai 3 −1 1 + exp (5.5− Ke) −1 . (2.30)

Since all types of transmembrane currents, i.e., also pumps lead to changes of the mem-brane potential, the net pump current Ip must be added in the rate Eq. (2.16) for the

membrane potential. The dynamics of the ion concentrations is now the result of HH currents through the sodium, potassium and chloride channels and the Na+/K+–pumps:

dNai dt = − γ ωi(INa+ 3Ip ) , (2.31) dKi dt = − γ ωi (IK− 2Ip) , (2.32) dCli dt = + γ ωi ICl . (2.33)

The factor γ in the above equations converts currents to ion fluxes and depends on the membrane surface areaAm and Faraday’s constant F :

γ = Am

F (2.34)

Dividing the ion fluxes by the volume ωi of the ICS gives the change rates for the ion

concentrations in the ICS.

We remark that in Ref. [30] and also in a later study by Barreto and Cressman [31] no dynamics of chloride is included, even though the chloride current enters the rate equa-tion of the membrane potential. This is physically inconsistent and, in fact, changes the phase space dramatically as we will discuss below. Hence we modify the model to include Eq. (2.33) to our list of rate equations. Furthermore in Refs. [30, 31] intra-cellular potassium is not modeled as an independent dynamical variable, but instead approximated by

Ki = 140.0 + (18.0− Nai) . (2.35)

We do not impose such a constraint from the start, but we will see below that a very similar relation for the intracellular ion concentrations follows from the time scales of the system.

Until now we only consider ion dynamics that is mediated by transmembrane currents— either through ion channels or because of pumping intracellular sodium and extracellular potassium. This means that the system that is modeled is closed and mass conservation holds. Consequently ion concentrations in the ECS can be computed from those in the ICS:

ione= ion0e+

ωi

ωe

(ion0_{i − ion}i) , (2.36)

So altogether the dynamics for the ion–based model are given by the rate Eqs. (2.16), (2.28), and (2.31)–(2.33). They are complemented by the gating constraints Eqs. (2.26) and (2.27) and the mass conservation constraints

Nae = Na0e+ ωi ωe (Na0_{i − Na}i), (2.37) Ke = Ke0+ ωi ωe (K_i0_{− K}i) , (2.38) Cle = Cl0e+ ωi ωe (Cl0_{i − Cl}i) . (2.39)

Dynamic ion concentrations imply that the Nernst potentials in the current Eqs. (2.17)– (2.19) are no longer model parameters, but are now also dynamic (cf. the definition of the Nernst potential in Eq. (2.15)). In addition to the parameters of membrane dynamics that are already contained in a classical HH–like model, i.e., conductancesgg/l_ionof gated and leak channels, and the membrane capacitance Cm, also morphological parameters

must now be specified. Ion dynamics in our model depends on the membrane surface area Am and compartmental volumes ωi/e which are all given in Tab. 2.2.

The numerical values for Am and ωi are taken from the rather physiologically detailed

model by Kager et al. [23], but depending on the specific neuron type these parame-ters vary. The extracellular volume was then chosen to be in agreement with recent experimental data, so it differs a slightly from Ref. [23].

In cortical ion–based models, it is customary to use the extracellular volume fraction

fe=

ωe

ωe+ωi

(2.40)

for a description of volumes. In computational models its value ranges from 13% in Ref. [23] to 33% in Ref. [34]. In experimental studies, the extracellular volume fraction is usually defined with respect to the whole tissue, i.e., a system comprised of neurons, glia cells and the extracellular space, which is a different definition than Eq. (2.40). A common experimental value is about 20% which can, however, increase to 27% in a certain type of focal cortical dysplasias (cf. Ref. [37]) or even to 32% during sleep, if we transfer the observations from mouse data to humans (cf. Ref. [38]). Assuming equally sized neuronal and glial volume fractions of 40% each, an experimentally measured value of 20% would in our model, which does not include the volume of the glial syncytium, correspond to fe = 0.33 or 33%. By setting fe to 25% we have chosen an intermediate

value of the here cited range of volume fractions. The influence of this choice is addressed in a specific analysis of the dependence of the phase space structure onfe (cf. Sect. 3.4).

In Tab. 2.2 these morphological parameters are given in the commonly used units which are appropriate to their order of magnitude. We prefer this to unifying all parameters, so the cell volume is given in µm3, while ion concentrations are typically given in ‘mil-limolar’. A ‘molar’ is the amount of ions measured in moles divided by one liter, i.e., concentrations are given in units of mM = mmol/l. All parameters of the ion–based model are given in Tab. 2.2. Due to the choice of units γ from Tab. 2.2 must be multi-plied by a factor 10 in the computer code to correctly convert current densities (given in µA/cm2) to change rates of ion concentrations. All other numerical values of parameters in Tab. 2.2 can be directly used.

We will address the time scales of ion–based neuron models more explicitly in Sect. 4.1.1 of Chapt. 4, but it is noteworthy already at this point that the extremely small value of γ/ωi makes the ion dynamics five orders of magnitude slower than the membrane

dynamics, i.e., the dynamics of V and n. Much slower ion dynamics should, of course, be expected since in the original HH model ion concentrations are even assumed to be constant.

Table 2.2: Parameters for reduced ion–based model

Name Value & unit Description

Cm 1µF/cm2 membrane capacitance

ϕ 3/msec gating timescale parameter gl_Na 0.0175mS/cm2 sodium leak conductance

gg_Na 100mS/cm2 max. gated sodium conductance gl

K 0.05mS/cm

2 _{potassium leak conductance}

gg_K 40mS/cm2 max. gated potassium conductance gl_Cl 0.05mS/cm2 chloride leak conductance

ωi 2,160µm3 volume of ICS

ωe 720µm3 volume of ECS

F 96485C/mol Faraday’s constant Am 922µm2 membrane surface

γ 9.556e–3µm2_Cmol conversion factor ρ 5.25µA/cm2 max. pump current P_Nal 0.0264µm/sec leak sodium permeability P_Nag 150.77µm/sec gated sodium permeability Pl

K 0.0169µm/sec leak potassium permeability

P_Kg 13.488µm/sec gated potassium permeability P_Cll 0.0521µm/sec leak chloride permeability

The model that we have derived in this section has a physiological equilibrium that is comparable to that of the classical HH model. The values of the dynamical variables and derived quantities like Nernst potentials and currents are listed in Tab. 2.3. Unless otherwise stated we will use these physiological resting state values as initial conditions in our simulations. We remark that the main task of the ion pumps is to compensate for the leak currentsI_Nal and I_Kl , while the gated currentsI_Nag andI_Kg practically vanish

in the normal resting state (see Tab. 2.3). The leak and the gated currents have the obvious definitions: I_Nal = gl_{Na· (V − E}Na) , (2.41) I_Nag = gg_Nam3h_{· (V − E}Na), (2.42) I_Kl = gl_K· (V − EK) , (2.43) I_Kg = gg_Kn4· (V − EK) . (2.44)

Table 2.3: Physiological resting conditions for reduced ion–based model

Name Value & unit Description

V −68mV membrane potential n 0.065 sodium activation

Nai 27mM ECS sodium concentration

Nae 120mM ICS sodium concentration

Ki 131mM ECS potassium concentration

Ke 4mM ICS potassium concentration

Cli 9.7mM ECS chloride concentration

Cle 124mM ICS chloride concentration

EN a 39.7mV sodium Nernst potential

EK −92.9mV potassium Nernst potential

ECl −68mV chloride Nernst potential

Il

Na −1.886µA/cm

2 _{sodium leak current}

I_Nag −0.011µA/cm2 gated sodium current I_Kl 1.247µA/cm2 potassium leak current I_Kg 0.018µA/cm2 gated potassium current Ip 0.632µA/cm2 pump current

We will investigate the dynamics of this model by means of direct numerical simulations of time series, but we will also apply bifurcation analyses to systematically scan the dynamical regimes. However, prior to a bifurcation analysis we need to discuss a con-servation law (symmetry) of the rate Eqs. (2.16), (2.28) and (2.31)–(2.33). The direct extension of the membrane model to include ion dynamics as presented above naturally leads to a linear dependence of the dynamical variables. In our case this is reflected by the following relation

Cm dV dt = ωi γ  dNai dt + dKi dt − dCli dt  , (2.45)

forIapp= 0. As a consequence, the determinant of the Jacobian is always zero and the

system is nowhere hyperbolic. For the continuation techniques used by software tools like AUTO, however, the Jacobian of the system must be invertible, so they cannot be applied to the system unless this degeneracy is resolved. Furthermore the phase

space structure of such nonhyperbolic systems can be changed with arbitrarily small perturbations, which is why they are called structurally unstable [39].

The physiological view on the structural instability should be as follows. Assume that the system is in its physiological resting state and apply a constant current Iapp to the

voltage rate Eq. (2.16). Then the ion rate Eqs. (2.31)–(2.33) and Eq. (2.16) imply that the equilibrium conditions ˙V = 0 and ˙Ki = ˙Nai = ˙Cli = 0 are contradictory, and hence

the equilibrium will vanish even for arbitrarily small currents. In fact for any constant, positiveIapp, the system will evolve in a highly non–physiological manner withKi, Nai

and Cle slowly tending to zero.

To avoid such unphysiological behavior we will now use Eq. (2.45) to reduce the system and thereby make it structurally stable. We can, for example, eliminate V and express it in terms of the ICS ion concentrations rather than treating V as an independent dynamical variable: d dt  V − ωi Cmγ(Nai+Ki− Cli)  = 0 ⇒ V = V(0)+ ωi CmγNai− Na 0 i +Ki− Ki0− Cli+ Cl0i  (2.46)

This was also done in Ref. [36]. The physiological meaning of this reduction is simply that the possibility of unspecified applied currents is ruled out. For instance, a perturbation on the sodium rate Eq. (2.31) should be interpreted as a sodium current. The above constraint describes the simultaneous effect on V . It would be equivalent to apply perturbations to the rate equations for sodium and the membrane potential consistently to model the full effect of an applied sodium current. The constraint Eq. (2.46) should hence be seen as a consistency condition.

This consistency rule does not at all change the dynamics unless unspecified currents are applied, and even then it practically does not change the dynamics, because any deviation in ion concentrations scales with the extremely small factorγ/ωi and is hence

negligible unless the current is constantly applied. The structural instability is thus a rather formal feature of the degenerate model and we remark that its physiological equilibrium is nevertheless stationary. Instabilities that lead to an unphysiological drift of ion concentrations for very long simulation times have been reported and resolved in cardiac cell models [15, 24]. Our case is different though, because the physiological state is a stationary one and the response to moderate stimulation is physiologically realistic. For the simulations and bifurcation analyses presented in this paper we have eliminated Nai rather than V . This is for numerical reasons, because V , n and ioni have all

different time scales and the stiff solvers we are using are much more reliable when no fast dynamics is hidden in the constraints. So in our reduction we have replaced the

sodium rate Eq. (2.31) by the following constraint:

Nai = Na0i − Ki+Ki0+ Cli− Cl0i +

Cmγ

ωi

(V − V0) . (2.47)

This relates to the approximation Eq. (2.35) for intracellular potassium that was used in Refs. [30, 31]. First note that the contribution of the membrane potential V to Nai

scales withγ/ωiand can hence be ommitted. Then, if no chloride dynamics is included,

the above condition can be simplified to

Nai = Na0i − Ki+Ki0 . (2.48)

With the right choice of initial conditions this is equivalent to Eq. (2.35). In particular the sum of intracellular sodium ion and potassium ion concentrations is then constant. Since the chloride concentration is set to be constant, the entire intracellular concen-tration of electrical charge is conserved and ion fluxes across the membrane are elec-troneutral. However, we consider dynamics of chloride, too, and will use the constraint Eq. (2.47). Electroneutrality also holds in this case as we will see in Sect. 3.1.

The non–degenerate model is now defined by the rate Eqs. (2.16), (2.28), (2.32) and (2.33), and the constraint Eqs. (2.26), (2.27), (2.37)–(2.39) and (2.47). So due to gating approximations, mass conservation, and degeneracy the initial set of dynamical vari-ables has been reduced to a total of four: the membrane potential V and the sodium activationn represent the membrane dynamics, and Ki and Cli are the two independent

ion concentrations. It is noteworthy that already the non–reduced HH model has as much as four dynamical variables, while the here presented model covers both, action potential and ion dynamics with different but not more independent variables.

2.3

Reduced Model Variants

In the previous section we have presented the derivation of one particular reduced ion– based neuron model. There is a great degree of freedom to this extension and we will cover certain model variants so that we can confirm the robustness of the results obtained for the initial model from Sect. 2.2. This robustness test will be done in Sect. 3.4. One important point is to compare Nernst models to Goldman–Hodgkin–Katz (GHK) models. In Sect. 2.2 the ion dynamics makes Nernst potentials time–dependent. For a more accurate description of the simultaneous effect of a diffusive and an electrical force acting on a solution of ions the GHK equation should be used though. Yet Nernst cur-rents are preferable in another respect: Nernst curcur-rents are defined by the conductances of the ion channels which are the common empirically estimated parameters, while for

the GHK model data about the permeabilities of the specific ion channels is needed, which is much harder to find. Nevertheless we should try to get a qualitative idea of the effect of our choice to use the Nernst model.

For a systematic comparison of the current models, we will convert the Nernst currents in Eqs. (2.17)–(2.19) to GHK currents. The natural way to do this is to assume that the ion concentrations in the physiological resting state (see Tab. 2.3) are the same. To make sure that these ion concentrations are also stable for the GHK model we require the GHK ion currents to equal their Nernstian counterparts under physiological resting conditions (see also Tab. 2.3) and choose the permeabilities accordingly. For example, the GHK version of the sodium current is

I_NaGHK = (P_Nal +P_Nag m3h)_{· F ·} z

2_V

26.64 ·

Nai− Naeexp(−zV/26.64)

1_{− exp(−zV/26.64)} (2.49) with membrane permeabilities Pl

Na and P

g

Na instead of conductances. The valence z

equals one for sodium. To compute the permeabilities we set the above expression equal to its Nernstian version

INa= (gNal +gNag m3h)· (V − ENa) (2.50)

for the physioligical resting conditions. This leads to a common conversion factor from conductancesgl/g_Na to permeabilities P_Nag/l. With this ansatz we obtain conversion factors for the three different ion species that lead to the permeabilities in Tab. 2.2. We remark that both current models are used in the computational neuroscience literature.

Another ambiguity in ion–based models is the role of the chloride leak current. Many ion–based models use a physically inconsistent leak current with fixed ion concentrations that lead to a fixed Nernst potential. We have resolved this inconsistency by adding chloride dynamics to the model. This is not in accordance with Ref. [30] that we have followed mostly in Sect. 2.2, but was, for example, done in a very similar model by Zandt et al. that has been used to study ion dynamics during anoxia [34]. In the classical HH model the leak current takes the role of stabilizing the physiological fixed point. We test if this is also the effect of leak currents for ion dynamics, and therefore simply compare models that do and do not contain the chloride leak.

There is also a certain freedom in the choice of a pump model. It is a general feature of Na+/K+–pumps that their activity is enhanced by the elevation of ECS potassium and ICS sodium. Still different models exist, and to investigate the role of the particular pump model we replace the pump by Cressman et al. from Eq. (2.30), now referred to