The dynamics of memory retrieval in hierarchical networks

Abstract

Memory retrieval is of central importance to a wide variety of brain functions. To understand the dynamic nature of memory retrieval and its underlying neurophysiological mechanisms, we develop a biologically plausible spiking neural circuit model, and demonstrate that free memory retrieval of sequences of events naturally arises from the model under the condition of excitation-inhibition (E/I) balance. Using the mean-field model of the spiking circuit, we gain further theoretical insights into how such memory retrieval emerges. We show that the spiking neural circuit model quantitatively reproduces several salient features of free memory retrieval, including its semantic proximity effect and log-normal distributions of inter-retrieval intervals. In addition, we demonstrate that our model can serve as a platform to examine memory retrieval deficits observed in neuropsychiatric diseases such as Parkinson’s and Alzheimer’s diseases. Furthermore, our model allows us to make novel and experimentally testable predictions, such as the prediction that there are long-range correlations in the sequences of retrieved items.

Appendices

Appendix A: Mean-field analysis methods

1.1 A.1 Derivation of the mean-field dynamics

In our model, the ensemble of neurons can be divided into several sub-populations, including the excitatory neural modules $a, b {\dots } h$ and the inhibitory population, within each of which the neurons share the same dynamic properties and, in the statistical sense, the same structure of connections with the rest neurons in the network. In the mean-field approach, it is often the collective behaviour of such neural sub-populations (for instance, the steady-state mean firing rates) instead of the behaviour of individual neurons that is of interest, and the steady-state mean firing rate of any sub-population can be obtained by studying a typical neuron in it as we will show. For this purpose, we simplify our notations in Eqs. (1)–(4) by dropping the subscript i, which gives the dynamic equation of a typical neuron in the sub-population α

$$ C\frac{dV^{\alpha}}{dt} = -g_{L}(V^{\alpha}-V_{L}) + I^{\alpha}(t) + I^{\alpha}_{ext}, $$

(15)

where the recurrent current is

$$\begin{array}{@{}rcl@{}} I^{\alpha}(t) &=& \sum\limits_{\beta}I^{\alpha\beta}(t) \end{array} $$

(16)

$$\begin{array}{@{}rcl@{}} &=&\sum\limits_{\beta}\left[-{\sum}_{j=1}^{N^{\beta}}a^{\alpha\beta}J^{\alpha\beta}s_{j}^{\beta}(t)\left( V^{\alpha}-V_{rev}^{\beta}\right)\right] . \end{array} $$

(17)

Note that we have also dropped the subscript i and superscript α in $s_{ij}^{\alpha \beta }(t)$, this is done by neglecting the conduction delay in our analysis ($d_{ij}^{\alpha \beta } \approx 0$), which is the only dependence of $s_{ij}^{\alpha \beta }(t)$ on i and α in Eq. (5). The recurrent current I ^αβ(t) can be expressed as the sum of its mean $\bar {I}^{\alpha \beta }$ and fluctuation δ I ^αβ(t). This is done by expressing the gating variable in Eq. (17) as the sum of its mean and fluctuation as $s_{j}^{\beta }(t) = \bar {s}^{\beta } + \delta {s}_{j}^{\beta }(t)$, which gives

$$\begin{array}{@{}rcl@{}} I^{\alpha\beta}(t) &=& \bar{I}^{\alpha\beta} + \delta I^{\alpha\beta}(t) \end{array} $$

(18)

$$\begin{array}{@{}rcl@{}} &=& -\bar{C}^{\alpha\beta}J^{\alpha\beta}\bar{s}^{\beta}\left( V^{\alpha}-V_{rev}^{\beta}\right)\\ &&-\sum\limits_{j=1}^{N^{\beta}}a^{\alpha\beta}J^{\alpha\beta}\delta s_{j}^{\beta}(t)\left( V^{\alpha}-V_{rev}^{\beta}\right), \end{array} $$

(19)

where in the last step we have introduced the average number of connections $\bar {C}^{\alpha \beta }=\left \langle {\sum }_{j=1}^{N^{\beta }}a^{\alpha \beta } \right \rangle =N^{\beta }\Pr (a^{\alpha \beta } = 1)$. Note that the actual number of connections to different neurons are different, following a binomial distribution (Erdős-Rényi random graph). By only considering the average number of connections $\bar {C}^{\alpha \beta }$ in the current mean-field approach, the moments of the binomial distribution higher than one are neglected. Substituting Eqs. (19) into Eq. (15) and, after some rearrangement, we have

