User:Markus A. Dahlem/Proposed/Models of cortical spreading depression

• Markus A. Dahlem, Department of Physics, Humboldt-University of Berlin, Berlin, Germany

• Mr. Niklas Hübel, Technische Universität Berlin, Berlin, Germany

Dr. Markus A. Dahlem accepted the invitation on 25 February 2009 (self-imposed deadline: 31. January 2013).

Models of cortical spreading depression refers to mathematical descriptions of cortical spreading depression. Cortical spreading depression is a massive but temporary perturbation in the cortical ionic homeostasis leading to a depression of neuronal activity that spreads (hence the name) through the cortex and other gray matter regions in the brain. Cortical spreading depression is an emergent phenomenon arising from the interaction local nonlinear processes with different volume transmission phenomena in cortical tissue. Neural and glial compartments, the extracellular space, and the vasculature play a role in cortical spreading depression (SD).

Introduction: The pathophysiological basis for modeling SD

This is a breif review focusing mainyl on simple cellular and generic macroscopic models of SD. Note that there are also experimental models of this phenomenon referring to a particular experimental procedure, animal species, and/or brain region in which SD is studied.

Terminology

In this article, SD is used as the abbreviation for spreading depression. Omitting the "C" is advisable as some experimental data are not obtained from the cortex. Moreover, this phenomenon is in recent years frequently called spreading depolarization (Dreier, 2011), a term with the same abbreviation but that reflects what is thought to reflect better the mechanism, namely the nearly complete depolarization of the membrane potential.

Distinctive features of SD

During SD a massive redistribution of ions across cell membranes and nearly eliminates the large gradients across the cell membrane. This redistribution peaks after several seconds in a nearly complete neuronal depolarization with a depression of neural activity, followed by a much slower recovery process taking up to minutes during which ion gradients are re-established towards their physiological values. This maximal ionic perturbation clearly distinguishes SD from all other brain states such as epileptic seizure activity, functional activation or the physiological resting state. The sequence of ionic perturbation and its recovery spreads at a pace of about 3~mm/min over cortical regions.

Clinical relevance of SD

Cortical spreading depression is closely related to migraine with aura and several of of stroke and brain injuries (Dreier, 2011). A pathway was recently found explaining how SD triggers headache in migraine by activating large neuronal transmembrane channels connecting the intracellular and extracellular space (Karatas, 2013). In stroke, re-entrant SD waves are believed to have the potential to worsen outcome in incremental steps with each wave circling near the infarct core (anatomical block) but could also have some bene�cial component by stimulating blood ow in the penumbra zone far from the infarct tissue (Strong, 2007; Dahlem et al, 2010).

Early modeling attempts

Figure 1: First page from Wiener and Rosenblueth paper (1941). Copy with annotations made by the late Art Winfree.

Figure 2: Spiral-shaped retinal SD, overlay.

In 1946, Norbert Wiener and Arturo Rosenblueth proposed a mathematical framework to describe wave patterns explaining certain types of cardiac arrhythmia. The authors made also the connection to brain dynamics, stating that ``Nervous elements and cardiac [...] fibers are excitable. [...] The laws which apply to the muscle fibers are also applicable to the nerve fibers``. Wiener and Rosenblueth knew very well, in fact, that waves similar to the those in cardiac arrhythmia occur in the cortex. Arturo Rosenblueth, in collaboration with Hallowell Davis, supervised a few years earlier a young Brazilian PhD student, Aristides A. P. Leao who was working in the department of Physiology at Harvard Medical School and discovered the phenomenon of SD in 1944.

Leão stated in (1945) "Much has been written about vascular phenomena both in clinical epilepsy and the presumably related condition of migraine. The latter disease with the marked dilatation of major blood vessels and the slow march of scotomata in the visual or somatic sensory sphere is suggestively similar to the experimental phenomenon here described, in spite of the fact that known scotomata are still felt to be vasoconstrictor in nature."

