Mathematical approaches to modeling of cortical spreading depression.

Mathematical approaches to modeling of cortical spreading depression Robert M. Miura, Huaxiong Huang, and Jonathan J. Wylie Citation: Chaos: An Interdisciplinary Journal of Nonlinear Science 23, 046103 (2013); doi: 10.1063/1.4821955 View online: http://dx.doi.org/10.1063/1.4821955 View Table of Contents: http://scitation.aip.org/content/aip/journal/chaos/23/4?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Transcranial magnetic stimulation of mouse brain using high-resolution anatomical models J. Appl. Phys. 115, 17B303 (2014); 10.1063/1.4862217 Modeling oscillatory dynamics in brain microcircuits as a way to help uncover neurological disease mechanisms: A proposal Chaos 23, 046108 (2013); 10.1063/1.4829620 Randomness switches the dynamics in a biophysical model for Parkinson Disease AIP Conf. Proc. 1479, 1434 (2012); 10.1063/1.4756429 Mathematical models of brain dynamics: Progress and perspectives AIP Conf. Proc. 1450, 18 (2012); 10.1063/1.4724111 Stochastic diffusion model of multistep activation in a voltage-dependent K channel J. Chem. Phys. 132, 145101 (2010); 10.1063/1.3368602

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 138.251.14.35 On: Sun, 21 Dec 2014 00:38:28

CHAOS 23, 046103 (2013)

Mathematical approaches to modeling of cortical spreading depression Robert M. Miura,1,a) Huaxiong Huang,2,b) and Jonathan J. Wylie3,c) 1

Department of Mathematical Sciences and Center for Applied Mathematics and Statistics, New Jersey Institute of Technology, Newark, New Jersey 07102 USA 2 Department of Mathematics and Statistics, York University, Toronto, Ontario M3J 1P3, Canada 3 Department of Mathematics, City University of Hong Kong, Kowloon, Hong Kong

(Received 16 April 2013; accepted 9 September 2013; published online 10 October 2013) Migraine with aura (MwA) is a debilitating disease that afflicts about 25%–30% of migraine sufferers. During MwA, a visual illusion propagates in the visual field, then disappears, and is followed by a sustained headache. MwA was conjectured by Lashley to be related to some neurological phenomenon. A few years later, Le~ao observed electrophysiological waves in the brain that are now known as cortical spreading depression (CSD). CSD waves were soon conjectured to be the neurological phenomenon underlying MwA that had been suggested by Lashley. However, the confirmation of the link between MwA and CSD was not made until 2001 by Hadjikhani et al. [Proc. Natl. Acad. Sci. U.S.A. 98, 4687–4692 (2001)] using functional MRI techniques. Despite the fact that CSD has been studied continuously since its discovery in 1944, our detailed understandings of the interactions between the mechanisms underlying CSD waves have remained elusive. The connection between MwA and CSD makes the understanding of CSD even more compelling and urgent. In addition to all of the information gleaned from the many experimental studies on CSD since its discovery, mathematical modeling studies provide a general and in some sense more precise alternative method for exploring a variety of mechanisms, which may be important to develop a comprehensive picture of the diverse mechanisms leading to CSD wave instigation and propagation. Some of the mechanisms that are believed to be important include ion diffusion, membrane ionic currents, osmotic effects, spatial buffering, neurotransmitter substances, gap junctions, metabolic pumps, and synaptic connections. Discrete and continuum models of CSD consist of coupled nonlinear differential equations for the ion concentrations. In this review of the current quantitative understanding of CSD, we focus on these modeling paradigms and various C 2013 AIP Publishing LLC. mechanisms that are felt to be important for CSD. V [http://dx.doi.org/10.1063/1.4821955]

Cortical spreading depression (CSD) has been the focus of intensive experimental and theoretical study since its discovery by Le~ ao in 1944. The connection between CSD waves and migraine with aura (MwA) has made the mathematical study of CSD even more compelling. The quantitative study of CSD has led to mathematical models that incorporate many of the mechanisms known to affect CSD speed, shape, and duration. CSD involves the movement of the major ions (potassium, sodium, chloride, and calcium) in the brain. These ions are transported by diffusion and can go through cell membranes as a result of ionic membrane channels and energy consuming pumps that restore the ionic concentrations to homeostasis. The massive variations in ion concentrations associated with CSD waves also give rise to significant osmotic effects and puts intense stress on the neurovascular coupling that regulates oxygen supply to the brain. The mathematical models of CSD described here vary in their levels of detail and the objectives they aim to achieve. A definitive model has yet to be formulated.

a)

E-mail: [email protected] E-mail: [email protected] c) E-mail: [email protected] b)

1054-1500/2013/23(4)/046103/12/$30.00

I. INTRODUCTION

Chronic migraine is a genetic neurological disease and is a debilitating headache disorder that is ranked number 19 by the World Health Organization among all diseases worldwide causing disability. It is estimated to have a global prevalence of about 14%.53 Also, it has been found to occur in women three times more often than in men.5 MwA occurs in approximately 30% of people35 who suffer from migraine and the aura (causing visual, sensory, and/or speech symptoms) either precedes or occurs during the headache. In the visual case, the effects of the migraine headache are compounded because of the illusions in the visual field that are created within the brain from the visual cortex. These visual effects vary with multifarious forms: blind spots, fortification patterns, moving patterns, etc. (For more detailed discussion of various forms of migraine and especially migraine with aura, the readers are referred to the chapter by Dahlem11 in this issue.) Although these auras and their connections with migraine have been recognized as a clinical condition for over three centuries,16 it was not until 1941 that Lashley34 conjectured that visual aura was connected with a propagating cerebral phenomenon. In 1944, almost 70 years ago, Le~ao from Brazil, then a physiology doctoral student

23, 046103-1

C 2013 AIP Publishing LLC V


046103-2

Miura, Huang, and Wylie

studying at the Harvard Medical School, was trying to understand the “cortical electrogram” data collected from experiments on epilepsy in rabbits.38 He noticed that sometimes the epileptic wave would be followed by a slow wave emanating from the location of the initial stimulus that reduced the neuronal excitability in the cortex of the brain.38 Le~ao initially described the phenomenon as a “spreading depression” and it is now referred to as “CSD” and identifies the depression and spread of the electroencephalogram record. It had been conjectured for many years that the propagating cerebral phenomenon referred to by Lashley34 was CSD. However, during the aura, reliable data confirming this link between CSD and MwA did not appear until 2001 in the work of Hadjikhani et al.,24 who measured blood oxygenation level-dependent signal changes with functional MRI which showed many of the characteristics of CSD. Neurophysiological experiments in vivo and in vitro have provided an enormous amount of data on CSD, and they have delineated many of the mechanisms involved. However, it has proved to be an onerous task to determine the detailed causal effects of CSD on MwA given the difficulty in attempting to capture the aura events when they occur in humans24 and having the lack of access to animal observations of the aura. In this survey, we focus our attention on the mathematical modeling of CSD. A. Cortical spreading depression (CSD)

Waves of CSD are slow chemical waves that have been observed in the cortices of many brain structures in laboratory animals. These waves reduce excitability of neuronal tissue and lead to massive changes in some of the extracellular and intracellular ionic concentrations.6,39,50,51 These ions include sodium, potassium, chloride, and calcium; the latter being an important signaling ion and a second messenger in the cytosol. Waves of CSD can be instigated using a variety of stimuli, including applications of potassium chloride (KCl) on the cortical surface, mechanical impact, electrical stimulation, as well as other means.6,50 However, some of the precise mechanisms involved in the instigation of CSD waves under these various stimuli are not involved in the propagation of CSD waves. The focus in this review is on mathematical models concerned mainly with propagation of CSD waves. CSD waves also lead to shrinkage of the extracellular space (ECS) and increased metabolism that depletes local oxygen in brain tissue. Most studies of CSD have focussed on the propagation of CSD waves, which are extremely slow with speeds of 1–10 mm/min. Since diffusion coefficients of the various ions are assumed to be constant, an explanation of the variability in the speed of CSD involves tissue cell types and structures, volume fraction of the extracellular space, and membrane ionic channel distributions and currents, as well as other possible mechanisms. CSD wave propagation basically occurs in the cortex of a brain structure (effectively a two-dimensional sheet of grey matter in three dimensions). However, a CSD wave, assumed to be propagating in one dimension, has many qualitative properties in common with action potential propagation in nerve axons. Both waves appear to be solitary waves,