$$ \tilde{\tau}_{m}^{\alpha}\frac{dV^{\alpha}}{dt} = - \left( V^{\alpha}- V^{\alpha}_{ss}\right) + \frac{\delta \tilde{I}^{\alpha}(t)} {\tilde{g}_{L}^{\alpha}}, $$

(20)

where the effective leaky conductance is $\tilde {g}_{L}^{\alpha } = g_{L} + {\sum }_{\beta }\bar {C}^{\alpha \beta }J^{\alpha \beta } \bar {s}^{\beta }$, the effective membrane time constant is $\tilde {\tau }_{m}^{\alpha } = C/\tilde {g}_{L}^{\alpha }$, and the steady-state membrane potential in the absence of threshold is $V^{\alpha }_{ss} = (g_{L}V_{L} + {\sum }_{\beta }\bar {C}^{\alpha \beta }J^{\alpha \beta } \bar {s}^{\beta }V^{\beta }_{rev} + I^{\alpha }_{ext}) / \tilde {g}_{L}^{\alpha }$ (Renart et al. 2003). Note that to avoid dealing with the multiplicative dependence of the fluctuation in the recurrent current δ I ^αβ(t) on the membrane potential V ^α(t), the membrane potential V ^α(t) in Eq. (19) is approximated by its mean value $\bar {V}^{\alpha }$ (Renart et al. 2003), that is,

$$\begin{array}{@{}rcl@{}} \delta I^{\alpha\beta}(t) &\approx& \delta \tilde{I}^{\alpha\beta}(t) \end{array} $$

(21)

$$\begin{array}{@{}rcl@{}} &=& - \sum\limits_{j=1}^{N^{\beta}}a^{\alpha\beta}J^{\alpha\beta}\delta s_{j}^{\beta}(t)\left( \bar{V}^{\alpha}-V_{rev}^{\beta}\right). \end{array} $$

(22)

The mean membrane potential $\bar {V}^{\alpha }$ is different from the steady-state membrane potential in the absence of threshold $V^{\alpha }_{ss}$ due to the mechanisms of resetting membrane potential to V _{r t} for an absolute refractory period τ _{r e f} immediately after firing in the integrate-and-fire neuron model. An approximate expression of $\bar {V}^{\alpha }$ is Renart et al. (2003)

$$ \bar{V}^{\alpha} = V_{ss}^{\alpha} - (V_{th} - V_{rt})\tilde{\tau}_{m}v^{\alpha} - \left( V_{ss}^{\alpha} - V_{rt}\right)\tau_{ref}v^{\alpha}, $$

(23)

where v ^α is the average firing rate of the sub-population α.

To simplify the dynamics of the gating variable described by Eq. (5), we neglect the conduction delay $\left (d_{ij}^{\alpha \beta } \approx 0\right )$ and the saturation contributed by two consecutive spikes $t^{\beta }_{j,k}$ and $t^{\beta }_{j,k+1}$, since the firing rate of interest is much lower than $\left (\tau _{d}^{\beta }\right )^{-1} \approx 200$ Hz. And we further approximate the time-course of the neurotransmitter concentration as a delta function, that is

$$ h\left( t - t^{\beta}_{j,k}\right) \approx A^{\beta} \delta\left( t-t^{\beta}_{j,k}\right), $$

(24)

where A ^β is a normalization factor to keep ${\int }_{-\infty }^{+\infty } \hat {s}_{j}^{\beta }(t) dt$ unchanged under the simplification Eq. (24) and $\hat {s}_{j}^{\beta }(t)$ is the time course of the gating variable $s_{j}^{\beta }(t)$ in response to a single spike coming at t = 0 with the initial condition $s_{j}^{\beta }(-\infty ) = 0$. This normalization factor is given by,

$$ A^{\beta} = \frac{\left( \tau_{r}^{\beta}\right)^{2} + \tau_{r}^{\beta}\tau_{d}^{\beta} + \left( \tau_{r}^{\beta}\right)^{2}}{\tau_{r}^{\beta} + \tau_{d}^{\beta}}\left( 1 - e^{-1-\tau_{r}^{\beta}/\tau_{d}^{\beta}}\right), $$

(25)

Substituting these simplifications into Eq. (5) gives

$$ \frac{ds_{j}^{\beta}}{dt} = -\frac{s_{j}^{\beta}}{\tau_{d}^{\beta}}+A^{\beta}\sum\limits_{t^{\beta}_{j}}\delta\left( t-t^{\beta}_{j}\right). $$

(26)

