Mapping input noise to escape noise in integrate-and-fire neurons: a level-crossing approach

Abstract

Noise in spiking neurons is commonly modeled by a noisy input current or by generating output spikes stochastically with a voltage-dependent hazard rate (“escape noise”). While input noise lends itself to modeling biophysical noise processes, the phenomenological escape noise is mathematically more tractable. Using the level-crossing theory for differentiable Gaussian processes, we derive an approximate mapping between colored input noise and escape noise in leaky integrate-and-fire neurons. This mapping requires the first-passage-time (FPT) density of an overdamped Brownian particle driven by colored noise with respect to an arbitrarily moving boundary. Starting from the Wiener–Rice series for the FPT density, we apply the second-order decoupling approximation of Stratonovich to the case of moving boundaries and derive a simplified hazard-rate representation that is local in time and numerically efficient. This simplification requires the calculation of the non-stationary auto-correlation function of the level-crossing process: For exponentially correlated input noise (Ornstein–Uhlenbeck process), we obtain an exact formula for the zero-lag auto-correlation as a function of noise parameters, mean membrane potential and its speed, as well as an exponential approximation of the full auto-correlation function. The theory well predicts the FPT and interspike interval densities as well as the population activities obtained from simulations with colored input noise and time-dependent stimulus or boundary. The agreement with simulations is strongly enhanced across the sub- and suprathreshold firing regime compared to a first-order decoupling approximation that neglects correlations between level crossings. The second-order approximation also improves upon a previously proposed theory in the subthreshold regime. Depending on a simplicity-accuracy trade-off, all considered approximations represent useful mappings from colored input noise to escape noise, enabling progress in the theory of neuronal population dynamics.

1 Introduction

Neurons in the brain must operate under highly non-stationary conditions. In fact, most behaviorally relevant sensory stimuli as well as internal signals are rarely constant in time but may change rapidly. In the presence of noise, such dynamic stimuli can be reliably encoded in the time-dependent population activity of a large population of spiking neurons (Gerstner et al. 2014). The time-dependent population activity also provides a concise coarse-grained description of the collective dynamics of interacting spiking neurons. Therefore, theories that predict the population activity in response to a time-dependent signal have been of fundamental interest in theoretical neuroscience (Knight 1972; Gerstner 2000; Augustin et al. 2017; Schwalger et al. 2017).

The population activity of noisy spiking neurons can be mathematically described by population density equations (Nykamp and Tranchina 2000; Chizhov 2017). The form of the population density equation depends on the noise model. Two popular ways to model neuronal noise consist of modeling noise either in the input or in the output of the neuron (Gerstner et al. 2014). In the first model class (input noise), noise enters the dynamical equations of the membrane potential, currents or conductances leading to stochastic differential equations. If the noise is Gaussian white noise, the subthreshold dynamics becomes a diffusion process and the input noise is also called diffusive noise (Gerstner 2000). The corresponding population density equation is a Fokker–Planck equation and the population activity can be obtained as the probability flux across the threshold (Abbott and van Vreeswijk 1993; Brunel and Hakim 1999; Fourcaud and Brunel 2002; Nykamp and Tranchina 2000; Richardson 2008; Augustin et al. 2017). Models based on diffusive noise naturally appear as the result of modeling biophysical processes such as synaptic shot-noise or ion channel noise. In particular, a frequently considered source of noise is background synaptic input modeled as external Poisson processes (Brunel 2000; Potjans and Diesmann 2014). The fluctuating part of this external shot noise leads, via a diffusion approximation (Gerstner et al. 2014), to Gaussian white noise driving the synaptic input current or conductance. Besides its biophysical interpretability, input noise has the advantage that it permits modeling both temporal (Fourcaud and Brunel 2002; Schwalger et al. 2015) and spatial (Lindner et al. 2005) correlations of synaptic inputs and enables mean-field theories for recurrent networks of sparsely connected integrate-and-fire neurons (Brunel and Hakim 1999; Brunel 2000).

In the second model class (called output noise or escape noise (Gerstner 2000)), the dynamical equations for the state variables are deterministic, while spikes (“output”) are generated stochastically through a hazard rate or conditional intensity (Gerstner 2000; Paninski 2004; Truccolo et al. 2005; Pillow et al. 2008; Pillow and Latham 2008; Naud and Gerstner 2012; Brea et al. 2013; Galves and Löcherbach 2016; Gerhard et al. 2017; Raad et al. 2020). This hazard rate depends on the state variables via a link function. For example, it may be given as ${{\hat{\lambda }}}(t)=\varPsi (u(t),{\hat{t}}(t))$, where u(t) is the membrane potential and ${\hat{t}}(t)$ is the last spike time of the neuron at time t. If the neuron model is a non-homogeneous renewal or quasi-renewal (Naud and Gerstner 2012) process, the corresponding population density equation is a renewal integral equation or, equivalently, a refractory density equation (Gerstner et al. 2014; Gerstner 2000; Naud and Gerstner 2012; Chizhov and Graham 2007, 2008; Dumont et al. 2016; Schwalger and Chizhov 2019). Although output noise is of phenomenological nature without a quantitative link to biophysical mechanisms, it has several advantages (Schwalger and Chizhov 2019) owing to its simpler mathematical tractability: First, the refractory density or integral equation admits an extension to finite numbers of neurons (Schwalger et al. 2017; Schwalger and Chizhov 2019; Schmutz et al. 2020, 2021). This extension allows to account for finite-size fluctuations of the population activity at the mesoscopic scale. Second, models with output noise provide analytical expressions for the likelihood function, and thus, model parameters can be efficiently fitted to experimental data of single neuron recordings (Paninski 2004; Truccolo et al. 2005; Pillow et al. 2008; Gerhard et al. 2017; Mensi et al. 2012; Pozzorini et al. 2015; Teeter et al. 2018). And third, the state space for models with output noise remains approximately one-dimensional even for multi-dimensional conductance-based neuron models (Chizhov and Graham 2007). The one-dimensional description permits highly efficient numerical solutions, in contrast to Fokker–Planck equations (Apfaltrer et al. 2006), which become intractable and computationally inefficient for several state variables.

