Switching Time Statistics for Driven Neuron Models: Analytic Expressions versus Numerics

0.3 MB·English

by  Michael Schindler

Firing Time Statistics for Driven Neuron Models: Analytic Expressions versus
Numerics
Michael Schindler, Peter Talkner, and Peter Ha¨nggi
Institut fu¨r Physik, Universit¨at Augsburg, Universit¨atsstraße 1, D–86135 Augsburg, Germany
(Dated: February 4, 2008)
Analyticalexpressionsareputforwardtoinvestigatetheforcedspikingactivityofabstractneuron
models suchas thedrivenleaky integrate-and-firemodel. Themethod isvalid in awide parameter
regimebeyondtherestraininglimitsofweakdriving(linearresponse)and/orweaknoise. Thenovel
approximationisbasedon adiscretestateMarkovianmodelingofthefulllong-timedynamicswith
time-dependentrates. TheschemeyieldsexcellentagreementwithnumericalLangevinandFokker-
6
Plancksimulationsofthefullnon-stationarydynamics,notonlyforthefirst-passagetimestatistics,
0
butalso for theimportant interspikeinterval (residence time) distribution.
0
2
PACSnumbers: 87.19.La,05.40.-a,87.18.Sn,89.75.Hc

n
a
J The detailed modeling of neural behavior presents a reading

