\section*{Introduction}


Familial hemiplegic migraine (FHM) is a rare monogenic, autosomal dominantly inherited syndrome with hemiparesis during the aura phase of migraine. Three distinct genetic mutations for FHM have been identified, in the CACNA1A calcium channel gene (FHM1), in the ATP1A2 Na,K-ATPase gene (FHM2), and in the SCN1A sodium channel gene (FHM3).  It has been proposed that all three phenotypes reflect hyperexcitability in the form of increased susceptibility for  spreading depression (SD). However, the functional connection between the molecular findings and a facilitated   generation of SD is unclear.  

To determine the electrophysiological consequences of such a genetic defect, we integrate a mutation of FHM3 into three types of computational models of neuronal dynamics. This allows us to bridge the gap between genotype and phenotype. A similar approach was used by \cite{CLA99}. We use a standard Hodgkin-Huxley model for action potentials (AP) \citep{HOD52} and a model of SD \citep{HUE14} to evaluate the change in the threshold of generating SD by tolerating various brief intervals of transient ischemic attacks. Moreover, we use a model for anoxic deporlarization (AD) \citep{ZAN11} that is derived form a seizure model \citep{CRE09} as a test of the robustness of our results.     

The paper is organized as follows. In the Sec.\ Methods we introduce three computational models and our method to incorporate measured tail currents in FHM3 \citep{DIC05,VAN07} into the Hodgkin-Huxley framework. In Sec.\ Results we present simulations and analysis of the wild-type and mutant models. We end with Sec.\ Discussion were we focus on three topics: (i) the appropriateness of the terms hypoexcitable vs. hyperexcitable, (ii) the seemingly paradoxically increased susceptibility to SD in the mutant model if one considers the firing rate, a measure that is usually used to quantify slow neural dynamics, and (iii) the inadequate concept of a threshold as a quantity measured by a single value.   

\section*{Methods}\label{sec:methods}

All three models are based on Hodgkin-Huxley type dynamics with different degree of complexity from the classical model to a second generation with time-dependent ion concentrations. 

\subsection*{Hodgkin-Huxley model}\label{sec:HH}

\begin{table}[b]
\caption{\label{tab:1} Model parameters for Hodgkin-Huxley model.}
\begin{tabular}{|l|l|l|}
\hline
Name 			& Value \& unit 		& Description \\ 
\hline
$C_m$			& 1 $\mu$F/cm$^2$ 		& membrane capacitance \\
$\bar{g}_{Na}$		& 120 mS/cm$^2$ 		& max. sodium conductance \\
$\bar{g}_K$			& 36 mS/cm$^2$ 			& max. potassium leak conductance \\
$g_l$			& 0.3 mS/cm$^2$ 		& leak conductance \\
$E_{Na}$		& 50 mV					& sodium reversal potential \\
$E_K$			& -77 mV 				& potassium reversal potential \\
$E_{leak}$			& -54.402 mV 			& leak reversal potential \\
\hline
\end{tabular}
\end{table}


The Hodgkin-Huxley (HH) model is one of the most widely used computational models in neurscience. It is a conductance-based neuron model \citep{HOD52} 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 \in \{K,Na,leak\}$. Note that  the conductance  $g_j$ for voltage-gated channels, i.e., $j \in \{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}$. 
\begin{eqnarray}
\frac{\mathrm{d}V}{\mathrm{d}t} &=& -\frac{1}{C_m}(I_{Na^+}+I_{K^+} + I_{leak} - I_{app}), 	\label{eq:1}\\
I_{Na^+} &=& \bar{g}_{Na} m^3 h \left ( V - E_{Na} \right ),										\label{eq:2}\\
I_{K^+} &=& \bar{g}_K n^4 \left ( V - E_K \right ), 												\label{eq:3}\\
I_{leak} &=& g_{l} \left ( V - E_{leak} \right ). 												\label{eq:4}
\end{eqnarray}
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:
\begin{eqnarray}
\frac{\mathrm{d}x}{\mathrm{d}t}				&=& 	\frac{x_\infty - x}{\tau_x } \quad\mathrm{with} 		\label{eq:5}\\
x_\infty 	   								&=& 	\frac{\alpha_x}{\alpha_x + \beta_x}\quad\mathrm{and}	\label{eq:6}\\
\tau_x 						   				&=& 	\frac{1}{\alpha_x+\beta_x}\quad\mathrm{for}\ x \in \{n,\ m,\ h\}.\label{eq:7}
\end{eqnarray}
$x_\infty$ describes the steady-state of the gating variables and $\tau_x$ is the time constant. 