In view of the wide use of biologically interpretable input noise and the mathematical advantages of output noise, an intriguing question is whether input noise can be approximately mapped to output noise, so as to take full advantage of both noise models. Mathematically, such a map requires the specification of the hazard rate ${{\hat{\lambda }}}(t)$ in terms of a link function $\varPsi$, which depends on some dynamical variables and defines the escape-noise model. Unfortunately, a standard method to derive such a link function does not exist. To see this, let us consider the example of non-homogeneous renewal processes as a popular class of neuron models. In these models, the probability density $P(t|{\hat{t}})$ to fire the next spike at time t given a spike at time ${\hat{t}}$, ${\hat{t}}<t$, does not depend on the state of the model before time ${\hat{t}}$, i.e., the memory of renewal neurons only reaches back to its last spike. An important example of non-homogeneous renewal models in neuroscience are one-dimensional integrate-and-fire neurons driven by white input noise (Gerstner et al. 2014). For this model class one can formally construct the hazard rate via the formula $\lambda (t|{\hat{t}})=P(t|{\hat{t}})/\bigl [1-\int _{{\hat{t}}}^tP(s|{\hat{t}})\,\hbox {d}s\bigr ]$ (Gerstner et al. 2014). However, there are two obstacles: first, in order to apply this formula, the “interspike interval (ISI) density” $P(t|{\hat{t}})$ would be needed in analytical form for arbitrary, time-dependent input currents $\{I(t')\}_{t'\in ({\hat{t}},t)}$ that occurred since the last spike. However, the calculation of the ISI density for time-dependent inputs is equivalent to a first-passage-time (FPT) problem with time-dependent parameters or boundary. The solution of this FPT problem requires the solution of the Fokker–Planck equation with moving absorbing boundary, which is known to be a hard theoretical problem (Bulsara et al. 1996; Schindler et al. 2004; Lindner 2004b). Second, even if one succeeds to derive an approximate formula for the hazard rate $\lambda (t|{\hat{t}})$, it is still challenging to represent the hazard rate in the form of a link function $\varPsi \bigl (u(t),\{z(t)\},{\hat{t}}\bigr )$ that depends on some voltage-like variable u(t), the last spike time ${\hat{t}}$ and possibly further dynamical variables $\{z(t)\}$ locally in time (as opposed to a “non-local” functional of $\{u(t'),z(t')\}_{t'\in ({\hat{t}},t)}$).

Several theoretical studies have suggested approximate local hazard rates for leaky integrate-and fire (LIF) models driven by white (Plesser and Gerstner 2000; Herrmann and Gerstner 2001; Chizhov and Graham 2007) or exponentially correlated (Chizhov and Graham 2008) Gaussian noise, or quasi-static (frozen) noise (Goedeke and Diesmann 2008). In this paper, we explore an alternative approach to the hazard rate and the first-passage-time density based on the theory of level crossings (Verechtchaguina et al. 2006). In Sect. 2, we introduce the LIF model with time-dependent driving and constant threshold and map this process an equivalent model with constant input and moving barrier. In Sect. 3, we consider the level crossing statistics with respect to this moving barrier and use the Wiener–Rice series and approximations thereof to provide formal expressions for the FPT density. These expressions form the starting point for deriving approximate hazard rates that are local in time. This derivation reveals some unexpected results concerning the correlations of level-crossings of Gaussian processes at small time lags (Sect. 3.4). Then, we turn to the LIF model and the problem of mapping input noise to escape noise (Sect. 4) and apply this map to predict the time-dependent population activity of LIF neurons with colored input noise (Sect. 5). Each of the Sects. 3,  4 and 5 closes with a comparison of the level-crossing theory with simulations and a previous theory by Chizhov and Graham (2008). Detailed derivations are provided in Appendix.

2 Leaky integrate-and-fire models and the associated first-passage-time problem

As a spiking neuron model with input noise, we consider a leaky integrate-and-fire model driven by synaptically filtered (“colored”) noise (Schwalger and Schimansky-Geier 2008; Gerstner et al. 2014; Schuecker et al. 2015). In this model, spikes are emitted whenever the membrane potential V(t) reaches a threshold $V_\text {T}$. The subthreshold dynamics for $V<V_\text {T}$ can be written as

$$\begin{aligned} \tau _\text {m}\dot{V}&=-V+\mu (t)+\eta (t), \end{aligned}$$

(1a)

$$\begin{aligned} \tau _\text {s}{{\dot{\eta }}}&=-\eta +\sqrt{2\tau _\text {s}}\sigma _\eta \xi (t), \end{aligned}$$

(1b)

where $\tau _\text {m}$ is the membrane time constant and $\mu (t)=V_\text {rest}+RI(t)$ is the mean neuronal drive consisting of a constant resting potential $V_\text {rest}$ and a time-dependent input current I(t) (R denotes the membrane resistance). Furthermore, $\eta (t)$ is a colored noise modeled as a one-dimensional Ornstein–Uhlenbeck process with correlation time $\tau _\text {s}$ and variance $\sigma _\eta ^2$, and $\xi (t)$ is a zero-mean Gaussian white noise with auto-correlation function $\langle \xi (t)\xi (t')\rangle =\delta (t-t')$. The colored noise captures the effect of various intrinsic and extrinsic noise sources, such as fluctuations of synaptic background activity in vivo (shot noise due to random spike arrival from background neurons). After threshold crossing and spike emission, V(t) is reset to a reset potential $V_R$, $V_R<V_T$, and the subthreshold dynamics Eq. (1) resumes after an absolute refractory period of length $t_\text {ref}$ following the reset.

We are seeking a corresponding spiking neuron model with escape-noise (Gerstner 2000) given by a hazard rate (conditional intensity) of the form $\varPsi \bigl (u(t),\dot{u}(t),\{z_i(t)\},t-{\hat{t}}\bigr )$. Here, u(t) is a membrane-potential variable that obeys the noiseless membrane dynamics of the LIF model between spikes:

$$\begin{aligned} \tau _\text {m}\dot{u}=-u+\mu (t). \end{aligned}$$

(2a)

Furthermore, we allow an explicit dependence on the speed of the membrane potential $\dot{u}(t)$ (in accordance with previous studies (Plesser and Gerstner 2000; Herrmann and Gerstner 2001; Chizhov and Graham 2007; Goedeke and Diesmann 2008)), the time since the last spike $t-{\hat{t}}$, and possibly further auxiliary variables $\{z_i\}$ whose dynamics between spikes is given by ordinary differential equations. Given these variables at time t, a spike is fired independently in the next time step with probability

$$\begin{aligned}&\text {Pr}\left( \text {spike in }(t,t+\hbox {d}t)|u(t),\dot{u}(t),\{z_i(t)\},t-{\hat{t}}\right) \nonumber \\&\quad =\varPsi \bigl (u(t),\dot{u}(t),{\{z_i(t)\}},t-{\hat{t}}\bigr )\hbox {d}t \end{aligned}$$

(2b)

where dt is a small step size. This probabilistic firing rule is the counterpart of the firing rule with a hard threshold in the LIF model with input noise. After a spike, u(t) is reset to $V_\text {R}$ and the auxiliary variables $\{z_i\}$ are also reset to some suitable fixed reset value. During an absolute refractory period of length $t_\text {ref}$, the variables are clamped to their reset values and the hazard rate is set to zero. Because all memory is erased upon resetting, the escape-noise model is a non-homogeneous renewal process .

The main goal is to map the model with colored input noise, Eq. (1) to the model with escape noise, Eq. (2). Strictly speaking, mapping the two models is an ill-posed problem because the model with input noise is a non-renewal process, whereas the escape-noise model is a (non-homogeneous) renewal process. In fact, the temporal correlations of the colored noise in Eq. (1) introduces memory that is not erased upon spiking. This memory leads to correlations between interspike intervals (ISIs) (Lindner 2004a; Schwalger and Schimansky-Geier 2008; Schwalger et al. 2015). However, if the correlation time $\tau _\text {s}$ of the colored noise is much smaller than the mean interspike interval, these correlations will be small and the model with input noise can be regarded as approximately renewal. In this case, it is sufficient to match the ISI densities of the two models in order to obtain an approximate mapping. Therefore, our goal of mapping the two models can be phrased more modestly as follows: Can we find a link function $\varPsi$ of the escape-noise model such that for an arbitrary given stimulus $\mu (t)$ the time-dependent ISI densities $P(t|{\hat{t}})$ of the two models approximately match for all t and ${\hat{t}}<t$? We emphasize that this definition of the mapping rests on the assumption of sufficiently small correlation times of the colored input noise. Biologically, this assumption seems to be reasonable given that typical time scales of excitatory and inhibitory postsynaptic currents are often only on the order of a few milliseconds (Gerstner et al. 2014).

To derive the link function $\varPsi$ that maps input to output noise, one needs to solve a first-passage-time (FPT) problem: As mentioned in Introduction, the hazard rate can be obtained from the ISI density of the model with input noise, Eq. (1). In this model, the interspike interval is determined by the “first-passage time” that is needed for the membrane potential to travel from the reset potential to the threshold. Thus, the ISI density $P(t|{\hat{t}})$ is equivalent to the FPT density (apart from a time shift due to the deterministic absolute refractory period). To compute the FPT density, one needs to choose suitable initial conditions for the colored noise $\eta (t)$. The ISI starting at the last spike time ${\hat{t}}$ is composed of the initial absolute refractory period of length $t_\text {ref}$ and the stochastic FPT $t^*$. We thus need the initial value $\eta ({\hat{t}}+t_\text {ref})$ of the noise at the starting time ${\hat{t}}+t_\text {ref}$ of the stochastic motion. At the firing time ${\hat{t}}$, the distribution of the noise $p_\text {fire}(\eta ,{\hat{t}})$ is biased toward positive values of $\eta$ (Lindner 2004a; Schwalger and Schimansky-Geier 2008; Schwalger 2013; Schwalger et al. 2015), in contrast to the stationary distribution $p_\text {st}(\eta )$ of the Ornstein–Uhlenbeck noise, which has zero mean. During the absolute refractory period, the noise distribution relaxes toward the stationary distribution. Even though the noise at time ${\hat{t}}+t_\text {ref}$ may not be fully stationary yet, it is reasonable to assume stationary initial conditions, where $\eta ({\hat{t}}+t_\text {ref})\sim {\mathcal {N}}(0,\sigma _\eta ^2)$ is drawn from a normal distribution with variance $\sigma _\eta ^2$. This initial condition is justified because the noise correlation time $\tau _\text {s}$ has been assumed to be much smaller than the mean ISI; hence, we do not expect that the precise shape of the initial noise distribution has a significant effect on the FPT density.

Because in the following we focus on the FPT starting at ${\hat{t}}+t_\text {ref}$, we will conveniently choose the time origin such that ${\hat{t}}+t_\text {ref}=0$. Furthermore, since we are only interested in the first threshold crossing after time $t=0$, we can omit the voltage resetting for $t>0$ without changing the FPT statistics. The resulting non-resetting process ${{\hat{V}}}(t)$ is the freely evolving solution of Eq. (1) without reset and with initial conditions ${{{\hat{V}}}}(0)=V_\text {R}$, $\eta (0)\sim {\mathcal {N}}(0,\sigma _\eta ^2)$ (Fig. 1a). This non-resetting process will be useful for the level-crossing approach below.

Fig. 1

First-passage time of an integrate-and-fire neuron model and an equivalent model with moving boundary. a At time $t=0$, different realizations of the non-resetting membrane potential ${{\hat{V}}}(t)$ (colored thin lines) are released from the reset potential $V_\mathrm{R}$. The non-resetting membrane potential follows a Gaussian process with time-dependent mean $\langle {{\hat{V}}}(t)\rangle$ (gray thick line). Shown are three realizations (green, red, blue lines) that have an identical threshold crossing at time $t=t^*$ (blue circle), which is not necessarily the first crossing (indicated by an arrow). b Transformation to an equivalent time-homogeneous process x(t) with moving boundary b(t), in which the positions of threshold crossings are preserved. Parameters: $\tau _\text {s}=4$ ms, $\tau _\text {m}=10$ ms, $\sigma _V{{:}{=}}\sigma _x(\infty )=0.25(V_\text {T}-V_\text {R})$

For mathematical convenience, we will now reformulate the FPT problem in terms of a time-homogeneous process x(t) and a moving boundary b(t), so as to eliminate the time-dependent parameter $\mu (t)$ in Eq. (1) (Fig. 1b). This is achieved by subtracting the mean non-resetting membrane potential $\langle {{{\hat{V}}}}(t)\rangle =u(t)$:

$$\begin{aligned} x(t)&={{{\hat{V}}}}(t)-u(t) \end{aligned}$$

(3)

$$\begin{aligned} b(t)&=V_T-u(t), \end{aligned}$$

(4)

where u(t) is given by Eq. (2a) with initial condition $u(0)=V_\text {R}$. Furthermore, setting $y=\eta /\tau _\text {m}$, $\gamma =1/\tau _\text {m}$, $\tau _y=\tau _\text {s}$ and $D=\tau _\text {s}\sigma _\eta ^2/\tau _\text {m}^2$, we find the Langevin equation

$$\begin{aligned} \dot{x}&=-\gamma x+y \end{aligned}$$

(5a)

$$\begin{aligned} \tau _y\dot{y}&=-y+\sqrt{2D}\xi (t) \end{aligned}$$

(5b)

with initial conditions

$$\begin{aligned} x(0)=0,\qquad y(0)\sim {\mathcal {N}}(0,\sigma _y^2) \end{aligned}$$

(6)

The dynamics of x(t) can be interpreted as an overdamped motion of Brownian particle in a parabolic potential subject to a colored noise y(t) (Ornstein–Uhlenbeck process). Here, D and $\tau _y$ are intensity and the correlation time of the noise, respectively, and $\gamma$ is the friction coefficient. As before, $\xi (t)$ is a zero-mean Gaussian white noise. At time $t=0$, the random initial condition for the colored noise y corresponds to a stationary Gaussian distribution with mean zero and variance $\sigma _y^2=D/\tau _y$. By construction, the domain of the particle is bounded from above by the time-dependent boundary b(t), where $b(0)>0$ and b(t) is a differentiable function of time. The FPT $t^*$ is defined as the time when x(t) exits the domain, i.e., when it reaches the boundary, for the first time. The FPT density will be denoted by P(t), i.e., $P(t)\hbox {d}t=\text {Prob}(t^*\in [t,t+\hbox {d}t))$ for an infinitesimal time interval of length dt. We emphasize again that the FPT density of the Brownian particle x(t) with moving boundary b(t) is the same as the FPT density of the membrane potential V(t) with respect to the constant threshold $V_\text {T}$.

Beyond neuroscience, the escape of the doubly low-pass filtered process, Eq. (5), from a domain with moving boundary b(t) may serve as a simple archetypal model for non-stationary FPT problems. One prominent example is reaction times of bimolecular chemical reactions (Hänggi et al. 1990). If x(t) is interpreted as a reaction coordinate and the domain $x<b(t)$ corresponds to the reactant state, the boundary b(t) can be interpreted as a time-dependent energy barrier that needs to be surpassed to reach the product state. Accordingly, the first-passage time can be interpreted as the reaction time.

3 Level-crossing theory for a moving barrier

3.1 Hazard-rate representation of first-passage-time density

To find approximations to the FPT density from approximate hazard rates, we use concepts from renewal theory, especially the notion of hazard rate and survival probability (Cox 1962). Because the process Eq. (5) starts at time 0, the hazard rate $\lambda (t)$ is defined here as the conditional probability per small time interval $\hbox {d}t$ to find a boundary crossing in the interval $(t,t+\hbox {d}t)$ given the absence of crossings in the interval (0, t). On the other hand, the survival probability S(t) is defined as the probability of an absence of crossings in (0, t). The two definitions imply that $S(t+\hbox {d}t)=S(t)(1-\lambda (t)\hbox {d}t)$, hence $\hbox {d}S(t)/\hbox {d}t=-\lambda (t)S(t)$. Because the survival probability is unity at time $t=0$, we thus obtain $S(t)=\exp \left( -\int _{0}^t\lambda (s)\,\hbox {d}s \right)$ for $t>0$. The probability to find the first crossing after time 0 in the interval $(t,t+\hbox {d}t)$ is equal to the probability to find a crossing in $(t,t+\hbox {d}t)$ and to have no crossings in (0, t). Hence, the FPT density is given by the product $P(t)=\lambda (t)S(t)$, or

$$\begin{aligned} P(t)=\lambda (t)\exp \left( -\int _{0}^t\lambda (s)\,\hbox {d}s \right) . \end{aligned}$$

(7)

Given the hazard rate $\lambda (t)$ for $t>0$, Eq. (7) provides a simple formula for the FPT density. An advantage of this representation is that the exponential factor can be turned into a first-order differential equation,

$$\begin{aligned} P(t)=\lambda (t)S(t),\qquad \frac{\mathrm {d}S}{\mathrm {d}t}=-\lambda (t)S(t),\qquad S(0)=1. \end{aligned}$$

(8)

Thus, if the hazard rate $\lambda (t)$ can be efficiently computed for $t>0$, this representation permits an efficient numerical integration of the first-passage-time density forward in time. Therefore, the main strategy in this paper is to derive computationally efficient approximations for the hazard rate.

In general, the calculation of the hazard rate is as difficult as the calculation of the FPT density itself. However, finding approximations for $\lambda (t)$ has several advantages over direct approximations of P(t). Firstly, as a probability density, P(t) must satisfy the normalization to unity. Thus, the value of the FPT density at different times cannot be calculated independently. In particular, the value of P(t) strongly depends on the values for $t'\in (0,t)$. By contrast, $\lambda (t)$ is not a probability density and can thus, in principle, be arbitrary as long as it is nonnegative and $S(t)=\exp \left( -\int _{0}^t\lambda (s)\,\hbox {d}s \right)$ converges to zero as $t\rightarrow \infty$. Thus, if we are able to find any approximation for $\lambda (t)$, the normalization of P(t) is guaranteed by Eq. (7).

Secondly, the character of the hazard rate is more local in time than the FPT density, and thus, we expect more efficient approximations for the hazard rate. The non-local character of P(t) has been already mentioned above. Moreover, the non-locality becomes particularly evident by the integral in Eq. (7), which accumulates the history of hazard rates. The exponential factor S(t) shaped by this integral thus contributes a trivial history-dependence of the FPT density P(t), which is present already for time-homogeneous processes. By contrast, this trivial history-dependence is divided out in the hazard rate $\lambda (t)=P(t)/S(t)$. The remaining time-dependence of the hazard rate singles out effects of non-stationarity and explicit time-dependence of the system, which can be captured by local variables. Thirdly, because of the locality in time, time-dependent rates are interesting in its own right as they are often the natural choice to model escape processes in terms of a Markovian dynamics and master equations.

From the above considerations, it becomes clear that the hazard rate representation, Eq. (7), is only useful if we succeed to derive approximations for $\lambda (t)$ that are local in time. This means that we are seeking an approximation of the hazard rate in the form

$$\begin{aligned} \lambda (t)\approx \varPhi \bigl (b(t),\dot{b}(t),\dotsc , \{z_i(t)\},t\bigr ), \end{aligned}$$

(9)

which may depend on time explicitly and through a few variables such as the value and its derivative of the time-dependent boundary, b(t) and $\dot{b}(t)$, respectively, and possibly through a few auxiliary variables $z_i(t)$ that obey simple ordinary differential equations. Note that we use the notations $\varPhi$ for the boundary-dependent hazard rate of the model Eq. (5) and $\varPsi$ for the voltage-dependent hazard rate of the model Eq. (2b). The two functions are related in a simple way, see Sect. 4.1.

3.2 Wiener–Rice series

Our approach to tackle the time-dependent FPT problem is to employ the level-crossing statistics of a Gaussian process (Rice 1945; Ricciardi and Sato 1983; Verechtchaguina et al. 2006; Braun and Thul 2017; Azaïs and Wschebor 2009). To this end, let us consider the sub-set of all realizations of x(t) that cross the barrier b(t) from below in the time interval $(t^*,t^*+\varDelta t)$, a so-called “up-crossing” (Fig. 1b). The up-crossing at time $t^*$ is not necessarily the first one but could be the second, third (and so on) up-crossing (e.g., green and red lines in Fig. 1b). To compute the density of the first up-crossing, one can make use of the statistics of repeated up-crossing events. These events form a point process in the time interval $[0,t^*]$

$$\begin{aligned} s(t)=\sum _{i=1}^{N(t^*)}\delta (t-{\hat{t}}_i), \end{aligned}$$

(10)

where $N(t^*)$ denotes the (random) number of up-crossings in that interval, $\{{\hat{t}}_i\}_{i=1,\dotsc ,N(t^*)}$ are the up-crossing times and $\delta (\cdot )$ is the Dirac $\delta$-function. The statistics of the point process can be fully described by the set of moment functions $f_k(t_1,\dotsc ,t_k)=\langle s(t_1)\cdots s(t_k)\rangle$, for $k=1,2,\dotsc$ and non-coinciding time arguments $t_i$ (Stratonovich 1967a; van Kampen 1992). The moment functions can be interpreted such that for a small time step $\varDelta t$ the quantity $f_k(t_1,\dotsc ,t_k)\varDelta t^k$ yields the probability to find an up-crossing events in each of the non-overlapping intervals $(t_1,t_1+\varDelta t)$, ..., $(t_k,t_k+\varDelta t)$. For instance, $f_1(t)$ yields the rate of up-crossings at time t, and $f_2(t_1,t_2)/f(t_1)$ is the conditional rate of an upcrossing at time $t_2$ given an up-crossing at time $t_1$. For level-crossings of Gaussian processes, the distribution functions $f_k$ can be calculated explicitly, for both stationary and non-stationary processes (see Appendix, Sect. 3).

The distribution functions $f_k$ allow for an exact series expression of the FPT density, sometimes called Wiener–Rice series (Verechtchaguina et al. 2006; Braun and Thul 2017):

$$\begin{aligned} P(t)=\sum _{k=0}^\infty \frac{(-1)^k}{k!}\int _{0}^t\hbox {d}t_1\cdots \hbox {d}t_k\,f_{k+1}(t_1,\dotsc ,t_k,t) \end{aligned}$$

(11)

A detailed explanation of this formula is given in reference (Verechtchaguina et al. 2006). In brief, it counts—for a large ensemble of trajectories—the number of those trajectories that have a crossing in $[t,t+\hbox {d}t)$ but no crossing in (0, t). Starting with the fraction $f_1(t)\hbox {d}t$ of all trajectories that cross the boundary at time t ($k=0$ term), the fraction with no previous crossing can be computed by subtracting those trajectories that crossed the boundary before time t. The second term $\int _0^tf_2(t_1,t)\hbox {d}t_1$ in Eq. (11) accounts for these trajectories but overestimates their number because some trajectories are counted multiply. This corresponds to trajectories that cross the boundary more than once before time t (e.g., red line in Fig. 1). To correct for the excessive subtraction, one needs to add the fraction of trajectories with two or more crossings before t. This is taken into account by the third term $\frac{1}{2}\int _0^t\int _0^tf_3(t_1,t_2,t)\hbox {d}t_1\hbox {d}t_2$ which computes the mean number of crossing pairs $\{{\hat{t}}_1,{\hat{t}}_2\}$ per trajectory (e.g., in Fig. 1, the blue and green curve contributes zero and the red curve contributes one such pair; the factor $\frac{1}{2}$ accounts for permutations of ${\hat{t}}_1$ and ${\hat{t}}_2$). Again, this term overestimates the fraction of trajectories with double crossing events because trajectories with more than two crossings are multiply counted (e.g., a trajectory with three crossings gives rise to three pairs $\{{\hat{t}}_1,{\hat{t}}_2\}$, $\{{\hat{t}}_1,{\hat{t}}_3\}$, $\{{\hat{t}}_2,{\hat{t}}_3\}$). Continuing this correction procedure for trajectories with arbitrary number of crossings leads to the infinite series expression Eq. (11).

An alternative statistical description of the point process s(t) is given by the k-th-order cumulant functions $g_k(t_1,\dotsc ,t_k)$ (see (Stratonovich 1967a; van Kampen 1992) and Sect. 1), which remove the dependence on lower-order moment functions: for instance, $g_1(t)=f_1(t)$ and $g_2(t_1,t_2)=f_2(t_1,t_2)-f_1(t_1)f_1(t_2)$. The probability to find no event in the interval (0, t) (i.e., the survival probability) is related to the cumulant functions by (Stratonovich 1967a; van Kampen 1992)

$$\begin{aligned} S(t)=\exp \left( \sum _{k=1}^\infty \frac{(-1)^k}{k!}\int _{0}^t\mathrm {d}t_1\cdots \mathrm {d}t_k\,g_k(t_1,\dotsc ,t_k)\right) .\nonumber \\ \end{aligned}$$

(12)

From this expression, the Wiener–Rice series for the FPT density, Eq. (11) is recovered by $P(t)=-\hbox {d}S(t)/\hbox {d}t$. Similarly, the hazard rate can be obtained by $\lambda (t)=-\hbox {d}(\ln S)/\hbox {d}t$. As infinite series expressions, Eq. (11) and Eq. (12) are of no practical use for direct computations of the FPT density. However, these formal expressions are used as a starting point for further approximations.

3.3 Decoupling approximations

The series expression for the survival probability, Eq. (12), simplifies considerably if higher-order cumulant functions $g_k$ are approximated in terms of lower-order cumulant functions, thereby neglecting higher-order dependencies between up-crossings. In this section, we review two approximations based on such a decoupling of (temporal) interactions between events (Stratonovich 1967a): a first-order decoupling approximation, where all up-crossing events are assumed to be independent, and a second-order decoupling approximation, in which higher-order interactions are modeled in terms of pairwise interactions. While the first-order approximation readily results in local hazard rates, the more accurate pairwise interaction approximation is highly non-local and therefore not useful for practical calculations. However, as we shall show in Sect. 3.5, the pairwise interaction model can be used as a starting point for deriving an efficient local approximation of the hazard rate (second-order decoupling approximation) that accounts for higher-order interactions between up-crossings.

3.3.1 Independent upcrossings

If the correlation time of the process x(t) is much smaller than the (typical) intervals between upcrossings, up-crossing events can be regarded as independent, i.e., the series of up-crossing events is an inhomogeneous Poisson process with rate $f_1(t)$. Mathematically, this corresponds to neglecting higher-order cumulants except for the first one: $g_1(t)=f_1(t)$ and $g_k\approx 0$ for all $k\ge 2$ (Stratonovich 1967a). In this case, Eq. (12) reduces to $S(t)=\exp \left( -\int _{0}^{t}f_1(\tau )\,\hbox {d}\tau \right)$, and hence the FPT density reads

$$\begin{aligned} P(t)\approx f_1(t)\exp \left\{ -\int _{0}^{t}f_1(\tau )\,\hbox {d}\tau \right\} . \end{aligned}$$

(13)

From this expression, we see that the hazard rate is simply given by the upcrossing rate of the freely evolving process x(t): $\lambda (t)\approx f_1(t)$. The upcrossing rate $f_1(t)$ can be calculated analytically in terms of the current value of the boundary b(t) and its derivative $\dot{b}(t)$ (see Appendix 3 and 4). The result is the first-order decoupling approximation:

$$\begin{aligned} \lambda (t)\approx f_1(t)= & {} \varPhi _1\bigl (b(t),\dot{b}(t),t\bigr )\nonumber \\:= & {} \frac{\sqrt{\sigma _x^2\sigma _y^2-\sigma _{xy}^2}}{2\pi \sigma _x^2}H\left( \frac{(\gamma \sigma _x^2-\sigma _{xy})b+\sigma _x^2\dot{b}}{\sqrt{2(\sigma _x^2\sigma _y^2-\sigma _{xy}^2)}\sigma _x}\right) e^{-B(b,\dot{b},t)}, \nonumber \\ \end{aligned}$$

(14)

where $H(x)=1-\sqrt{\pi }xe^{x^2}\text {erfc}(x)$ and

$$\begin{aligned} B(b,\dot{b},t)&=\frac{\left( \gamma ^{2}\sigma _x^2-2\gamma \sigma _{xy}+\sigma _y^2\right) b^2+2\left( \gamma \sigma _x^2-\sigma _{xy}\right) b\dot{b}+\sigma _x^2\dot{b}^2}{2\left( \sigma _x^2\sigma _y^2-\sigma _{xy}^2\right) }. \end{aligned}$$

(15)

In these equations, the time-dependent moments $\sigma _{xy}(t)=\langle x(t)y(t)\rangle$ and $\sigma _x^2(t)=\langle x^2(t)\rangle$ are given by

$$\begin{aligned} \sigma _{xy}(t)&={\tilde{\tau }}\sigma _y^2\left( 1-e^{-t/{\tilde{\tau }}}\right) , \end{aligned}$$

(16a)

$$\begin{aligned} \sigma _x^2(t)&=\frac{{\tilde{\tau }}\sigma _y^2}{\gamma }\left( 1-e^{-2\gamma t}\right) +\frac{2{\tilde{\tau }}\sigma _y^2}{2\gamma -{\tilde{\tau }}^{-1}}\left( e^{-2\gamma t}-e^{-\frac{t}{{\tilde{\tau }}}}\right) \end{aligned}$$

(16b)

with $\sigma _y^2=D/\tau _y$ and ${\tilde{\tau }}^{-1}=\gamma +\tau _y^{-1}$ (see Sect. 2, esp. Eq. (59) for a numerically stable ODE representation of the moments).

Similar expressions for the level-crossing density in the time-inhomogeneous case have been derived in previous studies (Ricciardi and Sato 1983; Badel 2011).

3.3.2 Upcrossings correlated in pairs

If the average time between upcrossings $1/f_1(t)$ is on the order of or smaller than the correlation time of x(t) given by $\tau _\mathrm{cor}=\gamma ^{-1}+\tau _y$, upcrossing events cannot be regarded as being independent anymore. To account for correlations between upcrossings, we follow a decoupling approximation (DA) of higher-order correlation functions $g_k(t_1,\dotsc ,t_k)$, $k\ge 3$, proposed by Stratonovich (Stratonovich 1967a, b). This approximation assumes that higher-order correlations are governed by the same time scales as pair-wise correlations and can therefore be expressed in terms of the first two correlation functions $f_1(t)$ and $g_2(t_1,t_2)$. Specifically, correlation functions with $k\ge 2$ are approximated by the ansatz (Stratonovich 1967a, b)

$$\begin{aligned} g_k(t_1,\dotsc ,t_k)=(k-1)!f_1(t_1)\cdots f_1(t_k)\{R(t_2,t_1)\cdots R(t_k,t_1)\}_\mathrm{sym}.\nonumber \\ \end{aligned}$$

(17)

here the function $R(t,t')$ describes the pairwise interactions between events at time t and $t'$, and $\{\cdots \}_\mathrm{sym}$ denotes the operation of symmetrization (i.e., the arithmetic mean of all permutations of the time arguments). As suggested in (Stratonovich 1967a, b), we choose $R(t,t')$ as the normalized auto-correlation function

$$\begin{aligned} R(t,t')=\frac{f_2(t,t')}{f_1(t)f_1(t')}-1, \end{aligned}$$

(18)

which makes the ansatz Eq. (17) exact for $k=2$. Note that compared to (Stratonovich 1967a, b), we use an opposite sign in the definition of R for mathematical convenience. The auto-correlation function $R(t,t')$ can be interpreted as the conditional probability density of an event at time $t'$ given an event at time t normalized by the unconditional probability density $f_1(t')$ and shifted by the mean such that $R(t,t')=0$ if events at time t and $t'$ are independent. For stationary point processes, $R(t,t')=R(|t-t'|)$ only depends on the time difference. In analogy to the common use for spatial point processes, $R(t-t')$ will be called pair correlation function in this case.

We expect the following behavior of the auto-correlation function: firstly, if events are far apart, $|t-t'|\gg \tau _\text {cor}$, they occur independently, hence $f_2(t,t')\approx f_1(t)f_1(t')$. This implies a vanishing auto-correlation function $R(t,t')\approx 0$. Secondly, the behavior when t and $t'$ are close depends on the correlations between events: if close events occur independently as in the case of an inhomogeneous Poisson process, $R(t,t')$ vanishes. In contrast, a positive pair correlation function $R(t,t')>0$ at small time lag indicates that events are attractive and tend to cluster. Conversely, for a negative pair correlation function $R(t,t')<0$ at small time lag, events are repulsive, i.e., the occurrence of close events is less frequent than expected for a Poisson process. In particular, if a point process exhibits a refractory period after each event (“hardcore interaction”), we find that $f_2(t,t')=0$ and hence $R(t,t')=-1$ if t and $t'$ fall within a refractory period. Similarly, non-approaching random points (Stratonovich 1967a) are characterized by $R(t,t)=-1$ in the limit of vanishing time lag. Interestingly, it has been assumed by some authors that level crossings of differentiable processes are non-approaching events with $R(t,t)=-1$ (Verechtchaguina et al. 2006; Puelma Touzel and Wolf 2016). In Sect. 3.4, we shall investigate this assumption in more detail.

While the decoupling approximation (DA), Eq. (17), is exact for $k=2$ by construction, it must be considered as a physically motivated, heuristic ansatz for $k\ge 3$, which in general is not expected to be exact. Nevertheless, the ansatz and the above-described behavior of $R(t,t')$ ensure some important properties of the higher-order correlation functions $g_k$: first, the DA is exact for an inhomogeneous Poisson process because in this case $R(t,t')\equiv 0$ and thus Eq. (17) recovers the expected result $g_k\equiv 0$ for all $k\ge 2$. Second, $g_k$ does not depend on the order of the time arguments because of the symmetrization operation in Eq. (17). Third, $g_k(t_1,\dotsc ,t_k)\approx 0$ if the time difference of two arguments is much larger than $\tau _\text {cor}$ because their pair correlation vanishes. And forth, it is known that for a system of non-approaching random points $g_k(t,\dotsc ,t)=(-1)^k(k-1)!f_1^k(t)$ (Stratonovich 1967b), which is consistent with Eq. (17) and $R(t,t)=-1$.

Substituting the DA, Eq. (17), into the general expression for the survival probability, Eq. (12), yields (Stratonovich 1967a, b; Verechtchaguina et al. 2006; van Meegen and van Albada 2019)

$$\begin{aligned} S(t)\approx \exp \left\{ -\int _{0}^{t}f_1(\tau )\frac{\ln \left[ 1+q(t,\tau )\right] }{q(t,\tau )}\,d\tau \right\} , \end{aligned}$$

(19)

where

$$\begin{aligned} q(t,\tau )= & {} \int _{0}^tR(\tau ,\tau ')f_1(\tau ')\,\hbox {d}\tau '\nonumber \\= & {} \frac{1}{f_1(\tau )}\int _0^t\left[ f_2(\tau ,\tau ')-f_1(\tau )f_1(\tau ')\right] \,\hbox {d}\tau ' \end{aligned}$$

(20)

is a measure of upcrossing correlations on the time scale t. The formula Eq. (19) has been termed Stratonovich approximation (Verechtchaguina et al. 2006). Comparing the Stratonovich approximation with the first-order decoupling approximation, Eq. (13), we observe that the upcrossing rate $f_1(\tau )$ is multiplied by a correction factor $\ln (1+q)/q$. However, this correction factor depends explicitly on time t, which precludes a direct interpretation of the integrand in Eq. (19) as the hazard rate (but see (van Meegen and van Albada 2019) for a hazard rate approximation of the integrand in the time-homogeneous case). For the Stratonovich approximation to be applicable, one has to require that

$$\begin{aligned} q(t,\tau )>-1 \end{aligned}$$

(21)

for all t and $\tau$ so as to keep the argument of the logarithm positive (Verechtchaguina et al. 2006).

In practice, Eq. (19) is not useful as a computational tool. A numerical evaluation is highly inefficient because Eq. (19) contains nested integrals on three levels: for each $\tau$ of the outer integral, the integral $q(t,\tau )$ needs to be evaluated independently for each time t. Furthermore, the numerical integration of $q(t,\tau )$ is itself computationally complex because $R(\tau ,\tau ')$ involves a further integration (taking already into account that one of the two integrals in the definition of $f_2$, Eq. (81), Sect. 5, can be evaluated analytically (Ricciardi and Sato 1983; Verechtchaguina et al. 2006); we note that $f_2$ can also be expressed in terms of Owen’s T function (van Meegen and van Albada 2019)). Therefore, we will further simplify Eq. (19) by deriving a local approximation of the hazard rate.

3.4 The auto-correlation function of level crossings for small time lags

We now proceed with calculating the auto-correlation function $R(t,t+\tau )$ in the limit of small time lags $\tau$. Based on the zero-lag limit we then propose a rough estimation of the temporal correlation structure for $\tau >0$, which will be required for the simplification of the Stratonovich approximation in the next section. While the rate of level-crossings has been studied extensively (e.g., (Rice 1945; Stratonovich 1967b; Verechtchaguina et al. 2006; Tchumatchenko et al. 2010b)), the calculation of second-order statistics such as the auto-correlation function has not received much attention. To the best of our knowledge, closed-form analytical formulas for the auto-correlation function of non-stationary level crossings have not been published previously. In Appendix Sect. 2, we also provide formulas for the auto-correlation function of general Gaussian level-crossing processes in the stationary state (see also (Burak et al. 2009) for special cases and (Jung 1994) for the related but distinct result for the stationary auto-correlation function of the two-state process triggered by level crossings).

According to Eq. (18), the auto-correlation function at zero time lag is given by

$$\begin{aligned} R_0(t)=\frac{f_2(t,t)}{f_1^2(t)}-1, \end{aligned}$$

(22)

where $f_2(t,t)\equiv \lim _{\tau \rightarrow 0}f_2(t,t+\tau )$ is defined through the limiting procedure $\tau \rightarrow 0$. This corresponds to the continuous part of the auto-correlation function, i.e., $f_2(t,t)$ excludes the singular self-correlation of points given by $f_1(t)\delta (\tau )$. The correlations between upcrossing in the limit of vanishing lag can be calculated within a saddle-point approximation (see Appendix, Sect. 5). The result is

$$\begin{aligned} f_2(t,t)&=\frac{3\sqrt{3}-\pi }{36\pi ^2}\frac{\sigma _y^2/\tau _y}{\sqrt{\sigma _x^2\sigma _y^2-\sigma _{xy}^2}}e^{-B(b,\dot{b},t)} \end{aligned}$$

(23)

$$\begin{aligned}&=:{{\hat{f}}}_2\bigl ( b(t),\dot{b}(t),t \bigr ) \end{aligned}$$

(24)

It is instructive to discuss the stationary case, $\dot{b}=0$ and $t\rightarrow \infty$, in which the pair correlation function $R(\tau )$ for vanishing time lag $\tau$ obtains the simple form

$$\begin{aligned} R_0=\beta \frac{1+\gamma \tau _y}{\sqrt{\gamma \tau _y}}\exp \left( \frac{b^2}{2\sigma _x^2} \right) -1 \end{aligned}$$

(25)

with the numerical constant $\beta =(3\sqrt{3}-\pi )/9\approx 0.228284$. For any fixed value of $\gamma \tau _y$, this expression becomes minimal at $b=0$ (Fig. 2c, blue dashed line). From this, we infer that $R_0$ is always positive if $\gamma \tau _y<0.0583757$ (“white noise regime”) or $\gamma \tau _y>17.1304$ (strong friction or large noise correlation time). In this case, upcrossings tend to cluster. In the wide intermediate range $0.0583757<\gamma \tau _y<17.1304$, the sign of $R_0$ depends on the ratio $|b|/\sigma _x$ of barrier height to standard deviation of x(t). For vanishing or low barrier height such that $|b|/\sigma _x$ is below the critical value

$$\begin{aligned} \frac{b_\mathrm{crit}}{\sigma _x}=\sqrt{2\ln \left( \frac{\sqrt{\gamma \tau _y}}{\beta (1+\gamma \tau _y)} \right) }, \end{aligned}$$

(26)

the pair correlation function will be negative at small time lags, i.e., upcrossings tend to repel each other. Closer inspection of Eq. (25) shows that for any barrier height b, $R_0$ becomes minimal (i.e., most negative) if $\gamma \tau _y=1$. The absolute achievable minimum is found as $R_0=-0.543431$. Therefore, the value $R_0=-1$ expected for non-approaching points is never realized for level crossings of a doubly low-pass-filtered white noise such as Eq. (1) and (5) for the membrane potential and overdamped Brownian particle driven by a one-dimensional Ornstein–Uhlenbeck noise, respectively. This result is in marked contrast to the assumption of non-approaching level crossings made in previous studies (Verechtchaguina et al. 2006; Puelma Touzel and Wolf 2016).

Fig. 2

Correlations of level crossings of a stationary process x(t). a Normalized auto-correlation function $R(\tau )\equiv R(t,t+\tau )$ as a function on the time lag $\tau$ (in units of $\tau _1{\mathop {=}\limits ^{\text {\tiny def}}}\tau _\text {m}=\gamma ^{-1}$, $\tau \ne 0$) for constant barriers b (as indicated on top) and small time constant $\tau _y=0.4\tau _\text {m}$. The solid magenta lines show the exact semi-analytical result obtained from numerical integration of Eq. (81), and the blue dashed lines show the exponential approximation, Eq. (30), respectively. b Same as a but with $\tau _y=2.5\tau _\text {m}$. c Correlations in the limit of vanishing time lag, $R(0)=\lim _{\tau \rightarrow 0}R(\tau )$, as a function of the time scale ratio $\tau _2/\tau _1=\tau _y/\tau _\text {m}$ for three different constant ($\dot{b}=0$) threshold levels b (as indicated). d Correlations for vanishing time lag as a function of the instantaneous threshold level b(t) for three different slopes $\dot{b}(t)$ (at $\tau _y=0.4\tau _\text {m}$): decreasing thresholds lower probability of observing two infinitesimally close level crossings (blue dashed line), whereas increasing threshold increase this probability (finely dashed red line) compared to constant thresholds (solid green line). In all panels, black dotted lines indicate the zero baseline corresponding to a Poisson statistics

On the other hand, for high barriers such that $|b|>b_\mathrm{crit}$, the pair correlation function is positive at small time lags, implying that upcrossing events tend to cluster. Intuitively, upcrossings are mediated by large fluctuations of x(t) in order to reach the high barrier. Once the barrier is reached, x(t) persists at high values for some period because values of x(t) are positively correlated at short time lags. During this period, the probability to cross the barrier for a second time is strongly increased. That is, upcrossings tend to cluster in periods on the order of the correlation time of x(t). This clustering corresponds to a positive pair correlation $R_0$

3.5 Local hazard function

From the Stratonovich approximation, Eq. (19), we obtain the corresponding hazard rate by differentiating $-\ln S(t)$ with respect to t. Using Eq. (20), the result can be written as

$$\begin{aligned} \lambda (t)=f_1(t)\left\{ F\bigl (q(t,t)\bigr )+\int _{0}^t f_1(\tau )F'\bigl (q(t,\tau )\bigr )R(t,\tau )\,\mathrm {d}\tau \right\} \nonumber \\ \end{aligned}$$

(27)

where $F(q)=\ln (1+q)/q$. Because of the integral in Eq. (27), the hazard rate is still non-local in time. In order to obtain a local approximation, we make two ad hoc approximations. First, Eq. (27) can be considerably simplified if $F'(q(t,\tau ))$ only weakly depends on $\tau$ such that we can pull this function out of the integral. Under this assumption and using again Eq. (20), the hazard rate reduces to the particularly simple form

$$\begin{aligned} \lambda (t)=\frac{f_1(t)}{1+q(t)}, \end{aligned}$$

(28)

where we used the short-hand notation

$$\begin{aligned} q(t)\equiv q(t,t)=\int _{0}^tR(t,t')f_1(t')\,\hbox {d}t'. \end{aligned}$$

(29)

The above ad hoc approximation seems plausible because the pair-correlation function $R(t,\tau )$ is different from zero only in a region of width $|\tau -t|\sim \tau _\mathrm{corr}$ around its peak at the integration boundary $\tau =t$, where $\tau _\mathrm{corr}$ is the correlation time defined in Eq. (31) (Fig. 2a, b). On this time scale, $q(t,\tau )$ represents indeed a slowly varying function of $\tau$ since it results from an integration over R (cf. Eq. (20)). Note that an alternative approximation has been suggested in (van Meegen and van Albada 2019), which neglects the second term in Eq. (27).

The formula Eq. (28) reveals a simple relation between the upcrossing rate and the hazard rate, which is the relevant quantity for the FPT: In the absence of correlations between upcrossings, $q=0$, the two rates are equal, while negative correlations (repulsion of up-crossings) increase the hazard rate and positive correlations (attraction or clustering of up-crossings) decreases the hazard rate compared to the up-crossing rate $f_1$.

Second, to find a local estimation of q(t) we need to turn the integral in Eq. (29) into a differential equation for q. A simple way to achieve this is to use an exponential approximation for the pair correlation function

$$\begin{aligned} R(t,t')\approx R_0(t)\exp \left( -\frac{|t-t'|}{\tau _\mathrm{corr}} \right) , \end{aligned}$$

(30)

where $R_0(t)=f_2(t,t)/f_1^2(t)-1$ is the limit of vanishing time lag $\tau \rightarrow 0$. Accordingly, the function $f_2(t,t)$ has to be understood as the limit $\lim _{\tau \rightarrow 0}f_2(t,t+\tau )$, which has been calculated analytically in the previous section. Furthermore, $\tau _{corr}$ is the typical correlation time with which correlations between upcrossings decay as function of their temporal distance. As a rough approximation, this correlation time is given by the correlation time of the stationary process x(t) itself:

$$\begin{aligned} \tau _\mathrm{corr}=\int _0^\infty \frac{C_{xx}(\tau )}{C_{xx}(0)}\,\mathrm {d}\tau =\tau _\text {m}+\tau _\text {s}. \end{aligned}$$

(31)

here $C_{xx}(\tau )$ is the auto-correlation function of x(t) in the stationary state. In fact, comparison of the exponential approximation with numerical evaluation of the exact quadrature formula of the correlation function confirms our choice of $\tau _\mathrm{corr}$ and also shows that that the exponential ansatz is reasonable as long as $R_0$ is significantly different from zero (Fig. 2 a,b, left and right panels). In the crossover region from negative to positive $R_0$ when the barrier height b is increased, the auto-correlation function has both positive and negative phases that are not captured by an exponential function (Fig. 2 a,b, middle panels). However, these deviations are less significant because absolute correlations are small in this case.

Inserting the exponential ansatz Eq. (30) into Eq. (29), we can pull $R_0(t)$ in front of the integral and obtain:

$$\begin{aligned} q(t)\approx R_0(t)z(t), \end{aligned}$$

(32)

where $z(t)=\int _0^t\exp \left[ -(t-t')/\tau _\mathrm{corr} \right] f_1(t')$ defines a new auxiliary variable that satisfies the differential equation

$$\begin{aligned} \frac{\mathrm {d}z}{\mathrm {d}t}=-\frac{1}{\tau _\mathrm{corr}}z+f_1(t) \end{aligned}$$

(33)

with $z(0)=0$. We note that the slightly different ansatz $R(t,t')\approx R_0(t')\exp \left( -\frac{t-t'}{\tau _\mathrm{corr}} \right)$ yields slightly different equations with similarly good results. In Sect. 5.2, we will thus only show the results for the above ansatz, Eq. (30).

We note that in the limit of vanishing correlations between upcrossings, $R_0(t)\equiv 0$, the first-order DA $\lambda (t)\approx f_1(t)$ is recovered from Eq. (28). Thus, the first-order approximation, Eq. (14), is expected to be valid if

$$\begin{aligned} |q(t,t)|\ll 1 \end{aligned}$$

(34)

for all $t>0$.

In summary, the local hazard rate in the second-order DA is given by

$$\begin{aligned} \lambda (t)\approx \varPhi _2\bigl (b(t),\dot{b}(t),z(t),t\bigr ){\mathop {=}\limits ^{\text {\tiny def}}}\frac{\varPhi _1\bigl ( b(t),\dot{b}(t),t\bigr )}{1+{{\hat{R}}}_0\bigl ( b(t),\dot{b}(t),t \bigr )z(t)}.\nonumber \\ \end{aligned}$$

(35)

here $\varPhi _1$ is given by Eq. (14) and

$$\begin{aligned} {\hat{R}}_0\bigl ( b,\dot{b},t \bigr )=\frac{{\hat{f}}_2\bigl ( b,\dot{b},t \bigr )}{\left[ \varPhi _1\bigl ( b,\dot{b},t \bigr ) \right] ^2}-1 \end{aligned}$$

(36)

is the zero-lag correlation between up-crossings, Eq. (22), where ${\hat{f}}_2$ is given by Eq. (23). In contrast to the first-order approximation $\varPhi _1$, the hazard rate $\varPhi _2$ depends on the additional local variable z that obeys

$$\begin{aligned} \frac{\mathrm {d}z}{\mathrm {d}t}=-\frac{1}{\tau _\mathrm{corr}}z+\varPhi _1\bigl ( b(t),\dot{b}(t),t \bigr ),\quad z(0)=0. \end{aligned}$$

(37)

Together with Eq. (8), this ordinary differential equation provides an update rule for the numerical evaluation of the FPT density P(t) forward in time.

3.6 First-passage-time densities

Fig. 3

First-passage-time density for periodically moving barrier. a Low amplitude $\alpha =0.25$ (subthreshold regime). Top: illustration of moving barrier (green dashed line) and a sample trajectory x(t) (black solid line). The shaded region indicates the mean $\langle x\rangle =0$ (horizontal dashed line) ± the standard deviation $\sigma _x(t)$. Bottom: first-passage-time density P(t) from simulations (gray circles) and theory (first- and second-order decoupling approximation—Eq. (13) (green dashed line) and Eq. (7) (blue solid line), respectively; and the Chizhov–Graham theory—Eqs. (96)–(101) (red thin line)). b Same with high amplitude $\alpha =1.2$ (suprathreshold regime). Parameters: $\sigma _x(\infty )=0.5$, $\tau _x=1$, $\tau _y=0.2$, $f=0.5$

Being equipped with local approximations of the hazard rate, the FPT density P(t) can be easily obtained from Eq. (8). To test the performance of our theory, we compare the first- and second-order decoupling approximations (DA) with simulations and an alternative hazard-rate theory proposed by Chizhov and Graham (2008). An extended variant of the Chizhov–Graham (C&G) theory is presented in Appendix 3, Eq. (101).

For concreteness, we consider a periodically moving boundary:

$$\begin{aligned} b(t)=1+\alpha \cos (2\pi ft) \end{aligned}$$

(38)

(Fig. 3, top panels). The case, where the amplitude of the oscillating boundary is smaller than unity, $\alpha <1$, corresponds to the subthreshold firing regime of LIF neurons. In this case, both the first- and second-order DA (Eq. (8) with $\lambda (t)$ given by Eq. (14) and (35), respectively) yield excellent agreements with simulations (Fig. 3a). In contrast, the C&G theory (Eq. (8) with $\lambda (t)$ given by Eq. (101)), shows clear deviations from simulations at the peaks of the FPT density and during the time spans when the boundary is increasing ($\dot{b}>0$), i.e., when the boundary moves away from zero. In these regions, the drift component, Eq. (96), of the C&G hazard rate is set to zero, leaving only diffusion as a source of threshold crossings. The rectification of the drift component also leads to a characteristic kink at the local extrema of the boundary ($\dot{b}(t)=0$).

The case of large amplitude oscillations of the boundary ($\alpha >1$) is equivalent to a LIF model that is periodically driven into the supra-threshold regime. In this case, the first-order DA performs significantly worse than the second-order approximation and the C&G theory, which both agree well with simulation results (Fig. 3b). In particular, the first peak in the FPT density (green dotted line in Fig. 3b) is underestimated if correlations between upcrossings are neglected. The underestimation is caused by a reduced hazard rate, which can be understood from the simple formula Eq. (28): in the first order approximation, the hazard rate is given by the level-crossing rate $\lambda (t)\approx f_1(t)$, while in the second-order approximation $\lambda (t)\approx f_1(t)/[1+q(t)]$ with $q(t)=R_0(t)z(t)$. The factor $1/(1+q)$ accounts for the correlations between upcrossings. At the peak, the boundary b(t) is close to zero. In this case, the zero lag pair correlation $R_0$ is negative representing the reduced probability of nearby crossings (“repulsion,” Fig. 2, left panels). Since z is positive, we have $-1<q<0$ and thus the factor $1/[1+q]$ is larger than unity (note that $q>-1$ by the assumption Eq. (21)). Therefore, correlations between upcrossings lead to an increased hazard rate and thus a stronger first peak of the FPT density.

4 Mapping colored input noise to escape noise in the leaky integrate-and-fire model

4.1 Link function

We now come back to our initial motivation to map colored noise in the input to escape noise in the output of a LIF neuron. Having derived the hazard rate $\varPhi$ for the FPT with moving boundary b(t), it is easy to formulate the link function $\varPsi$ in Eq. (2) that provides the escape-noise model corresponding to the LIF model with input noise Eq. (1). To this end, we only need to shift time such that the FPT starts at time ${\hat{t}}+t_\text {ref}$ instead of $t=0$, enforce a zero hazard rate during the absolute refractory period, and express the moving threshold b(t) in terms of the mean membrane potential u(t) for $t>{\hat{t}}+t_\text {ref}$ using Eq. (4). Accordingly, we also replace the temporal derivative of the moving boundary by

$$\begin{aligned} \dot{b}(t)=-\dot{u}(t)=\frac{u(t)-\mu (t)}{\tau _\text {m}} \end{aligned}$$

(39)

for $t>{\hat{t}}+t_\text {ref}$. The last expression shows that, instead of the two functions u(t) and $\dot{u}(t)$, one can also use the two functions u(t) and $\mu (t)$ if the stimulus $\mu (t)$ is known.

With these changes, we obtain the link function in the first-order DA as

$$\begin{aligned} \varPsi _1\bigl ( u,\dot{u},\tau \bigr )=\theta (\tau -t_\text {ref})\varPhi _1(V_\text {T}-u,-\dot{u},\tau -t_\text {ref}). \end{aligned}$$

(40)

here $\theta (t)=\mathbf {1}_{t\ge 0}$ is the Heaviside step function and $\varPhi _1$ is given by Eq. (14). Note that in the first-order DA, the link function $\varPsi (u,\dot{u},z,\tau )=\varPsi _1(u,\dot{u},\tau )$ does not depend on an auxiliary variable z. In contrast, the second-order DA exhibits an additional auxiliary variable z. Taking the last spike time and the absolute refractory period into account, its dynamics reads

$$\begin{aligned} \dot{z}=-\frac{z}{\tau _\text {m}+\tau _\text {s}}+\varPsi _1(u,-\dot{u},t-{\hat{t}}) \end{aligned}$$

(41)

with initial condition $z({\hat{t}})=0$. We can now write the link function $\varPsi$ in the second-order DA as

$$\begin{aligned} \varPsi _2\bigl ( u,\dot{u},z,\tau \bigr )=\theta (\tau -t_\text {ref})\varPhi _2(V_\text {T}-u,-\dot{u},z,\tau -t_\text {ref}),\nonumber \\ \end{aligned}$$

(42)

where $\varPhi _2$ is given by Eq. (35).

Fig. 4

First-passage-time density, survivor function and hazard rate under non-stationary driving of a neuron that fired its last spike at time ${\hat{t}}=0$. a Weak subthreshold stimulus $\mu (t)$ (top panel) leads to a mean membrane potential response u(t|0) below threshold at $V_\text {T}=1$ (second panel). The first-passage-time density P(t|0) for the first threshold crossing of ${{\hat{V}}}(t)$ is shown in the third panel (gray circles: MC simulations of $10^6$ trials; red solid line: Chizhov–Graham theory, Eq. (7), (101); blue dashed line: first-order decoupling approximation (independent up-crossings), Eq. (45), (43); blue solid line: second-order decoupling approximation (correlated upcrossings), Eq. (46), (43). The survival probability $S(t|0)=-dP(t|0)/\hbox {d}t$ and the corresponding hazard rate $\lambda (t|0)$ are shown in the two bottom panels, respectively. For MC simulations, the hazard rate is computed from the ratio $\lambda (t|0)=P(t|0)/S(t|0)$. b The same for a suprathreshold stimulus, for which the mean membrane potential u(t|0) reaches the threshold. In both figures, $\tau _\text {s}=4$ ms, $\tau _\text {m}=10$ ms and $\sigma _\eta$ is such that the standard deviation of ${{\hat{V}}}$ is $\sigma _V=0.25$

4.2 Comparison with simulation and C&G theory

To judge the performance of the level-crossing theory given by the link functions $\varPsi _1$ and $\varPsi _2$, we compared ISI densities, survival functions and hazard rates with Monte-Carlo simulations of the LIF model with colored input noise, Eq. (1), and the C&G theory. These functions are obtained from the link functions as

$$\begin{aligned} P(t|{\hat{t}})&=\lambda (t|{\hat{t}})S(t|{\hat{t}}), \end{aligned}$$

(43)

$$\begin{aligned} S(t|{\hat{t}})&=\exp \left( -\int _{{\hat{t}}}^t\lambda (s|{\hat{t}})\,\hbox {d}s \right) , \end{aligned}$$

(44)

where for the first-order decoupling approximation(DA)

$$\begin{aligned} \lambda (t|{\hat{t}})\approx \varPsi _1\bigl ( u(t|{\hat{t}}),\dot{u}(t|{\hat{t}}),t-{\hat{t}} \bigr ), \end{aligned}$$

(45)

and for the second-order DA,

$$\begin{aligned} \lambda (t|{\hat{t}})&\approx \varPsi _2\bigl ( u(t|{\hat{t}}),\dot{u}(t|{\hat{t}}),z(t|{\hat{t}}),t-{\hat{t}} \bigr ) \end{aligned}$$

(46)

with $\varPsi _1$ and $\varPsi _2$ given by Eq. (40) and Eq. (42), respectively. In Eq. (45) and (46), we have introduced the membrane potential and the auxiliary variable as deterministic functions of t and ${\hat{t}}$. For $t>{\hat{t}}+t_\text {ref}$, these functions obey the first-order dynamics

$$\begin{aligned} \dot{u}(t|{\hat{t}})&=-\frac{u(t|{\hat{t}})+\mu (t)}{\tau _\text {m}}, \end{aligned}$$

(47)

$$\begin{aligned} \dot{z}(t|{\hat{t}})&=-\frac{z(t|{\hat{t}})}{\tau _\text {m}+\tau _\text {s}}+\varPsi _1\Bigl (u(t|{\hat{t}}),-\dot{u}(t|{\hat{t}}),t-{\hat{t}}\Bigr ). \end{aligned}$$

(48)

with initial conditions $u({\hat{t}}+t_\text {ref}|{\hat{t}})=V_\text {R}$ and $z({\hat{t}}+t_\text {ref}|{\hat{t}})=0$.

The time-dependent stimulus $\mu (t)$, shown in Fig. 4 (top panels), was obtained as $\mu (t)=\mu _0+\mu _1(t)$, where $\mu _1(t)$ is a fixed realization of a band-limited white-noise process with a cut-off frequency of $100\text { Hz}$. Without loss of generality, we also choose the last spike time as the time origin, ${\hat{t}}=0$. The membrane potential u(t|0) is shown in Fig. 4 (second panel from top). Note that in simulations and figures, we measured voltages in units of $V_\text {T}-V_\text {R}$ and chose the arbitrary reference potential such that $V_\text {R}=0$, and hence $V_\text {T}=1$. For subthreshold stimuli (Fig. 4A), $u(t|0)<V_\text {T}$, both the first- and second-order decoupling approximations agree well with the interval distribution obtained from simulations of the model with colored input noise. As in the case of periodic subthreshold driving (Fig. 3a), the C&G theory exhibits again clear deviations at the peaks of the ISI density and in periods where the slope of the mean membrane potential is negative, $\dot{u}(t|0)<0$ (Fig. 4a, middle panel). The overall performance is better visible in the survival function (Fig. 4a, second panel from bottom), which is related to the cumulative ISI distribution via $S(t|{\hat{t}})=1-\int _{{\hat{t}}}^tP(s|{\hat{t}})\,\hbox {d}s$. It confirms the excellent performance of both decoupling approximations in the subthreshold regime. For completeness, we also compared the hazard rates (Fig. 4a, bottom panel). Note that the initial transient of u(t|0) from reset to resting potential $\mu _0$ realizes a relative refractory period, where the probability to fire is low.

For suprathreshold stimuli, where the mean membrane potential exceeds the threshold, the first-order DA deviates significantly from simulation results (Fig. 4b). This is because level crossings occur more frequently when u is close to the threshold and thus exhibit stronger (negative) correlations. In this case, the assumption of independent upcrossing is no longer valid. Again, the underestimation of the first peak in the ISI density and the hazard rate (dotted lines in Fig. 4b, middle and bottom panel) if correlations are neglected can be understood from the simple formula Eq. (28): under the assumption of independent upcrossings, the hazard rate is given by the level-crossing rate $\lambda (t|0)\approx f_1(t)$, while correlations between upcrossings are accounted for in the second-order approximation as $\lambda (t|0)\approx f_1(t)/[1+R_0(t)z(t|0)]$. We have seen that if u is close to the threshold (corresponding to $b=0$), the zero lag pair correlation $R_0$ is negative representing the reduced probability of nearby crossings (“repulsion,” Fig. 2, left panels). Since z is positive, the factor $1/[1+R_0z]$ is larger than unity (Note that $q\equiv R_0z>-1$ by assumption Eq. (21) for the applicability of the Stratonovich approximation). Therefore, correlations between upcrossings lead to an increased hazard rate (2nd-order DA) as compared to the theory with independent upcrossings (1st-order DA) (blue solid vs. blue dotted line in Fig. 4b, bottom).

Fig. 5

Error of the theoretical approximations for different stimulus properties. The error is measured as the Kolmogorov–Smirnov distance ${\mathcal {D}}$ between the theoretical and simulated ISI distribution. The stimulus $\mu (t)$ driving the LIF model is sampled from an Ornstein–Uhlenbeck process with mean ${{\bar{\mu }}}$, standard deviation $\sqrt{1+\tau _\text {m}/\tau _\mu }{{\bar{\sigma }}}$ and correlation time $\tau _\mu$. a Color-coded value of ${\mathcal {D}}$ as a function of ${{\bar{\mu }}}$ and ${{\bar{\sigma }}}$ for a rapidly varying stimulus, $\tau _\mu =1$ ms (left: 1st-order DA , middle: 2nd-order DA, right: Chizhov–Graham theory). b Same as a but for a moderately fast stimulus, $\tau _\mu =10$ ms. c Same as a but for a slow stimulus, $\tau _\mu =100$ ms. Other parameters as in Fig. 4

To characterize the error of the theoretical approximations more systematically, we compare theory and simulations as a function of the stimulus properties (Fig. 5). To this end, we model $\mu (t)$ as a complex stimulus sampled from a stationary Ornstein–Uhlenbeck process with correlation time $\tau _\mu$, mean ${{\bar{\mu }}}$ and variance $(1+\tau _\text {m}/\tau _\mu ){{\bar{\sigma }}}^2$. This parametrization has been chosen such that the non-resetting membrane potential ${{\hat{V}}}$ has mean ${{\bar{\mu }}}$ and standard deviation ${{\bar{\sigma }}}$ in the stationary state. For a given realization $\mu (t)$, we quantify the deviation of the theoretical ISI distribution $P_\mu (t|0)$ from the simulated one ${{\hat{P}}}_\mu (t|0)$ using the Kolmogorov–Smirnov (KS) statistics (Press et al. 1992). This statistics is then averaged over the stimulus ensemble (the subscript $\mu$ indicates the dependence on a given realization $\mu (t)$). Explicitly, the mean KS statistics is defined as

$$\begin{aligned} {\mathcal {D}}&=\left\langle \max _{t>0}\left| \int _0^t P_\mu (s|0)\,\hbox {d}s-\int _0^t {{\hat{P}}}_\mu (s|0)\,\hbox {d}s\right| \right\rangle _\mu \nonumber \\&=\left\langle \max _{t>0}\left| S_\mu (t|0)-{{\hat{S}}}_\mu (t|0)\right| \right\rangle _\mu , \end{aligned}$$

(49)

where $\langle \cdot \rangle _\mu$ denotes the ensemble average over realizations $\mu (t)$. Thus, the KS statistics can also be interpreted as the largest absolute difference between the survival function $S_\mu (t|0)$ and the simulated survival function ${{\hat{S}}}_\mu (t)$ (see Fig. 4, second panels from bottom).

The analysis confirms our previous observations that the decoupling approximations perform best in the subthreshold regime (${{\bar{\mu }}}<1$) at small stimulus variations ${{\bar{\sigma }}}$ (Fig. 5); they both become worse in the tonically firing regime (${{\bar{\mu }}}>1$). Although the qualitative dependence on the stimulus parameters is similar between the first- and second-order DA, the error is considerably smaller for the second-order DA throughout stimulus parameters. On the other hand, the Chizhov–Graham (C&G) theory has an opposite dependence; it generally performs well in the tonically firing regime (${{\bar{\mu }}}>1$) and shows small weaknesses in the subthreshold regime (Fig. 5b, ${{\bar{\mu }}}<1$), but it exhibits a good overall performance. For all three approximations, the error is larger for a rapidly changing stimulus (Fig. 5a). Interestingly, in the strongly mean-driven regime (${{\bar{\mu }}}>1$), a constant or weakly fluctuating stimulus (${{\bar{\sigma }}}\ll 1$) turns out to more difficult for the second-order DA than a more strongly fluctuating stimulus (Fig. 5b, c).

5 Population activity of LIF neurons (time-dependent firing rate)

5.1 Integral equation

As an application of the noise mapping, we consider the dynamics of the time-dependent firing rate, or equivalently the population activity of LIF neurons with colored input noise. Being in possession of an approximate hazard rate, it is straightforward to use the renewal integral equation (Gerstner 2000; Gerstner et al. 2014) (or equivalently, the refractory density equation (Chizhov and Graham 2007, 2008; Dumont et al. 2016; Schwalger and Chizhov 2019; Pietras et al. 2020)) to compute the population activity forward in time. To this end, let us consider a population of N uncoupled LIF neurons with colored input noise, Eq. (1). The spike train $X_i(t)$ of a given neuron i, $i=1,\dotsc ,N$ is defined as $X_i(t)=\sum _{k}\delta (t-t_{i,k})$, where $\{t_{i,k}\}_{k\in {\mathbb {Z}}}$ are the spike times of that neuron. The population activity is defined as the total number of spikes in a small time bin $(t,t+\varDelta t)$ divided by the number of neurons N and the time step $\varDelta t$. In the limit of infinitely many neurons and infinitesimally small time steps, we obtain the deterministic population activity

$$\begin{aligned} A(t)=\lim _{\varDelta t\rightarrow \infty }\lim _{N\rightarrow \infty }\frac{1}{N\varDelta t}\sum _{i=1}^N\int _t^{t+\varDelta t}X_i(t')\,\hbox {d}t'. \end{aligned}$$

(50)

Note that this expression can also be interpreted as an ensemble or trial average of a single neuron spike train, i.e., A(t) is equivalent to the time-dependent firing rate of a single neuron measured over many trials or realizations of a statistical ensemble. The evolution of the population activity is given by the renewal equation (Cox 1962; Gerstner et al. 2014)

$$\begin{aligned} A(t)=P\left( t|t_0\right) +\int _{t_0^+}^t P\left( t|{\hat{t}}\right) A\left( {\hat{t}}\right) \,\hbox {d}{\hat{t}}, \end{aligned}$$

(51)

where $P(t|{\hat{t}})$ is given by Eq. (43) and $t_0^+$ denotes the right-sided limit. In Eq. (51), we assumed that the population is initialized with a spike of each neuron at time $t_0$ (“synchronized initial condition”). The first term $P(t|t_0)$ represents the contribution from neurons that fire at time t for the first time after the initial spike at $t_0$. The integral equation (51) can be efficiently solved numerically (Gerstner and Kistler 2002). In particular, for numerical solutions, it is useful to turn the exponential factor into a differential equation as in Eq. (8):

$$\begin{aligned} \frac{\mathrm {d}S(t|{\hat{t}})}{\mathrm {d}t}=-\lambda (t|{\hat{t}})S(t|{\hat{t}}),\qquad S({\hat{t}}|{\hat{t}})=1 \end{aligned}$$

(52)

for all ${\hat{t}}<t$.

Fig. 6

Macroscopic population activity of non-adapting neurons under non-stationary driving. a Weak subthreshold stimulus $\mu (t)$ (i) as in Fig. 4a leads to a mean membrane potential response $u(t|t_0)$ below threshold at $V_\text {T}=1$ (ii). The resulting population activity A(t) is shown in (iii) and (iv) for strong ($\sigma _V=0.25$) and weak ($\sigma _V=0.1$) background noise, respectively. Gray circles: MC simulations of $10^6$ trials; red solid line: Chizhov–Graham theory, Eq. (7), (101); blue dashed line: level-crossing theory with independent upcrossings (first-order decoupling approximation), Eq. (45), (43); blue solid line: level-crossing theory with correlated upcrossings (second-order decoupling approximation), Eq. (46), (43). b The same for a suprathreshold stimulus as in Fig. 4b, for which the mean membrane potential u(t) reaches the threshold. In both panels, $\tau _\text {s}=4$ ms, $\tau _\text {m}=10$ ms, $t_\text {ref}=4$ ms and the population was initialized at time $t_0=-25$ ms (initial transient not shown)

5.2 Comparison with simulations and C&G theory

As an example, we studied the response of the population activity to the complex stimulus $\mu (t)$ shown in Fig. 6 Ai and Bi. In the subthreshold regime, where the membrane potential remains below threshold (Fig. 6a), the level-crossing theory well predicts the population activity obtained from simulations, while the C&G prediction exhibits small deviations as expected from the deviations of the ISI density in the subthreshold regime discussed above (Figs. 3 and 4). The agreement is good for both strong and weak noise.

For suprathreshold stimuli, where the membrane potential exceeds the threshold, the first-order decoupling approximation shows clear deviations (Fig. 6b). However, accounting for correlations between level-crossings in the second-order approximation recovers the population activity of simulated neurons for both strong and weak noise. Similarly, the C&G theory shows an excellent agreement with simulations.

6 Discussion

We developed a level-crossing theory for the hazard rate of a leaky integrate-and-fire neuron driven by colored input noise. To this end, we generalized the Stratonovich approximation for the first-passage-time (FPT) density (Stratonovich 1967b; Verechtchaguina et al. 2006; van Meegen and van Albada 2019) to the time-inhomogeneous case, where the stimulus or boundary is time-dependent, and derived a simplification that is local in time. Because higher-order correlations between upcrossings are approximated through their pair-wise correlations, we referred to this theory as the second-order decoupling approximation (DA). Besides the mean membrane potential u(t), the simplified hazard rate depends on the speed $\dot{u}$ and one additional variable z(t), which accounts for correlations between level crossings. Therefore, the escape-noise model defined by this hazard rate consists of only one extra first-order differential equation, Eq. (41), besides the dynamics of u, Eq. (2a). Our simulation results for the time-dependent interspike-interval (ISI) density and population activity show that the mapped LIF model with escape-noise well matches the LIF model with colored input noise. Thus, the hazard rate in the second-order DA (link function Eq. (42) and dynamics of z, Eq. (41)) provides a novel map from input noise to escape noise. We note that the dependence on the speed $\dot{u}$ is important and qualitatively differs from commonly used escape-noise models, where the link function only depends on the mean membrane potential u. Given the extensive theoretical literature on population models with simple link functions $\varPsi (u)$ (Gerstner 2000; Cormier et al. 2020; Schmutz et al. 2021), it will be an interesting question for further studies how the mean-field dynamics is influenced by an additional dependence on the membrane potential speed $\dot{u}$.

The map based on the second-order DA should be compared to the 1st-order DA, which neglects any correlations between upcrossings and represents a time-dependent generalization of the Hertz approximation (Verechtchaguina et al. 2006), and the previously proposed map by Chizhov and Graham (C&G) (Chizhov and Graham 2008). The generalized Hertz approximation (first-order DA) involves less ad hoc approximations compared to the second-order DA (cf. Eqs. (28) and (30)), and performs well in the fluctuation-driven (subthreshold) firing regime at low firing rates. On the other hand, its region of validity, Eq. (34), is rather limited, especially transiently large firing rates and mean-driven (suprathreshold) firing are not well described by the first-order approximation. Furthermore, the gain in numerical efficiency compared to the second-order DA is minor, e.g., simulating the firing rate trajectory of 200ms in Fig. 6b (middle) took 134ms for the first-order DA versus 165 ms for the second-order DA (Julia code run on an Intel(R) Core(TM) i7-8550U CPU @ 1.80 GHz).

On the contrary, the C&G map exhibits some weaknesses in the fluctuation-driven regime, while it has an excellent performance for short, mean-driven firing-rate transients. This behavior is expected because the theory represents an interpolation between two limit cases, where the theory is expected to work well: strong positive drift toward the threshold without diffusion (cf. also (Goedeke and Diesmann 2008)) and pure diffusion without drift. During short mean-driven periods the drift-induced firing dominates and diffusion effects can be safely neglected. An advantage of the C&G hazard rate, Eq. (103), is its simpler mathematical form and thus easier numerical implementation than the hazard rates based on the level-crossing theory (first- and second-order DA). Furthermore, the C&G theory permits to take the white-noise limit, $\tau _\text {s}\rightarrow 0$, whereas the level-crossing theory is not well defined in this limit: for $\tau _\text {s}\rightarrow 0$, the upcrossing rate $f_1$ diverges (Rice 1945; Stratonovich 1967b) (cf. Eq. (87)). Despite the divergence in the white-noise limit, we found in simulations that the 2nd-order DA performs well in the physiologically relevant range of synaptic time scales including synaptic time constants as small as $\tau _\text {s}=1$ ms (relative to $\tau _\text {m}=10$ ms, data not shown). On the other hand, the numerical efficiency of the C&G and the second-order DA are comparable (e.g., 175 ms and 165 ms run time, respectively, for the 200ms firing rate trajectory in Fig. 6b, middle). Overall, the C&G theory represents a good compromise between simplicity and accuracy.

Apart from the mapping of input noise to escape noise, the analysis performed in this paper also provided some analytical insights into the Stratonovich approximation. First, the ansatz of Stratonovich, Eq. (17), has been originally proposed for a system of “non-approaching” random points (level crossings) (Stratonovich 1967a, b). In our terminology, this means that the pair correlation function at zero time lag is $R(t,t)=-1$. Put differently, the conditional rate $\nu _{\text {cond}}(t,\tau )=f_2(t,t+\tau )/f_1(t)$ of an upcrossing to occur a time lag $\tau$ after a crossing at time t vanishes for $\tau \rightarrow 0$ if upcrossings are non-approaching. However, we found that in our case of the membrane potential driven by an exponentially correlated Gaussian noise, i.e., a doubly low-pass-filtered white noise (cf. Eq. (1) or (5)), the upcrossings do not form a system of non-approaching points. The conditional rate $\nu _\text {cond}$ at zero time lag has a non-vanishing minimum (corresponding to a reduced probability of close upcrossings, $\nu _\text {cond}<f_1$) and can even exceed the stationary upcrossing rate, $\nu _\text {cond}>f_1$, (the probability of an upcrossing is increased by an immediately preceding upcrossing, as already noted by (Burak et al. 2009) for stationary level-crossings). Given the excellent agreement of the second-order DA with simulations, the ansatz Eq. (17) seems to be more general and not limited to systems of non-approaching random points.

Based on the assumption of non-approaching level crossings, the threshold-crossing process has been frequently used as an analytically tractable model of neural spike generation. Examples include the calculation of information rates (DeWeese and Bialek 1995), pairwise correlations and synchronization of neurons due to shared inputs (Tchumatchenko et al. 2010b, a; Burak et al. 2009) and stochastic resonance (Jung 1995). The intuition behind this assumption is that level crossings exhibit refractoriness (Puelma Touzel and Wolf 2016) or a silence period (Tchumatchenko et al. 2010b) because it takes some time for a trajectory to re-cross the threshold from below. While this intuition is true for sufficiently smooth Gaussian processes (Tchumatchenko et al. 2010b, a) (auto-correlation function must be at least four times differentiable at 0), it fails if the velocity of the process is not differentiable (third derivative of auto-correlation at 0 does not exist), as in the present study and in (Verechtchaguina et al. 2006; Tchumatchenko and Wolf 2011; Badel 2011; Puelma Touzel and Wolf 2016). Because neurons exhibit some degree of refractoriness, the Gaussian processes of threshold-crossing neurons should be sufficiently smooth to be useful as a spiking neuron model.

By mapping input noise to escape noise we could apply the renewal integral equation to predict the time-dependent population activity of infinitely many LIF neurons with colored input noise. This detour via an approximate escape-noise model allowed us to circumvent the direct numerical solution of the two-dimensional Fokker–Planck equation associated with the LIF model Eq. (1). An intriguing question is whether the same indirect approach could be used to solve the important problem of finitely many neurons with input noise. Neural circuits in the brain are often modeled by networks of integrate-and-fire neurons driven by Poissonian input noise (e.g., (Diesmann et al. 1999; Potjans and Diesmann 2014; Donoso et al. 2018)). In these network models, the number of neurons per cell type range from about hundred to a few thousand cells, consistent with experimental estimations (Lefort et al. 2009). On this mesoscopic scale, finite-size fluctuations of the population activity cannot be neglected. It is, however, unknown how to generalize the Fokker–Planck equation to a stochastic population equation in the case of finite neuron numbers, so as to account for finite-size fluctuations. On the other hand, the problem of finite-size neural population equations has been recently solved for LIF neurons with escape noise in the form of a stochastic integral equation (Schwalger et al. 2017; Schmutz et al. 2021). In the original paper (Schwalger et al. 2017), we have applied the stochastic integral equation to the cortical microcircuit model of (Potjans and Diesmann 2014) by roughly fitting an escape-noise model with the simple link function $\varPsi (u)=ce^{\beta u}$ to match mean population activities of simulation data. However, with the map derived in this paper, where $\varPsi$ depends on u and $\dot{u}$, it should be possible to directly use the stochastic integral equation as a new mesoscopic population model for finite-size populations of LIF neurons driven by colored input noise.

Appendices

Appendix

FPT density from level-crossing statistics

1.1 General expression for survivor function

The sequence of upward crossings of the freely evolving, non-resetting membrane potential across the threshold, or shortly the set of “upcrossings,” forms a point process $\{t_1,t_2,\dotsc \}$ in time with $t_i>0$. Thus, the upcrossing times are defined by ${{{\hat{V}}}}(t_i)=V_\text {T}$ and ${\dot{{\hat{V}}}}(t_i)>0$. As any point process, the upcrossing times for $t>0$ can be fully characterized by the joint distribution functions $f_1(t_1)$, $f_2(t_1,t_2)$, $f_3(t_1,t_2,t_3)$, ... (see, e.g., (Stratonovich 1967a; van Kampen 1992)). These functions are defined such that

$$\begin{aligned} f_k(t_1,\dotsc ,t_k)\mathrm {d}t_1\cdots \mathrm {d}t_k+{\mathcal {O}}(\mathrm {d}t) \end{aligned}$$

(53)

is the probability to find an upcrossing in each of the non-overlapping intervals $[t_1,t_1+\mathrm {d}t_1)$, ..., $[t_k,t_k+\mathrm {d}t_k)$, with sufficiently small intervals $\mathrm {d}t_1,\dotsc ,\mathrm {d}t_k<\mathrm {d}t$ and non-coinciding arguments $t_i\ne t_j$ for all $i\ne j$. In the case of coinciding arguments $t_i=t_j$ for some $i\ne j$, the function $f_k$ is understood to be its limit value for $t_i\rightarrow t_j$.

For our purpose, it will be more convenient to use the correlation functions $g_1(t_1)$, $g_2(t_1,t_2)$, $g_3(t_1,t_2,t_3)$, ... (see, e.g., (Stratonovich 1967a; van Kampen 1992)). Similar to the joint distribution functions $\{f_k\}$, the system of correlation functions $\{g_k\}$ completely characterizes the statistics of the upcrossing times. To define the correlation functions, we first introduce the generating functional for the $f_k$ given by

$$\begin{aligned} L[v]\equiv \left\langle \prod _{t_i>0}\left( 1+v(t_i) \right) \right\rangle , \end{aligned}$$

(54)

where v(t) is a test function (Stratonovich 1967a; van Kampen 1992). It can be shown that expanding the generating functional in powers of v(t) yields

$$\begin{aligned} L[v]= & {} 1+\sum _{k=1}^\infty \frac{1}{k!}\int _{0}^\infty \!\!\cdots \!\!\int _{0}^\infty f_k(t_1,\dotsc ,t_k)\nonumber \\&\times v(t_1)\cdots v(t_k)\,\mathrm {d}t_1\cdots \mathrm {d}t_k, \end{aligned}$$

(55)

i.e., the functions $f_k$ are the expansion coefficients of the generating functional. Therefore, the joint distribution functions $f_k$ can be uniquely generated by functional differentiation of L[v]. In analogy to the cumulants of a random variable that are generated from the logarithm of the moment generating function, the correlation functions $g_k$ can be obtained from $\ln L$ as follows:

$$\begin{aligned} g_k(t_1,\dotsc ,t_k)=\left. \frac{\delta ^k\ln L[v(t)]}{\delta v(t_1)\cdots \delta v(t_k)}\right| _{v(t)\equiv 0} \end{aligned}$$

(56)

In particular, the first two correlation functions read

$$\begin{aligned} g_1(t)&=f_1(t), \end{aligned}$$

(57)

$$\begin{aligned} g_2(t_1,t_2)&=f_2(t_1,t_2)-f_1(t_1)f_1(t_2). \end{aligned}$$

(58)

By means of the correlation functions, the survivor function S(t), i.e., the probability for having no upcrossing in the interval [0, t), can be expressed as Eq. (12).

1.2 Moments and correlation functions of the Gaussian process

In contrast to the vanishing first moments $\langle x\rangle =\langle y\rangle =0$ and the stationary variance $\sigma _y^2=\langle y^2(t)\rangle$, the second moments $\sigma _x^2(t)=\langle x^2(t)\rangle$ and $\sigma _{xy}(t)=\langle x(t)y(t)\rangle$ are time-dependent. They obey the differential equation (Risken 1984)

$$\begin{aligned} \frac{\mathrm {d}(\sigma _x^2)}{\mathrm {d}t}&=-2\left( \gamma \sigma _x^2-\sigma _{xy}\right) , \end{aligned}$$

(59)

$$\begin{aligned} \frac{\mathrm {d}\sigma _{xy}}{\mathrm {d}t}&=-{\tilde{\tau }}^{-1}\sigma _{xy}+\sigma _y^2. \end{aligned}$$

(60)

with $\sigma _y^2=D/\tau _y$, ${\tilde{\tau }}^{-1}=\gamma +\tau _y^{-1}$ and $\sigma _x^2(0)=\sigma _{xy}(0)=0$. The explicit solution is

$$\begin{aligned} \sigma _{xy}(t)&={\tilde{\tau }}\sigma _y^2\left( 1-e^{-t/{\tilde{\tau }}}\right) , \end{aligned}$$

(61a)

$$\begin{aligned} \sigma _x^2(t)&=\frac{{\tilde{\tau }}\sigma _y^2}{\gamma }\left( 1-e^{-2\gamma t}\right) \nonumber \\&\quad +\frac{2{\tilde{\tau }}\sigma _y^2}{2\gamma -{\tilde{\tau }}^{-1}}\left( e^{-2\gamma t}-e^{-\frac{t}{{\tilde{\tau }}}}\right) \end{aligned}$$

(61b)

For large t, the process [x(t), y(t)] becomes stationary with the following constant moments

$$\begin{aligned} \sigma _x^2=\frac{1}{\gamma }\sigma _{xy}=\frac{{\tilde{\tau }}}{\gamma }\sigma _y^2. \end{aligned}$$

(62)

1.3 Joint distribution functions for upcrossings

Let us denote the point process of the upcrossings by $\{{\hat{t}}_i\}_{i=1,2,\dotsc }$. Using the properties of the Dirac-delta function with respect to a change of variable, the corresponding spike train can be written as

$$\begin{aligned} s(t)&=\sum _{i=1}^\infty \delta (t-{\hat{t}}_i),\nonumber \\&=\left[ \dot{x}(t)-\dot{b}(t) \right] \delta \bigl (x(t)-b(t)\bigr )\theta \bigl (\dot{x}(t)-\dot{b}(t)\bigr ). \end{aligned}$$

(63)

here the Heaviside function $\theta \bigl (\dot{x}(t)-\dot{b}(t)\bigr )$ accounts for the fact that the velocity of the process x(t) must be larger than the velocity of the boundary b(t) at an upcrossing, i.e., it ensures that x(t) crosses the moving boundary b(t) from below. Note that Eq. (63) can be seen as an extension of the Kac–Rice formula (Azaïs and Wschebor 2009) to moving boundaries. The joint distribution function is defined as

$$\begin{aligned} f_k(t_1,\dotsc ,t_k)=\left\langle s(t_1)\cdots s(t_k) \right\rangle \end{aligned}$$

(64)

(for $t_i\ne t_j$ for $i,j=1,\dotsc ,k$, $i\ne j$). Substituting Eq. (63) into Eq. (64) and taking the average yields

$$\begin{aligned}&f_k(t_1,\dotsc ,t_k)=\int _{\dot{b}_1}^\infty \!\!\cdots \!\!\int _{\dot{b}_k}^\infty \mathrm {d}\dot{x}_1\cdots \mathrm {d}\dot{x}_k\,\nonumber \\&\quad \times \left( \dot{x}_1-\dot{b}_1 \right) \cdots \left( \dot{x}_k-\dot{b}_k \right) p_{2k}^{(x,\dot{x})}(b_1,\dotsc ,b_k,\dot{x}_1,\dotsc ,\dot{x}_k),\nonumber \\ \end{aligned}$$

(65)

where $b_i$ and $\dot{b}_i$ is short-hand for $b(t_i)$ and $\dot{b}(t_i)$, respectively. Furthermore, $p_{2k}^{(x,\dot{x})}(x_1,\dotsc ,x_k,\dot{x}_1,\dotsc ,\dot{x}_k)$ is the joint probability density for the variables $x_i=x(t_i)$ and $\dot{x}_i=\dot{x}(t_i)$. In our case of the two-dimensional Ornstein–Uhlenbeck process, Eq. (5), $p_{2k}^{(x,\dot{x})}$ can be simply expressed by the joint probability density $p_{2k}(x_1,\dotsc ,x_k,y_1,\dotsc ,y_k)$ of the variables $x_i=x(t_i)$ and $y_i=y(t_i)$:

$$\begin{aligned}&p_{2k}^{(x,\dot{x})}(b_1,\dotsc ,b_k,\dot{x}_1,\dotsc ,\dot{x}_k)\nonumber \\&\quad s=p_{2k}(b_1,\dotsc ,b_k,\gamma b_1+\dot{x}_1,\dotsc ,\gamma b_k+\dot{x}_k). \end{aligned}$$

(66)

Inserting this expression into Eq. (65) yields

$$\begin{aligned}&f_k(t_1,\dotsc ,t_k)=\int _{0}^\infty \!\!\cdots \!\!\int _{0}^\infty \mathrm {d}w_1\cdots \mathrm {d}w_k\,w_1\cdots w_k\nonumber \\&\quad \times p_{2k}(b_1,\dotsc ,b_k,\gamma b_1+\dot{b}_1+w_1,\dotsc ,\gamma b_k+\dot{b}_k+w_k),\nonumber \\ \end{aligned}$$

(67)

where we made the substitution $\dot{x}_i=\dot{b}_i+w_i$ with new integration variables $w_i$. We note, however, that for higher-dimensional models, it is generally more convenient to directly compute the density $p_{2k}^{(x,\dot{x})}$ and use Eq. (65). For example, for a $(n+1)$-dimensional Gaussian process $\mathbf {x}(t)=[x(t),y_1(t),\dotsc ,y_n(t)]^T$, this density is determined by the time-dependent correlation functions $\langle x(t)x(t+\tau )\rangle$, $\langle x(t)\dot{x}(t+\tau )\rangle$, $\langle \dot{x}(t) x(t+\tau )\rangle$ and $\langle \dot{x}(t)\dot{x}(t+\tau )\rangle$, which can be obtained from the time-dependent covariance matrix of $\mathbf {x}(t)$ in a straightforward manner.

1.4 Upcrossing rate $f_1(t)$

Using the moments $\sigma _x^2(t)$, $\sigma _{xy}(t)$ and $\sigma _y^2$ derived in Sect. 2, the joint probability density of x and y is given by the bivariate Gaussian distribution

$$\begin{aligned} p_2(x,y,t)=\frac{1}{2\pi \sqrt{|C_2|}}\exp \left( -\frac{\sigma _y^2x^2-2\sigma _{xy}xy+\sigma _x^2y^2}{2|C_2|}\right) \nonumber \\ \end{aligned}$$

(68)

with $|C_2|=\sigma _x^2\sigma _y^2-\sigma _{xy}^2$. This allows us to calculate the upcrossing rate $f_1(t)$ from Eq. (67). The integration can be performed analytically resulting in the formula Eq. (14).

1.5 Correlations between upcrossings for small time lag

Here, we are interested in the probability that two upcrossings occur very close to each other. More precisely, we want to calculate the probability density $f_2(t,t+\tau )$ in the limit when the distance $\tau$ between upcrossings goes to zero.

1.5.1 Time-dependent boundary

To this end, we need the probability density of the four-dimensional vector $z=[x(t),x(t+\tau ),y(t),y(t+\tau )]^T$, which is given by the multivariate Gaussian distribution

$$\begin{aligned} p_4(z)=\frac{1}{4\pi ^2\sqrt{|C_4|}}\exp \left( -\frac{1}{2}z^TC_4^{-1}z \right) . \end{aligned}$$

(69)

This distribution is determined by the correlation matrix $C_4$ with elements $\left( C_4 \right) _{ij}=\left\langle z_iz_j \right\rangle$:

$$\begin{aligned} C_4= \begin{pmatrix} \sigma _x^2(t)&{}\quad C_{xx}(t,\tau )&{}\quad \sigma _{xy}(t)&{}\quad C_{xy}(t,\tau )\\ C_{xx}(t,\tau )&{}\quad \sigma _x^2(t+\tau )&{}\quad C_{yx}(t,\tau )&{}\quad \sigma _{xy}(t+\tau )\\ \sigma _{xy}(t)&{}\quad C_{yx}(t,\tau )&{}\quad \sigma _{y}^2&{}\quad C_{yy}(\tau )\\ C_{xy}(t,\tau )&{}\quad \sigma _{xy}(t+\tau )&{}\quad C_{yy}(\tau )&{}\quad \sigma _{y}^2 \end{pmatrix},\nonumber \\ \end{aligned}$$

(70)

where we used the notations $C_{xx}(t,\tau )\equiv \left\langle x(t)x(t+\tau ) \right\rangle$, $C_{xy}(t,\tau )\equiv \left\langle x(t)y(t+\tau ) \right\rangle$ and $C_{yx}(t,\tau )\equiv \left\langle y(t)x(t+\tau ) \right\rangle$. Furthermore, note that $\sigma _y^2=\left\langle y^2(t) \right\rangle$ and $C_{yy}(\tau )\equiv \left\langle y(t)y(t+\tau ) \right\rangle$ do not depend on time because of the stationarity of y(t). The correlation functions for $\tau >0$ can be computed from the regression theorem (Risken 1984):

$$\begin{aligned} C_{xx}(t,\tau )&=G_{xx}(\tau )\sigma _x^2(t)+G_{xy}(\tau )\sigma _{xy}(t), \end{aligned}$$

(71)

$$\begin{aligned} C_{xy}(t,\tau )&=G_{yx}(\tau )\sigma _x^2(t)+G_{yy}(\tau )\sigma _{xy}(t), \end{aligned}$$

(72)

$$\begin{aligned} C_{yx}(t,\tau )&=G_{xx}(\tau )\sigma _{xy}(t)+G_{xy}(\tau )\sigma _{y}^2(t), \end{aligned}$$

(73)

$$\begin{aligned} C_{yy}(\tau )&=G_{yy}(\tau )\sigma _y^2(t), \end{aligned}$$

(74)

where we used the elements of the Green’s function

$$\begin{aligned} G(\tau )= \begin{pmatrix} G_{xx}(\tau )&{}G_{xy}(\tau )\\ 0&{}G_{yy}(\tau ) \end{pmatrix} \end{aligned}$$

(75)

of the Ornstein–Uhlenbeck process Eq. (5). Using the negative drift matrix of the Ornstein–Uhlenbeck process $A=\bigl ( {\begin{matrix} \gamma &{}-1\\ 0&{}1/\tau _y \end{matrix}}\bigr )$, the Green’s function is obtained from $G(\tau )=e^{-A\tau }$:

$$\begin{aligned} G_{xx}(\tau )=e^{-\gamma \tau },\qquad G_{yy}(\tau )=e^{-\tau /\tau _y}, \end{aligned}$$

(76)

$$\begin{aligned} G_{xy}(\tau )=\frac{\tau _y}{1-\gamma \tau _y}\left( e^{-\gamma \tau }-e^{-\tau /\tau _y} \right) . \end{aligned}$$

(77)

In Eq. (70), we also need the time-shifted moments $\sigma _x^2(t+\tau )$ and $\sigma _{xy}(t+\tau )$. These can be obtained from $\sigma _x^2(t)$ and $\sigma _{xy}(t)$ by propagating Eq. (59). This yields

$$\begin{aligned} \sigma _{xy}(t+\tau )&=e^{-\frac{\tau }{{\tilde{\tau }}}}\sigma _{xy}+{\tilde{\tau }}\left( 1-e^{-\frac{\tau }{{\tilde{\tau }}}} \right) \sigma _{y}^2 \end{aligned}$$

(78)

$$\begin{aligned} \sigma _x^2(t+\tau )&=e^{-\frac{2\tau }{\tau _\text {m}}}\sigma _{x}^2+\frac{\sigma _y^2{\tilde{\tau }}}{\gamma }\bigl (1-e^{-2\gamma \tau }\bigr )\nonumber \\&\quad +\frac{2{\tilde{\tau }}(\sigma _{xy}-\sigma _y^2{\tilde{\tau }})}{1-2\gamma {\tilde{\tau }}}\bigl (e^{-2\gamma \tau }-e^{-\tau /{\tilde{\tau }}}\bigr ). \end{aligned}$$

(79)

Because we are interested in the limit $\tau \rightarrow 0$, we can expand the moving threshold at time t to linear order such that

$$\begin{aligned} b(t+\tau )=b(t)+\dot{b}(t)\tau +{\mathcal {O}}(\tau ^2). \end{aligned}$$

(80)

The two-point joint density follows from Eq. (67) and (69):

$$\begin{aligned} f_2(t,t+\tau )= & {} \frac{1}{4\pi ^2\sqrt{|C_4|}}\iint _{0}^\infty \mathrm {d}w_1\mathrm {d}w_2\,w_1w_2\nonumber \\&\times \exp \left( -\frac{1}{2}z^TC_4^{-1}z \right) \end{aligned}$$

(81)

with $z=[b,b+\dot{b} \tau ,\gamma b+\dot{b}+w_1,\gamma (b+\dot{b}\tau )+\dot{b}+w_2]^T$. A straightforward but lengthy series expansion of the exponent ${{\tilde{B}}}=-\frac{1}{2}z^TC_4^{-1}z$ for small $\tau$ yields

$$\begin{aligned} \tilde{B}=-\frac{w_1^2+w_1w_2+w_2^2}{(\sigma _y^2/\tau _y)\tau }-\tilde{B}_0(b,\dot{b},w_1,w_2)+{\mathcal {O}}(\tau ),\nonumber \\ \end{aligned}$$

(82)

where

$$\begin{aligned}&{\tilde{B}}_0(b,\dot{b},w_1,w_2)=\frac{1}{4\sigma _y^2|C_2|}\Big \{\left( \gamma ^2\sigma _x^2-2\gamma \sigma _{xy}+\sigma _y^2\right) b^2 \nonumber \\&\quad +2[\gamma \left( \sigma _{xy}^2+\sigma _x^2\sigma _y^2\right) -2\sigma _{xy}\sigma _y^2]b(\dot{b}+w_1)\nonumber \\&\quad -2\gamma \left( \sigma _{xy}^2-\sigma _x^2\sigma _y^2\right) b\left( \dot{b}+w_2\right) \nonumber \\&\quad +\left[ \left( 1-\gamma \tau _y\right) \sigma _x^2\sigma _y^2+\left( 1+\gamma \tau _y\right) \sigma _{xy}^2\right] (\dot{b}+w_1)^2\nonumber \\&\quad -(1+\gamma \tau _y)(\sigma _{xy}^2-\sigma _x^2\sigma _y^2)(\dot{b}+w_2)^2\Big \} \end{aligned}$$

(83)

is an ${\mathcal {O}}(1)$ term. The first term of ${\tilde{B}}$ is of order $1/\tau$ and has a maximum at the lower integration boundary $w_1=w_2=0$. Therefore, the neighborhood of the point $w_1=w_2=0$ dominates the integral in the limit $\tau \rightarrow 0$. At this point the term ${{\tilde{B}}}_0(b,\dot{b},w_1,w_2)$ coincides with $B(b,\dot{b})$ in Eq. (15). Thus, we can write

$$\begin{aligned} f_2(t,t+\tau )\sim \frac{1}{4\pi ^2\sqrt{|C_4|}}e^{-B(b,\dot{b})}I(\tau ),\qquad \tau \rightarrow 0 \end{aligned}$$

(84)

with the Gaussian integral

$$\begin{aligned} I(\tau )&=\iint _{0}^\infty \mathrm {d}w_1\mathrm {d}w_2\,w_1w_2\exp \left( -\frac{w_1^2+w_1w_2+w_2^2}{(\sigma _y^2/\tau _y)\tau } \right) \nonumber \\&=\frac{9-\sqrt{3}\pi }{27}\frac{\sigma _y^4}{\tau _y^2}\tau ^2. \end{aligned}$$

(85)

As a last step, we expand the determinant $|C_4|$ to lowest order in $\tau$:

$$\begin{aligned} |C_4|=\frac{\left( \sigma _x^2\sigma _y^2-\sigma _{xy}^2\right) \sigma _y^4}{3\tau _y^2}\tau ^4+{\mathcal {O}}(\tau ^5). \end{aligned}$$

Combining all factors yields the two-point upcrossing density in the limit of zero lag given by Eq. (23).

In the stationary case, $\dot{b}=0$ and $t\rightarrow \infty$, the formula for $f_2(t,t)$ reduces to

$$\begin{aligned} f_2^{(s)}(t,t)=\frac{\left( 3\sqrt{3}-\pi \right) \tau _y}{18\pi \left( 1+\gamma \tau _y\right) }f_1 \end{aligned}$$

(86)

with the stationary upcrossing rate

$$\begin{aligned} f_1=\frac{1}{2\pi }\sqrt{\frac{\gamma }{\tau _y}}\exp \left( -\frac{b^2}{2\sigma _x^2} \right) . \end{aligned}$$

(87)

This expression results in the pair correlation function Eq. (25).

1.5.2 Auto-correlation function of up-crossings for stationary, differentiable Gaussian processes

In the stationary case, the rate of upcrossings $f_1$ is constant and the second order distribution function as well as the auto-correlation function of x only depend on the time difference, hence $f_2(t,t+\tau )=f_2(\tau )$ and $C_{xx}(t,t+\tau )=C_{xx}(\tau )$. A classical result for the upcrossing rate is (Rice 1945)

$$\begin{aligned} f_1=\frac{\sqrt{|C_{xx}''(0)|/\sigma _x^2}}{2\pi }\exp \left( -\frac{b^2}{2\sigma _x^2} \right) . \end{aligned}$$

(88)

Here, we derive the asymptotic behavior of $f_2(\tau )$ for small time lag $\tau$. To this end, we expand $C_{xx}(\tau )$

$$\begin{aligned} C_{xx}(\tau )=c_0+\sum _{k=2}^\infty \frac{c_k}{k!}|\tau |^k. \end{aligned}$$

(89)

where $c_k=C_{xx}^{(k)}(0)$ denotes the k-th right-sided derivative of the correlation function at zero time lag. Here, we have taken into account that the auto-correlation function is an even function. Furthermore, we have not included the first-order term $c_1|\tau |$ because the derivative $C_{xx}'(0)=C_{x\dot{x}}(0)$ must be zero for differentiable processes x(t), i.e., for velocities $\dot{x}$ with finite variance. For example, the one-dimensional Ornstein–Uhlenbeck process is excluded because it exhibits a kink in its auto-correlation function $C_{xx}(\tau )\sim e^{-|\tau |/\tau _{cor}}$ at zero lag (i.e., $c_1<0$) implying an infinite variance of the velocity, $\sigma _{\dot{x}}^2=-c_2=\infty$, and hence a diverging up-crossing rate, Eq. (88). This divergence arises for any one-dimensional Langevin dynamics, for which the velocity $\dot{x}$ exhibits a white noise component, and reflects the fractal nature of Markovian diffusion processes (Jung 1994). In the following, we distinguish three cases, all of which have occurred in previous studies: (i) $c_3\ne 0$ corresponding to a kink in the velocity correlation function $C_{\dot{x} \dot{x}}(\tau )=C_{xx}''(\tau )$. This case is considered in the present work as well as in previous models (Verechtchaguina et al. 2006; Tchumatchenko and Wolf 2011; Badel 2011; Puelma Touzel and Wolf 2016). (ii) $c_3=0$ and $c_5\ne 0$ corresponding to a kink in the correlation function of the acceleration $\ddot{x}(t)$, as in (Verechtchaguina et al. 2006). And (iii), $c_3=0$ and $c_5=0$ which occurs, e.g., for smooth correlation functions as used in (Tchumatchenko et al. 2010b, a).

In the first case, $c_3\ne 0$, i.e., when $C_{xx}''(\tau )$ has a kink at zero lag and thus the acceleration $\ddot{x}$ has infinite variance as in our model Eq. (1), we find in lowest order in $\tau$

$$\begin{aligned} f_2(\tau )\sim \frac{3\sqrt{3}-\pi }{18\pi }f_1\left| \frac{c_3}{c_2}\right| ,\qquad \tau \rightarrow 0. \end{aligned}$$

(90)

This expression recovers a previous result obtained in (Burak et al. 2009). Furthermore, the case $c_3=0$, $c_5\ne 0$, yields the following lowest-order behavior

$$\begin{aligned} f_2(\tau )\sim \frac{c_5^{3/2}\exp \left( -\frac{b^2}{2\left( c_0-c_2^2/c_4\right) } \right) }{90\sqrt{15}\pi ^2\sqrt{|c_2|^3+c_0c_2c_4}}\tau ^{\frac{5}{2}},\qquad \tau \rightarrow 0. \end{aligned}$$

(91)

To the best of our knowledge, this expression is a novel result. Finally, the third case $c_3=0$ and $c_5=0$, yields in lowest order

$$\begin{aligned} f_2(\tau )\sim \frac{\left| c_4^2-c_2c_6\right| \exp \left( -\frac{b^2}{2\left( c_0-c_2^2/c_4\right) } \right) }{1296\pi ^2\sqrt{c_2^2-c_0c_4}}\tau ^4,\qquad \tau \rightarrow 0,\nonumber \\ \end{aligned}$$

(92)

which has been reported before (Burak et al. 2009). In the derivation of Eqs. (90)–(92), we have used the Gaussian integral

$$\begin{aligned} \iint _0^\infty \hbox {d}\dot{x}_1\hbox {d}\dot{x}_2\,\dot{x}_1\dot{x}_2\exp \left( -\frac{(\dot{x}_1+\dot{x}_2)^2}{\beta } \right) =\frac{\beta ^2}{12}. \end{aligned}$$

(93)

Chizhov–Graham theory

An elegant approximation for the FPT problem has been put forward by Chizhov and Graham (2007, 2008), which we will state here without proof. The idea is to construct the hazard function from two limit cases: First, for an excitatory current that is much faster than the diffusion time, the probability flux across the threshold is dominated by the deterministic positive drift, whereas the noise can be treated as frozen. For a monotonic movement of the mean membrane potential toward the threshold ($\dot{u}(t)>0$), one can simply shift the Gaussian probability density along its time-dependent center and calculate the survival probability as the total probability mass that is below the threshold at time t:

$$\begin{aligned} S_\text {drift}(t)&=\int _{-\infty }^{b(t)}\frac{\mathrm {d}x}{\sqrt{2\pi }\sigma _x(t)}\exp \left( -\frac{x^2(t)}{2\sigma _x^2(t)} \right) , \end{aligned}$$

(94)

$$\begin{aligned}&=\frac{1}{2}\left[ 1+\text {erf}\left( \frac{b(t)}{\sqrt{2}\sigma _x(t)} \right) \right] \end{aligned}$$

(95)

In contrast, for negative movement of the center of mass, i.e., downward and away from the threshold, the survival probability is kept constant. The hazard rate corresponding to the deterministic drift is given by $-\frac{\hbox {d}}{\hbox {d}t}\ln (S_\text {drift})$ resulting in

$$\begin{aligned} \varPhi _\text {drift}(b,\dot{b},t)=\frac{2}{\sqrt{\pi }}\left[ -\dot{T}\right] _+\frac{\exp (-T^2)}{\text {erfc}(T)}. \end{aligned}$$

(96)

Following (Chizhov and Graham 2007, 2008), we introduced the dimensionless quantity

$$\begin{aligned} T(t)=\frac{b(t)}{\sqrt{2}\sigma _x(t)}, \end{aligned}$$

(97)

the temporal derivative of which is given by

$$\begin{aligned} \dot{T}(t)=\frac{1}{\sqrt{2}\sigma _x}\left( \dot{b}+\frac{b\left( \sigma _x^2-\tau _\text {m}\sigma _{xy}\right) }{\tau _\text {m}\sigma _x^2} \right) . \end{aligned}$$

(98)

The moments $\sigma _x^2(t)$ and $\sigma _{xy}(t)$ have been derived above (Sect. 2) and $\tau _\text {m}=1/\gamma$. Note that the second term in Eq. (98) accounts for the non-stationarity of the variance $\sigma _x^2(t)$. This term is absent in the original formula in (Chizhov and Graham 2007, 2008), which assumed stationary fluctuations with $\sigma _x(t)=\text {const}.$. This version with stationary fluctuations has also been derived in (Goedeke and Diesmann 2008).

Second, the effect of diffusion can be captured in the quasi-stationary limit case of slow driving. In this case, the survival probability can be calculated analytically, resulting in the corresponding hazard rate

$$\begin{aligned} \varPhi _\text {diff}(b)&=\varPhi _\text {diff}^\text {wn}(T)\left[ 1-\left( 1+\frac{\tau _\text {m}}{\tau _\text {s}}\right) ^{-0.71+0.0825(T+3)}\right] , \end{aligned}$$

(99)

$$\begin{aligned} \varPhi _\text {diff}^\text {wn}(T)&=\tau _\text {m}^{-1}\exp \bigl (6.1\cdot 10^{-3}-1.12T-0.25T^2\nonumber \\&\quad -0.072T^3-0.0117T^4\bigr ). \end{aligned}$$

(100)

here the numerical coefficients have been fitted to the exact solution (Chizhov and Graham 2008). The total hazard rate is simply given by the sum of the two limit cases:

$$\begin{aligned} \lambda (t)=\varPhi _\text {drift}(b,\dot{b},t)+\varPhi _\text {diff}(b,t). \end{aligned}$$

(101)

Thus, we obtain for the hazard rate of the LIF neuron with absolute refractory period $t_\text {ref}$ and given last spike time ${\hat{t}}$

$$\begin{aligned} \lambda (t|{\hat{t}})\approx \varPsi _\text {CG}\bigl ( u(t|{\hat{t}}),\dot{u}(t|{\hat{t}}),t-{\hat{t}} \bigr ) \end{aligned}$$

(102)

with the Chizhov–Graham link function

$$\begin{aligned} \varPsi _\text {CG}\bigl ( u,\dot{u},\tau \bigr )= & {} \theta (\tau -t_\text {ref}) \bigl [\varPhi _\text {drift}(V_\text {T}-u,-\dot{u},\tau -t_\text {ref})\nonumber \\&+\varPhi _\text {diff}(V_\text {T}-u,\tau -t_\text {ref})\bigr ]. \end{aligned}$$

(103)