First passage times in integrate-and-fire neurons with stochastic thresholds

We consider a leaky integrate-and-fire neuron with deterministic subthreshold dynamics and a firing threshold that evolves as an Ornstein-Uhlenbeck process. The formulation of this minimal model is motivated by the experimentally observed widespread variation of neural firing thresholds. We show numerically that the mean first passage time can depend non-monotonically on the noise amplitude. For sufficiently large values of the correlation time of the stochastic threshold the mean first passage time is maximal for non-vanishing noise. We provide an explanation for this effect by analytically transforming the original model into a first passage time problem for Brownian motion. This transformation also allows for a perturbative calculation of the first passage time histograms. In turn this provides quantitative insights into the mechanisms that lead to the non-monotonic behaviour of the mean first passage time. The perturbation expansion is in excellent agreement with direct numerical simulations. The approach developed here can be applied to any deterministic subthreshold dynamics and any Gauss-Markov processes for the firing threshold. This opens up the possibility to incorporate biophysically detailed components into the subthreshold dynamics, rendering our approach a powerful framework that sits between traditional integrate-and-fire models and complex mechanistic descriptions of neural dynamics.


I. INTRODUCTION
In vivo recordings of the membrane voltage in many types of neurons display stereotypical upstrokes known as spikes [1]. The timing of these spikes has been shown to exhibit large variability. While fluctuations are often considered as obscuring biological function, the constructive role of stochasticity in neural information processing is now well established [2][3][4].
The random nature of spike generation has attracted considerable attention in the field of mathematical neuroscience. Among the many models that have been suggested to describe neural stochasticity, the class of integrate-and-fire (IF) models [5] has been used to great effect [6]. At the heart of IF models is the notion of neural excitability. A voltage spike is elicited when the membrane potential reaches a threshold. Traditionally, the threshold is considered as either constant or a given function of time that can depend on the spike history, while the membrane potential is described by stochastic differential equations [7][8][9][10][11][12][13][14].
In the present study we investigate the timing of spikes that are generated when a deterministic subthreshold voltage crosses a stochastic firing threshold. A fluctuating threshold reflects experimental findings that the membrane voltage at which a spike is elicited varies. This has been demonstrated e.g. in cortical neurons and hippocampal pyramidal cells [15,16]. Numerous mechanisms have been suggested that could give rise to a variable threshold including adaptation, channel noise and dynamic modulation of the axon initial segment [17- * pmxwb1@nottingham.ac.uk 20]. The physiological relevance of a variable threshold was recently demonstrated in cortical neurons where synchrony detection was significantly improved [15]. In [21] a leaky IF model coupled to a threshold that evolves as a Gaussian process successfully described spiking behaviour of regularly firing stellate neurons within the ventral cochlear nucleus. A stochastic threshold also captures inherent uncertainty in both the detection of spiking thresholds and the spike generation mechanism. The former arises from non-standardised methods to determine threshold values from experimental records [22] and the fact that spikes are generated at the axon initial segment, but are often recorded at the soma. The latter results from still incomplete knowledge of the molecular components that trigger a spike [20].
IF models with a fluctuating threshold have been investigated in the past. However, the approach presented here significantly differs from these studies. Often a distribution of threshold values has been assumed at each time point [23][24][25][26]. As such, successive threshold values are uncorrelated, and in principle arbitrarily large changes of the threshold can occur due to the unbounded support of the chosen threshold distributions. In [27] the concept of hazard functions is used to determine optimal parameters of the firing probability. The tested hazard functions have unbounded support, and the analysis assumes an inhomogeneous Poisson process as the spike generation mechanism [28]. A more general approach was taken in [29] where a subthreshold stochastic process intersects with a frozen stochastic threshold giving rise to a phase transition in the distribution of first passage times.
In principle it is always possible to transform an IF model with a deterministic firing threshold and stochastic subthreshold dynamics with additive noise into a model with a deterministic description of the membrane voltage and a fluctuating firing threshold (see e.g. [30]). Indeed our model is closely related to a leaky IF model driven by white noise and subject to an exponentially decaying firing threshold [31]. The time-dependent firing threshold can be transformed into a constant one giving rise to an exponential drive in the subthreshold dynamics [31]. Generally such explicitly time-dependent stochastic subthreshold dynamics can only be studied numerically. Perturbative expressions for the first and second moment of the first passage time (FPT) have been obtained in the limit of a weak exponential drive [32]. In contrast, we derive analytical expressions of the FPT distribution which do not rely on a small parameter. Our findings provide quantitative insights into non-trivial effects such as the non-monotonic behaviour of the the mean FPT (MFPT) as a function of the noise amplitude.

