Linking a genetic defect in migraine to spreading depression in a computational model

Hodgkin–Huxley model

The Hodgkin–Huxley (HH) model is one of the most widely used computational models in neuroscience. It is a conductance-based neuron model (Hodgkin & Huxley, 1952) and consists of four differential equations describing the membrane potential V and three gating variables m, n and h that determine the conductances of potassium and sodium channels. The change in membrane potential is proportional to the current that is flowing across the membrane with the proportionality constant given by the capacitance of the membrane C_m. The individual currents are modeled as the conductance g_i of the respective channel times the driving force, which is given by the difference between the membrane potential and the respective ion’s reversal potential E_i, where i∈{K, Na, leak}. Note that the conductance g_j for voltage-gated channels, i.e., j∈{K, Na}, is given by the maximal conductance ${\bar{g}}_j$ times the respective gating variables as introduced below. The model takes into account a sodium current I_Na⁺, a potassium current I_K⁺, a leak current I_leak that is carried by unspecified ions, and an applied current I_app. (1)${\displaystyle \frac{\mathrm{d}V}{\mathrm{d}t}=-\frac{1}{C_m}\left(I_{{\mathrm{Na}}^{+}}+I_{{\mathrm{K}}^{+}}+I_{\mathrm{leak}}-I_{\mathrm{app}}\right),}$ (2)${\displaystyle I_{{\mathrm{Na}}^{+}}={\bar{g}}_{\mathrm{Na}}{m}^{3}h\left(V-E_{\mathrm{Na}}\right),}$ (3)${\displaystyle I_{{\mathrm{K}}^{+}}={\bar{g}}_{\mathrm{K}}{n}^4\left(V-E_{\mathrm{K}}\right),}$ (4)${\displaystyle I_{\mathrm{leak}}=g_l\left(V-E_{\mathrm{leak}}\right).}$ In the HH model the potassium current is modeled as a delayed rectifier current with activation gate n while the sodium current is described by a transient current with an activation gate m and an inactivation gate h. All gating variables are voltage dependent and are given by the following equations: (5)${\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}=\frac{x_{\infty }-x}{{\tau }_x}\mathrm{with}}$ (6)${\displaystyle x_{\infty }=\frac{{\alpha }_x}{{\alpha }_x+{\beta }_x}\mathrm{and}}$ (7)${\displaystyle {\tau }_x=\frac{1}{{\alpha }_x+{\beta }_x}\mathrm{for}x\in \left\{n,m,h\right\}.}$ x_∞ describes the steady-state of the gating variables and τ_x is the time constant.

The rate equations for α_x and β_x are voltage-dependent and given by (8)${\displaystyle {\alpha }_m=\frac{0.1\left(V+40\right)}{1-exp\left(-\left(V+40\right)/10\right)},}$ (9)${\displaystyle {\beta }_m=4exp\left(-\left(V+65\right)/18\right),}$ (10)${\displaystyle {\alpha }_n=\frac{0.01\left(V+55\right)}{1-exp\left(-\left(V+55\right)/10\right)},}$ (11)${\displaystyle {\beta }_n=0.125exp\left(-\left(V+65\right)/80\right),}$ (12)${\displaystyle {\alpha }_h=0.07exp\left(-\left(V+65\right)/20\right),}$ (13)${\displaystyle {\beta }_h=\frac{1}{1+exp\left(-0.1\left(V+35\right)\right)}.}$ This model is capable of producing action potentials in response to depolarizations of the membrane caused by an appropriate externally applied current I_app. All model parameters that were used in the simulations of the HH model can be found in Table 1. It is interesting to remark that trying to study the effect of the mutation in a reduced two-dimensional model in the phase plane did not lead to promising results because the mutation quickly led to bistability, which is consistant with our results of a prolonged plateau of action potential and early depolarization block in the form of bistability.

Spreading depression model