In Figure 1, the title page of this seminal paper by Wiener and Rosenblueth is shown with many annotations made 1970 by Art Winfree, in particular, mentioning rotors, that is, a special concept of re-entrant waves that are not pinned by an anatomical block but rotate freely. These waves are also called spiral waves. Winfree was a leading theoretical biologist who made major contributions to the field of biological oscillations. In 1972, he published a paper on Spiral Waves of Chemical Activity observed in the Belousov-Zhabotinsky (BZ) reaction. Three years later, Reshodko and Bures (1975) published the first computer simulation of re-entrant SD waves in a network of cell automata.

These early modeling attempts that focus on the two-dimensional propagation pattern of SD on a centimeter scale have recently regained much interest, for monitoring the human brain provided evidence of re-entrant SD patterns during stroke. Stroke outcome could critically depend on SD as these events pose dramatic metabolic stress to the tissue. Whether or not the loss of potentially salvageable tissue, i.e., tissue at risk of infarction, is increased by the number of SD is a current research focus. Computational model of SD complement traditional research methodologies (Ruppin et al. 1999).

Neurons and also glial cells in brain tissue are in a non-equilibrium state stationary state. During SD, the brain tissue approaches its equilibrium state, which would eventually lead to the loss of this tissue, if the recovery to the non-equilibrium state stationary is not achieved within minutes. In a nutshell for physicists, SD is thus a pulse-like solitary wave in which the tissue temporarily approaches nearly an equilibrium state and is then driven back far from it.

Ion-based model

In a nutshell, physiologically realistic models of SD are based on a Hodgkin-Huxley type description, that is, they use Kirchhoff's current law the describe electrophysiological events across the neuronal cell membrane in terms of voltage-gated currents. Such models are also refered to conductance-based models, because the other dynamical varibabes---beside the transmembrane Potential---are gating-variables that determine the conductance of the ions.

As an extention to the original formalism the Hodgkin-Huxley model of action potentials, models of SD must also consider surrounding compartments to account for changes in the intra- and extracelluar ion concentrations. We call these models ion-based models. In the literature of myocard action potential, ion-based models are also called second-generation Hodgkin-Huxley model (Arce et al. 2000, Endersen and Skarland 2000). It is often mentioned in this literature that such second-generation models can show generacy of equilibrium and show very slow long-term drifts of ion concentrations. A detailed comparison between ion-based models in neurons and cardiac myocytes is beyond the scope of this review, however.

Under physiological conditions, the chemical gradients due to the different ion concentrations inside and outside the cells are rather stable batteries and can therefore be considered as being constant in most neural models. By contrast, the characteristic event in SD is the break-down of the ionic homeostasis. To keep track of the changes in the ion concentration three phenomena are usually accounted for in models of SD:

• the Goldman-Hodgkin-Katz (GHK) current equation,
• ion exchange pumps,
• glial buffering,
• diffusion from and into the blood and
• cell swelling.

Obtaining the spatial continuum limit of such discrete microscopic cellular models of SD is by no means straightforward, although efforts can be made in this direction. The continuum limit is necessary to model the spread, which can under conditions of migraine and stroke extend over several centimeters in the human cortex. Continuum limit requires appropriate conditions on the spatial coupling, that is, on the communication pathways between cortical and glial elements and the neurovascular coupling. Communication between neurons is usually classified in the cortex by two schemes:

• wiring transmission and
• volume transmission.

Volume transmission is mainly characterized by diffusion of chemical signals in the extracellular space and ephaptic coupling through extracellular electric fields. Wiring transmission is based on the structural and functional cortical connectivity, and other intercellular communication occurring through a well-defined network. To approach the continuum limit in a mathematical model, an effective medium theory is needed. This approach can become increasingly complicated, in particular because the neural and glial elements swell considerably under conditions of SD. Already the 3D flow of a liquid through a porous medium with only absorbing and swelling elements but without any nonlinear reactions is a hard problem (Fasano and Mikelic, 2002).