II. MODEL AND GOVERNING EQUATIONS
In the subthreshold regime the membrane voltage v(t) obeys the ordinary differential equation of a leaky integrator with α and β positive constants. A spike occurs when v(t) hits the threshold h(t) from below, at which point v is discontinuously reset to a constant v r = v(0), i.e.
We hence define the firing times T n by The fluctuating threshold h is given by where h > v r denotes the mean of the threshold and ε > 0 measures the coupling strength to a stochastic process X(t). We model X(t) as an Ornstein-Uhlenbeck process (OUP) whose dynamics reads as [33] Here, γ is the inverse of the correlation time of the OUP, D is a positive constant, and dW is the increment of a standard Wiener process W . We choose the initial condition X(0) = 0 such that h(0) = h and reset the OUP at each threshold crossing to X(0). A fixed value of X(0) renders the above model a renewal process. Without the discontinuous reset of the OUP, consecutive interspike intervals (ISI) are not independent identically distributed and hence do not describe a renewal process [34]. Fig. 1 illustrates the generation of ISIs by the deterministic subthreshold dynamics of v(t) and the OUP in the renewal regime. Note that the full non-stationary correlation function of the OUP enters our analysis, i.e.
for T n ≤ t, t ≤ T n+1 . It is worth noting that ε in Eq. (4) does not need to be small. As both v and X are reset at a threshold crossing such that v r <h, the threshold never reaches negative values although the OUP in general takes on all values on the real line. Throughout this paper we set v r = 0.