Chaos 23, 046103 (2013)

although multiple waves are seen experimentally. These waves are all or none, i.e., either a stimulus generates a sustained wave or any initially generated disturbance eventually decays. They exhibit refractoriness, i.e., a line stimulus parallel to and following a one-dimensional CSD wavefront will not generate a second wave if it is too close to the refractory part of the wave. On the other hand, numerical experiments have been carried out whereby a line stimulus following a wave can be located essentially where the tissue is refractory on one side of the stimulus and no longer refractory on the other side. This will result in a single wave front propagating in a direction opposite to the first wave. Finally, two onedimensional CSD waves will annihilate each other when they collide head on. (For a discussion of two-dimensional wave patterns of CSD waves and their clinical importance, see Dahlem12). What distinguishes CSD waves from action potentials, however, is the difference in time scale, namely, CSD evolves on a time scale of minutes whereas action potentials occur over a time scale of milliseconds. Also, CSD waves propagate at very slow speeds, on the order of mm/min, relative to action potentials, which travel on the order of m/ms. The literature on cortical spreading depression is large, and a number of excellent reviews have been written.6,19,39,42,49,50 There have been numerous studies on ionic and chemical correlates of the phenomenon. Furthermore, the complexity of CSD waves elevates the role of mathematical modelling, through which one can alter, control, and focus on specific mechanisms. There are numerous references in the literature to putative relationships between CSD and a variety of pathological states of the brain, including migraine with aura,2,10,19,23,36,43,44 stroke,37 cerebral ischemia-infarctions,19,20 and transient global amnesia.19 As mentioned above, the relationship between CSD and MwA has been established.24 There also are suggestions of the protective role that CSD could play in reducing ischemic brain injury.33 These are important biomedical correlates of CSD, and each one deserves further study from the modeling point-of-view. Our focus here is to identify the relevant mechanisms that appear involved in CSD and to include them in a comprehensive model. In spite of knowing many of the basic mechanisms involved in CSD, we still do not understand the relative importance of these mechanisms and how they conspire to produce the observed wave phenomena. We can greatly extend our basic knowledge of how the brain functions by understanding what happens and causes the “massive failure of ion homeostasis” during CSD.32 Furthermore, CSD serves as a paradigm for our understanding of the basic mechanisms in the brain that are important in the study of neuroscience. It is extremely difficult to perform detailed experiments that record the multifarious variables that are changing during CSD. Therefore, a good understanding of the complex interactions between these mechanisms can only come from a detailed mathematical modeling viewpoint where one can turn mechanisms on and off. Cortical spreading depression remains an enigma today, and identification of how the precise mechanisms involved lead to the instigation and propagation of CSD waves has remained inconclusive. Cortical


046103-3


spreading depression is a wonderful research tool to help us understand basic neurophysiological mechanisms for the normal functioning brain. The study of CSD continues to challenge experimental and mathematical neuroscientists in the 21st century. B. Mathematical models of CSD

In spite of all the experimental data and observations, the precise synthesis of these data and putative mechanisms into a coherent picture of the instigation and propagation of CSD waves remains difficult to capture. Additional tools that have been used are mathematical models that provide a precise and generalized method for exploring a variety of mechanisms that may be important for CSD wave instigation and propagation. In this survey of mathematical approaches that have been used for understanding CSD, we describe a synthesis of discrete and continuum modeling methods combined with various mechanisms that can be applied to CSD. We note that the earliest attempts at modeling used cellular automata representations based on cardiac modeling46,59 or on simplified neuronal models, see Grafstein21,22 that utilized the simplified Hodgkin-Huxley equations,26 e.g., the FitzHugh-Nagumo model.15,22,40 By discrete modeling of CSD, we mean models that treat each of the individual neurons as discrete entities in space and then couple them together electrically and chemically. A single neuron model has been studied extensively in the context of CSD by Kager et al.28,29 using the computer program, NEURON.25 Simplified versions of these neurons, e.g., only soma or soma and dendrites, were coupled in a network by Huang et al.27 to simulate the instigation and propagation of CSD waves. By continuum modeling, we mean models in which the temporal evolution and spatial dependence of the dependent variables are continuous. Thus, these models are naturally described by ordinary and partial differential equations. These models come from basic principles and are derived by taking suitable averages over the ECS and the intracellular space (ICS) of two cell types, namely, neurons and glia. The concept of the modeling in this case is to treat each of the spaces, ECS and ICS, as continuum spaces, i.e., the independent variables of time and space and the dependent variables are all continuously varying quantities.

Chaos 23, 046103 (2013)

CSD. As a result, the diffusion of ions plays an important role in the modeling of CSD.1 In fact, in the 1950s, Grafstein21,22 postulated the importance of extracellular potassium and its diffusion in CSD. There are two aspects of diffusion that are important in the modeling of CSD. In the ECS, the diffusion of ions is, in general, more restrictive than that in free space due to the presence of cells and in some cases, the vascular network. At the scale of CSD, this slowing of the effective diffusive spreading can be modeled by introducing a tortuosity, which can be estimated using simple mathematical models. For example, Dai and Miura14 used a lattice Boltzmann method to estimate the value of tortuosity by carrying out numerical simulations of a diffusive substance spreading through a cellular medium. In the brain-cell-microenvironment, the value of tortuosity is estimated to be around 1.6. In other words, the diffusion in a tortuous space such as the extracellular space is slowed by a factor of 1.6. On the length scale of CSD, the diffusion within the ICS is negligible inside the neurons. However, for glial cells, diffusion in the ICS is significant since glia often form a connected syncytial net which allows ions to move within the network. Therefore, an effective diffusion coefficient can be introduced to capture the network nature of the intracellular transport of ions inside the glia syncytium via gap junctions. A typical value is given in Chen and Nicholson.9 B. Membrane currents and pumps

The membrane currents during CSD are quite different from those under normal conditions as the ions channels are operating outside their normal range. In addition, not all ion channels are important and necessary for the massive fluxes of ions that lead to the depolarization of membrane potential. However, identifying the channels that are important for CSD proves to be a non-trivial task. In addition to the passive membrane channels, there are energy consuming pumps that move ions against their electrodiffusive gradients. These pumps consume a significant amount of energy and are essential to maintain brain homeostasis. During CSD, the pumps are put under much greater stress than in homeostasis as they try to restore the resting ion concentrations and membrane potential. 1. Ions channels

II. MECHANISMS INVOLVED IN CSD

One of the main objectives in mathematically modeling CSD is to determine which of the multitude of mechanisms involved in CSD are essential for the phenomenon to occur, and to what extent they interact to produce the emergent phenomenon of CSD wave propagation. Here, we list some of the physical and physiological mechanisms that are believed to be important for the occurrence of CSD and its various properties. A. Ion diffusion

Relative to the normally functioning brain, there is a massive redistribution of ions in the ECS and ICS during

The most relevant ion channels for CSD have been identified while some otherwise important channels during normal activity are found to be not needed for the instigation and spreading of CSD. For example, the fast sodium channels, which are essential for generating action potentials, are not important for the instigation and spreading of CSD. Both experimental and theoretical studies have confirmed that blocking these channels does not prevent CSD.50,54,56 On the other hand, it is generally believed that NMDA (N-methyl-Daspartate) channels are of critical importance due to the observation that blocking NMDA receptors can prevent CSD. A simple mathematical model that describes these channel currents for each ion is the Hodgkin-Huxley (HH) formula