The classical HH model neglects the time-dependency of ion concentrations caused by spiking dynamics. Ions accumulate very slowly but also progressively due to the fluxes across the neuronal membrane. Therefore, changes in concentrations become significant either in the course of many rapid action potentials or under metabolic stress with insufficient ion pump activity, such as during transient ischemic attacks. Hence both the onset of spiking and also the response to reduced ion pump activity are of interest. These can be modeled by the spreading depression model described in more detail by Hübel, Schöll & Dahlem (2014).

This model is also based on HH dynamics, but uses several changes and extensions. Instead of an unspecified leak current, a combined Na⁺–K⁺-leak current is used. The equations for sodium and potassium currents, including a pump current I_p that is introduced below, therefore change to (14)${\displaystyle I_{{\mathrm{Na}}^{+}}=\left(g_{\mathrm{Na}}^l+{\bar{g}}_{\mathrm{Na}}^g{m}^{3}h\right)\cdot \left(V-E_{\mathrm{Na}}\right)+3I_p,}$ (15)${\displaystyle I_{{\mathrm{K}}^{+}}=\left(g_{\mathrm{K}}^l+{\bar{g}}_{\mathrm{K}}^g{n}^4\right)\cdot \left(V-E_{\mathrm{K}}\right)-2I_p.}$ Furthermore, the SD model uses dynamic ion concentrations to be able to model the breakdown of the ion gradients that is observed during SD. The intracellular potassium concentration K_i and extracellular potassium concentration K_e are modeled explicitly as dynamical variables, while the intra- and extracellular sodium concentrations (Na_i and Na_e) are computed from the potassium concentration due to the constraint of electroneutrality (16)${\displaystyle \frac{\mathrm{d}{\mathrm{K}}_i}{\mathrm{d}t}=-\frac{\gamma }{{\omega }_i}{I}_{{\mathrm{K}}^{+}},}$ (17)${\displaystyle \frac{\mathrm{d}{\mathrm{K}}_e}{\mathrm{d}t}=\frac{\gamma }{{\omega }_e}{I}_{{\mathrm{K}}^{+}}+J_{\mathrm{diff}}\left({\mathrm{K}}_e\right)}$ (18)${\displaystyle {\mathrm{Na}}_i={\mathrm{Na}}_i^{\left(0\right)}-{\mathrm{K}}_i+{\mathrm{K}}_i^{\left(0\right)},}$ (19)${\displaystyle {\mathrm{Na}}_e=\frac{{\omega }_i}{{\omega }_e}\left({\mathrm{Na}}_i^{\left(0\right)}-{\mathrm{Na}}_i\right)+{\mathrm{Na}}_e^{\left(0\right)}.}$ The factor γ converts currents to ion fluxes and depends on the membrane surface A_m and Faraday’s constant F: (20)${\displaystyle \gamma =\frac{A_m}{F},}$ ω_i and ω_e are constants describing the intra- and extracellular volume, respectively, and the buffer flux J_diff is (21)${\displaystyle J_{\mathrm{diff}}=F_{\mathrm{diff}}\left({\mathrm{K}}_{\mathrm{bath}}-{\mathrm{K}}_e\right).}$ An overview of all constants and the values that were used in the simulations can be found in Table 2.

If ion concentrations are time-dependent, they actually change drastically during neuronal activity. To still maintain homeostasis an ion pump has to be included that pumps Na⁺ ions out of and K⁺ ions into the cell at a 3/2 ratio. The pump current thus depends on the extracellular potassium and the intracellular sodium concentration. The pump is modeled according to Barreto & Cressman (2011) (22)${\displaystyle I_p\left({\mathrm{Na}}_i,{\mathrm{K}}_e\right)=\rho {\left(1+exp\left(\frac{25-{\mathrm{Na}}_i}{3}\right)\right)}^{-1}{\left(1+exp\left(5.5-{\mathrm{K}}_e\right)\right)}^{-1},}$ with ρ being the pump current strength. Note that the pump current also shows up in the equations for Na⁺- and K⁺-currents (Eqs. (14) and (15)).