III. MEAN FIRST PASSAGE TIMES
We computed the MFPT T n+1 − T n as a function of γ and ε when β/α > h. Angular brackets indicate the average of the firing times T n+1 −T n for the renewal OUP. For notational convenience, we will denote the MFPT by T throughout the paper. In this case firing occurs even for a deterministic threshold at h. Our results are shown in Fig. 2. For small values of γ, the MFPT first increases as a function of ε before decreasing. As we increase γ, the range over which T grows monotonically as well as the amplitude both decrease. The amplitude is defined as the maximal increase of T with respect to the noise-free constant firing time. The latter always exists for the parameter regime under investigation and can be computed from Eqs. (1)-(4) with ε = 0 as Here, H denotes the Heaviside step function, i.e. H(x) = 1 for x ≥ 0 and zero otherwise. For γ = 0.5 the MFPT already decays monotonically. As the OUP can be simulated exactly [35] and the solution for the subthreshold dynamics is known in closed form, the main error for computing ISIs results from the determination of the threshold crossing. We employed two approaches for this. One is based on linear interpolation [36], while the other makes use of a Brownian bridge [37]. Both methods give almost identical results (green and red diamonds in Fig. 2). We used more than 10 6 independent realisations for the computation of the MPFT at each data point. Another approach to compute the MFPT consists of solving the partial differential equation (PDE) which describes the MFPT T = T (v 0 , h 0 ) in dependence on general initial values of the voltage and the threshold, i.e. v(0) = v 0 and h(0) = h 0 . In the present study we are interested in the case h 0 =h and v 0 = 0. Equation (8) can be derived by applying the Feynman-Kac formula [38,39] to the stochastic system given by Eqs.
(1)- (5). In the derivation we used the stochastic differential equation for h, which follows from Eqs. (4)-(5) and Itō's formula [33] as We solve Eq. (8) on a two-dimensional rectangular domain bounded by the line v 0 = h 0 , on which T = 0. We use no-flux boundary conditions on the remaining three sides of the rectangle. The size of the rectangle is chosen such that these three sides are sufficiently far away from the point of interest (v r , h) eliminating any impact of the no-flux boundary conditions on the solution T (v r , h) [40]. Insights into the non-monotonic behaviour of the MFPT can be obtained by introducing a new stochastic variable g = e γt X, where the dynamics for X is given by Eq. (5). Using Itō's formula, the dynamics of g obeys and the time-dependent variance of g follows as We can interpret Eq. (11) as the definition of a new time s in terms of the original time t. We then define a new stochastic process V (s) = g(t(s)) − g(0) in time s where t(s) denotes the inverse of s(t), which always exists since s(t) is strictly monotonically increasing. The variance of V is V 2 (s) − V (s) 2 = s, and together with V (0) = 0 this renders V a standard Brownian motion. This transformation is discussed extensively in [41] and was first used in [42]. From the threshold crossing condition in the new time s, v(t(s)) = h(t(s)) = h + εX(t(s)), we find where we have used which follows from the definitions of g and V . Therefore, the crossing of a voltage v(t) through a threshold h(t) that is described by an OUP is equivalent to the crossing of a new subthreshold voltage v(s) through a standard Brownian motion V (s). Figure 3 shows the new subthrehold voltage v(s) for different values of γ and ε. We observe that for a fixed value of γ the zero crossing of v is independent of ε. It can be shown that this crossing occurs at a value s 0 that corresponds to the deterministic first passage time (Eq. (7)). Note that the value of s 0 grows as we increase γ. For a fixed s < s 0 the value of v(s) increases as we increase ε. Therefore a given realisation of the standard Brownian motion V will intersect with the subthreshold voltage v earlier for larger ε than for smaller ε for times s < s 0 . For γ = 1 almost all crossings of v and V occur before s 0 , which is reflected in the monotonically decreasing behaviour of the MFPT T as a function of ε (Fig. 2). For smaller values of γ there is a significant probability that no threshold crossing occurs before s 0 . Since v(s) is a decreasing function of ε for fixed s > s 0 , a given realisation of V that has not intersected with v before s 0 will cross v earlier for smaller values of than for larger values. This realisation of V has a larger FPT for growing values of ε. The above arguments hold for single realisations of V and hence individual FPT.
To elucidate the increase of the MFPT for growing ε we need to differentiate between different processes. On the one hand more realisations of V could intersect with v at times larger than s 0 . On the other hand the spread of FPTs that are larger than s 0 grows. To determine which factor -or the combination of both -is relevant for the current observations we compute numerically the FPT distributions.

