Approximation of the first passage time density of a Wiener process to an exponentially decaying threshold by two-piecewise linear threshold. Application to neuronal spiking activity

The first passage time density of a diffusion process to a time varying threshold is of primary interest in different fields. Here we consider a Brownian motion in presence of an exponentially decaying threshold to model the neuronal spiking activity. Since analytical expressions of the first passage time density are not available, we propose to approximate the curved boundary by means of a continuous two-piecewise linear threshold. Explicit expressions for the first passage time density towards the new boundary are provided. Then we introduce different approximating linear threshold and describe how to choose the optimal one minimizing the distance to the curved boundary and hence the error in the corresponding passage time density. Theoretical means, variances and coefficients of variation given by our method are then compared with empirical quantities from simulated data as well as other firing statistics derived under the assumption of a small amplitude of the time-dependent change in the threshold. Finally maximum likelihood and moment estimators of the parameters of the Wiener process are also proposed and compared on simulated data.


1.
Introduction. Stochastic models have been extensively used in theoretical neuroscience since the pioneer work by Gerstein and Mandelbrot in 1964 [12]. There they considered a Wiener process (also known as Brownian motion or Perfect-Integrate-and-Fire model) to model the voltage across the membrane. An action potential, also known as spike, is generated whenever the membrane potential reaches a certain constant threshold. After that, the membrane voltage is reset to its resting value and the evolution restarts. From a mathematical point of view, a spike is the first passage time (FPT) of a stochastic process to a constant threshold. The collection of spike epochs of a neuron, called spike train, defines a renewal process, with independent and identically distributed inter-spike intervals (ISIs). Despite the excellent fit with some experimental data, the Gerstein-Mandelbrot model was criticized because it disregards features involved in neuronal coding.

MASSIMILIANO TAMBORRINO
A first extension, combining both mathematical tractability and biological realism, is represented by Leaky-Integrate-and-Fire (LIF) models [28,36]. Despite some criticisms on the lack of fit of experimental data [16,31], these models are still largely used.
Another common generalization is represented by Wiener processes (or more generally LIF models) with time-dependent threshold [35,37]. These models can be chosen to reproduce biological features such as the afterhyperpolarization in neurons. For exponentially decaying thresholds, these processes can be used to model a neuron with an exponential time-dependent drift, as shown by Lindner and Longtin [19]. They investigated the effect of an exponentially decaying threshold on the firing statistics of a stochastic integrate-and-fire neuron [19]. Using a perturbation method [18], they derived analytical expressions of the firing statistics under the assumption that the amplitude of the time-dependent change in the threshold is small. These statistics are useful to characterize the spontaneous neural activity and to investigate the neuronal signal transmission. In particular, they can suggest under which conditions a decaying threshold may facilitate or deteriorate signal processing by stochastic neurons. For a Wiener process, these quantities can also be obtained using the approach in [38]. Also this method assumes a small amplitude , but it has the advantage of providing an explicit approximation of the FPT density.
Here we consider a Wiener process with exponentially decaying threshold. The first aim of the paper is to provide an alternative method to approximate the firing statistics and the FPT density for any possible amplitude , extending the results in [19,38]. Different estimators are proposed, as mentioned in Section 1.1 and discussed in Section 4.2. Means, variances, coefficients of variation (CVs) and distributions of the FPTs are compared on simulated data and the most suitable are recommended. A comparison with the results in [19,38] under the assumption of a small amplitude is also performed. The second aim of this work is the estimation of drift and diffusion coefficients of the Wiener process. Maximum likelihood and moment estimators are derived and evaluated on simulated data. Our results show a good approximation of both firing statistics and parameters of the underlying model. Although the considered model generates a renewal process, the proposed method can also be applied to non-renewal processes, e.g. adaptive threshold models [8,17]. Recently, an increasing interest arose towards these models, interest motivated by the excellent fit of the firing statistics of electrosensory neurons [7,9]. The novelty of these models is that the threshold has a jump immediately after a spike. Since the boundary depends on the previous firing epochs, the ISIs are not independent anymore. However, the distribution between two consecutive spikes, conditioned on the initial position of the threshold, is the same of that studied here. Hence, our results may represent a first step towards an understanding of the more complicated adapting-threshold models.
1.1. Mathematical background. FPTs of diffusion processes to constant or timedependent thresholds have been extensively studied in the literature. Explicit expressions for constant thresholds are available for the Wiener process [10,11], for a special case of the Ornstein Uhlenbeck (OU) process [26], for the Cox-Ingersoll-Ross process [6], and for those processes which can be obtained from the previous through suitable measure or space-time transformations, see e.g. [2,6,25]. For most of the processes arising from applications and for time-varying thresholds, analytical expressions are not available. Numerical algorithms based on solving integral equations have been proposed in [4,5,27,29,32,34], while approximations based on Monte-Carlo path-simulation methods in [13,14,20].
A different approach to tackle the FPT problem consists in focusing directly on the two-sided boundary crossing probability (BCP), i.e. the probability that a process is constrained to be between two boundaries. If one of the boundary is set to −∞, the resulting one-sided BCP equals the survival probability of the FPT to the other boundary [39]. Explicit formulas for the BCP of a standard Wiener process for continuous and piecewise-linear thresholds are known (see [3,22,23,39,40]). In general, the BCP of a diffusion process through an exponential decaying threshold is available only for those processes which can be expressed as a piecewise monotone functional of a standard Brownian motion. Examples are the OU process or the geometric Brownian motion with time dependent drift for specific parameter values [40]. The simple but powerful idea is to approximate both one and two-sided curvilinear BCPs by similar probabilities for close boundaries of simpler form, namely n piecewise-linear thresholds, whose computation of the BCP for Wiener is feasible. Under some mild assumptions, the approximated two-sided BCP converges to the original one when n → ∞ [40], with rate of convergence given in [3,39].
For the exponential decaying threshold considered in this paper, the convergence can be obtained by choosing piecewise linear thresholds approximating the curved boundary from above and below, with approximation accuracy given by their distance [40]. However, all the available formulas for the BCPs require either Monte-Carlo simulation methods or heavy numerical approximations.
Here we consider a two-piecewise linear threshold as an approximation of the curvilinear boundary. Since n = 2, the asymptotic convergence of the BCPs does not hold. However, we can derive analytical expression for the FPT density to the two-piecewise linear boundary, and use it as an approximation of the unknown FPT density. Four possible piecewise thresholds are proposed and optimized to minimize the distance to the original threshold.