The detailed ion-based model contains a neural (somatic and/or dentritic) and a glial compartment plus an extracellular space. There is a neural and a glia membrane seperating the respective compartments from the extracellular space. The ion species considered in the model are sodium, potassium and chloride with concetrations \([\mathit{Na}]_{cpt}\), \([K]_{cpt}\) and \([\mathit{Cl}]_{cpt}\), where the subscript "cpt" specifies the compartment ("so" for soma, "gl" for glia and "ex" for extracellular space). The soma membrane contains active ("gated") sodium and potassium channels and additional passive ("leak") channels for all ion species. We do not include in this first descrition a dentritic compartment. This modeling scheme can easily be extended, however. Also, for example, ion fluxes \(Ca^{2+}\) can likewise be computed and their feedback effect on the neuron can be taken into account. Furthermore, for both \(Na^+\) and \(K^+\) various currents can be implemented with possibly different distribution in the various neural comaprtments (see Fig. ). For \(Na^+\) a fast transient sodium current but not the persistent sodium current. For \(K^+\) we consider only a simple rectifying potassium current.

The glia membrane has passive channels only.

The dynamical variables of the model are both membrane potentials, all ion concentrations, the gating variables and the volumes of the compartments.

Some of these 16 variables are not independent though. For example we will assume the volume of the entire system to be constant. Also some ion concentrations follow from conservation of the total number of ions of the species.

Model dynamics are given by standard rate equations for the membrane potentials

\( \dot{V}_{so} =\frac{1}{C_{so}}\cdot(-I_{Na}^{(so)}-I_K^{(so)}-I_{Cl}^{(so)}+I_{app}) \tag{1} \)

\( \dot{V}_{gl} =\frac{1}{C_{gl}}\cdot(-I_{Na}^{(gl)}-I_K^{(gl)}-I_{Cl}^{(gl)}) \tag{2} \)

Soma currents for sodium and potassium depend on gating variables, which obey Hogkin-Huxley-like dynamics.

\( \dot{n} =\phi\cdot[\alpha_n(V_{so})(1-n)-\beta_n(V_{so})n] \tag{3} \)

\( \dot{h} =\phi\cdot[\alpha_h(V_{so})(1-h)-\beta_h(V_{so})h] \tag{4} \)

Changes in ion concentrations are due to transmembrane currents and in case of extracellular potassium diffusion to an extracellular bath.

\( \dot{[\mathit{Na}]}_{so} =-[\mathit{Na}]_{so}\cdot\frac{\dot{\mathit{Vol}_{so}}}{\mathit{Vol}_{so}}+\gamma \cdot\frac{\mathit{Vol}_{so}^0}{\mathit{Vol}_{so}}\cdot-I_{Na}^{(so)} \tag{5} \)

\( \dot{[\mathit{Na}]}_{ex} =-[\mathit{Na}]_{ex}\cdot\frac{\dot{\mathit{Vol}_{ex}}}{\mathit{Vol}_{ex}}+\gamma \cdot\frac{\mathit{Vol}_{so}^0}{\mathit{Vol}_{ex}}\cdot(I_{Na}^{(so)}+I_{Na}^{(gl)}) \tag{6} \)

\( \dot{[\mathit{Na}]}_{gl}=-[\mathit{Na}]_{gl}\cdot\frac{\dot{\mathit{Vol}_{gl}}}{\mathit{Vol}_{gl}}+\gamma \cdot\frac{\mathit{Vol}_{so}^0}{\mathit{Vol}_{gl}}\cdot-I_{Na}^{(gl)} \tag{7} \)

\( \dot{[K]}_{so} =-[K]_{so}\cdot\frac{\dot{\mathit{Vol}_{so}}}{\mathit{Vol}_{so}} +\gamma \cdot\frac{\mathit{Vol}_{so}^0}{\mathit{Vol}_{so}}\cdot-I_K^{(so)} \tag{8} \)

\( \dot{[K]}_{ex} =-[K]_{ex}\cdot\frac{\dot{\mathit{Vol}_{ex}}}{\mathit{Vol}_{ex}} +\gamma \cdot\frac{\mathit{Vol}_{so}^0}{\mathit{Vol}_{ex}}\cdot(I_K^{(so)}+I_K^{(gl)}) + J_{diff} \tag{9} \)

\( \dot{[K]}_{gl} =-[K]_{gl}\cdot\frac{\dot{\mathit{Vol}_{gl}}}{\mathit{Vol}_{gl}} +\gamma \cdot\frac{\mathit{Vol}_{so}^0}{\mathit{Vol}_{gl}}\cdot-I_K^{(gl)} \tag{10} \)