046103-4


Chaos 23, 046103 (2013)

Iion;HH ¼ gion;HH ðEm Eion Þ;

(1)

where gion,HH is the conductance for each ion (sodium, potassium, chloride, and other ions), Em is the membrane potential, and Eion is the Nernst potential for each ion. A more detailed model (that includes the experimentally observed phenomena that channels can be activated and deactivated and that the channel currents depend on the ion concentrations in the ECS and ICS) is given by the GoldmanHodgkin-Katz (GHK) formula Em FPion;GHK Em ½ioni exp ½ione / ; Iion;GHK ¼ mpion nqion Em / 1 exp / (2) where mion and nion are the activation and inactivation gating variables, respectively, p and q are integers that depend on the channel type, F is the Faraday constant, / ¼ RT=F, and Pion,GHK is the effective permeability. Here, R and T are the universal gas constant and absolute temperature. The variables [ion]i and [ion]e are the ICS and ECS ionic concentrations, respectively. For a more detailed description, we refer the readers to Bennett et al.4, Huang et al.27, and Yao et al.60 2. Ion pumps

There are many different pumps and exchangers. One of the most important pumps is the potassium-sodium exchange pump that moves two potassium ions into a cell for every three sodium ions removed from the cell. The mathematical form of this pump current has been measured in experiments and a good empirical approximation is given by an expression with Hill-type dependence on extracellular potassium and intracellular sodium concentrations, namely, IK;pump ¼ rNaK

½Ke ½Ke þ ae

!2

½Nai ½Nai þ ai

3 INa;pump ¼ IK;pump ; 2

!3 f ðEm Þ;

(3)

(4)

where [K]e and [Na]i are the ECS potassium and ICS sodium concentrations, respectively; rNaK, ae, and ai are constants; and f(Em) is a function of the membrane potential Em. The values of the parameters and the exact functional form for f can be found in the supplementary material of Ref. 4, along with other types of pumps and exchangers. C. Spatial buffering

First introduced by Kuffler and colleagues in their studies on Necturus glial cells,45 spatial buffering is the direct consequence of basic biophysical properties. In a glial syncytium connected via gap junctions, the membrane of connected cells tends to remain isopotential. An increased extracellular potassium concentration causes a local depolarization across the glial cell membrane, which by the cable

properties of the membrane, spreads electrotonically down the cells to more remote regions of the syncytium. Local currents are formed by the unbalanced spatial distribution of membrane potential differences and completed through the ECS. Since glial cell membranes exhibit a high potassium conductance, the asymmetric potential differences result in an influx of ECS sodium at the site of local depolarization and an efflux of potassium into the ECS at more distal regions where its concentration is still near resting levels. The potassium is then passively transported from an ECS location of high concentration and dispersed to the regions with low concentration. The circuit current loop is closed by intracellular and extracellular ionic current flows, primarily those of sodium and chloride. It has been shown that this passive buffering mechanism is energy-independent and generally more efficient than diffusion for transporting potassium through the interstitial space.17,18 Mathematical models have been constructed to study spatial buffering, using either a continuum approach in one space dimension9 or a discrete lattice-Boltzmann technique in a two-dimensional setting.52 D. Osmosis and cell swelling

Under normal operating conditions, the cell volume is regulated and the brain-cell-microenvironment is under osmotic balance that maintains isotonicity in the ICS and ECS. During CSD, on the other hand, the isotonicity is no longer maintained and the massive redistribution of ions induces an osmotic pressure that drives water into the cells. The osmotic pressure is given in Keener and Sneyd,31 Eq. (2.83) p ¼ RTwNc;

(5)

where w is the osmotic coefficient, which indicates the deviation of the solvent from ideal behavior, c is the solute molar concentration, and N is the number of ions the compound dissociates into. When there is solute on both sides of the membrane, the difference of osmotic pressures provides the driving force for the water flux across the membrane. When CSD passes through a region of the brain, the influx of the ions is greater than the efflux. As a result, the osmotic pressure is greater in the ECS. This drives water into the cells and they swell. According to Shapiro,47 the rate of change of ECS volume fraction, f^, from the flow of water across the membrane due to osmotic pressure differences in a single ion can be modeled as Pf VW S @ f^ ¼ ð½ioni ½ione Þ; @t V

(6)

with f^ ¼ Ve/V, where Ve is the volume of the ECS, Vi is the volume of the ICS, V ¼ Ve þ Vi is the total fixed volume, Pf is the osmotic water permeability of the membrane, VW is the partial molar volume of water, and S is the surface area of the membrane. In the current setting, the osmotic pressure differences due to the three ions combine additively and the relationship is modified by assuming that the imbalance is offset by immobile ions that cannot cross the membrane. The new equation is given by


046103-5


Pf VW S @ f^ ¼ ½Nai þ ½Ki ½Cli ½Nae @t V Vi ½Ke þ ½Cle þ 0 A ; Vi

Chaos 23, 046103 (2013)

III. MATHEMATICAL MODELS FOR CSD

(7)

where Vi,0 is the resting value of Vi and A ¼ [Na]e,0 þ [K]e,0 [Cl]e,0 [Na]i,0 [K]i,0 þ [Cl]i,0 represents the sum of the resting values of the ionic concentrations that add up to the concentration of the immobile anions trapped inside the cells.

Historically, there have been several attempts to model CSD from a mathematical point of view. Above, we mentioned earlier usage of mathematical models to study CSD. Note that there are two aspects of CSD that should be treated differently. They are instigation of CSD and propagation of CSD. The first attempt to incorporate detailed physiological mechanisms into CSD propagation models was due to Tuckwell and Miura.56 A. Tuckwell-Miura model

E. Other mechanisms 1. Glutamate

Glutamate as a possible agent for CSD was proposed by van Harreveld in 1959.57 In 1978, he proposed a dual hypothesis with two separate mechanisms for CSD.58 The basic hypothesis is that the release of ATP from glia (astrocytes), which propagates independently of the calcium wave, leads to the release of glutamate from astrocytes. This glutamate acts on the NMDA receptors of neurons to trigger a large depolarization of the membrane potential (via massive redistribution of ions) and release of glutamate from the neurons, which leads to additional releases of ATP from astrocytes that complete the cycle.4 While there is strong evidence in favor of the glutamate hypothesis, other experimental work also supports the potassium mechanism even though the details of the original theory may not be correct.50 It is very likely that both glutamate and potassium are important and there is evidence in favor of the dual hypothesis. For a detailed discussion, we refer the readers to the book by Somjen.50 2. Vascular effects

CSD is associated with profound vascular changes that may be a significant factor in the clinical response to CSD events such as migraine. Recent experimental studies7 using a combination of optical intrinsic signal imaging, electrophysiology, potassium sensitive electrodes, and spectroscopy have revealed a more detailed picture of the vascular response to CSD. Two distinct phases of altered neurovascular function have been identified, one during the propagating CSD wave and a second much longer phase after the passage of the wave. The direct current shift associated with the wave was accompanied by marked arterial constriction and desaturation of cortical haemoglobin. After recovery from the initial wave, a second phase of prolonged, negative direct current shift, arterial constriction, and haemoglobin desaturation occurred. Persistent disruption of neurovascular coupling was demonstrated by a loss of coherence between electrophysiological activity and perfusion. ECS potassium concentration increased during the CSD wave, but recovered and remained at baseline after passage of the wave, consistent with different mechanisms underlying the first and second phases of neurovascular dysfunction. The interactions between the vascular effects and CSD have just been modeled for the first time.8