As a result of the dynamic ion concentrations also the reversal potentials become dynamic (23)${\displaystyle E_{\mathrm{ion}}=\frac{26.64}{z_{\mathrm{ion}}}ln\left({\left[\mathrm{ion}\right]}_e/{\left[\mathrm{ion}\right]}_i\right).}$

The fast gating dynamics of the m-gate is modeled adiabatically as (24)${\displaystyle m=m_{\infty }\left(V\right).}$ Note that in this model shifted versions of the rate equations are used (Cressman et al., 2009; Cressman et al., 2011) (25)${\displaystyle {\alpha }_m=\frac{0.1\left(V+30\right)}{1-exp\left(-\left(V+30\right)/10\right)},}$ (26)${\displaystyle {\beta }_m=4exp\left(-\left(V+55\right)/18\right),}$ (27)${\displaystyle {\alpha }_n=\frac{0.01\left(V+34\right)}{1-exp\left(-\left(V+34\right)/10\right)},}$ (28)${\displaystyle {\beta }_n=0.125exp\left(-\left(V+44\right)/80\right),}$ (29)${\displaystyle {\alpha }_h=0.07exp\left(-\left(V+44\right)/20\right),}$ (30)${\displaystyle {\beta }_h=\frac{1}{1+exp\left(-0.1\left(V+14\right)\right)}.}$ Furthermore, the time constants are scaled by a factor ϕ (31)${\displaystyle {\tau }_x=\frac{1}{\varphi \left({\alpha }_x+{\beta }_x\right)}.}$ In contrast to Hübel, Schöll & Dahlem (2014) we did not reduce the dimension of the model further by assuming a linear or sigmoidal relation between n and h. Instead, h was kept dynamic since the changes caused by the mutation affect the h-gate.

Anoxia model

As a test of the robustness of our results we investigate the effects of FHM3 also in a mutant model of anoxia (Zandt et al., 2011). In fact, migraine with aura has been linked to a higher risk of ischemic stroke (Kurth & Diener, 2012). For furthere details on the rationale, see the Sec. Results.

The anoxia model is similar to the SD model, but uses five more dynamic variables, in particular, it also models chloride ion dynamics. The other dimensions are due to explicitly modeling intra- and extracellular ion concentrations and not assuming mass conservation, and also electroneutrality is not assumed in this model.

Therefore, in addition to Na⁺- and K⁺-currents as in Eqs. (14) and (15) a chloride (Cl⁻) channel is included, which contributes to the leak current (32)${\displaystyle \frac{\mathrm{d}V}{\mathrm{d}t}=-\frac{1}{C_m}\left(I_{{\mathrm{Na}}^{+}}+I_{{\mathrm{K}}^{+}}+I_{\mathrm{Cl}}\right)}$ (33)${\displaystyle I_{{\mathrm{Cl}}^{-}}=g_{\mathrm{Cl}}^l\left(V-E_{\mathrm{Cl}}\right).}$

Intra- and extracellular ion concentrations are dynamic and modeled as (34)${\displaystyle \frac{\mathrm{d}{\mathrm{Na}}_i}{\mathrm{d}t}=-\frac{A}{VF}{I}_{{\mathrm{Na}}^{+}}}$ (35)${\displaystyle \frac{\mathrm{d}{\mathrm{Na}}_e}{\mathrm{d}t}=\frac{\beta A}{VF}{I}_{{\mathrm{Na}}^{+}}}$ (36)${\displaystyle \frac{\mathrm{d}{\mathrm{Cl}}_i}{\mathrm{d}t}=-\frac{A}{VF}{I}_{{\mathrm{Cl}}^{-}}}$ (37)${\displaystyle \frac{\mathrm{d}{\mathrm{Cl}}_e}{\mathrm{d}t}=\frac{\beta A}{VF}{I}_{{\mathrm{Cl}}^{-}}}$ (38)${\displaystyle \frac{\mathrm{d}{\mathrm{K}}_i}{\mathrm{d}t}=-\frac{A}{VF}{I}_{{\mathrm{K}}^{+}}}$ (39)${\displaystyle \frac{\mathrm{d}{\mathrm{K}}_e}{\mathrm{d}t}=\frac{\beta A}{VF}{I}_{{\mathrm{K}}^{+}}-I_g-I_d.}$