9 purnodmerisnteanntdcinhgalloefngneeuornalthceodinintrgigpuriningcivpilsetsa. toTwhaerdlseatkhye ∂tρ=L(t)ρ=∂x U′(x,t)ρ +D∂x2ρ, (2)
] (cid:0) (cid:1)
C integrate-and-fire (LIF) model, whose deterministic for- with the absorbing boundary and initial conditions
N mulation has been introduced long ago [1], likely is one
. of the most studied abstract neuron models [2]. It is ρ(a,t|x0,s)=0 for all t, s, and x0 (3)
o characterized by marked simplicity and lack of memory: ρ(x,s x0,s)=δ(x x0). (4)
i | −
b Whenever the neuron has been excited to fire a pulse it After firing the process immediately restarts at the in-
- is reset to a predefined state. The beneficial rˆole of an
q stantaneous minimum of the potential.
appreciabledoseofnoisehasprovedtobestowakeypart
[ The set of eqs. (1–4) defines our starting point for ob-
to the interspike statistics of neurons [2]. There exist by
taining the firing statistics of this driven neuron model.
2
now numerous studies and important generalizations of
v Thisisaratherintricateproblembecausethepresenceof
realistic synaptic models, mostly of a numerical nature,
5 non-stationarityandmultiple time-scalesfordrivingand
1 that demonstrate the rich behavior of the renewal firing relaxation, in combination with the absorbing boundary
0 probability, e.g. see in Ref. [3, 4, 5].
condition prohibits an analytical exact solution [7]. Our
1 In presence of a time-dependent input-stimulation the main objective is, nevertheless, to develop a most accu-
0
stochastic firing process becomes non-stationary, which rate analytical approximation that supersedes all prior
4
0 in turn significantly complicates the stochastic firing attempts known to us. Those attempts, in fact, all in-
/ statistics. Nevertheless, the signal transmission and its volve the use of either of the following limiting approxi-
o
i detection can exhibit a remarkable improvement via the mationschemessuchasthelimitoflinearresponsetheory
b phenomenon of Stochastic Resonance [6]. The dynamics (i.e.aweakstimulusf(t))[8],thelimitofasymptotically
q- of the neuronal firing probability emerges due to a large weaknoise[6,9]ortheuseofthemethodofimageswhich
: bombardment of synaptic spike events; consequently, it appears to present an uncontrollable approximation for
v
i iscustomarytoemployadiffusionapproximationforthe the case with λ = 0 [4, 5]. A most appealing numeri-
X stochasticdynamicsofthemembranepotentialx(t). The cal approach is b6ased on an exact integral equation for
r complexity of the driven abstract LIF model thus as- the first-passage time density of time-dependent Gauss-
a
sumes the archetype, non-stationary Langevin dynamics Markovprocesseswith an absorbing boundary [10]. Our
schemedetailedbelowyieldsnovelanalyticandtractable
expressions beyond the linear response and weak noise
x˙(t)= λx(t)+µ+f(t)+√2Dξ(t) (1) limit; these are limited solely by the use of a discrete,
−
Markovianstochasticdynamics forthe populationofthe
where the process starts at a time s at x(s) = x0 attracting domain and slowly varying (in comparison to
and fires when it reaches the threshold voltage x = a. intra-wellrelaxationtime-scale) stimuli f(t). As demon-
Here, f(t) presents a general, time-dependent stimu- strated below, this novelscheme indeed provides analyt-
lus which, for example, can be chosen to be oscilla- ical formulae that compare very favorably with precise
tory, and ξ(t) is white Gaussian noise. The dynamics numerical results of the full dynamics in eqs. (1, 2–4).
ofthe processx(t) is equivalently describedbya Fokker- Different from other approaches, we obtain the distribu-
Planck(FP)equationfortheconditionalprobabilityden- tion not only of the first-passage time but also of the
sityρ(x,t x0,s)inatime-dependentquadraticpotential, residence time, which is the more interesting variable,
| 2
U(x,t)=λ x xmin(t) /2withxmin(t)= µ+f(t) /λ, concerning neurons.
−
(cid:2) (cid:3) (cid:0) (cid:1)
2
To start, we approximate the solution to eqs. (2–4) error function leads to the well-known weak noise result
in the regime where the statistics of times at which the for the time-dependent rate [6, 7], i.e.
threshold is reached can be characterized by a time-
dependent firing rate κ(t) [7, 11]. κwn(t)=λ ∆U(t)/(πD) exp( ∆U(t)/D). (10)
−
This rate then follows from a time-scale separation in p
the full FP dynamics (2) with boundary condition (3). Firingtimedistributions.—Withtheexpressionforthe
Afterafew timesλ−1 ofthefastrelaxationtheprobabil- exit rate κ(t) in (9) we can calculate the properties of
ity density ρ(x,t) assumes a slowly varying pattern that interest, namely the densities for the first-passage time
decays with the rate κ(t). As in the time-independent andtheresidencetime[12]oftheattracting”integrating”
case, this slowly varying part of ρ(x,t) can be expressed state that covers the domain <x(t)<a.
−∞
by a product of the (normalized) instantaneous station- Thefirst-passagetimedistributionisgivenbytheneg-
ary solution ρ0(x,t) exp U(x,t)/D to the FP ative rate of change of probability finding the process at
equation, satisfying L(∝t)ρ0(x,(cid:8)t−) = 0, and(cid:9)a form func- time t in the ”integrating” state, i.e.
tion ζ(x,t). Our ansatz thus reads a
g(t s)= ∂t ρ(x,t x0,s)dx (11)
(t−s)>λ−1 t ′ ′ | − Z−∞ |
ρ(x,t|x0,s) ≃ ζ(x,t)ρ0(x,t)exp(cid:16)−Rsκ(s)ds(cid:17). (5) =κ(t) exp tκ(s′)ds′ , (12)
Theinitialtimesentersonlythroughtheexponentialfac- (cid:16)−Rs (cid:17)
tor. Thedependenceoftheconditionalprobabilityonthe Here, the integral over the spatial part of ρ(x,t) yields
initialvaluex0 decaysexponentiallyonthetimescaleλ−1 unity and the exponential factor in (5) is obtained. It
and therefore can be neglected for long times. gives the probability for the process to stay in the ”inte-
Deep inside the attracting well, i.e. for x a there grating” state from time s until t without interruption.
≪
is no sensible difference between ρ(x,t) and ρ0(x,t), and The distribution of the residence times h(τ), also
consequently ζ(x,t) approaches one. At the absorbing termedtheinterspikeintervaldensity,followsastheaver-
boundary, however,ρ(a,t) and consequently ζ(a,t) both ageofthefirst-passagetimedensityoverthedensityofre-
must vanish. The quantitative form of ζ(x,t) in the setting times s,whichcoincideswiththe firing rateκ(s).
crossoverregion follows from It thus reads:
L+(t)ζ(x,t)=0, (6) T g(τ +s s)κ(s)ds
h(τ)= lim −T | . (13)
where the potential in the backward operator L+(t) = T→∞R −TT κ(s)ds
U′(x,t)∂x+D∂x2 canbe linearizedabout the threshold R
−x=a. Eq.(6)thenyieldsfortheformfunctiontheresult Equations (12) and (13), together with the expression
for the rate (9) constitute the main results of this work.
Their quantitative validity for an extended parameter
′
U (a,t)
regime will be checked next.
ζ(x,t)=1 exp (x a) . (7)
− n − D o Numerical comparison.—We have employed three
methods for the numerical analysis. A first one is based
The rate κ(t)is determined by multiplying the FP equa-
on the Langevin equation (1) where the position x(t) is
tion (2) in the long-time limit (5) by the form func-
updated sequentially. For the second we have solved the
tionζ(x,t)andintegratingoverxfrom tothethresh-
−∞ FP equation(2) numerically,using a Chebychev colloca-
oldvoltagea. Indoingso,weaccountforprominentfinite
tion method to reduce the problem to a coupled system
barrier corrections, yielding
ofordinarydifferentialequations,seealso[13]. Thethird
−a∞dxζ(x,t)L(t)ζ(x,t)ρ0(x,t) method solves an integral equation for the first-passage
κ(t)=−R −a∞dxζ2(x,t)ρ0(x,t) . (8) tmimetehpordosbhaabvileitpyrdoevnidsietdyparnadcitsicdaelslycriidbeendtiinca[l1r0e].suAltlsl.three
R
In the prior stimulating work [4] it has been left open
Upon insertion of eq. (7) for the form function one can
inwhatrelevantparameterregimetheemployedapprox-
exactlyperformthe integrationsandobtainsfor the rate
imation possesses validity [5]. Here, using a periodic
modulationf(t)=Acos(ωt) we like to determine a min-
∆U(t) 1 erf ∆U(t)/D imal set of relevant parameters that can be taken for
κ(t)=λ − , (9)
D 1 exp(cid:0)(p∆U(t)/D)(cid:1) comparison with physiological measurements. The most
− − generalLangevin equation (1) considered has seven con-
where ∆U(t) denotes the instantaneous potential height stantparameters,includingthethresholdaandthereset-
atthethresholdasseenfromtheminimum, anderf(z)is position x0. Three of them, (λ,µ,a), can be chosen to
the errorfunction. For verysmallD anexpansionofthe be(1,0,1)bytransformingtodimensionlesscoordinates,
3
using the time-unit λ−1, the space-unit (a µ/λ), and where U¯+ and U¯− are the maximum and the minimum
the coordinate-origin at µ/λ. The resulting−parameters ofU¯(1,t¯)duringawholeperiodofmodulation. Theseso
thus read (the bars indicate dimensionless coordinates) chosenparametershave the benefit that they providean
on-hand estimate for the validity of our approximations
A¯=A/(λa µ), ω¯ =ω/λ, D¯ =D/ λ(a µλ)2 . (14) and can be evaluated directly from (1).
− −
(cid:0) (cid:1) The equation we have used in our simulations thus
For our purposes it is advantageous to use instead the reads (omitting the overbars)
equivalent set
x˙(t)= x(t)+A cos(ωt)+√2D ξ(t) (16)
−
U¯+/D¯, U¯−/D¯, and ω¯, (15)
withthethresholdlocatedatx=1. Forobtainingtheres-
idencetime,xhasbeenresetintotheminimumofU(x,t)
immediately after firing. The Figs. 1 and 2 depict the
probability densities of the first-passage and the resi-
FIG.1: First-passagetimedensityforparametersU+/D=8,
U−/D = 5, ω = 0.05. The jagged line depicts the his-
togram obtained from iterations of the Langevin equation FIG. 3: Testing extreme limits. First-passage time density
in (16). Note that fluctuations in the histogram depend on foran extremelysmall lowerbarrier U−/D=3. Theremain-
the total number of events and the width of the histogram ing parameters are as in Fig. 1. Langevin simulation results
bins. These fluctuations stay completely within their ex- (jagged), analytical result in (12) with eq. (9) (solid), like-
pected range which is indicated as the gray shaded area in wise, with eq. (10) (dashed). The inset compares the rate
the inset. The height of this area is twice the expected stan- κ(t) obtained from simulations of (16) with the analytic re-
dard deviation of the histogram levels. The solid line shows sults from eq. (9) and eq. (10), respectively.
theanalyticfirst-passagetimedensityfromeq.(12),withthe
rate used in (9). Both lines are in excellent agreement.
FIG. 4: Testing extreme limits. Probability density of the
first-passage time for a very fast driving, ω = 0.5. The
other parameters are as in Fig. 1. The analytic approxi-
mation (eqs. (12) and (9), dashed line) still depicts maxima
FIG. 2: Residencetime density vs. time for thesame param- that are approximately located at 1,2,... while the numeri-
eters as in Fig. 1. The jagged line shows the histogram ob- calresultsareshiftedtolatertimes. Thejagged linepresents
tainedfromiterationsof(16). Again,thenumericspractically the Langevin-iterations from (16). The solid curve presents
coincides within the line-width with the analytic expression within its linewidth both the numerical solution to the FP
in (13) evaluated with therate in (9) (solid line). equation (2) and to the integral equation from [10].
4
dence times, respectively. The residence time distribu- cannot be expected a priori.
tion exhibits a less pronounced modulation in compari-
Our method, primarily aimed at describing first-
sonwiththefirst-passagetimedistribution. Bothanalyt-
passage time and residence time probabilities of driven
ical expressions in (12) and (13) compare very favorably
dynamical systems, is also readily extended to more re-
with the numerical results obtained by iteration of the
alistic neuron models such as e.g. the two-dimensional
Langevinequation(16). Theremainingdeviationsinthe
driven FitzHugh-Nagumo model [14] for neuronal spik-
two figures are of purely statistical nature (see the inset
ing activity, whose multiple attractors may be consid-
in Fig. 1) and can be diminished further by increasing
ered as discrete states. Likewise, the scheme can also
the number of events in the simulations.
be employed to study yet other time-dependent switch-
In order to further test the range of validity of our
ing dynamics and synchronization phenomena such as
novelapproximationschemewe–onpurpose–havecho- the paradigmofStochastic Resonance [6]anddiscreteor
sen extreme values for the lower barrier height U−/D continuous Brownian motor transport [15].
andangulardrivingfrequencyω,respectively,seeFigs.3
ThisworkhasbeensupportedbytheDeutscheForsch-
and 4. Here, deviations from the numerical results are
ungsgemeinschaft via Projects No. HA1517/13-4 and
nottheresultofstatisticsbutaresystematic. Forthelow
SFB-486, Projects No. A10 and No. B13.
potential barrier in Fig. 3, the time-scales in the process
are not separated sufficiently. The fast intra-well fluctu-
ations begin to influence the behavior of the modulated
firing dynamics. Moreover, the difference between the
moderate noise result for the time-dependent rate κ(t) [1] L. Lapicque, J. Physiol. (Paris) 9, 620 (1907).
in (9) and its weak noise approximationin (10) becomes [2] H. C. Tuckwell, Stochastic Processes in the Neurosci-
visiblyincreased,asexpected. Figure4depictstheother ences, (SIAM,Philadelphia, 1989).
extreme situation with a modulation time-scale that is [3] B.W.Knight,J. Gen.Phys. 59,734(1972); M.Stemm-
not slow enough. Because the system cannot follow the ler, Network-Comp. Neural 7, 687 (1996); P. Lansky,
Phys.Rev.E55,2040(1997);M.T.GiraudoandL.Sac-
driving instantaneously we find a shift in the maxima
erdote, BioSystems 48, 77 (1998); T. Shimokawa et al.,
of the first-passage time density. This shift is not repro- Phys.Rev.E59,3461(1999);H.E.PlesserandT.Geisel,
ducedbyourapproximationin(5)and(12);nevertheless, Phys. Rev. E 59, 7008 (1999); M. J. Chacron et al.,
ourschemeyieldsamazinglygoodresultsevenwithinthis Phys. Rev.Lett. 85 1576 (2000); N. Brunel et al.,Phys.
extreme parameter regime. The results basedonthe nu- Rev. Lett. 86, 2186 (2001); L. Sacerdote and P. Lansky,
mericalevaluationoftheFPequation(2)andtheintegral BioSystems67,213(2002);J.W.Middletonetal.,Phys.
equation [10] virtually collapse into one curve and per- Rev. E 68, 021920 (2003); B. Lindner, A. Longtin and
A. Bulsara, Neural Comp. 15, 1761 (2003).
fectly coincide with the Langevinsimulations. The same
[4] A.R.Bulsara,T.C.Elston,C.R.Doering,S.B.Lowen,
result was obtained for the other parameter values. and K.Lindenberg, Phys. Rev.E 53, 3958 (1996).
Conclusions.—The precise theoretical modeling of the [5] H.E.PlesserandW.Gerstner,Neurocomputing 32–33,
neuronal spiking activity under external time-dependent 219(2000);H.E.PlesserandT.Geisel,Phys.Rev.E59,
driving presents a challenge of considerable importance 7008 (1999).
in neurophysiology and physics. Due to the presence of [6] L. Gammaitoni, P. H¨anggi, P. Jung, and F. Marchesoni,
Rev. Mod. Phys. 70, 223 (1998).
non-stationarity, absorbing boundary conditions of the
[7] P.H¨anggi,P.Talkner,andM.Borkovec,Rev.Mod.Phys.
underlyingfirst-passageproblemanddifferingtime-scales
62, 251 (1990).
the task of obtaining reliable analytical estimates for
[8] B. Lindner and L. Schimansky-Geier, Phys. Rev. Lett.
the firing statistics is anything but trivial. By refer- 86, 2934 (2001); N. Fourcaud and N. Brunel, Neural
ence to a discrete Markovian dynamics for the corre- Comp. 14, 2057 (2002); B. Lindner et al., Phys. Rep.
sponding full space-continuous stochastic dynamics we 392, 321 (2004).
succeededinobtaininganalyticalapproximationsforthe [9] P. Jung and P. H¨anggi, Phys. Rev. A 44, 8032 (1991);
V. A. Shneidman et al., Phys. Rev. Lett. 72, 2682
time-dependentfirst-passagetimeandtheresidencetime
(1994); N. G. Stocks, Nuovo Cimento D 17, 925 (1995);
statistics that are valid beyond the restraining limits of
J. Lehmann, P. Reimann, and P. H¨anggi, Phys. Rev.
linear response and asymptotically weak noise. We have Lett. 84, 1639 (2000); A. Nikitin, N. G. Stocks and A.
tested our findings for the case of a periodically driven R. Bulsara, Phys. Rev. E 68, 016103 (2003); J. Casado-
LIF model. The obtained agreement with precise nu- Pascual et al.,Phys.Rev. Lett.91, 210601 (2003).
merical simulations of either the Langevin type in (16) [10] E. Di Nardo, A.G. Nobile, E. Pirozzi, L. M. Ricciardi,
or,equivalently,oftheFPtypein(2)turnsouttobevery Adv. Appl. Probab. 33, 453 (2001); A. Buonocore, et
al., Adv.Appl. Probab. 19, 784 (1987).
good. Ourmethodisnotrestrictedtoanoscillatoryforc-
[11] P. Talkner and J. Luczka,arXiv:condmat/0307498.
ing but applies as well to arbitrary drive functions f(t)
[12] R. L¨ofstedt and S. N. Coppersmith, Phys. Rev. E 49,
such as an exponentially decaying drive (e.g. simulat-
4821 (1994); M. H. Choi, R. F. Fox, and P. Jung, Phys.
ing a decaying threshold). Our scheme even yields good Rev. E 57, 6335 (1998); P. Talkner, Physica A 325, 124
results in extreme parameter regimes where agreement (2003).
5
[13] J. Lehmann, P. Reimann, and P. H¨anggi, Phys. Rev. E [15] R. D. Astumian and P. H¨anggi, Phys. Today 55 No.11,
62, 6282 (2000). 33 (2002); P. Reimann, Phys.Rep. 361, 57 (2002).
[14] R. A. FitzHugh, Biophys. J. 1, 445 (1961); J. Nagumo,
et al.,Proc. IRE50, 2061 (1962).