Now we want to cast the above dynamic equation of gating variable into a dynamic equation of the fluctuations in recurrent current $\delta \tilde {I}^{\alpha }(t)$, because it is $\delta \tilde {I}^{\alpha }(t)$, rather than the gating variable $s^{\beta }_{j}(t)$, that directly affects the membrane potential as shown in Eq. (1). For this purpose, we substitute $s_{j}^{\beta }(t) = \bar {s}^{\beta } + \delta {s}_{j}^{\beta }(t)$ into Eq. (26), multiply both sides by $\tau _{d}^{\beta }a^{\alpha \beta }J^{\alpha \beta }$ and sum over j; this yields

$$\begin{array}{@{}rcl@{}} \tau_{d}^{\beta}\frac{d}{dt} \left( \sum\limits_{j=1}^{N^{\beta}}a^{\alpha\beta}J^{\alpha\beta} \delta s^{\beta}_{j} \right)= &-&\bar{C}^{\alpha\beta}J^{\alpha\beta}\bar{s}^{\beta} - \sum\limits_{j=1}^{N^{\beta}}a^{\alpha\beta}J^{\alpha\beta}\delta s_{j}^{\beta}\\ &+&\!\tau_{d}^{\beta}J^{\alpha\beta}A^{\beta}{\sum}_{j=1}^{N^{\beta}}a^{\alpha\beta}\sum\limits_{k}\delta\!\left( t\!-t^{\beta}_{j,k}\right). \end{array} $$

(27)

To simply Eq. (27), firstly note that the number of synapses into one neuron is large ($\bar {C}^{\alpha \beta } \gg 1$) and that the individual synaptic couplings are weak, by which we mean the peak amplitude of an unitary EPSP or IPSP due to recurrent input (<1m V) is much less than what is required to excite the post-synaptic neuron from resting potential to fire (20m V). With the above two conditions, we can apply the diffusion approximation to simplify the dynamics with an assumption that the presynaptic neurons in the sub-population β fire approximately as independent Poisson processes (Stein 1967; Tuckwell 2005; Richardson and Gerstner 2005). The diffusion approximation simplifies the last term on the right hand side of Eq. (27) into Gaussian white noise,

$$\begin{array}{@{}rcl@{}} \tau_{d}^{\beta}J^{\alpha\beta}A^{\beta}\sum\limits_{j=1}^{N^{\beta}}a^{\alpha\beta}\sum\limits_{k}\delta\left( t-t^{\beta}_{j,k}\right) &=& \tau_{d}^{\beta}J^{\alpha\beta}A^{\beta}\left[\bar{C}^{\alpha\beta}v^{\beta}\right.\\ &&\,+\, \left.\sqrt{\bar{C}^{\alpha\beta}v^{\beta}}\eta^{\beta}\!(t)\!\right ]\!, \end{array} $$

(28)

where η(t)^β is zero-mean Gaussian white noise with unit variance. Substituting the diffusion approximation Eqs. (28) into Eq. (27) gives,

$$\begin{array}{@{}rcl@{}} \tau_{d}^{\beta}\frac{d}{dt} \left( \sum\limits_{j=1}^{N^{\beta}}a^{\alpha\beta}J^{\alpha\beta} \delta s^{\beta}_{j} \right) &=& \bar{C}^{\alpha\beta}J^{\alpha\beta}\left( A^{\beta}v^{\beta}\tau_{d}^{\beta}-\bar{s}^{\beta}\right)\\ &&- \sum\limits_{j=1}^{N^{\beta}}a^{\alpha\beta}J^{\alpha\beta}\delta s_{j}^{\beta}\\ &&+ \tau_{d}^{\beta}J^{\alpha\beta}A^{\beta}\sqrt{\bar{C}^{\alpha\beta}v^{\beta}}\eta^{\beta}(t), \end{array} $$

(29)

where letting the first term on the right hand side be zero gives a self-consistent solution for the stationary mean of the gating variable, $\bar {s}^{\beta } = A^{\beta }v^{\beta }\tau _{d}^{\beta }$.

Multiplying both sides of Eq. (29) by $-\left (\bar {V}^{\alpha }-V_{rev}^{\beta }\right )$ gives the dynamic equation governing the fluctuation of each individual recurrent current $\delta \tilde {I}^{\alpha \beta }$

$$ \tau_{d}^{\beta} \frac{d}{dt} (\delta \tilde{I}^{\alpha\beta}) = - \delta \tilde{I}^{\alpha\beta} + \sigma_{I}^{\alpha\beta} \eta^{\beta}(t), $$

(30)