\( \dot{[\mathit{Cl}]}_{so}=-[\mathit{Cl}]_{so}\cdot\frac{\dot{\mathit{Vol}_{so}}}{\mathit{Vol}_{so}}+\gamma \cdot\frac{\mathit{Vol}_{so}^0}{\mathit{Vol}_{so}}\cdot-I_{Cl}^{(so)} \tag{11} \)

\( \dot{[\mathit{Cl}]}_{ex}=-[\mathit{Cl}]_{ex}\cdot\frac{\dot{\mathit{Vol}_{ex}}}{\mathit{Vol}_{ex}}+\gamma \cdot\frac{\mathit{Vol}_{so}^0}{\mathit{Vol}_{ex}}\cdot(I_{Cl}^{(so)}+I_{Cl}^{(gl)}) \tag{12} \)

\( \dot{[\mathit{Cl}]}_{gl}=-[\mathit{Cl}]_{gl}\cdot\frac{\dot{\mathit{Vol}_{gl}}}{\mathit{Vol}_{gl}}+\gamma \cdot\frac{\mathit{Vol}_{so}^0}{\mathit{Vol}_{gl}}\cdot-I_{Cl}^{(gl)} \tag{13} \)

Compartment volumes change according to osmotic imbalances. Dynamics are given by a first order process with asymptotic values depending on ion concentrations in neighbouring compartments (see below).

\( \dot{\mathit{Vol}_{so}} =\frac{\mathit{Vol}_{so}^\infty-\mathit{Vol}_{so}}{250} \tag{14} \)

\(\dot{\mathit{Vol}_{ex}} =\frac{\mathit{Vol}_{ex}^\infty-\mathit{Vol}_{ex}}{250} \tag{15}\)

\( \dot{\mathit{Vol}_{gl}} =\frac{\mathit{Vol}_{gl}^\infty-\mathit{Vol}_{gl}}{250}\tag{16} \)

Quantities occuring these equations are listed in the table below. The factor \(\gamma\) converts currents to ion fluxes in the appropriate equations at initial volume size \( \gamma = \frac{\mathit{Surf}_{so}}{\mathit{Vol}_{so}^0F} \), with initial soma volume \(\mathit{Vol}_{so}^0\), Faraday's constant \(F\) and mebrane surface \(\mathit{Surf}_{so}\). Mutliplication with appropiate volume fractions give the right conversion for each compartment.

Currents

Transmembrane currents are computed as the product of some (potentially voltage-dependent) conductance and a Nernst-like voltage term. Conductances are

\( g_{Na}^{(so)} = g_{Na,l}^{(so)}+g_{Na,g}^{(so)}m_\infty^3h \tag{17} , \)

\( g_{K}^{(so)} = g_{K,l}^{(so)}+g_{K,g}^{(so)}n^4 \tag{18} , \)

\( g_{Cl}^{(so)} = g_{Cl,l}^{(so)} \tag{19} , \)

where \(m_\infty\) is the equilibrium value of the fast sodium activator. Subscripts \(l\) and \(g\) denote leak and gated channel conductances, respectively. Glia conductances \(g_{Na}^{(gl)}\), \(g_K^{(gl)}\) and \(g_{Cl}^{(gl)}\) are only leak conductances. Reversal potentials \(E_{ion}^{(cpt)}\) have to be computed for every timestep, because ion concentrations are changing.

\( E_{Na}^{(so)} =+26{.}64\ \ln([\mathit{Na}]_{ex}/[\mathit{Na}]_{so}), \tag{20}\\ \)

\( E_{K}^{(so)} =+26{.}64\ \ln([K]_{ex}/[K]_{so}) \tag{21} \\ \)

\( E_{Cl}^{(so)} =-26{.}64\ \ln([\mathit{Cl}]_{ex}/[\mathit{Cl}]_{so})\tag{22}\\ \)

\( E_{Na}^{(gl)} =+26{.}64\ \ln([\mathit{Na}]_{ex}/[\mathit{Na}]_{gl})\tag{23}\\ \)

\( E_{K}^{(gl)} =+26{.}64\ \ln([K]_{ex}/[K]_{gl}) \tag{24} \\ \)

\( E_{Cl}^{(gl)} =-26{.}64\ \ln([\mathit{Cl}]_{ex}/[\mathit{Cl}]_{gl})\tag{25} \)

