∎

11institutetext: T. Schwalger 22institutetext: Institute of Mathematics,
Technical University Berlin, 10623 Berlin, Germany
22email<EMAIL_ADDRESS>33institutetext: T. Schwalger
44institutetext: Bernstein Center for Computational Neuroscience Berlin, 10115
Berlin, Germany

# Mapping Input Noise to Escape Noise in Integrate-and-fire neurons: A Level-
Crossing Approach

Tilo Schwalger

(Received: date / Accepted: date)

###### 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.

###### Keywords:

Integrate-and-fire neuron Interspike interval density First-passage-time
density Colored noise Escape noise Hazard rate Threshold-crossing statistics
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 GerKis14 . 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 Kni72 ; Ger00 ; AugLad17 ; SchDeg17 .

The population activity of noisy spiking neurons can be mathematically
described by population density equations NykTra00 ; Chi17 . 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 GerKis14 . 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_ Ger00 . The corresponding
population density equation is a Fokker-Planck equation and the population
activity can be obtained as the probability flux across the threshold AbbVre93
; BruHak99 ; FouBru02 ; NykTra00 ; Ric08 ; AugLad17 . 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 Bru00 ; PotDie14 . The fluctuating part of this
external shot noise leads, via a diffusion approximation GerKis14 , 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 FouBru02 ; SchDro15  and spatial LinDoi05
correlations of synaptic inputs and enables mean-field theories for recurrent
networks of sparsely-connected integrate-and-fire neurons BruHak99 ; Bru00 .

In the second model class (called _output noise_ or _escape noise_ Ger00 ),
the dynamical equations for the state variables are deterministic while spikes
(“output”) are generated stochastically through a hazard rate or conditional
intensity Ger00 ; Pan04 ; TruEde05 ; PilShl08 ; PilLat08 ; NauGer12 ; BreSen13
; GalLoe16 ; GerDeg17 ; RaaDit20 . This hazard rate depends on the state
variables via a link function. For example, it may be given as
$\hat{\lambda}(t)=\Psi(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 NauGer12
process, the corresponding population density equation is a renewal integral
equation or, equivalently, a refractory density equation GerKis14 ; Ger00 ;
NauGer12 ; ChiGra07 ; ChiGra08 ; DumHen16 ; SchChi19 . Although output noise
is of phenomenological nature without a quantitative link to biophysical
mechanisms, it has several advantages SchChi19  owing to its simpler
mathematical tractability: First, the refractory density or integral equation
admits an extension to finite numbers of neurons SchDeg17 ; SchChi19 ;
SchGer20 ; SchLoe21_arxiv . 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 Pan04 ; TruEde05 ; PilShl08 ; GerDeg17 ;
MenNau12 ; PozMen15 ; TeeIye18 . And third, the state space for models with
output noise remains approximately one-dimensional even for multi-dimensional
conductance-based neuron models ChiGra07 . The one-dimensional description
permits highly efficient numerical solutions, in contrast to Fokker-Planck
equations ApfLy06 , 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 $\Psi$, 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
nonhomogeneous 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 nonhomogeneous renewal
models in neuroscience are one-dimensional integrate-and-fire neurons driven
by white input noise GerKis14 . 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}}^{t}P(s|\hat{t})\,ds\bigr{]}$
GerKis14 . 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^{\prime})\\}_{t^{\prime}\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 BulEls96 ; SchTal04 ; Lin04b . 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 $\Psi\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^{\prime}),z(t^{\prime})\\}_{t^{\prime}\in(\hat{t},t)}$).

Several theoretical studies have suggested approximate local hazard rates for
leaky integrate-and fire (LIF) models driven by white PleGer2000 ; HerGer01 ;
ChiGra07  or exponentially-correlated ChiGra08  Gaussian noise, or quasi-
static (frozen) noise GoeDie08 . 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 VSS2006b . In Sec. 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 Sec. 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 (Sec. 3.4). Then, we turn to the LIF
model and the problem of mapping input noise to escape noise (Secs. 4) and
apply this map to predict the time-dependent population activity of LIF
neurons with colored input noise (Sec. 5). Each of the sections 3, 4 and 5
closes with a comparison of the level-crossing theory with simulations and a
previous theory by Chizhov and Graham ChiGra08 . Detailed derivations are
provided in the 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 SchLSG08 ;
GerKis14 ; SchDie15 . 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

$\displaystyle\tau_{\text{m}}\dot{V}$ $\displaystyle=-V+\mu(t)+\eta(t),$ (1a)
$\displaystyle\tau_{\text{s}}\dot{\eta}$
$\displaystyle=-\eta+\sqrt{2\tau_{\text{s}}}\sigma_{\eta}\xi(t),$ (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^{\prime})\rangle=\delta(t-t^{\prime})$. 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 Ger00
given by a hazard rate (conditional intensity) of the form
$\Psi\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:

$\tau_{\text{m}}\dot{u}=-u+\mu(t).$ (2a) Furthermore, we allow an explicit
dependence on the speed of the membrane potential $\dot{u}(t)$ (in accordance
with previous studies PleGer2000 ; HerGer01 ; ChiGra07 ; GoeDie08 ), 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 $\text{Pr}(\text{spike in
$(t,t+dt)$}|u(t),\dot{u}(t),\\{z_{i}(t)\\},t-\hat{t})\\\
=\Psi\bigl{(}u(t),\dot{u}(t),{\color[rgb]{0,0,0}\\{z_{i}(t)\\}},t-\hat{t}\bigr{)}dt$
(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) Lin04 ; SchLSG08 ; SchDro15 .
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
$\Psi$ 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 GerKis14 .

To derive the link function $\Psi$ that maps input to output noise, one needs
to solve a first-passage-time (FPT) problem: As mentioned in the 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 towards positive values of $\eta$
Lin04 ; SchLSG08 ; Sch13 ; SchDro15 , 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 towards 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
${\color[rgb]{0,0,0}\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.

Figure 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_{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{\color[rgb]{0,0,0}\hat{V}}(t)\rangle=u(t)$:

$\displaystyle x(t)$ $\displaystyle={\color[rgb]{0,0,0}\hat{V}}(t)-u(t)$ (3)
$\displaystyle b(t)$ $\displaystyle=V_{T}-u(t),$ (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

$\displaystyle\dot{x}$ $\displaystyle=-\gamma x+y$ (5a)
$\displaystyle\tau_{y}\dot{y}$ $\displaystyle=-y+\sqrt{2D}\xi(t)$ (5b)

with initial conditions

$x(0)=0,\qquad y(0)\sim\mathcal{N}(0,\sigma_{y}^{2})$ (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)dt=\text{Prob}(t^{*}\in[t,t+dt))$ 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 are
reaction times of bimolecular chemical reactions HanTal90 . 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 Cox62 .  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 $dt$ to find a boundary crossing in the interval
$(t,t+dt)$ 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+dt)=S(t)(1-\lambda(t)dt)$, hence $dS(t)/dt=-\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)\,ds\right)$ for $t>0$. The probability
to find the first crossing  after time $0$ in the interval $(t,t+dt)$ is equal
to the probability to find a crossing in $(t,t+dt)$ and to have no crossings
in $(0,t)$. Hence, the FPT density is given by the product
$P(t)=\lambda(t)S(t)$, or