The same pump current as in the SD model is used (Eq. (22)). While the total amount of sodium and chloride is constant, the extracellular potassium concentration can be buffered by glial cells (I_g) and diffuse into and out of the blood (I_d) (40)${\displaystyle I_g=G{\left(1+exp\left(\frac{18-{\mathrm{K}}_e}{2.5}\right)\right)}^{-1},}$ (41)${\displaystyle I_d=ϵ\left({\mathrm{K}}_e-k_{\infty }\right),}$ h and n are dynamic and given by Eqs. (5), (6) and (25)–(31). The sodium activation gate m is adiabatically modeled as in Eq. (24). For parameter values see Table 3.

Under physiological conditions this model behaves normally, as it responds with a single action potential to a short current pulse and with periodic firing when a larger current of 1.5 mA/cm² or more is injected (not shown). This model is also able to show seizure activity (Cressman et al., 2009; Cressman et al., 2011).

Modified time constant function based on tail currents

The three models introduced above are given in their ‘wild-type’ formulation. The ‘mutant’ formulation has only a single difference, a modified I_Na current, as described in the following and illustrated in Fig. 1.

From experimental data we know that the mutation leads to a two- to four-fold faster deinactivation (Dichgans et al., 2005) and to a two- to four-fold slower inactivation (Vanmolkot et al., 2007). We checked the robustness of our simulations within this range. The simulations presented here, however, were performed at an intermediate value of a three-fold change.

To change the responsiveness of inactivation and deinactivation accordingly, we need to modify the time constant τ_h of the gating variable h. In the mutant model this time constant is replaced by (42)${\displaystyle {\tau }_h^{\ast }\left(V\right)={\tau }_h\left(V\right)\cdot \left({\kappa }_1\cdot tanh\left(\sigma \cdot \left(\mathrm{V}-{\mathrm{V}}_{\max }\right)\right)+{\kappa }_2\right).}$ The parameter V_max shifts the sigmoidal tanh-function to the position of the maximum of the time constant function τ_h(V). The slope factor of the sigmoidal tanh-function is σ = 0.1 to ensure sufficiently rapid convergence to the limit of a three-fold change. The other parameters are κ₁ = 1.335 and κ₂ = 1.665. These parameters result from the two constraints κ₁ + κ₂ = f and κ₂−κ₁ = 1/f for an f-fold change. We chose f = 3.

To test the mutant time constant ${\tau }_h^{\ast }$, we simulated the experimental protocol performed by Dichgans et al. (2005) in the computational model. The membrane voltage is clamped to a holding potential of −120 mV and then stepped to a potential of −10 mV. At −120 mV the h-gate is completely deinactivated, i.e., open. The step to −10 mV causes the h-gate to inactivate. Therefore, we can measure the time constant of inactivation with this protocol (see Fig. 1B). In contrast, holding the membrane potential at −10 mV and then stepping back to −120 mV allows us to measure the time constant of the process of deinactivation. At −10 mV the h-gate is completely inactivated, i.e., closed, and the step to −120 mV causes the gate to deinactivate again, i.e., the gate reopens. An illustration of this protocol can be found in Fig. 1C. By using this procedure and measuring the two different time constants, it was assured that the chosen parameters lead to a 3-fold slower inactivation and a 3-fold faster deinactivation. The main part of Fig. 1 shows the inactivation time constants τ_h (black line) and ${\tau }_h^{\ast }$ (green line) for the wild-type and mutant model, respectively, as a function of the membrane potential V.