We also include ATPase as a pump current that is added to the sodium and potassium current in each channel. It is given by

\( I_P^{(cpt)}=\frac{I_{max}^{(cpt)}}{\left(1+\exp[(\mathit{Na}_{P}^{(cpt)}-[\mathit{Na}]_{cpt})/3]\right)\cdot\left(1+\exp(K_P^{(cpt)}-[K]_{ex})\right)}, \tag{26} \)

with pump constants \(\mathit{Na}_{P}^{(so)}\), \(\mathit{Na}_{P}^{(gl)}\), \(K_P^{(so)}\), \(K_P^{(gl)}\) refering to the considered mebrane and maximal pump rates \(I_{max}^{(so)}\) and \(I_{max}^{(gl)}\). We introduce an overall factor \(f_{glia}\) to adjust the glial strength. Explicit formulae for the currents occuring in the membrane rate equations ((1), (2)) are

\( I_{Na}^{(so)} =g_{Na}^{(so)}\cdot(V_{so}-E_{Na}^{(so)})+3I_P^{(so)} \tag{27} \\ \)

\( I_K^{(so)} =g_{K}^{(so)} \cdot(V_{so}-E_{K}^{(so)} )-2I_P^{(so)} \tag{28} \\ \)

\( I_{Cl}^{(so)}=g_{Cl}^{(so)}\cdot(V_{so}-E_{Cl}^{(so)}) \tag{29} \\ \)

\( I_{Na}^{(gl)}=f_{glia}\cdot[g_{Na}^{(gl)}\cdot(V_{gl}-E_{Na}^{(gl)})+3I_P^{(gl)}]\tag{30} \\ \)

\( I_K^{(gl)} =f_{glia}\cdot[g_{K}^{(gl)} \cdot(V_{gl}-E_{K}^{(gl)} )-2I_P^{(gl)}]\tag{31} \\ \)

\( I_{Cl}^{(gl)}=f_{glia}\cdot[g_{Cl}^{(gl)}\cdot(V_{gl}-E_{Cl}^{(gl)})] \tag{32} \)

Gating Variables

The gating variables \(m\) (see adiabatic value \(m_\infty\) in eq.\ ((17))), \(h\) and \(n\) are refered to as sodium activator, sodium inactivator and potassium activator, respectively. They obey Hodgkin-Huxely-like dynamics defined by the following exponential functions\[ \alpha_m(V_{so})=0{.}1\cdot\frac{V_{so}+30}{1-\exp[-0{.}1(V_{so}+30)]} \ ,\tag{33} \]

\( \beta_m(V_{so}) =4\cdot\exp[-(V_{so}+55)/18] \ ,\tag{34} \)

\( \alpha_n(V_{so})=0{.}01\cdot\frac{V_{so}+34}{1-\exp[-0{.}1(V_{so}+34)]}\ ,\tag{35} \)

\( \beta_n(V_{so}) =0{.}125\cdot\exp[-(V_{so}+44)/80] \ ,\tag{36} \)

\( \alpha_h(V_{so})=0{.}07\cdot\exp[-(V_{so}+44) /20] \ ,\tag{37} \)

\( \beta_h(V_{so}) =\frac{1}{1+\exp[-0{.}1(V_{so}+14)]} \ .\tag{38} \)

Analogously to eq.\ ((3), (4)) it is understood that

\( \dot{m}=\phi\cdot[\alpha_m(V_{so})(1-m)-\beta_m(V_{so})m]\ .\tag{39} \)

Equilibrium values follow immediately from the gating kinetics\[ m_\infty(V_{so})=\frac{\alpha_m(V_{so})}{\alpha_m(V_{so})+\beta_m(V_{so})}\ ,\tag{40} \\ \]

\( h_\infty(V_{so})=\frac{\alpha_h(V_{so})}{\alpha_h(V_{so})+\beta_h(V_{so})}\ ,\tag{41}\\ \)

\( n_\infty(V_{so})=\frac{\alpha_n(V_{so})}{\alpha_n(V_{so})+\beta_n(V_{so})}\ .\tag{42} \)

Changes in Ion Concentrations