IV. DISTRIBUTION OF FIRST PASSAGE TIMES
Our calculation of FPT distributions is based on an approach developed by Durbin and Williams [43]. They derived an alternating series for the FPT density p of a Brownian motion through a curved boundary. Using the notation of Sec. III for the transformed dynamics in time s, we have with Here, we introduced where we formally set v(s i ) = s i = 0. f i (s) denotes the joint probability density of V (s), V (s 1 ), . . . , V (s i−1 ) at the boundary values v(s), v(s 1 ), . . . , v(s i−1 ). In Eq. (14) the term r k (s) represents the error that is made by truncating the infinite series after k terms. Note that Eq. (16) differs from the original expression since we study downcrossings of a Brownian motion instead of up-crossings as in [43]. For practical purposes it is convenient to introduce F k (s) = k i=1 (−1) i−1 q i (s). We can transform results obtained in the new time s to the original time t via F k (t) = F k (s(t))ds(t)/dt. In Fig. 4 we compare F 2 (t) with direct numerical simulations for different values of γ and ε. For γ = 0.1 we observe that for increasing ε the maximum of the distribution slightly shifts to the left, while the tail of the distribution becomes significantly longer. This wider spread of long FPTs is the reason for the non-monotonic behaviour of the MFPT as it dominates the increase of trajectories with shorter FPTs. grow as we increase ε, but the shift of the distribution to shorter times is more pronounced. The latter leads to the monotonic decrease of the MFPT for increasing ε. The occurrence of long FPTs for small values of γ can be attributed to the correlation time of the OUP τ OU = 1/γ. For small values of γ, changes in the OUP occur slowly. Since for times larger than 1/α, the voltage has almost reached its equilibrium β/α, the OUP can then spend a considerable time in the vicinity of the voltage without crossing it. It is worth noting that F 2 (t) captures the full FPT distribution extremely well for most of the probability mass. F 2 (t) only fails at correctly predicting the tail of the FPT distribution. The extent to which F k (t) agrees with the true distribution depends on the k. As we illustrate in Fig. 5 F k (t) approximates results from direct simulations better as we increase k. In their  Fig. 3 shows, however, that the intercept is negative for small values of s and positive for large values. Our results suggest that convergence of the FPT distribution also occurs if the above assumption is violated, but a rigorous proof is still missing.

V. CUMULATIVE DENSITY FUNCTIONS
We next studied the emergence of long tails in the FPT distributions in more detail. In Fig. 6 we show the cumulative distribution function (CDF) of the FPTs for different values of γ. We employed the method developed by Wang and Pötzelberger [44]. The CDF is expressed as an expectation value g using a piecewise linear approximation of the boundary v. An advantage of this method is that each instance of g is known in closed form and only the average needs to be computed numerically. We observe that for times smaller than 2 the two CDFs agree well. This entails that the FPT distributions almost co- incide since they follow from the corresponding CDFs by differentiation with respect to time. The CDF for γ = 1 levels off at significantly earlier times than the CDF for γ = 0.1 illustrating the more pronounced tail of the FPT distribution for γ = 0.1. For comparison we also include results from direct numerical simulations, which lie on top of the quasi-analytical results for the CDF. The findings in Fig. 6 provide an independent test for the non-monotonic behaviour of the MFPT. Moreover, convergence of the CDF as we increase the number of linear segments to approximate v is guaranteed.
The shift of the FPT distributions to smaller times can also be understood by computing the probability c that a given realisation of V (s) intersects with v(s) at times s < s 0 . We computed c based on the approach in [44] and plot results as a function of ε for different values of γ in Fig. 7. We observe that c is an increasing function of ε. Hence a growing number of realisations of V have FPTs that are smaller than s 0 . As we increase γ for fixed ε we find that c grows monotonically demonstrating that more realisations of V reach v before s 0 . The rate of increase of c with ε is more pronounced for larger values of γ. Therefore, the relative increase in realisations of V that have FPTs smaller than s 0 is stronger for larger values of γ compared with smaller values. Taken together the results in Figs. 6 and 7 demonstrate the subtle interplay between the emergence of long tails in the FPT distribution and the increase in the probability mass of FPTs smaller than T det in the generation of non-monotonic MFPTs.