Note that in the Hodgkin–Huxley formalism, the gating subunits of a channel are assumed to be identical and the inactivation and deinactivation as being independent. Therefore this formalism cannot represent certain dependencies in a straightforward manner in the kinetic states. For example, the inactivation of the Na⁺ channel (represented by the h-subunit) has a greater probability of occurring when all subunits are open, therefore the inactivation depends on activation (represented by the three m-subunits). This violates the assumption of independent gating. Because of this independence in the HH formulation, the dynamics of the h-gate is only described by a single time constant function τ_h. An alternative ansatz is to use a Markov model to compute the occupancy of the channel in its various kinetic states as done by Clancy & Rudy (1999).

Mutant AP with marked plateau, increased responsiveness, delayed excitation block, and firing onset unchanged

We first consider the shape of APs. The AP is rather directly affected by FHM3 through altered voltage gating in h. In other words, the results are consistent with the measured tail currents and therefore the results for a mutant AP are even to some degree predictable. This situation will change, when we model dynamics separated three orders of magnitude from AP dynamics.

For a single AP stimulated by a transient applied current I_app(t) of 3 ms duration and 3 µA cm⁻² amplitude (labeled ‘excitatory’ in Fig. 2), we observe that the mutant model compared to wild-type model leads to a prolonged AP with a marked plateau. This is consistent with the larger inactivation time constant ${\tau }_h^{\ast }\left(V^{\mathrm{dep}}\right)$ of the mutant as compared to the wild-type inactivation time scale τ_h(V^dep), cf. tail currents in Fig. 1B. Note that we omitted before the explicit voltage dependency of the time constants, but now we make the dependency explicit because the mutant time constant function ${\tau }_h^{\ast }\left(V\right)$ is in FHM3 increased only for the regime of the membrane potential V being depolarised. This voltage regime is indicated by the superscript “dep” and it corresponds to an inactivation of h (closed h gates).

Furthermore and a bit more subtle to observe, the mutant dynamics reacts faster to a sudden brief stimulation. The mutant model fires an AP that reaches its maximal amplitude just below 2 ms after the I_app is turned off again, while in the wild-type model the maximal amplitude is reached only after about 3 ms. Again, this is also consistent with the defect in the time constant function ${\tau }_h^{\ast }\left(V\right)$. In this case it is explained by the decreased and therefore faster regime τ_h(V^pol) compared to the wild-type. The mutant time constant function ${\tau }_h^{\ast }\left(V\right)$ is decreased for V being in the polarised resting state indicated by the superscript “pol”, cf. tail currents in Fig. 1A. This is the regime of deinactivation (open h gates).

The modified AP profile is also observed during spiking, i.e., in the oscillatory regime, when a constant I_app larger than—by definition (see below)—the rheobase current I_rh is applied. Individual APs in the spike train show this plateau (labeled ‘oscillatory’ in Fig. 2). As a result the spiking frequency is reduced in the mutant model, despite the overall increased responsiveness (Fig. 3). This decreased spiking frequency can be associated with hypoexcitability as the neural response is usually characterized by the firing-rate function.

To get some further quantitative measures of the effects of FHM3 with regard to excitability, we investigated the change of stability in the resting state by varying the input current I_app. This is a bifurcation analysis (Fig. 2). The determined two so-called bifurcation points mark the beginning and end of the oscillatory spiking regime. The first Hopf bifurcation point (HB1) is the onset of oscillation at a minimal value of I_app, which is the definition of the rheobase current I_rh. For the wild-type model the first Hopf bifurcation (HB1) is at $I_{\mathrm{app}}^{\mathrm{HB1}}\equiv I_{\mathrm{rh}}=9.78\mu \mathrm{A}\hspace{0.17em}{\mathrm{cm}}^{-2}$ and the second Hopf bifurcation (HB2) at $I_{\mathrm{app}}^{\mathrm{HB2}}=154.5\mu \mathrm{A}\hspace{0.17em}{\mathrm{cm}}^{-2}$, which determines the excitation block as the oscillation ceases at this point. For the mutant model these Hopf bifurcations occur at I_rh = 9.72 µA cm⁻² and $I_{\mathrm{app}}^{\mathrm{HB2}}=175.0\mu \mathrm{A}\hspace{0.17em}{\mathrm{cm}}^{-2}$. The first Hopf bifurcations (HB1) are subcritical, while HB2 are both supercritical. This means that if the I_app is not slowly ramped towards the rheobase current I_rh, one can observe the oscillatory regime even before the I_rh. Hence the two firing-rate functions in Fig. 3 start slightly before the values given here for HB1, with the mutant model starting again earlier.