The rate equations for $\alpha_x$ and $\beta_x$ are voltage-dependent and given by
\begin{eqnarray}
\alpha_m	&=& \frac{0{.}1(V+40)}{1 - \exp(-(V+40)/10)}\ , 	\label{eq:8}\\
\beta_m		&=& 4 \exp(-(V+65)/18)\ , 							\label{eq:9}\\
\alpha_n	&=& \frac{0{.}01(V+55)}{1 - \exp(-(V+55)/10)}\ , 	\label{eq:10}\\
\beta_n		&=& 0{.}125 \exp(-(V+65)/80)\ , 					\label{eq:11}\\
\alpha_h	&=& 0{.}07 \exp(-(V + 65)/20))\ ,					\label{eq:12}\\
\beta_h 	&=& \frac{1}{1 + \exp(-0{.}1(V + 35))}\ .			\label{eq:13}
\end{eqnarray}
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 \ref{tab: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 bistabilty, which is consistant with our results of a prolognoed plateau of action potential and early depolarization block in form of bistability.

\subsection*{Spreading depression model} \label{sec:sd}

\begin{table}[b!]
\caption{\label{tab:2} Model parameters for the SD model}
\begin{tabular}{|l|l|l|}
\hline
Name 			& Value \& unit 												& Description \\
\hline
$C_m$			& 1 $\mu$F/cm$^2$ 		& membrane capacitance \\
$g_{Na}^l$		& 0{.}0175 mS/cm$^2$ 	& sodium leak conductance \\
$g_{Na}^g$		& 100 mS/cm$^2$ 		& max.\ gated sodium conductance \\
$g_K^l$			& 0{.}05 mS/cm$^2$ 		& potassium leak conductance \\
$g_K^g$			& 40 mS/cm$^2$ 			& max.\ gated potassium conductance \\
$\mathit{Na}_i$	& 27 mMol/$l$ 			& ECS sodium concentration \\
$\mathit{Na}_e$	& 120 mMol/$l$ 			& ICS sodium concentration \\
$K_i$			& 130{.}99 mMol/$l$		& ECS potassium concentration \\
$K_e$			& 4 mMol/$l$			& ICS potassium concentration \\
$E_{Na}$		& 39{.}74 mV			& sodium reversal potential \\
$E_K$			& -92{.}94 mV 			& potassium reversal potential \\
$\omega_i$ 		& 2{.}16 $\mu$m$^3$ 											& volume of ICS \\ 
$\omega_e$		& 0{.}72 $\mu$m$^3$ 											& volume of ECS \\
$F$ 			& 96485 C/Mol	 												& Faraday's constant \\ 
$A_m$ 			& 0{.}922 $\mu$m$^2$ 											& membrane surface \\ 
$\gamma$ 		& 9{.}556e-6 $\frac{\mu\mathrm{m}^2\mathrm{Mol}}{\mathrm{C}}$ 	& conversion factor \\ 
$\rho$ 			& 5{.}25 $\mu$A/cm$^2$ 											& max.\ pump current \\
$\phi$ 			& 3/msec 														& gating timescale parameter\\ 
$F_{diff}$		& 6{.}66e-6/msec													& diffusion parameter		\\
$K_{bath}$		& 4 mMol/$l$													& potassium bath concentration \\
\hline
\end{tabular}
\end{table}


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 boththe 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 \cite{HUE14}.

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 
\begin{eqnarray}
  I_{Na^+} &=& (g_{Na}^l + \bar{g}_{Na}^g m^3h) \cdot (V - E_{Na}) + 3 I_p\ , 	\label{eq:14} \\
  I_{K^+}  &=& (g_K^l + \bar{g}_K^g n^4) \cdot (V - E_K) - 2 I_p.\ 				\label{eq:15}
\end{eqnarray}
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
\begin{eqnarray}
\frac{\mathrm{d}K_i}{\mathrm{d}t}	&=&	-\frac{\gamma}{\omega_i}I_{K^+} , \label{eq:16}\\
\frac{\mathrm{d}K_e}{\mathrm{d}t}			&=& \frac{\gamma}{\omega_e}I_{K^+} + J_{diff}(K_e) 	\label{eq:21}\\
\mathit{Na}_i 						&=& \mathit{Na}_i^{(0)}-K_i+K_i^{(0)} , \label{eq:17}\\
\mathit{Na}_e 						&=& \frac{\omega_i}{\omega_e}(\mathit{Na}_i^{(0)}- \mathit{Na}_i) + \mathit{Na}_e^{(0)}. \label{eq:18}
\end{eqnarray}
The factor $\gamma$ converts currents to ion fluxes and depends on the membrane surface $A_m$ and Faraday's constant $F$:
\begin{eqnarray}
\gamma = \frac{A_m}{F},	\label{eq:19}
\end{eqnarray}
$\omega_i$ and $\omega_e$ are constants describing the intra- and extracellular volume, respectively, and the buffer flux $J_{diff}$ is 
\begin{eqnarray}
J_{diff} &=& F_{diff} (K_{bath}-K_e)  		\label{eq:22}	
\end{eqnarray}
An overview of all constants and the values that were used in the simulations can be found in Table \ref{tab: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 \cite{CRE09,BAR11}
\begin{eqnarray}
I_p(\mathit{Na}_i,K_e) &=& \rho\bigg(1+\exp{\bigg(\frac{25-\mathit{Na}_i}{3}\bigg)}\bigg)^{-1}  \bigg(1+\exp{(5.5-K_e)}\bigg)^{-1}\ , \label{eq:20}
\end{eqnarray}
with $\rho$ being the pump current strength. Note that the pump current also shows up in the equations for $Na^+$- and $Na^+$-currents (Eqs.~(\ref{eq:14}) and (\ref{eq:15})).


As a result of the dynamic ion concentrations also the reversal potentials become dynamic
\begin{eqnarray}
E_{ion}=\frac{26{.}64}{z_{ion}}\ln(\mathit{[ion]}_e/\mathit{[ion]}_i).	\label{eq:23}
\end{eqnarray}

The fast gating dynamics of the $m$-gate is modeled adiabatically as
\begin{eqnarray}
m=m_\infty(V)\ . \label{eq:24}
\end{eqnarray}
Note that in this model shifted versions of the rate equations are used \cite{CRE09}
\begin{eqnarray}
\alpha_m	&=& \frac{0{.}1(V+30)}{1 - \exp(-(V+30)/10)}\ , 	\label{eq:25}\\
\beta_m		&=& 4 \exp(-(V+55)/18)\ , 							\label{eq:26}\\
\alpha_n	&=& \frac{0{.}01(V+34)}{1 - \exp(-(V+34)/10)}\ , 	\label{eq:27}\\
\beta_n		&=& 0{.}125 \exp(-(V+44)/80)\ , 					\label{eq:28}\\
\alpha_h	&=& 0{.}07 \exp(-(V + 44)/20))\ ,					\label{eq:29}\\
\beta_h 	&=& \frac{1}{1 + \exp(-0{.}1(V + 14))}\ .			\label{eq:30}
\end{eqnarray}
Furthermore, the time constants are scaled by a factor $\phi$
\begin{eqnarray}
\tau_x 		&=& 	\frac{1}{\phi(\alpha_x+\beta_x)}.			\label{eq:31}
\end{eqnarray}
In contrast to \cite{HUE14} 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. 

\subsection*{Anoxia model}\label{sec:waveofdeath}


\begin{table}[b!]
\caption{\label{tab:3} Model parameters for anoxia model}
\begin{tabular}{|l|l|l|}
\hline
Name 			& Value \& unit 	    & Description \\
\hline
$C_m$			& 1 $\mu$F/cm$^2$ 		& membrane capacitance \\
$g_{Na}^l$		& 0{.}0175 mS/cm$^2$ 	& sodium leak conductance \\
$g_{Na}^g$		& 100 mS/cm$^2$ 		& max.\ gated sodium conductance \\
$g_K^l$			& 0{.}05 mS/cm$^2$ 		& potassium leak conductance \\
$g_K^g$			& 40 mS/cm$^2$ 			& max.\ gated potassium conductance \\
$g_{Cl}^l$		& 0{.}05 mS/cm$^2$ 		& chloride leak conductance \\
$\mathit{Na}_i$	& 27 mMol/$l$ 			& ECS sodium concentration \\
$\mathit{Na}_e$	& 120 mMol/$l$ 			& ICS sodium concentration \\
$K_i$			& 130{.}99 mMol/$l$		& ECS potassium concentration \\
$K_e$			& 4 mMol/$l$			& ICS potassium concentration \\
$E_{Na}$		& 39{.}74 mV			& sodium reversal potential \\
$E_K$			& -92{.}94 mV 			& potassium reversal potential \\
$\phi$ 			& 3/msec 	 			& gating timescale parameter\\ 
A/VF			& 0{.}044 $\frac{mMol}{s} / (\frac{mA}{cm^2})$		& conversion factor\\
$\beta$			& 2{.}0					& ratio ICS/ECS \\
$\rho$			& 28{.}1 $\mu$A/cm$^2$  & Na-K-Pump rate \\
G				& 66 mMol/s				& glial buffering rate for K$^+$ \\
$\epsilon$		& 1{.}3 s$^{-1}$ 		& diffusion rate \\
$k_\infty$		& 4{.}0 mMol			& concentration K$^+$ in blood \\
T				& 310 K					& absolute temperature \\
\hline
\end{tabular}
\end{table}


As a test of the robustness of our results we investigate the effects of FHN3 also in a mutant model of anoxia \citep{ZAN11}. In fact, migraine with aura has been linked to a higher risk of ischemic stroke \citep{KUR12}. For furthere details on the rational, 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.~(\ref{eq:14}) and (\ref{eq:15}) a chloride ($Cl^-$) channel is included, which contributes to the leak current

\begin{eqnarray}
\frac{\mathrm{d}V}{\mathrm{d}t} &=& 	-\frac{1}{C_m}(I_{Na^+} + I_{K^+} + I_{Cl}) 	\label{eq:32}\\
I_{Cl^-} 						&=& 	g_{Cl}^l( V - E_{Cl}). \						\label{eq:33}
\end{eqnarray}

Intra- and extracellular ion concentrations are dynamic and modeled as
\begin{eqnarray}
\frac{\mathrm{d}Na_i}{\mathrm{d}t}	&=&	-\frac{A}{VF}I_{Na^+}\ 					\label{eq:34}\\
\frac{\mathrm{d}Na_e}{\mathrm{d}t}	&=&	\frac{\beta A}{VF}I_{Na^+}\ 			\label{eq:35}\\
\frac{\mathrm{d}Cl_i}{\mathrm{d}t}	&=&	-\frac{A}{VF}I_{Cl^-}\ 					\label{eq:36}\\
\frac{\mathrm{d}Cl_e}{\mathrm{d}t}	&=&	\frac{\beta A}{VF}I_{Cl^-}\ 			\label{eq:37}\\
\frac{\mathrm{d}K_i}{\mathrm{d}t}	&=&	-\frac{A}{VF}I_{K^+}\ 					\label{eq:38}\\
\frac{\mathrm{d}K_e}{\mathrm{d}t}	&=&	\frac{\beta A}{VF}I_{K^+} - I_g - I_d\ 	\label{eq:39}
\end{eqnarray}

The same pump current as in the SD model is used (Eq.~(\ref{eq:20})). 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$) 
\begin{eqnarray}
I_g		&=& G \bigg(1+\exp{\bigg(\frac{18-\mathit{K}_e}{2.5}\bigg)}\bigg)^{-1}, 	\label{eq:40}\\
I_d 	&=& \epsilon \left(\mathit{K}_e - k_\infty \right),  					\label{eq:41}
\end{eqnarray}
$h$ and $n$ are dynamic and given by Eqs.~(\ref{eq:5}), (\ref{eq:6}), (\ref{eq:25})-(\ref{eq:31}). The sodium activation gate $m$ is adiabatically modeled as in Eq.~(\ref{eq:24}). 

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^2$ or more is injected (not shown). This model is also able to show seizure activity \citep{CRE09}.


