Escape from an attractor generated by recurrent exit
EEscape from an attractor due to recurrent exit increases the escape time in Kramer’stheory of activation
Lou Zonca , and David Holcman Group of Computational Biology and Applied Mathematics,Ecole Normale Sup´erieure-PSL, 75005 Paris, France. Sorbonne University, Pierre et Marie Curie Campus,75005 Paris, France.
Kramer’s theory of activation over a potential barrier consists in computing the mean exit timefrom the boundary of a basin of attraction of a randomly perturbed dynamical system. Here wereport that for some systems, crossing the boundary is not enough, because stochastic trajectoriesreturn inside the basin with a high probability a certain number of times before escape far away.This situation is due to a shallow potential. We compute the mean and distribution of escape timesand show how this result explain the large distribution of interburst durations in neuronal networks.
Kramers theory [1–4] is the classical framework tostudy the escape time over a potential barrier: it consistsin computing the mean first passage time of a dynamicalsystem perturbed by a small noise to the boundary of abasin of attraction. The mean first passage time mea-sures the stability and provides great insight of the back-ward binding rate in chemistry [5], loss of lock for phasecontrollers in communication theory [6], escape of recep-tors from the post-synaptic density at neuronal synapse[7] and is used to evaluate future derivatives in the finan-cial market [8].In the limit of small noise, a trajectory escapes a basin ofattraction with probability one [9], but the escape timeis exponentially long depending on the topology of thenoiseless dynamics [10] and its behavior at the boundary.In addition, the distribution of exit points peaks at a dis-tance O ( √ σ ) from a saddle-point, where σ is the noiseamplitude [2, 11]. Interestingly, when a focus attractoris located near the boundary of the basin of attraction,the escape time deviates from an exponential distributionbecause trajectories oscillate inside the attractor beforeescape [12–16].In these previous examples, the escape ends at the firsttime a trajectory crosses the separatrix that delimits thebasin of attraction. We introduce here a class of shallowtwo-dimensional dynamical systems for which trajecto-ries first exit the basin of attraction, then make excur-sions outside before coming back inside the domain, abehavior that occurs several times before eventually es-caping far away. This situation is peculiar and theserecurrent entries need to be taken into account in com-puting the final escape time. This letter reports suchclass of systems and their escape time. We derive for-mulas for the mean and distribution of escape times andwe show that these recurrent reentries inside the basin ofattraction can multiply the escape time by a factor be-tween two and three. Finally, we discuss how this resultexplains long interbursts in neuronal networks. Appearance of a recurrent escape pattern.
We in- x h A B x h C P r obab ili t y
022 1h4 0 1.5x 10.5-1 0
A S x10 -2 D x h A C ГГ A SExit
A S Г FIG. 1.
Emergence of a of recurrent escape patternA.
Escaping trajectories reach the separatrix Γ for the firsttime (step 1, black) and cross the separatrix several timesgoing inside and outside of the basin of attraction (step 2,green, cyan, blue) before they escape (pink). B. Stochastictrajectories doing one (yellow) and two (orange) RT beforeescape. C. Distributions of successive exit points on Γ (500runs). D. Outer boundary layer C computed as the convexhull of all trajectories reentering the basin of attraction (red). troduce the 2D model˙ h = − αh + x + σ ˙ ω ˙ x = (cid:26) h − γx for h ≥ − γx for h ≤ , (1)where α ∈ ]0 , γ ∈ ]0 , α [, ˙ ω is a Gaussian white noiseand σ its amplitude. This system has two critical points:one attractor A = (0 ,
0) (fig. 1A yellow star) and onesaddle-point S = ( γ α, γα ) (fig. 1A cyan star) and theseparatrix (Γ) delimits the basin of attraction of A (fig.1A solid black).Stochastic escape of the basin of attraction occurs intwo steps. 1) A trajectory starting at A reaches Γ for thefirst time (fig. 1A black trajectory between A and the a r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p first exit point Exit S has only realeigenvalues λ ± = − (cid:16) − ( α + γ ) ± (cid:112) ( α + γ ) + 4 αγ (cid:17) , λ + ≈ . , λ − ≈ − . − ˜ p ) that we willcompute below. Characterization of the escape time
The escape timecan be decomposed into the time to reach the separatrixΓ for the first time plus the time spent for the trajec-tory to go back and forth around Γ before a final escape.Using Baye’s law and conditionning the time by the RTnumber, the mean escape time can be decomposed as (cid:104) τ esc (cid:105) = ∞ (cid:88) k =0 (cid:104) τ | k (cid:105) P RT ( k ) , (2)where (cid:104) τ | k (cid:105) (resp. P RT ( k )) is the mean time (resp. prob-ability) to return k times inside the basin of attraction.A trajectory has terminated its escape when it reaches asecond boundary C that delimits the region of the phase-space where trajectories still have a high probability ofreentering the basin of attraction, meaning that the es-cape process is not complete until the trajectory reachesthis new boundary. We approximate C as the convex hullof all the trajectories that have not escaped yet (fig. 1Dred, 500 runs). Each RT can be considered independentof the previous ones, thus the probability to escape after k RT is given by P RT ( k ) = ˜ p (1 − ˜ p ) k − , (3)and the mean escape time is (cid:104) τ esc (cid:105) = (cid:104) τ (cid:105) + ( (cid:104) τ ext (cid:105) + (cid:104) τ int (cid:105) )˜ p (cid:80) ∞ k =1 k (1 − ˜ p ) k − = (cid:104) τ (cid:105) + (cid:104) τ ext (cid:105) + (cid:104) τ int (cid:105) ˜ p . (4)where (cid:104) τ (cid:105) is the mean time to reach the separatrixfor the first time and (cid:104) τ ext (cid:105) (resp. (cid:104) τ int (cid:105) ) is the timespent on the outside (resp. inside) of the basin of at-traction of A at each RT (fig. 2A) Note that, if the BA τ RT = τ int +τ ext P r obab ili t y RT D u r a t i on ( s ) 〈 τ int 〉〈 τ ext 〉 RT P r obab ili t y
10 20 30 40 50 6000.050.10.15
Noise amplitude σ0.3 0.5 0.7 0.92468 R T 〈⟩ Number of RT before escapeRT durations
Time (s) D i s t r i bu t i on o f e x i t t i m e s Exit times for 0 RTExit times for 1 RT
Time (s) D i s t r i bu t i on o f e x i t t i m e s C D σ =
FIG. 2.
Distribution of RT and escape times. A.
Dis-tributions of the duration of a RT τ RT,k = τ ext,k + τ int,k for k ∈ [1 , (cid:104) τ ext,k (cid:105) (resp. (cid:104) τ int,k (cid:105) ) with respectto the RT number k . B. Distributions of the RT numberaround the separatrix before a trajectory escapes definitelyfor various values of σ (with γ = 0 . α = 1), 500 runsfor each value of σ . Inset: mean RT number with respect tothe noise amplitude σ . C. Distributions f (upper), resp. f (lower), of escape times for trajectories doing no RT, resp. 1RT, with the fit (8). D. Distribution of exit times with thecontribution of each RT number overlayed with the analyticalexit time distribution (equation 7). probability to escape directly ˜ p tends to zeros, the es-cape time tends to infinity which corresponds to thecase where the trajectory would be trapped going in-side and outside of the basin of attraction forever. Herewe obtain from our numerical simulations ˜ p ≈ .
12 thus (cid:104) τ esc (cid:105) ≈ (cid:104) τ (cid:105) + 8 . (cid:104) τ ext (cid:105) + (cid:104) τ int (cid:105) ). Note that with ourparameters (cid:104) τ (cid:105) ≈ . s and (cid:104) τ ext (cid:105) + (cid:104) τ int (cid:105) ≈ s meaningthat the escape time is increased 2.6 times. Interestinglythe noise amplitude does not influence the number of RTdone before escaping (fig. 2B), and for our value of theparameter γ = 0 .
6, trajectories do 8 RT on average (fig.2B, inset) which is coherent with our calculation of themean escape time. These results indicate that the noiseamplitude does not directly influence the probability toescape from the region contained between Γ and C , how-ever it could modify the distance between both curves.We now determine the distribution of escape times P ( τ esc < t ) = ∞ (cid:88) k =0 P ( τ k < t | k ) P RT ( k ) , (5) P ( τ k < t | k ) is the conditional distribution probability ofescape times after k RT. Given that the RT are i.i.d, thisprobability is the k -th convolution of the distribution oftimes of a single RT f ( t ) with the distribution of escapetimes without RT f ( t ) P ( τ k < t | k ) = f ( t ) ∗ f ( t ) ∗ k , (6)where f ( t ) ∗ k = f ( t ) ∗ f ( t ) ∗ ... ∗ f ( t ), k times. Thus thepdf of exit times is given by f ( t ) = ∞ (cid:88) k =0 f ( t ) ∗ f ( t ) ∗ k ˜ p (1 − ˜ p ) k − . (7)To compare this formula to the results of our numericalsimulations, we approximate the distributions f and f by a function of the form f i ( t ) = c i (cid:32) erf (cid:32) t − a i b i (cid:33)(cid:33) e − λ i t, for i = 0 , , (8)where erf ( x ) = 2 √ π (cid:90) x e − u du is the error function. Wefitted the distributions obtained from the numerical sim-ulations of the trajectories that escape without doingany RT ( f fig. 2C, upper) and after one single RT ( f fig. 2C, lower) conditioned to λ ≥ λ . We obtained c = 1 . c = 1 . λ = 0 . λ = 0 . a = − . a = − . b = − . b = − .
69. We then com-puted each term of the sum (7) and we could compare itto the corresponding parts of the distribution of escapetimes obtained from our numerical simulations (fig. 2D).
Application to interburst durations for a firingneuronal network
Burst and interburst durations arefundamental in neuronal rhythm generation. Interspikevariability has been modeled using a general class of neu-ronal models [17]. Network synchronization, which is fun-damental for burst generation, depends on such intervals[18]. However, the mechanisms leading to long interburstintervals are still under investigation. Long interburstshave been modeled using a two-state synaptic depression[19], or by accounting for the refractory period induced byAHP [20]. Here we show that recurrent escape patternscan also explain long interbursts intervals. We apply theprevious results to the depression-facilitation short-termsynaptic plasticity model of network neuronal bursting[21, 22]. This is a mean-field model which consists ofthree equations for the mean voltage h , the depression y ,and the facilitation x : τ ˙ h = − h + Jxyh + + √ τ σ ˙ ω ˙ x = X − xt f + K (1 − x ) h + (9)˙ y = 1 − yt r − Lxyh + , where h + = max ( h,
0) is a linear threshold functionof the synaptic current that gives the average popula-tion firing rate [23]. The mean number of connections (synapses) per neuron is accounted for by the parame-ter J [24] and the term Jxy represents the effect of theshort-term synaptic plasticity on the network activity.The parameters K and L describe how the firing rate istransformed into molecular events that are changing theduration and probability of vesicular release. The timescales t f and t r define the recovery of a synapse from thenetwork activity. Finally, ˙ ω is an additive Gaussian noiseand σ its amplitude, it represents fluctuations in the fir-ing rate.This system has 3 critical points, one attractor and twosaddles. Around the attractor A = (0 , X,
1) the dynam-ics are very anisotropic ( | λ | = 12 . (cid:29) | λ | = 1 . (cid:29)| λ | = 0 .
34) and we can project it on a 2D-plan y =constant:˙ y = 0 = 1 − yτ r − Lxyh + = 0 ⇐⇒ y = 11 + τ r Lxh + (10)the 2D simplified dynamics is˙ h = h ( Jx − − τ r Lxh + ) τ (1 + τ r Lxh + ) + √ τ σ ˙ ω ˙ x = X − xτ f + K (1 − x ) h + (11)This 2D deterministic system (for σ = 0) has 3 criticalpoints, two attractors and one saddle-point . Attractor A A first equilibrium point is given by h =0 and x = X . The Jacobian at this point is J A = − JXτ K (1 − X ) − τ f . (12)With our parameters (Table I) the eigenvalues λ = JX − τ ≈ − . λ = − τ f ≈ − .
11 are both nega-tive confirming A is an attractor. Saddle-point S The second critical-point is S ( h ≈ . x ≈ . λ ≈ − .
73 and λ ≈ .
43. It is a saddle-point.
Attractor A The third critical-point is A ( h ≈ . x ≈ . λ ≈ − . λ ≈ − .
33. It is another attractor. The two attractorsare separated by the 1D stable manifold of the saddle-point S (fig. 3A, solid black curve).The phase-space of system (11), restricted to { x ≤ . h ≤ } has the same topological properties as sys-tem (1): one attractor and one saddle-point, the sep-aratrix delimiting the basin of attraction is the stablemanifold of S (fig. 3A). The escaping trajectories exitsand re enters the basin of attraction several times beforeeventually escaping (fig. 3A, orange).Interburst intervals correspond to the exit times of thetrajectories of system (11) from the basin of attraction, A BC
Facilitation x F i r i ng r a t e h Time (s) D i s t r i bu t i on o f e x i t t i m e s Number of RT before escape P r obab ili t y R T 〈⟩ RT D u r a t i on ( s ) 〈 τ int 〉〈 τ ext 〉 D A S Γ FIG. 3.
Application to the dynamical system (11) . A.
2D phase-space restricted to { x ≤ . h ≤ } . The basinof attraction of A (yellow star) is delimited by the stablemanifold of S (solid black curve Γ) with an exiting trajectorydoing 1 RT (orange) around the separatrix before escape . B. Distribution of exit times with the contribution of the tra-jectories per RT number before escape with the analytical fit(equations 7, 13 and 14). C. Distribution of the RT numberfor σ ∈ [4 ,
7] and mean RT number with respect to noise (in-set). D. Values of (cid:104) τ ext (cid:105) (red) and (cid:104) τ int (cid:105) (black) with respectto the RT number. for trajectories starting close to the attractor. We canthus use formula (7) to fit the distribution of exit timesobtained. In this case, we have ˜ p ≈ .
13 and we obtain(fig. 3B) f ( t ) = 0 .
23 exp( − . t ) (cid:32) erf (cid:32) t − . . (cid:33)(cid:33) (13)and f ( t ) = 0 .
19 exp( − . t ) (cid:32) erf (cid:32) t + 15 . . (cid:33)(cid:33) . (14)Finally, we note that, as for the generic system (1), theRT number before escape does not depend on the noiseamplitude (fig. 3C), the trajectories do on average 8 RTbefore escape (inset). We use formula (4) to determinethe mean escape time (cid:104) τ esc (cid:105) ≈ (cid:104) τ (cid:105) + 7 . (cid:104) τ ext (cid:105) + (cid:104) τ int (cid:105) )where (cid:104) τ (cid:105) ≈ . s and (cid:104) τ ext (cid:105) + (cid:104) τ int (cid:105) ≈ . s (fig. 3D)thus multiplying the escape time by a factor 2.2. Conclusion and discussion
We described a new es-cape pattern for which reaching the boundary of the de-terministic basin of attraction is not sufficient to escape.We derived analytical formulas for both the mean escapetime and the distribution of escape times taking into ac-count the excursions inside and outside of the basin of at-traction before the final escape. We show that this type
Parameters Values τ Time constant for h J Synaptic connectivity 4.21 [20] K Facilitation rate 0.037Hz [20] X Facilitation resting value 0.08825 [20] L Depression rate 0.028Hz [20] τ r Facilitation time rate 2.9s [20] τ f Depression time rate 0.9s [20] T Depolarization parameter 0TABLE I. Model (9) parameters of escape pattern can apply to neuronal networks andexplain the long interburst durations observed in somecases. [1] H. A. Kramers, Physica , 284 (1940).[2] Z. Schuss, Theory and Applications of Stochastic Differ-ential Equations (Wiley, 1980).[3] Z. Schuss,
Theory and Applications of Stochastic Pro-cesses: An Analytical Approach. (Springer New York,2010).[4] W. C. Gardiner and A. Burcat,
Combustion chemistry (Springer, 1984).[5] A. Nitzan,
Chemical dynamics in condensed phases: re-laxation, transfer and reactions in condensed molecularsystems (Oxford university press, 2006).[6] Z. Schuss,
Nonlinear filtering and optimal phase tracking ,Vol. 180 (Springer Science & Business Media, 2011).[7] D. Holcman and Z. Schuss,
Asymptotics of Elliptic andParabolic PDEs (Springer, Cham, 2019).[8] J.-P. Fouque, G. Papanicolaou, and K. R. Sircar,
Deriva-tives in financial markets with stochastic volatility (Cam-bridge University Press, 2000).[9] B. J. Matkowsky and Z. Schuss, SIAM Journal on Ap-plied Mathematics , 365 (1977).[10] M. I. Freidlin and A. D. Wentzell, in Random perturba-tions of dynamical systems (Springer, 1998) pp. 15–43.[11] B. Bobrovsky and Z. Schuss, SIAM Journal on AppliedMathematics , 174 (1982).[12] T. Verechtchaguina, I. M. Sokolov, and L. Schimansky-Geier, Physical Review E , 031108 (2006).[13] T. Verechtchaguina, I. Sokolov, and L. Schimansky-Geier, EPL (Europhysics Letters) , 691 (2006).[14] T. Verechtchaguina, I. Sokolov, and L. Schimansky-Geier, Biosystems , 63 (2007).[15] K. Dao Duc, Z. Schuss, and D. Holcman, MultiscaleModeling & Simulation , 772 (2016).[16] K. D. Duc, Z. Schuss, and D. Holcman, Physical ReviewE , 030101 (2014).[17] B. S. Gutkin and G. B. Ermentrout, Neural computation , 1047 (1998).[18] B. Ermentrout, M. Pascal, and B. Gutkin, Neural com-putation , 1285 (2001).[19] C. Guerrier, J. A. Hayes, G. Fortin, and D. Holcman,Proceedings of the National Academy of Sciences ,9728 (2015).[20] L. Zonca and D. Holcman, arXiv preprintarXiv:2001.11432 (2020). [21] M. V. Tsodyks and H. Markram, Proc. Natl. Acad. Sci.USA , 719 (1997).[22] K. Dao Duc, C.-Y. Lee, P. Parutto, D. Cohen, M. Segal,N. Rouach, and D. Holcman, Plos one , e0124694 (2015).[23] D. Holcman and M. Tsodyks, PLoS Computational Bi-ology , 174 (2006).[24] E. Bart, S. Bao, and D. Holcman, Journal of Computa-tional Neuroscience19