The rate equations for the different ion concentrations in each compartment contain a contribution from ion flux due to transmembrane currents and a term accounting for concentration changes due to volume swelling or shrinking (cp.\ second and first term in eq.\ ((5)--(13)), respectively). We neglect changes of membrane surface size, but only consider how the transformation from current to ion flux depends on the compartment volumes. In addition extracellular potassium is regulated by a diffusion-like ion flux \(J_{diff}\) mimicing an extracellular potassium bath. It is defined by the diffusion strength \(f_{diff}\) and the potassium concentration \(K_{bath}\) of the bath\[ J_{diff}=f_{diff}\cdot([K]_{ex}-K_{bath})\tag{43} \]

Volume Changes

Coming soon \ldots

Stimulation

We investigate two different types of stimulation in this model. One option is to apply some unspecified current \(I_{app}\) to the soma (cp.\ eq.\ ((1))). This current is usally a steppulse specified by its amplitude and duration, which is applied at a certain time. We observe and analyse the type of phase space excursion such stimulus causes. Sometimes we also apply a constant current. The other option is to increase the potassium bath concetration \(K_{bath}\) and drive the system from a stable regime to different types of oscillatory motion.

Parameters

The model contains a large number of parameters, most of which are fixed. Default values are listed below.

Figure 3: Phase space excursion during SD

Macroscopic and other approaches

An alternative to the physiologically realistic ion-based models---with an yet unknown to otain the continuum limit on equal on equal footage of realistic detail---is a macroscopic model utilizing a top-down approach to describe patterns within the cortical tissue.

This seems particularly attractive for spatio-temporal patterns of SD on a space scale of several centimeters lasting minutes to hours. These space and time scales are in striking contrast to the millisecond and micrometer range of the physiologically realistic models and these scales are highly relevant to various clinical questions.

In particular, two-dimensional SD wave patterns that retract, reenter or are stationary waves have been attributed to migraine, stroke, and forms in-between, such as persistent migraine without infarction, migrainous infarction, and ischemia-induced migraine, respectively (Dahlem et al. 2010).

In the case of an macroscopic model, the mechanism of SD is described as a generic reaction-diffusion system. The macroscopic patterns obtained form these reaction-diffusion models might also help to understand what is called the migraine-stroke continuum. At least an approach from the other end, complementary to the microscopic models, is needed to cover the wide spectrum of phenomenons related to SD.

Furthermore, there are cellular automaton models of SD and kinematic descriptions of the SD wave front, but these models play a minor role to date, though there is some historical interest. Finally, there are hybrid models, that is, models using a macroscopic reaction-diffusion description in combination with cellular neural network models. Such a model model was first suggested to explain visual hallucinations caused by SD during the course of a migraine attack (Reggia and Montgomery, 1996), but recent developments focus the attention now on the neural circuity as an augmented transmission scheme that provides weak nonlocal coupling and thereby it may modulate cortical susceptibility to SD.

The articles is organized as follows. We start with ion-based models of SD. These models ground on a solid biophysical understanding. They are considered as physiologically realistic although, of course, only a limited number of dynamical variables and compartments can be considered.

In the second part, we take the top-down macroscopic approach. In this part, we will chronologically present the material, starting with a brief review of some early attempts to model SD. We end with a discussion.

Hodgkin-Grafstein model of SD

[Pending]

Canonical reaction-diffusion model of activator inhibitor type

The diversity of the behavior of traveling waves in two spatial dimensions was studied in canonical models depending on two generic parameters \(\beta\) and \(\varepsilon\) in Eqs.\,(\ref{eq:fhn1})-(\ref{eq:fhn2}), which determine the parameter plane of excitability (Winfree 1991) (in the system without mean field feedback control). The excitability of media is not a scalar quantity but rather a classification of bifurcation-separated regimes in this plane in which the same kind of traveling wave solutions exist. However, in a loose way the concept of excitability can be treated as if it were a scalar quantity in the following sense. Usually, patterns of discontinuous waves (with open ends) are used to probe excitability. For instance, the propagation of discontinuous waves usually evolve to spiral waves. They can rotate rigidly (without changing their shape), or---by changing either parameter in a direction that might be considered as increasing the excitability---the rotation frequency can destabilize in a Hopf bifurcation such that meandering or drifting spiral waves bifurcate (Barkely 1991). By increasing excitability further, the motion of the spiral center becomes more complex (Winfree 1991), a behavior that as been observed in spiral-shaped SD waves in an {\em in vitro} animal model (Dahlem & Müller 1997). The spiral wave can also break up, either in the core or in the far field by some instabilities (San 2007). The signature of break-up in the core field was also observed in this {\em in vitro} animal model (Dahlem & Müller 1997), so that in principal, the generic parameters \(\beta\) and \(\varepsilon\) can be identified by comparing the predicted patterns with experimental data.