\subsection*{Modified time constant function based on tail currents}
\label{sec:mutation}

\begin{figure}[t!]
\begin{center}
\includegraphics[width=\linewidth]{fhm3Fig1}
\end{center}
\caption{Voltage-dependent time constant for mutation ($\tau_h^*$) and wild-type ($\tau_h$). Insets show the response of $h$ to a voltage-clamp protocol. Left inset shows the deinactivation (i.e., recovery from inactivation) process as a response to a step in voltage from $-10$mV to $-120$mV. Right inset shows the inactivation process by stepping the voltage from $-10$mV to $-120$mV. The intersections of the $h$-curve with the 1/e- and (1-1/e)-lines, respectively, show the actual time constants. In the left inset $\tau_h^*$ is three-fold smaller than $\tau_h$. In the right inset $\tau_h^*$ is three-fold lager than $\tau_h$.\label{fig:methods}}
\end{figure}

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.~\ref{fig:methods}. 

From experimental data we know that the mutation leads to a two- to four-fold faster deinactivation \citep{DIC05} and to a two- to four-fold slower inactivation \cite{VAN07}. 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 $\tau_h$ of the gating variable $h$. In the mutant model this time constant is replaced by 
\begin{eqnarray}
\tau_h^*(V) &=& \tau_h(V) \cdot (\kappa_1 \cdot \tanh(\sigma \cdot (\mathrm{V}-\mathrm{V_{max}})) + \kappa_2) \label{eq:42}
\end{eqnarray}
The parameter $V_{max}$ shifts the  sigmoidal $\tanh$-function to the position of the maximum of the of time constant function $\tau_h$(V). The slope factor of the sigmoidal $\tanh$-function is $\sigma=0.1$ to ensure sufficiently rapid convergence to the limit of a three-fold change.  The other parameters are $\kappa_1=1{.}335$ and $\kappa_2=1{.}665$. These parameters results from the two constrains $\kappa_1+\kappa_2=f$ and $\kappa_2-\kappa_1=1/f$ for an $f$-fold change. We chose $f=3$. 