Preliminary studies by Sugaya et al.54 of CSD ionic mechanisms showed that application of tetrodotoxin (TTX) to a region of brain tissue, which eliminates neuronal action potentials, did not prevent the propagation of CSD through the treated region. Thus, the need for action potentials, which had been assumed to cause the increase of ECS potassium, was unnecessary, i.e., the fast transient sodium membrane currents could be eliminated from the model. This led Tuckwell and Miura56 to postulate the neurotransmitter hypothesis. In particular, increased ECS potassium depolarizes surrounding cell membranes. As a result, presynaptic boutons will depolarize and allow calcium to enter and cause exocytosis of neurotransmitter molecules from vesicles. The neurotransmitter molecules attach to potassium channel receptors and open up potassium channels allowing further release of potassium into the ECS. The novel modeling aspect of the Tuckwell-Miura model was the way in which they treated brain tissue mathematically. The difficulty with brain tissue is that neurons are isolated and not contiguous. A propagating CSD wave takes on the order of a minute to traverse a spatial position in the cortex, thus, the capacitive currents can be neglected for the bulk of the CSD wave. Furthermore, the space scale of the CSD wave is long compared to cell soma size. This means that the local variables, e.g., ionic concentrations, vary little as one moves from one location in the ECS to a nearby location in the ECS. The same holds for ion concentrations in the ICS of neighboring neurons or glial cells. Note that a CSD wave travels in the cortex parallel to the cortex surface and perpendicular to the general direction of the dendritic trees of neurons (cf. Fig. 1). The mathematical model created by Tuckwell and Miura56 treats the tissue as overlapping continuum spaces, namely, the ECS occupies the same physical space as the ICS for the neurons and glia cells (cf. Fig. 2). Because of the near uniform values of different variables in the ECS or in the different ICS for neurons and glial cells, if an ionselective microelectrode is inserted into the tissue at a given location, and the tip measurement records the ECS value of the ionic concentration, then that is considered the ECS measurement at that location. On the other hand, a small movement of the microelectrode could result in the tip measurement recording the ICS value of the ion concentration in a neuron, so that is considered the ICS measurement at that point in the ICS of the continuum neuronal space. Mathematically, the ECS is connected and ions can move unimpeded. The actual ICS for the neurons and glia cells are


046103-6


Chaos 23, 046103 (2013)

FIG. 1. Schematic of the cortex illustrating the direction of propagation of CSD waves.

not connected, but the ICS is treated mathematically as a connected continuum space as well. What distinguishes the ICS from the ECS is that ions do not diffuse in the ICS. In order for an ICS neuronal ion to move, it must first cross the membrane into the ECS, then diffuse and then enter a neighboring neuron. This concept of overlapping continuum spaces allows the model to consist of a system of coupled differential equations with no interior boundaries. For the Tuckwell-Miura56 model, which represents a minimal continuum model in one dimension, only two ions were taken into account, potassium, the most important ECS ion, and calcium, responsible for release of neurotransmitter molecules. The corresponding ECS and ICS model equations account for diffusion of ECS potassium and calcium and for mechanisms at the membrane level, including diffusion of ions, ionic membrane currents, and electrogenic pumps, namely, @½Ke @ 2 ½Ke ¼ DK þ q1 ðIK þ IK;pump Þ; @t @x2

(8)

FIG. 2. Simplified continuum model of CSD propagation through the cortical tissue.

@½Ki a ¼ q ðIK þ IK;pump Þ; @t 1a 1

(9)

@½Cae @ 2 ½Cae ¼ DCa þ q2 ðICa þ ICa;pump Þ; @t @x2

(10)

@½Cai a ¼ þ q2 ðICa þ ICa;pump Þ; @t 1a

(11)

where [K]e,[Ca]e and [K]i, [Ca]i are the ionic concentrations of potassium and calcium in the ECS and ICS, respectively, DK and DCa are the ion diffusion coefficients for potassium and calcium, respectively, q1 and q2 are constant coefficients signifying channel distributions and cell membrane area, IK, IK,pump and ICa, ICa,pump correspond to the membrane ionic currents as modeled by Eq. (1) with voltage and time dependent conductances and the energy-consuming pumps for potassium and calcium, respectively, and a is the fraction of the total volume occupied by the ECS. We note that there may be small variations of the effective diffusion coefficients with concentrations, temperature, and other factors, but since these variations tend to be difficult to measure in the complex environment of the brain, it was assumed that they are well-approximated by constant values. Hodgkin-Huxleytype currents26 with a voltage-calcium-dependent potassium conductance and a voltage-dependent calcium conductance were assumed and pump expressions were chosen to return the ECS potassium and ICS calcium back to their resting values. A result of this model was a wave speed of approximately 1 mm/min (taking into account tortuosity) that is within the range of experimentally measured values. B. Shapiro model

Shapiro47,48 proposed a model that is in the same spirit as the Tuckwell-Miura model in the sense that it is a continuum model that describes ionic concentrations in compartments for the ECS and ICS that are overlapping everywhere. Following Tuckwell-Miura, the model consists of a set of coupled reaction-diffusion equations that describe the fluxes of ions. However, whereas the Tuckwell-Miura model can be viewed as a minimal continuum model for CSD wave propagation that accounts for the mechanisms at the cellular level, the Shapiro model strives for a much more detailed picture of CSD propagation. The model includes two major features that are not included in the Tuckwell-Miura model. First, Shapiro hypothesized that diffusion of ions in the ICS through neuronal gap junctions is important in the propagation of CSD. Due to the absence of detailed measurements of the gating properties of neuronal gap junctions, he adopted an ad hoc diffusion model (motivated by diffusion observed in cardiac cells) with a diffusion coefficient that depends on the voltage drop across gap junctions and ionic concentrations of protons and calcium. Second, he assumed that the extremely high ionic concentrations arising when CSD propagates will give rise to significant osmotic forces that drive cell swelling as described in Sec. II D. In addition, Shapiro’s model provides a much more detailed view of the ionic fluxes that occur during CSD. It contains 20 different channels and pump currents and describes the concentrations of


046103-7


Chaos 23, 046103 (2013)

nine ionic species. Eight of these are given by potassium, sodium, chloride, and calcium in the ECS and ICS, respectively, and the ninth is buffered calcium. Each of these ionic concentrations is supplemented by Hodgkin-Huxley activation and inactivation variables. The model also includes IP3, inositol triphosphate, in the ECS and is used in signal transduction within the cells where it is assumed to be important in activating calcium channels on the endoplasmic reticulum membrane that can act as a calcium buffer. IP3 production was assumed to be induced by a membrane stretching mechanism that was modeled in a simple way. IP3 was also assumed to degrade following the model of Keener and Sneyd.31 Surprisingly, the results from the model showed that diffusion in the ECS was not necessary for CSD wave propagation. However, the effective diffusion of ions through the gap junctions provides an alternative way to transport ions and it was found that blocking the gap junctions prevents CSD waves from propagating. Recent experiments with carbenoxolone (CBX), known to block gap junctions, showed that CBX can accelerate the initiation and propagation of CSD,55 which shows that blocking gap junctions will not prevent CSD. The role of osmotic cell swelling was also examined. If cell swelling was not included in the model, then CSD waves were found not to propagate. However, other modeling approaches60 found that neglecting the effects of cell swelling in the model did not affect CSD wave propagation significantly. C. Kager-Wadman-Somjen model

In Refs. 28 and 29, Somjen and collaborators investigated CSD-like membrane depolarizations of a single neuron using the computer model NEURON.25 In their model, an isolated neuron is immersed in an ECS, forming a system of two compartments. Cross-membrane currents of the major ions (Na, K, and Cl, but not Ca) are included. For leak currents, they used the HH formula (1). For other currents such as the fast sodium currents or NMDA currents, they used the GHK formula (2). For energetic pump currents, they focus on the Na/K exchange pump where the currents are given by Eqs. (3) and (4). Both ICS and ECS ion concentrations are allowed to change due to the cross-membrane ion fluxes d½ioni Iion Scell ¼ ; dt FVi

(12)

d½ione Iion Scell ¼ ; dt FVe

(13)

