Asymptotics of extreme statistics of escape time in 1,2 and 3-dimensional diffusions
Kanishka Basnayake, Claire Guerrier, Zeev Schuss, David Holcman
aa r X i v : . [ q - b i o . S C ] O c t Asymptotics of extreme statistics of escape time in 1,2and 3-dimensional diffusions
K. Basnayake ∗ , C. Guerrier † , Z. Schuss ‡ , D. Holcman ∗ Abstract
The first of N identical independently distributed (i.i.d.) Brownian trajectoriesthat arrives to a small target, sets the time scale of activation, which in generalis much faster than the arrival to the target of only a single trajectory. Analyti-cal asymptotic expressions for the minimal time is notoriously difficult to computein general geometries. We derive here asymptotic laws for the probability densityfunction of the first and second arrival times of a large number of i.i.d. Browniantrajectories to a small target in 1,2, and 3 dimensions and study their range of va-lidity by stochastic simulations. The results are applied to activation of biochemicalpathways in cellular transduction. Fast activation of biochemical pathways in cell biology is often initiated by the first arrivalof a particle to a small target. This is the case of calcium activation in synapses of neuronalcells [1–3], fast photoresponse in rods, cones and fly photoreceptors [4–6], and many more.However, the time scale underlying these fast activations is not very well understood. Wepropose here that these mechanisms are initiated by the first arrival of one or more of themany identical independently distributed (i.i.d.) Brownian particles to small receptors(such as the influx of many calcium ions inside a synapse to receptors).In general, one or several particles are required to initiate a cascade of chemical reac-tions, such as the opening of a protein channel [7] of a cellular membrane, which amplifiesthe inflow of ions to an avalanche of thousand or more ions, resulting from binding ofcouple of ions. These statistic of the minimal arrival times are referred to in the statis-tical physics literature as extreme statistics [8]. Despite great efforts [8–16], there are noexplicit expressions for the probability distributions of arrival times of the first trajectory, ∗ Applied Mathematics & Computational Biology, Ecole Normale Sup´erieure, 46 rue d’Ulm 75005Paris, France. † The University of British Columbia Mathematics Department, Vancouver, B.C. Canada V6T 1Z2 ‡ Department of Applied Mathematics, Tel-Aviv University, Tel-Aviv 69978, Israel. N particles. Asymptotic expressions for the mean escape time of asingle Brownian path, the so called narrow escape time, computed in the narrow escapetheory [17, 18], depends on global and local geometric properties of the bounding domainand its boundary, such as the surface area (in 2 dimensions or volume (in 3 dimensions),and the local geometry near the absorbing window (mean curvature of the boundary atthe small window, the window’s shape, and relative size). The number and distribution ofabsorbing windows can influence drastically the narrow escape time. As shown below, theescape of the fastest particles selects trajectories that are very different from the typicalones, which determine the mean narrow escape time (NET).Moreover, the asymptotics of the expected first arrival time are not the same as ofthe mean first passage time (MFPT) of a single Brownian path to a small window. Theanalysis relies on the time-dependent solution of the Fokker-Planck equation and theshort-time asymptotics of the survival probability. Previous studies of the short-timeasymptotics of the diffusion equation concern the asymptotics of the trace of the heatkernel analysis [19, 20]. Here, an estimate is needed of the survival probability, whichrequires different analysis. Our method is based on the construction of the asymptoticsfrom Green’s function of the Helmholtz equation.The main results are explicit expressions for the statistics of the first arrival time seeattempt in [14]) in 1,2, and 3 dimensions and a formula for the expected shortest exit timefrom a neuronal spine with and without returns. The manuscript is organized as follows.First we present the general framework for the computation of the pdf of the first arrivaland conditional second arrival, given that the first one has already arrived, in a populationof N Brownian particles in a ray and in an interval. We then study the difference betweenPoissonian and diffusion escape time approximations and in particular, we consider thecase of a bulk domain with a window connected to a narrow cylinder (dendritic spineshape [17, 18]). We then compute the pdf of the extreme escape time in dimensions 2and 3 through small windows. Finally, we discuss applications to activation in cellularbiology.
The narrow escape problem (NEP) for the shortest arrival time of N non-interacting i.i.d.Brownian trajectories (ions) in a bounded domain Ω to a binding site is defined as follows.Denote by t i the arrival times and by τ the shortest one, τ = min( t , . . . , t N ) , (1)2here t i are the i.i.d. arrival times of the N ions in the medium. The NEP is to find thePDF and the MFPT of τ . The complementary PDF of τ is given byPr { τ > t } = Pr N { t > t } , (2)where Pr { t > t } is the survival probability of a single particle prior to binding at thetarget. This probability can be computed from the following boundary value problem. As-suming that the boundary ∂ Ω contains N R binding sites ∂ Ω i ⊂ ∂ Ω ( ∂ Ω a = N R S i =1 ∂ Ω i , ∂ Ω r = ∂ Ω − ∂ Ω a ), the pdf of a Brownian trajectory is the solution of the initial boundary valueproblem (IBVP) ∂p ( x , t ) ∂t = D ∆ p ( x , t ) for x ∈ Ω , t > p ( x ,
0) = p ( x ) for x ∈ Ω ∂p ( x , t ) ∂ n =0 for x ∈ ∂ Ω r p ( x , t ) =0 for x ∈ ∂ Ω a . The survival probability is Pr { t > t } = Z Ω p ( x , t ) d x , (4)so that Pr { τ = t } = ddt Pr { τ < t } = N (Pr { t > t } ) N − Pr { t = t } , (5)where Pr { t = t } = I ∂ Ω a ∂p ( x , t ) ∂ n dS x . (6)Pr { t = t } = N R I ∂ Ω ∂p ( x , t ) ∂ n dS x . (7)Putting all the above together results in the pdfPr { τ = t } = N N R Z Ω p ( x , t ) d x N − I ∂ Ω ∂p ( x , t ) ∂ n dS x . (8)The first arrival time is computed from the survival probability of a particle and the fluxthrough the target. Obtaining an explicit or asymptotic expression is not possible ingeneral. 3 .1 The pdf of the first arrival time in an interval To obtain an analytic expression for the pdf of the first arrival time (8) of a particleinside a narrow neck that could represent the dendritic spine neck, we model the narrowspine neck as a segment of length L , with a reflecting boundary at x = 0 and absorbingboundary at x = L . Then the diffusion boundary value problem (3) becomes ∂p∂t = D ∂ p∂x for 0 < x < L, t > p ( x,
0) = δ ( x ) for 0 < x < L (10) p ( L, t ) = ∂p (0 , t ) ∂x = 0 for t > , (11)where the initial condition corresponds to a particle initially at the origin. The generalsolution is given by the eigenfunction expansion p ( x, t ) = 2 ∞ X n =0 e − Dλ n t cos λ n x, (12)where the eigenvalues are λ n = πL (cid:18) n + 12 (cid:19) . (13)The survival probability (4) of a particle is thus given byPr { t > t } = L Z p ( x, t ) dx = 2 ∞ X n =0 ( − n λ n e − Dλ n t . (14)The pdf of the arrival time to L of a single Brownian trajectory is the probability effluxat the absorbing boundary ∂ Ω a , given by − I Ω a ∂p ( x , t ) ∂ n dS x = − ∂p ( L, t ) ∂x = 2 ∞ X n =0 ( − n λ n e − Dλ n t . (15)Therefore, the pdf of the first arrival time in an ensemble of N particles to one of N R independent absorbers is given byPr { τ (1) = t } = 2 N N R ∞ X n =0 ( − n λ n e − Dλ n t ! N − ∞ X n =0 ( − n λ n e − Dλ n t . (16)For numerical purposes, we approximate (16) by the sum of n terms,Pr { τ (1) = t } ≈ f n ( t ) = N N R n X n =0 ( − n λ n e − Dλ n t ! N − n X n =0 ( − n λ n e − Dλ n t . (17)4igs.1A-B show the pdf of the first arrival time for N = 5 and N = 500 Brownian particleswith diffusion coefficient D = 1, which start at x = 0 at time 0 and exit the interval at x = 1. These figures confirm the validity of the analytical approximation (16) with only n = 100 terms in the lowly converging alternating series.Figure 1: Histograms of the arrival times to the boundary of the fastest particle, obtained fromBrownian simulations with Euler’s scheme. The number of Brownian particles is N = 5 in A and N = 500 in B . The analytical solution (red curves) is obtained by setting n = 100 in (17). Next, we turn to the computation of the conditional pdf of the arrival time τ (2) of thesecond particle, which is that of the minimum of the shortest arrival time in the ensembleof N − τ (1) . The time τ (1) + τ (2) is that of arrival of the first two particles at reach the target. Theconditional distribution of the arrival time τ (2) of the second particle, given the positions ofthe N − τ (1) , can be computed from their joint probability distributionat positions ( x , . . . , x N − ) and the first particle has already arrived at time τ (1) = s ,Pr { τ (2) = t } (18)= Z t Z Ω .. Z Ω Pr { τ (2) = t, τ = s, x ( s ) = x , . . . , x N − ( s ) = x N − } dx · · · dx N ds and Pr { τ (2) = t, τ = s, x ( s ) = x , . . . , x N − ( s ) = x N − } = Pr { τ (2) = t | τ = s, x ( s ) = x , . . . , x N − ( s ) = x N − }× Pr { τ = s } Pr { x ( s ) = x , . . . , x N − ( s ) = x N − } . Because all particles are independent,Pr { x ( s ) = x , . . . , x N − ( s ) = x N − } = N − Y i =1 Pr { x i ( s ) = x i }
5o thatPr { τ (2) = t } = t Z Pr { τ (2) = t | τ = s } L Z Pr { x ( s ) = x } dx N Pr { τ = s } ds. (19) The pdf (19) can be evaluated under some additional assumptions. For example, ifthe Brownian trajectories escape from a deep potential well, the escape process is wellapproximated by a Poisson process with rate equal the reciprocal of the mean escape timefrom the well [21]. Also, when Brownian particles escape a domain Ω = B ∪ C , whichconsists of a bulk B and a narrow cylindrical neck C , the escape process from Ω can beapproximated by a Poisson process, according to the narrow scape theory [17]. Here themotion in the narrow cylinder C is approximated by one-dimensional Brownian motionin an interval of length L .Consequently, under the Poisson approximation, the arrival of the first particle is muchfaster than the escape of the second one from the bulk compartment B , thus we can usethe approximation that all particles are still in the bulk B after the arrival of the firstone. The bulk is represented by the position x = 0 in an approximate one-dimensionalmodel. This assumption simplifies (19) to L Z Pr { x ( s ) = x } dx N ≈ , (20)so that Pr { τ (2) = t } = t Z Pr { τ (2) = t | τ = s } Pr { τ = s } ds. (21)The Markovian property of the Poisson process givesPr { τ (2) = t | τ = s } = Pr { τ (2) = t − s, } (22)so that τ (2) has the same pdf as τ with N − N , that is, Pr { τ (2) = t } = t Z f ( t − s ) f ( s ) ds, (23)where (recall (16)) f ( s ) = Pr { τ = s } = N N R g ( t ) N h ( t ) . (24)6n one dimension, g ( t ) = P N t n =0 ( − n λ n e − Dλ n t and h ( t ) = P N t n =0 ( − n λ n e − Dλ n t . It followsthat Pr { τ (2) = t } ≈ N N R t Z g ( s ) N − h ( s ) g ( t − s ) N − h ( t − s ) ds. (25) Pr { τ (2) } of N Brownian i.i.d. trajectories in a segment
As in the first paragraph of section 3, equations (18) and (19) are valid with Ω replacedby the segment [0 , L ]. That is,Pr { τ (2) = t } = t Z L Z · · · L Z Pr { τ (2) = t, τ = s, x ( s ) = x , . . . , x N ( s ) = x N } dx · · · dx N ds. and Pr { τ (2) = t, τ = s, x ( s ) = x , . . . , x N − ( s ) = x N } = Pr { τ (2) = t | τ = s, x ( s ) = x , . . . , x N − ( s ) = x N }× Pr { τ = s } Pr { x ( s ) = x , . . . , x N − ( s ) = x N } . Because all particles are independent,Pr { x ( s ) = x , . . . , x N ( s ) = x N − } = N − Y i =1 Pr { x i +1 ( s ) = x i } , (26)hence,Pr { τ (2) = t } = t Z Pr { τ (2) = t | τ = s } L Z Pr { x ( s ) = x } dx N − Pr { τ = s } ds. (27)To compute the survival probability S ( s ) = L Z Pr { x ( s ) = x } dx , (28)we use the short-time asymptotics of the one-dimensional diffusion equation. Modifyingequation (9) for short-time diffusion of a particle starting at 0, we get ∂p ( x, t ) ∂t = D ∂ p ( x, t ) ∂x for x > , t > p ( x,
0) = δ ( x ) for x > , p ( L, t ) = 0 for t > . (29)7he short-time diffusion is well approximated by the fundamental solution (except at theboundary, where the error is exponentially small in 1 /t ) p ( x, t ) = 1 + o ( t ) √ Dπt exp (cid:26) − x Dt (cid:27) . (30)Thus the survival probability at short time t is S ( t ) = L Z o ( t ) √ Dπt exp (cid:26) − x Dt (cid:27) dx. (31)The short-time asymptotic expansion (42) (see below) and the change of variable x = u √ Dt in the integral (31), give S ( t ) = 1 − √ π ∞ Z L/ √ Dt (cid:2) exp (cid:8) − u (cid:9)(cid:3) du (32) ≈ − √ Dt exp n − ( L/ √ Dt ) o √ πL (cid:18) − DtL + O (cid:18) t L (cid:19)(cid:19) . (33)It follows from (27) that the pdf of the second arrival time isPr { τ (2) = t } (34)= [1 + o (1)] t Z Pr { τ = s } Pr { τ = t − s } − √ Ds exp n − ( L √ Ds ) o √ πL N − ds. Figure 2 compares results of the stochastic simulations with the analytical formula (25)for the second fastest arrival time τ (2) to the boundary 1 of the interval [0 ,
1] among 20particles. We use the analytical formula (25) (no correction) and (34), which contains theshift correction due to the distribution of the particles in the interval at time τ (1) , whenthe first particle arrives at x = L been absorbed.This figure shows how the corrected formula gives a better agreement with the Brow-nian simulations, thus proving that the distribution of the Brownian particles inside theinterval contributes to the decrease of the arrival time of the second particle. The fit tosimulation data is based on the eigenfunction expansionPr { τ (2) = t } = [1 + o (1)] t Z Pr { τ = s } Pr { τ = t − s } ∞ X n =0 ( − n λ n e − Dλ n t ! N − ds, (35)which is equivalent to (34). Formula (24) is used for Pr { τ = s } . Note that the alternatingseries contains an even number of terms. 8 r { τ ( ) = t } N=20 corrected with gN=20 uncorrected
Time ( a.u.)
Figure 2:
Histogram of the arrival time of the second fastest particle, obtained from Browniansimulations with Euler’s scheme. The fastest is computed for N = 20 in B. The analyticalsolution with no correction is given by (22) (blue) and compared to (34) with the correction(red). There are n = 6 terms in the series (17). To conclude, the internal distribution of particle is given by − √ Ds exp n − ( L √ Ds ) o √ πL N − and causes the faster arrival of the second particle relative to the first one. This exampleshows the deviation from the purely Poissonian approximation. ¯ τ N i.i.d. Brownian paths is given by¯ τ = ∞ Z Pr { τ > t } dt = ∞ Z [Pr { t > t } ] N dt, (36)where t is the arrival time of a single Brownian path. Writing the last integral in (36) as¯ τ = ∞ Z e N ln g ( t ) dt, (37)9t can be expanded for N ≫ g ( t ) = ∞ X n =0 ( − n λ n e − Dλ n t (38)(see (14)). Consider the case L = ∞ and the IBVP ∂p ( x, t ) ∂t = D ∂ p ( x, t ) ∂x for x > , t > p ( x,
0) = δ ( x − a ) for x > , p (0 , t ) = 0 for t > , (39)whose solution is p ( x, t ) = 1 √ Dπt (cid:20) exp (cid:26) − ( x − a ) Dt (cid:27) − exp (cid:26) − ( x + a ) Dt (cid:27)(cid:21) . (40)The survival probability with D = 1 isPr { t > t } = ∞ Z p ( x, t ) dx = 1 − √ π ∞ Z a/ √ t e − u du. (41)To compute the MFPT in (36), we use the expansion of the complementary error function2 √ π ∞ Z x e − u du = e − x x √ π (cid:18) − x + O ( x − ) (cid:19) for x ≫ , (42)which gives I N ≡ ∞ Z [Pr { t > t } ] N dt ≈ ∞ Z exp ( N ln − e − ( a/ √ t ) ( a/ √ t ) √ π !) dt, (43)and with the approximation I N ≈ ∞ Z exp ( − N √ te − a t a √ π ) dt = a ∞ Z exp ( − N √ ue − u √ π ) du. (44)To evaluate the integral (44), we make the monotone change of variable w = w ( t ) = √ te − /t , w ′ ( t ) = √ te − t (cid:18) t + 1 t (cid:19) . (45)10ote that for small t , w ′ ( t ) ≈ w t (46)and ln w ≈ − /t . Thus, w ′ ( t ) ≈ w (ln w ) . (47)Breaking with N ′ = N √ π I N ≈ a ∞ Z exp {− N ′ w } dwdt dw ≈ a δ Z exp {− N ′ w } a w (ln( w )) dw + ∞ Z δ exp {− N ′ w } dwdt dw for some 0 < δ <
1, the second integral turns out to be exponentially small in N and isthus negligible relative to the first one. Integrating by parts, I N ≈ a δ Z exp {− N ′ w } w (ln( w )) dw ≈ O (exp( − aN )) + a N ′ δ Z exp {− N ′ w } | w | dw and changing the variable to u = N ′ w , we obtain N ′ δ Z exp {− N ′ w } a | w | dw = N ′ δ Z a exp {− u } | ln u/N ′ | du. Expanding 1 | ln u/N ′ | = 1ln N ′ | ln u | ln N + O (cid:18) | ln u | ln N ′ (cid:19) ! for u > ε >
0, we obtain, N ′ δ Z exp {− N ′ w } a w dw ≈ N ′ δ Z exp {− u } a | ln N ′ | (cid:18) | ln u | ln N ′ (cid:19) du. Thus, breaking the integral into two parts, from [0 , ε ] (which is negligible) and [ ε, ∞ [, weget ¯ τ ≈ a D ln N √ π for N ≫ . (48)11 .2 Escape for the second fastest from half a line Equation (23) and the approximation (20) givePr { τ (2) = t } = t Z f ( t − s ) f ( s ) ds, (49)where f ( s ) = Pr { τ = s } = N N R g ( t ) N h ( t ) . (50)According to (41), g ( t ) = Pr { t > t } = ∞ Z p ( x, t ) dx = 1 − √ π ∞ Z a/ √ t e − u du ≈ − e − ( a/ √ t ) ( a/ √ t ) √ π (51)and h ( t ) = − d Pr { t > t } dt = 2 √ π a √ t e − ( a/ √ t ) . (52)The MFPT of the second arrival is obtained from (49) as¯ τ (2) = ∞ Z t Pr { τ (2) = t } dt = ∞ Z t t Z f ( t − s ) f ( s ) dsdt (53)= 2 ∞ Z sf ( s ) ds ∞ Z f ( t ) dt = 2¯ τ (1) ∞ Z f ( t ) dt = 2¯ τ (1) . (54) [0 , a ] We follow the steps of the previous section, where Green’s function for the homogenousIBVP is now given by the infinite sum p ( x, t | y ) = 1 √ Dπt ∞ X n = −∞ (cid:20) exp (cid:26) − ( x − y + 2 na ) t (cid:27) − exp (cid:26) − ( x + y + 2 na ) t (cid:27)(cid:21) . (55)12he conditional survival probability isPr { t > t | y } = a Z p ( x, t | y ) dx (56)= 1 √ Dπt ∞ X n = −∞ a Z (cid:20) exp (cid:26) − ( x − y + 2 na ) t (cid:27) − exp (cid:26) − ( x + y + 2 na ) t (cid:27)(cid:21) dx = a Z √ Dπt (cid:20) exp (cid:26) − ( x − y ) t (cid:27) − exp (cid:26) − ( x + y ) t (cid:27)(cid:21) dx + S ( y, t ) − S ( y, t ) , where S = 1 √ Dπt ∞ X n =1 a Z (cid:20) exp (cid:26) − ( x + y + 2 na ) t (cid:27) − exp (cid:26) − ( x − y + 2 na ) t (cid:27)(cid:21) dxS = 1 √ Dπt ∞ X n =1 a Z (cid:20) exp (cid:26) − ( x + y − na ) t (cid:27) − exp (cid:26) − ( x − y − na ) t (cid:27)(cid:21) dx. Note that the integrand in the third line of (29), denoted p ( x, t | y ), satisfies the initialcondition p ( x, | y ) = δ ( x − y ) and the boundary condition p (0 , t | y ) = p ( x, t |
0) = 0,but p ( a, t | y ) = 0 and p ( x, t | a ) = 0. However, with the first correction, p ( x, t | y ) = 1 √ Dπt (cid:20) exp (cid:26) − ( x − y ) t (cid:27) − exp (cid:26) − ( x + y ) t (cid:27) (57)+ exp (cid:26) − ( x − y − a ) t (cid:27) − exp (cid:26) − ( x + y − a ) t (cid:27) + exp (cid:26) − ( x − y + 2 a ) t (cid:27) − exp (cid:26) − ( x + y + 2 a ) t (cid:27)(cid:21) it satisfies the same initial condition for x and y in the interval, and the boundary condi-tions p ( x, t | a ) = 0 , p ( x, t |
0) = 1 √ Dπt (cid:20) exp (cid:26) − ( x + 2 a ) t (cid:27) − exp (cid:26) − ( x − a ) t (cid:27)(cid:21) . (58)Higher-order approximations correct the one boundary condition and corrupt the other,though the error decreases at higher exponential rates.The first line of (57) gives the approximation a Z √ πt (cid:20) exp (cid:26) − ( x + y ) t (cid:27) − exp (cid:26) − ( x − y ) t (cid:27)(cid:21) dx = 1 √ π y/ √ t Z ( y − a ) / √ t e − u du ∼ − max 2 √ t √ π " e − y / t y , e − ( a − y ) / t a − y as t → , (59)13igure 3: A. Plot of Pr { τ (1) = t } (escape from an interval) for N = 3 , , , n = 100 terms in the series of (17). B. Decay of the expected arrival time of the fastest particlevs N (red points). The plot of the asymptotic formula (65) is (blue) with parameter . N . where the maximum occurs at min[ y, a − y ] for 0 < y < a (the shortest ray from y to theboundary). Starting at x = a/
2, this givesPr { t > t } = a Z √ πt (cid:20) exp (cid:26) − ( x − a/ t (cid:27) − exp (cid:26) − ( x + a/ t (cid:27)(cid:21) dx as t → , so changing x + a/ z √ t in the first integral and x − a/ z √ t in the second, weget1 √ π a/ √ t Z − a/ √ t e − z dz − √ π a/ √ t Z a/ √ t e − z dz ≈ − √ te − a / t a √ π − √ te − a / t a √ π − √ te − a / t a √ π = 1 − √ te − a / t a √ π − √ te − a / t a √ π . (60)The second integral in the second line of (57) is I / = − a Z √ πt (cid:20) exp (cid:26) − ( x − a/ t (cid:27)(cid:21) dx. (61)Set x − a/ − z √ t , then (61) becomes I / = − √ π a/ √ t Z a/ √ t exp (cid:8) − z (cid:9) dz = − √ π ∞ Z a/ √ t exp (cid:8) − z (cid:9) dz + 1 √ π ∞ Z a/ √ t exp (cid:8) − z (cid:9) dz. − a Z √ πt (cid:20) exp (cid:26) − ( x − a/ t (cid:27) − exp (cid:26) − ( x − a/ t (cid:27)(cid:21) dx ≈ − √ ta √ π e − a / t + 2 √ t a √ π e − a / t , (62)and in the third line of (57), we get a Z √ πt (cid:20) exp (cid:26) − ( x + 3 a/ t (cid:27) − exp (cid:26) − ( x + 5 a/ t (cid:27)(cid:21) dx ≈ √ t a √ π e − a / t − √ t a √ π e − a / t , (63)hence ∞ Z [Pr { t > t } ] N dt ≈ ∞ Z exp (cid:26) N ln (cid:18) − √ ta √ π e − a / t (cid:19)(cid:27) dt (64)and the expected time of the fastest particle that starts at the center of the interval is(48) with a replaced by a/ N replaced by 2 N . That is,¯ τ ≈ a D ln N √ π for N ≫ . (65)Figure 3A shows plot of the pdf analytical approximation of shortest arrival time (17)with n = 100 terms, D = 1 and L = 1 for N = 4 ,
6, and 10. As the number of particlesincreases, the mean first arrival time decreases (Fig.3B) and according to equation (65),the asymptotic behavior is given by C/ log N, where C is a constant. We show below thepdf of the fastest Brownian particle. R , N i.i.d. Brownian particles in a boundeddomain Ω ⊂ R , , we assume that the particles are initially injected at a point y ∈ Ω andthey can escape through a single small absorbing window ∂ Ω a in the boundary ∂ Ω of thedomain. The pdf of the fist passage time to ∂ Ω a is given by (8).15 .1 Asymptotics in dimension 3 To determine the short-time asymptotics of the pdf, we use the Laplace transform ofthe IBVP (3) and solve the resulting elliptic mixed Neumann-Dirichlet BVP [22]. TheDirichlet part of the boundary consists of N well-separated small absorbing windows, ∂ Ω a = S Nj =1 ∂ Ω j and the reflecting (Neumann) part is ∂ Ω r = ∂ Ω − ∂ Ω a , so that the IBVP(3) has the form ∂p ( x , t | y ) ∂t = D ∆ p ( x , t | y ) (66) p ( x , | y ) = δ ( x − y ) for x , y ∈ Ω ∂p ( x , t | y ) ∂ n =0 for x ∈ ∂ Ω r p ( x , t | y ) =0 for t > , x ∈ ∂ Ω a . We consider the case N = 1. The Laplace transform of (3),ˆ p ( x , q | y ) = ∞ Z p ( x , t | y ) e − pt dt, (67)gives the BVP − δ ( x − y ) + q ˆ p ( x , q | y ) = D ∆ˆ p ( x , q | y ) for x , y ∈ Ω ∂ ˆ p ( x , q | y ) ∂ n =0 for x ∈ ∂ Ω r ˆ p ( x , q | y ) =0 for x ∈ ∂ Ω a . Green’s function for the Neumann problem in Ω is the solution of − ∆ x ˆ G ( x , q | y ) + q ˆ G ( x , q | y ) = δ ( x − y ) for x , y ∈ Ω , (68) ∂ ˆ G q ( x , q | y ) ∂n x =0 for x , y ∈ ∂ Ω . The asymptotic solution of (68) in R is given byˆ G ( x , q | y ) = e −√ q | x − y | π || x − y || + R q ( x , y ) , (69)where R q ( x , y ) is more regular than the first term. When x (or y ) is on the boundary,ˆ G ( x , q | y ) = e −√ q | x − y | (cid:18) π || x − y || + H ( x )2 π log | x − y | + R ( x , y ) (cid:19) , (70)16here R ( x , y ) is more regular than the logarithmic term and H ( x ) is a geometric factor[22]. Using Green’ identity, we obtain that Z Ω h ˆ p ( x , q | y )∆ x ˆ G ( x , q | y ) − ∆ x ˆ p ( x , q | y ) ˆ G ( x , q | y ) i d y = Z ∂ Ω " ˆ p ( x , q | y ) ∂ ˆ G ( x , q | y ) ∂n x − ∂ ˆ p ( x , q | y ) ∂n x ˆ G ( x , q | y ) dS y , hence ˆ p ( x , q | y ) = ˆ G ( x , q | y ) − Z ∂ Ω a ∂ ˆ p ( x , q | y ′ ) ∂n x ˆ G ( x , q | y ′ ) dS y ′ . (71)If the absorbing window ∂ Ω a is centered at x = A , then, for x ∈ ∂ Ω a ,0 = ˆ G ( A , q | y ) − Z ∂ Ω a ∂ ˆ p ( A , q | y ′ ) ∂n x ˆ G ( A , q | y ′ ) dS y ′ . (72)This is a Helmoltz equation and the solution is given by [22]ˆ p ( A , q | y ) = C √ a − r , (73)where r = | A − y | . Thus, C is computed from0 = G q ( A , y ) − Z ∂ Ω a ∂ ˆ p ( A , q, y | y ) ∂n x G q ( A , y ) dS y (74)and to leading order, G q ( A , y ) = Z ∂ Ω a Ce −√ q | y − A | √ a − r (cid:18) π || y − A || + H ( x )2 π log | y − A | + R ( y , A ) (cid:19) dS y . (75)If ∂ Ω a is a disk of radius a , then G q ( A , y ) ≈ C Z ∂ Ω a e −√ qr √ a − r πr πrdr = π I ( √ qa ) − L ( √ qa )) C, (76)where I is the modified Bessel function of the first kind and L the Struve function.Thus,ˆ p ( x , q | y ) = G q ( x , y ) − G q ( A , y ) 2 π (cid:0) I ( √ qa ) − L ( √ qa ) (cid:1) Z ∂ Ω a G q ( x , y ) dS y √ a − r . (77)17or | A − x | ≫ a and q large, I ( √ qa ) − L ( √ qa ) ≈ π √ qa ˆ p ( x , q | y ) ≈ G q ( x , y ) − G q ( A , y ) G q ( A , x ) 2 π ( I ( √ qa ) − L ( √ qa )) Z ∂ Ω a dS y √ a − r , hence, for a small circular window of radius a ˆ p ( x , q | y ) ≈ G q ( x , y ) − π √ qa G q ( A , y ) G q ( A , x ) + o ( a ) for a ≪ . (78)To leading order in small t and x , y ∈ Ω, we obtain for the first term using the leadingorder term in expression 69, L − ( G q ( x , y )) ≈ πt ) / e − | x − y | t . (79)We will now use the inverse Laplace transform [23][p.1026; 29.3.87] L − ( √ q e −√ q | x − y | | x − y | ) = 14 √ πt e − | x − y | t H ( | x − y | √ t ) , (80)where H ( x ) = 4 x − A ∈ ∂ Ω, the imagecharge (for the Dirichlet boundary) leads to a factor 1 / G q ( A , y ) G q ( A , x )2 π √ qa = √ qa e −√ q ( | A − y | + | A − x | )2 π | A − y || A − x | (81)and L − ( G q ( A , y ) G q ( A , x ) √ qa ) = a (4 πt ) e − ( | A − y | + | A − x | ) t | A − y || A − x | H ( | A − y | + | A − x | √ t ) . (82)Finally, L − (ˆ p ( x , q | y )) = 1 p (4 πt ) e − | x − y | t − a | A − y || A − x | e − ( | A − y | + | A − x | ) t H ( | A − y | + | A − x | √ t ) . (83) The short-time asymptotics of the survival probability with δ = | A − y | , S ( t ) ≈ Z Ω p t ( x , y ) d x (84)= 1 p (4 πt ) Z Ω e − | x − y | t − a | A − y || A − x | e − ( | A − y | + | A − x | ) t H ( | A − y | + | A − x | √ t ) d x = I ( t ) − I ( t ) − I ( t ) − I ( t ) , (85) H ( | A − y | + | A − x | √ t ) ≈ ( | A − y | + | A − x | ) t , (86)and I ( t ) = 1 p (4 πt ) Z Ω e − | x − y | t d x (87) I ( t ) = 1 p (4 πt ) a δt Z Ω | A − x | e − ( δ + | A − x | ) t d x (88) I ( t ) = 1 p (4 πt ) a t Z Ω e − ( δ + | A − x | ) t d x . (89) I ( t ) = 1 p (4 πt ) a δt Z Ω | A − x | e − ( δ + | A − x | ) t d x . (90)Each integral is evaluated in the short-time limit. I ( t ) = 1 p (4 πt ) Z Ω e − | x − y | t d x ≈ − √ π ∞ Z R a / √ t e − u du (91) ≈ − √ t e − ( R a / √ t ) R a √ π (cid:18) O (cid:18)(cid:16) R a / √ t (cid:17) (cid:19)(cid:19) , (92)where R a is the radius of the maximal ball inscribed in Ω. The integral I is evaluated bythe change of variables z = x − A and then η = ( δ + r ) / √ t , where r = | z | (recall that A is in Ω a ), I ( t ) = 1 p (4 πt ) a δt Z Ω | A − x | e − ( δ + | A − x | ) t d x (93)= 1 p (4 πt ) a δt Z Ω+ A | z | e − ( δ + | z | ) t π | z | d | z | (94) ≈ π p (4 πt ) a δt δ + R √ t Z δ √ t e − η ( √ tη − δ ) √ tdη, (95)19here R is the radius of the largest half-ball centered at A ∈ Ω a and inscribed in Ω. Forshort time, δ + R √ t Z δ √ t e − η dη ≈ √ tδ e − (cid:18) δ √ t (cid:19) (cid:18) − t δ + 12 t δ (cid:19) (96) δ + R √ t Z δ √ t ηe − η dη = 12 e − (cid:18) δ √ t (cid:19) − e − (cid:18) δ + R √ t (cid:19) ≈ e − (cid:18) δ √ t (cid:19) . (97)Therefore, I ( t ) ≈ p (4 πt ) πa δt te − (cid:18) δ √ t (cid:19) − δ √ t √ tδ e − ( δ √ t ) (1 − t δ + 12 t δ ) (98) ≈ a δ √ π √ t (1 − tδ ) e − (cid:18) δ √ t (cid:19) . (99) Now, I ( t ) = 1 p (4 πt ) a t Z Ω e − ( δ + | A − x | ) t d x (100)= 2 π p (4 πt ) a t δ + R √ t Z δ √ t e − η ( √ tη − δ ) √ tdη (101)= 2 π p (4 πt ) a t δ + R √ t Z δ √ t e − η (4 tη − √ tηδ + δ ) √ tdη, (102)that we write as I ( t ) = I (1)3 ( t ) + I (2)3 ( t ) + I (3)3 ( t ). The approximation δ + R √ t Z δ √ t e − η η dη ≈ δ √ t e − (cid:18) δ √ t (cid:19) + 14 √ tδ e − (cid:18) δ √ t (cid:19) (cid:18) − t δ + o ( t ) (cid:19) , (103)20ives I (1)3 ( t ) = 2 πδπ / a t √ t e − { δ √ t } + 2 a tπ / δ √ te − ( δ √ t ) (cid:18) − t δ (cid:19) I (2)3 ( t ) = − πa δt √ tπ / e − ( δ √ t ) I (3)3 ( t ) = 2 πa δt √ tπ / e − (cid:18) δ √ t (cid:19) (cid:18) − t δ + o ( t ) (cid:19) . Summing the three contributions, the leading order terms cancel and I ( t ) = 4 a √ tπ / δ e − (cid:18) δ √ t (cid:19) . (104)To compute I , we decompose it into 4 pieces: I ( t ) = 1 p (4 πt ) a δt Z Ω | A − x | e − ( δ + | A − x | ) t d x . (105)= 2 π p (4 πt ) a δt δ + R √ t Z δ √ t e − η ( √ tη − δ ) √ tdη (106)= J ( t ) + J ( t ) + J ( t ) + J ( t ) . (107)Direct computations give: J ( t ) = 2 π (4 t ) p (4 πt ) a δt δ + R √ t Z δ √ t e − η η dη = 4 a δ √ π (4 t ) / (1 + 4 tδ ) e − (cid:18) δ √ t (cid:19) , (108)where we used δ + R √ t Z δ √ t η e − η dη ≈ (cid:18) δ t (cid:19) e − (cid:18) δ √ t (cid:19) . (109)Next, J ( t ) = − π p ( π ) a δt δ + R √ t Z δ √ t e − η η δdη = − a δ √ π (4 t ) / (cid:18) tδ (1 − tδ + 12 t δ + o ( t )) (cid:19) e − (cid:18) δ √ t (cid:19) , (110) δ + R √ t Z δ √ t η e − η dη ≈ δ √ t + √ t δ (cid:18) − t δ + 12 t δ + o ( t ) (cid:19)! e − (cid:18) δ √ t (cid:19) (111)Using relation 96, J ( t ) = 2 π p (4 tπ ) a δt δ + R √ t Z δ √ t e − η ηδ dη = 12 a δ √ π (4 t ) / e − δ t . (112)Finally, J ( t ) = − π t p ( π ) a δt δ + R √ t Z δ √ t e − η δ dη = − a δ (4 t ) / √ π e − δ t (cid:18) − tδ + 12 t δ + o ( t ) (cid:19) . (113)Direct computations show that the terms in t − / and t − / cancels out in the computationof I from the four terms J , ..J and it remains only the term in t / I ( t ) = − a √ πδ t / e {− δ / t } . (114)Summing (91)-(104)-(105), we get S ( t ) = Z Ω p t ( x , y ) d x = 1 − √ t e − ( R a / √ t ) R a √ π − a δπ / √ t e − δ / t + o t / e − (cid:18) δ √ t (cid:19) ≈ − a δπ / √ t e − δ / t .
22t follows that in three dimensions, the expected shortest arrival time to a small circularwindow of radius a , the expected shortest time ¯ τ is given by¯ τ = ∞ Z [Pr { t > t } ] N dt ≈ ∞ Z exp N log (cid:18) − a δπ / √ t e − δ / t (cid:19) dt ≈ ∞ Z exp (cid:18) − N a/δ ) δπ / √ t e − δ / t (cid:19) dt ≈ δ ∞ Z exp (cid:18) − N ′ √ u e − / u (cid:19) du, where N ′ = N a π / δ Using the method develop in section 4.1 with the change of variable, w = w ( t ) = 1 √ t e − /t , w ′ ( t ) = 1 √ t e − t (cid:18) − t + 14 t / (cid:19) . (115)We have with w ′ = 4 w (log( w )) / ¯ τ ≈ δ ∞ Z exp( − N ′ w )4 w (log( w )) / du When the diffusion coefficient is D , the formula changes to¯ τ ≈ δ D s log (cid:18) N a π / δ (cid:19) . (116)The next term in the expansion can be obtained by accounting for the logarithmic singu-larity in the expansion of Green’s function. When there are p windows, whose distancesfrom the initial position of the Brownian particle are d k = dist ( P , P k ), formula (116)changes to ¯ τ ≈ δ m D s log (cid:18) N a π / δ (cid:19) , (117)where δ m = min( d , ..d p ). The asymptotic formula (116) is compared with results ofBrownian simulations and shows very good agreement (Fig. 4). When absorbing windowsare ellipses, the Green’s function approach, based on Narrow Escape methodology, canbe applied as well [17]. 23 P a ×P end B A ¯ (cid:1) ( ) ( a . u . ) ¯ (cid:0) (1) = (cid:2) log N − Figure 4:
Extreme statistics of the narrow escape time through a small window in three dimen-sions. A. The geometry of the NEP for the fastest particle. In the simulation, the sphere hasa radius 5 µm , the absorbing window ∂ S a , has a radius ε = 0 . µm and the diffusion coefficientis D = 0 . µm s − . The trajectory starts at point P (cross), and ends at point P end . B. Plotof the MFPT of the fastest particle versus the number of particle N . We simulated 2000 runs.The asymptotic solution (red curve) is A/ log( N + B ). .2 Asymptotics in dimension 2 We consider the diffusion of N Brownian i.i.d. particles in a two-dimensional domain Ωwith a small absorbing arc ∂ Ω a of length 2 ε on the otherwise reflecting boundary ∂ Ω. Tocompute the pdf of the shortest arrival time to the arc, we follow the steps of the analysisin dimension 3, presented in the previous subsection 6.1.The Neumann-Green function (68) in two dimensions is the solution of the BVP − ∆ x ˆ G ( x , q | y ) + q ˆ G ( x , q | y ) = δ ( x − y ) for x , y ∈ Ω , (118) ∂ ˆ G q ( x , q | y ) ∂n x =0 for x , y ∈ ∂ Ω , (119)is given for x , y ∈ ∂ Ω by [24, p.51]ˆ G ( x , q | y ) = 1 π K ( √ q | x − y | ) + R ( x , y ) , (120)where R ( x , y ) is its regular part. For a disk, the analytical expression is given by theseries R ( x , y ) = 1 π ∞ X σ n cos( n ( ψ − ψ )) K ′ n ( √ q ) I ′ n ( √ q ) I n ( r √ q ) I n ( r √ q ) , (121)where σ = 1 , σ n = 2 for n ≥ x = re iψ , y = r e iψ . The integral representation (71)of the solution isˆ p ( x , q | y ) = ˆ G ( x , q | y ) − Z ∂ Ω a ∂ ˆ p ( x , q | y ′ ) ∂n x ˆ G ( x , q | y ′ ) dS y ′ , (122)so choosing x ∈ ∂ Ω a ,0 = ˆ G ( x , q | y ) − Z ∂ Ω a ∂ ˆ p ( x , q | y ′ ) ∂n x ˆ G ( x , q | y ′ ) dS y ′ . (123)This Helmholtz equation has the constant solution [17] ∂ ˆ p ( x , q | y ′ ) ∂n x = C for all x = A ∈ ∂ Ω a . (124)To leading order, we getˆ G ( A , q | y ) = Cπ Z ∂ Ω a K ( √ q | A − y | ) ds y , (125)25here ds y is arclength element in ∂ Ω a . When | x − y | ≤ ε and √ qε ≪ q expansion of Green’s function, K ( √ q | x − y | ) = − log( √ q | x − y | ) + log 2 − γ + o (1) , (126)we obtain ˆ G ( A , q | y ) = Cπ Z ∂ Ω a [ − log( √ q | A − y | ) + log 2 − γ + o (1)] ds y ; (127)that is, ˆ G ( A , q | y ) = 2 Cπ ε Z [ − log( √ qr ) + log 2 − γ + o (1)] dr. (128)Therefore the leading order approximation of C is C = π ˆ G ( A , q | y )2 ε [ − log( √ qε ) + O ( ε )] . (129)Finally, (122) gives for | A − x | ≫ ε ˆ p ( x , q | y ) ≈ ˆ G ( x , q | y ) + π ˆ G ( A , q | y ) ˆ G ( A , q | x )log( √ qε ) + O ( ε ) . (130)The inversion formula [23, p.1028] for k > L − ( K ( k √ q ) = 12 t e − k t , (131)gives L − ˆ G ( x , q | y ) = 14 πt e − | x − y | t . (132)For an initial point far from the boundary layer near the window, the expansion [23, p.378] K ( z ) = r π z e − z (cid:18) O (cid:18) z (cid:19)(cid:19) for z ≫ , (133)gives in (120)ˆ G ( A , q | y ) ˆ G ( A , q | x ) = 12 π r qs s e −√ q ( s + s ) (cid:0) O ( q − / ) (cid:1) , (134)26here s = | A − y | and s = | A − x | , π ˆ G ( A , q | y ) ˆ G ( A , q | x ) − log( √ qε ) + O ( ε ) = 1 − √ qε ) r qs s e −√ q ( s + s ) (cid:0) O ( q − / ) (cid:1) , ≈ − ε ) + O (1) r qs s e −√ q ( s + s ) (cid:0) O ( q − / ) (cid:1) . The inversion formula L − (cid:18) √ q e − k √ q (cid:19) = 1 √ πt e − k t (135)gives π L − ( ˆ G ( A , q | y ) ˆ G ( A , q | x ) − log( √ qε ) + O ( ε ) ) = 1 − ε ) + O (1) 1 √ πts s e − ( s + s ) t . (136)Hence, we obtain the short-time asymptotics of the survival probability S ( t ) ≈ Z Ω p t ( x , y ) d x (137)= 14 πt Z Ω e − | x − y | t d x (138) − − ε ) √ s + O (1) 1 √ πt Z Ω s | A − x | e − ( | A − x | + s ) t d x (139)= R ( t ) + R ( t ) (140)where R ( t ) = 14 πt Z Ω e − | x − y | t d x ≈ − e − ( R a / √ t ) , (141)and R a is the radius of the maximal disk inscribed in Ω. The second term is R ( t ) = − − ε ) √ s + O (1) 1 √ πt Z Ω s | A − x | e − ( | A − x | + s ) t d x ≈ − − ε ) √ s + O (1) r πt R a Z e − ( r + s ) t √ r dr. (142)27he small t Laplace expansion and the two successive changes of variable u = s t r and v = u / give R a Z e − ( r + s ) t √ rdr ≈ (cid:18) ts (cid:19) / e − s t ∞ Z e − v / dv. (143)Thus, with I = ∞ Z e − v / dv = 3 √ π , (144)we obtain R ( t ) − − ε ) + O (1) √ πts e − s t . (145)We conclude therefore that the survival probability (137) is approximately S ( t ) ≈ −
12 log( ε ) √ πts e − s t , (146)where the contribution of (141) is negligible. Thus the MFPT of the fastest particle isgiven by¯ τ = ∞ Z [Pr { t > t } ] N dt ≈ ∞ Z exp N log −
12 log( ε ) √ πts e − s t dt ≈ ∞ Z exp − N
12 log( ε ) √ πts e − s t dt. The computation of the last integral follows the steps described in subsection 5. Thechange of variable w = te − s t leads to the asymptotic formula with diffusion coefficient D ¯ τ ≈ s D log π √ N (cid:0) ε (cid:1) ! . where s = | x − A | and x is the position of injection A the center of the absorbingwindow. This formula is compared with Brownian simulations in Fig. 5. However, thedependence on the window size is log( ε ). 28
000 2000 3000 4000 5000 6000N2022242628 (cid:3) ( ) ( a . u . ) (cid:4) (1) = (cid:5) A B ∂ a × P Figure 5:
Escape through a narrow opening in a planar disk. A. The geometry of the NEP forthe fastest particle. B. Plot of the MFPT of the fastest particle versus the number of particles N . The asymptotic solution (red curve) is of the form α log( N )+ β . The geometry of a dendritic spine (see Fig.6) is composed of a head with a small holeopening, connected to a cylindrical neck. Initially, all Brownian particles, which representcalcium ions that are uniformly distributed in the spine head at the time of their release.This geometry implies that the mean time τ to reach the base of the neck is the meantime τ to reach the small window for the first time plus the mean time τ spent in thespine neck, with no possible returns: we assume here that when a particle enters the neckcylinder, it cannot return to the head.We compute now the distribution of the arrival time for the fastest and second fastestBrownian particle in a dendritic spine geometry. The pdf of a particle arriving at thebase of the dendrite within the time τ = τ + τ is computed as follows, when the escapetime form the head is Poissonian:Pr { τ = τ + τ = t } = t Z Pr { τ = t − s | τ = s } Pr { τ = s } ds. (147)The the Markovian property implies thatPr { τ + τ = t } = t Z Pr { τ = t − s } Pr { τ = s } ds. (148)Using the narrow escape theory [17], the distribution of arrival time of a Brownian particleat the entrance of the dendritic neck is Poissonian,Pr { τ = s } = γe − γs , (149)29here γ − = | Ω | aD (cid:20) L ( ) + N ( )2 π a log a + o ( a log a ) (cid:21) , with | Ω | the volume of the spherical head, while a is the radius of the cylindrical neck [17]and L (0) and N (0) are the principal mean curvature. After the first particles reaches thecylinder (spine neck), we approximate its Brownian motion in the cylindrical domain byone-dimensional motion (1D). By substituting N = 1 for the first arriving particle, thisis given by (16): Pr { τ = t − s } = N R ∞ X n =0 ( − n λ n e − Dλ n ( t − s ) . (150)Hence, Pr { τ + τ = t } = γN R t Z e − γs ∞ X n =0 ( − n λ n e − Dλ n ( t − s ) ds (151)= γN R ∞ X n =0 ( − n " e − Dλ n t − e − γt γ − Dλ n . (152)This result represents the pdf of the arrival time of a Brownian particle at the baseof a spine, or a process with two time scales: one dictated by diffusion and the otherPoissonian.The result (151) for the arrival time is derived for the study of the statistics of a singleparticle. The expression for flux,Φ( t ) = Pr { τ + τ = t } , (153)gives f min ( t ) = Pr { τ = min( t , . . . , t N ) = t } = N − t Z Φ( s ) ds N − Φ( t ) (154)= N [ J ( t )] N − Φ( t ) , (155)and using (151), we get J ( t ) = 1 − t Z Φ( s ) ds = 1 − γ ∞ X n =0 ( − n ( γ − Dλ n ) " − e − Dλ n t Dλ n − − e − γt γ . (156)30n the Poissonian approximation, the pdf of the arrival time τ (2) of the second fastestparticle is given by Pr { τ (2) = t } = N t Z f min ( t − s ) f min ( s ) ds, (157)where f min ( s ) = J ( s ) N − Φ( s ). Thus,Pr { τ (2) = t } = t Z [ J ( t − s )] N Φ( t − s )[ J ( s )] N Φ( s ) ds. (158)Expressions (154) and (158) represent the distributions of arrival times of the first andsecond Brownian particles, initially injected in the spine head and escape at the end ofthe cylindrical neck . These expressions are computed under the assumption of no return:when a particle enters the neck cylinder, it cannot return to the head. Consider Brownian particles that escape a dendritic spine (Fig. 6) into a dendrite withany number of returns to the head after crossing into the neck. Recrossing is defined tothe stochastic separatrix [25], the position of which is not known exactly. This effect islikely to impact the first arrival time. The pdf of no return is given by (151). To computethe pdf of the shortest escape time τ a with returns, we use Bayes’ law for the escapedensity, conditioned on any number of returns, given byPr { τ a = t } = ∞ X k =0 Pr { τ a = t | k } Pr { k } , (159)where Pr { k } = k is the probability that the particle returns k times to the head. Theparticle hits the stochastic separatrix [21] and then returns to the head, before reachingthe dendrite. The probability of the escape, conditioned on k returns, Pr { τ a = t | k } , canbe computed from the successive arrivals times to the stochastic separatrix, τ , ..τ k , sothat Pr { τ a = t | k } = Pr { τ + .. + τ k = t } . (160)Assuming that the arrival time to the stochastic separatrix is Poissonian with rate λ S [17],we obtain that Pr { τ + .. + τ k = t } = λ S t Z ( λ S s ) n − ( n − f ( t − s ) ds, (161)31here f ( t ) is the pdf of no return (151). ThereforePr { τ a = t } = 12 f ( t ) + ∞ X n =1 t Z λ S ( λ S s ) n − ( n − f ( t − s ) ds k . (162)Finally, Pr { τ a = t } = 12 f ( t ) + t Z exp( − λ S s/ f ( t − s ) ds, (163)Expression (151) with λ S = γ gives the final expression for the pdf of the escape time f return ( t ) = Pr { τ a = t } (164)= 12 f ( t ) + γN R ∞ X n =0 ( − n λ n γ γ − Dλ n ) " e − γt/ − e − γt γ/ − e − γ/ t − e − Dλ n t Dλ n − γ . The maximum of f return is achieved at the point t max ≈ γ log 2. The pdfs of the first andsecond arrivals are computed as f (1) min ( t ) = Pr { τ = min ( t , . . . , t N ) = t } = N − t Z f return ( s ) ds N − f return ( t ) , and as in (157), f (2) min ( t ) = Pr { τ (2) = t } = N t Z f (1) min ( t − s ) f (1) min ( s ) ds, (165)The pdfs of the fastest and second fastest arrival times are computed from (164), as inthe previous sections. Fig. 6 shows the pdf of the arrival time τ (2) . Note that a correctionis needed in (158), because the second particle is not necessarily located inside the headwhen the first one arrives at the base. Finally, the arrival time formula 165 can be furthercorrect by adding the distribution of the particle inside the head. The correction is similarto the one we obtain in dimension one (formula 34 and 35). It can be written here usingthe survival probability 156 as f (2) further ( t ) = Pr { τ (2) = t } = N t Z f (1) min ( t − s ) f (1) min ( s ) J ( t ) N − ds. (166)32 * ini al posi on of all par cles B C P r { τ ( ) = t } Time ( a.u.)
UncorrectedCorrected for returnfurther corrected D E = + * posi on of the second par cle when the (cid:6) rstone arrives * * * * Small targets (2) = 3.41ms τ (1) N=500
Time ( a.u. ) P r { τ ( ) = t } UncorrectedCorrected for return F further corrected (2) = 2.08ms Time ( a.u.) (cid:7) r { τ ( ) = t } (cid:8) r { τ ( ) = t } G H
Time ( a.u. ) τ (1) when analytical solution analytical solution τ (1) N=1000 τ (1) when τ (2) N=500 τ ( ) when τ (2) N=1000 τ ( ) when Figure 6:
A-B-C-D:
The geometry of the spine is a spherical head and cylindrical neck.Particles are released at the center of the head and must first reach the top of the neck and thendiffuse through the neck to reach the base. Time taken for each process is represented with thenotation τ and τ , making the total time to be τ = τ + τ . C shows a trajectory that canreturn to the head. Note that in D, when the first particle arrives, the second one, which couldhave return to the head, is now located inside the neck. E-F.
Plot of Pr { τ (1) = t } for values N = 500 , G-H.
Plot of Pr { τ (2) = t } for values N = 500 , D = 600 µm /s Conclusion and applications of extreme statisticsto fast time scale activation in cell biology
We derived here new asymptotics for the expected arrival time of fastest Brownian parti-cles in several geometries: half a line, a segment, a bounded domain in dimension two andthree that contains a small window, and spine-shaped geometry (a ball connected to thincylinder). We found that the geometry is involved and explored by the fastest particleand the pdf is defined by the shortest ray (with reflections if the is an obstacle [20]) fromthe source to the target, in contrast with the narrow escape problem [17], where the maingeometrical feature is the size of the window and the surface or volume of the domain. Wederived here new laws for the first arrival time of Brownian particles to a target, whichcan be summarized as ¯ τ d = δ min D ln( N √ π ) , in dim 1 (167)¯ τ d ≈ s D log π √ N (cid:0) ε (cid:1) ! , in dim 2 (168)¯ τ d ≈ δ m D s log (cid:18) N a π / δ (cid:19) , in dim 3 , (169)where δ min is the shortest ray from the source to the window, D is the diffusion coefficientand N is the number of particles, s = | x − A | and x is the position of injection and thecenter of the window is A . Formula (167) is very different from the classical NET, whichinvolves volume or surface area and mean curvature (in dimension 3). When the windowis located at the end of a cusp, the asymptotics for the fastest particle are yet unresolved.We further found that the rate of arrival cannot be approximated as Poissonian, be-cause the fastest particle can arrive at a time scale that can fall into the short-timeasymptotic. We further studied here the arrival of a second particle. The mean arrivaltime of the second can be influence by the distribution of all particles, especially on asegment, because the distribution of particles at the time the first particle’s arrival is notDirac’s delta function (see section 3). However, the case of a spine geometry is interesting,because the first particle may have already arrived, but all other particles can still be inthe head (due to the narrow opening at the neck-head connection). We further computedthe pdf of the arrival time when particle can return to the head after sojourn in the neck(see section 7.2).The present asymptotics have several important applications: activation of molecu-lar processes are often triggered by the arrival of the first particles (ions or molecules)to target-binding sites. The simplest model of the motion of calcium ions in cell biol-ogy, such as neurons or a dendritic spine (neglecting electrostatic interactions) is that of34ndependent Brownian particles in a bounded domain. The first two calcium ions thatarrive at channels (such as TRP) can trigger the first step of biochemical amplificationleading to the photoresponse in fly photoreceptor. Another example is the activation of aRyanodine receptor (RyaR), mediated by the arrival of two calcium ions to the receptorbinding sites, which form small targets. Ryanodin receptors are located at the base of thedendritic spine. Computing the distribution of arrival times of Brownian particles at thebase, when they are released at the center of the spine head, is a model of calcium releaseduring synaptic activation. Computing the distribution of arrival time reveals that thefastest ions can generate a fast calcium response following synaptic activity. Thus thefastest two calcium ions can cross a sub-cellular structure, thus setting the time scale ofactivation, which can be much shorter than the time defined by the classical forward rate,usually computed as the steady-state Brownian flux into the target, or by the narrowescape time [17]. Acknowledgements
This research was supported by the Fondation pour la Recherche M´edicale - ´Equipes FRM2016 grant DEQ20160334882.