where $\sigma _{I}^{\alpha \beta } = \tau _{d}^{\beta }J^{\alpha \beta }A^{\beta }\sqrt {\bar {C}^{\alpha \beta }v^{\beta }}|\bar {V}^{\alpha }-V_{rev}^{\beta }|$ and any extra minus sign is absorbed into the Gaussian white noise η ^β(t). Note that Eq. (30) describes an Ornstein-Uhlenbeck (O-U) process (Uhlenbeck and Ornstein 1930) for the fluctuations in each recurrent current component $\delta \tilde {I}^{\alpha \beta }(t)$. The last step necessary for obtaining the dynamic equation governing the fluctuation in total recurrent current $\delta \tilde {I}^{\alpha }(t) ={\sum }_{\beta } \delta \tilde {I}^{\alpha \beta }(t)$ is to introduce an effective synaptic time constant $\tilde {\tau }_{d}^{\alpha }$ to replace $\tau _{d}^{\beta }$ in Eq. (30) (Fourcaud and Brunel 2002) so that we can sum Eq. (30) over β as

$$\begin{array}{@{}rcl@{}} \tilde{\tau}_{d}^{\alpha} \sum\limits_{\beta} \frac{d}{dt} (\delta \tilde{I}^{\alpha\beta}) &=& - \sum\limits_{\beta} \delta \tilde{I}^{\alpha\beta} + \sum\limits_{\beta} \sigma_{I}^{\alpha\beta} \eta^{\beta}(t), \end{array} $$

(31)

$$\begin{array}{@{}rcl@{}} \tilde{\tau}_{d}^{\alpha} \frac{d}{dt} (\delta \tilde{I}^{\alpha}) &=& - \delta \tilde{I}^{\alpha} + \sigma_{I}^{\alpha} \eta^{\alpha}(t), \end{array} $$

(32)

where in the last step we have used the fact that the sum of several independent zero mean Gaussian processes is still a zero mean Gaussian process, that is, ${\sum }_{\beta } \sigma _{I}^{\alpha \beta } \eta ^{\beta }(t) = \sigma _{I}^{\alpha } \eta ^{\alpha }(t)$ with a new variance $\left (\sigma _{I}^{\alpha }\right )^{2} = {\sum }_{\beta } \left (\sigma _{I}^{\alpha \beta }\right )^{2}$. The effective synaptic time constant $\tilde {\tau }_{d}^{\alpha }$ is chosen to preserve the stationary variance of $\delta \tilde {I}^{\alpha }$ (Fourcaud and Brunel 2002). Note that the stationary solution of the O-U process described by Eq. (30) is a zero-mean Gaussian process with a variance of $\left (\sigma _{I}^{\alpha \beta }\right )^{2}/\left (2\tau _{d}^{\beta }\right )$ (Uhlenbeck and Ornstein 1930). Therefore, preserving the total stationary variance summed over the sub-populations in the effective O-U process of Eq. (32) gives the equation for $\tilde {\tau }_{d}^{\alpha }$

$$ \sum\limits_{\beta} \frac{ \left( \sigma_{I}^{\alpha\beta}\right)^{2}}{2\tau_{d}^{\beta}} = \frac{{\sum}_{\beta} \left( \sigma_{I}^{\alpha\beta}\right)^{2}}{2\tilde{\tau}_{d}^{\alpha}}. $$

(33)

The effective synaptic time constant approximation holds valid when the difference among the individual time constants and variances are within an order of magnitude (Fourcaud and Brunel 2002), which is the case in our model.

1.2 A.2 Numerical solution to the stationary mean firing rate

For a typical neuron in the sub-population α, Eq. (20) and (32) together describe a bi-variate Markov process (Gillespie 1991; Risken 1984). To obtain the mean firing rate v ^α of the sub-population α is equivalent to find the expected firing rate of this Markov process. In the case where $k = \sqrt {\tilde {\tau }_{d}^{\alpha }/\tilde {\tau }_{m}^{\alpha }} \ll 1$, an analytical solution to the expected firing rate of the Markov process can be obtained via perturbation theory to the first order of k (Fourcaud and Brunel 2002). However, when k is comparable to 1, which is quite often the case in the conductance-based model, no analytical solution has been found, although extrapolation to the second order in k and extrapolation between small and large k limits have been attempted in similar models (Brunel and Sergi 1998; Brunel and Latham 2003). Therefore, to obtain the accurate mean firing rate when k≈1, we resort to numerical solution via efficient direct simulation of Eq. (20) and (32). Firstly, note that these two equations can be readily casted into a single stochastic harmonic oscillator by differentiating Eq. (20) once and eliminating the variable $\delta \tilde {I}^{\alpha }$, which gives