where Iion includes the cross-membrane current of the ion, Scell is the cell surface area, and Vi and Ve are the ICS and ECS volumes associated with this cell. Additional features of the model are a detailed treatment using NEURON of the dendritic structure and ion movement in different parts of the dendrites due to concentration differences. The control of ECS potassium accumulation is handled by a buffer. Osmotic forces are also considered. In a more recent paper by the same authors,30 the glial compartment and calcium are added to examine the effects of buffering of potassium in more detail.

The NEURON simulations show that CSD-like depolarization can be induced either by an injected current over a short period of time or by shutting down the Na/K exchange pump in the neuronal membrane (but not in the glial cell membrane). Even though not all the cross-membrane ion currents are included in the model, essential features of the CSD time course are captured. On the other hand, the propagation of CSD cannot be investigated using a single neuron model. D. Bennett-Farnell-Gibson model

In Ref. 4, it was proposed that glutamate released by astrocytes is the main agent responsible for the propagation of CSD, along with the ATP released by neurons. The model consists of a one-dimensional array of astrocytes and neurons, which are coupled via a two-way communication. The astrocyte-to-neuron communication is through the NMDA receptors on the neuronal membranes, which are activated by the glutamate released by astrocytes. On the other hand, the neuron-to-astrocyte communication goes through the metabotropic receptors on the astrocytes, which trigger the release of ATP from the astrocytes. The ATP then diffuses through the ECS. The concentrations of the major ions and the membrane potential undergo changes due to changes in ion channels including those activated by NMDA receptors. In this model, NMDA currents are the key currents and the gating variables for NMDA currents are modeled as being glutamate-sensitive but do not depend on potassium explicitly Em ½ione gion;GHK FEm ½ioni exp / ; INMDA;ion ¼ Popen Em / 1 exp / (14) where / ¼ RT/F and Popen is the glutamate dependent channel open probability, given by dPopen ¼ r1 ½GluA ð1 Popen Þ r2 Popen ; dt

(15)

where r1 and r2 are constants and [Glu]A is the concentration of glutamate released by astrocytes. The production of glutamate results from several biochemical pathways that involve astrocytes and neurons, stimulated by ATP as @½GluA VA ½ATP0:98 ¼ kA ½GluA ; @t 7 þ ½ATP0:98

(16)

where VA and kA are constants and [ATP] is the local ATP concentration (released by astrocytes). The first term on the right-hand side represents the production of glutamate due to ATP and the second term represents the loss of free glutamate due to absorption, e.g., by the receptors on the neuronal membranes. The concentration of glutamate released by neurons [Glu]N is modeled by


046103-8

Miura, Huang, and Wylie 2 @½GluN ¼ VN e0:0044ðEm 8:66Þ kN ½GluN ; @t

Chaos 23, 046103 (2013)

(17)

where VN and kN are constants. Note that the release of glutamate from neurons depends on the membrane potential while that by astrocytes depends on ATP. The mechanism of ATP release from astrocytes has not been established. It was proposed in Ref. 3 that IP3 acts as the intermediate agent that triggers the release of ATP into the ECS. The intermediate compound IP3 is released from the internal stores of astrocytes but its concentration decays due to degradation and spatial diffusion (via gap-junctions within the astrocyte network) d½IP3 j dt

¼ rh ðGATP þ GGlu Þj kdeg ½IP3 j þ cð½IP3 j1 ½IP3 j Þ;

(18)

where the concentration is measured at the neighbouring boundaries of astrocytes (j 1 and j) and rh, kdeg, and c are constants. The fraction of G-protein activated by ATP is given by GATP ¼

qATP þ dAT ½ATP ; ; qATP ¼ KATP þ qATP þ dATP KRAT P þ ½ATP (19)

where dATP is a constant. The additional fraction of G-protein can be activated by glutamate released by neurons as GGlu ¼

½GluN qGlu ; qGlu ¼ : KGlu þ qGlu KRGlu þ ½GluN

“localized transient” waves. Dahlem and Isele13 used a generic reaction-diffusion model that generalized the model by Grafstein, Hodgkin, and Huxley @u u3 ¼ u v þ r2 u; 3 @t ð @v ¼ u þ b þ K HðuÞdxdy: @t e

(22) (23)

Here, u (activator) and v (inhibitor) are both variables that qualitatively represent all of the excitatory and inhibitory processes, H is the Heaviside function, K (not potassium) is the strength of the mean field inhibitory coupling (representing the vascular feedback), and b is a control parameter that determines the dynamics of the model. Without the mean field inhibitory coupling, i.e., K ¼ 0, the system is bistable. The homogeneous state u ¼ 0 represents the normal resting state. The non-homogeneous state u > 0 (associated with CSD) is also stable. When a critical mass is reached, termed “critical nucleation,” CSD is ignited and it engulfs all the other regions (tissue) after initial ignition. When K > 0, this critical nucleation is controlled by a widely spread inhibitory feedback. The travelling wave is no longer global and becomes localized, and in many cases, it disappears completely. From the dynamical systems point of view, this represents a saddle-node bifurcation, which can be used to explain the transient nature of the CSD wave associated with migraine. For more detailed discussion, the readers are referred to the paper by Dahlem in this issue.11

(20) F. Huang-Miura-Yao model

Finally, the extracellular potential VSD is computed using the following equation: 1 @½ione X @VSD RT X ¼ zion Dion zion Dion ½ione (21) @x @x F ion ion where Zion is the valence of each ion and the extracellular ion concentration [ion]e is calculated using interpolation since diffusion is not included in the model. E. Dahlem model

While the models we have discussed so far attempted to model individual ions and mediators, such as glutamate and ATP, and their collective behavior using the physiologically based conductance for membrane currents, phenomenological mean field models have been put forward by Dahlem and his collaborators. Their main focus is to provide an understanding of the role played by CSD in migraine from the dynamical systems point of view. The key observation is that CSD associated with migraine is not a regular wave of depolarization that propagates outward along the cortex from the location of initiation. Instead, these CSD waves pop up and die out in various locations in the cortex as revealed by fMRI images of patients undergoing migraine attacks.24 In other words, these are

Huang, Miura and Yao27 constructed a simplified neuronal model by extending the work of Kager et al.28 They simplified the NEURON model used in Kager et al.28 by replacing the soma-dendrite system with a point neuron while retaining the basic cross-membrane currents for sodium and potassium ions. The equation for the membrane potential is given by dEm ¼ I; (24) Cm dt where t is the time, Cm is the membrane capacitance per unit surface area, and I is the total cross-membrane ionic current per unit surface area. The current due to sodium ions is given by INa ¼ INa;T þ INa;P þ INa;Leak þ INa;Pump where INa;T is the transient sodium current and represents the fast sodium current involved in action potentials, INa;P is the persistent sodium current, INa;Leak is the sodium leak current, and INa;Pump is the sodium exchange pump current, responsible for restoring the ionic sodium concentration back to its homeostatic state. Also involved in the total cross-membrane current is the potassium current, IK ¼ IK;DR þ IK;A þ IK;Leak þ IK;Pump , where IK;DR is the delayed rectifier potassium current, IK;A is the potassium A current, IK;Leak is the potassium leak current, and IK;Pump is the potassium exchange pump current. The time evolution of the concentrations of potassium and sodium ions are determined by the cross-membrane currents as follows:


046103-9


Chaos 23, 046103 (2013)

d½ioni S ¼ Iion ; dt FVi

(25)

d½ione S ¼ Iion ; dt FVe

(26)

where ion ¼ Na and K and S, Vi, and Ve are the soma surface area and intra- and extracellular volumes, respectively. Their simplified point neuron model reproduces the essential features of the more complex NEURON model and the CSD-like membrane depolarization shown in Fig. 3. The spreading of the membrane depolarization is modeled by a small network of simplified point neurons. These neurons are considered to be in a one-dimensional array, with a properly defined spacing between them so that the volume fractions for the ICS and ECS are consistent with the values found in the brain. The cross-membrane currents remain the same and the sodium and potassium ions move in the ECS due to molecular diffusion. The governing equations are given by