With regard to the rheobase current, the values for the wild-type and mutant differ by less then 0.6%, with the mutant value being smaller, which, at least in principal, corresponds to hyperexcitability, though due to the small magnitude this seems negligible for all practical purposes. However, the excitation block observed at the second critical transition HB2 occurs at larger values of I_app for the mutant model. The mutant channels tolerate an increased maximal $I_{\mathrm{app}}^{\mathrm{HB2}}$ by 13% compared to the wild-type. This means that the mutant neurons exhibit oscillatory behavior in a larger range of applied currents. Therefore, this shift establishes a gain-of-function, which indicates hyperexcitability.

To summarize, while the reduced firing frequency indicates hypoexcitability, increased responsiveness and delayed excitation block indicate hyperexcitability.

Mutant more vulnerable to SD

We now focus on effects of FHM3 upon cellular functioning that occurs in the same neural substrate that generates APs but on time scales at least three orders of magnitude separated from AP dynamics, that is, effects that occur during several tens of seconds up to minutes. This is the time scale of SD. It is therefore relevant for pathological conditions, for instance, in migraine with aura. In accordance with this pathophysiological context, we select the stimulations of SD in the wild-type and mutant model as rather large perturbations to neural homeostasis such as a compromised energy supply during focal hypoperfusion that induces and occurs in conjunction with migraine aura symptoms (Olesen et al., 1993; Friberg et al., 1994).

In particular, we investigate the effect of a breakdown of the Na⁺-K⁺-pump upon the membrane potential V and reversal potentials E_Na and E_K. For this purpose the maximum pump rate ρ is linearly down-regulated to 20% of its physiological value within 10 s, then ρ is kept at 20% for a variable time window, and finally ρ is linearly up-regulated back to 100% within 5 s. The stimulation trace of ρ is shown in Fig. 4 with the dashed-dotted line. The specific choice of the variable time window is additionally marked for the wild-type and mutant stimulation trace by an annotated two-headed arrow. Let us remark that in our studies we also used two other perturbations, namely a transient increase in extracellular K⁺ concentration (by increasing K_bath) and a large current pulse I_app, with basically the same results (not shown).

We determined the minimal duration of the variable time window with reduced pump rate (20%) that is just no longer tolerated and results in a long lasting but transient breakdown of the reversal potentials E_Na and E_K characteristic for SD. For this purpose we increased the variable time window by 0.1 s steps. While the wild-type model could not tolerate a period of 13.6 s of reduced pump rate at 20%, the mutant model was less robust and could not tolerate a period of 7.2 s of reduced pump activity (Fig. 4). Therefore, the mutant model is approximately only half (53%) as robust to periods of reduced ion pump activity as the wild-type model is.

Shorter stimulation periods did not lead to full blown SD signals. In this case, the spiking ceased about a second after the interval began that increased the pump rate back from 20% to 100% (this interval lasts 5 s) and, more importantly, both membrane potential V and reversal potentials E_Na and E_K recovered within only a few seconds back to physiological values (not shown). Thus, SD profiles of these potentials, which followed longer stimulation periods, are clearly distinguished by a all-or-none phenomenon. Not only do membrane potential V and reversal potentials E_Na and E_K change dramatically after the stimulation is off, but also full recovery from SD to the initial physiological values takes very long. Of course recovery reaches the resting state only asymptotically. For up to one to two hours the changes in particular in E_Na are observable, while the signals in Fig. 4 are shown only for 100 s. It is noteworthy that the neuronal state is already back to basic functioning emitting APs if stimulated after the repolarization, that is, even if the resting state is not fully recovered. Similar dynamics is described in other computational models of SD by Kager, Wadman & Somjen (2000), Yao, Huang & Miura (2011) and Hübel, Schöll & Dahlem (2014).