VI. DISCUSSION
In this paper, we have formulated and investigated an IF model for neural firing with a stochastic process controlling the firing threshold, motivated by experimental studies showing apparent variability or randomness in the spiking threshold. The model is designed to be as simple as possible while allowing fluctuations of the threshold around a 'preferred' valueh. We therefore use an OUP reverting to a meanh for the threshold, and a leaky IF model for the sub-threshold dynamics v. The simplicity of the model permits the use of techniques from stochastic calculus [38].
We have demonstrated that the MFPT may depend non-monotonically on the noise amplitude , finding a maximum of the MFPT as the noise strength is increased. The occurrence of a minimum in spiking activity at non-vanishing noise strength has been termed inverse stochastic resonance and has been described in stochastic Hodgkin-Huxley models [45]. The emergence of inverse stochastic resonance in these models is directly linked to the phase-space structure of the corresponding deterministic models. In contrast the deterministic LIF model does not predict the complex behaviour of its stochastic counterpart. Inverse stochastic resonance has been observed in vivo [46] and may even have therapeutic applications [47].
Local extrema of the MFPT have been reported before. In [32] escape from a linear and parabolic potential in the presence of a weak exponential driving force is considered, while a Kramers problem with a fluctuating barrier is investigated in [48]. In both studies extrema of the MFPT emerge when a time scale is varied, which is the decay rate of the weak exponential forcing in [32] and the correlation time of the random barrier in [48]. The dependence of the MFPT as a function of the noise amplitude is always monotonic in these studies. This sets our results apart from these earlier findings. Moreover, transforming our model to one with a constant threshold and an exponential drive in the subthreshold voltage results in an amplitude of the time-dependent drive that is not small for the parameter values that we have considered. This precludes us from applying the perturbative results in [31,32] for computing the first two moments of the MFPT. An exponential drive of sufficient amplitude and decay time is essential for observing the non-monotonic behaviour of the MFPT. The non-monotonicity results from the competition between small passage times, which arise when the OUP crosses v during the rising phase of the subthreshold voltage, and large passage times, which occur when the OUP intersects with the almost stationary value of v.
The dependence of the MFPT on the noise strength was found both by direct Monte Carlo simulations and by solving a backward Kolmogorov-type PDE, with virtually identical results. The increase of the MFPT with noise strength can be explained quantitatively by inspection of the FPT distributions. As the noise becomes stronger, the peak of the FPT distribution moves towards shorter firing times, but the tail of the distribution strengthens significantly, so that rare long times between firings become more frequent. Further understanding of this effect was achieved by transforming the dynamics of the system to that of a standard Brownian motion. This allows a series method [43] to be used to obtain approximations to the FPT distribution analytically; this series appears to converge rapidly, showing good agreement with numerical results.
A number of previous studies on LIF models with constant threshold employed different noise intensities of the OUP. For example, a noise scaling of √ 2Dγ was used in [49], which entails that the variance of the OUP is given by D and hence does not depend on the correlation time 1/γ. By setting the noise intensity to √ 2Dγ the authors of [50] ensure that the variance of the subthreshold voltage remains almost constant in the limit of long correlation times of the OUP, for which the variance is given by Dγ. Since the emergence of the maximum of the MFPT depends on the interplay of long correlation times and sufficient variance of the OUP we tested both of the above noise intensities. Our results show that the new noise scalings of the OUP have no effect on the nonmonotonicity of the MFPT. The maximum still exists but is shifted to larger values of ε for a given value of γ.
An interesting question is how broadly applicable the transformation of a stochastic threshold to Brownian motion is, particularly since an IF model with additive noise subthreshold dynamics and constant firing threshold can always be transformed into a model where the subthreshold voltage evolves deterministically and hits a stochastic threshold. If the stochastic threshold is given by a Gauss-Markov process a transformation to Brownian motion always exists [41]. A wider class of practically relevant stochastic threshold processes are general diffusion processes [33]. Here a transformation to Brownian motion is only feasible for specific forms of the drift and diffusion functions. Wang and Pötzelberger [51] provide a practical PDE that these functions need to satisfy to derive an equivalent Brownian motion. For the most general conditions to recast a diffusion process as Brownian motion we refer the reader to [52].
Future work will consider two variations of the set-up of the problem investigated here. We have considered a 'renewal' model in which the threshold is reset to a constant value after firing, but an alternative procedure is to set it to a value selected from the statistical distribution of the OUP. A second possibility is to use a full spike train, in which the threshold OUP is not reset after each spike, but evolves freely. This will give rise to serial correlations of the interspike intervals, and it will be of interest to investigate whether the serial correlations display stronger or weaker correlations than the OUP threshold process itself.