FIG. 3. CSD-like evolution of the membrane and Nernst potentials with 200 ms electrical stimulation: (a) the first 300 ms and (b) the first 30 s. Also plotted are ICS and ECS concentrations of potassium and sodium and their Nernst potentials.

d½Nae;j dt d½Ke;j dt

¼

¼

S INa;j þ cNa ð½Nae;jþ1 þ ½Nae;j1 2½Nae;j Þ; FVe (27) S IK;j þ cK ð½Ke;jþ1 þ ½Ke;j1 2½Ke;j Þ; FVe

(28)

for neuron j and the ECS associated with it. Here, cNa and cK are coefficients that are obtained from the molecular diffusion coefficients, DNa and DK, scaled by the square of the average distance (d) between the neurons. In other words, cion ¼ Dion/d2.47 For the membrane potential and intracellular concentrations, the equations are given by Eqs. (24) and (25) evaluated at each neuron. When an electrical current stimulus is applied for 200 ms to three neurons, all three neurons depolarize after a prolonged period of action potentials. This triggers a significant buildup of the ECS potassium concentration and influx of sodium ions. As a consequence of the diffusion of potassium in the ECS to the neighboring neurons, action potentials are induced in the neighboring neurons, resulting in a cascade of depolarization of the neuronal membrane potentials one after the other. This represents a spreading of CSD-like depolarization, cf. Fig. 4. The small network of simplified neurons was used to explore the spreading of the membrane depolarization by application of KCl when INa,T is removed. The application of KCl was assumed to be confined to the ECS of the first three neurons. While the first three neurons are depolarized almost instantly after the application of KCl, the other neurons undergo depolarization one-by-one, without the firing of action potentials (results not shown). They also extended their model to investigate the role of NMDA currents. Since the NMDA currents only occur in dendrites, they include the dendrites in their model by adding another compartment. Their simulations indicate that once a wave of CSD has been instigated, it propagates in a similar way to the wave that occurs in the model without dendrites. On the other hand, blocking the NMDA channels prevents CSD instigation and spreading since membrane potentials in the neuronal soma membrane are significantly altered. Most interestingly, if the NMDA channels are blocked in the dendrites of neurons away from the region on which KCl is applied, then the CSD wave is not blocked. We note that as the wave passes through this region, there are changes in the level of depolarizations in the soma. Another important feature of the model is its ability to make quantitative predictions. They have computed the propagation speed of the CSD wave and found that it travels at approximately 8 and 4 mm/min when the fast transient sodium channels are included and removed, respectively. Both are within the range of experimentally observed values. While the discrete models reveal some interesting phenomena, they become computationally inefficient when applied to the large scale phenomena observed in experiments. In these cases, continuum models are favored. In Sec. III G, we review one of the recent models.


046103-10


Chaos 23, 046103 (2013)

@½Nae ½Nae @ f^ 1 @ ^ @½Nae S ¼ ; INa þ f DNa @t @x FV f^ f^ @t f^@x @½Nai S ½Nai @ f^ ¼ INa þ ; @t FVð1 f^Þ 1 f^ @t

FIG. 4. Evolution of the membrane potential and concentration of ECS potassium. In order to instigate a CSD wave, it was necessary to stimulate at least three neurons. The black curves represent the third neuron being stimulated, and the blue and red curves represent neighboring neurons, starting from the nearest ones and alternating in color for ease of visualization. (a) Membrane potentials of individual neurons. (b) Concentrations of potassium in the ECS of neurons.

G. Yao-Huang-Miura model

In Ref. 60, a continuum formulation of the discrete model described in Sec. III F (Ref. 27) was proposed. The length scale of a CSD wave is large compared to the neuron size, thus, they formulated a one-dimensional continuum model with mathematical structure similar to that of Tuckwell and Miura.56 They followed Ref. 56 in assuming that diffusion only occurs in the ECS, but there are several important differences. First of all, their model allows for intracellular and extracellular volume changes. Second, and more importantly, they do not assume any special mechanisms such as the use of neurotransmitters or gap-junctions. Third, a fuller repertoire of membrane ionic currents as well as more realistic ionic pumps, namely, those included in Kager et al.,28,29 were included. The membrane potential Em is governed by Eq. (24) where dEm/dt is replaced by the partial derivative @Em/@t, the concentration of sodium in the ECS and ICS is governed by

(29)

(30)

where f^ ¼ Ve/V is defined after Eq. (6). The first terms on the right-hand sides of Eqs. (29) and (30) are the membrane currents of sodium, which also account for the differences between the ECS and ICS volume fractions. The second terms account for ECS and ICS volume changes. The third term in Eq. (29) is the diffusive flux of sodium in the ECS. In deriving the diffusive flux, one needs to know the ECS fraction of the surface that bounds a control volume. It was assumed that this surface fraction can be well approximated by the ECS volume fraction. This is the reason that the volume fraction, f^, appears in the diffusive flux. Similar equations govern the concentrations of potassium and chloride in the ECS and ICS ð½Ke ; ½Ki ; ½Cle ; ½Cli Þ. During cell swelling, @ f^=@t, given by Eq. (7), is negative and the total intracellular concentration of the ions decreases, while the total extracellular concentration increases. Numerical simulations were carried out to investigate the instigation and spreading of SD waves under different conditions. First, the cell volume and chloride concentrations were held constant in time (no swelling). They showed that the speed of the CSD waves decreases as f^ increases, consistent with an expanded extracellular space leading to a slower temporal buildup of extracellular potassium concentration. The estimated speeds are within the experimentally observed range. Finally, they also investigated the effect of chloride ions and cell swelling by using the full model. As Cl is the only extracellular anion in their simulations, its initial extracellular concentration, [Cl]e,o, was determined by the sum of both cations: [Cl]e,o ¼ [Na]e,o þ [K]e,o, to ensure electroneutrality. They chose the intra-cellular concentration of chloride, [Cl]i, so that isotonicity is achieved. Table I lists the estimated CSD propagation speed under different conditions. It can be seen that including Cl as a variable, similar to sodium and potassium rather than as simply a passive membrane current term, affects the speed of CSD while the effect of cell swelling, which occurs mainly behind the CSD wavefront, is negligible. H. Neurovascular coupling model

As noted at the end of Sec. II E 2, the interactions of the vascular effects and CSD have been modeled for the first time.8 None of the models above for CSD include the TABLE I. CSD propagation speed with different conditions. Condition Without volume change nor C1 With C1 only With volume change only with both C1 and volume change

Speed (mm/min) 9.68 10.12 9.61 10.10


046103-11


important limitations that the vascular system places on the supply of oxygen and other chemicals that are needed by brain tissue during CSD. Energy consuming pumps utilize such chemicals and are responsible for returning the ionic concentrations in the ECS and the ICS back to their homeostatic states. Chang et al.8 (including the authors) focus on modeling important features essential to understanding the implications of neurovascular coupling during CSD. In this model, the sodium-potassium-ATPase, mainly responsible for ionic homeostasis and active during CSD, operates at a rate that is dependent on the supply of oxygen. The supply of oxygen is determined by modeling blood flow through a lumped vascular tree with an effective local vessel radius that is controlled by the extracellular potassium concentration. This model replicates the qualitative and quantitative behavior of CSD—vasoconstriction, oxygen depletion, potassium elevation, prolonged depolarization—found in experimental studies, and elucidates the effect of oxygen deprivation on CSD recovery. It is shown that during CSD, the metabolic demands of the cortex exceed the physiological limits placed on oxygen delivery, regardless of vascular constriction or dilation. However, vasoconstriction and vasodilation play important roles in determining the susceptibility of the cortical tissue to CSD and its recovery. This study incorporates relevant perfusion and metabolic factors into a model of CSD, and in doing so, helps to explain this phenomenon in vitro and in vivo. IV. CONCLUDING REMARKS

