Skip to main content

Fokker–Planck and Fortet Equation-Based Parameter Estimation for a Leaky Integrate-and-Fire Model with Sinusoidal and Stochastic Forcing

Abstract

Analysis of sinusoidal noisy leaky integrate-and-fire models and comparison with experimental data are important to understand the neural code and neural synchronization and rhythms. In this paper, we propose two methods to estimate input parameters using interspike interval data only. One is based on numerical solutions of the Fokker–Planck equation, and the other is based on an integral equation, which is fulfilled by the interspike interval probability density. This generalizes previous methods tailored to stationary data to the case of time-dependent input. The main contribution is a binning method to circumvent the problems of nonstationarity, and an easy-to-implement initializer for the numerical procedures. The methods are compared on simulated data.

List of Abbreviations

LIF: Leaky integrate-and-fire

ISI: Interspike interval

SDE: Stochastic differential equation

PDE: Partial differential equation

1 Introduction

Information processing in the nervous system is carried out by spike timings in neurons. To study the neural code in such a complicated system, a first step is to understand signal processing and transmission in single neurons. Stochastic leaky integrate-and-fire (LIF) neuronal models are a good compromise between biophysical realism and mathematical tractability, and are commonly applied as theoretical tools to study properties of real neuronal systems. A central issue is then to perform statistical inference from experimental data and estimate model parameters. Many electrophysiological experiments on neurons, namely extra-cellular recordings, are only capable of detecting the time of the spike and not the detailed voltage trajectory leading up to the spike. Estimating the parameters of the LIF model from this type of data is equivalent to estimating the parameters of a stochastic model from the statistics of the first-passage times only. A common assumption is that the data are well described by a renewal process, thus basing the statistical inference on the interspike intervals (ISIs), assuming these are realizations of independent and identically distributed random variables. Since only partial information about the process is available, the statistical problem becomes more difficult, and no explicit expression for the likelihood is available.

Different methods have been proposed. In the seminal paper [1], a point process approach is proposed. The spike trains of a collection of neurons are represented as counting processes. Time is discretized and the point processes approximated by 0–1 time series. Then the probability of firing in the next time interval is modeled as a function of the spike history. In this way, maximum likelihood estimation is feasible. External stimuli are not considered. In [2], a numerically involved moment method is developed. It uses the first two moments of the first-passage times of the Ornstein–Uhlenbeck process to a constant threshold, which are given as series expressions, and equates them to their empirical counterparts. In [3, 4], certain explicit moment relations derived from the Laplace transform of the first-passage time distribution are applied, but these are only valid under stimulation (supra-threshold regime). In [5], inference is based on numerical inversion of the Laplace transform. In [6], a functional of a three-dimensional Bessel bridge is applied to obtain a maximum likelihood estimator. None of these methods are feasible to extend to the time-inhomogeneous case, which is of our interest. In [7, 8], an integral equation is used to derive an estimator in the time-homogeneous setting. This approach is readily extended to time varying input, which we will explore in this paper. Some of the above methods are compared in [9]. Finally, a review of estimation methods is provided in [10].

Many sensory stimuli, like sound, contain an oscillatory component [11, 12]. Such inputs will cause oscillating membrane potentials in the neuron, generating rhythmic spiking patterns. The oscillation frequency determines the basic rhythm of spiking, and is considered to be significant for neuronal information processing. The dynamics of periodically forced neuron models have been extensively studied; see [1318] and references therein. Even so, attempts to solve the estimation problem in these nonstationary settings have been rare. One problem is that the ISIs are no longer independent nor identically distributed. In [19], a more complicated model with linear filters is considered, allowing also for the spike history to influence the membrane potential dynamics. The estimation problem is solved through numerical solutions to the Fokker–Planck equation, and it is shown that the log-likelihood is concave, thus ensuring a global maximum; see also [20, 21]. Because their model is more involved, some approximations to the solution of the Fokker–Planck equation are applied, to ensure acceptable computing times. We will apply the full Fokker–Planck equation to solve our estimation problem, since the computing time is always lower than 2 seconds for a sample size of 1000 spikes.

In this paper, we thus describe and discuss two methods to estimate parameters of LIF models with the added complexity of a time-varying input current. We assume that the time-varying current is a sinusoidal wave, but we believe that the approaches generalize to an arbitrary periodic forcing with known frequency. One approach relies on the Fortet integral equation, which is readily extended to the time-inhomogeneous case. An advantage of this approach is that if the transition density of the diffusion in the LIF model is known, as is the case for the Ornstein–Uhlenbeck and the Feller model, the computational burden is limited. A second approach involves numerical solution of the Fokker–Planck equation, where the time-dependence is explicitly accounted for. After a numerical differentiation, the likelihood function can be calculated providing the maximum likelihood estimator. Nevertheless, we chose an alternative loss function which seems marginally more robust, directly comparing the survival function provided by the solution of the Fokker–Planck equation with its empirical counterpart. The two approaches give similar results and they are more carefully compared in the supplementary online material. Both methods need sensible starting values for the optimization algorithms, and we provide an easy-to-implement initializer. The estimation procedures are compared on simulated data and we find that both algorithms are able to find estimates close to the true values for several different dynamical regimes. We find that for small sample sizes the Fokker–Planck algorithm can be considered marginally preferable, whereas for larger sample sizes the Fortet algorithm becomes marginally superior. Moreover, at high frequencies of the sinusoidal forcing, the Fortet is better at identifying the parameters, though in general there is less information in the data to distinguish between a constant input and the amplitude of the periodic forcing.

2 Model

The time evolution of the voltage of a spiking neuron is modelled by a stochastic process, V, given as solution to the following stochastic differential equation (SDE):