2.
Model. We describe the membrane potential evolution of a single neuron by a Wiener process X(t), starting at some initial value x 0 . We assume X(t) given as the solution to a stochastic differential equation where W (t) is a standard (driftless) Wiener process. The drift µ > 0 and the diffusion coefficient σ > 0 represent input and noise intensities, respectively. A spike occurs when the membrane potential X(t) exceeds the exponentially decaying for the first time. Here, δ k denotes the time of the kth spike for k > 0, and can be interpreted as a relative refractory period. We set δ 0 to be the starting time of the process, i.e. δ 0 = t 0 . The term λ represents the decay rate of the threshold, while is interpreted as the amplitude of the time-dependent change in the boundary. After a spike, the membrane potential is reset to its resting position x 0 < b 0 + , and its evolution is restarted, as illustrated in Fig. 1. The presence of δ k in (2) ensures that the ISIs are independent and identically distributed. Denote by b(t) the threshold b * (t) for k = 0, i.e.
Then, all ISIs are distributed as the FPT of X to b(t), namely Quantities of interest are the probability density function (pdf) and the cumulative distribution function (cdf) of T b , denoted by f T b and F T b , respectively. Another relevant quantity is the two-sided BCP given by Here τ > t 0 is fixed, boundaries a(t) and c(t) are real functions satisfying a(t) < c(t) for all t 0 < t ≤ τ and a(t 0 ) < x 0 < c(t 0 ). Setting a(t) = −∞ yields the one-sided BCP P X (−∞, c, τ ) = P(T c > τ ) = 1 − F Tc (τ ), which corresponds to the survival probability of T c . For a standard Wiener process W , Wang and Pötzelberger [40] showed that, if the sequences of piecewise linear functions a n and c n converge uniformly to a(t) and c(t) on [t 0 , τ ] respectively, then, for the continuity property of probability measure, it holds lim n→∞ P W (a n , c n , τ ) = P W (a, c, τ ).
When a(t) = −∞ and c(t) = b(t), the convergence of P(T bn > τ ) to P(T b > τ ) can be obtained by choosing piecewise linear thresholds approximating b(t) from above, denoted by b + n (t), or from below, b − n (t). That is, b + n (t) ≥ b(t) and b − n (t) ≤ b(t), ∀t ∈ [t 0 , τ ], respectively. Since the considered curved boundary is convex, we have [23] i.e.