To test the mutant  time constant $\tau^*_h$, we simulated the experimental protocol performed by \cite{DIC05} in the computational model. The membrane voltage is clamped to  a holding potential of $-120\mathrm{mV}$ and then stepped to a potential of $-10\mathrm{mV}$. At $-120\mathrm{mV}$ the $h$-gate is completely deinactivated, i.e., open. The step to $-10\mathrm{mV}$ causes the $h$-gate to inactivate. Therefore, we can measure the time constant of \textit{inactivation} with this protocol (see right inset of Fig.~\ref{fig:methods}).
In contrast, holding the membrane potential at $-10\mathrm{mV}$ and then stepping back to $-120\mathrm{mV}$ allows us to measure the time constant of the process of \textit{deinactivation}. At $-10\mathrm{mV}$ the $h$-gate is completely inactivated, i.e., closed, and the step to $-120\mathrm{mV}$ causes the gate to deinactivate again, i.e., the gate reopens. An illustration of this protocol can be found in the left inset of Figure \ref{fig:methods}. 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. 



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 $\tau_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 \cite{CLA99}. 



\section*{Results} 
\label{sec:results}

Three different models are investigated, a model of action potentials (AP), a model of spreading depression (SD), and a model of anoxic depolarization (AD). These models describe normal cell functions in terms of the dynamic repertoire either without genetic defect (three wild-type models) or with altered cell functions in FHM3 (three mutant models). The three  mutant models  (AP, SD, and AD) are the same as the wild-type models except that the $I_{Na}$ current has a different voltage-gating mechanism in the fast gating variable $h$. This is described in the wild-type model by the time constant $\tau_h$ and in the mutant model by $\tau^*_h$ (see Sec.~Methods). The observed functional consequences of FHM3 occur on time scales ranging from milliseconds to several tens of seconds. 