In the other direction, that is, lowering excitability, the non-excited circular core of rigidly rotating spirals enlarges. At a critical excitability, the spiral becomes a half plan wave and its core can be considered as being infinitely large. In other words, at this critical excitability, called the rotor boundary \(\partial R_{\infty}\), discontinuous waves do not curl-in and re-enter the same medium. Beyond the rotor boundary, discontinuous waves start to retract at their open ends and any discontinuous waves will disappear (Mikhailov & Zykov 1991, Hakim & Karma 1999). The border <math>\partial R_{\infty}<math> is a saddle-node bifurcation at which discontinuous waves collide with their corresponding nucleation solution. This leads to the key idea of our model. The linear mean field feedback control, as we will describe in the following, moves this saddle-node bifurcation towards distinct localized wave segments with a characteristic form (shape, size) and behind this bifurcation these waves become transient objects.

Hybrid models

A continuous reaction-diffusion SD model forcing a discrete neural network is what we consider a hybrid model. Although there is only a single paper on this, we feel this deserves mentioning, because in the future one should think about replacing such discrete neural network with continuous neural field models. It is only for historical reasons that SD was studied in conductance-based HH models and their extensions to ion-based models, while rate-based or activity-based neural field models may as well be extended toward incooperating the depletion ion gradients.

Figure 4: A model of hallucinatory zigzag caused by SD wave .

Figure 5: Simulation.

In the mid 1990s, Reggia and Montgomery (1996) proposed a SD model to explain the hallucinatory zigzag pattern that is thought to be caused by the initial excitation phase of SD. To test this, at that time controversial, hypothesis, they extended a reaction-diffusion (RD) model based on potassium dynamics, which was similar to the Hodgkin-Grafstein model only with a quartic polynomial, by including a second dynamic variable describing a recovery phase, identical to the FitzHugh-Nagumo inhibitor rate equation, and connected this RD model to a neural network. The neural network was initially not developed to study migraine aura but used to study cortical dynamics and sensory map reorganization.

The distinguishing and novel feature of the Reggia and Montgomery (1996) SD model is that it connected two types of models to become a hybrid model. The key result of their simulations is that on the leading edge of the potassium wave, the elsewhere largely uniform neural activity was replaced by a pattern of small, irregular patches and lines of highly active elements. In 1971, Richards suggested in a Scientific American article about visual migraine aura that the zigzag seen in a visual aura is caused by the spatial layout of a specific type of neurons in visual cortex. Therefore, SD by evoking migraine hallucinations provides a route to the neural mechanisms of normal brain function.

It was argued that the approaching SD wave initially affects cortical cells located in a special layer of the cortex (IVc) that possess the highest spontaneous activity. In layer IVc, the afferent inputs to the cortex terminate, and cells located there posses a concentric circular receptive field structure. Cells in the surrounding laminae possess fairly low spontaneous discharge rate, and presumably are strongly inhibited by intracortical inhibition. Therefore, those cells which converge to orientations which are perpendicular to the SD front will be inhibited first. Their excitatory input to the surrounding simple cells will be attenuated, and the activity of these surrounding cells will be depressed. Correspondingly, their neighbors on either side will be disinhibited, and will experience an elevation in firing rate. These columns correspond to orientations which are in the range of 30° to 60° from perpendicular to the advancing front, and thus would be perceived as the jagged fortification.

Neural network models have become sophisticated enough to constrain and validate possible an SD attack underlying cortical circuitry. Simulations with SD driven population codes interpreting cortical feature maps an transforming them into a percept in the visual field have lead to simulating of the aura that allow the patient to virtually re-experience their disorders (see Figure 4). Preliminary results indicate that such networks provide realistic simulations that can be matched by patients against their visual symptoms.