$$ \tilde{\tau}_{m}^{\alpha} \frac{d^{2}}{dt^{2}} \hat{V}^{\alpha} + \left( 1+\frac{\tilde{\tau}_{m}^{\alpha}} {\tilde{\tau}_{d}^{\alpha}}\right) \frac{d}{dt} \hat{V}^{\alpha} + \frac{1}{\tilde{\tau}_{d}^{\alpha}} \hat{V}^{\alpha} = \frac{\tilde{g}^{\alpha}_{L}\tilde{\tau}_{d}^{\alpha}}{\sigma_{I}^{\alpha}} \eta^{\alpha}(t), $$

(34)

where $\hat {V}^{\alpha } = V^{\alpha } - V_{ss}^{\alpha }$. An exact simulation scheme has been developed for such harmonic oscillator driven by Gaussian white noise (Nørrelykke and Flyvbjerg 2011). After removing the numerical singularity, the simulation scheme can be applied to Eq. (34) to obtain the mean firing rates (see Appendix B). This numerical solution essentially provides a stochastic non-linear input-output relation v _{o u t} = Φ(v _{i n} ), where v _{i n} is the vector of the input firing rates to the neuron sub-populations and v _{o u t} is the resultant vector of the output firing rates. This is conceptually equivalent to the input-output relations used in firing-rate models, such as the simple sigmoid function (Wilson and Cowan 1972). To solve the stationary mean firing rates of the sub-populations in the recurrent network is to find the roots of the self-consistent equation v = v _{i n} = v _{o u t} , or, equivalently,

$$ \Phi (\boldsymbol{v}) - \boldsymbol{v} = 0, $$

(35)

where v = [v ^α]. Solving the above self-consistent equation can be casted as a stochastic optimization problem with an objective function |Φ(v)−v| to be minimized, which can be readily solved by stochastic optimization algorithms such as the pattern search algorithm (Hooke and Jeeves 1961). Starting from different initial values for the optimization algorithm, different local minimums are found and the residue errors are calculated, which is defined as |Φ(v)−v|/|v|; only the local minimums with a residue error less than 20 % are accepted as solutions to Eq. (35). The solutions are summarized in Table 2.

Approximate non-stationary firing rate dynamics of the neuron sub-populations are commonly used for stability analysis in current based integrate-and-fire neuron circuit models to bypass the difficulty of analyzing the exact non-stationary dynamics (Amit and Brunel 1997; Mascaro and Amit 1999; Renart et al. 2003). However, in our conductance based model, the time-dependency of the system parameters $\tilde {\tau }_{m}^{\alpha }$, $\tilde {\tau }_{d}^{\alpha }$ and $\tilde {g}^{\alpha }_{L}$ in Eq. (20) renders such approximate method invalid, because the approximate dynamics assume that such system parameters are constant. Nonetheless, without stability analysis, the above mean-field results for fixed points, when compared with results from simulating the full spiking neural circuit model, still provide valuable insights into the dynamics of the latter.

Qualitatively, from the perspective of attractor dynamics, the IRIs may be understood as the escape times of a “particle” in a simplified multi-well potential being driven by noise, as suggested by Litwin-Kumar and Doiron (2012). Based on Kramer’s classic reaction-rate theory (Hänggi et al. 1990), the time it takes for the particle to escape from one well into another is mainly controlled by the statistical properties of the noise and the depths of the two wells. As shown by the mean-field analysis (Eqs. (30) and (33)), the noise can be approximated as Gaussian white noise, and its intensity depends on several factors, including the synaptic time constants, coupling strengths, presynaptic firing rates, number of synapses (connection probabilities and network size) and the mean membrane potential. And the potential well depths may also depend on such a broad range of system parameters. Especially, among all the parameters, the coupling strengths and hierarchical connection probabilities play an important role in determining the IRI timescales, as demonstrated in Figs. 6 and 7, where the switching/transition frequencies are, by definition, inversely correlated with the IRI timescale. However, Kramer’s reaction-rate theory predicts an exponential tail of the escape time distribution in the low noise intensity limit (Hänggi et al. 1990), which is qualitatively different from the log-normal distributions observed in our model as well as in a recent firing-rate network model (Stern et al. 2014), and it remains an open problem to explain the generative mechanisms of such a log-normal distribution.