\subsection*{Mutant AP with marked plateau, increased responsiveness, delayed excitation block, and firing onset unchanged}
\label{sec:resultsAP}

\begin{figure}[t!]
\begin{center}
\includegraphics[width=0.8\linewidth]{fhm3Fig2}
\end{center}
\caption{Comparison of wild-type (upper panel) and mutant (lower panel) spiking behavior. The main plots show bifurcation diagrams by varying the external current $I_{app}$. For the wild-type model Hopf bifurcations can be found at $I_{app}=9.77994$ and $I_{app}=154.527\mu\mathrm{A\,cm}^{-2}$ . For the mutant model Hopf bifurcations occur at $I_{app}=9.72266\mu\mathrm{A\,cm}^{-2}$ and $I_{app}=175.027\mathrm{A\,cm}^{-2}$. The right insets show the response of the models in the excitatory regime to a $3$ms long current pulse with amplitude $3\mathrm{A\,cm}^{-2}$. The left insets show behavior in the oscillatory regime as a response to a constant input current of $12\mathrm{A\,cm}^{-2}\mu$.  \label{fig:ap}}
\end{figure}

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 3ms duration and $3\mu\mathrm{A\,cm}^{-2}$ amplitude (see inset labeled `excitatory' in Fig.~\ref{fig:ap}), 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^*(V^{\mathrm{dep}})$ of the mutant as compared to the wild-type inactivation time scale $\tau_h(V^{\mathrm{dep}})$, cf.\ tail currents in the right inset of Fig.~\ref{fig:methods}. 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 {\it function} $\tau_h^*(V)$ 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 2ms after the $I_{app}$ is turned off again, while in the wild-type model the maximal amplitude is reached only after about 3ms. Again, this is also consistent with the defect in the time constant function $\tau_h^*(V)$. In this case it is explained by the decreased and therefore faster regime  $\tau_h(V^{\mathrm{pol}})$ compared to the wild-type.  The mutant  time constant function $\tau_h^*(V)$ is decreased for $V$ being in the polarised resting state indicated by the superscript ``pol'', cf.\ tail currents left inset in Fig.~\ref{fig:methods}. 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 (inset labeled `oscillatory' in Fig.~\ref{fig:methods}). As a result the spiking frequency is reduced in the mutant model, despite the overall increased responsiveness (Fig.~\ref{fig:rate}). 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.~\ref{fig:ap}). 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^{HB1}_{app}\equiv I_{rh}=9.78\mu\textrm{A\,cm}^{-2}$ and the second Hopf bifurcation (HB2) at $I^{HB2}_{app}=154.5\mu\textrm{A\,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\mu\textrm{A\,cm}^{-2}$ and $I^{HB2}_{app}=175.0\mu\textrm{A\,cm}^{-2}$. These Hopf bifurcations are subcritical. 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.~\ref{fig:rate} start slightly before the values given here for HB1, with the mutant model starting again earlier.  

\begin{figure}[t!]
\begin{center}
\includegraphics[width=0.8\linewidth]{fhm3Fig3}
\end{center}
\caption{Nonlinear firing-rate function $\mathrm{F}(I_{app})$ for wild-type model (black, solid) and mutant model (green, dashed).  \label{fig:rate}}
\end{figure}

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^{HB2}_{app}$ 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.

\subsection*{Mutant more vulnerable to SD}
\label{sec:resultsSD}

\begin{figure}[b!]
\begin{center}
\includegraphics[width=\linewidth]{fhm3Fig4}
\end{center}
\caption{Development of SD in wild-type (upper panel) and mutant (lower panel) models. A SD is elicited by down-regulating the pump current to 20\% of its maximal value for $13.6$s (wild-type) and $7.2$s (mutant), respectively (see blacked dashed line). The red and blue dashed lines show the temporal development of the sodium and potassium reversal potentials.\label{fig:2}}
\end{figure}

We now focus on effects of FHM3 upon cellular functioning that occur 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 \citep{OLE93,FRI94}. 

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 $\rho$ is linearly down-regulated to $20\%$ of its physiological value within 10s, then $\rho$ is kept at 20\% for a variable time window, and finally $\rho$ is linearly up-regulated back to 100\% within 5s. The stimulation trace of $\rho$ is shown in Fig.~\ref{fig:2} 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 break-down of the reversal potentials $E_{Na}$ and $E_K$ characteristic for SD. For this purpose we increased the variable time window by 0.1s 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.~\ref{fig:2}). 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 5s) 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.~\ref{fig:2} are shown only for $100s$. 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 \cite{KAG00,YAO11,HUE14}. 