Figure 6: Flow chart of feedback loops and pathways in SD.

References

• Arce H., Xu A., González H., Guevara M.R., Alternans and higher-order rhythms in an ionic model of a sheet of ischemic ventricular muscle. Chaos 10,411 (2000).

• Barkley D., Euclidean symmetry and the dynamics of rotating spiral waves. Phys. Rev. Lett. 72,164(1994).

• Dreier, J. P., The role of spreading depression, spreading depolarization and spreading ischemia in neurological disease, Nat. Med. 17, 439, 2011.

• Dahlem M.A., Müller S.C., Self-induced splitting of spiral-shaped spreading depression waves in chicken retina. Exp. Brain Res. 115,319 (1997).

• Dahlem, M. A. , Graf, R. , Strong, A. J. , Dreier, J. P. , Dahlem, Y. A., Sieber, M., Hanke, W. , Podoll, K. and Schöll, E., Two-dimensional wave patterns of spreading depolarization: retracting, re-entrant, and stationary waves, Physica D 239, 889, 2010.

• Dahlem M.A., Schneider F.M., Schöll E., Failure of feedback as a putative common mechanism of spreading depolarizations in migraine and stroke. Chaos 18,026110 (2008).

• Endresen L.P. and Skarland N., Limit cycle oscillations in pacemaker cells. IEEE Trans Biomed Eng 47,1134 (2000)

• Fasano, A. and Mikelic, A., The 3D flow of a liquid through a porous medium with absorbing and swelling granules, Interfaces and Free Boundaries 4, 239, 2002.

• Hakim V., Karma A. Theory of spiral wave dynamics in weakly excitable media: asymptotic reduction to a kinematic model and applications. Phys Rev E 60,5073 (1999).

• Kager, H., 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).

• Karatas, H. , Erdener, S. E. , Gursoy-Ozdemir, Y. , Lule, S. , Eren-Kocak, E. , Sen, Z. D. and Dalkara, T. , Spreading depression triggers headache by activating neuronal Panx1 channels, Science 339, 6123, 1092 (2013).

• Leão A.A.P., Spreading depression of activity in the cerebral cortex. J. Neurophysiol 7,359–390, 1944.

• Leão A.A.P., Morrison R.S., Propagation of spreading cortical depression. J. Neurophysiol 8:33-46, 1945.

• Mikhailov A.S., Zykov V.S. Kinematical theory of spiral waves in excitable media: comparison with numerical simulations. Physica D 52,379 (1991).

• Reggia, J. A. and Montgomery, D. , A computational model of visual hallucinations in migraine, Comput. Biol. Med. 26, 133-139 (1996).

• Ruppin, E. , Revett, K. , Ofer, E. , Goodall, S. and Reggia, J. A. , Penumbral tissue damage following acute stroke: a computational investigation, Prog. Brain Res. 121, 243-255 (1999)

• Reshodko, L. V. and Bures J., Computer simulation of reverberating spreading depression in a network of cell automata, Biol. Cybern. 18, 181-189 (1975).

• Sandstede B., Scheel A., Period-doubling of spiral waves and defects. SIAM Journal on Applied Dynamical Systems 6,494 (2007).

• A. J. Strong, P. J. Anderson, H. R. Watts, D. J. Virley, A. Lloyd, E. A. Irving, T. Nagafuji,

M. Ninomiya, H. Nakamura, A. K. Dunn, and R. Graf, Peri-infarct depolarizations lead to loss of perfusion in ischaemic gyrencephalic cerebral cortex, Brain 130, 995-1008 (2007).

• Winfree, A. T., Spiral Waves of Chemical Activity, Science 175, 634-36 (1972).

Recommended reading

• Somjen,G.G., Ions in the brain: normal function, seizures, and stroke, Oxford University Press (2004)

Internal references

• Gregoire Nicolis and Anne De Wit (2007) Reaction-diffusion systems. Scholarpedia, 2(9):1475.
• Eugene M. Izhikevich and Richard FitzHugh (2006) FitzHugh-Nagumo model. Scholarpedia, 1(9):1349.