d V ( t ) = ( μ V ( t ) τ + A sin ( ω t ) ) d t + σ d W ( t ) , t 0 = 0 ; V ( t 0 ) = v 0 , t n = inf { t > t n 1 : V ( t ) = v th } for  n 1 , { V ( t n + ) = v 0 , J n = t n t n 1 .
(1)

Here, μ is a bias current acting on the cell, τ is the decay time, A and ω are the amplitude and (angular) frequency of the sinusoidal current acting on the cell, σ is the strength of the stochastic fluctuations, W= { W t } t 0 is a standard Wiener process, and t n + denotes the right limit taken at t n . A spike occurs when the membrane voltage V(t) crosses a voltage threshold, v th , and then V(t) is instantaneously reset to the resting potential v 0 . The difference between subsequent spike times, J n = t n t n 1 , is called the interspike interval (ISI).

We will assume that τ is known (but see Sect. 6 for a discussion of the alternative) and nondimensionalize Eq. (1) as follows:

s = t τ , X s = V ( t ) v 0 v th v 0 , W s = W ( t ) τ , x th = 1 , α = μ τ v th v 0 , β = σ τ v th v 0 , γ = A τ v th v 0 , Ω = ω τ

to obtain

d X s = ( α X s + γ sin ( Ω s ) ) d s + β d W s , s 0 = 0 ; X s 0 = 0 , s n = inf { s > s n 1 : X s = x th = 1 } for  n 1 , { X s n + = 0 , I n = s n s n 1 ,
(2)

where we have defined I n = J n /τ. We can also write the dynamics between two spike times s n and s n + 1 in terms of elapsed time since the last spike, s =s s n , s < I n + 1 ,

d X s = ( α X s + γ sin ( Ω ( s + ϕ n ) ) ) d s + β d W s , { s = s s n , ϕ n = s n mod 2 π Ω .
(3)

This form of the dynamics highlights that this is not a renewal process since different trajectories between spikes have different phase shifts ϕ n = s n modulo 2π/Ω. This will be important in the following discussion. The shape of the ISI distribution depends on the model parameters, and it is natural to divide the parameter space in different regimes characterized by their qualitative behavior. Four distinct parameter regimes will be considered; supra-threshold, critical, subthreshold, and supersinusoidal. To understand the reasoning behind the regime names, observe that in the absence of noise, β=0, the deterministic model will produce spikes if and only if

α+ γ 1 + Ω 2 >1,

see the discussion in [14], which can be directly inferred from the solution in Eq. (12) below. In both the supra-threshold and supersinusoidal regimes, α+γ/ 1 + Ω 2 >1. The difference between the two is that in the supra-threshold regime the constant bias current alone is sufficient for spikes to occur, also in absence of noise, that is, α>1. In the supersinusoidal regime, the sinusoidal current is necessary for spikes to occur in absence of noise, that is, α+γ/ 1 + Ω 2 >1 and α1. In the critical regime, the sum of the two terms is just barely enough to guarantee deterministic spiking, that is α+γ/ 1 + Ω 2 1. Finally, in the subthreshold regime, there would be no spikes without the noise, α+γ/ 1 + Ω 2 <1.

Table 1 tabulates examples of corresponding parameter values for each regime, while Fig. 1 shows examples of individual voltage trajectories and their associated spike trains. Figures 2 and 3 illustrate how each regime behaves for selected ϕ’s by plotting the survivor distribution, G ¯ ϕ (t), and the probability density, g ϕ (t), both defined in Eq. (4) below.

Fig. 1
figure 1

Example trajectories from Eq. (2) for the four different parameter regimes using the parameter values given in Table 1. a supra-threshold, b supersinusoidal, c critical, d subthreshold. In the supra-threshold regime spikes occur regularly and often; in the supersinusoidal regime spikes cluster near the peak of the sine wave; in the critical regime they occur less often; and in the subthreshold regime, spikes occur rarely. For all regimes, Ω=1

Fig. 2
figure 2

The four different parameter regimes using the parameter values given in Table 1. Illustrated are the probability density functions, g ϕ m (t), for representative ϕ m =2π/Ω×{0,0.25,0.5,0.75}. Varying ϕ m has, for the most part, the effect of shifting the curves laterally, while varying α, β, γ changes their characteristic form. For all regimes, Ω=1. a supra-threshold, b supersinusoidal, c critical, d subthreshold

Fig. 3
figure 3

The four different parameter regimes using the parameter values given in Table 1. Illustrated are the survivor distribution functions, G ¯ ϕ m (t), for representative ϕ m =2π/Ω×{0,0.25,0.5,0.75}. Varying ϕ m has, for the most part, the effect of shifting the curves laterally, while varying α, β, γ changes their characteristic form. a supra-threshold, b supersinusoidal, c critical, d subthreshold

Table 1 Example of α, β, γ parameter values for the different regimes, given Ω=1

With regards to Figs. 2 and 3, it is worth noting explicitly that combinations of noise and sinusoidal forcing can cause firing patterns in which spikes are phase locked, but skip a certain number of cycles. This leads to multimodal ISI densities. There are many different dynamical mechanisms that can yield such patterns, and the particular correlations between the ISIs will depend on the underlying voltage dynamics (which, in our case, we assume to be given by Eq. (1)); in particular, it may be difficult to distinguish whether the dynamics are subthreshold or supra-threshold, since both can show similar ISI densities; see [22].

2.1 Basic ISI Probability Density Functions

Here, we introduce the notation for the probability density, distribution and survival functions of I n , an ISI arising from a trajectory produced by Eq. (3),

g ϕ ( τ ) d τ : = P ( I n + 1 [ τ , τ + d τ ) | ϕ n = ϕ ) (probability density) , G ϕ ( t ) : = P ( I n + 1 t | ϕ n = ϕ ) = 0 t g ϕ ( τ ) d τ (cumulative distribution) , G ¯ ϕ ( t ) : = P ( I n + 1 > t | ϕ n = ϕ ) = 1 G ϕ ( t ) (survivor distribution) .
(4)

The subscript ϕ is to stress that g, G, and G ¯ depend on the value of ϕ n in Eq. (3). This is the formal statement that in a sinusoidally-driven neuron, the interspike intervals are not identically distributed, and are only independent conditioned on the sinusoidal phase at an interval’s onset. Knowing these distributions would provide the likelihood function, offering estimation by the preferred method of choice, the maximum likelihood estimator. Unfortunately, explicit expressions for the ISI distribution are not available except for the special case of γ=0 and α=1; see [3]. Different representations of the likelihood function are available though, see [23], one of which we will use below.

2.2 Fokker–Planck Equation with Absorbing Boundaries

The Fokker–Planck equation is a partial differential equation (PDE) describing the evolution of the probability density, f(x,t), of X t . For the sinusoidally-forced Ornstein–Uhlenbeck process, Eq. (3), with the threshold x th =1, the PDE is

t f ( ϕ ) ( x , t ) = x [ ( α x + γ sin ( Ω ( t + ϕ ) ) ) f ( ϕ ) ] + x 2 [ β 2 2 f ( ϕ ) ] , x ( , 1 ) .
(5)

Due to the reset, we have that at time t=0, X t =0 and so for the initial conditions we can write

f ( ϕ ) (x,t=0)=δ(x),
(6)

where δ() is the Dirac delta function. The spike is represented as a zero boundary condition for f at x=1

f(1,t)=0.

The natural way of using the Fokker–Planck equation in first-hitting-times problems is as follows. Denote the integral of f ( ϕ ) by F ( ϕ ) (x,t)= ξ x f ( ϕ ) (ξ,t)dξ. F ( ϕ ) (x,t) can be related to the ISI’s survivor distribution function, G ¯ ϕ (t), by

G ¯ ϕ (t)= F ( ϕ ) (1,t).
(7)

Equation (7) forms the basis of one of the methods below for estimating the structural parameters from the observed data.

Since Eq. (5) has to be solved numerically, we will need to truncate its domain from below. The most natural way to do this, given the dynamics, is to impose reflecting boundary conditions at some x= x (αγ/ 1 + Ω 2 ) where the probability mass is very small. For the left (lower) limit of the computational domain, we use the formula

x =min ( α γ / 1 + Ω 2 mean 2 β / 2 s t d . d e v . , 0.25 ) .

This choice requires some explanation. In the t limit, the distribution of X t in Eq. (3) without thresholding is Gaussian with mean given by Eq. (12) (below) and variance equal to β 2 /2. Thus, to truncate the computational domain for the thresholded process from below, we take the lowest value of the asymptotic mean, αγ/ 1 + Ω 2 , then from this we subtract two standard deviations, 2β/ 2 and set the result to be the lower bound, x . Finally, if this value for x happens to be larger than −0.25, we enforce that x 0.25.

Numerical considerations lead us to solve for F, instead of f, since delta functions are difficult to represent in floating point, while the initial conditions for F, the Heaviside step function, H(x), faces no such difficulties [24]. The Heaviside step function is defined to be equal to 0 for x<0 and to be equal to 1 for x0. At this point, we need to derive the PDE for the distribution F, starting from the PDE for the density, f, Eq. (5).

First, at the lower boundary, it is intuitive that the distribution should be zero, F( x ,t)=0, while f(1,t)=0 implies that at the upper boundary x F(1,t)=0. Inside the domain, the PDE itself reformulates as

t f(x,t)= x [ 1 2 x [ β 2 f ] ( α x + γ sin ( Ω ( t ϕ ) ) ) f ]

so that

x t F(x,t)= x [ β 2 2 x 2 F ( α x + γ sin ( Ω ( t + ϕ ) ) ) x F ] .

Integrating with respect to x then gives

t F(x,t)= β 2 2 x 2 F ( α x + γ sin ( Ω ( t + ϕ ) ) ) x F+C(t),

where C(t) is a constant of integration depending on t. Now consider the lower boundary condition, x= x . Here, F( x ,t)=0 implies that t F=0 and so

C(t)= [ β 2 2 x 2 F ( α x + γ sin ( Ω ( t + ϕ ) ) ) x F ] .
(8)

The right-hand side in Eq. (8) is precisely the reflecting boundary condition on f once we recall that x F=f. Therefore, C(t)0.

Thus, the fully specified PDE for F, which we will be solving frequently in what follows, is

t F ( ϕ ) ( x , t ) = β 2 2 x 2 F ( ϕ ) ( α x + γ sin ( Ω ( t + ϕ ) ) ) x F ( ϕ ) , { F ( ϕ ) ( x , 0 ) = H ( x ) , F ( ϕ ) ( x , t ) | x = x 0 , x F ( ϕ ) ( x , t ) | x = 1 0 .
(9)

Numerical solutions for Eq. (9) are shown in Fig. 4. We have used the standard Crank–Nicholson finite-difference algorithm (central-differences in space with equally weighted implicit-explicit terms in time, see [25]).

Fig. 4
figure 4

Example solution to Eq. (9) for (α,β,γ)=(0.5,0.3,0.5 2 ); Ω=1, ϕ=π/2. In a, b, c, we show the full solution in space–time F(x,t). In d we show the time solution at the upper boundary, F(1,t)

2.3 Fortet Equation

Consider a general form of Eq. (3),

d Y t =b(t, Y t )dt+σ(t, Y t )d W t .

Let Φ(y,t| y 0 , t 0 ):=P[ Y t y| Y t 0 = y 0 ] be the transition cumulative distribution of Y. Note that this is the distribution of Y t in absence of a threshold, different from the distribution given in Eq. (7), which is the distribution of the process constrained to be below the threshold. Now consider an arbitrary time-dependent threshold v th (t). The Fortet equation (see [26]) convolves the first-hitting time probabilities, g(t), with the transition density of the process. Integrating over (, v th (t)), we obtain

1Φ ( v th ( t ) , t | v 0 , 0 ) = 0 t g(τ) [ 1 Φ ( v th ( t ) , t | v th ( τ ) , τ ) ] dτ.
(10)

The left-hand side is simply the probability of exceeding v th at time t starting at v 0 at time 0. This can also be written as the probability of hitting v th for the first time at time τ<t and then exceeding v th at time t starting at v th at time τ, integrated over all τ.

The Fortet equation is particularly appealing to use when we have an analytical expression for Φ. For the problem at hand, Φ is complicated due to the time-dependent forcing. However, the following transformation yields a time-homogeneous Y for which Φ will be tractable, along with an associated moving threshold, v th (t). This makes feasible the use of the Fortet equation. To obtain this transformation, cf. [27], consider the deterministic version of the SDE in Eq. (3)

d v ( t ) = ( α v + γ sin ( Ω ( t + ϕ ) ) ) d t , v ( 0 ) = 0
(11)

with solution

v ( t ) = α ( 1 exp ( t ) ) + γ 1 + Ω 2 [ sin ( Ω ( t + ϕ ) ψ ) exp ( t ) sin ( ϕ Ω ψ ) ] ; ψ = arctan ( Ω ) .
(12)

Now take X t , the solution to Eq. (3) and v(t), Eq. (12), and let Y t = X t v(t). Then

d Y t = Y t dt+βdW,
(13)

which has the time and parameter dependent threshold

v th { α , γ ; ϕ } (t)= v th v(t).
(14)

That is, X t hits the constant threshold v th if and only if Y t hits the moving threshold v th { α , γ ; ϕ } (t), where the subindex indicates the dependence on α, γ and ϕ. Therefore, the ISIs produced by X and Y are the same and so are their distributions. Thus, g ϕ (τ) satisfies

1 Φ { β } ( v th { α , γ ; ϕ } ( t ) , t | 0 , 0 ) = 0 t g ϕ ( τ ) [ 1 Φ { β } ( v th { α , γ ; ϕ } ( t ) , t | v th { α , γ ; ϕ } ( τ ) , τ ) ] d τ ,
(15)

where

Φ { β } (y,t| y 0 , t 0 )= 1 π β 2 ( 1 e 2 ( t t 0 ) ) y exp ( ( x y 0 e ( t t 0 ) ) 2 β 2 ( 1 e 2 ( t t 0 ) ) ) dx

is the conditional cumulative distribution function of Y t defined in Eq. (13).

3 Parameter Estimation Algorithms

The unknown parameters in Eq. (3) are α, β, and γ, while we assume Ω known. The reason why the amplitude, γ, is often unknown while the frequency, Ω, is known is that one can usually observe the sinusoidal input and thus its frequency. Further, the encoding of the input into neuronal firing patterns often involves phase locking to the sinusoidal component. However, the actual forcing amplitude at the level of the neuron is usually modified by various synaptic and other filtering processes, unless the cell receives direct sinusoidal current injection.

Our goal is to estimate the structural parameters (α, β, γ) from a sample of spike time data, { i 1 ,, i N }. There are several algorithms for estimating the parameters for the simpler and more common case of γ=0. One such algorithm relies on the Fortet equation (see [7, 8]), which we extend to the presence of a time-varying current. A more basic approach is to directly solve the Fokker–Planck equation for the probability density of X t , [1921], from which one can derive the survival distribution of I n , and use this to compare against the empirical survival distribution of I n obtained from data. An approximate maximum likelihood approach is also possible by numerical differentiation. The relation between Fokker–Planck equations and the first-passage time problem is discussed in most introductory books on stochastic analysis; see, for example, [28]. A recent review of this approach for the simple γ=0 case in neuronal modeling can be found in [21], wherein the first passage problem is discussed at great lengths in the context of spiking neurons. We will use this in Sect. 2.2. A more elaborate approach using the Fokker–Planck equation to approximate the hitting time distribution is given in [29]. The techniques in [29] avoid the need to compute the Fokker–Planck PDE numerically, instead approximating it with analytically known solutions. This approach might offer significant computational savings, but since this would at most amount to a computational speed-up of our algorithm, we have left this unexplored for now.

The immediate problem in generalizing the aforementioned approaches to the case of γ0 is that the I n ’s are no longer identically distributed since the phase ϕ n 1 of the n th interval I n depends on t n 1 , the time the previous spike occurred. The I n ’s are also dependent, but conditionally independent given ϕ n 1 . So the trajectories in each interval are parameterized by the value of ϕ n 1 at the time of the last spike/reset. We overcome this obstacle by splitting the I n ’s in groups, and approximating the I n ’s within groups as coming from identically distributed trajectories in a sense to be specified below. This approximation which solves the challenge of dependent and non-identically distributed ISIs is the primary contribution of this paper.

3.1 ϕ-Binning

Before we can use Eq. (9) or (15), we need to deal with the fact that ϕ is not fixed, but instead each I n starts with a distinct ϕ n . Our approach is to partition the interval [0,2π/Ω] into M bins, where MN, and represent each bin by the midpoint of the bin, ϕ m . Then we approximate the N observed ϕ n ’s by the closest ϕ m and pretend that any observed I n was not produced by a trajectory of the form in Eq. (3) with ϕ= ϕ n , but with ϕ= ϕ m . Our hope is that for a judicious choice of M, we can balance the error of ϕ n ϕ m with having enough data points in each bin in order to obtain a useful estimate from Eq. (9) or (15).

There is clearly much freedom in how one sets up these bins, but we will do the simplest thing and make them all of equal width, δϕ=2π/(ΩM). Each ϕ n will belong to one and only one of the bins [ ϕ m δ ϕ / 2 , ϕ m + δ ϕ / 2 ) m = 1 M , with center points ϕ m =δϕ/2+(m1)δϕ, for m=1,,M. Thus, given an empirically observed I n with associated ϕ n , we will pretend that it was produced by the process

d X s =(α X s )ds+γsin ( Ω ( s + ϕ m ( n ) ) ) ds+βd W s ,

where

ϕ m (n)= arg min ϕ m | ϕ n ϕ m |.

This binning is illustrated in Fig. 5.

Fig. 5
figure 5

The raw ( i n , ϕ n ) pairs (a) are binned into a set of M bins with a representative ϕ m (b) and the ISIs within each bin are treated as a renewal process. In this illustration, M=8, Ω=1 while the parameters α, β, γ are taken from the supra-threshold regime

While we have no rigorous approach to determine the value of M, our limited experience suggests that given N=1000 ISIs, M=10, or M=20 gives satisfactory results for very different parameter regimes. In general, choosing M is a balancing act. For M too high, the resulting bins will have too few data points to approximate G ¯ (I) accurately. Therefore, M is forced to be small when sample size is not large. For M too low, the approximation of the phase shifts will be poor, leading to a biased estimate of G ¯ (I). We illustrate the effect of increasing M in Fig. 6. Generally, as long as there are sufficient data points, as M increases, the approximation of using the survival distribution with ϕ m instead of ϕ n improves since ϕ m (n) ϕ n as M. In the sequel, we will use M=20 for sample sizes of N=1000 and M=8 for sample sizes of N=100.

Fig. 6
figure 6

Effect of M, the number of bins, on the approximate survival distribution. The full-drawn blue curve is the true survivor distribution given in Eq. (9), the red points are the approximation given in Eq. (16). In the figures, the least populous (above) and most populous (below) bin for each M is shown. The width of the bins is δϕ=2π/(ΩM). We have used a, eM=5; b, fM=10; c, gM=20; d, iM=40. As M increases, the approximation of using the survival distribution using ϕ m instead of ϕ n improves since ϕ m (n) ϕ n as M. The data are generated using parameter values from the supersinusoidal regime and N=1000. For this particular data set the largest generated ISI was 6.55 time units

3.2 Fokker–Planck Algorithm

Within each bin it is clear how to apply Eq. (7). In the m th bin, for a given ϕ m , we approximate G ¯ ϕ (t) by

G ¯ ˆ ϕ m (t)= # [ i n > t | ϕ n 1 [ ϕ m δ ϕ / 2 , ϕ m + δ ϕ / 2 ) ] N m ,
(16)

where N m is the number of ISIs in bin m. Using Eq. (7), we define the loss function

L(α,β,γ)= ϕ m N m { sup t > 0 | G ¯ ˆ ϕ m ( t ) F α , β , γ ϕ m ( x th , t ) | } .
(17)

The weight N m is included so that bins with larger sample sizes have a larger influence on the estimates.

To evaluate the supremum in Eq. (17), we spline interpolate the empirically discrete G ¯ ˆ for each ϕ m , sample at the time nodes of the PDE discretization and finally take the maximum amongst the sampled values. We then minimize L using an optimization algorithm (see below, Sect. 4) and take our estimates α ˆ , β ˆ , γ ˆ to be

α ˆ , β ˆ , γ ˆ = arg min α , β , γ L(α,β,γ).

Note that the relation between the spike time survival density, G ¯ ϕ and the transition distribution, F ϕ , in Eq. (7) could also allow for an approximate maximum likelihood estimator (MLE), based on maximizing

L MLE (α,β,γ)= n log ( g ϕ n 1 ( i n ) ) = n log [ t F α , β , γ ϕ n 1 ( x th , t ) ] | t = i n ,

where the derivative has to be approximated by finite differences. We can then again use binning to avoid having to compute the PDE separately for each ( i n , ϕ n 1 ). Our experience with the MLE approach has been that the quality of the estimates provided are similar to those obtained by minimizing Eq. (17) and that the associated computing times are on the same order. Due to this similarity and in order to keep the paper concise, we include details of the MLE estimates only in the supplementary online material.

3.3 Fortet Algorithm

An alternative approach is to form a loss function from Eq. (15). This is similar to what is done in [7, 8] for the simpler case of a constant threshold. Noting that 0 t g(τ)[1Φ]dτ=E[(1Φ) 1 I t ] where the expectation is taken with respect to the distribution of the random variable I, we can use the fact that the ISIs are approximately independent and invoke the law of large numbers to estimate the integral as

0 t g ϕ m ( τ ) [ 1 Φ { β } ( ϕ ) ( v th { α , γ ; ϕ } ( t ) , t | v th { α , γ ; ϕ } ( τ ) , τ ) ] d τ 1 N m i n < t [ 1 Φ { β } ( ϕ ) ( v th { α , γ ; ϕ } ( t ) , t | v th { α , γ ; ϕ } ( i n ) , i n ) ] .

We then define the loss function

L ( α , β , γ ) = ϕ m N m { sup s > 0 | 1 Φ { β } ( ϕ m ) ( v th { α , γ ; ϕ } ( s ) , s | 0 , 0 ) 1 N m i n < s [ 1 Φ { β } ( ϕ m ) ( v th { α , γ ; ϕ } ( s ) , s | v th { α , γ ; ϕ } ( i n ) , i n ) ] | / ω ( ϕ m ; α , β , γ ) } .
(18)

We divide each inner term by ω( ϕ m ;α,β,γ)= sup s > 0 |1 Φ α , β , γ ( ϕ m ) ( v th (s),s| v 0 )|, following the suggestion in [8]. This scaling ensures that Eq. (15) divided by ω(α,β,γ) will vary between 0 and 1 for all parameter values thus giving sense to the measure defined by the loss function. Since we can solve in closed form for Φ, we have all we need given an observed spike train of i n ’s. We evaluate the sup by sampling at K=500 uniformly spaced points in (0, I max +ϵ] and taking the maximum of the sampled values.

3.4 Initialization of the Algorithms

The parameter search can be initialized in a simple way using the fact that the Fokker–Planck PDE is almost an “advection-diffusion” equation whose solution is almost a Gaussian. Then G ¯ (t) can be approximated by the amount of probability mass of a Gaussian to the left of the threshold at time t. The idea is as follows. Suppose we are solving the following PDE:

t ρ=U x [ρ]+ β 2 2 x 2 [ρ].
(19)

Its solution given an initial condition ρ(x,0)=δ(x) will be a Gaussian bell moving to the right with speed U and standard deviation σ=β t .

The survivor function G ¯ (t) can be thought of as the amount of area that has passed the threshold (from the left moving to the right). We can then invert the information about G ¯ to estimate U and β. In particular, a Gaussian bell has ≈0.158 of its mass more than one standard deviation to the right of its mean. Thus, at time t 1 such that G ¯ ( t 1 )=0.842, the right tail of more than one standard deviation of the Gaussian bell has crossed the threshold. The threshold is at x=1 and we obtain the following equation:

U t 1 +β t 1 =1.
(20)

Similarly, at time t 2 such that G ¯ ( t 2 )=10.842, the Gaussian bell has crossed the threshold except for the left tail and we have

U t 2 β t 2 =1.
(21)

If U and β were constant, then Eqs. (20) and (21) provide two equations in two unknowns. However, U=U(x,t)=(αx+γsin(Ω(t+ϕ))) is not constant and we approximate U as

U(x,t)α0.5+γ 1 t 0 t sin ( Ω ( τ + ϕ ) ) dτ,
(22)

i.e., we approximate the space-dependent term, x, with the mid-point between the reset value, v 0 =0, and the threshold, v th =1, and we approximate the time-dependent term, sin(Ω(τ+ϕ)), by its time-average value between 0 and t. If we use the 0th, 1st, and 2nd standard deviation points, we can form 5 equations in 3 unknowns as follows:

α t 1 + γ s ( t 1 ) + 2 β t 1 = 1 + 0.5 t 1 , α t 2 + γ s ( t 2 ) + β t 2 = 1 + 0.5 t 2 , α t 3 + γ s ( t 3 ) + 0 β = 1 + 0.5 t 3 , α t 4 + γ s ( t 4 ) 1 β t 4 = 1 + 0.5 t 4 , α t 5 + γ s ( t 5 ) 2 β t 5 = 1 + 0.5 t 5

with the time-average weighting function s(t)=(cos(Ωϕ)cos(Ω(t+ϕ)))/Ω. However, the approximation is best for earlier times, when the solution is closer to a Gaussian bell that is approaching the threshold, but less correct for later times, since it neglects the loss of probability mass and thus overestimates the backward probability current. Indeed, we have found it to be best to use only t 1 and t 2 . In the following we use only these equations:

α t 1 + γ s ( t 1 ) + 2 β t 1 = 1 + 0.5 t 1 , α t 2 + γ s ( t 2 ) + β t 2 = 1 + 0.5 t 2

for the initializer. We can form these equations separately for each ϕ m bin, thus resulting in M×2 equations for the unknowns α, β, and γ. Since we have more equations than unknowns, we use least-squares estimates in a regression to pick out unique α, β, and γ estimates.

The proposed initialization procedure has two advantages. First, it is automatic, i.e., it requires only the data and no input or guidance from the user. Second, it is extremely fast. While the precise effect of the initializer is shown in Sect. 4, it is intuitively clear that it will work best in the supra-threshold parameter regime when the bell curve is truly moving past the threshold as a whole and less so for subthreshold regimes, when only the diffusive force serves to propel the process to reach v th . The behavior of the initializer in the different regimes is illustrated in Fig. 7. What we show in Fig. 7 is the following: First, we show the survival distribution for a given regime and ϕ m fixed. Then using data generated from such a regime and with ϕ n in the m th bin, the initializer tries to find the best approximation by the motion of a Gaussian bell which will fit these data, in the sense of solving for α, β, γ as previously described. Once this is done, we then show in red the amount of area under this Gaussian bell to the left of the threshold. Of course, the interpretation of the survival distribution for an ISI as a fraction of the area under a moving bell with conserved total area is wrong, but the assumption is useful in automatically generating initial values for the more appropriate approximations to start their work.

Fig. 7
figure 7

The blue curves are the numerically obtained survivor distributions G ¯ ϕ for the exact parameters in the four regimes (as in Table 1) and Ω=1. The red curves are obtained in the following manner: Simulations using the true parameters were used to generate sample spikes. Using these samples, the initializer algorithm was used to generate estimates for α, β, γ. Using these estimates, the bell curve discussed in Sect. 3.4 was formed and evolved in time. Thus, the red curve drawn in the figures measures the area under this bell that is to the left of the threshold at time t. a supra-threshold, b supersinusoidal, c critical, d subthreshold

4 Method Comparison on Simulated Data

We will now use our algorithms on spike trains simulated from the four different regimes: the supra-threshold, the critical, the supersinusoidal and the subthreshold. We have used 100 sample spike trains per regime, with N=100 as well as N=1000 spikes per train. In order to perform the numerical minimization of Eqs. (17) and (18), we have used an implementation of the Nelder–Mead algorithm from the SciPy library [30]. The Nelder–Mead algorithm is a non-linear minimization routine which uses a bounding-polygon method to zero-in on the minimum and thus avoids the need to provide the gradient of the loss function. It is the standard non-gradient minimization algorithm.

The estimation results are shown in Figs. 8, 9, 10, 11, where we plot box plots for the estimates, α ˆ , β ˆ , γ ˆ in the four regimes. We also tabulate the average and the empirical 95 % confidence intervals of the estimates in Tables 2 and 3. Conclusions that can be drawn from these results are as follows. The initializer method is effective for the supra-threshold regime and gives reasonable ballpark estimates for all regimes, though the error can be substantial for the supersinusoidal regime. In general, both the Fortet and Fokker–Planck algorithm estimate the parameters well in the supra-threshold, critical and supersinusoidal regimes. The estimators’ variance is especially low in the supra-threshold regime, while it is higher for the critical and supersinusoidal regimes. In the supersinusoidal regime, the two algorithms give accurate estimates even though the initializer can be quite off. On the other hand, in the subthreshold regime, the initializer has a performance comparable to that of the two more involved methods. It seems that distinguishing between the constant bias and the sinusoidal current is difficult if their sum is not sufficient to generate spikes without noise.

Fig. 8
figure 8

Boxplots of parameter estimates for the supra-threshold regime. The upper plots (a, b, c) show estimates using N=100 sample spikes per estimation, while the lower plots (d, e, f) use N=1000. The dashed line indicates the true parameter value, while the red line inside the boxes indicates the median of the estimates. The boxes contain the central 50 % of the estimates. The bars indicate the range of the estimates, except for outliers given by the points outside the bars, and defined to be more than 1.5 times the interquantile range (the height of the box) from the box

Fig. 9
figure 9

Boxplots of parameter estimates for the supersinusoidal regime. The upper plots (a, b, c) show estimates using N=100 sample spikes per estimation, while the lower plots (d, e, f) use N=1000. The dashed line indicates the true parameter value, while the red line inside the boxes indicates the median of the estimates. The boxes contain the central 50 % of the estimates. The bars indicate the range of the estimates, except for outliers given by the points outside the bars, and defined to be more than 1.5 times the interquantile range (the height of the box) from the box

Fig. 10
figure 10

Boxplots of parameter estimates for the critical regime. The upper plots (a, b, c) show estimates using N=100 sample spikes per estimation, while the lower plots (d, e, f) use N=1000. The dashed line indicates the true parameter value, while the red line inside the boxes indicates the median of the estimates. The boxes contain the central 50 % of the estimates. The bars indicate the range of the estimates, except for outliers given by the points outside the bars, and defined to be more than 1.5 times the interquantile range (the height of the box) from the box

Fig. 11
figure 11

Boxplots of parameter estimates for the subthreshold regime. The upper plots (a, b, c) show estimates using N=100 sample spikes per estimation, while the lower plots (d, e, f) use N=1000. The dashed line indicates the true parameter value, while the red line inside the boxes indicates the median of the estimates. The boxes contain the central 50 % of the estimates. The bars indicate the range of the estimates, except for outliers given by the points outside the bars, and defined to be more than 1.5 times the interquantile range (the height of the box) from the box

Table 2 Averages and empirical 95 % confidence intervals of the estimates for N=100 spikes per train
Table 3 Averages and empirical 95 % confidence intervals of the estimates for N=1000 spikes per train

The Fokker–Planck method has a larger bias but a smaller spread than the Fortet method for N=100, Table 2. However for N=1000, the two methods have comparable spreads, while the Fortet method retains a smaller bias, see Table 3. More precisely, for N=1000, the Fokker–Planck method has a smaller spread in the subthreshold regime, while the Fortet method has a smaller spread in the supersinusoidal regime. As such, at least in the supersinusoidal regime, the Fortet method seems superior. A detailed comparison between the Fortet and the Fokker–Planck estimators for each parameter in each regime can be seen in Fig. 12 for N=100 and in Fig. 13 for N=1000.

Fig. 12
figure 12

Estimates based on samples of N=100 spikes obtained from the Fokker–Planck algorithm against the Fortet algorithm for the four different parameter regimes, with parameter values given in Table 1, fixing Ω=1. Each row corresponds to one regime and one set of simulations. Each column corresponds to a parameter, with the specific value indicated above each plot. a, b, c supra-threshold; d, e, f supersinusoidal; g, h, i critical; j, k, l subthreshold

Fig. 13
figure 13

Estimates based on samples of N=1000 spikes obtained from the Fokker–Planck algorithm against the Fortet algorithm for the four different parameter regimes, with parameter values given in Table 1, fixing Ω=1. Each row corresponds to one regime and one set of simulations. Each column corresponds to a parameter, with the specific value indicated above each plot. a, b, c supra-threshold; d, e, f supersinusoidal; g, h, i critical; j, k, l subthreshold

The two algorithms are numerically intensive. For N=100 and N=1000 spikes, we show the times taken for the estimation in Table 4. While we have done most of our numerical work in Python/SciPy [30], we have implemented the critical components of both algorithms in C. That is, we solve the inner part of Eq. (18) and the Fokker–Planck PDE, Eq. (9), in C using the GSL libraries [31]. From Table 4, we can verify that the computing time for the Fortet algorithm scales proportionally with the number of spikes. This is to be expected, since the Fortet equation has a term of the form i n which in turn has N terms and this forms the bulk of the computing time for the Fortet equation. The Fokker–Planck algorithm, on the other hand, scales less-than-linearly with N, since the dependency on N is in forming the approximation, G ¯ ˆ to the survivor function and that is not computationally intensive (solving the PDE is).

Table 4 Average times ± standard deviations in seconds for the algorithm in various regimes. Left: N=100 spikes; right: N=1000 spikes

5 The Effect of Ω

So far, we have held Ω constant and equal to 1. We now investigate the effect of varying Ω on the quality of estimates. To narrow the scope, we focus on increasing Ω while keeping the parameters in the critical regime such that α+γ/ 1 + Ω 2 =1 and α=0.5. This amounts to increasing γ with Ω. We do the estimations for four values of Ω=[1,5,10,20]. Similarly to the previous section, we use 100 sample spike trains per parameter set, with each spike train consisting of N=1000 ISIs.

We show box plots of the estimates for each Ω in Fig. 14. We then directly compare the two algorithms, Fortet vs. Fokker–Planck, in Fig. 15. The immediate observation is that the Fokker–Planck algorithm fails to keep up at the higher frequencies and consistently underestimates γ. The Fortet algorithm does better, but still underestimates γ. In general, this underestimation of γ is accompanied by an overestimation of α. This is exacerbated at higher Ω. We illustrate the relation between estimates for α vs. γ in Fig. 16, where it is quite clear that an underestimation of γ is proportional to the overestimation of α.

Fig. 14
figure 14

Boxplots of parameter estimates for varying Ω across [1,5,10,20] while holding γ/ 1 + Ω 2 constant as to keep the parameters in the critical regime. acΩ=1, dfΩ=5, giΩ=10, jlΩ=20. The boxes contain the central 50 % of the estimates. The bars indicate the range of the estimates, except for outliers given by the points outside the bars, and defined to be more than 1.5 times the interquantile range (the height of the box) from the box

Fig. 15
figure 15

Estimates based on samples of N=1000 spikes obtained from the Fokker–Planck algorithm against the Fortet algorithm for a parameter set in the critical regime, while varying Ω across [1,5,10,20] and holding γ/ 1 + Ω 2 and α constant. a, b, cΩ=1; d, e, fΩ=5; g, h, iΩ=10; j, k, lΩ=20

Fig. 16
figure 16

Comparison of α ˆ vs. γ ˆ parameter estimates while varying Ω across [1,5,10,20], holding γ/ 1 + Ω 2 constant as to keep the parameters in the critical regime. a, b, cΩ=1; d, e, fΩ=5; g, h, iΩ=10; j, k, lΩ=20

For completeness, we also include the estimates’ average and empirical 95 % confidence intervals in Table 5.

Table 5 Averages and empirical 95 % confidence intervals of estimates for N=1000 spikes per train in the critical regime for varying Ω across [1,5,10,20]. Note that the upper subtable corresponds to the third subtable in Table 3; numbers differ slightly due to statistical fluctuations in the simulations

6 Discussion and Outlook

We have shown two methods to estimate parameters in Eq. (2) from ISI data. Our methods are based on binning the spikes into bins with representative phase shifts. We have devised a constructive procedure to automatically initialize the methods from the data.

Our computational results suggest that for low frequencies the Fortet algorithm is superior for large sample sizes, especially in the supersinusoidal regime, while the Fokker–Planck algorithm has a comparable accuracy and a lower variance for small sample sizes. Both algorithms find sensible estimates most of the time, although they seem less effective in the subthreshold regime. Their performance can be partially attributed to the ability of the initializer algorithm to supply good guesses for starting the optimization iterations.

The Fokker–Planck equation allows for approximate maximum likelihood estimation. We chose an alternative loss function, though, because it marginally appeared more robust, possibly because a numerical derivation step is avoided. This is further investigated by simulations in the supplementary online material. The simulations suggest that the distribution of the maximum likelihood estimates in the supersinusoidal regime appears bimodal, which is not the case for the alternative loss function, Eq. (17).

We have also made a preliminary exploration of the effect of Ω on the quality of the estimates. Our results show that an increase in Ω makes the parameters α and γ more difficult to estimate accurately and at high Ω, γ is underestimated, while α is over-estimated. We find that in this scenario, the Fortet algorithm does a markedly more accurate job then the Fokker–Planck algorithm.

We have assumed the time-constant τ of the leak term to be known. In most experiments that is not realistic, and it would be preferable to estimate τ alongside the other parameters. However, it is difficult to estimate [32]. When we tried to estimate it together with the other parameters, we usually obtained results which were not accurate. The obtained estimates resulted in ISIs that very well matched the data, no worse than the ISIs obtained from the true parameters. This leads us to believe that the simultaneous estimation of τ along with α, β, γ using only ISI data suffers from identifiability problems. In [5], they were able to estimate τ in the simpler nonsinusoidally-driven model, but concluded that adding τ as an unknown dramatically reduced the accuracy in the estimation of the other unknown parameters. The reason is that if τ is also estimated from a single dataset alongside the other parameters, then a reasonable fit can be found to the data for various combinations of α, β, γ, and τ, but the so-obtained parameter values can be far from the true values.

Our model is relatively simple and ignores neurophysiological realism, such as the fact that the spiking threshold is likely nonconstant, with a time-dependent functional form that would involve further unknown parameters. A recent paper attempting the parameter estimation in such a model, but without sinusoidal forcing, is [20]. Furthermore, intracellular recordings suggest that a hard threshold is a rough approximation and an exponential voltage-dependent spiking intensity is more realistic [33].

While our work has used a very specific form of the periodic forcing term, namely γsin(Ωt), it is clear how to apply the approach to an arbitrary periodic function. This can be done as long as one knows where in the period of oscillation a spike has occurred. If that is the case, then the binning procedure can be applied and the estimation methods proposed can be attempted.

References

  1. Brillinger DR: Maximum likelihood analysis of spike trains of interacting nerve cells. Biol Cybern 1988, 59: 189–200. 10.1007/BF00318010

    Article  MATH  Google Scholar 

  2. Inoue J, Sato S, Ricciardi L: On the parameter estimation for diffusion models of single neuron’s activities. Biol Cybern 1995, 73: 209–221. 10.1007/BF00201423

    Article  MATH  Google Scholar 

  3. Ditlevsen S, Lánský P: Estimation of the input parameters in the Ornstein–Uhlenbeck neuronal model. Phys Rev E 2005., 71: Article ID 011907 Article ID 011907

    Google Scholar 

  4. Ditlevsen S, Lánský P: Estimation of the input parameters in the Feller neuronal model. Phys Rev E 2006., 73: Article ID 061910 Article ID 061910

    Google Scholar 

  5. Mullowney P, Iyengar S: Parameter estimation for a leaky integrate-and-fire neuronal model from ISI data. J Comput Neurosci 2008, 24(2):179–194. 10.1007/s10827-007-0047-5

    Article  MathSciNet  Google Scholar 

  6. Zhang X, You G, Chen T, Feng J: Maximum likelihood decoding of neuronal inputs from an interspike interval distribution. Neural Comput 2009, 21(11):3079–3105. 10.1162/neco.2009.06-08-807

    Article  MATH  MathSciNet  Google Scholar 

  7. Ditlevsen S, Ditlevsen O: Parameter estimation from observations of first-passage times of the Ornstein–Uhlenbeck process and the Feller process. Probab Eng Mech 2008, 23: 170–179. 10.1016/j.probengmech.2007.12.024

    Article  Google Scholar 

  8. Ditlevsen S, Lánský P: Parameters of stochastic diffusion processes estimated from observations of first-hitting times: application to the leaky integrate-and-fire neuronal model. Phys Rev E 2007., 76: Article ID 041906 Article ID 041906

    Google Scholar 

  9. Ditlevsen S, Lánský P: Comparison of statistical methods for estimation of the input parameters in the Ornstein–Uhlenbeck neuronal model from first-passage times data. American Institute of Physics Proceedings Series CP1028. In Collective Dynamics: Topics on Competition and Cooperation in the Biosciences. Edited by: Ricciardi LM, Buonocore A, Pirozzi E. Am Inst of Phys, Melville; 2008:171–185.

    Google Scholar 

  10. Lánský P, Ditlevsen S: A review of the methods for signal estimation in stochastic diffusion leaky integrate-and-fire neuronal models. Biol Cybern 2008, 99(4–5):253–262. 10.1007/s00422-008-0237-x

    Article  MATH  Google Scholar 

  11. Braun HA, Wissing H, Schafer K, Hirsch MC: Oscillation and noise determine signal-transduction in shark multimodal sensory cells. Nature 1994, 367(6460):270–273. 10.1038/367270a0

    Article  Google Scholar 

  12. Chacron MJ, Longtin A, St-Hilaire M, Maler L: Suprathreshold stochastic firing dynamics with memory in P-type electroreceptors. Phys Rev Lett 2000, 85(7):1576–1579. 10.1103/PhysRevLett.85.1576

    Article  Google Scholar 

  13. Bulsara AR, Elston TC, Doering CR, Lowen SB, Lindenberg K: Cooperative behavior in periodically driven noisy integrate-fire models of neuronal dynamics. Phys Rev E 1996, 53(4):3958–3969. 10.1103/PhysRevE.53.3958

    Article  Google Scholar 

  14. Burkitt AN: A review of the integrate-and-fire neuron model: II. Inhomogeneous synaptic input and network properties. Biol Cybern 2006, 95(2):97–112. 10.1007/s00422-006-0082-8

    Article  MATH  MathSciNet  Google Scholar 

  15. Lánský P: Sources of periodical force in noisy integrate-and-fire models of neuronal dynamics. Phys Rev E 1997, 55(2):2040–2043. 10.1103/PhysRevE.55.2040

    Article  Google Scholar 

  16. Longtin A, Bulsara A, Pierson D, Moss F: Bistability and the dynamics of periodically forced sensory neurons. Biol Cybern 1994, 70(6):569–578. 10.1007/BF00198810

    Article  MATH  Google Scholar 

  17. Sacerdote L, Giraudo MT: Stochastic integrate and fire models: a review on mathematical methods and their applications. Lecture Notes in Mathematics 2058. In Stochastic Biomathematical Models with Applications to Neuronal Modeling. Edited by: Bachar M, Batzel J, Ditlevsen S. Springer, Berlin; 2013.

    Google Scholar 

  18. Shimokawa T, Pakdaman K, Takahata T, Tanabe S, Sato S: A first-passage-time analysis of the periodically forced noisy leaky integrate-and-fire model. Biol Cybern 2000, 83(4):327–340. 10.1007/s004220000156

    Article  MATH  Google Scholar 

  19. Paninski L, Pillow J, Simoncelli E: Maximum likelihood estimation of a stochastic integrate-and-fire neural encoding model. Neural Comput 2004, 16(12):2533–2561. 10.1162/0899766042321797

    Article  MATH  Google Scholar 

  20. Dong Y, Mihalas S, Russell A, Etienne-Cummings R, Niebur E: Estimating parameters of generalized integrate-and-fire neurons from the maximum likelihood of spike trains. Neural Comput 2011, 23(11):2833–2867. 10.1162/NECO_a_00196

    Article  MATH  Google Scholar 

  21. Sirovich L, Knight B: Spiking neurons and the first passage problem. Neural Comput 2011, 23(7):1675–1703. 10.1162/NECO_a_00139

    Article  MATH  MathSciNet  Google Scholar 

  22. Longtin A: Mechanisms of stochastic phase locking. Chaos, Interdiscip J Nonlinear Sci 1995, 1(5):209–215.

    Article  Google Scholar 

  23. Alili L, Patie P, Pedersen J: Representations of the first hitting time density of an Ornstein–Uhlenbeck process. Stoch Models 2005, 21: 967–980. 10.1080/15326340500294702

    Article  MATH  MathSciNet  Google Scholar 

  24. Hurn A, Jeisman J, Lindsay K: ML estimation of the parameters of SDE’s by numerical solution of the Fokker–Planck equation. MODSIM 2005: International Congress on Modelling and Simulation: Advances and Applications for Management and Decision Making 2005, 849–855.

    Google Scholar 

  25. Karniadakis GE, Kirby RM: Parallel Scientific Computing in C++ and MPI. Cambridge University Press, Cambridge; 2003.

    Book  MATH  Google Scholar 

  26. Fortet R: Les fonctions aléatoires du type de Markoff associées à certaines équations linéaires aux dérivées partielles du type parabolique. J Math Pures Appl (9) 1943, 22: 177–243.

    MATH  MathSciNet  Google Scholar 

  27. Shimokawa T, Pakdaman K, Sato S: Time-scale matching in the response of a leaky integrate-and-fire neuron model to periodic stimulus with additive noise. Phys Rev E 1999, 59(3):3427–3443. 10.1103/PhysRevE.59.3427

    Article  Google Scholar 

  28. Jacobs K: Stochastic Processes for Physicists: Understanding Noisy Systems. Cambridge University Press, Cambridge; 2010.

    Book  Google Scholar 

  29. Lo C, Chung TK: First passage time problem for the Ornstein–Uhlenbeck neuronal model. 13th International Conference on Neural Information Processing 2006, 324–333.

    Chapter  Google Scholar 

  30. Jones E, Oliphant T, Peterson P, et al: SciPy: open source scientific tools for Python; 2001.

    Google Scholar 

  31. Galasi M: GNU Scientific Library Reference Manual. 3rd edition. Network Theory Ltd, Godalming; 2009.

    Google Scholar 

  32. Ditlevsen S, Lánský P: Only through perturbation can relaxation times be estimated. Phys Rev E 2012., 86(5): Article ID 050102 Article ID 050102

    Google Scholar 

  33. Jahn P, Berg RW, Hounsgaard J, Ditlevsen S: Motoneuron membrane potentials follow a time inhomogeneous jump diffusion process. J Comput Neurosci 2011, 31: 563–579. 10.1007/s10827-011-0326-z

    Article  MathSciNet  Google Scholar 

Download references

Acknowledgements

SD is supported by the Danish Council for Independent Research|Natural Sciences. AL acknowledges support from NSERC Canada. AI is supported by a PGS-D scholarship from NSERC Canada. The work is part of the Dynamical Systems Interdisciplinary Network, University of Copenhagen.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Alexandre Iolov.

Additional information

Competing Interests

The authors declare that they have no competing interests.

Authors’ Contributions

All authors contributed equally to the writing of this paper. All authors read and approved the final manuscript.

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Authors’ original file for figure 3

Authors’ original file for figure 4

Authors’ original file for figure 5

Authors’ original file for figure 6

Authors’ original file for figure 7

Authors’ original file for figure 8

Authors’ original file for figure 9

Authors’ original file for figure 10

Authors’ original file for figure 11

Authors’ original file for figure 12

Authors’ original file for figure 13

Authors’ original file for figure 14

Authors’ original file for figure 15

Authors’ original file for figure 16

Authors’ original file for figure 17

Authors’ original file for figure 18

Authors’ original file for figure 19

Authors’ original file for figure 20

Authors’ original file for figure 21

Authors’ original file for figure 22

Authors’ original file for figure 23

Authors’ original file for figure 24

Authors’ original file for figure 25

Authors’ original file for figure 26

Authors’ original file for figure 27

Authors’ original file for figure 28

Authors’ original file for figure 29

Authors’ original file for figure 30

Authors’ original file for figure 31

Authors’ original file for figure 32

Authors’ original file for figure 33

Authors’ original file for figure 34

Authors’ original file for figure 35

Authors’ original file for figure 36

Authors’ original file for figure 37

Authors’ original file for figure 38

Authors’ original file for figure 39

Authors’ original file for figure 40

Authors’ original file for figure 41

Authors’ original file for figure 42

Authors’ original file for figure 43

Authors’ original file for figure 44

Authors’ original file for figure 45

Authors’ original file for figure 46

Authors’ original file for figure 47

Authors’ original file for figure 48

Authors’ original file for figure 49

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Iolov, A., Ditlevsen, S. & Longtin, A. Fokker–Planck and Fortet Equation-Based Parameter Estimation for a Leaky Integrate-and-Fire Model with Sinusoidal and Stochastic Forcing. J. Math. Neurosc. 4, 4 (2014). https://doi.org/10.1186/2190-8567-4-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/2190-8567-4-4

Keywords