\begin{figure}[b!]
\begin{center}
\includegraphics[width=0.75\linewidth]{fhm3Fig5}
\end{center}
\caption{Response of wild-type (upper panel) and mutant (lower panel) membrane potential to a complete breakdown of	 pump, glial and diffusion currents at $t=5s$ (black dashed vertical line). Red and blue dotted lines show the $Na^+$ and $K^+$ reversal potentials over time. The time from the onset of spiking until the beginning of the excitation block is approximately $6.7s$ without and $2.7s$ with mutation. \label{fig:3}}
\end{figure}



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.~\ref{fig:rate}). This will be further discussed in Sec. Discussion.

 


\subsection*{Effects in anoxia model consistent with SD model}
\label{sec:resultsAD}

Last, we study  a model of AD \citep{ZAN11}; the AD model shares many features with the SD model but is more detailed (see Sec.\ 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 \citep{ZAN11}. 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.~\ref{fig:3}. In fact, the upper panel of  Fig.~\ref{fig:3} with the wild-type model is a reproduction of the simulations performed by ~\cite{ZAN11}.

Note that patients with migraine with aura are at greater risk for stroke \citep{KUR12}. Thus there is a rational 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 \cite{ZAN11} is derived from a model suggested by \cite{CRE09}. 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 Sec. Introduction). Investigating the underlying ion homeostasis in the three conditions of epilepsy, migraine, and stroke may yield interesting results in future investigtions 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.~\ref{fig:3} 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 \citep{HUE14}. The similarity supports the robustness of our results, as this model is an established model showing anoxia \citep{ZAN11} and seizure activity \citep{CRE09}. 

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.

 






