- Email: [email protected]
On the evaluation of firing densities for periodically driven neuron models
Available online at www.sciencedirect.com
Mathematical Biosciences 214 (2008) 122–133 www.elsevier.com/locate/mbs
On the evaluation of firing densities for periodically driven neuron models q Aniello Buonocore a,*, Luigia Caputo b, Enrica Pirozzi a a
Dipartimento di Matematica e Applicazioni, Universita` di Napoli Federico II, Via Cintia, 80126 Napoli, Italy b Dipartimento di Matematica, Universita` di Torino Via Carlo Alberto 10, 10123 Torino, Italy Received 29 December 2007; received in revised form 1 February 2008; accepted 9 February 2008 Available online 23 February 2008
Abstract The leaky integrate-and-fire model for neuronal spiking events driven by a periodic stimulus is studied by using the Fokker–Planck formulation. To this purpose, an essential use is made of the asymptotic behavior of the first-passage-time probability density function of a time homogeneous diffusion process through an asymptotically periodic threshold. Numerical comparisons with some recently published results derived by a different approach are performed. Use of a new asymptotic approximation is then made in order to design a numerical algorithm of predictor–corrector type to solve the integral equation in the unknown first-passage-time probability density function. Such algorithm, characterized by a reduced (linear) computation time, is seen to provide a high computation accuracy. Finally, it is shown that such an approach yields excellent approximations to the firing probability density function for a wide range of parameters, including the case of high stimulus frequencies. Ó 2008 Elsevier Inc. All rights reserved. Keywords: Leaky integrate-and-fire model; Fokker–Planck equation; First-passage-time; Integral equations
1. Introduction Since the classical paper by Gerstein and Mandelbrot [1], numerous attempts have been made to formulate stochastic models for single neuron’s activity able to reproduce some of the essential features of the behavior exhibited by living neural cells under spontaneous or stimulated conditions (see [2–4] and reference therein). While in complex deterministic models, such as the Hodgkin– Huxley or the FitzHugh–Nagumo models, the generation of action potentials is automatically included, for stochastic models a firing threshold has to be introduced. The generation times of action potentials (spikes in the following) are then identified with the instants when the membrane depolarization, which is modeled on the analogy of a q Work performed under partial support of MIUR (PRIN 2005), GNCS-INdAM and Campania Region. * Corresponding author. E-mail address: [email protected] (A. Buonocore).
0025-5564/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2008.02.003
Brownian particle, attains the firing threshold. These instants are uniquely specified for deterministic models, whereas in the case of stochastic models, similarly for what occurs for many types of real neurons, the time interval that elapses between a spike and the next one is at least partially random. In [1], such a time interval was modeled as a random variable to be mathematically investigated by studying the related first-passage-time (FPT) problem. Within such a framework, the time course of the depolarization of the neuronal membrane potential was viewed as a diffusion process characterized by constant drift and infinitesimal variance, namely as a Wiener process with drift. In such a case, the probability density function (pdf) of the FPT through the firing threshold can be written in closed form, and the statistics of the firing times can be quantitatively studied. However, despite the excellent fittings of a variety of data, the neuronal model based on the Wiener process is undoubtedly too crude a simplification of the physiological reality as, for instance, it does not include the well-known spontaneous decay of the
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
membrane depolarization. A diffusion model that includes this specific feature has therefore been successively proposed and rigorously derived (see, for instance, [5–7]). This is a stochastic counterpart of the so called phenomenological leaky-integrate-and-fire (LIF) deterministic model. It identifies with the Ornstein–Uhlenbeck (OU) diffusion process having infinitesimal moments x A1 ðxÞ ¼ þ A h
A2 ¼ 2D
where h denotes the time constant of the neuronal membrane depolarization x and A 2 R and D > 0 are constants describing the input properties. The diffusion model of membrane depolarization, henceforth denoted as the LIF model, ultimately rests on the assumption that the inputs to the neuron are summed up over time until the firing threshold is crossed. At that instant a spike is generated and the accumulation starts anew. A brief review of the related essential literature will now follow in order to provide a better characterization of our results. As is well-known, the LIF model has been the object of numerous investigations aiming to elucidate some of its hidden features. In particular, differently from the case of Wiener neuronal model, in the LIF case the FPT pdf through a constant firing threshold of physiological relevance is not known (see [8]). Nevertheless, some information on the moments of the firing time can be obtained, though at the cost of very cumbersome numerical computations (see, for instance, [9,10]). In this respect, a conceptually interesting and practically relevant question was raised in [11], where the determination of the neuronal firing pdf under a sine-like varying input was approached within the LIF model. There, for suitable choices of the numerous parameters appearing in the model, two approximations to the FPT pdf were obtained. The first approximation was constructed by using the method of images, whereas the second one was based on a linearization of the quadratic potential well (the opposite of the integral of the drift) from the resting potential to the firing threshold. A thorough discussion of the motivations for considering periodic stimuli on the grounds of neuro-physiological evidence can be found in [12]. There, the origin of this kind of stimuli is also discussed with reference to intensity or frequency features of post-synaptic excitatory and inhibitory potentials. In addition, this seminal paper also includes the fundamental question as to when and why the FPT pdf may be taken as representative of the interspike pdf. In short, two alternatives have to be considered: (i) after each firing, the membrane potential is instantly reset to a fixed initial value x0 and an analogous reset occurs for the periodic stimulus as well; (ii) the resetting of the membrane potential to the initial value x0 is not accompanied by a similar reset of the stimulus that keeps evolving without being affected by the firing occurrence. In [12], case (i) is referred to as to an ‘endogenous’ stimulus, while (ii) depicts
123
a situation where the involved stimulus is of an ‘exogenous’ nature. Clearly, for endogenous stimuli, the spike train, as experimentally recorded, is to be viewed as a renewal process. Hence, the FPT pdf through the firing threshold of the diffusion process (originating at a fixed state) that is used to model the membrane potential dynamics, also represents the pdf of the inter-arrival time between two successive spikes. Quite different is the situation for exogenous stimuli. Indeed, in such a case the pdfs of the successive interarrival times differ from the first one, as well as among themselves, because of the changes of initial phases of the periodic stimulus. The determination of the pdfs of the inter-arrival times is then a new challenging problem. In [13] this problem is approached by averaging in ð0; 2pÞ, by means of a weight function h1 ðuÞ, the FPT pdf with respect to the initial phase u of a periodic sine-like stimulus. This weight function is obtained as the fixed point of an operator named ‘Stochastic Phase Transition Operator’ (SPTO) by the authors. Such an approach is then extended in [14] to more general P-periodic stimuli. The determination of the FPT pdf, synthetizing the most relevant information both for endogenous and for exogenous stimuli, is then carried in [13,14] by solving first or second order integral equations by numerical techniques. An interesting approach adopted in [15,16] to arise to the pdf of inter-arrival times should also be recalled. This is based on use of the ‘Phase Transition Matrix’ that, in some way, plays the role of the above SPTO. In addition to the above numerical techniques, here the time-dependent exponential approximation to the FTP pdf proposed in [17] is also used. The time-dependent rate, interpreted as the firing ‘risk’, depends both on previous firing time and on the stimulus. However, such rate is specified by resorting to physico-chemical models, such as Arrhenius, sigmoidal and Tuckwell models. Hence, the exploited exponential approximation must be viewed as a purely phenomenological tool, since no constructive derivation of it has been performed. The LIF model in presence of a periodic sine-like stimulus has been more recently re-considered in [18]. The specific aim was to prove that an exponential approximation to FPT pdf holds, with a time-dependent rate kðtÞ whose period equals that of the stimulus. This is achieved under a twofold assumption, namely low stimulus frequency and small noise intensity: use of Floquet eigenvalues and eigenfunctions of the backward Fokker–Planck operator associated to the LIF stochastic differential equation could therefore be implemented. It is not difficult to realize that such a method may also be used for P-periodic stimuli of more general types. In [19] it is finally shown that under exogenous stimulus condition the function h1 ðuÞ of reference R[13,14] can be replaced by the normalized rate P kðtÞ= 0 kðsÞ ds. In the sequel, we shall take up the problem considered in [18] with the aim of relaxing the above assumptions. We shall indeed prove, by a novel approach, that an exponential approximation to the FPT pdf indeed holds under the
124
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
sole assumption of small noise intensity, whatever the stimulus frequency is. We shall then compare our results with those obtained in [18], both analytically and via numerical computations. These are based on the well-known procedure of [20], also implemented in [14,18], as well as on a new method leading to a predictor-corrector type algorithm. 2. The LIF model
With the aim to describe the neuronal dynamics modeled by Eq. (1) and to obtain an approximation to the instantaneous firing rate, we start from the Fokker–Planck (FP) equation: o f ðx; tjx0 ; t0 Þ ¼ LðtÞf ðx;tjx0 ;t0 Þ ot i o hx o Acosðxt þ /Þ þ D f ðx;tjx0 ;t0 Þ ox h ox ð2Þ
We consider the non-stationary LIF neuronal model whose dynamics is described by the following Langevin equation: pffiffiffiffiffiffi x x_ ¼ þ A cos ðxt þ /Þ þ 2DnðtÞ ð1Þ h where xðtÞ stands for the electric potential difference across the neuron’s somatic membrane, namely the so-called membrane potential. Here, we have shifted the potential scale in such a way that the physiological resting potential is zero. Assuming xðt0 Þ ¼ x0 , the neuron fires when xðtÞ attains for the first time the firing threshold, that here is assumed to be a constant, say S. In Eq. (1) the parameter h represents the deterministic relaxation time towards the resting potential x ¼ 0 while the constant D summarizes the contributions to the membrane potential changes induced by the post-synaptic potentials. In addition, nðtÞ is assumed to be a standard white Gaussian noise whereas A cosðxt þ /Þ is viewed as a P-periodic deterministic stimulus ðP ¼ 2p=xÞ. In Fig. 1 we have plotted sample paths obtained via a discretization of Eq. (1) for a small driving frequency ðx ¼ 0:05=hÞ and for a fourfold larger value ðx ¼ 0:2=hÞ. Other parameters are specified in the caption.
where f ðx; tjx0 ; t0 Þ is the transition probability density function (pdf) of a diffusion process X ðtÞ with state space ðr1 ; r2 Þ ð1; þ1Þ, both end points being natural barriers and where LðtÞ denotes the Fokker–Planck operator. We now perform the transformation 0 0 x0 ¼ x0 sðt0 Þ x ¼ x sðtÞ ð3Þ 0 t00 ¼ t0 t ¼t ox f 0 x0 ; t0 jx00 ; t00 ¼ f ðx; tjx0 ; t0 Þ 0 ox where we have set sðtÞ ¼
Ah cos ðxt þ /Þ þ xh sin ðxt þ /Þ et=h : 2 2 1þx h
Hence Eq. (2) becomes
o 0 0 o x0 o 0 0 0 f x ; tjx0 ; t0 ¼ 0 þ D 0 f x ; tjx00 ; t0 ot ox h ox
ð4Þ
that is the FP equation for the transition pdf f 0 x0 ; tjx00 ; t0 of a stationary Ornstein–Uhlenbeck process. Eq. (4) is the equivalent FP formulation of the transformed LIF model:
Fig. 1. Sample paths obtained via discretization of Eq. (1) ending at the instant of first attainment of threshold S ¼ 1, for the two indicated values of the driving frequencies. The periodic stimuli (times h) are also shown. For an easier comparison, time scales, measured in units of h, have been chosen in a way to yield FPTs at roughly the same time ðt 300Þ. Other parameters have been chosen as follows: x0 ¼ 0:2; h ¼ 1, A ¼ 0:1; / ¼ 0 and D ¼ 0:1.
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
x_ 0 ¼
x0 pffiffiffiffiffiffi þ 2DnðtÞ: h
ð5Þ
Indeed, this can be viewed as describing the dynamics of the neuronal voltage in the presence of the asymptotically periodic firing threshold S 0 ðtÞ ¼ S
Ah cosðxt þ /Þ þ xh sinðxt þ /Þ et=h ; 2 2 1þx h ð6Þ
with the fixed starting voltage x0 ðt0 Þ ¼ x00 x0 sðt0 Þ: Fig. 2 shows two sample paths obtained via Eq. (5). Here, the indicated values of x refer to the angular frequencies of the asymptotically periodic transformed thresholds (6).
x_ ¼
125
oU ðx; tÞ pffiffiffiffiffiffi þ 2DnðtÞ; ox
ð8Þ
in which the potential 2
U ðx; tÞ ¼
½x Ah cosðxt þ /Þ 2h
explicitly appears. Note that minx
By making use of Floquet theory, the authors were able to show that under the two conditions DU ðtÞ=D > 4
We recall that, for a regular stochastic process fX ðtÞ; t 2 ½0; þ1Þg; the first-passage-time (FPT) through a time varying threshold SðtÞ is the random variable T ¼ inf ft P 0 : X ðtÞ > SðtÞg;
X ðt0 Þ ¼ x0 :
The pdf of T, gðtÞ gðSðtÞ; tjx0 ; t0 Þ ¼
d ProbðT 6 tÞ; dt
ð9Þ
xh1
3. FPT and neuronal firing density
ð7Þ
is the relevant quantity in our context. Recently, an equivalent formulation of the non-stationary OU process for the description of the neuronal firing properties has been given in [18]. There, the authors consider the LIF model (1) written in the form:
the following approximation to the instantaneous firing rate kðtÞ holds: qffiffiffiffiffiffiffiffi
DU ðtÞ 1 erf D DU ðtÞ h i: kðtÞ k1 ðtÞ ’ ð10Þ hD 1 exp DU ðtÞ D Here k1 ðtÞ is the first instantaneous eigenvalue of the backward operator related to LðtÞ as defined in Eq. (2). Furthermore, as described in [18,19], the waiting-time probability (i.e. the interspike interval pdf) can be expressed as Z t P ðtÞ ’ exp k1 ðsÞ ds : ð11Þ 0
Hence, from (10) and (11) one obtains:
Fig. 2. Sample paths obtained via discretization of Eq. (5) ending at the instant of first attainment of threshold S 0 ðtÞ given in Eq. (6) with S ¼ 1, for the two indicated values of the driving frequencies. Time scales, measured in units of h, have been chosen in a way to yield FPTs at roughly the same time ðt 300Þ. Other parameters have been chosen as follows: x0 ¼ 0:2; h ¼ 1; A ¼ 0:1; / ¼ 0 and D ¼ 0:1.
126
gS ðtÞ ¼
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
dP ðtÞ ¼ kðtÞ e dt
Rt 0
kðsÞ ds
ð12Þ
where gs ðtÞ denotes the firing pdf. As the authors state, their approach yields analytic and tractable expression not only in the case of the linear response (i.e., for weak stimuli, A 1) or of weak-noise limits, as far as conditions (9) are satisfied. The first of conditions (9) involves the value of the potential at the boundary: specifically, the difference between the potential value on the threshold and its minimum must always be larger than a few units of noise intensity D. This assumption implies an exponential timescale separation between the average (very long) time in which the threshold is reached from the minimum of the potential and the (very short) time of the deterministic relaxation towards the potential minimum. The second condition in (9) requires that the driving frequency x be much smaller than the relaxation rate appearing in the parabolic potential. It must be emphasized that approximation in Eq. (12) is valid for t h. In order to provide a quantitative evaluation of the goodness of the approximation to firing pdf discussed in the foregoing, as well as for some different approximations that will be given in the sequel, we outline hereafter a procedure able to estimate the FPT pdf through any C 2 ð½0; þ1ÞÞ threshold that holds for any time-homogeneous diffusion process for which the transition pdf is known. These conditions are satisfied by the OU process, which is the solution of Eq. (5), and by the threshold S 0 ðtÞ defined in (6). Let fOU ðx; tjy; sÞ be the transition pdf obtained as the solution of FP Eq. (4) satisfying the associated initial delta-condition. This is a Gaussian pdf with mean y expfðt sÞ=hg and variance Dh½1 expf2ðt sÞ=hg. Hence, setting t0 ¼ 0 and SðtÞ SðtÞ SðsÞeðtsÞ=h _ fOU ½SðtÞ;tjSðsÞ;s wðt;sÞ ¼ SðtÞ þ 2 h½1 e2ðtsÞ=h h t=h _ þ SðtÞ 2 SðtÞ x0 e wx0 ðtÞ ¼ SðtÞ fOU ½SðtÞ; tjx0 ; 0; 2t=h h½1 e h ð13Þ 0
the FPT pdf gðtÞ of the OU process through threshold S ðtÞ, modeling the neuronal firing pdf, satisfies the following Volterra integral equation of second kind ([20]): Z t gðtÞ ¼ wx0 ðtÞ þ wðt; sÞgðsÞ ds ð14Þ 0
whose kernel is non-singular since lims"t wðt; sÞ ¼ 0 as seen from the first of Eq. (13). As is well-known, the class of time dependent thresholds for which Eq. (14) leads to closed form solutions does not include threshold S 0 ðtÞ given in (6). Hence, one is obliged to resort to numerical evaluations that, however, can provide reliable approximations to the unknown firing density. A very simple and efficient numerical method, based on the composite trapezoidal rule, was implemented in [20], where an analysis of its convergence features was also given. This method, exhibiting a quadratic time-complexity, will be
now implemented to test the goodness of approximation (12), as well as of our forthcoming analytical and numerical approximations. Fig. 3 shows the FPT pdfs obtained by the method of [20] and by approximations (12) for different values of the driving frequencies and fixed other parameters. The former is recognized as it originates at zero. Note that approximations (12), as well as the phase locking to the input frequencies, improve as the driving frequencies decrease. 4. A different asymptotic approximation In [21], making use of the integral Eq. (14), the asymptotic behavior (for large times and large average threshold values) of the function gðtÞ in the case of asymptotically periodic thresholds, was determined for a class of onedimensional diffusion processes possessing a steady-state density. Along similar lines, we shall now determine the asymptotic behavior of the FPT pdf for the OU process through the time dependent threshold S 0 ðtÞ. To this purpose, let W OU ðxÞ denote the OU steady state pdf (that is Gaussian with zero mean and variance Dh) and let 2p 0 V ðtÞ ¼ lim S t þ n n!1 x ¼ S Ah
cos ðxt þ /Þ þ xh sin ðxt þ /Þ : 1 þ x2 h2
It is easy to prove that 2p lim w t þ n ; s ¼ hðtÞ n!1 x
and
ð15Þ
2p lim wx0 t þ n n!1 x
¼ hðtÞ where
V ðtÞ hðtÞ ¼ V_ ðtÞ W OU ½V ðtÞ: h
ð16Þ
As for the function hðtÞ, we note that V ðtÞ _ hðtÞ > 0; 8t > 0 () V ðtÞ > 0; 8t > 0 h 1 x2 h2 cosðxt þ /Þ þ 2xh sinðxt þ /Þ S () A > 0; 8t > 0 h 1 þ x2 h2 1 x2 h2 cosðxt þ /Þ þ 2xh sinðxt þ /Þ () S > Ahmax t>0 1 þ x2 h 2 S ð17Þ () A < : h
In addition, there holds: lim hðtÞ ¼ 0;
S!þ1
8t > 0:
ð18Þ
Under condition (17), after setting Z 2p=x x k¼ hðsÞ ds; 2p 0 in [21] it was proved that there exists a non-negative strictly increasing function /ðtÞ such that
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
127
Fig. 3. Plots of the numerically evaluated FPT pdfs and of the corresponding asymptotic approximations (12) for different values of the driving frequency x. Other parameters have been chosen as follows: S ¼ 1; x0 ¼ 0:2; h ¼ 1; A ¼ 0:1, / ¼ 0. Such values imply that in order for the second of conditions (12) to hold, it must be D < 0:10125. Hence, here we have taken the intermediate value D ¼ 0:05.
Z
/ðtÞ
hðsÞ ds ¼ kt;
/ðt þ np=xÞ ¼ /ðtÞ þ np=x:
0
The boundness of the periodic functions in the right-hand side of (16) implies:
t=h lim h / ¼ 0: ð19Þ S!þ1 k This, in turn, yields t=h ¼ þ1: lim / S!þ1 k
ð20Þ
Finally, making use of (20) we obtain: h i h i wx0 / t=h w / t=h ; / s=h k k k h i ¼ lim h i lim ¼ 1: t=h S!þ1 S!þ1 h / t=h h / k k ð21Þ
By virtue of (19) and (21), for large S the following asymptotic behavior of the FPT pdf follows: Rt hðsÞ ds ; t h: ð22Þ g0 ðtÞ ’ hðtÞ e 0 Due to the identification of g0 ðtÞ with an approximation to the neuronal firing pdf, the function hðtÞ takes the meaning of neuronal firing rate. The excellent agreement of the asymptotic FPT pdfs of the OU process through constant thresholds with exponenpffiffiffiffiffiffiffiffi ffi tial densities, whenever the thresholds exceed 2 2Dh, was pinpointed in [22,23]. Here we now look for a specific condition such as to insure that an asymptotic exponential trend of gðtÞ holds for our case of the transformed timedependent threshold S 0 ðtÞ. To this purpose, we note that Ah min V ðtÞ ¼ S pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t>0 1 þ x2 h2
128
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
so that inf min V ðtÞ ¼ S Ah:
modeling), the condition (23) can be re-expressed as follows:
x>0 t>0
Hence, 8x > 0, we conclude that approximation (22) holds for all firing thresholds S such that pffiffiffi S Ah pffiffiffiffiffiffi > 2 2: ð23Þ Dh Indeed, condition (23) implies A < S=h which, in turn, has been seen to insure hðtÞ > 0; 8t > 0. Therefore, asymptotic approximation (22) holds. On the other hand, if one assumes that S is assigned (which is customary in neuronal
D<
ðS AhÞ2 : 8h
ð24Þ
This condition, with S > Ah, is suggestive of the effects of the global level of input randomness on neuronal membrane. In Fig. 4 the FPT pdfs obtained as numerical solutions of integral Eq. (14) and the corresponding asymptotic approximations (22) are plotted for four different angular frequencies and for a fixed set of other parameters. The
Fig. 4. Plots of the numerically evaluated FPT pdfs and of the corresponding asymptotic approximations (22) for different values of the driving frequency x. Other parameters have been chosen as in Fig. 3.
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
excellent agreement is evident, apart from a slight overestimation of the latter around some of the peaks. The approximation (22) improves as t grows larger. Note the evident phase-locking with the threshold and hence, ultimately, with the driving frequency. We explicitly remark that this method can be extended to any periodic stimulus IðtÞ as far as sðtÞ of Eq. (3) can be obtained as solution of equation s_ ¼ s=h þ IðtÞ. For a further improvement of the accuracy of our approximation, hereafter we shall outline a numerical method for solving Eq. (14) that makes essential use of the asymptotic form (22) of gðtÞ. 5. A predictor–corrector type procedure In this Section we refer to the LIF model (5) and (6) and to the related FPT pdf gðtÞ both as the solution of integral Eq. (14) and as its asymptotic approximation (22).
129
Let a > 0 and b > a be real numbers (whose meaning and specification will soon follow) and set: ( g~ðtÞ 0; 0 < t < ah; Rt gðtÞ ¼ ð25Þ hðsÞ ds ; t P ah g^ðtÞ hðtÞ e ah ( ~ x ðtÞ w ðtÞ; 0 < t < bh w x0 0 x ðtÞ ¼ w ð26Þ 0 ^ wx0 ðtÞ hðtÞ; t P bh; ( ~ sÞ wðt; sÞ; 0 < t s < bh wðt; wðt; sÞ ¼ ð27Þ ^ sÞ hðtÞ; t s P bh: wðt; To obtain an approximation g1 ðtÞ of gðtÞ, we perform the following substitutions on the right-hand side of Eq. (14): gðtÞ ! gðtÞ;
x ! w ; w x0 0
sÞ ! wðt; sÞ wðt;
and consider the following partition:
Fig. 5. Plots of the numerically evaluated FPT pdfs and of the corresponding asymptotic approximations (28) for different values of the driving frequency x. Here a ¼ 0:5 and b ¼ 5. Other parameters have been chosen as in Fig. 3.
130
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
0 6 t < ah;
bh 6 t < ða þ bÞh;
ah 6 t < bh;
t P ða þ bÞh: For the sake of conciseness, we shall now refer to case t P ða þ bÞh as, mutatis mutandis, this will also indicate how to proceed for all other intervals. Making use of (25)–(27) one can see that the following holds: gðtÞ ’ wx0 ðtÞ þ þ
Z
Z
ah
wðt;sÞ gðsÞds þ
Z
0
tbh
ah
t
sÞ ^ x ðtÞ þ wðt; gðsÞds ¼ w 0
Z
t
~ sÞ^ wðt; gðsÞds ¼ hðtÞ hðtÞ
tbh t
þ
Z
tbh
Hence,
Z
tbh
^ sÞ^ wðt; gðsÞds
ah
tbh
þ
sÞ wðt; gðsÞds
wðt; sÞ^ gðsÞds ¼ hðtÞe
Z
R tbh ah
g1 ðtÞ ¼ 8 wx0 ðtÞ; > > > Rt > > < wx0 ðtÞ þ ah wðt; sÞ^ gðsÞ ds; Rt gðsÞ ds; hðtÞ þ ah wðt; sÞ^ > > > R tbh > > : hðtÞ e ah hðsÞ ds þ R t wðt; sÞ^ gðsÞ ds; tbh
0 6 t < ah ah 6 t < bh bh 6 t < ða þ bÞh t P ða þ bÞh: ð28Þ
We point out that the computational complexity of the scheme (28) is linear in time due to the predictor-corrector type of the implemented numerical procedure. Hence, computations are now one order of magnitude faster than those of [20]. As for a and b (h fixed) we take:
tbh
g^ðsÞds
ah hðsÞ ds
þ
Z
t
tbh
wðt; sÞ^ gðsÞds:
(i) a such that ah is very close to (but less then) the first maximum of wx0 ðtÞ; (ii) b such that the product bh represents the time span after which the process can be reasonably assumed to have entered the stationary regime.
Fig. 6. Numerically evaluated FPT pdfs and corresponding asymptotic approximations (28) for different values of noise intensity D. Here x ¼ p=30; a ¼ 0:5 and b ¼ 5. Other parameters have been chosen as in Fig. 3.
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
131
Fig. 7. As in Fig. 6 except that now x ¼ p=5.
6. Some comparisons We report some numerical comparisons between the improved evaluations of firing pdf g1 ðtÞ obtained via the predictor-corrector procedure (28) and the asymptotic estimations given by (22). Comparing Fig. 5 with Fig. 4, we note that the estimation g1 ðtÞ of gðtÞ obtained by means of the numerical procedure (28) is better then that given by the asymptotic approximation g0 ðtÞ defined in Eq. (22) for large t. In addition, the evaluation (28) appears to be in good agreement with the numerical solution of (14) also in the case of small values of t. The role of noise intensity D emerges from inspection of Figs. 6 and 7 referring to angular frequencies x ¼ p=30 (low frequency) and x ¼ p=5 (high frequency), respectively. Indeed, it is evident that the goodness of the approximations drops as D approaches its maximum allowed value specified by condition (24),
that equals 0.10125 in the case of the chosen parameters. Looking at Fig. 8, one sees that a plays a significant role when the period of the asymptotically periodic threshold largely exceeds the relaxation time h, to become less and less effective as the noise intensity D increases. In conclusion, while the first condition of (9) coincides with (24), the second of (9) is not required by our approaches. Hence, differently from the case of [18], our methods for neuronal firing pdfs evaluations appear to be valid for any frequency x. As far as the asymptotic approximations are concerned, the method in [18], when implementable, does not yield relevant discrepancies with respect to ours, the relative errors being in both cases within the range of 7%. By virtue of our second method these errors are further reduced under a still acceptable computation time that turns out to be linear. Finally, the goodness of such an approximation is noteworthy not only asymptotically, but even for small times.
132
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
Fig. 8. Numerically evaluated FPT pdfs and corresponding asymptotic approximations (28) for different values of noise intensity D and for x ¼ 4p. Here b ¼ 5. Other parameters have been chosen as in Fig. 3. The values of a have been selected in a way to obtain best fit of gðtÞ by g1 ðtÞ.
References [1] G.L. Gerstein, B. Mandelbrot, A random walk models for the Spyke activity of a single neuron, Biophys. J. 4 (1964) 41. [2] L.M. Ricciardi, A. Di Crescenzo, V. Giorno, A.G. Nobile, An outline of theoretical and algorithmic approaches to first passage time problems with applications to biological modeling, Sci. Mat. Jpn. 50 (1999) 247. [3] L.M. Ricciardi, P. La`nsky`, Diffusion models and neural activity, in: L. Nadel (Ed.), Encyclopedia of Cognitive Sciences I, Nature Publishing Group, London, 2002, p. 968. [4] M.A. Arbib, The Handbook of Brain Theory and Neural Networks, MIT, Cambridge, 2002. [5] R.B. Stein, A theoretical analysis of neuronal variability, Biophys. J. 5 (1965) 173. [6] B. Gluss, A model for neuron firing with exponential decay of potential resulting in diffusion equations for probability density, Bull. Math. Biol. 29 (2) (1967) 233. [7] L.M. Ricciardi, Diffusion processes and related topics in biology, Lecture Notes in Biomathematics, vol. 14, Springer, Berlin, 1977.
[8] R. Capocelli, L.M. Ricciardi, Diffusion approximation and first passage time problem for a model neuron, Biol. Cyber. 8 (2) (1971) 214. [9] L.M. Ricciardi, S. Sato, First-passage-time density and moments of Ornstein–Uhlenbeck process, J. Appl. Probab. 25 (1988) 43. [10] J. Inoue, S. Sato, L.M. Ricciardi, A note on the moments of the firstpassage-time of the Ornstein–Uhlenbeck process with a reflecting boundary, Ricerche Mat. XLVI (1) (1997) 87. [11] A.R. Bulsara, T.C. Elston, C.R. Doering, S.B. Lowen, K. Linderberg, Cooperative behavior in periodically driven noisy integrate-fire models of neuronal dynamics, Phys. Rev. E 53 (4) (1996) 3958. [12] P. La´nsky´, Sources of periodical force in noisy integrate-and-fire models of neuronal dynamics, Phys. Rev. E 55 (2) (1997) 2040. [13] T. Shimokawa, K. Pakdaman, S. Sato, Time-scale matching in the response of a leaky integrate-and-fire neuron model to periodic stimulus with additive noise, Phys. Rev. E 59 (3) (1999) 3427. [14] T. Shimokawa, K. Pakdaman, S. Sato, A first-passage-time analysis of the periodically forced noisy leaky integrate-and-fire model, Biol. Cyber. 83 (4) (2000) 327. [15] H.E. Plesser, T. Geisel, Markov analysis of stochastic resonance in a periodically driven integrate-and-fire neuron, Phys. Rev. E 59 (6) (1999) 7008.
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133 [16] H.E. Plesser, T. Geisel, Stochastic resonance in neuron models: endogenous stimulation revisited, Phys. Rev. E 63 (3) (2001) 031916-1. [17] H.E. Plesser, W. Gerstner, Noise in integrate-and-fire neurons: from stochastic input to escape rates, Neural Computat. 12 (2) (2000) 367. [18] M. Schindler, P. Talkner, P. Ha¨nggi, Escape rates in periodically driven Markov processes, Physica A 351 (2005) 40. [19] M. Kostur, M. Schindler, P. Talkner, P. Ha¨nggi, Neuron firing in driven non linear integrate-and-fire models, Math. Biosci 207 (2006) 302. [20] A. Buonocore, A.G. Nobile, L.M. Ricciardi, A new integral equation for the evaluation of first-passage-time probability densities, Adv. Appl. Probab. 19 (1987) 784.
133
[21] V. Giorno, A.G. Nobile, L.M. Ricciardi, On the asymptotic behavior of first-passage-time densities for one-dimensional diffusion processes and varying, Adv. Appl. Probab. 22 (4) (1990) 883. [22] A.G. Nobile, L.M. Ricciardi, L. Sacerdote, Exponential trends of Ornstein–Uhlenbeck first-passage-time densities, J. Appl. Probab. 22 (1985) 360. [23] A.G. Nobile, L.M. Ricciardi, L. Sacerdote, Exponential trends of first-passage-time densities for a class of diffusion processes with steady-state distribution, J. Appl. Probab. 22 (1985) 611.
Mathematical Biosciences 214 (2008) 122–133 www.elsevier.com/locate/mbs
On the evaluation of firing densities for periodically driven neuron models q Aniello Buonocore a,*, Luigia Caputo b, Enrica Pirozzi a a
Dipartimento di Matematica e Applicazioni, Universita` di Napoli Federico II, Via Cintia, 80126 Napoli, Italy b Dipartimento di Matematica, Universita` di Torino Via Carlo Alberto 10, 10123 Torino, Italy Received 29 December 2007; received in revised form 1 February 2008; accepted 9 February 2008 Available online 23 February 2008
Abstract The leaky integrate-and-fire model for neuronal spiking events driven by a periodic stimulus is studied by using the Fokker–Planck formulation. To this purpose, an essential use is made of the asymptotic behavior of the first-passage-time probability density function of a time homogeneous diffusion process through an asymptotically periodic threshold. Numerical comparisons with some recently published results derived by a different approach are performed. Use of a new asymptotic approximation is then made in order to design a numerical algorithm of predictor–corrector type to solve the integral equation in the unknown first-passage-time probability density function. Such algorithm, characterized by a reduced (linear) computation time, is seen to provide a high computation accuracy. Finally, it is shown that such an approach yields excellent approximations to the firing probability density function for a wide range of parameters, including the case of high stimulus frequencies. Ó 2008 Elsevier Inc. All rights reserved. Keywords: Leaky integrate-and-fire model; Fokker–Planck equation; First-passage-time; Integral equations
1. Introduction Since the classical paper by Gerstein and Mandelbrot [1], numerous attempts have been made to formulate stochastic models for single neuron’s activity able to reproduce some of the essential features of the behavior exhibited by living neural cells under spontaneous or stimulated conditions (see [2–4] and reference therein). While in complex deterministic models, such as the Hodgkin– Huxley or the FitzHugh–Nagumo models, the generation of action potentials is automatically included, for stochastic models a firing threshold has to be introduced. The generation times of action potentials (spikes in the following) are then identified with the instants when the membrane depolarization, which is modeled on the analogy of a q Work performed under partial support of MIUR (PRIN 2005), GNCS-INdAM and Campania Region. * Corresponding author. E-mail address: [email protected] (A. Buonocore).
0025-5564/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2008.02.003
Brownian particle, attains the firing threshold. These instants are uniquely specified for deterministic models, whereas in the case of stochastic models, similarly for what occurs for many types of real neurons, the time interval that elapses between a spike and the next one is at least partially random. In [1], such a time interval was modeled as a random variable to be mathematically investigated by studying the related first-passage-time (FPT) problem. Within such a framework, the time course of the depolarization of the neuronal membrane potential was viewed as a diffusion process characterized by constant drift and infinitesimal variance, namely as a Wiener process with drift. In such a case, the probability density function (pdf) of the FPT through the firing threshold can be written in closed form, and the statistics of the firing times can be quantitatively studied. However, despite the excellent fittings of a variety of data, the neuronal model based on the Wiener process is undoubtedly too crude a simplification of the physiological reality as, for instance, it does not include the well-known spontaneous decay of the
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
membrane depolarization. A diffusion model that includes this specific feature has therefore been successively proposed and rigorously derived (see, for instance, [5–7]). This is a stochastic counterpart of the so called phenomenological leaky-integrate-and-fire (LIF) deterministic model. It identifies with the Ornstein–Uhlenbeck (OU) diffusion process having infinitesimal moments x A1 ðxÞ ¼ þ A h
A2 ¼ 2D
where h denotes the time constant of the neuronal membrane depolarization x and A 2 R and D > 0 are constants describing the input properties. The diffusion model of membrane depolarization, henceforth denoted as the LIF model, ultimately rests on the assumption that the inputs to the neuron are summed up over time until the firing threshold is crossed. At that instant a spike is generated and the accumulation starts anew. A brief review of the related essential literature will now follow in order to provide a better characterization of our results. As is well-known, the LIF model has been the object of numerous investigations aiming to elucidate some of its hidden features. In particular, differently from the case of Wiener neuronal model, in the LIF case the FPT pdf through a constant firing threshold of physiological relevance is not known (see [8]). Nevertheless, some information on the moments of the firing time can be obtained, though at the cost of very cumbersome numerical computations (see, for instance, [9,10]). In this respect, a conceptually interesting and practically relevant question was raised in [11], where the determination of the neuronal firing pdf under a sine-like varying input was approached within the LIF model. There, for suitable choices of the numerous parameters appearing in the model, two approximations to the FPT pdf were obtained. The first approximation was constructed by using the method of images, whereas the second one was based on a linearization of the quadratic potential well (the opposite of the integral of the drift) from the resting potential to the firing threshold. A thorough discussion of the motivations for considering periodic stimuli on the grounds of neuro-physiological evidence can be found in [12]. There, the origin of this kind of stimuli is also discussed with reference to intensity or frequency features of post-synaptic excitatory and inhibitory potentials. In addition, this seminal paper also includes the fundamental question as to when and why the FPT pdf may be taken as representative of the interspike pdf. In short, two alternatives have to be considered: (i) after each firing, the membrane potential is instantly reset to a fixed initial value x0 and an analogous reset occurs for the periodic stimulus as well; (ii) the resetting of the membrane potential to the initial value x0 is not accompanied by a similar reset of the stimulus that keeps evolving without being affected by the firing occurrence. In [12], case (i) is referred to as to an ‘endogenous’ stimulus, while (ii) depicts
123
a situation where the involved stimulus is of an ‘exogenous’ nature. Clearly, for endogenous stimuli, the spike train, as experimentally recorded, is to be viewed as a renewal process. Hence, the FPT pdf through the firing threshold of the diffusion process (originating at a fixed state) that is used to model the membrane potential dynamics, also represents the pdf of the inter-arrival time between two successive spikes. Quite different is the situation for exogenous stimuli. Indeed, in such a case the pdfs of the successive interarrival times differ from the first one, as well as among themselves, because of the changes of initial phases of the periodic stimulus. The determination of the pdfs of the inter-arrival times is then a new challenging problem. In [13] this problem is approached by averaging in ð0; 2pÞ, by means of a weight function h1 ðuÞ, the FPT pdf with respect to the initial phase u of a periodic sine-like stimulus. This weight function is obtained as the fixed point of an operator named ‘Stochastic Phase Transition Operator’ (SPTO) by the authors. Such an approach is then extended in [14] to more general P-periodic stimuli. The determination of the FPT pdf, synthetizing the most relevant information both for endogenous and for exogenous stimuli, is then carried in [13,14] by solving first or second order integral equations by numerical techniques. An interesting approach adopted in [15,16] to arise to the pdf of inter-arrival times should also be recalled. This is based on use of the ‘Phase Transition Matrix’ that, in some way, plays the role of the above SPTO. In addition to the above numerical techniques, here the time-dependent exponential approximation to the FTP pdf proposed in [17] is also used. The time-dependent rate, interpreted as the firing ‘risk’, depends both on previous firing time and on the stimulus. However, such rate is specified by resorting to physico-chemical models, such as Arrhenius, sigmoidal and Tuckwell models. Hence, the exploited exponential approximation must be viewed as a purely phenomenological tool, since no constructive derivation of it has been performed. The LIF model in presence of a periodic sine-like stimulus has been more recently re-considered in [18]. The specific aim was to prove that an exponential approximation to FPT pdf holds, with a time-dependent rate kðtÞ whose period equals that of the stimulus. This is achieved under a twofold assumption, namely low stimulus frequency and small noise intensity: use of Floquet eigenvalues and eigenfunctions of the backward Fokker–Planck operator associated to the LIF stochastic differential equation could therefore be implemented. It is not difficult to realize that such a method may also be used for P-periodic stimuli of more general types. In [19] it is finally shown that under exogenous stimulus condition the function h1 ðuÞ of reference R[13,14] can be replaced by the normalized rate P kðtÞ= 0 kðsÞ ds. In the sequel, we shall take up the problem considered in [18] with the aim of relaxing the above assumptions. We shall indeed prove, by a novel approach, that an exponential approximation to the FPT pdf indeed holds under the
124
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
sole assumption of small noise intensity, whatever the stimulus frequency is. We shall then compare our results with those obtained in [18], both analytically and via numerical computations. These are based on the well-known procedure of [20], also implemented in [14,18], as well as on a new method leading to a predictor-corrector type algorithm. 2. The LIF model
With the aim to describe the neuronal dynamics modeled by Eq. (1) and to obtain an approximation to the instantaneous firing rate, we start from the Fokker–Planck (FP) equation: o f ðx; tjx0 ; t0 Þ ¼ LðtÞf ðx;tjx0 ;t0 Þ ot i o hx o Acosðxt þ /Þ þ D f ðx;tjx0 ;t0 Þ ox h ox ð2Þ
We consider the non-stationary LIF neuronal model whose dynamics is described by the following Langevin equation: pffiffiffiffiffiffi x x_ ¼ þ A cos ðxt þ /Þ þ 2DnðtÞ ð1Þ h where xðtÞ stands for the electric potential difference across the neuron’s somatic membrane, namely the so-called membrane potential. Here, we have shifted the potential scale in such a way that the physiological resting potential is zero. Assuming xðt0 Þ ¼ x0 , the neuron fires when xðtÞ attains for the first time the firing threshold, that here is assumed to be a constant, say S. In Eq. (1) the parameter h represents the deterministic relaxation time towards the resting potential x ¼ 0 while the constant D summarizes the contributions to the membrane potential changes induced by the post-synaptic potentials. In addition, nðtÞ is assumed to be a standard white Gaussian noise whereas A cosðxt þ /Þ is viewed as a P-periodic deterministic stimulus ðP ¼ 2p=xÞ. In Fig. 1 we have plotted sample paths obtained via a discretization of Eq. (1) for a small driving frequency ðx ¼ 0:05=hÞ and for a fourfold larger value ðx ¼ 0:2=hÞ. Other parameters are specified in the caption.
where f ðx; tjx0 ; t0 Þ is the transition probability density function (pdf) of a diffusion process X ðtÞ with state space ðr1 ; r2 Þ ð1; þ1Þ, both end points being natural barriers and where LðtÞ denotes the Fokker–Planck operator. We now perform the transformation 0 0 x0 ¼ x0 sðt0 Þ x ¼ x sðtÞ ð3Þ 0 t00 ¼ t0 t ¼t ox f 0 x0 ; t0 jx00 ; t00 ¼ f ðx; tjx0 ; t0 Þ 0 ox where we have set sðtÞ ¼
Ah cos ðxt þ /Þ þ xh sin ðxt þ /Þ et=h : 2 2 1þx h
Hence Eq. (2) becomes
o 0 0 o x0 o 0 0 0 f x ; tjx0 ; t0 ¼ 0 þ D 0 f x ; tjx00 ; t0 ot ox h ox
ð4Þ
that is the FP equation for the transition pdf f 0 x0 ; tjx00 ; t0 of a stationary Ornstein–Uhlenbeck process. Eq. (4) is the equivalent FP formulation of the transformed LIF model:
Fig. 1. Sample paths obtained via discretization of Eq. (1) ending at the instant of first attainment of threshold S ¼ 1, for the two indicated values of the driving frequencies. The periodic stimuli (times h) are also shown. For an easier comparison, time scales, measured in units of h, have been chosen in a way to yield FPTs at roughly the same time ðt 300Þ. Other parameters have been chosen as follows: x0 ¼ 0:2; h ¼ 1, A ¼ 0:1; / ¼ 0 and D ¼ 0:1.
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
x_ 0 ¼
x0 pffiffiffiffiffiffi þ 2DnðtÞ: h
ð5Þ
Indeed, this can be viewed as describing the dynamics of the neuronal voltage in the presence of the asymptotically periodic firing threshold S 0 ðtÞ ¼ S
Ah cosðxt þ /Þ þ xh sinðxt þ /Þ et=h ; 2 2 1þx h ð6Þ
with the fixed starting voltage x0 ðt0 Þ ¼ x00 x0 sðt0 Þ: Fig. 2 shows two sample paths obtained via Eq. (5). Here, the indicated values of x refer to the angular frequencies of the asymptotically periodic transformed thresholds (6).
x_ ¼
125
oU ðx; tÞ pffiffiffiffiffiffi þ 2DnðtÞ; ox
ð8Þ
in which the potential 2
U ðx; tÞ ¼
½x Ah cosðxt þ /Þ 2h
explicitly appears. Note that minx
By making use of Floquet theory, the authors were able to show that under the two conditions DU ðtÞ=D > 4
We recall that, for a regular stochastic process fX ðtÞ; t 2 ½0; þ1Þg; the first-passage-time (FPT) through a time varying threshold SðtÞ is the random variable T ¼ inf ft P 0 : X ðtÞ > SðtÞg;
X ðt0 Þ ¼ x0 :
The pdf of T, gðtÞ gðSðtÞ; tjx0 ; t0 Þ ¼
d ProbðT 6 tÞ; dt
ð9Þ
xh1
3. FPT and neuronal firing density
ð7Þ
is the relevant quantity in our context. Recently, an equivalent formulation of the non-stationary OU process for the description of the neuronal firing properties has been given in [18]. There, the authors consider the LIF model (1) written in the form:
the following approximation to the instantaneous firing rate kðtÞ holds: qffiffiffiffiffiffiffiffi
DU ðtÞ 1 erf D DU ðtÞ h i: kðtÞ k1 ðtÞ ’ ð10Þ hD 1 exp DU ðtÞ D Here k1 ðtÞ is the first instantaneous eigenvalue of the backward operator related to LðtÞ as defined in Eq. (2). Furthermore, as described in [18,19], the waiting-time probability (i.e. the interspike interval pdf) can be expressed as Z t P ðtÞ ’ exp k1 ðsÞ ds : ð11Þ 0
Hence, from (10) and (11) one obtains:
Fig. 2. Sample paths obtained via discretization of Eq. (5) ending at the instant of first attainment of threshold S 0 ðtÞ given in Eq. (6) with S ¼ 1, for the two indicated values of the driving frequencies. Time scales, measured in units of h, have been chosen in a way to yield FPTs at roughly the same time ðt 300Þ. Other parameters have been chosen as follows: x0 ¼ 0:2; h ¼ 1; A ¼ 0:1; / ¼ 0 and D ¼ 0:1.
126
gS ðtÞ ¼
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
dP ðtÞ ¼ kðtÞ e dt
Rt 0
kðsÞ ds
ð12Þ
where gs ðtÞ denotes the firing pdf. As the authors state, their approach yields analytic and tractable expression not only in the case of the linear response (i.e., for weak stimuli, A 1) or of weak-noise limits, as far as conditions (9) are satisfied. The first of conditions (9) involves the value of the potential at the boundary: specifically, the difference between the potential value on the threshold and its minimum must always be larger than a few units of noise intensity D. This assumption implies an exponential timescale separation between the average (very long) time in which the threshold is reached from the minimum of the potential and the (very short) time of the deterministic relaxation towards the potential minimum. The second condition in (9) requires that the driving frequency x be much smaller than the relaxation rate appearing in the parabolic potential. It must be emphasized that approximation in Eq. (12) is valid for t h. In order to provide a quantitative evaluation of the goodness of the approximation to firing pdf discussed in the foregoing, as well as for some different approximations that will be given in the sequel, we outline hereafter a procedure able to estimate the FPT pdf through any C 2 ð½0; þ1ÞÞ threshold that holds for any time-homogeneous diffusion process for which the transition pdf is known. These conditions are satisfied by the OU process, which is the solution of Eq. (5), and by the threshold S 0 ðtÞ defined in (6). Let fOU ðx; tjy; sÞ be the transition pdf obtained as the solution of FP Eq. (4) satisfying the associated initial delta-condition. This is a Gaussian pdf with mean y expfðt sÞ=hg and variance Dh½1 expf2ðt sÞ=hg. Hence, setting t0 ¼ 0 and SðtÞ SðtÞ SðsÞeðtsÞ=h _ fOU ½SðtÞ;tjSðsÞ;s wðt;sÞ ¼ SðtÞ þ 2 h½1 e2ðtsÞ=h h t=h _ þ SðtÞ 2 SðtÞ x0 e wx0 ðtÞ ¼ SðtÞ fOU ½SðtÞ; tjx0 ; 0; 2t=h h½1 e h ð13Þ 0
the FPT pdf gðtÞ of the OU process through threshold S ðtÞ, modeling the neuronal firing pdf, satisfies the following Volterra integral equation of second kind ([20]): Z t gðtÞ ¼ wx0 ðtÞ þ wðt; sÞgðsÞ ds ð14Þ 0
whose kernel is non-singular since lims"t wðt; sÞ ¼ 0 as seen from the first of Eq. (13). As is well-known, the class of time dependent thresholds for which Eq. (14) leads to closed form solutions does not include threshold S 0 ðtÞ given in (6). Hence, one is obliged to resort to numerical evaluations that, however, can provide reliable approximations to the unknown firing density. A very simple and efficient numerical method, based on the composite trapezoidal rule, was implemented in [20], where an analysis of its convergence features was also given. This method, exhibiting a quadratic time-complexity, will be
now implemented to test the goodness of approximation (12), as well as of our forthcoming analytical and numerical approximations. Fig. 3 shows the FPT pdfs obtained by the method of [20] and by approximations (12) for different values of the driving frequencies and fixed other parameters. The former is recognized as it originates at zero. Note that approximations (12), as well as the phase locking to the input frequencies, improve as the driving frequencies decrease. 4. A different asymptotic approximation In [21], making use of the integral Eq. (14), the asymptotic behavior (for large times and large average threshold values) of the function gðtÞ in the case of asymptotically periodic thresholds, was determined for a class of onedimensional diffusion processes possessing a steady-state density. Along similar lines, we shall now determine the asymptotic behavior of the FPT pdf for the OU process through the time dependent threshold S 0 ðtÞ. To this purpose, let W OU ðxÞ denote the OU steady state pdf (that is Gaussian with zero mean and variance Dh) and let 2p 0 V ðtÞ ¼ lim S t þ n n!1 x ¼ S Ah
cos ðxt þ /Þ þ xh sin ðxt þ /Þ : 1 þ x2 h2
It is easy to prove that 2p lim w t þ n ; s ¼ hðtÞ n!1 x
and
ð15Þ
2p lim wx0 t þ n n!1 x
¼ hðtÞ where
V ðtÞ hðtÞ ¼ V_ ðtÞ W OU ½V ðtÞ: h
ð16Þ
As for the function hðtÞ, we note that V ðtÞ _ hðtÞ > 0; 8t > 0 () V ðtÞ > 0; 8t > 0 h 1 x2 h2 cosðxt þ /Þ þ 2xh sinðxt þ /Þ S () A > 0; 8t > 0 h 1 þ x2 h2 1 x2 h2 cosðxt þ /Þ þ 2xh sinðxt þ /Þ () S > Ahmax t>0 1 þ x2 h 2 S ð17Þ () A < : h
In addition, there holds: lim hðtÞ ¼ 0;
S!þ1
8t > 0:
ð18Þ
Under condition (17), after setting Z 2p=x x k¼ hðsÞ ds; 2p 0 in [21] it was proved that there exists a non-negative strictly increasing function /ðtÞ such that
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
127
Fig. 3. Plots of the numerically evaluated FPT pdfs and of the corresponding asymptotic approximations (12) for different values of the driving frequency x. Other parameters have been chosen as follows: S ¼ 1; x0 ¼ 0:2; h ¼ 1; A ¼ 0:1, / ¼ 0. Such values imply that in order for the second of conditions (12) to hold, it must be D < 0:10125. Hence, here we have taken the intermediate value D ¼ 0:05.
Z
/ðtÞ
hðsÞ ds ¼ kt;
/ðt þ np=xÞ ¼ /ðtÞ þ np=x:
0
The boundness of the periodic functions in the right-hand side of (16) implies:
t=h lim h / ¼ 0: ð19Þ S!þ1 k This, in turn, yields t=h ¼ þ1: lim / S!þ1 k
ð20Þ
Finally, making use of (20) we obtain: h i h i wx0 / t=h w / t=h ; / s=h k k k h i ¼ lim h i lim ¼ 1: t=h S!þ1 S!þ1 h / t=h h / k k ð21Þ
By virtue of (19) and (21), for large S the following asymptotic behavior of the FPT pdf follows: Rt hðsÞ ds ; t h: ð22Þ g0 ðtÞ ’ hðtÞ e 0 Due to the identification of g0 ðtÞ with an approximation to the neuronal firing pdf, the function hðtÞ takes the meaning of neuronal firing rate. The excellent agreement of the asymptotic FPT pdfs of the OU process through constant thresholds with exponenpffiffiffiffiffiffiffiffi ffi tial densities, whenever the thresholds exceed 2 2Dh, was pinpointed in [22,23]. Here we now look for a specific condition such as to insure that an asymptotic exponential trend of gðtÞ holds for our case of the transformed timedependent threshold S 0 ðtÞ. To this purpose, we note that Ah min V ðtÞ ¼ S pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t>0 1 þ x2 h2
128
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
so that inf min V ðtÞ ¼ S Ah:
modeling), the condition (23) can be re-expressed as follows:
x>0 t>0
Hence, 8x > 0, we conclude that approximation (22) holds for all firing thresholds S such that pffiffiffi S Ah pffiffiffiffiffiffi > 2 2: ð23Þ Dh Indeed, condition (23) implies A < S=h which, in turn, has been seen to insure hðtÞ > 0; 8t > 0. Therefore, asymptotic approximation (22) holds. On the other hand, if one assumes that S is assigned (which is customary in neuronal
D<
ðS AhÞ2 : 8h
ð24Þ
This condition, with S > Ah, is suggestive of the effects of the global level of input randomness on neuronal membrane. In Fig. 4 the FPT pdfs obtained as numerical solutions of integral Eq. (14) and the corresponding asymptotic approximations (22) are plotted for four different angular frequencies and for a fixed set of other parameters. The
Fig. 4. Plots of the numerically evaluated FPT pdfs and of the corresponding asymptotic approximations (22) for different values of the driving frequency x. Other parameters have been chosen as in Fig. 3.
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
excellent agreement is evident, apart from a slight overestimation of the latter around some of the peaks. The approximation (22) improves as t grows larger. Note the evident phase-locking with the threshold and hence, ultimately, with the driving frequency. We explicitly remark that this method can be extended to any periodic stimulus IðtÞ as far as sðtÞ of Eq. (3) can be obtained as solution of equation s_ ¼ s=h þ IðtÞ. For a further improvement of the accuracy of our approximation, hereafter we shall outline a numerical method for solving Eq. (14) that makes essential use of the asymptotic form (22) of gðtÞ. 5. A predictor–corrector type procedure In this Section we refer to the LIF model (5) and (6) and to the related FPT pdf gðtÞ both as the solution of integral Eq. (14) and as its asymptotic approximation (22).
129
Let a > 0 and b > a be real numbers (whose meaning and specification will soon follow) and set: ( g~ðtÞ 0; 0 < t < ah; Rt gðtÞ ¼ ð25Þ hðsÞ ds ; t P ah g^ðtÞ hðtÞ e ah ( ~ x ðtÞ w ðtÞ; 0 < t < bh w x0 0 x ðtÞ ¼ w ð26Þ 0 ^ wx0 ðtÞ hðtÞ; t P bh; ( ~ sÞ wðt; sÞ; 0 < t s < bh wðt; wðt; sÞ ¼ ð27Þ ^ sÞ hðtÞ; t s P bh: wðt; To obtain an approximation g1 ðtÞ of gðtÞ, we perform the following substitutions on the right-hand side of Eq. (14): gðtÞ ! gðtÞ;
x ! w ; w x0 0
sÞ ! wðt; sÞ wðt;
and consider the following partition:
Fig. 5. Plots of the numerically evaluated FPT pdfs and of the corresponding asymptotic approximations (28) for different values of the driving frequency x. Here a ¼ 0:5 and b ¼ 5. Other parameters have been chosen as in Fig. 3.
130
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
0 6 t < ah;
bh 6 t < ða þ bÞh;
ah 6 t < bh;
t P ða þ bÞh: For the sake of conciseness, we shall now refer to case t P ða þ bÞh as, mutatis mutandis, this will also indicate how to proceed for all other intervals. Making use of (25)–(27) one can see that the following holds: gðtÞ ’ wx0 ðtÞ þ þ
Z
Z
ah
wðt;sÞ gðsÞds þ
Z
0
tbh
ah
t
sÞ ^ x ðtÞ þ wðt; gðsÞds ¼ w 0
Z
t
~ sÞ^ wðt; gðsÞds ¼ hðtÞ hðtÞ
tbh t
þ
Z
tbh
Hence,
Z
tbh
^ sÞ^ wðt; gðsÞds
ah
tbh
þ
sÞ wðt; gðsÞds
wðt; sÞ^ gðsÞds ¼ hðtÞe
Z
R tbh ah
g1 ðtÞ ¼ 8 wx0 ðtÞ; > > > Rt > > < wx0 ðtÞ þ ah wðt; sÞ^ gðsÞ ds; Rt gðsÞ ds; hðtÞ þ ah wðt; sÞ^ > > > R tbh > > : hðtÞ e ah hðsÞ ds þ R t wðt; sÞ^ gðsÞ ds; tbh
0 6 t < ah ah 6 t < bh bh 6 t < ða þ bÞh t P ða þ bÞh: ð28Þ
We point out that the computational complexity of the scheme (28) is linear in time due to the predictor-corrector type of the implemented numerical procedure. Hence, computations are now one order of magnitude faster than those of [20]. As for a and b (h fixed) we take:
tbh
g^ðsÞds
ah hðsÞ ds
þ
Z
t
tbh
wðt; sÞ^ gðsÞds:
(i) a such that ah is very close to (but less then) the first maximum of wx0 ðtÞ; (ii) b such that the product bh represents the time span after which the process can be reasonably assumed to have entered the stationary regime.
Fig. 6. Numerically evaluated FPT pdfs and corresponding asymptotic approximations (28) for different values of noise intensity D. Here x ¼ p=30; a ¼ 0:5 and b ¼ 5. Other parameters have been chosen as in Fig. 3.
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
131
Fig. 7. As in Fig. 6 except that now x ¼ p=5.
6. Some comparisons We report some numerical comparisons between the improved evaluations of firing pdf g1 ðtÞ obtained via the predictor-corrector procedure (28) and the asymptotic estimations given by (22). Comparing Fig. 5 with Fig. 4, we note that the estimation g1 ðtÞ of gðtÞ obtained by means of the numerical procedure (28) is better then that given by the asymptotic approximation g0 ðtÞ defined in Eq. (22) for large t. In addition, the evaluation (28) appears to be in good agreement with the numerical solution of (14) also in the case of small values of t. The role of noise intensity D emerges from inspection of Figs. 6 and 7 referring to angular frequencies x ¼ p=30 (low frequency) and x ¼ p=5 (high frequency), respectively. Indeed, it is evident that the goodness of the approximations drops as D approaches its maximum allowed value specified by condition (24),
that equals 0.10125 in the case of the chosen parameters. Looking at Fig. 8, one sees that a plays a significant role when the period of the asymptotically periodic threshold largely exceeds the relaxation time h, to become less and less effective as the noise intensity D increases. In conclusion, while the first condition of (9) coincides with (24), the second of (9) is not required by our approaches. Hence, differently from the case of [18], our methods for neuronal firing pdfs evaluations appear to be valid for any frequency x. As far as the asymptotic approximations are concerned, the method in [18], when implementable, does not yield relevant discrepancies with respect to ours, the relative errors being in both cases within the range of 7%. By virtue of our second method these errors are further reduced under a still acceptable computation time that turns out to be linear. Finally, the goodness of such an approximation is noteworthy not only asymptotically, but even for small times.
132
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133
Fig. 8. Numerically evaluated FPT pdfs and corresponding asymptotic approximations (28) for different values of noise intensity D and for x ¼ 4p. Here b ¼ 5. Other parameters have been chosen as in Fig. 3. The values of a have been selected in a way to obtain best fit of gðtÞ by g1 ðtÞ.
References [1] G.L. Gerstein, B. Mandelbrot, A random walk models for the Spyke activity of a single neuron, Biophys. J. 4 (1964) 41. [2] L.M. Ricciardi, A. Di Crescenzo, V. Giorno, A.G. Nobile, An outline of theoretical and algorithmic approaches to first passage time problems with applications to biological modeling, Sci. Mat. Jpn. 50 (1999) 247. [3] L.M. Ricciardi, P. La`nsky`, Diffusion models and neural activity, in: L. Nadel (Ed.), Encyclopedia of Cognitive Sciences I, Nature Publishing Group, London, 2002, p. 968. [4] M.A. Arbib, The Handbook of Brain Theory and Neural Networks, MIT, Cambridge, 2002. [5] R.B. Stein, A theoretical analysis of neuronal variability, Biophys. J. 5 (1965) 173. [6] B. Gluss, A model for neuron firing with exponential decay of potential resulting in diffusion equations for probability density, Bull. Math. Biol. 29 (2) (1967) 233. [7] L.M. Ricciardi, Diffusion processes and related topics in biology, Lecture Notes in Biomathematics, vol. 14, Springer, Berlin, 1977.
[8] R. Capocelli, L.M. Ricciardi, Diffusion approximation and first passage time problem for a model neuron, Biol. Cyber. 8 (2) (1971) 214. [9] L.M. Ricciardi, S. Sato, First-passage-time density and moments of Ornstein–Uhlenbeck process, J. Appl. Probab. 25 (1988) 43. [10] J. Inoue, S. Sato, L.M. Ricciardi, A note on the moments of the firstpassage-time of the Ornstein–Uhlenbeck process with a reflecting boundary, Ricerche Mat. XLVI (1) (1997) 87. [11] A.R. Bulsara, T.C. Elston, C.R. Doering, S.B. Lowen, K. Linderberg, Cooperative behavior in periodically driven noisy integrate-fire models of neuronal dynamics, Phys. Rev. E 53 (4) (1996) 3958. [12] P. La´nsky´, Sources of periodical force in noisy integrate-and-fire models of neuronal dynamics, Phys. Rev. E 55 (2) (1997) 2040. [13] T. Shimokawa, K. Pakdaman, S. Sato, Time-scale matching in the response of a leaky integrate-and-fire neuron model to periodic stimulus with additive noise, Phys. Rev. E 59 (3) (1999) 3427. [14] T. Shimokawa, K. Pakdaman, S. Sato, A first-passage-time analysis of the periodically forced noisy leaky integrate-and-fire model, Biol. Cyber. 83 (4) (2000) 327. [15] H.E. Plesser, T. Geisel, Markov analysis of stochastic resonance in a periodically driven integrate-and-fire neuron, Phys. Rev. E 59 (6) (1999) 7008.
A. Buonocore et al. / Mathematical Biosciences 214 (2008) 122–133 [16] H.E. Plesser, T. Geisel, Stochastic resonance in neuron models: endogenous stimulation revisited, Phys. Rev. E 63 (3) (2001) 031916-1. [17] H.E. Plesser, W. Gerstner, Noise in integrate-and-fire neurons: from stochastic input to escape rates, Neural Computat. 12 (2) (2000) 367. [18] M. Schindler, P. Talkner, P. Ha¨nggi, Escape rates in periodically driven Markov processes, Physica A 351 (2005) 40. [19] M. Kostur, M. Schindler, P. Talkner, P. Ha¨nggi, Neuron firing in driven non linear integrate-and-fire models, Math. Biosci 207 (2006) 302. [20] A. Buonocore, A.G. Nobile, L.M. Ricciardi, A new integral equation for the evaluation of first-passage-time probability densities, Adv. Appl. Probab. 19 (1987) 784.
133
[21] V. Giorno, A.G. Nobile, L.M. Ricciardi, On the asymptotic behavior of first-passage-time densities for one-dimensional diffusion processes and varying, Adv. Appl. Probab. 22 (4) (1990) 883. [22] A.G. Nobile, L.M. Ricciardi, L. Sacerdote, Exponential trends of Ornstein–Uhlenbeck first-passage-time densities, J. Appl. Probab. 22 (1985) 360. [23] A.G. Nobile, L.M. Ricciardi, L. Sacerdote, Exponential trends of first-passage-time densities for a class of diffusion processes with steady-state distribution, J. Appl. Probab. 22 (1985) 611.