Table Of ContentEscape rates in periodically driven Markov processes
6 Michael Schindler, Peter Talkner , Peter Ha¨nggi
∗
0
0
2 Institutfu¨rPhysik,Universita¨tAugsburg, D–86135 Augsburg, Germany
n
a
J
9 Abstract
]
n We present an approximate analytical expression for the escape rate of time-dependent
a
drivenstochastic processes withanabsorbing boundary suchasthedrivenleakyintegrate-
-
a and-fire model for neural spiking. The novel approximation is based on a discrete state
t
a Markovian modeling of the full long-time dynamics with time-dependent rates. It is valid
d
inawideparameterregimebeyondtherestraining limitsofweakdriving(linearresponse)
.
s
and/orweaknoise.Theschemeiscarefullytestedandyieldsexcellentagreementwiththree
c
i different numerical methods based on the Langevin equation, the Fokker-Planck equation
s
y andanintegralequation.
h
p
[
Keywords: absorbing boundary, non-stationary Markovprocesses, Ornstein-Uhlenbeck
1 process, rateprocess, periodic driving, drivenneuronmodels,
v
PACS:05.40.-a, 05.10.Gg,02.70.Rr,82.20.Uv,89.75.Hc
5
4
0
1
0
6
0
/ 1 Introduction
s
c
i
s
y
Although the solution of the stationary and unbounded Ornstein-Uhlenbeck pro-
h
p cess has been found long ago, it is not yet possible to give an analytic exact ex-
v: pressionthatincludestime-dependentdrivingand absorbingboundaries[1,2]. Yet,
i such processes with a linear restoring force an a periodic driving which termi-
X
nate at a prescribed threshold are widely used as models for numerous physical
r
a
effects. Examples range from rupturing experiments on molecules [3] where the
time-dependence is introduced as linear movement of the absorbing boundary up
tototallydifferentmodelsliketheleakyintegrate-and-fire(LIF)modelforneuronal
spikingevents[4,5,6,7,8].Thelatteristheapplicationweprimarilythinkofinthis
paper.Thestochasticvariablestandsforthecellsoma’selectricpotentialx(t)thatis
changing due to a great many incoming signals from other neurons. It is thus cus-
tomary to employ a diffusion approximation for the stochastic dynamics of x(t).
[email protected](PeterTalkner)
∗
PreprintsubmittedtoElsevierScience 2February2008
The driven abstract LIF model assumes the non-stationary Langevin dynamics (in
dimensionlesscoordinates)
x˙(t) = x(t)+Acos(ωt+φ)+√2Dξ(t) (1)
−
where the process starts at time t = 0 at x(0) = x and fires when it reaches the
0
threshold voltage x = a 1. ξ(t) is white Gaussian noise. Here, a sinusoidal
≡
stimulushas been chosen for thesake ofconvenience. Thefollowinganalysismay
easily beextended to general periodicstimuli.Thedynamicsof theprocess x(t) is
equivalentlydescribedbyaFokker-Planck(FP)equationfortheconditionalproba-
bilitydensityfunction(PDF)ρ(x,t x ,0)inatime-dependentquadraticpotential,
0
2 |
U(x,t) = x Acos(ωt+φ) /2,reading
(cid:16) − (cid:17)
∂ ∂ ∂2
′
ρ = L(t)ρ = U (x,t)ρ +D ρ, (2)
∂t ∂x (cid:16) (cid:17) ∂x2
withtheabsorbingboundaryand initialconditions
ρ(a,t x ,0) = 0 forall tand x (3)
0 0
|
ρ(x,0 x ,0) = δ(x x ). (4)
0 0
| −
After firing the process immediately restarts at the instantaneous minimum of the
potential.
The set of eqs. (1–4) defines our starting point for obtaining the firing statistics of
this driven neuron model. Our main objective is to develop an accurate analytical
approximationthat avoids certain restrictiveassumptionsof prior attempts.Those,
infact,allinvolvetheuseofeitherofthefollowinglimitingapproximationschemes
such as the limit of linear response theory (i.e. a weak stimulus A 1) [7,9]
≪
or the limit of asymptotically weak noise [10,11,12,13,14]. Our scheme detailed
below yields novel analytic and tractable expressions beyond the linear response
and weak noise limit; as will be demonstrated, this novel scheme indeed provides
analytical formulae that compare very favorably with precise numerical results of
the full dynamics in eqs. (1, 2–4). The arguments given for the agreement of the
first-passage time distribution also hold for the residence-time [15] which is not
furtherconsidered here.
2 Reduction to adiscrete model
TheperiodicityoftheexternaldrivingwiththeperiodT = 2π/ωallowsonetorep-
resent the time-dependent solution ρ(x,t) of the Fokker Planck equation in terms
of Floquet eigenfunctions and eigenvalues of the Fokker-Planck operator, vi(x,t)
2
and µi,respectively,[10,16]
∂
vi(x,t)+L(t)vi(x,t) = µivi(x,t), (5)
−∂t
where the eigenfunctions are periodic in time, integrable in x from to a, and
−∞
fulfilltheabsorbingboundaryconditionat x = a
vi(a,t) = 0. (6)
The time-dependent PDF can be written as a weighted sum of the Floquet eigen-
functions
ρ(x,t) = civi(x,t) exp(µit) (7)
Xi
wherethecoefficientsci aredeterminedbytheinitialPDF.Notethatbecauseofthe
absorbing boundary condition at x = a the total probability is not conserved and
thereforeall Floqueteigenvalueshaveanon-vanishingnegativereal part.
The first main assumption that we impose concerns the value of the potential at
the boundary: The minimum of the potential must always belong to the “allowed”
regionleftofthethreshold,and,moreover,thepotentialdifferencebetweenthresh-
oldand minimum,denotedby ∆U(t), mustalways belarger thanat least afew D,
i.e. ∆U(t)/D > 4. This assumption implies an exponential time-scale separation
between the average time τκ in which the threshold is reached from the minimum
ofthepotentialcompared to thetimeτr ofthedeterministicrelaxation towardsthe
potential minimum. In the dimensionless units used here τr = 1. For the Floquet
spectrum this implies the presence of a large gap between the first eigenvalue µ
1
which is of the same order as 1/τκ and the higher ones which are of the order
−
1orsmaller.Afterashortinitialtimeoftheorder1,allcontributionsfromhigher
−
Floqueteigenvaluescan beneglectedand onlythefirst onesurvives:
ρ(x,t) v (x,t) exp(µ t) (8)
1 1
≈
In general, the Floquet eigenfunctionsand the correspondingeigenvalues are diffi-
culttodetermine.Aformalexpansionintermsoftheinstantaneouseigenfunctions
ψi(x,t) ofL(t) fulfilling
L(t)ψk(x,t) = λk(t)ψk(x,t) (9)
isalways possiblethoughnotalways helpful
vi(x,t) = dik(t)ψk(x,t). (10)
Xk
The periodicity of vi(x,t) and ψk(x,t) implies that the coefficients dik(t) also are
periodicfunctionsof time.Expansion(10), togetherwith theFloquet equation (5),
yieldsacoupledsetofordinarydifferentialequationsforthecoefficientsdik(t)[17]
3
∂
d˙ik(t) (λk(t) µi)dik(t) = dil(t) ϕk(t),ψl(t) , (11)
− − (cid:28)∂t (cid:29)
Xl
where ϕk(x,t) denotes the instantaneous eigenfunction of the backward opera-
torL+(t) belongingto theeigenvalueλk(t)
L+(t)ϕk(x,t) = λk(t)ϕk(x,t). (12)
Theeigenfunctionsψk(x,t)andϕk(x,t)constituteabi-orthogonalsetoffunctions
thatalways can benormalizedsuch that
ϕl(t),ψk(t) = δkl. (13)
h i
Here, the scalar product f,g is defined as theintegral overthereal axis up to the
h i
threshold:
a
f,g = dxf(x)g(x) (14)
h i Z−∞
Withoursecond assumptionwerequirethat thedrivingfrequency ω is smallcom-
pared to the relaxation rate in the parabolic potential. Under this condition, the
matrix elements ∂ϕk(t)/∂t,ψl(t) that are proportional to the frequency ω are
h i
also small and may be neglected to lowest order in the equations for the coeffi-
cientsdik(t)[17].Theresultingequationsareuncoupledandreadilysolvedtoyield
withtheperiodicboundaryconditions
t
′ ′
d1k(t) δ1kexp dt λ1(t) µ1t , (15)
≈ (cid:18)Z0 − (cid:19)
where µ = 1 T λ (t)dt follows from the periodicity of d (t). Together with
1 T 0 1 11
eqs. (8)and(10)Rweobtainforthelong-timebehaviorofthePDF
t
′ ′
ρ(x,t) exp dt λ (t) ψ (x,t). (16)
1 1
≈ (cid:18)Z0 (cid:19)
Note,thatthefirstFloqueteigenvaluehascanceled.Thelowestinstantaneouseigen-
functionsψ (x,t) and ϕ (x,t) are related by
1 1
ψ (x,t) = ϕ (x,t)ρ (x,t), (17)
1 1 0
where
ρ (x,t) exp U(x,t)/D . (18)
0
∝ (cid:16)− (cid:17)
Forthecorrespondingeigenvalueλ (t) wefind from(9)
1
a
dxϕ (x,t)L(t)ϕ (x,t)ρ (x,t)
−∞ 1 1 0
λ (t) = . (19)
1 R a dxϕ2(x,t)ρ (x,t)
−∞ 1 0
R
An explicit expression, valid for high potential differences, can be given after lin-
earization ofU about a
′
ϕ (x,t) = 1 exp (x a)U (a,t)/D (20)
1
− (cid:16) − (cid:17)
4
whichgivesforλ (t)
1
∆U(t) 1 erf ∆U(t)/D
λ1(t) = − (cid:16)q (cid:17) . (21)
− D 1 exp( ∆U(t)/D)
− −
whereerf(z) istheerrorfunction.
Thewaiting-timeprobability[18]can beexpressedas
a t
′ ′
P(t) = dxϕ (x,t)ρ(x,t) = exp dt λ (t) . (22)
1 1
Z−∞ (cid:16)Z0 (cid:17)
Therefore, the eigenvalue λ (t) coincides with the negative of the time-dependent
1
escaperate κ(t).
With the expression(21) for the escape rate we can calculate the property of inter-
est, namely thePDF forthefirst-passagetime(FPT) oftheattracting ”integrating”
state that covers the domain < x(t) < a. The FPT-PDF is given by the nega-
−∞
tiverate ofchangeofthewaitingtimeprobability,i.e.
dP(t) t
′ ′
g(t) = = κ(t) exp κ(t)dt , (23)
− dt (cid:18)−Z0 (cid:19)
The quantitative validity of these expressions for an extended parameter regime
willbechecked next.
3 Numerical analysis
We implemented three different numerical methods to obtain both the FPT-PDF
and therate in order to have a reliablebasis for comparisonwith the analytical ex-
pression (21). The first method performs explicit time-steps of the Langevin equa-
tion (1). We used an elaborate technique for the time-integrationof the fluctuating
force ξ(t). For points away from the threshold a it is sufficient to take a normal
distributed random variable for the displacement due to ξ(t). Quite the contrary in
thevicinityoftheabsorbingboundary.Here,theintegralofξ(t)ratherbehaveslike
aWienerprocesswithabsorbingboundary,as illustratedinFig.1.Theappropriate
transition distribution, is known analytically as the weighted difference between
twonormal distributions[1]
p(xi+1,t+δt xi,t) = N1(xi+1,δt xi,0) N2(xi+1,δt xi,0)
| | − | (24)
=: N1(xi+1,δt xi,0) 1 Pout(xi+1,xi,δt)
| (cid:16) − (cid:17)
The multiplication on the right-hand side stands for a logical AND that leads to a
correction step in the algorithm. First, a new position xi+1 is proposed according
5
to the normal distribution density N1. With the probability Pout(xi+1,xi,δt) the
trajectoryhasalreadycrossedtheboundaryduringthistime-stepδtfromxi toxi+1
and, therefore, isto beended. TheexplicitformsofN and N give
1 2
N 1
2
Pout(xi+1,xi,δt) = = exp (a xi+1)(a xi) . (25)
N (cid:18)−Dδt − − (cid:19)
1
Thesameformulahas been givenby[19]butwith adifferentreasoning.
In order to get the correctly normalized FPT-PDF g(t) we counted the number of
trajectories hittingthe absorbing boundary withintheinterval [t,t+δt). The FPT-
PDF is then estimated by this number divided by δt and by the total number of
trajectories.Therate isgivenby
κ(t) = g(t)/P(t), (26)
xi a
Fig. 1. Thetransition probabilities p(xi+1,t+δt xi,t) (black line) and N1(xi+1,t+δt
| |
xi,t)(dashedline)fromxiforasingletime-step,withandwithouttheabsorbingboundary,
respectively. Theverticallineindicates theboundary.
0.1 0.1
g(t)/ω g(t)/ω
0.01 0.01
0.001 0.001
0.2
κ(t)/ω 0.1 κ(t)/ω
0.15
0.075
0.1
0.05
0.05 0.025
0 0
0 2 4 6 8 10 0 2 4 6 8 10
t/T t/T
Fig. 2. First-passage time density g(t) (upper) and rate κ(t) (lower plots) as functions
of time. Displayed are all three numerical methods we used for testing (solid lines)
and the approximation based on (21) (dashed lines). The parameters in the left plots,
∆U(t)/D [5,8] and ω = 0.05, are chosen to yield a very good approximation of the
∈
rate by (21). The right plots display extreme parameters, ∆U(t)/D [3,8] and ω = 0.5,
∈
wherenogoodapproximation oftheratecanbeexpected. Inbothcasesφ = 0.
6
where 1 P(t) is estimated by the number of trajectories that have escaped up to
−
timet, dividedby thetotalnumberoftrajectories.
For the second numerical method we have solved the FP equation (2) using a
Chebychev collocation method to reduce the problem to a coupled system of or-
dinarydifferentialequations[13,20].ThisgivesP(t)astheintegralofρ(x,t)from
to a. The FPT-PDF in thefigures is then calculated according to eq. (23), and
−∞
therateagainby (26).
The third method solves Ricciardi’s integral equation for the FPT-PDF and is de-
tailedin[21,22].Foremployinghisalgorithmtheprocessmustbetransformedinto
astationaryOrnstein-Uhlenbeckprocesswith amovingabsorbingboundary
A
−t
S(t) = a cos(ωt+φ)+ωsin(ωt+φ) e . (27)
− 1+ω2(cid:20) − (cid:21)
Allthreemethodsprovidepracticallyidenticalresultsascan beseeninFigs.2and
3. The results for the FPT-PDF and for the rate all collapse into one single line.
Differences between the numerical methods, e.g. fluctuations in the histogram of
the Langevin equation method are visible only in the plots of the relative errors
(Fig. 3 middleand lowerrows).
Figure 2 shows that the FPT-PDF is extremely well approximated by expression
(21) for the rate κ(t). In the left plots we used quite a high barrier with quite slow
driving compared to the time-scale τr of the process. Good agreement is thus ex-
pected.Intherightplotsweshowthesituationwithextremeparameters.Thelower
barrierheight∆U /D goesdownto3where arate-descriptionisunlikelytosuf-
min
fice. Moreover,thedrivingisfaster,ω = 0.5.Thesystemcannotfollowthedriving
instantaneously,andwefindashiftinthemaximumoftheFPT-PDFtolatertimes.
Under these conditions it is impressing how good the novel approximation still
works.
A moredelicatemeasure fortheerrors oftheapproximationare therate κ(t) itself
and its relative deviation from the three numerically calculated rates. Both can be
seen in Fig. 3. The upper row of plots shows the approximation error of the rate
for the same two parameter sets as in Fig. 2. Especially at the maximum the rate
is over-estimated. This leads to a faster decay of the FPT-PDF which is scarcely
visible in Fig. 2. Also, the shift of the maxima (indicated by vertical lines) can be
observed.It isnegligiblysmallfor ω = 0.05but morepronouncedfor ω = 0.5.
In the middle row of Fig. 3 a systematic error of the approximation becomes vis-
ible. The relative error with respect to the numerical results behaves roughly si-
nusoidally with a phase-shift of π/2 relative to the driving and with an additional
constant offset. For the instantaneous rate expression (19) to be valid it is neces-
sary that the driving signal is sufficiently slow. If this assumption is violated then
a rate can still be defined if the barrier is sufficiently high. But in addition to the
7
leading term d (t) in (10) the higher instantaneous eigenfunctions must be taken
11
into account. The coupling to the coefficients d1k(t) is induced by the matrix ele-
ments φ˙k(t),ψl(t) ,see eq. (11), containinga timederivativethat introducesnon-
h i
adiabaticcorrectionsto therateand, consequently,to thestatisticsoftheFPT.
It is quiteastonishing,that the huge relativeerror in the right middleplot of Fig. 3
leads to such a good result in Fig. 2. The explanation for this is that the FPT-
PDF (23)uses the time-integratedrate. Therefore, errors are importantonly where
therateislarge.Acloserlookontheplotshowsthataroundthemaximaoftherate
therelativeerror is comparably small. Because theerrors are linear in timearound
the rate’s maxima they cancel out when integrated over time in (23). The same is
2 2.5 3 3.5 4 2 2.5 3 3.5 4
0.2 0.1
κ(t)/ω 0.08 κ(t)/ω
0.15
0.06
0.1
0.04
0.05
0.02
0 0
0.1 κα−κ 1 κα−κ
κα κα
0.05 0.5
0 0
0.05 0.5
− −
0.1 1
− −
0.03 0.03
0.02 κα−κβ 0.02 κα−κβ
κα+κβ κα+κβ
0.01 0.01
0 0
0.01 0.01
− −
0.02 0.02
− −
0.03 0.03
− 2 2.5 3 3.5 4 − 2 2.5 3 3.5 4
t/T t/T
Fig.3.Comparison ofthe numerical rates andthenovel approximation from eq.(21). The
respective parameters in the left/right plots are the same as in Fig. 2. Upper plots: The
numerically determined rates are displayed as solid lines: Langevin equation simulations
(black); Fokker-Planck equation (blue); Ricciardi’s integral equation (green). Thetheoret-
ical approximation κ(t) from eq. (21) is displayed as the red dashed line. Middle plots:
Relative error of the approximation κ(t) with respect to each numerical rate κα(t) (with
thesamecolorcodingasabove).Lowerplots:Errorsofthenumericalrateswithrespectto
eachother. Thethinverticallinesindicate thepositions ofthenumericalrates’maxima.
8
0 0
0.025 0.025
− −
0.05 0.05
− −
0.075 κRIC κ 0.075 κRIC κ
− − − −
κRIC κRIC
0.1 0.1
− −
0.008 0.08
0.006 0.06
0.004 0.04
Differencebetweenthepositions Differencebetweenthepositions
0.002 ofthemaximaofκRIC andκ 0.02 ofthemaximaofκRIC andκ
inunitsofthedrivingperiodT inunitsofthedrivingperiodT
0 0
0 5 10 15 20 25 30 0 5 10 15 20 25
∆Umin/D ∆Umin/D
Fig.4.Relativeerrorandrelativetime-shiftoftherate(21).Asthebasisofcomparisonwe
used the rate κRIC obtained by solving Ricciardi’s integral equation [21] at times t > τr
whereithasbecomeperiodic. Upperplots:Theerrorrelativetoκ evaluatedatthemax-
RIC
ima of κ . Shown are data for (∆U /D ∆U /D) (0.1,1,2,3,5), from top to
RIC max min
− ∈
bottom with the colors (red,green,blue,cyan,magenta), and the phase φ (0,π/2,π)
∈
with the symbols (×, , +). In the left panel the driving is slow, ω = 0.05, in the right
◦
it is fast, ω = 0.5. Note that the relative error is of the same order of magnitude for slow
and for fast driving. A dependence on the phase φ cannot be observed. Lower plots: The
differenceofthemaxima’spositionofκ andrate(21)inunitsoftheperiodT,againfor
RIC
ω = 0.05(leftpanel)andω = 0.5(rightpanel).Colorandsymbolcodingsarethesameas
above.
0
−0.2 0.15 DofiffthereemncaexbimetawoefenκRthICeapnodsitκions
inunitsofthedrivingperiodT
0.4 0.1
−
κRIC κ
−
0.6 κRIC 0.05
−
0.8 0
− 0.01 0.1 1 0.01 0.1 1
ω ω
Fig. 5. Relative error and relative time-shift of the maxima as a function of ω. For all
data points ∆U /D = 5. The color and symbol codings for ∆U /D and φ are the
min max
same asinFig. 4. Forfast driving the relative error starts depending on theinterval length
of∆U(t)/D.
9
validfortheresidence timewhosePDF alsocontainsintegralsoftherate[15,23].
Figure4showsthisrelativeerrorofκ(t)atthemaximaofthenumericallyobtained
rateasafunctionofthebarrierheight.Again,twodifferentdrivingfrequenciesare
given. In both cases the relative error has the same order of magnitude, and thus
explainswhyboth parametersetsinFig. 2 yieldgoodapproximations.
Finally,wewouldliketo pointthereader’s attentiontothelimitationsofthelinear
response approximation. For linear response the parameter ratio A/D needs to be
small. In our validating example in Fig. 2 (left plots) it takes on the value A/D =
1.5. Thus,ourapproximationschemeisvalidbeyondthelinearresponselimit.
The time-scale of the driving force is mainly restricted by the relaxation time-
scale τr and much less by the magnitude of the rate itself. There is no restriction
ontherelativemagnitudesofκandT.Instead,bothω andκhavetobesufficiently
small.Fig.5indicatesthatboththerelativeerrorandthetime-shiftofthemaxima’s
positionsaremodestforω < 0.1.
4 Conclusions
By reference to a discrete Markovian dynamics for the corresponding full space-
continuous stochasticprocess we succeeded in obtaining an analytical approxima-
tion for the time-dependent escape rate which can be used for calculating first-
passage time statistics. This result is valid beyond the restraining limits of linear
responseorasymptoticallyweak noiseand ofadiabaticallyslowdriving.
We checked our findings using simulations of the Langevin equation (1) and nu-
merical solutions of the equivalent FP equation in (2) and of the integral equation
in [21]. We found an impressiveagreement for thefirst-passagetimedensityand a
goodmatch fortherate whichisthemoredelicateproperty forcomparison.
Finally, we note that our method is not restricted to a periodic forcing but applies
also to arbitrary drive functions. However, in the oscillatory case some of the ap-
proximation errors cancel out. This leads to useful results in extreme parameter
regimeswhereagreement cannot beexpected a priori.
ThisworkhasbeensupportedbytheDeutscheForschungsgemeinschaftviaproject
HA1517/13-4and SFB-486, projectsA10 andB13.
10