The approximation accuracy is given by
, with bounds given in [3]. Obviously, the accuracy in the BCP increases when the distance between the two thresholds decreases.
3. FPT to continuous piecewise linear threshold. The transition density function of a standard Brownian motion in x 1 , x 2 , . . . , x n at time t 1 , t 2 , . . . , t n , constrained to be below the absorbing threshold c(t) defined by n piecewise-linear threshold over [t 0 , t n ], is given in [39]. Extending that result to the case of a Brownian motion with drift µ and diffusion coefficient σ, starting in x 0 < c(t 0 ) = c 0 at time t 0 , we obtain where c i = c(t i ) and HITTING TIME TO EXPONENTIALLY DECAYING THRESHOLDS where δ k denotes the kth spike. The membrane potential starts in x 0 = 0 at time δ 0 := t 0 = 0, and it evolves until it hits the boundary for the first time. Then, a spike is generated, the voltage X(t) is reset to its resting potential x 0 , the threshold is reset to b 0 + and the evolution restarts. For the considered spiking generation rule, all ISIs are independent and identically distributed. Here the parameters are µ = 1, σ 2 = 1, b 0 = 1, λ = 1 and = 5.
Since we approximate b(t) by means of a continuous two-piecewise linear threshold, we denote byb the linear threshold c(t) when n = 2. We havẽ where 1 A denotes the indicator function of the set A and α 1 , α 2 , β 1 , β 2 ∈ R.
Throughout the paper, we set α 2 = α 1 + β 1 (t 1 − t 0 ) to guarantee the continuity ofb(t). This allows to provide analytical expressions of (4) and (6), which we use as an approximation of f T b . In particular, we have P(Tb < t) = P(Tb < t, Tb < t 1 ) + P(Tb < t, Tb > t 1 ) with pb 1 (x 1 , t 1 ) given by (4) for n = 1. Mimicking [33], by taking the derivative of (9) with respect to t, and plugging (7) in it for proper values of α and β, we get This result extends that for a driftless Brownian motion, see e.g. [1,30]. As expected, setting α 1 = α 2 = α and β 1 = β 2 = β yields the pdf of the FPT of a Wiener process to a linear threshold c(t) = α + β(t − t 0 ). By definition, the first two moments and variance of Tb are given by and can be numerically computed.

4.1.
Parameter estimation of the piecewise-linear threshold. The primary aim of this paper is the approximation of the FPT distribution (and relevant statistics) for a curved boundary b(t), by means of the FPT distribution for a continuous two-piecewise linear thresholdb(t). As discussed in Section 2, the quality of the approximation improves when the distance betweenb and b decreases.
Denote byθ + ,θ − ,θ betw andθ free the estimators of θ from the boundaries b + , b − , b betw and b free , respectively. From (3), it follows that the best approximation of P X (−∞, b, τ ) is obtained when the distance between b + and b − is minimized. For this reason, we defineθ + andθ − as the estimators minimizing the area of the squared distance between the two boundaries, i.e.
is chosen to avoid possible numerical issues in the optimization procedure.
Once b + and b − have been computed, the estimatorθ betw is defined aŝ i.e. it is the equidistant line from b + (t) and b − (t).
Finally, the estimatorθ free is the one minimizing the area of the squared distance between the piecewise-line and the curved boundary, i.e.
Note that the estimation of θ does not depend on the observations of the FPTs but it is performed theoretically under the assumption that the parameters λ, and b 0 of b(t) are known. The proposed estimators and their assumptions are summarized in Table 1. All the minimizations have been performed in the computing environment R [24]. Since the parameter values need to fulfil some conditions (cf. Table 1), minimizing the areas is a constrained optimization problem. We perform it by means of the built-in R function optim, penalizing those parameter values not fulfilling the conditions by returning 10 10 .

4.2.
Parameter estimation of the process. The second aim of the paper is the estimation of the parameters µ and σ 2 of the Wiener process from a sample {r i } n i=1 of n independent observations of T b . That is, we want to estimate φ = (µ, σ 2 ) under the assumption that the parameters of the threshold are known.

4.2.1.
Maximum likelihood estimator of φ. First, we derive the parameters θ of the thresholdb(t), as described in Section 4.1. Then, we use maximum likelihood estimator (MLE) as follows. Since the observations are independent and identically distributed, the log-likelihood function is given by where f Tb is the pdf given by (10) with θ replaced byθ. Then, the log-likelihood function can be maximized numerically to obtain the unknown parameter φ. Since the parameter values of µ and σ need to be positive, when minimizing l r (φ) by means of the function optim, we penalize negative values of µ and σ by returning 10 10 . We denote byφ MLE (θ) the MLE of φ derived from the thresholdb with parametersθ.