To summarize, in terms of susceptibility to SD the mutant model is hyperexcitable. This seems to be in contrast to the major effect of the mutant upon the AP firing frequency that indicates that the mutant model is hypoexcitable (Fig. 3). This will be further discussed in the Discussion.

Effects in anoxia model consistent with SD model

Last, we study a model of AD (Zandt et al., 2011); the AD model shares many features with the SD model but is more detailed (see Methods) and hence effects obtained with this model serve as control to compare them with effects obtained from the SD model. The model was first published to study slow waves after decapitation in a computational model (Zandt et al., 2011). By repeating this with a mutant version of this model, our focus is set very similar to the previous section. In the decapitation study, anoxia is modeled by completely switching off all pump, glial, and diffusion currents, see Fig. 5. In fact, Fig. 5A with the wild-type model is a reproduction of the simulations performed by Zandt et al. (2011).

Note that patients with migraine with aura are at greater risk for stroke (Kurth & Diener, 2012). Thus there is a rationale to perform this comparison beyond the mere confirmation of plausibility of our results obtained above with the SD model. However, the multiplicity of potential links include not only common genetic risk factors but also indirect links like common triggers outside the brain, e.g., microemboli caused by cardiac shunts. Furthermore, the model investigated by Zandt et al. (2011) is derived from a model suggested by Cressman et al. (2009). This model exhibits periodic bursting similar to seizure activity. Both migraine and epilepsy have genetically based forms caused by various mutations in genes, while the mutation in FHM3 differs markedly within the several mutations in SCN1A therein that it is not associated with epilepsy (see Introduction). Investigating the underlying ion homeostasis in the three conditions of epilepsy, migraine, and stroke may yield interesting results in future investigations of computational models that can unify certain dynamical aspects and link disease genotype to phenotype. However, this is clearly beyond the scope of this study.

In this study, let us only refer to the dynamics resulting from switching off all pump, glial, and diffusion currents until the excitation block and compare the wild-type and mutant model. After a gradual rise of the membrane potential that lasts in either case about 30 s (note that the simulated ‘decapitation’ occurs in Fig. 5 at t = 5 s), the membrane potential reaches the AP threshold, subsequently resulting in a final burst of spiking. These initial, less than a minute lasting, phases in the wild-type and mutant model are indeed very similar to the initial phases in the SD model following a transient energy failure. A minor difference is that the gradual rise is overall slower, but this is explained by a slightly different geometry (larger extracellular space) and by the chloride ion dynamics (Hübel, Schöll & Dahlem, 2014). The similarity supports the robustness of our results, as this model is an established model showing anoxia (Zandt et al., 2011) and seizure activity (Cressman et al., 2009; Cressman et al., 2011).

To summarize, also for AD the slow gradual fall of the potentials does not significantly differ during the initial leak phase in the wild-type and mutant model, while once the model is spiking the excitation block occurs about 2.5-times faster, corresponding to a faster breakdown of ion gradients due to spiking, in the mutant model.

Hypoexcitable vs. hyperexcitable

The increased susceptiblility to SD does not contradict the reduced firing frequency for a given stimulation current I_app, although this change in firing frequency indicates that the mutant model is hypoexcitable.

Firing a single action potential (AP) is a form of cellular excitability manifested as a transmembrane voltage jump without significant changes in ion concentrations. SD is a form of cellular excitability manifested by massive changes in ion concentrations. There is not necessarily a direct relation between the two excitable systems, not even with regard to merely classifying terms such as hypoexcitable and hyperexcitable. Rather, AP and SD can be viewed as largely independent phenomena, because while sharing the same neural substrate, AP and SD are separated by time scales differing in three orders of magnitude (see below). Notwithstanding, the massive breakdown of ion gradients in SD is, of course, mediated by APs that occur on the fast time scale.