Mathematical modeling of CSD has progressed significantly since its original discovery by Le~ao38 in 1944. Detailed mechanisms have been identified and more and more pieces of the complete puzzle are being added. Although this survey of the mathematical models is not exhaustive, we have discussed some of the historically relevant models that have been proposed. However, in spite of the extensive experiments and modeling attempts to date that have been developed to explain how this complicated phenomenon works, major successes in this direction have remained elusive. Models have been created weaving a plethora of mechanisms into the mix. How these mechanisms qualitatively and quantitatively interact to form a CSD wave remains an open question. In fact, it is not even clear what will constitute an explanation of CSD. Is it sufficient to reproduce the observed experimental phenomena or is it adequate to identify targets for medical treatment of migraine with aura in humans? In the long run, understanding the brain is always our objective, and as Nicholson states,41 an understanding of how the brain works requires an understanding of CSD. ACKNOWLEDGMENTS

The authors thank Horacio Rotstein for the invitation to write this chapter and we thank our colleagues, K. C. Brennan, Josh Chang, Greg Lewis, Yoichiro Mori, Chris Sotak (deceased), Kazuyasu Sugiyama, Shu Takagi, Louis Tao, Philip Wilson, and Wei Yao, who were involved in

Chaos 23, 046103 (2013)

earlier discussions on some of these models. The authors also acknowledge financial support from the National Science Foundation (R.M.M., H.H., J.J.W.), the Natural Sciences and Engineering Research Council of Canada (H.H.), and the Research Grants Council of the Hong Kong Special Administrative Region, China [CityU 104211] (J.J.W.).

1

Even though electro-neutrality is maintained under normal conditions, it may not be the case under CSD. Therefore, it may be necessary to include electro-diffusion. However, in this chapter, we will only consider the diffusive component of extracellular ion transport. 2 T. R. Anderson and R. D. Andrew, “Spreading depression: Imaging and blockade in the rat neo-cortical brain slice,” J. Neurophysiol. 88, 2713–2725 (2002). 3 M. R. Bennett, L. Farnell, and W. G. Gibson, “A quantitative model of purinergic junctional transmission of calcium waves in astrocyte networks,” Biophys. J. 89, 2235–2250 (2005). 4 M. R. Bennett, L. Farnell, and W. G. Gibson, “A quantitative model of cortical spreading depression due to purinergic and gap-junction transmission in astrocyte networks,” Biophys. J. 95, 5648–5660 (2008). 5 K. C. Brennan, M. R. Reyes, H. E. Lopez-Valdes, A. P. Arnold, and A. C. Charles, “Reduced threshold for cortical spreading depression in female mice,” Ann. Neurol. 61, 603–606 (2007). 6 J. Bures, O. Buresova, and J. Krivanek, The Mechanisms and Applications of Le~ ao’s Spreading Depression of Electrical Activity (Academia, Prague, 1974). 7 J. C. Chang, L. L. Shook, J. Biag, E. N. Nguyen, A. W. Toga, A. C. Charles, and K. C. Brennan, “Biphasic direct current shift, haemoglobin desaturation and neurovascular uncoupling in cortical spreading depression,” Brain 133, 996–1012 (2010). 8 J. C. Chang, K. C. Brennan, D. He, H. Huang, R. M. Miura, P. L. Wilson, and J. J. Wylie, “A mathematical model of the metabolic and perfusion effects on cortical spreading depression,” PLoS ONE 8, e70469 (2013). 9 K. C. Chen and C. Nicholson, “Spatial buffering of potassium ions in brain extracellular space,” Biophys. J. 78, 2776–2797 (2000). 10 R. Choudhuri, L. Cui, C. Yong, S. Bowyer, R. M. Klein, K. M. A. Welch, and N. E. J. Berman, “Cortical spreading depression and gene regulation: Relevance to migraine,” Ann. Neurol. 51, 499–506 (2002). 11 M. A. Dahlem, “Migraine generator network and spreading depression dynamics as neuromodulation targets in episodic migraine,” Chaos 23, 046101 (2013). 12 M. A. Dahlem, R. Graf, A. J. Strong, J. P. Dreier, Y. A. Dahlem, M. Sieber, W. Hanke, K. Podoll, and E. Sch€ oll, “Two-dimensional wave patterns of spreading depolarization: Retracting, re-entrant, and stationary waves,” Physica D 239, 889–903 (2010). 13 M. A. Dahlem and T. M. Isele, “Transient localized wave patterns and their application to migraine,” J. Math. Neurosci. 3, 7 (2013). 14 L. X. Dai and R. M. Miura, “A lattice cellular automata model for ion diffusion in the braincell microenvironment and determination of tortuosity and volume fraction,” SIAM J. Appl. Math. 59, 2247–2273 (1999). 15 R. FitzHugh, “Impulses and physiological states in theoretical models of nerve membrane,” Biophys. J. 1, 445–466 (1961). 16 J. Fothergill, “Remarks on that complaint commonly known under the name of sick-headache,” Med. Obs. Inq. 6, 103–137 (1784). 17 A. R. Gardner-Medwin, “A study of the mechanisms by which potassium moves through brain tissues in the rat,” J. Physiol. (London) 335, 353–374 (1983). 18 A. R. Gardner-Medwin and C. Nicholson, “Changes of extracellular potassium activity induced by electric current through brain tissue,” J. Physiol. (London) 335, 375–392 (1983). 19 A. Gorji, “Spreading depression: A review of the clinical relevance,” Brain Res. Rev. 38, 33–60 (2001). 20 A. Gorji, D. Scheller, H. Straub, F. Tegtmeier, R. Koliling, J.-M. Hohling, I. Tuxhorn, A. Ebner, P. Wolf, H. W. Panneck, F. Oppel, and E.-J. Speckmann, “Spreading depression in human neocortical slices,” Brain Res. 906, 74–83 (2001). 21 B. Grafstein, “Mechanism of spreading cortical depression,” J. Neurophysiol. 19, 154–171 (1956). 22 B. Grafstein, “Neuronal release of potassium during spreading depression,” in Brain Function. Cortical Excitability and Steady Potentials, edited by M. A. B. Brazier (University of California Press, 1963), pp. 87–124.


046103-12 23