The observation that the transition probability decreases while the transition duration (IRI) increases as the connection probability decreases (Fig. 3) might also be explained from the perspective of attractor dynamics. As strongly indicated by the simulation results in Fig. 2 and the mean-field analytical result in Table 2, a transition from one clustered firing state to another (mean-field states of the first kind) is achieved by the system via going through one or multiple intermediate disordered states (mean-field states of the second kind). For example, in Fig. 2b, the level-3 transition happening between 302 and 305 seconds consists of obviously more intermediate disordered states than the previous transitions at lower levels. This suggests that, when two excitatory neural modules are connected with a smaller probability (i.e., at a higher hierarchy level), their corresponding clustered firing states may be further apart in the attractor landscape, and as a result the transition between them will require a longer duration to go through the intermediate attractors and have a smaller probability to happen.

Appendix B: Numerical methods for simulating stochastic harmonic oscillator

Here we describe in details the method used to simulate the stochastic process described by Eq. (34) in the main text, which is modified from the methods in Nørrelykke and Flyvbjerg (2011) to remove the numerical singularity at k = 1 in the original expressions. Firstly, we introduce the definitions of the parameters that will be used. The cyclic frequency of the undamped oscillator is $\omega ^{\alpha }_{0}=1/\sqrt {\tilde {\tau }_{d}^{\alpha }\tilde {\tau }_{m}^{\alpha }}$ and the characteristic time τ ^α of the exponential decrease with time that the momentum of the “particle” undergoes in the absence of all but friction forces is

$$ \tau^{\alpha} = \frac{\tilde{\tau}_{d}^{\alpha}\tilde{\tau}_{m}^{\alpha}}{\tilde{\tau}_{d}^{\alpha} + \tilde{\tau}_{m}^{\alpha}}. $$

(36)

The diffusion constant is defined as

$$ D^{\alpha} = \frac{\left( \sigma_{I}^{\alpha}\right)^{2}}{2\tilde{g}_{L}^{2}\left( \tilde{\tau}_{d}^{\alpha} + \tilde{\tau}_{m}^{\alpha}\right)^{2}} $$

(37)

And the cyclic frequency of the damped oscillator is $\omega ^{\alpha } = \sqrt {(\omega ^{\alpha }_{0})^{2}-1/(2\tau ^{\alpha })^{2}}$, which is real for under-damped oscillator. However, substituting Eq. (36) into the expression shows that ω ^α is zero when $\tilde {\tau }_{d}^{\alpha } = \tilde {\tau }_{m}^{\alpha }$ and always imaginary otherwise, indicating an over-damped oscillator. For mathematical convenience, we keep the imaginary frequency ω ^α. Note that the difficulty in obtaining analytical solution to the mean first passage time of the bi-variate Markov process described by Eqs. (20) and (32) in the main text is consistent with the difficulty in solving it for an over-damped stochastic harmonic oscillator (Verechtchaguina et al. 2006).

With the above definitions, together with suitable time step Δt, t _j = jΔt, $\hat {V}^{\alpha }_{j} = \hat {V}^{\alpha }(t_{j})$, $d\hat {V}^{\alpha }/dt = Y^{\alpha }$, and $Y^{\alpha }_{j} = Y^{\alpha }(t_{j})$, the recursive relation for simulation is given by

$$\begin{array}{@{}rcl@{}} \left( \begin{array}{l} {\hat{V}}_{j+1}^{\alpha}\\ {Y}_{j+1}^{\alpha} \end{array}\right) &=& e^{-\mathbf{M}^{\alpha}\Delta t} \left( \begin{array}{l} {\hat{V}}_{j}^{\alpha}\\ {Y}_{j}^{\alpha} \end{array}\right)\\ &&+ \left( \begin{array}{cc} {\sigma}_{vv}^{\alpha} & 0\\ \left( {\sigma}_{vy}^{\alpha}\right)^{2} / {\sigma}_{vv}^{\alpha} & \sqrt{ \left( {\sigma}_{yy}^{\alpha}\right)^{2} - \left( {\sigma}_{vy}^{\alpha}\right)^{4}/ \left( {\sigma}_{vv}^{\alpha}\right)^{2}} \end{array}\right)\\ &&\times \left( \begin{array}{l} {\xi}_{j}^{\alpha} \\ {\zeta}_{j}^{\alpha}, \end{array}\right) \end{array} $$

(38)

where