\section*{Discussion}
\label{sec:discussion}

Our main result is that the mutant model is more susceptible to spreading depression (SD). With our computational model, we bridge the gap between the tail currents measured by \cite{DIC05} and altered cell function that constitutes the phenotype of migraine with aura. Importantly, in a computational model we can follow in all needed detail how the complex interactions of channel dynamics lead to altered cell function. A similar approach was taken, for instance, to link a genetic defect to its cellular phenotype in a cardiac arrhythmia by \cite{CLA99}.

In the discussion, we mainly highlight aspects of hypoexcitable vs.\ hyperexcitable and the concept of a threshold. 



\subsection*{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 times scales differing in three orders of magnitude (see below). Notwithstanding, the massive break-down 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 \cite{SHE01,COP07} 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.~\ref{fig:rate}). 




\subsection*{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.~\ref{fig:ap}) 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 break-down 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 larger threshold exists.

To show the actual situation in  Fig.~\ref{fig:2}, 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 \cite{ZAN11} 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 transfered across the membrane during spiking. This, in particular the intracellular ion concentration, cannot easily be measured even in an {\it in vitro} setup. In Fig. \ref{fig:2}, we see this by the much steeper slope of the reversal potential $E_{K}$ in the mutant model. 

\subsection*{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 \citep{KAG00}, 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 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 \cite{FIT55} and as discussed in a modern style by \cite{MIT13} and applied to migraine by \cite{DAH13}. 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. 

 





