\label{sec:monte-carlo-simul}
Monte-Carlo simulations are generated using an extension of the algorithm presented in \cite{keener11a}. Instead of using the Gillespie algorithm as in \cite{keener11a}, we use the next reaction method along the lines of \cite{bokes13a}. The algorithm is exact in the sense that the transition times can be approximated to any desired precision. The simulations were coded in C (using the GNU Scientific Library for random number generators) and carried out in Python, using the Scipy package. In between each jump in the number of open channels, the voltage is evolved according to the deterministic dynamics \[\label{eq:13} \frac{dv}{dt} = \frac{n}{N}f_{{\rm Na}}(v) + \frac{m}{M}f_{{\rm K}}(v) + f_{{\mathrm{leak}}}(v) + I_{{\mathrm{app}}},\] The solution provides the relationship between voltage and time, \[\label{eq:14} v(t) = \left(v(t_{0}) - \frac{c_{2}}{c_{1}}\right)e^{-c_{1}(t-t_{0})} + \frac{c_{2}}{c_{1}},\] where \[\begin{aligned} \label{eq:15} c_{1} &= \frac{n}{N}g_{{\rm Na}} + \frac{m}{M}g_{{\rm K}} + g_{{\mathrm{leak}}}, \\ c_{2} &= \frac{n}{N}g_{{\rm Na}}v_{{\rm Na}} + \frac{m}{M}g_{{\rm K}}v_{{\rm K}} + g_{{\mathrm{leak}}}v_{{\mathrm{leak}}} + I_{{\mathrm{app}}}.\end{aligned}\] To compute the next jump time, we compute four random jump times for each of the four possible transitions: \(n\to n\pm 1\) and \(m \to m \pm 1\). Each transition time is distributed according to \[\label{eq:16} W^{-}_{{\rm Na}}(t) =1 - e^{-\beta_{{\rm Na}}n (t - t_{0})}, \quad W^{+}_{{\rm Na}}(t) = 1- {\exp\left[ -\beta_{{\rm Na}}\int_{t_{0}}^{t}\Omega^{+}_{{\rm Na}}(v(\tau))d\tau \right]},\] \[W^{\pm}_{{\rm K}}(t) = 1- {\exp\left[ -\beta_{{\rm K}}\int_{t_{0}}^{t}\Omega^{\pm}_{{\rm K}}(v(\tau))d\tau \right]}.\] After integrating the voltage dependent transition rates we obtain for \((i=+, j = {\rm Na})\) and \((i=\pm, j = {\rm K})\), \[\label{eq:17} \int_{t_{0}}^{t}\Omega^{i}_{j}(v(\tau))d\tau = \frac{1}{c_{1}}\Omega^{i}_{j}(\frac{c_{2}}{c_{1}})(E_{\rm i}(z^{i}_{j}e^{-c_{1}(t-t_{0})}) - E_{\rm i}(z^{i}_{j})),\] where \[\label{eq:18} z^{+}_{{\rm Na}} = 4\gamma_{{\rm Na}}\left(v(t_{0}) - \frac{c_{2}}{c_{1}}\right), \quad z^{\pm}_{{\rm K}} = \pm \gamma_{{\rm K}}\left(v(t_{0}) - \frac{c_{2}}{c_{1}}\right),\] and \(E_{\rm i}\) is the exponential integral function defined as the Cauchy principal value integral, \[\label{eq:19} E_{\rm i}(x) = \int_{-\infty}^{x}t^{-1}e^{t}dt, \quad x\neq 0.\] Denote the jump times by \(t^{i}_{j}\), \(i = \pm\) and \(j = {\rm Na}, {\rm K}\), and let \(U\) be a uniform random variable. The jump times are given by the solution to \(W^{i}_{j}(t^{i}_{j}) = U\). There is one voltage independent jump time, \[\label{eq:21} t^{-}_{{\rm Na}} = -\frac{\log(U)}{n\beta_{{\rm Na}}}.\] Because three of the transition rates depend on voltage, and therefore time, the distributions for the jump times are not explicitly invertible. Hence, the next jump times are given implicitly by \[\frac{1}{c_{1}}\Omega^{+}_{{\rm Na}}(\frac{c_{2}}{c_{1}})(E_{\rm i}(z^{+}_{{\rm Na}}e^{-c_{1}(t^{+}_{{\rm Na}}-t_{0})}) - E_{\rm i}(z^{+}_{{\rm Na}})) = -\frac{\log(U)}{\beta_{{\rm Na}}},\] \[\label{eq:22} \frac{1}{c_{1}}\Omega^{\pm}_{{\rm K}}(\frac{c_{2}}{c_{1}})(E_{\rm i}(z^{\pm}_{{\rm K}}e^{-c_{1}(t^{\pm}_{{\rm K}}-t_{0})}) - E_{\rm i}(z^{\pm}_{{\rm K}})) =-\frac{\log(U)}{\beta_{{\rm K}}}.\] To generate the voltage dependent jump times, a root finding algorithm is applied to \eqref{eq:22} with a tolerance of \( 10^{-8}\). Once all four transition times have been computed, the next transition time is \(t^{i_{*}}_{j_{*}} = \min_{i=\pm,j={\rm Na},{\rm K}}\{t_{j}^{i}\}\). The global time is updated with \(t \leftarrow t + t^{i_{*}}_{j_{*}}\). The state is updated with \(v \leftarrow v(t^{i_{*}}_{j_{*}})\) (where \(v(t)\) is given by \eqref{eq:14} with \(t_{0}\) the time of the previous jump), \(n \leftarrow n + i_{*}\) if \(j_{*} = {\rm Na}\), and \(m \leftarrow m + i_{*}\) if \(j_{*} = {\rm K}\).