Moment estimator.
A different approach consists in equating the theoretical moments of Tb, given by Eq. (11), with the empirical moments of T b . In particular, we numerically solve a system of two equations (given by the first two moments) in the two unknown parameters φ = (µ, σ 2 ). We denote byφ ME the moment estimator (ME) of φ.
Estimator Assumption onb(t) on [τ 0 , τ * ] Unknown parameters Parameter conditionŝ Table 1. Proposed estimators of the parameters of the piecewise linear thresholds b + , b − , b betw and b free given in Section 4.1 under different assumptions.
When is small, approximated mean and variance of T b are available [19,38]. In particular, for a general parameter b 0 > x 0 , we have We denote byφ ME the moment estimator of φ obtained from (13) and (14) when is small.

Monte Carlo simulations.
We simulate FPTs of the Wiener process X to b(t) as described in [19,28]. Applying the Euler-Maruyama scheme to the stochastic differential equation (1), we generate realizations of X, denoted by x i := X(s i ), at discrete times s i = i∆s, i ≥ 1. We set X 0 = x 0 = 0 and ∆s = 0.001 as time step. To avoid the risk of not detecting a crossing of the boundary due to the discretization of the sample path, at each iteration step we compute the probability that the bridge process time s i and constrained to be in x i+1 < b(s i+1 ) at time s i+1 , crosses the threshold in between s i and s i+1 . For a Wiener process, this probability is given by [15] P( A FPT is observed if x i hits/exceeds the threshold b at time s i , i.e.  [19,38]. Finally, for each value of σ 2 , and λ, we repeat simulation of data set 1000 times, obtaining 1000 statistically indistinguishable and independent trials.

Set up.
In the simulations we are mainly concerned to illustrate the performance of our method under the assumption that the threshold b(t) is known, i.e. b 0 , the rate λ and the amplitude are given. Two scenarios are considered: both µ and σ 2 are known; no information about the parameter of the Wiener process is given. In the first case it is of interest to evaluate the error in the estimation of mean, variance, CV and cdf of T b by comparing theoretical (11) and empirical firing statistics. When is small, a further comparison with (13) and (14) is carried out.
To measure the error in the estimation of F T b , we consider the relative integrate absolute error (R IAE ), defined as We replace the unknown quantities F T b and E[T b ] by their empirical estimators, defined by F n (t) = 1 T bi /n, respectively. Both empirical quantities are based on n = 1000000 simulated FPTs, ensuring the closeness to the theoretical counterparts by the law of large numbers. This first scenario is meant to understand the goodness of our approximation through simulations.
Another relevant question is the performance of the MLEs and MEs of µ and σ 2 , as described in Section 4.2. To compare different estimators, we use the relative mean error R ME to evaluate the bias and the relative mean square error R MSE , which incorporates both the variance and the bias. They are defined as the average over the 1000 repetitions of the quantities and likewise for σ 2 .

5.3.
Theoretical results for cdf and firing statistics of T b . In Fig. 3 are reported theoretical and empirical means, variances and CVs of T b as a function of the rate λ, for small values of the amplitude and for σ 2 = 0.2, 0.4 and 1. The given theoretical quantities are obtained from (11) for the piecewise linear threshold b free . Note how the mean of T b does not depend on σ 2 , as it also happens for a linear threshold, while both variance and CV increase with growing σ 2 . We refer to [19] for a detailed discussion on other qualitative features of the firing statistics, e.g. monotonic decrease on the mean with growing λ, existence of a minimum value for the variance, etc. What is relevant to emphasize is the excellent fit of the firing statistics provided by our method for any λ, and for both small (cf. Fig. 3) and large (cf. Fig. 4) values of . When is small, our theoretical firing statistics are at least as good as those in [19,38], outperforming them when grows. The firing statistics of T b betw are almost identical to those of T b free , while those of T b+ and T b− are slightly different for increasing . This can be seen in Fig. 5, left panel, where the percentages of the R IAE (F T ) for the four proposed estimators are given. As expected, the best approximation of F T b is provided by F T b free , since b free is the only threshold whose parameters are obtained from a non-constrained optimization problem. The performance of the estimators gets worse for large σ 2 and . The highest error is observed for the value of λ that minimizes the variance of T b . However, all errors are smaller than 2%, confirming the good performance of the proposed estimators.

5.4.
Parameter estimation of (µ, σ 2 ). We have seen that T b free yields the best approximation of T b in terms of both cdf and firing statistics. For this reason, we limit our study to the estimatorsφ based on b free . In Fig. 6 the R ME and the R MSE ofμ andσ 2 are reported. As expected, the MLE provides the best estimate of φ, while both MEs are acceptable only for small values of . The performance ofμ is highly satisfactory, with R ME (μ) smaller than 0.5%, and R MSE (μ) < 0.2%. Larger but still good R ME and R MSE are observed forσ 2 . The performance ofφ MLE gets worse for growing σ 2 , as shown in Fig. 7. However, except the R ME (σ 2 ) for large values of , all errors are between 0 and 2 − 3%. Two last remarks should be done: Top panels: σ 2 = 0.2. Central panels: σ 2 = 0.4. Bottom panels: σ 2 = 1. Empirical quantities from simulations (symbols), theoretical quantities given by (11) for the piecewise linear threshold b free (solid lines), and theoretical quantities (13) and (14) when is small (solid gray lines). Also shown are the firing statistics of T b when = 0 (horizontal dashed lines).
first, the R MSE ofμ for small values of approaches the corresponding values of σ 2 . Second, R MSE (σ 2 ) seems not to depend on λ, and σ 2 , but to be equal to 2%. This error decreases when increasing the sample size. For example, the R MSE (σ 2 ) ≈ 1% when n = 200 (results not shown).
6. Discussion. As a consequence of the recent increasing interest towards adaptingthreshold models for the description of the neuronal spiking activity, a need of suitable mathematical tools to deal with hitting times of diffusion processes to   (13) and (14) when is small (solid gray lines). Also shown are the firing statistics of T b when = 0 (horizontal dashed lines). R ME (μ) in percentage q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q λ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q ) in percentage q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q λ q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q ) in percentage Figure 6. Dependence of R ME (μ), R MSE (μ), R ME (σ 2 ) and R MSE (σ 2 ) (average over 1000 simulations) on λ and when X is a Wiener process with µ = 1 and σ 2 = 0.2. Different estimators of φ = (µ, σ 2 ) are considered: maximum likelihood estimatorφ MLE (solid lines with triangles), moment estimatorφ ME (dashed lines with circles) and moment estimator from (13) and (14) when is small,φ M E (gray solid lines with gray circles). The values of between consecutive vertical dotted lines are fixed and equal to 0.05, 0.1, 0.2, 1, 5, 10, while λ varies between 0.02 and 10.
represented by the work of Wang and Pötzelberger, who provide an explicit expression which should then be evaluated through Monte-Carlo simulations. The idea behind the works of Lindner and Longtin and of Urdapilleta, was to simplify some mathematical difficult equations arising from the study of the FPT by linearizing them in , the amplitude of the decaying threshold. As a consequence, the quality of the approximation rapidly decreases when increases. R ME (μ) in percentage q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q  . Dependence of R ME (μ), R MSE (μ), R ME (σ 2 ) and R MSE (σ 2 ) (average over 1000 simulations) on λ, and σ 2 when X is a Wiener process with µ = 1 and σ 2 equal to 0.2 (solid lines with circles), 0.4 (dashed lines with triangles) and 1 (gray solid lines with gray circles). Here only the maximum likelihood estimatorφ MLE of φ = (µ, σ 2 ) is considered. The values of between consecutive vertical dotted lines are fixed and equal to 0.05, 0.1, 0.2, 1, 5, 10, while λ varies between 0.02 and 10.
The method proposed here has no restriction on the parameter of the thresholds and it is based on the simple idea of replacing the boundary by a continuous twopiecewise linear threshold. This allows us to derive the analytical expression of the FPT density to the two-piecewise threshold, and to use it to approximate the desired distribution. To some extent, the presence of two linear thresholds can be considered as a second order approximation of the problem.
Numerical simulations show a good performance of the proposed method both when computing the main firing statistics, such as means, variances and CVs, and when calculating the FPT distribution. Different approximating thresholds have been proposed. We suggest choosing the one minimizing the distance with the curvilinear threshold and to restrict the interval where to perform the optimization as described in the paper. Among the estimators of the drift and diffusion coefficients of the Wiener process, we suggest applying MLE which always estimates the parameters reasonably well.
The method proposed here may yield several interesting developments. First of all, it can be used to characterize the firing statistics of the Wiener process to the exponential decaying threshold, extending the previous considerations obtained for small values of . Then, it may be extended to the case of a Wiener process with an adapting decaying threshold, as suggested in the introduction. Finally, our results may also be applied to all those processes which can be expressed as a piecewise monotone functional of a standard Brownian motion [40], as well as to Wiener processes with time-varying drift [19,21].