$$\begin{array}{@{}rcl@{}} e^{-\mathbf{M}^{\alpha}\Delta t} &=& e^{-\frac{\Delta t}{2\tau^{\alpha}}} \cos(\omega^{\alpha} \Delta t) \left( \begin{array}{ll} 1 & 0 \\ 0 & 1 \end{array}\right)\\ &&+ e^{-\frac{\Delta t}{2\tau^{\alpha}}} \text{sinc}(\omega^{\alpha} \Delta t)\Delta t\! \left( \begin{array}{cc} \frac{1}{2\tau^{\alpha}} & 1 \\ -\left( {\omega}_{0}^{\alpha}\right)^{2} & -\frac{1}{2\tau^{\alpha}} \end{array}\right), \end{array} $$

(39)

$$\begin{array}{@{}rcl@{}} \left( {\sigma}_{vv}^{\alpha}\right)^{2} &=& \frac{D^{\alpha}}{\left( {\omega}_{0}^{\alpha}\right)^{2}\tau^{\alpha}} - \frac{D^{\alpha} e^{-\Delta t/ \tau^{\alpha}}}{2\left( \omega^{\alpha}_{0}\right)^{2}(\tau^{\alpha})^{3}}\\ &&\times \left\{ \left[\Delta t \text{sinc}(\omega^{\alpha} \Delta t) + \tau^{\alpha}\right]^{2} + (\tau^{\alpha})^{2} \right\}, \end{array} $$

(40)

$$\begin{array}{@{}rcl@{}} \left( {\sigma}_{yy}^{\alpha}\right)^{2} &=& \frac{D^{\alpha}}{\tau^{\alpha}} - \frac{D^{\alpha} e^{-\Delta t/ \tau^{\alpha}}}{2(\tau^{\alpha})^{3}}\\ &&\times \left\{ \left[{\Delta} t \text{sinc}(\omega^{\alpha} {\Delta} t) - \tau^{\alpha}\right]^{2} + (\tau^{\alpha})^{2} \right\}, \end{array} $$

(41)

$$ \left( {\sigma}_{vy}^{\alpha}\right)^{2} = \frac{D^{\alpha}\Delta t^{2}}{(\tau^{\alpha})^{2}} e^{-\Delta t/ \tau^{\alpha}} \text{sinc}^{2}(\omega^{\alpha} \Delta t), $$

(42)

where ${\xi }_{j}^{\alpha }$ and ${\zeta }_{j}^{\alpha }$ are independent zero-mean Gaussian white noise with unit variance, and $\text {sinc}(x) = \sin (x) / x$ is known as the unnormalized cardinal sine function. Note that the above simulation scheme is exact for time steps of arbitrary size (Nørrelykke and Flyvbjerg 2011). In practice, a time step of Δt = 1 ms is chosen for a balance between performance and temporal resolution. Analytical solution to special cases when k≪1 (Fourcaud and Brunel 2002) is used for code testing.

With the above simulation scheme, the mean first passage time of the dynamic system described by Eq. (34) in the main text can be obtained with the following procedure. Firstly the initial value of the membrane potential is set to be at the reset potential, or, equivalently, ${\hat {V}}_{0}^{\alpha } = V_{rt} - {V}_{ss}^{\alpha }$. And the initial value ${Y}_{0}^{\alpha }$ is self-consistently chosen from the stationary solution of ${Y}_{t}^{\alpha }$, which is a zero-mean Gaussian distribution with the variance given by Chandrasekhar (1943)

$$ \sigma^{2} \left( {Y}_{\infty}^{\alpha}\right) = \left( \frac{1}{{\tilde{\tau}}_{d}^{\alpha}} + \frac{1}{{\tilde{\tau}}_{m}^{\alpha}}\right) D^{\alpha}. $$

(43)

Then, for each trial, the time ΔT ^α it takes for the system to reach the threshold ${\hat {V}}_{th}^{\alpha } = V_{th} - {V}_{ss}^{\alpha }$ for the first time is recorded (first passage time). Finally, the sum of the refractory period τ _{r e f} and the mean first passage time averaged over trials 〈ΔT ^α〉 gives an estimate on the inverse of the mean firing rate, that is

$$ v^{\alpha} = (\tau_{ref} + \langle \Delta T^{\alpha} \rangle)^{-1}. $$

(44)

In practice, 5000 trails are used for each estimation.

Appendix C: Bayesian model comparison

Here we describe the detailed procedures of the Bayesian comparison between the ex-Gaussian model and the log-normal model of the IRI distributions. Firstly, the distribution of IRIs data is extracted from Fig. 1 in study by Rohrer and Wixted (1994) and the data is in the form of a histogram. We denote the bin width of the histogram as ΔB and the i ^th bin centre as B _i = B ₁+(i−1)ΔB, where $i = 1, 2 {\dots } M$ and M is the total number of bins. And we denote the number of entries (IRIs) in the i ^th bin as d _i . In the data, ΔB = 1 second, B ₁ = 0.5 second, M = 20 and the sample size $N = {\sum }_{i = 1}^{M} d_{i} = 1186$.