In our view, “hypoexcitable” vs. “hyperexcitable” is in any case too simple a classification scheme even considering AP and SD in isolation on their respective time scale. To support this criticism of classifying neural dynamics in migraine, let us mention that this problem was also addressed in the psychophysical and clinical contexts, see studies by Shepherd (2001) and Coppola, Pierelli & Schoenen (2007) and references therein; further support comes from the mathematical picture (below)—which are two sides of the same coin.

To illustrate this with only a single example, consider, as already mentioned above, that the mutant channels exhibit an increased range of spiking activity with a delayed excitation block by 13% compared to the wild-type. We argued that this larger spiking range establishes a gain-of-function. Consider further the increased responsiveness of the mutant model. Both indicate a form of hyperexcitability with regard to AP. In contrast, the change in firing frequency of AP indicates at the same time that the mutant model is hypoexcitable (Fig. 3).

SD susceptibility

How do these three diverse effects observed for APs (delayed excitation block, increased responsiveness, and lower firing frequency) manifest on the longer time scale under the condition of SD?

In terms of susceptibility to SD, the shifted excitation block (see HB2 in Fig. 2) might misleadingly suggest that the mutant model is less susceptible to SD. This is similar to the lower firing frequency that we considered above. Since the characteristic sustained breakdown of the reversal potentials E_Na and E_K is ignited in our model only if the system is driven by any stimulation into the excitation block, its delay in the mutant model seems to suggest that a longer stimulation may be needed and therefore a higher threshold exists.

To show the actual situation in Fig. 4, we highlighted a critical time window by a gray shade. This critical time window opens with start of the reduced pump rate recovery (from 20% back to 100%) and it closes with the beginning of the excitation block. Considering only the delay of the excitation block and the low frequency, it may seem surprising at first, this critical period lasts 3.4 s in the wild-type model and only 2.5 s in the mutant model. Note that this ‘paradox’ can also be observed in the overall shorter duration of the whole initial firing pattern in the mutant SD model. Our attention should be on signals that can actually be measured in a clinical setting, hence our focus is on these signals also in the presentation of the computational model, where we can “measure” everything. The reduced pump rate corresponds to hypoperfusion signals. The excitation block in SD corresponds to the first peak in an electroencephalography (EEG) signal, cf. the work by Zandt et al. (2011) where the simulated membrane potential is averaged and high-pass filtered, cut-off at 0.1 Hz, to estimate the EEG—although this EEG might only be observable intracranially.

That the mutant model is more susceptible to spreading depression (SD) is exclusively explained by the much larger amount of ions transferred across the membrane during spiking. This, in particular the intracellular ion concentration, cannot easily be measured even in an in vitro setup. In Fig. 4, we see this by the much steeper slope of the reversal potential E_K in the mutant model.

The multidimensional concept of thresholds

The complex question about susceptibility requires a deeper understanding of what a threshold is. In fact, the very reason why we have to get beyond the idea of “hypoexcitable” vs. “hyperexcitable” as a useful characterization of the system (see above) is that there is no one-dimensional ansatz to determine a threshold as a demarcation.

Before explaining this further, let us give one more explicit example. In other model variants of SD (Kager, Wadman & Somjen, 2000), a stimulation of SD may even stop before the excitation block is reached. In this case a sustained afterdischarge carries the system into the depolarisation block that then marks the start of the actual SD events. Clearly, in this case the depolarisation block cannot be considered being the actual threshold, because the system is ‘before’ this point when the stimulation is already off again.

In general, excitability or all-or-none phenomena do not possess a threshold in terms of single quantity, whether it is a particular membrane depolarisation that demarcates the all-or-none response in the form of an AP or a critical duration of hypoperfusion that demarcates the all-or-none response in form of SD. A detailed analysis of neural models shows that a threshold is a multidimensional surface (manifold) not a single number as first shown by FitzHugh (1955) and as discussed in a modern style by Mitry et al. (2013) and applied to migraine by Dahlem (2013). So the actual use of computational models goes far beyond numerical simulations. We gain a deeper understanding of the principal mechanisms in precise mathematical relationships, of which we can only give a very general overview in this paper.