K. Haerter, C. Ayata, and M. A. Moskowitz, “Cortical spreading depression: A model for understanding migraine biology and future drug targets,” Headache Curr. 2, 97–103 (2005). 24 N. Hadjikhani, M. S. del Rio, O. Wu, D. Schwartz, D. Bakker, B. Fischi, K. K. Kwaong, F. M. Cutrer, B. R. Rosen, R. B. H. Tootell, A. G. Sorensen, and M. A. Moskowitz, “Mechanisms of migraine aura revealed by functional MRI in human visual cortex,” Proc. Natl. Acad. Sci. U.S.A. 98, 4687–4692 (2001). 25 M. Hines and N. T. Carnevale, “The NEURON simulation environment,” Neural Comput. 9, 1179–1209 (1997). 26 A. L. Hodgkin and A. F. Huxley, “A quantitative description of membrane current and is application to conduction and excitation in nerve,” J. Physiol. (London) 117, 500–544 (1952). 27 H. Huang, R. M. Miura, and W. Yao, “A simplified neuronal model for the instigation and propagation of cortical spreading depression,” Adv. Appl. Math. Mech. 3, 759–773 (2011). 28 H. Kager, W. J. Wadman, and G. G. Somjen, “Simulated seizures and spreading depression in a neuron model incorporating interstitial space and ion concentrations,” J. Neurophysiol. 84, 495–512 (2000). 29 H. Kager, W. J. Wadman, and G. G. Somjen, “Conditions for the triggering of spreading depression studied with computer simulations,” J. Neurophysiol. 88, 2700–2712 (2002). 30 H. Kager, W. J. Wadman, and G. G. Somjen, “Computer simulations of neuron-glia interactions mediated by ion flux,” J. Comput. Neurosci. 25, 349–365 (2008). 31 J. Keener and J. Sneyd, Mathematical Physiology (Springer-Verlag, New York, 1998). 32 R. K€ohling, U. R. Koch, G. Hagemann, C. Redecker, H. Straub, and E.-J. Speckmann, “Differential sensitivity to induction of spreading depression by partial disinhibition in chronically epileptic human and rat as compared to native rat neocortical tissue,” Brain Res. 975, 129–134 (2003). 33 R. P. Kraig and P. E. Kunkler, “Spreading depression—A teleologic means for self protection from brain ischemia,” in Cerebrovascular Disease, edited by P. H. Chan (Cambridge University Press, Cambridge, 2002), pp. 142–157. 34 K. Lashley, “Patterns of cerebral integration indicated by scotomas of migraine,” Arch. Neurol. Psychiatr. 46, 331–339 (1941). 35 L. J. Launer, G. M. Terwindt, and M. D. Ferrari, “The prevalence and characteristics of migraine in a population-based cohort: The GEM study,” Neurology 53, 537–542 (1999). 36 M. Lauritzen, “Pathophysiology of the migraine aura. The spreading depression theory,” Brain 117, 199–210 (1994). 37 M. Lauritzen, J. P. Dreier, M. Fabricius, J. A. Hartings, R. Graf, and A. J. Strong, “Clinical relevance of cortical spreading depression in neurological disorders: Migraine, malignant stroke, subarachnoid and intracranial hemorrhage, and traumatic brain injury,” J. Cereb. Blood Flow Metab. 31, 17–35 (2011). 38 A. A. P. Le~ao, “Spreading depression of activity in the cerebral cortex,” J. Neurophysiol. 7, 359–390 (1944). 39 H. Martins-Ferreira, M. Nedergaard, and C. Nicholson, “Perspectives on spreading depression,” Brain Res. Rev. 32, 215–234 (2000). 40 J. Nagumo, S. Arimoto, and S. Yoshizawa, “An active pulse transmission line simulating nerve axon,” Proc. IRE 50, 2061–2070 (1962). 41 C. Nicholson, “Preface II,” in Migraine: Basic Mechanisms and Treatments, edited by A. Lehmenkuhler, K.-H. Grotemeyer, and F. Tegtmeier (Urban

Chaos 23, 046103 (2013) and Schwartzenberg, Munich, 1993), p. 2; and H. Martins-Ferreira, Brain Res. Rev. 32, 6–8 (2000). C. Nicholson, “Diffusion and related transport mechanisms in brain tissue,” Rep. Prog. Phys. 64, 815–884 (2001). 43 J. Olesen, B. Larsen, and M. Lauritzen, “Focal hyperemia followed by spreading oligemia and impaired activation of rCBF in classic migraine,” Ann. Neurol. 9, 344–352 (1981). 44 T. S. Olsen, L. Friberg, and N. A. Lassen, “Ischemia may be the primary cause of the neurologic deficits in classic migraine,” Arch. Neurol. 44, 156–161 (1987). 45 R. K. Orkand, J. G. Nicholls, and A. R. Kuffer, “Effect of nerve impulses on the membrane potential of glial cells in the central nervous system of amphibia,” J. Neurophysiol. 29, 788–806 (1966). 46 L. V. Reshodko and J. Bures, “Computer simulation of reverberating spreading depression in a network of cell automata,” Biol. Cybern. 18, 181–189 (1975). 47 B. E. Shapiro, “An electrophysiological model of gap-junction mediated cortical spreading depression including osmotic volume changes,” Ph.D. thesis (Department of Biomathematics, University of California, Los Angeles, 2000). 48 B. E. Shapiro, “Osmotic forces and gap junctions in spreading depression: A computational model,” J. Comput. Neurosci. 10, 99–120 (2001). 49 G. G. Somjen, “Mechanisms of spreading depression and hypoxic spreading depression-like depolarization,” Physiol. Rev. 81, 1065–1096 (2001). 50 G. G. Somjen, Ions in the Brain, Normal Function, Seizures, and Stroke (Oxford University Press, Oxford, UK, 2004). 51 G. G. Somjen, “Aristides Le~ao’s discovery of cortical spreading depression,” J. Neurophysiol. 94, 2–4 (2005). 52 B. Steinberg, Y. Q. Wang, H. Huang, and R. M. Miura, “Spatial buffering mechanism: Mathematical model and computer simulations,” Math. Biosci. Eng. 2, 675–702 (2005). 53 L. J. Stovner, K. Hagaen, R. Jensen, Z. Katsarava, R. Lipton, A. Scher, T. Steiner, and J. A. Zwart, “The global burden of headache: A documentation of headache prevalence and disability worldwide,” Cephalalgia 27, 193–210 (2007). 54 E. Sugaya, M. Takato, and Y. Noda, “Neuronal and glial activity during spreading depression in cerebral cortex of cat,” J. Neurophysiol. 38, 822–841 (1975). 55 K. Tamura, B. Alessandri, A. Heimann, and O. Kempski, “The effect of a gap-junction blocker, carbenoxolone, on ischemic brain injury and cortical spreading depression,” Neuroscience 194, 262–271 (2011). 56 H. C. Tuckwell and R. M. Miura, “A mathematical model for spreading cortical depression,” Biophys. J. 23, 257–276 (1978). 57 A. van Harreveld, “Compounds in brain extracts causing spreading depression of cerebral cortical activity and contraction of crustacean muscle,” J. Neurochem. 3, 300–315 (1959). 58 A. van Harreveld, “Two mechanisms for spreading depression in the chicken retina,” J. Neurobiol. 9, 419–431 (1978). 59 N. Wiener and A. Rosenblueth, “The mathematical formulation of the problem of conduction of impulses in a network of connected excitable elements, specifically in cardiac muscle,” Arch. Inst. Cardiol. Mex. 16, 205–265 (1946). 60 W. Yao, H. Huang, and R. M. Miura, “A continuum neuronal model for the instigation and propagation of cortical spreading depression,” Bull. Math. Biol. 73, 2773–2790 (2011). 42


Mathematical approaches to modeling development and reprogramming.

Regulation of cerebral metabolism during cortical spreading depression.

Cortical Variability and Challenges for Modeling Approaches.

Neurovascular coupling during cortical spreading depolarization and -depression.

Learned taste aversion induced by cortical spreading depression.

Cortical spreading depression modulates the caudate nucleus activity.

Cerebral hemodynamics during cortical spreading depression in rabbits.

Cortical spreading depression and familial hemiplegic migraine 2015.

Parenchymal spin-lock fMRI signals associated with cortical spreading depression.

Cortical blood flow changes during spreading depression in cats.

Cerebral blood flow and oxygen consumption in cortical spreading depression.

Cortical spreading depression and memory transfer: a methodological critique.

How does spreading depression spread? Physiology and modeling.

Effects of frontal cortical lesions and transverse cortical bisection on spreading depression-induced feeding.

Cortical involvement in tonic immobility ("animal hypnosis"): effect of spreading cortical depression.

Steady potential concomitants of cortical spreading depression: a reply to Huston.

Mathematical modeling of tumor growth and metastatic spreading: validation in tumor-bearing mice.

Habituation to intense acoustic stimulation during cortical spreading depression in rats.

Mathematical approaches to optimization of cancer chemotherapy.

Mathematical modeling.

Geometry Shapes Propagation: Assessing the Presence and Absence of Cortical Symmetries through a Computational Model of Cortical Spreading Depression.

Cortical Spreading Depression Causes Unique Dysregulation of Inflammatory Pathways in a Transgenic Mouse Model of Migraine.

Propagation of cortical spreading depression into the hippocampus: The role of the entorhinal cortex.

Modulation of intracranial meningeal nociceptor activity by cortical spreading depression: a reassessment.