The two models for the distribution of IRIs to be compared are the log-normal distribution and the ex-Gaussian distribution. The probability density function (PDF) for a log-normal distribution is

$$ f_{1}(x; \mu_{1}, \sigma_{1}) = \frac{1}{x\sigma_{1} \sqrt{2\pi}} e^{-\frac{(\ln x - \mu_{1})^{2}}{2{\sigma_{1}^{2}}}}, $$

(45)

where μ ₁ and σ ₁ are the distribution parameters. The PDF for a ex-Gaussian distribution is

$$\begin{array}{@{}rcl@{}} f_{2}(x; \lambda_{2}, \mu_{2}, \sigma_{2}) &=& \frac{\lambda_{2}}{2} e^{\frac{\lambda_{2}}{2} \left( 2\mu_{2} + \lambda_{2} {\sigma_{2}^{2}} - 2x\right)} \\ &&\cdot H\left( \frac{\mu_{2} + \lambda_{2} {\sigma_{2}^{2}} - x}{\sqrt{2} \sigma_{2}}\right), \end{array} $$

(46)

where λ ₂, μ ₂ and σ ₂ are the distribution parameters and $H(x) = 2/\sqrt {\pi } {\int }^{+\infty }_{x} e^{- t^{2}} dt$, which is the complementary error function. For convenience, we refer to the two PDF models as f _j (x;𝜃 _j ) or f _j (𝜃 _j ), where j = 1,2 is the model index and 𝜃 ₁ = (μ ₁, σ ₁) and 𝜃 ₂ = (λ ₂, μ ₂, σ ₂) are the model parameter vectors.

Given the sample size N, i.e. the total number of entries in all of the bins, the expected number of entries t _{i, j} in the i ^th bin given by the j ^th model is

$$ t_{i,j} = p_{i,j} N, $$

(47)

where p _{i, j} is the probability for an observed IRI to be within the i ^th bin given by the j ^th model, which can be obtained by

$$ p_{i,j} = \frac{ {\int}^{B_{i} + {\Delta} B /2}_{B_{i} - {\Delta} B /2} f_{j}(x; \boldsymbol{\theta}_{j} ) dx}{{\int}^{B_{M} + {\Delta} B /2}_{B_{1} - {\Delta} B /2} f_{j}(x; \boldsymbol{\theta}_{j} ) dx}. $$

(48)

Note that, in the above expression, the term on the denominator renormalizes the PDF over the range of bins, that is, the distribution is truncated after the last bin. With enough number of bins and a sample size large enough, the observed number of entries in the i ^th bin should follow a Poisson distribution with the parameter t _{i, j} . Therefore, the likelihood for t _{i, j} expected and d _i actually observed entry in the i ^th bin given by the j ^th model is $L_{i,j}(\boldsymbol {\theta }_{j}) = e^{-t_{i,j}}t_{i,j}^{d_{i}}/d_{i} !$. Summing this over all the bins, taking the logarithm and discarding model independent terms give the log-likelihood of the observed histogram under the j ^th model with parameter vector 𝜃 _j

$$ \ln L_{j}(\boldsymbol{\theta}_{j}) = \sum\limits_{i=1}^{M} (-t_{i,j} + d_{i} \ln t_{i,j} ). $$

(49)

With this expression, we can readily use the Schwarz criterion (Schwarz et al. 1978) to approximate the Bayes factor B ₁₂, that is

$$ \ln B_{12} \approx \ln L_{1}(\hat{\boldsymbol{\theta}}_{1}) - \ln L_{2}(\hat{\boldsymbol{\theta}}_{2}) - \frac{1}{2}(k_{1} - k_{2}) \ln N, $$

(50)

where $\hat {\boldsymbol {\theta }}_{1}$ and $\hat {\boldsymbol {\theta }}_{2}$ are the maximum likelihood estimators of the model parameter vectors that maximize the log-likelihood of the observed histogram as in Eq. (49), N is the sample size and k ₁ and k ₂ are the dimensions of the model parameter vectors 𝜃 ₁ and 𝜃 ₂, respectively. This approximation should provide a reasonable indication for the model comparison in large samples (N≫k ₁, k ₂) (Kass and Raftery 1995), and in our case N is around 10³ while k ₁ = 2, k ₂ = 3.