$P(t)=\lambda(t)\exp\left(-\int_{0}^{t}\lambda(s)\,ds\right).$ (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,

$P(t)=\lambda(t)S(t),\qquad\frac{\mathrm{d}S}{\mathrm{d}t}=-\lambda(t)S(t),\qquad
S(0)=1.$ (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^{\prime}\in(0,t)$. By contrast, $\lambda(t)$ is not a
probability density and can thus, in principle, be arbitrary as long as it is
non-negative and $S(t)=\exp\left(-\int_{0}^{t}\lambda(s)\,ds\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

$\lambda(t)\approx\Phi\bigl{(}b(t),\dot{b}(t),\dotsc,\\{z_{i}(t)\\},t\bigr{)},$
(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 $\Phi$ for the boundary-dependent hazard rate of the model Eq.
(5) and $\Psi$ for the voltage-dependent hazard rate of the model Eq. (2b).
The two functions are related in a simple way, see Sec. 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 Ric45 ; RicSat83 ; VSS2006b ;
BraThu17 ; AzaWsc09 . 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^{*}+\Delta 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^{*}]$

$s(t)=\sum_{i=1}^{N(t^{*})}\delta(t-\hat{t}_{i}),$ (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})\dotsb s(t_{k})\rangle$, for
$k=1,2,\dotsc$ and non-coinciding time arguments $t_{i}$ Str67I ; van92 . The
moment functions can be interpreted such that for a small time step $\Delta t$
the quantity $f_{k}(t_{1},\dotsc,t_{k})\Delta t^{k}$ yields the probability to
find an up-crossing events in each of the non-overlapping intervals
$(t_{1},t_{1}+\Delta t)$, …, $(t_{k},t_{k}+\Delta 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,
both for stationary and non-stationary processes (see appendix, Sec. A.3).

The distribution functions $f_{k}$ allow for an exact series expression of the
FPT density, sometimes called Wiener-Rice series VSS2006b ; BraThu17 :

$P(t)=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!}\int_{0}^{t}dt_{1}\dotsb
dt_{k}\,f_{k+1}(t_{1},\dotsc,t_{k},t)$ (11)

A detailed explanation of this formula is given in reference VSS2006b . In
brief, it counts – for a large ensemble of trajectories – the number of those
trajectories that have a crossing in $[t,t+dt)$ but no crossing in $(0,t)$.
Starting with the fraction $f_{1}(t)dt$ 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}^{t}f_{2}(t_{1},t)dt_{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}^{t}f_{3}(t_{1},t_{2},t)dt_{1}dt_{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 Str67I ;
van92  and Sec. A.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 Str67I ; van92

$S(t)=\exp\left(\sum_{k=1}^{\infty}\frac{(-1)^{k}}{k!}\int_{0}^{t}\mathrm{d}t_{1}\dotsb\mathrm{d}t_{k}\,g_{k}(t_{1},\dotsc,t_{k})\right).$
(12)

From this expression, the Wiener-Rice series for the FPT density, Eq. (11) is
recovered by $P(t)=-dS(t)/dt$. Similarly, the hazard rate can be obtained by
$\lambda(t)=-d(\ln S)/dt$. 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 Str67I : 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 Sec. 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\geq 2$ Str67I . In this case, Eq. (12) reduces to
$S(t)=\exp\left(-\int_{0}^{t}f_{1}(\tau)\,d\tau\right)$, and hence the FPT
density reads

$P(t)\approx f_{1}(t)\exp\left\\{-\int_{0}^{t}f_{1}(\tau)\,d\tau\right\\}.$
(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 A.3 and A.4). The result is  the first-order
decoupling approximation:

$\lambda(t)\approx f_{1}(t)=\Phi_{1}\bigl{(}b(t),\dot{b}(t),t\bigr{)}\\\
:=\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)},$
(14)

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

$\displaystyle B(b,\dot{b},t)$
$\displaystyle=\frac{(\gamma^{2}\sigma_{x}^{2}-2\gamma\sigma_{xy}+\sigma_{y}^{2})b^{2}+2(\gamma\sigma_{x}^{2}-\sigma_{xy})b\dot{b}+\sigma_{x}^{2}\dot{b}^{2}}{2(\sigma_{x}^{2}\sigma_{y}^{2}-\sigma_{xy}^{2})}.$
(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

$\displaystyle\sigma_{xy}(t)$
$\displaystyle=\tilde{\tau}\sigma_{y}^{2}\left(1-e^{-t/\tilde{\tau}}\right),$
(16a) $\displaystyle\sigma_{x}^{2}(t)$
$\displaystyle=\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)$ (16b)

with $\sigma_{y}^{2}=D/\tau_{y}$ and $\tilde{\tau}^{-1}=\gamma+\tau_{y}^{-1}$
(see Sec. A.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 RicSat83 ; Bad11 .

####  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_{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\geq 3$, proposed by Stratonovich Str67I ;
Str67II . 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\geq 2$ are
approximated by the ansatz  Str67I ; Str67II

$g_{k}(t_{1},\dotsc,t_{k})=(k-1)!f_{1}(t_{1})\dotsb
f_{1}(t_{k})\\{R(t_{2},t_{1})\dotsb R(t_{k},t_{1})\\}_{sym}.$ (17)

Here, the function $R(t,t^{\prime})$ describes the pairwise interactions
between events at time $t$ and $t^{\prime}$, and $\\{\dotsb\\}_{sym}$ denotes
the operation of symmetrization (i.e. the arithmetic mean of all permutations
of the time arguments).  As suggested in Str67I ; Str67II , we choose
$R(t,t^{\prime})$ as the normalized auto-correlation function

$R(t,t^{\prime})=\frac{f_{2}(t,t^{\prime})}{f_{1}(t)f_{1}(t^{\prime})}-1,$
(18)

which makes the ansatz Eq. (17) exact for $k=2$. Note that compared to Str67I
; Str67II , we use an opposite sign in the definition of $R$ for mathematical
convenience. The auto-correlation function $R(t,t^{\prime})$ can be
interpreted as the conditional probability density of an event at time
$t^{\prime}$ given an event at time $t$ normalized by the unconditional
probability density $f_{1}(t^{\prime})$ and shifted by the mean such that
$R(t,t^{\prime})=0$ if events at time $t$ and $t^{\prime}$ are independent.
For stationary point processes, $R(t,t^{\prime})=R(|t-t^{\prime}|)$ only
depends on the time difference. In analogy to the common use for spatial point
processes, $R(t-t^{\prime})$ 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^{\prime}|\gg\tau_{\text{cor}}$, they occur
independently, hence $f_{2}(t,t^{\prime})\approx f_{1}(t)f_{1}(t^{\prime})$.
This implies a vanishing auto-correlation function $R(t,t^{\prime})\approx 0$.
Secondly, the behavior when $t$ and $t^{\prime}$ 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^{\prime})$ vanishes. In
contrast, a positive pair correlation function $R(t,t^{\prime})>0$ at small
time lag indicates that events are attractive and tend to cluster. Conversely,
for a negative pair correlation function $R(t,t^{\prime})<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^{\prime})=0$ and hence
$R(t,t^{\prime})=-1$ if $t$ and $t^{\prime}$ fall within a refractory period.
Similarly, non-approaching random points Str67I  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$ VSS2006b ; PueWol16 . In Sec. 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\geq 3$, which in general is not expected to be exact.
Nevertheless, the ansatz and the above-described behavior of $R(t,t^{\prime})$
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^{\prime})\equiv 0$ and thus Eq. (17) recovers the expected
result $g_{k}\equiv 0$ for all $k\geq 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)$
Str67II , 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 Str67I ; Str67II ; VSS2006b ; MeeAlb21_arxiv

$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\\},$
(19)

where

$\displaystyle q(t,\tau)$
$\displaystyle=\int_{0}^{t}R(\tau,\tau^{\prime})f_{1}(\tau^{\prime})\,d\tau^{\prime}$
(20)
$\displaystyle=\frac{1}{f_{1}(\tau)}\int_{0}^{t}\left[f_{2}(\tau,\tau^{\prime})-f_{1}(\tau)f_{1}(\tau^{\prime})\right]\,d\tau^{\prime}$

is a measure of upcrossing correlations on the time scale $t$. The formula Eq.
(19) has been termed Stratonovich approximation VSS2006b . 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 MeeAlb21_arxiv  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

$q(t,\tau)>-1$ (21)

for all $t$ and $\tau$ so as to keep the argument of the logarithm positive
VSS2006b .

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^{\prime})$ involves a further integration (taking already into
account that one of the two integrals in the definition of $f_{2}$, Eq. (81),
Sec. A.5, can be evaluated analytically RicSat83 ; VSS2006b ; we note that
$f_{2}$ can also be expressed in terms of Owen’s T function MeeAlb21_arxiv ).
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

Figure 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}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\tau_{\text{m}}=\gamma^{-1}$, $\tau\neq 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 shows 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.

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. Ric45 ; Str67II ; VSS2006b ; TchMal10 ), 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 the Appendix Sec. A.5.2, we
also provide formulas for the auto-correlation function of general Gaussian
level-crossing processes in the stationary state (see also BurLew09  for
special cases and Jun94  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

$R_{0}(t)=\frac{f_{2}(t,t)}{f_{1}^{2}(t)}-1,$ (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, Sec. A.5). The
result is

$\displaystyle f_{2}(t,t)$
$\displaystyle=\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)}$
(23) $\displaystyle=:\hat{f}_{2}\bigl{(}b(t),\dot{b}(t),t\bigr{)}$ (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

$R_{0}=\beta\frac{1+\gamma\tau_{y}}{\sqrt{\gamma\tau_{y}}}\exp\left(\frac{b^{2}}{2\sigma_{x}^{2}}\right)-1$
(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

$\frac{b_{crit}}{\sigma_{x}}=\sqrt{2\ln\left(\frac{\sqrt{\gamma\tau_{y}}}{\beta(1+\gamma\tau_{y})}\right)},$
(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 Eq. (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
VSS2006b ; PueWol16 .

On the other hand, for high barriers such that $|b|>b_{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

$\lambda(t)=f_{1}(t)\left\\{F\bigl{(}q(t,t)\bigr{)}+\int_{0}^{t}f_{1}(\tau)F^{\prime}\bigl{(}q(t,\tau)\bigr{)}R(t,\tau)\,\mathrm{d}\tau\right\\}$
(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^{\prime}(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

$\lambda(t)=\frac{f_{1}(t)}{1+q(t)},$ (28)

where we used the short-hand notation

$q(t)\equiv q(t,t)=\int_{0}^{t}R(t,t^{\prime})f_{1}(t^{\prime})\,dt^{\prime}.$
(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_{corr}$ around its peak at the integration boundary
$\tau=t$, where $\tau_{corr}$ is the correlation time defined in Eq. (31)
below (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
MeeAlb21_arxiv , 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

$R(t,t^{\prime})\approx
R_{0}(t)\exp\left(-\frac{|t-t^{\prime}|}{\tau_{corr}}\right),$ (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:

$\tau_{corr}=\int_{0}^{\infty}\frac{C_{xx}(\tau)}{C_{xx}(0)}\,\mathrm{d}\tau=\tau_{\text{m}}+\tau_{\text{s}}.$
(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_{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:

$q(t)\approx R_{0}(t)z(t),$ (32)

where
$z(t)=\int_{0}^{t}\exp\left[-(t-t^{\prime})/\tau_{corr}\right]f_{1}(t^{\prime})$
defines a new auxiliary variable that satisfies the differential equation

$\frac{\mathrm{d}z}{\mathrm{d}t}=-\frac{1}{\tau_{corr}}z+f_{1}(t)$ (33)

with $z(0)=0$. We note that the slightly different ansatz
$R(t,t^{\prime})\approx
R_{0}(t^{\prime})\exp\left(-\frac{t-t^{\prime}}{\tau_{corr}}\right)$ yields
slightly different equations with similarly good results. In Sec. 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

$|q(t,t)|\ll 1$ (34)

for all $t>0$.

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

$\lambda(t)\approx\Phi_{2}\bigl{(}b(t),\dot{b}(t),z(t),t\bigr{)}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\frac{\Phi_{1}\bigl{(}b(t),\dot{b}(t),t\bigr{)}}{1+\hat{R}_{0}\bigl{(}b(t),\dot{b}(t),t\bigr{)}z(t)}.$
(35)

Here, $\Phi_{1}$ is given by Eq. (14) and

$\hat{R}_{0}\bigl{(}b,\dot{b},t\bigr{)}=\frac{\hat{f}_{2}\bigl{(}b,\dot{b},t\bigr{)}}{\left[\Phi_{1}\bigl{(}b,\dot{b},t\bigr{)}\right]^{2}}-1$
(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 $\Phi_{1}$, the hazard rate $\Phi_{2}$ depends on the additional
local variable $z$ that obeys

$\frac{\mathrm{d}z}{\mathrm{d}t}=-\frac{1}{\tau_{corr}}z+\Phi_{1}\bigl{(}b(t),\dot{b}(t),t\bigr{)},\quad
z(0)=0.$ (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

Figure 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) $\pm$ 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 ChiGra08 . An extended variant of the Chizhov-Graham (C&G) theory is
presented in Appendix B, Eq. (101).

For concreteness, we consider a periodically moving boundary:

$b(t)=1+\alpha\cos(2\pi ft)$ (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
$\Phi$ for the FPT with moving boundary $b(t)$, it is easy to formulate the
link function $\Psi$ 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

$\dot{b}(t)=-\dot{u}(t)=\frac{u(t)-\mu(t)}{\tau_{\text{m}}}$ (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

$\Psi_{1}\bigl{(}u,\dot{u},\tau\bigr{)}=\theta(\tau-
t_{\text{ref}})\Phi_{1}(V_{\text{T}}-u,-\dot{u},\tau-t_{\text{ref}}).$ (40)

Here, $\theta(t)=\mathbb{1}_{t\geq 0}$ is the Heaviside step function and
$\Phi_{1}$ is given by Eq. (14). Note that in the first-order DA, the link
function $\Psi(u,\dot{u},z,\tau)=\Psi_{1}(u,\dot{u},\tau)$ does not depend on
an auxiliary variable $z$. In contrast, the 2nd-order DA exhibits an
additional auxiliary variable $z$. Taking the last spike time and the absolute
refractory period into account, its dynamics reads

$\dot{z}=-\frac{z}{\tau_{\text{m}}+\tau_{\text{s}}}+\Psi_{1}(u,-\dot{u},t-\hat{t})$
(41)

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

$\Psi_{2}\bigl{(}u,\dot{u},z,\tau\bigr{)}=\theta(\tau-
t_{\text{ref}})\Phi_{2}(V_{\text{T}}-u,-\dot{u},z,\tau-t_{\text{ref}}),$ (42)

where $\Phi_{2}$ is given by Eq. (35).

###  4.2 Comparison with simulation and C&G theory

Figure 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)/dt$ 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$.

To judge the performance of the level-crossing theory given by the link
functions $\Psi_{1}$ and $\Psi_{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

$\displaystyle P(t|\hat{t})$ $\displaystyle=\lambda(t|\hat{t})S(t|\hat{t}),$
(43) $\displaystyle S(t|\hat{t})$
$\displaystyle=\exp\left(-\int_{\hat{t}}^{t}\lambda(s|\hat{t})\,ds\right),$
(44)

where for the first-order decoupling approximation(DA)

$\lambda(t|\hat{t})\approx\Psi_{1}\bigl{(}u(t|\hat{t}),\dot{u}(t|\hat{t}),t-\hat{t}\bigr{)},$
(45)

and for the second-order DA,

$\displaystyle\lambda(t|\hat{t})$
$\displaystyle\approx\Psi_{2}\bigl{(}u(t|\hat{t}),\dot{u}(t|\hat{t}),z(t|\hat{t}),t-\hat{t}\bigr{)}$
(46)

with $\Psi_{1}$ and $\Psi_{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

$\displaystyle\dot{u}(t|\hat{t})$
$\displaystyle=-\frac{u(t|\hat{t})+\mu(t)}{\tau_{\text{m}}},$ (47)
$\displaystyle\dot{z}(t|\hat{t})$
$\displaystyle=-\frac{z(t|\hat{t})}{\tau_{\text{m}}+\tau_{\text{s}}}+\Psi_{1}\Bigl{(}u(t|\hat{t}),-\dot{u}(t|\hat{t}),t-\hat{t}\Bigr{)}.$
(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}}^{t}P(s|\hat{t})\,ds$. 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 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_{0}z]$ is larger than unity (Note that $q\equiv R_{0}z>-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).

Figure 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 PreTeu92 . 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

$\displaystyle\mathcal{D}$
$\displaystyle=\left\langle\max_{t>0}\left|\int_{0}^{t}P_{\mu}(s|0)\,ds-\int_{0}^{t}\hat{P}_{\mu}(s|0)\,ds\right|\right\rangle_{\mu}$
$\displaystyle=\left\langle\max_{t>0}\left|S_{\mu}(t|0)-\hat{S}_{\mu}(t|0)\right|\right\rangle_{\mu},$
(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 1st- and 2nd-
order DA, the error is considerably smaller for the 2nd-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 2nd-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 Ger00 ; GerKis14
(or equivalently, the refractory density equation ChiGra07 ; ChiGra08 ;
DumHen16 ; SchChi19 ; PieGal20 ) 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+\Delta t)$ divided by the number of neurons $N$ and the time step
$\Delta t$. In the limit of infinitely many neurons and infinitesimally small
time steps, we obtain the deterministic population activity

$A(t)=\lim_{\Delta t\to\infty}\lim_{N\rightarrow\infty}\frac{1}{N\Delta
t}\sum_{i=1}^{N}\int_{t}^{t+\Delta t}X_{i}(t^{\prime})\,dt^{\prime}.$ (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 Cox62 ; GerKis14

$A(t)=P(t|t_{0})+\int_{t_{0}^{+}}^{t}P(t|\hat{t})A(\hat{t})\,d\hat{t},$ (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 GerKis02 . In particular,
for numerical solutions, it is useful to turn the exponential factor into a
differential equation as in Eq. (8):

$\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$ (52)

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

Figure 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 (1st-order decoupling
approximation), Eq. (45), (43); blue solid line: level-crossing theory with
correlated upcrossings (2nd-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. 6Ai 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 (Fig. 3 and Fig. 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 Str67II ;
VSS2006b ; MeeAlb21_arxiv  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 2nd-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
$\Psi(u)$ Ger00 ; CorTan20 ; SchLoe21_arxiv , 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 2nd-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 VSS2006b , and the
previously proposed map by Chizhov and Graham (C&G) ChiGra08 . The generalized
Hertz approximation (1st-order DA) involves less ad-hoc approximations
compared to the 2nd-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 2nd-order DA is minor: e.g.,
simulating the firing rate trajectory of 200ms in Fig. 6B (middle) took 134ms
for the 1st-order DA versus 165ms for the 2nd-order DA  (Julia code run on an
Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz).

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 towards the threshold without
diffusion (cf. also GoeDie08 ) 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 (1st- and 2nd-order DA).
Furthermore, the C&G theory permits to take the white-noise limit,
$\tau_{\text{s}}\to 0$, whereas the level-crossing theory is not well defined
in this limit: for $\tau_{\text{s}}\to 0$, the upcrossing rate $f_{1}$
diverges Ric45 ; Str67II  (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 2nd-order DA are comparable (e.g. 175ms and
165ms 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) Str67I ; Str67II . 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\to 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 BurLew09  for
stationary level-crossings). Given the excellent agreement of the 2nd-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 DeWBia , pairwise correlations and synchronization of neurons due to
shared inputs TchMal10 ; TchGei10 ; BurLew09  and stochastic resonance Jun95b
. The intuition behind this assumption is that level crossings exhibit
refractoriness PueWol16  or a silence period TchMal10  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 TchMal10 ;
TchGei10  (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 VSS2006b ; TchWol11 ; Bad11 ; PueWol16 . 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. DieGew99 ;
PotDie14 ; DonSch18 ). In these network models, the number of neurons per cell
type range from about hundred to a few thousand cells, consistent with
experimental estimations LefTom09 . 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 SchDeg17 ; SchLoe21_arxiv
. In the original paper SchDeg17 , we have applied the stochastic integral
equation to the cortical microcircuit model of PotDie14  by roughly fitting an
escape-noise model with the simple link function $\Psi(u)=ce^{\beta u}$ to
match mean population activities of simulation data. However, with the map
derived in this paper, where $\Psi$ 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.

##  Appendix A FPT density from level-crossing statistics

###  A.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
${\color[rgb]{0,0,0}\hat{V}}(t_{i})=V_{\text{T}}$ and
${\color[rgb]{0,0,0}\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. Str67I ; van92 ). These functions are
defined such that

$f_{k}(t_{1},\dotsc,t_{k})\mathrm{d}t_{1}\dotsb\mathrm{d}t_{k}+\mathcal{O}(\mathrm{d}t)$
(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}\neq t_{j}$ for all $i\neq j$. In the case of coinciding
arguments $t_{i}=t_{j}$ for some $i\neq 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.
Str67I ; van92 ). 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

$L[v]\equiv\left\langle\prod_{t_{i}>0}\left(1+v(t_{i})\right)\right\rangle,$
(54)

where $v(t)$ is a test function Str67I ; van92 . It can be shown that
expanding the generating functional in powers of $v(t)$ yields

$L[v]=1+\sum_{k=1}^{\infty}\frac{1}{k!}\int_{0}^{\infty}\\!\\!\dotsb\\!\\!\int_{0}^{\infty}f_{k}(t_{1},\dotsc,t_{k})\\\
\times v(t_{1})\dotsb v(t_{k})\,\mathrm{d}t_{1}\dotsb\mathrm{d}t_{k},$ (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:

$g_{k}(t_{1},\dotsc,t_{k})=\left.\frac{\delta^{k}\ln L[v(t)]}{\delta
v(t_{1})\dotsb\delta v(t_{k})}\right|_{v(t)\equiv 0}$ (56)

In particular, the first two correlation functions read

$\displaystyle g_{1}(t)$ $\displaystyle=f_{1}(t),$ (57) $\displaystyle
g_{2}(t_{1},t_{2})$
$\displaystyle=f_{2}(t_{1},t_{2})-f_{1}(t_{1})f_{1}(t_{2}).$ (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).

###  A.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 Ris84

$\displaystyle\frac{\mathrm{d}(\sigma_{x}^{2})}{\mathrm{d}t}$
$\displaystyle=-2\left(\gamma\sigma_{x}^{2}-\sigma_{xy}\right),$ (59)
$\displaystyle\frac{\mathrm{d}\sigma_{xy}}{\mathrm{d}t}$
$\displaystyle=-\tilde{\tau}^{-1}\sigma_{xy}+\sigma_{y}^{2}.$ (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

$\displaystyle\sigma_{xy}(t)$
$\displaystyle=\tilde{\tau}\sigma_{y}^{2}\left(1-e^{-t/\tilde{\tau}}\right),$
(61a) $\displaystyle\sigma_{x}^{2}(t)$
$\displaystyle=\frac{\tilde{\tau}\sigma_{y}^{2}}{\gamma}\left(1-e^{-2\gamma
t}\right)$
$\displaystyle\quad+\frac{2\tilde{\tau}\sigma_{y}^{2}}{2\gamma-\tilde{\tau}^{-1}}\left(e^{-2\gamma
t}-e^{-\frac{t}{\tilde{\tau}}}\right)$ (61b)

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

$\sigma_{x}^{2}=\frac{1}{\gamma}\sigma_{xy}=\frac{\tilde{\tau}}{\gamma}\sigma_{y}^{2}.$
(62)

###  A.3 Joint distribution functions for upcrossings

Let us denote the point process of the upcrossings by
$\\{\hat{t}_{i}\\}_{i=1,2,\dotsc}$. The corresponding spike train can be
written as

$\displaystyle s(t)$ $\displaystyle=\sum_{i=1}^{\infty}\delta(t-\hat{t}_{i}),$
$\displaystyle=\left[\dot{x}(t)-\dot{b}(t)\right]\delta\bigl{(}x(t)-b(t)\bigr{)}\theta\bigl{(}\dot{x}(t)-\dot{b}(t)\bigr{)}.$
(63)

Note that this equation can be seen as an extension of the Kac-Rice formula
AzaWsc09  to moving boundaries. The joint distribution function is defined as

$f_{k}(t_{1},\dotsc,t_{k})=\left\langle s(t_{1})\dotsb s(t_{k})\right\rangle$
(64)

(for $t_{i}\neq t_{j}$ for $i,j=1,\dotsc,k$, $i\neq j$). Substituting Eq. (63)
into Eq. (64) and taking the average yields

$f_{k}(t_{1},\dotsc,t_{k})=\int_{\dot{b}_{1}}^{\infty}\\!\\!\dotsb\\!\\!\int_{\dot{b}_{k}}^{\infty}\mathrm{d}\dot{x}_{1}\dotsb\mathrm{d}\dot{x}_{k}\,\\\
\times\left(\dot{x}_{1}-\dot{b}_{1}\right)\dotsb\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}),$
(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})$:

$p_{2k}^{(x,\dot{x})}(b_{1},\dotsc,b_{k},\dot{x}_{1},\dotsc,\dot{x}_{k})=p_{2k}(b_{1},\dotsc,b_{k},\gamma
b_{1}+\dot{x}_{1},\dotsc,\gamma b_{k}+\dot{x}_{k}).$ (66)

Inserting this expression into Eq. (65) yields

$f_{k}(t_{1},\dotsc,t_{k})=\int_{0}^{\infty}\\!\\!\dotsb\\!\\!\int_{0}^{\infty}\mathrm{d}w_{1}\dotsb\mathrm{d}w_{k}\,w_{1}\dotsb
w_{k}\\\ \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}),$ (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.

###  A.4 Uprossing rate $f_{1}(t)$

Using the moments $\sigma_{x}^{2}(t)$, $\sigma_{xy}(t)$ and $\sigma_{y}^{2}$
derived in Sec. A.2, the joint probability density of $x$ and $y$ is given by
the bivariate Gaussian distribution

$p_{2}(x,y,t)=\frac{1}{2\pi\sqrt{|C_{2}|}}\exp\left(-\frac{\sigma_{y}^{2}x^{2}-2\sigma_{xy}xy+\sigma_{x}^{2}y^{2}}{2|C_{2}|}\right)$
(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).

###  A.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.

####  A.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

$p_{4}(z)=\frac{1}{4\pi^{2}\sqrt{|C_{4}|}}\exp\left(-\frac{1}{2}z^{T}C_{4}^{-1}z\right).$
(69)

This distribution is determined by the correlation matrix $C_{4}$ with
elements $\left(C_{4}\right)_{ij}=\left\langle z_{i}z_{j}\right\rangle$:

$C_{4}=\begin{pmatrix}\sigma_{x}^{2}(t)&C_{xx}(t,\tau)&\sigma_{xy}(t)&C_{xy}(t,\tau)\\\
C_{xx}(t,\tau)&\sigma_{x}^{2}(t+\tau)&C_{yx}(t,\tau)&\sigma_{xy}(t+\tau)\\\
\sigma_{xy}(t)&C_{yx}(t,\tau)&\sigma_{y}^{2}&C_{yy}(\tau)\\\
C_{xy}(t,\tau)&\sigma_{xy}(t+\tau)&C_{yy}(\tau)&\sigma_{y}^{2}\end{pmatrix},$
(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 Ris84 :

$\displaystyle C_{xx}(t,\tau)$
$\displaystyle=G_{xx}(\tau)\sigma_{x}^{2}(t)+G_{xy}(\tau)\sigma_{xy}(t),$ (71)
$\displaystyle C_{xy}(t,\tau)$
$\displaystyle=G_{yx}(\tau)\sigma_{x}^{2}(t)+G_{yy}(\tau)\sigma_{xy}(t),$ (72)
$\displaystyle C_{yx}(t,\tau)$
$\displaystyle=G_{xx}(\tau)\sigma_{xy}(t)+G_{xy}(\tau)\sigma_{y}^{2}(t),$ (73)
$\displaystyle C_{yy}(\tau)$ $\displaystyle=G_{yy}(\tau)\sigma_{y}^{2}(t),$
(74)

where we used the elements of the Green’s function

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

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

$\displaystyle G_{xx}(\tau)=e^{-\gamma\tau},\qquad
G_{yy}(\tau)=e^{-\tau/\tau_{y}},$ (76) $\displaystyle
G_{xy}(\tau)=\frac{\tau_{y}}{1-\gamma\tau_{y}}\left(e^{-\gamma\tau}-e^{-\tau/\tau_{y}}\right).$
(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

$\displaystyle\sigma_{xy}(t+\tau)$
$\displaystyle=e^{-\frac{\tau}{\tilde{\tau}}}\sigma_{xy}+\tilde{\tau}\left(1-e^{-\frac{\tau}{\tilde{\tau}}}\right)\sigma_{y}^{2}$
(78) $\displaystyle\sigma_{x}^{2}(t+\tau)$
$\displaystyle=e^{-\frac{2\tau}{\tau_{\text{m}}}}\sigma_{x}^{2}+\frac{\sigma_{y}^{2}\tilde{\tau}}{\gamma}\bigl{(}1-e^{-2\gamma\tau}\bigr{)}$
$\displaystyle\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{)}.$
(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

$\displaystyle b(t+\tau)=b(t)+\dot{b}(t)\tau+\mathcal{O}(\tau^{2}).$ (80)

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

$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_{1}w_{2}\exp\left(-\frac{1}{2}z^{T}C_{4}^{-1}z\right)$
(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^{T}C_{4}^{-1}z$ for small $\tau$ yields

$\tilde{B}=-\frac{w_{1}^{2}+w_{1}w_{2}+w_{2}^{2}}{(\sigma_{y}^{2}/\tau_{y})\tau}-\tilde{B}_{0}(b,\dot{b},w_{1},w_{2})+\mathcal{O}(\tau),$
(82)

where

$\tilde{B}_{0}(b,\dot{b},w_{1},w_{2})=\frac{1}{4\sigma_{y}^{2}|C_{2}|}\left\\{(\gamma^{2}\sigma_{x}^{2}-2\gamma\sigma_{xy}+\sigma_{y}^{2})b^{2}\right.\\\
+2[\gamma(\sigma_{xy}^{2}+\sigma_{x}^{2}\sigma_{y}^{2})-2\sigma_{xy}\sigma_{y}^{2}]b(\dot{b}+w_{1})\\\
-2\gamma(\sigma_{xy}^{2}-\sigma_{x}^{2}\sigma_{y}^{2})b(\dot{b}+w_{2})\\\
+[(1-\gamma\tau_{y})\sigma_{x}^{2}\sigma_{y}^{2}+(1+\gamma\tau_{y})\sigma_{xy}^{2}](\dot{b}+w_{1})^{2}\\\
\left.-(1+\gamma\tau_{y})(\sigma_{xy}^{2}-\sigma_{x}^{2}\sigma_{y}^{2})(\dot{b}+w_{2})^{2}\right\\}$
(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

$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$ (84)

with the Gaussian integral

$\displaystyle I(\tau)$
$\displaystyle=\iint_{0}^{\infty}\mathrm{d}w_{1}\mathrm{d}w_{2}\,w_{1}w_{2}\exp\left(-\frac{w_{1}^{2}+w_{1}w_{2}+w_{2}^{2}}{(\sigma_{y}^{2}/\tau_{y})\tau}\right)$
$\displaystyle=\frac{9-\sqrt{3}\pi}{27}\frac{\sigma_{y}^{4}}{\tau_{y}^{2}}\tau^{2}.$
(85)

As a last step, we expand the determinant $|C_{4}|$ to lowest order in $\tau$:

$|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}).$

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

$f_{2}^{(s)}(t,t)=\frac{(3\sqrt{3}-\pi)\tau_{y}}{18\pi(1+\gamma\tau_{y})}f_{1}$
(86)

with the stationary upcrossing rate

$f_{1}=\frac{1}{2\pi}\sqrt{\frac{\gamma}{\tau_{y}}}\exp\left(-\frac{b^{2}}{2\sigma_{x}^{2}}\right).$
(87)

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

####  A.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 Ric45

$f_{1}=\frac{\sqrt{|C_{xx}^{\prime\prime}(0)|/\sigma_{x}^{2}}}{2\pi}\exp\left(-\frac{b^{2}}{2\sigma_{x}^{2}}\right).$
(88)

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

$C_{xx}(\tau)=c_{0}+\sum_{k=2}^{\infty}\frac{c_{k}}{k!}|\tau|^{k}.$ (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}^{\prime}(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 Jun94 . In the following, we distinguish three cases, all of which
have occurred in previous studies: (i) $c_{3}\neq 0$ corresponding to a kink
in the velocity correlation function
$C_{\dot{x}\dot{x}}(\tau)=C_{xx}^{\prime\prime}(\tau)$. This case is
considered in the present work as well as in previous models VSS2006b ;
TchWol11 ; Bad11 ; PueWol16 . (ii) $c_{3}=0$ and $c_{5}\neq 0$ corresponding
to a kink in the correlation function of the acceleration $\ddot{x}(t)$, as in
VSS2006b . And (iii), $c_{3}=0$ and $c_{5}=0$ which occurs, e.g., for smooth
correlation functions as used in TchMal10 ; TchGei10 .

In the first case, $c_{3}\neq 0$, i.e. when $C_{xx}^{\prime\prime}(\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$

$f_{2}(\tau)\sim\frac{3\sqrt{3}-\pi}{18\pi}f_{1}\left|\frac{c_{3}}{c_{2}}\right|,\qquad\tau\rightarrow
0.$ (90)

This expression recovers a previous result obtained in BurLew09 . Furthermore,
the case $c_{3}=0$, $c_{5}\neq 0$, yields the following lowest-order behavior

$f_{2}(\tau)\sim\frac{c_{5}^{3/2}\exp\left(-\frac{b^{2}}{2(c_{0}-c_{2}^{2}/c_{4})}\right)}{90\sqrt{15}\pi^{2}\sqrt{|c_{2}|^{3}+c_{0}c_{2}c_{4}}}\tau^{\frac{5}{2}},\qquad\tau\rightarrow
0.$ (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

$f_{2}(\tau)\sim\frac{|c_{4}^{2}-c_{2}c_{6}|\exp\left(-\frac{b^{2}}{2(c_{0}-c_{2}^{2}/c_{4})}\right)}{1296\pi^{2}\sqrt{c_{2}^{2}-c_{0}c_{4}}}\tau^{4},\qquad\tau\rightarrow
0,$ (92)

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

$\iint_{0}^{\infty}d\dot{x}_{1}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}.$
(93)

##  Appendix B Chizhov-Graham theory

An elegant approximation for the FPT problem has been put forward by Chizhov
and Graham ChiGra07 ; ChiGra08 , 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 towards 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$:

$\displaystyle S_{\text{drift}}(t)$
$\displaystyle=\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),$
(94)
$\displaystyle=\frac{1}{2}\left[1+\text{erf}\left(\frac{b(t)}{\sqrt{2}\sigma_{x}(t)}\right)\right]$
(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{d}{dt}\ln(S_{\text{drift}})$ resulting in

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

Following ChiGra07 ; ChiGra08 , we introduced the dimensionless quantity

$T(t)=\frac{b(t)}{\sqrt{2}\sigma_{x}(t)},$ (97)

the temporal derivative of which is given by

$\dot{T}(t)=\frac{1}{\sqrt{2}\sigma_{x}}\left(\dot{b}+\frac{b(\sigma_{x}^{2}-\tau_{\text{m}}\sigma_{xy})}{\tau_{\text{m}}\sigma_{x}^{2}}\right).$
(98)

The moments $\sigma_{x}^{2}(t)$ and $\sigma_{xy}(t)$ have been derived above,
Sec. A.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 ChiGra07 ; ChiGra08 , which assumed
stationary fluctuations with $\sigma_{x}(t)=\text{const}.$. This version with
stationary fluctuations has also been derived in GoeDie08 .

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

$\displaystyle\Phi_{\text{diff}}(b)$
$\displaystyle=\Phi_{\text{diff}}^{\text{wn}}(T)\left[1-\left(1+\frac{\tau_{\text{m}}}{\tau_{\text{s}}}\right)^{-0.71+0.0825(T+3)}\right],$
(99) $\displaystyle\Phi_{\text{diff}}^{\text{wn}}(T)$
$\displaystyle=\tau_{\text{m}}^{-1}\exp\bigl{(}6.1\cdot
10^{-3}-1.12T-0.25T^{2}$ $\displaystyle\quad-0.072T^{3}-0.0117T^{4}\bigr{)}.$
(100)

Here, the numerical coefficients have been fitted to the exact solution
ChiGra08 . The total hazard rate is simply given by the sum of the two limit
cases:

$\lambda(t)=\Phi_{\text{drift}}(b,\dot{b},t)+\Phi_{\text{diff}}(b,t).$ (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}$

$\lambda(t|\hat{t})\approx\Psi_{\text{CG}}\bigl{(}u(t|\hat{t}),\dot{u}(t|\hat{t}),t-\hat{t}\bigr{)}$
(102)

with the Chizhov-Graham link function

$\Psi_{\text{CG}}\bigl{(}u,\dot{u},\tau\bigr{)}=\theta(\tau-
t_{\text{ref}})\bigl{[}\Phi_{\text{drift}}(V_{\text{T}}-u,-\dot{u},\tau-
t_{\text{ref}})\\\ +\Phi_{\text{diff}}(V_{\text{T}}-u,\tau-
t_{\text{ref}})\bigr{]}.$ (103)

###### Acknowledgements.

I would like to thank Sven Goedeke and Markus Diesmann for numerous inspiring
and fruitful discussions throughout this project, especially about
applications to cortical synchronization dynamics, Alexander van Meegen for
interesting discussions on the Stratonovich approximation and for sharing his
unpublished manuscript, and Wulfram Gerstner for his support during part of
this project.

## Conflict of interest

The author declares that he has no conflict of interest.

## Code availability

The code will become available at the following GitHub link after publication:

    https://github.com/schwalger/LIF_hazard_levelcross

## References

  * (1) W. Gerstner, W.M. Kistler, R. Naud, L. Paninski, _Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition_ (Cambridge University Press, Cambridge, 2014)
  * (2) B.W. Knight, J. Gen. Physiol. 59, 734 (1972)
  * (3) W. Gerstner, Neural Comput. 12, 43 (2000)
  * (4) M. Augustin, J. Ladenbauer, F. Baumann, K. Obermayer, PLoS Comput. Biol. 13(6), e1005545 (2017).  DOI 10.1371/journal.pcbi.1005545
  * (5) T. Schwalger, M. Deger, W. Gerstner, PLoS Comput. Biol. 13(4), e1005507 (2017).  DOI 10.1371/journal.pcbi.1005507
  * (6) D.Q. Nykamp, D. Tranchina, J. Comput. Neurosci. 8(1), 19 (2000)
  * (7) A.V. Chizhov, Biol. Cybern. 111(5-6), 353 (2017)
  * (8) L.F. Abbott, C. van Vreeswijk, Phys. Rev. E 48, 1483 (1993)
  * (9) N. Brunel, V. Hakim, Neural Comput. 11, 1621 (1999)
  * (10) N. Fourcaud, N. Brunel, Neural Comput. 14, 2057 (2002)
  * (11) M.J.E. Richardson, Biol. Cybern. 99(4-5), 381 (2008)
  * (12) N. Brunel, J. Comput. Neurosci. 8, 183 (2000)
  * (13) T.C. Potjans, M. Diesmann, Cereb Cortex 24(3), 785 (2014)
  * (14) T. Schwalger, F. Droste, B. Lindner, J. Comput. Neurosci. 39(1), 29 (2015).  DOI 10.1007/s10827-015-0560-x
  * (15) B. Lindner, B. Doiron, A. Longtin, Phys. Rev. E 72(6), 061919 (2005)
  * (16) L. Paninski, Netw. Comput. Neural Syst. 15(4), 243 (2004)
  * (17) W. Truccolo, U.T. Eden, M.R. Fellows, J.P. Donoghue, E.N. Brown, J. Neurophysiol. 93(2), 1074 (2005)
  * (18) J.W. Pillow, J. Shlens, L. Paninski, A. Sher, A.M. Litke, E.J. Chichilnisky, E.P. Simoncelli, Nature 454(7207), 995 (2008)
  * (19) J.W. Pillow, P.E. Latham, in _Adv. Neural Inf. Process. Syst._ (2008), pp. 1161–1168
  * (20) R. Naud, W. Gerstner, PLoS Comput. Biol. 8(10) (2012)
  * (21) J. Brea, W. Senn, J.P. Pfister, J. Neurosci 33(23), 9565 (2013)
  * (22) A. Galves, E. Löcherbach, J. Soc. Fr. Stat. 157, 17 (2016)
  * (23) F. Gerhard, M. Deger, W. Truccolo, PLOS Computat. Biol. 13(2), e1005390 (2017)
  * (24) M.B. Raad, S. Ditlevsen, E. Löcherbach, Ann. Inst. H. Poincaré Probab. Statist. 56(3), 1958 (2020).  DOI 10.1214/19-AIHP1023.  URL https://doi.org/10.1214/19-AIHP1023
  * (25) A.V. Chizhov, L.J. Graham, Phys. Rev. E 75(1), 011924 (2007)
  * (26) A.V. Chizhov, L.J. Graham, Phys. Rev. E 77(1), 011910 (2008)
  * (27) G. Dumont, J. Henry, C.O. Tarniceriu, J. Theor. Biol. 406, 31 (2016).  DOI https://doi.org/10.1016/j.jtbi.2016.06.022
  * (28) T. Schwalger, A.V. Chizhov, Curr. Opin. Neurobiol. 58, 155 (2019)
  * (29) V. Schmutz, W. Gerstner, T. Schwalger, J. Math. Neurosc. 10(5) (2020)
  * (30) V. Schmutz, E. Löcherbach, T. Schwalger.  On a finite-size neuronal population equation (2021)
  * (31) S. Mensi, R. Naud, C. Pozzorini, M. Avermann, C.C.H. Petersen, W. Gerstner, J Neurophysiol (2012)
  * (32) C. Pozzorini, S. Mensi, O. Hagens, R. Naud, C. Koch, W. Gerstner, PLoS Comput Biol 11(6), e1004275 (2015)
  * (33) C. Teeter, R. Iyer, V. Menon, N. Gouwens, D. Feng, J. Berg, A. Szafer, N. Cain, H. Zeng, M. Hawrylycz, et al., Nat. Commun. 9(1), 709 (2018)
  * (34) F. Apfaltrer, C. Ly, D. Tranchina, Netw. Comput. Neural Syst. 17(4), 373 (2006)
  * (35) A. Bulsara, T.C. Elston, C.R. Doering, S.B. Lowen, K. Lindenberg, Phys. Rev. E 53, 3958 (1996)
  * (36) M. Schindler, P. Talkner, P. Hanggi, Physical Review Letters 93(4), 048102 (2004)
  * (37) B. Lindner, J. Stat. Phys. 117, 703 (2004)
  * (38) H.E. Plesser, W. Gerstner, Neural Comput. 12, 367 (2000)
  * (39) A. Herrmann, W. Gerstner, J Comp Neurosci 11(2), 135 (2001)
  * (40) S. Goedeke, M. Diesmann, New J. Phys. 10, 015007 (2008)
  * (41) T. Verechtchaguina, I. Sokolov, L. Schimansky-Geier, Phys. Rev. E 73, 031108 (2006)
  * (42) T. Schwalger, L. Schimansky-Geier, Phys. Rev. E 77, 031914 (2008).  DOI 10.1103/PhysRevE.77.031914
  * (43) J. Schuecker, M. Diesmann, M. Helias, Phys. Rev. E 92, 052119 (2015).  DOI 10.1103/PhysRevE.92.052119
  * (44) B. Lindner, Phys. Rev. E 69, 022901 (2004)
  * (45) T. Schwalger, The interspike-interval statistics of non-renewal neuron models.  Ph.D. thesis, Humboldt-Universität zu Berlin, Mathematisch-Naturwissenschaftliche Fakultät I (2013).  URL http://edoc.hu-berlin.de/docviews/abstract.php?id=40328
  * (46) P. Hänggi, P. Talkner, M. Borkovec, Rev. Mod. Phys. 62, 251 (1990)
  * (47) D.R. Cox, _Renewal Theory_ (Methuen, London, 1962)
  * (48) S.O. Rice, Bell Syst. Tech. J 24, 45 (1945)
  * (49) L.M. Ricciardi, S. Sato, IEEE Trans. Inf. Theory 29, 454 (1983)
  * (50) W. Braun, R. Thul, Phys. Rev. E 95, 012114 (2017).  DOI 10.1103/PhysRevE.95.012114
  * (51) J.M. Azaïs, M. Wschebor, _Level sets and extrema of random processes and fields_ (John Wiley & Sons, 2009)
  * (52) R.L. Stratonovich, _Topics in the Theory of Random Noise_ , vol. 1 (Gordon and Breach, New York, 1967)
  * (53) N.G. van Kampen, _Stochastic Processes in Physics and Chemistry_ (North-Holland, Amsterdam, 1992)
  * (54) L. Badel, Phys. Rev. E 84, 041919 (2011).  DOI 10.1103/PhysRevE.84.041919
  * (55) R.L. Stratonovich, _Topics in the Theory of Random Noise_ , vol. 2 (Gordon and Breach, New York, 1967)
  * (56) M. Puelma Touzel, F. Wolf, PLOS Comput. Biol. 11(12), 1 (2016).  DOI 10.1371/journal.pcbi.1004636.  URL https://doi.org/10.1371/journal.pcbi.1004636
  * (57) A. van Meegen, S.J. van Albada.  A microscopic theory of intrinsic timescales in spiking neural networks (2019)
  * (58) T. Tchumatchenko, A. Malyshev, T. Geisel, M. Volgushev, F. Wolf, Phys. Rev. Lett. 104(5), 058102 (2010)
  * (59) Y. Burak, S. Lewallen, H. Sompolinsky, Neural Comput. 21(8), 2269 (2009)
  * (60) P. Jung, Phys. Rev. E 50, 2513 (1994)
  * (61) W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, _Numerical Recipes in C_ , 2nd edn. (Cambridge University Press, Cambridge, USA, 1992)
  * (62) B. Pietras, N. Gallice, T. Schwalger, Phys. Rev. E 102, 022407 (2020).  DOI 10.1103/PhysRevE.102.022407
  * (63) W. Gerstner, W.M. Kistler, _Spiking Neuron Models: Single Neurons, Populations, Plasticity_ (Cambridge University Press, Cambridge, 2002)
  * (64) Q. Cormier, E. Tanré, R. Veltz, Stoch. Process. Their Appl. 130(5), 2553 (2020).  DOI https://doi.org/10.1016/j.spa.2019.07.010
  * (65) M. DeWeese, W. Bialek, Nuovo cimento D 17, 733 (1995)
  * (66) T. Tchumatchenko, T. Geisel, M. Volgushev, F. Wolf, Front Comput Neurosci. 4(1) (2010)
  * (67) P. Jung, Phys. Lett. A 207, 93 (1995)
  * (68) T. Tchumatchenko, F. Wolf, PLOS Comput. Biol. 7(10), 1 (2011).  DOI 10.1371/journal.pcbi.1002239.  URL https://doi.org/10.1371/journal.pcbi.1002239
  * (69) M. Diesmann, M.O. Gewaltig, A. Aertsen, Nature 402, 529 (1999)
  * (70) J.R. Donoso, D. Schmitz, N. Maier, R. Kempter, J. Neurosci. 38(12), 3124 (2018)
  * (71) S. Lefort, C. Tomm, J.C.F. Sarria, C.C.H. Petersen, Neuron 61(2), 301 (2009)
  * (72) H. Risken, _The Fokker-Planck Equation_ (Springer, Berlin, 1984)