Trapping of a run-and-tumble particle in an inhomogeneous domain: the weak noise limit
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b Trapping of a run-and-tumble particle in an inhomogeneous domain: the weak noiselimit
Paul C. Bressloff
Department of Mathematics, University of Utah, Salt Lake City, UT 84112 USA
A one-dimensional run-and-tumble particle (RTP) switches randomly between a left and rightmoving state of constant speed v . This type of motion arises in a wide range of applications in cellbiology, including the unbiased growth and shrinkage of microtubules or cytonemes, the bidirectionalmotion of molecular motors, and the “run-and-tumble” motion of bacteria such as E. coli . RTPsare also of more general interest within the non-equilibrium statistical physics community, both atthe single particle level and at the interacting population level, where it provides a simple exampleof active matter. In this paper we use asymptotic methods to calculate the mean first passagetime (MFPT) for a one-dimensional RTP to escape an effective trapping potential generated byspace-dependent switching rates. Such methods are part of a more general framework for studyingmetastability in so-called piecewise deterministic Markov processes (PDMPs), which include theRTP as a special case.
I. INTRODUCTION
Velocity jump processes, whereby a particle randomlyswitches between different velocity states, are finding agrowing number of applications in cell biology. The par-ticle could represent a bacterial cell such as
E. coli un-dergoing chemotaxis [1–3], a motor-cargo complex walk-ing along a cytoskeletal filament [4–9], the tip of a mi-crotubule undergoing alternating periods of growth andshrinkage (catastrophes) [10], or the tip of a cytonemefilament searching for a target cell during morphogenesis[11]. One of the simplest examples of a velocity jump pro-cess, is the so called one-dimensional run-and-tumble par-ticle (RTP), which switches between two velocity states ± v . (Within the context of E. coli , a run refers to a periodof almost constant ballistic motion, whereas tumbling isthe disordered local motion that selects a new randomdirection for the next run.) The run-and-tumble modelhas also attracted considerable recent attention withinthe non-equilibrium statistical physics community, bothat the single particle level and at the interacting popu-lation level, where it provides a simple example of activematter [12–14]. Studies at the single particle level includeproperties of the position density of a free RTP in oneand higher dimensions [15–18], first-passage time (FPT)properties [19–25], RTPs under stochastic resetting [26–28], and non-Boltzmann stationary states for an RTP ina confining potential [29–33].In the case of bacterial run-and-tumble, a chemotac-tic concentration gradient can bias the tumbling rate sothat the bacterium executes motion towards a source ofchemoattractant or away from a source of chemorepel-lant. In order to model one-dimensional chemotaxis us-ing the simplified RTP model, it is necessary to introducesome bias in the stochastic switching (tumbling) betweenthe velocity states ± v that depends on the extracellularconcentration gradient c [34, 35]. An alternative RTPmodeling paradigm is to assume that switching favors thenegative velocity state for x → ∞ and the positive veloc-ity state for x → −∞ . The spatially-dependent switching thus acts as an effective confining potential, which canlead to a non-Boltzmann-like stationary probability dis-tribution [36].The role of spatially dependent switching rates has alsobeen explored in a variety of intracellular transport mod-els, including both diffusive transport [37, 38] and activetransport. An example of the latter arises in the so-calledtug-of-war model of motor-driven bidirectional transportalong microtubules [39, 40]. Microtubules are polarizedpolymeric filaments with biophysically distinct ( + ) and (−) ends, and this polarity determines the preferred di-rection in which an individual molecular motor moves.In particular, kinesin motors move towards the (+) endwhereas dynein motors move towards the (−) end. Ifboth kinesin and dynein motors are attached to a vesic-ular cargo, then the velocity state will be determined byhow many of the kinesin and dynein motors are bound tothe microtubule at any one time. In addition, the switch-ing between different velocity states will depend on therates of binding and unbinding of individual motors tothe filament track. One mechanism for generating space-dependent transition rates involves microtubule associ-ated proteins (MAPs). These molecules bind to micro-tubules and effectively modify the free energy landscapeof motor-microtubule interactions. For example, tau is aMAP found in the axon of neurons and is known to bea key player in Alzheimer’s disease. Experiments haveshown that tau significantly alters the dynamics of ki-nesin; specifically, by reducing the rate at which kinesinbinds to the microtubule [41]. This can be interpretedas an effective space-dependent increase in the rate ofswitching to negative velocity states.The effect of local tau signaling on a tug-of-war modelhas been explored in terms of a multi-state velocity jumpprocess with space-dependent switching rates [7]. Anal-ogous to the more recent study of an RTP [36], a lo-cal increase in the tau concentration acts as an effec-tive confining potential for the motor complex. This canbe understood heuristically as follows. When a kinesindriven cargo encounters the MAP-coated trapping regionthe motors unbind at their usual rate and can’t rebind.Once the dynein motors are strong enough to pull the re-maining kinesin motors off the microtubule, the motor-complex quickly transitions to (−) end directed trans-port. After the dynein-driven cargo leaves the MAP-coated region, kinesin motors can then re-establish (+) end directed transport until the motor-complex returnsto the MAP-coated region. This back-and-forth motionrepeats until eventually the motor-complex is able tomove forward past the MAP-coated region. Interestingly,particle tracking experiments have observed oscillatorybehavior of motor-driven mRNA particles around synap-tic targets in the dendrites of neurons [42, 43]. This hasled to the hypothesis that local tau signaling enhancesthe probability of a motor complex delivering its vesic-ular cargo to a target [7]. The amount of time that themotor complex spends within the target domain can thenbe formulated as a mean FPT (MFPT) problem [8].One of the assumptions in Ref. [8] is that the switchingrates are fast relative to other dynamical processes (weaknoise assumption). This means that the escape from theeffective confining potential involves rare events that can-not be accurately captured using a diffusion approxima-tion of the velocity jump process. Instead, a combinationof the Wentzel-Kramers-Brillouin (WKB) method andmatched asymptotics are used to calculate the MFPT.Such methods have also been applied to a more generalclass of stochastic processes known as piecewise deter-ministic Markov processes (PDMPs), with particular ap-plications to stochastic ion channels [44–47], gene net-works [48–50] and stochastic neural networks [51–53]. APDMP involves the coupling between a discrete Markovchain N ( t ) ∈ { , , . . . , M } and a continuous process x ( t ) ∈ R d that evolves deterministically between jumpsin the discrete random variables [54]. That is, ˙ x = F n ( x ) when N ( t ) = n , where { F n ( x ) , n = , , . . . M − } is a setof vector fields. A velocity jump process is a special classof PDMP for which F n ( x ) = v n , where v n is the n -th ve-locity state, and a one-dimensional RTP corresponds tothe case M = v = v, v = − v .As far as we are aware, the connection between thestatistical physics of RTPs and the more general theoryof velocity jump processes and PDMPs has not been ex-plored in any detail. In this paper, we show how meth-ods developed to analyze metastability in PDMPs canbe used to study a one-dimensional RTP with an effec-tive trapping potential due to space-dependent switchingrates. In Sect. II we introduce the basic model and dis-cuss various choices for the transition rates. One of thesimplifying features of the model compared to more gen-eral PDMPs is that an exact solution for the stationarydistribution of the RTP position can be derived withoutrecourse to some approximation scheme such as WKB.The main part of the paper is developed in Sect. III,where we use the asymptotic analysis developed in Ref.[8] to calculate the MFPT for the RTP to escape the effec-tive trapping potential generated by the space-dependentswitching rates. II. RUN-AND-TUMBLE PARTICLE WITHSPACE-DEPENDENT SWITCHING
Consider an RTP that randomly switches between twoconstant velocity states labeled by n = , v = v and v = − v for some v >
0. The position X ( t ) of theparticle at time t evolves according to the velocity jumpprocess dXdt = v [ − n ( t )] , (2.1)where n ( t ) = ,
1. Furthermore, suppose that the particlereverses direction according to a two-state Markov chainwith space-dependent transition rates0 β ( x ) ÐÐ⇀↽ÐÐ α ( x ) . (2.2)Let p n ( x, t ) be the probability density of the RTP at po-sition x ∈ R at time t > n = ) and to the left ( n = ∂p ∂t = − v ∂p ∂x − β ( x ) p + α ( x ) p , (2.3a) ∂p ∂t = v ∂p ∂x + β ( x ) p − α ( x ) p . (2.3b)This is supplemented by the initial conditions x ( ) = x and n ( ) = n with probability ρ ,n such that ρ , + ρ , =
1. In matrix form, we can write ∂p n ∂t = − v n ∂p n ∂x + ∑ n = , Q nm ( x ) p m , (2.4)with v n = v ( − n ) and Q = ( − β ( x ) α ( x ) β ( x ) − α ( x ) ) . (2.5)The matrix version is easily generalizable to more thantwo velocity states.Let L be some characteristic distance, which could berelated to the space constant of a chemical concentra-tion gradient in chemotaxis, or the size of a target inmotor-driven cargo transport. This introduces a naturaltime scale T = L / v . Suppose that the transition rates α ( x ) , β ( x ) ≫ / T for all x ∈ R , so that we can we cantake α, β = O ( / ǫ ) on relevant length and time scales.Rescaling the transition rates in Eq. (2.4) thus gives ∂p n ∂t = − v n ∂p n ∂x + ǫ ∑ m = , Q nm ( x ) p m . (2.6)For a given x , define the average velocity V ( x ) = vρ ( x ) − vρ ( x ) , (2.7)where ρ ( x ) = α ( x ) α ( x ) + β ( x ) , ρ ( x ) = − ρ ( x ) (2.8)is the stationary probability distribution of the two-state Markov chain with generator Q ( x ) , that is, ∑ m = , Q nm ( x ) ρ m ( x ) =
0. Intuitively speaking, one ex-pects Eq. (2.1) to reduce to the deterministic dynamicalsystem dx ( t ) dt = V ( x ( t )) , x ( ) = x (2.9)in the fast switching or adiabatic limit ε →
0. That is,for sufficiently small ε , the Markov chain undergoes manyjumps over a small time interval ∆ t during which ∆ x ≈ m isapproximately ρ m ( x ) . This can be made precise in termsof a law of large numbers for velocity jump processes, aswell as more general PDMPs [55–58]. A. Concentration gradient
One possible source of space-dependent switching ortumbling rates is a chemical concentration gradient c ( x ) .For the sake of illustration, consider a simple phenomeno-logical model, in which the tumbling rates depend on thetime derivative of the concentration c ( t ) = c ( x ( t )) alongthe particle trajectory, where x ( t ) is the particle positionat time t [35]. Using the fact that ˙ c = ± vdc / dx , we take β ( x ) = k + k vc ′ ( x ) , α ( x ) = k − k vc ′ ( x ) . (2.10a)(For simplicity, switching depends on the instantaneousvalue of the concentration gradient rather than a timeaveraged change in concentration as is typical in bacte-rial chemotaxis [1].) The stationary probability densitiessatisfy the pair of equations v dp dx = − β ( x ) p ( x ) + α ( x ) p ( x ) , − v d p dx = β ( x ) p ( x ) − α ( x ) p ( x ) . Adding these two equations gives v dp dx − v dp dx = , which implies that the difference p ( x ) − p ( x ) = constant.Assuming that −∞ < x < ∞ , normalizability of the prob-ability densities requires this constant to be zero. Hence, p , ( x ) = p ( x )/ p ( x ) satisfying the single equation v dpdx = [ α ( x ) − β ( x )] p ( x ) = − k vc ′ ( x ) p ( x ) . This has the straightforward solution p ( x ) = N e − k c ( x ) , (2.11) where N is a normalization factor. If the signalingmolecules correspond to a chemoattractant then the rateof tumbling decreases in the direction for which ˙ c > k <
0, and maxima of the stationary solution(2.11) coincide with maxima of the concentration c ( x ) .Conversely, k > p ( x ) coincide with minima of the concentration. B. Localized trap
In this paper we are interested in a different form ofspace-dependent switching, namely one that traps theRTP within a local region which, without loss of gener-ality, we take to be in a neighborhood of the origin. Wewill consider two different examples of trapping mecha-nisms as illustrated in Fig. 1. The first mechanism favorsthe right-moving velocity when x < x >
0, Fig. 1(a). This can be imple- (a) β (x) α (x) (b) β (x) α (x) x = 0 s w i t c h i ng r a t e s w i t c h i ng r a t e κ κ κ κ FIG. 1. Space-dependent switching rates α ( x ) , β ( x ) for aone-dimensional RTP. (a) Switching rates (2.12). (b) Switch-ing rates (2.15). Thickness of the arrows indicates preferredvelocity direction. mented using the switching rates α ( x ) = κ + ( κ − κ )( + tanh ( x / γ )) , (2.12a) β ( x ) = κ + ( κ − κ )( + tanh ( x / γ )) . (2.12b)Clearly ( α ( x ) , β ( x )) → ( κ , κ ) as x → ∞ and ( α ( x ) , β ( x )) → ( κ , κ ) as x → −∞ . In addition, thesharpness of the transition is determined by γ such thatin the limit γ → β ( x ) = κ + ( κ − κ ) Θ ( x ) , (2.13a) α ( x ) = κ + ( κ − κ ) Θ ( x ) , (2.13b)where Θ ( x ) is the Heaviside function. The average ve-locity (2.7) is given by V ( x ) = κ − κ κ + κ tanh ( x / γ ) v. (2.14)Hence, if κ < κ then x = p ( x, t ) can be derived without restricting to the weaknoise regime [36].The second mechanism assumes that the left-movingstate is favored in a local region of the origin, whereasthe right-moving state is favored on either side of thisdomain, Fig. 1(b). The corresponding switching ratesare taken to be of the from α ( x ) = κ + ( κ − κ ) tanh ( x / γ )) , (2.15a) β ( x ) = κ + ( κ − κ )( − tanh ( x / γ )) . (2.15b)It can be seen that α ( ) = κ , β ( ) = κ , whereas α ( x ) → κ , β ( x ) → κ as ∣ x ∣ → ∞ . Moreover, the average velocity(2.7) is V ( x ) = v κ − κ κ + κ ( − ( x / γ )) . (2.16)Now there are two fixed points at x = ± ¯ x wheretanh ( ¯ x / γ ) = √ . (2.17)If κ < κ , then the fixed point − ¯ x is stable and the fixedpoint ¯ x is unstable. (The second mechanism is analo-gous to the trapping of a molecular motor complex by alocal region of enhanced tau concentration [7], which wasdescribed in the introduction.)The differences between the two cases becomes clearerby rewriting the deterministic Eq. (2.9) as the gradientsystem dxdt = − dU ( x ) dx , (2.18) -2 -1.5 -1 -0.5 0.5 1.5 210(a)(b) position xU γ = 0.1γ = 0.5γ = 1γ = 2 -0.8-0.6-0.4-0.200.20.40.60.8-3 -2 -1 1 320position xU γ = 0.5γ = 1γ = 2 FIG. 2. Deterministic potential U ( x ) corresponding to (a)the switching rates (2.12) and (b) the switching rates (2.15)for various gains γ . Other parameters are v = κ = . κ = with U ( x ) = vγ κ − κ κ + κ ln cosh ( x / γ ) , (2.19)for the switching rates (2.12) and U ( x ) = v κ − κ κ + κ [ − x + γ tanh ( x / γ )] (2.20)for the switching rates (2.15). As shown in Fig. 2(a),Eq. (2.19) corresponds to a global, symmetric potentialwell that has a unique minimum at x =
0, and sharpensas γ →
0. On the other hand, the potential of Eq. (2.20)is a cubic that is characterized by a potential well inthe domain ( −∞ , ¯ x ) with a minimum at − ¯ x and barrierheight U ( ¯ x ) − U ( − ¯ x ) , see Fig. 2(b). The deterministicpotential U ( x ) also determines the steady-state solutionof the CK Eq. (2.6). That is, the steady-state solutionis p ( x ) = p ( x ) = p ( x )/ p ( x ) = N e − Φ ( x )/ ǫ , (2.21)where N is a normalization factor andΦ ( x ) = − ∫ x α ( y ) − β ( y ) v dy = − κ + κ v ∫ x V ( y ) dy = κ + κ v U ( x ) . (2.22)One can identify Φ ( x ) as the so-called quasipotential. C. First passage time (FTP) problem
Given the potentials U ( x ) of Fig. 2, we would like todetermine the mean time for the RTP to escape a neigh-borhood of the origin in the weak noise regime. In thecase of the unimodal potential (2.19), we consider theMFPT for the particle to reach a location x ∗ ≫ x =
0. Onthe other hand, for the cubic potential (2.20), we con-sider the MFPT to reach the maximum ¯ x given that theparticle started at − ¯ x . These two escape problems areillustrated in Fig. 3. In both cases we use asymptoticmethods developed for general velocity jump processes[8]. In particular, we show that considerable simplifica-tion occurs in the case of an RTP where, for example, (a)(b) U x-x x0 x* _ _ FIG. 3. FPT problems for the (a) unimodal potential and (b)cubic potential. the exact stationary density (2.21) is known without anyrecourse to approximation schemes such as WKB. (Thetwo-state velocity jump process was not explicitly con-sidered in [8].) Moreover, we highlight a subtle featureof the asymptotic analysis of the cubic potential, arisingfrom the fact that the escape point is a maximum of thepotential, see also [49]. Note that one constraint on theuse of asymptotic methods is that the quasipotential istwice differentiable. Hence, space-dependent switchingrates such as Eq. (2.13) would need to be regularizedby replacing the Heaviside function with a sharp sigmoidfunction.
III. ASYMPTOTIC ANALYSIS OF THE MFPT
Consider the RTP with switching rates given by Eq.(2.12) or (2.15). In order to calculate the MFPT to es-cape a neighborhood of the origin, we supplement theCK equation (2.6) by the absorbing boundary condition p ( x ∗ , t ) = . (3.1)Note that the absorbing boundary condition is only im-posed on the component p associated with the negativevelocity, ensuring that once the RTP reaches x ∗ it cannever reenter the domain. In the case of the unimodalpotential we take x ∗ ≫
0, whereas for the cubic potentialwe set x ∗ = ¯ x , see Fig. 3. Let T denote the (stochastic)FPT for which the system first reaches x ∗ , given that itstarted at x = x . For the unimodal case x =
0, whereasfor the cubic case x = − ¯ x . The distribution of FPTs isrelated to the survival probability that the system hasn’tyet reached x ∗ , that is, P { t > T } = S ( t ) ≡ ∫ x ∗ −∞ ∑ n = , p n ( x, t ) dx. The FPT density is then f ( t ) = − dSdt = − ∫ x ∗ −∞ ∑ n = , ∂p n ( x, t ) ∂t dx. (3.2)Substituting for ∂p n / ∂t using the CK equation (2.6) andnoting that ∑ n Q nm ( x ) =
0, shows that f ( t ) = ∫ x ∗ −∞ ⎡⎢⎢⎢ ⎣ ∑ n = , v n ∂p n ( x, t ) ∂x ⎤⎥⎥⎥⎦ dx = ∑ n = , v n p n ( x ∗ , t ) = vp ( x ∗ , t ) ≡ J ( x ∗ , t ) , (3.3)where J ( x ∗ , t ) is the probability flux through the absorb-ing boundary. A. Quasistationary approximation
Consider an eigenfunction expansion of the time-dependent solution, p ( x, t ) = ∞ ∑ j = C j ( t ) φ j ( x ) , (3.4)where the eigenfunction φ j = ( φ j, , φ j, ) ⊺ satisfies thematrix operator equation L φ j ≡ diag ( v, − v ) ∂ φ j ( x ) ∂x − ǫ Q ( x ) φ j = λ j φ j , (3.5)together with the boundary condition φ j, ( x ∗ ) = . (3.6)Here diag ( a, b ) denotes the diagonal matrix with eigen-values a, b , Similarly, we define a set of eigenfunctions forthe adjoint operator L † given by L † ξ j ≡ − diag ( v, − v ) ∂ ξ j ( x ) ∂x − ǫ Q ⊺ ( x ) ξ j = λ j ξ j , (3.7)and the boundary condition ξ j, ( x ∗ ) = . (3.8)The two sets of eigenfunctions form a biorthogonal setaccording to the inner product rule ⟨ ξ j , φ k ⟩ ≡ ∫ x ∗ −∞ ∑ n = , ξ j,n ( x ) φ k,n ( x ) dx = δ j,k . (3.9)Substituting the eigenvalue expansion into Eq. (2.6)shows that the coefficients evolve according to the de-coupled equations dC j ( t ) dt = − λ j C j ( t ) . (3.10)If the absorbing boundary at x ∗ is replaced by a re-flecting boundary, then there is a single zero eigenvaluewhose corresponding eigenfunction is the stationary so-lution on the domain ( −∞ , x ∗ ) , and all other eigenvalueshave positive real parts. We can thus introduce the or-dering 0 = λ < Re [ λ ] ≤ Re [ λ ] ≤ . . . On the other hand, when there is an absorbing boundary,the stationary solution no longer exists due to an expo-nentially small probability flux leaving the system at x ∗ .(This assumes that 0 < ǫ ≪ λ is perturbed from zero,becoming an exponentially small positive principal eigen-value: 0 < λ ≪ Re [ λ ] ≤ Re [ λ ] ≤ . . . Hence, on intermediate time scales for which the proba-bility of escape is still negligible, contributions from alleigenvalues λ j , j ≥
1, have decayed to zero and we canmake the quasistationary approximation p n ( x, t ) ∼ C ( t ) φ ,n ( x ) , n = , . (3.11)In addition, φ ( x ) is almost identical to the stationarysolution (2.21) outside a neighborhood of x ∗ , so that φ , ( x ) = φ , ( x ) = φ ǫ ( x ) , φ ǫ ( x ) = e − Φ ( x )/ ǫ , (3.12) with Φ ( x ) given by Eq. (2.22). Clearly, the quasista-tionary solution breaks down around x ∗ since it does notsatisfy the absorbing boundary condition.It can be checked that under the quasistationary ap-proximation the solution C ( t ) = C ( ) e − λ t still holds.This follows from taking the inner product of Eq. (2.6)with the adjoint eigenvector ξ and substituting for p n ( x, t ) using the quasistationary approximation: ⟨ ξ , ∂ p ∂t ⟩ = − ⟨ ξ , L p ⟩ ⇒ ˙ C ⟨ ξ , φ ǫ ⟩ = − ⟨ L † ξ , φ ǫ ⟩ = − λ ⟨ ξ , φ ǫ ⟩ , (3.13)that is, ˙ C = − λ C . Hence, substituting the quasista-tionary solution into Eq. (3.2) gives f ( t ) ∼ C ( ) λ e − λ t ∫ x ∗ −∞ φ ǫ ( x ) dx. (3.14)The constant C ( ) can be determined from the initialcondition p n ( x, ) = δ ( x ) δ n,n , and the projection of theeigenfunction expansion onto the adjoint eigenfunction ξ : ⟨ ξ , C ( ) φ ǫ ( x )⟩ = ⟨ ξ , p ( x, )⟩ = ξ ,n ( ) . In the case of a reflecting boundary at x ∗ , the adjointeigenfunction ξ = ( , ) . This will still hold in the bulkof the domain for an absorbing boundary at x ∗ so thatwe can take C ( ) = [ ∫ x ∗ −∞ φ ǫ ( x ) dx ] − . (3.15)This establishes that under the quasistationary approxi-mation f ( t ) ∼ λ e − λ t , (3.16)and λ − can be identified as the MFPT to escape at x = x ∗ .In summary, the calculation of the MFPT reduces tothe problem of estimating the principal eigenvalue λ . Ifthe exact eigenfunctions φ and ξ were known then wecould use either of the inner product identities λ ⟨ ξ , φ ⟩ = ⟨ L † ξ , φ ⟩ , or λ ⟨ ξ , φ ⟩ = ⟨ ξ , L φ ⟩ . (3.17)On the other hand, simultaneously using the approxima-tions ξ = ( , ) and φ = φ ǫ yields λ =
0, reflecting thebreakdown of the quasistationary approximation at theboundary. Therefore, we only apply the quasistationaryapproximation to φ so that λ ∼ ⟨ L † ξ , φ ǫ ⟩⟨ ξ , φ ǫ ⟩ . (3.18)Substituting for L † , using integration by parts on thedomain ( −∞ , x ∗ ] , and using L φ ǫ =
0, shows that λ ∼ − vφ ǫ ( x ∗ )[ ξ , ( x ∗ ) − ξ , ( x ∗ )]⟨ ξ , φ ǫ ⟩ . (3.19)Following [8], the adjoint eigenfunction ξ ( x ) can be ap-proximated using singular perturbation methods. It isat this stage that escape from the unimodal and cubicpotentials have to be treated separately. B. Calculation of principal eigenvalue: unimodalpotential
In order to construct an approximate solution that alsosatisfies the absorbing boundary condition, we constructa boundary layer in a neighborhood of x ∗ by performingthe change of variables x = x ∗ − ǫz and setting A n ( z ) = ξ ,n ( x ∗ − ǫz ) . Eq. (3.7) for j = v n dA n ( z ) dz − ∑ m = , Q mn ( x ∗ ) A m ( z ) = , (3.20)together with the boundary condition A ( x ∗ ) = . (3.21)This inner solution has to be matched with the outersolution ξ = , which means thatlim z → ∞ A n ( z ) = , n = , . (3.22)Consider the eigenvalue equation ∑ m = , S m Q mn ( x ) v − n = µS n . (3.23)One solution is S = ( , ) and µ =
0, whereas theother is S ( x ) = ( β ( x ) , α ( x )) and µ = − Φ ′ ( x ) = ( α ( x ) − β ( x ))/ v . We now expand the solution A n ( z ) in terms ofthe pair of eigenfunctions at x = x ∗ : A n ( z ) = c + c S ,n ( x ∗ ) e − Φ ′ ( x ∗ ) z . (3.24)Since Φ ′ ( x ∗ ) > x ∗ > A n ( z ) → c as z → ∞ ,which implies c =
1. The constant c is then determinedfrom the boundary condition A ( ) = c = − α ( x ∗ ) . (3.25)It follows that ξ , ( x ∗ ) − ξ , ( x ∗ ) = α ( x ∗ ) − β ( x ∗ ) α ( x ∗ ) . (3.26)Substituting the expressions for φ ǫ ( x ∗ ) and ξ , ( x ∗ ) − ξ , ( x ∗ ) into Eq. (3.19) and simplifying the denominatorusing the outer solution ξ ,n ∼
1, we obtain the result λ ∼ N v β ( x ∗ ) − α ( x ∗ ) α ( x ∗ ) e − Φ ( x ∗ )/ ǫ , (3.27) x* l og E [ T ] γ = 0.1γ = 0.5γ = 1 FIG. 4. Plot of log E [ T ] for the unimodal potential as a func-tion of escape position x ∗ and various gains γ . The MFPT E [ T ] = λ − with λ given by Eq. (3.29). Other parametervalues are κ = . κ = v =
1, and ǫ = . where N = [ ∫ x ∗ −∞ exp ( − Φ ( x ) ǫ )] − . The latter can be approximated using Laplace’s methodto give
N ∼ √ Φ ′′ ( x ) πǫ exp ( Φ ( x ) ǫ ) . (3.28)Hence, we obtain the following expression for the inverseMFPT: λ ∼ v β ( x ∗ ) − α ( x ∗ ) α ( x ∗ ) √ Φ ′′ ( x ) πǫ e −( Φ ( x ∗ )− Φ ( x ))/ ǫ . (3.29)Setting x = E [ T ] = λ − as a function of the escapeposition x ∗ for various degrees of sharpness γ . The re-sults are shown in Fig. 4. C. Calculation of principal eigenvalue: cubicpotential
The boundary layer analysis of the unimodal potentialbreaks down in the case of the cubic potential due to thefact that Φ ′ ( x ∗ ) =
0; this is a consequence of the escapepoint x ∗ being a local maximum of the potential. Inparticular, the eigenfunction expansion (3.24) no longerholds since the zero eigenvalue is doubly degenerate at x = x ∗ . Hence, the solution needs to include a secularterm involving the generalized eigenvector ̂ S , ∑ n = , Q mn ( x ∗ ) ̂ S m ( x ∗ ) = − v n , (3.30) γ = 0.3 v = 0.5 v = 1 l og E [ T ] κ l og E [ T ] = 0.3 v = 0.5 v = 1 FIG. 5. Plot of log E [ T ] for the cubic potential as a functionof (a) the gain γ and (b) the rate κ . The MFPT is E [ T ] = λ − with λ given by Eq. (3.35). Baseline parameter values are γ = κ = . κ = v =
1, and ǫ = . which implies that ̂ S ( x ∗ ) − ̂ S ( x ∗ ) = vα ( x ∗ ) . (3.31)Note that the Fredholm alternative theorem ensuresthat ̂ S exists and is unique, since the stationary dis-tribution ρ m ( x ∗ ) is the right null vector of Q ( x ∗ ) and ∑ n = , ρ n ( x ∗ ) v n ≡ V ( x ∗ ) =
0; the latter reflects the factthat x ∗ = ¯ x is a fixed point of the deterministic equation(2.9). The solution for Q ( z ) is now A n ( z ) = c + c ( ̂ S n ( x ∗ ) − z ) . (3.32)The presence of the secular term means that the solutionis unbounded in the limit z → ∞ , which implies that theinner solution cannot be matched with the outer solu-tion. One way to remedy this situation is to introducean alternative scaling in the boundary layer of the form x = x ∗ + ǫ / z , as detailed in Ref. [48]. One can then elim-inate the secular term − c z and show that, see appendix A, c = − ¯ c √ π ∣ Φ ′′ ( x ∗ )∣ , c = √ ǫ ¯ c , (3.33)with ¯ c determined by imposing the boundary condition A ( ) =
0: ¯ c ∼ − √ ∣ Φ ′′ ( x ∗ )∣ π + O ( ǫ / ) , (3.34)Substituting the expressions for φ ǫ ( x ∗ ) and ξ , ( x ∗ ) − ξ , ( x ∗ ) into Eq. (3.19), simplifying the denominatorusing the outer solution ξ ,n ∼ λ ∼ π vα ( ¯ x ) √ Φ ′′ ( − ¯ x )∣ Φ ′′ ( ¯ x )∣ e − [ Φ ( ¯ x ) − Φ ( − ¯ x )/ ǫ , (3.35)We have also set x ∗ = ¯ x and x = − ¯ x with ¯ x determinedby Eq. (2.17). Example plots of the MFPT E [ T ] = λ − are shown in Fig. 5 for Φ given by Eqs. (2.20) and(2.22), and α ( x ) obtained from Eq. (2.15). The MFPTis a monotonically increasing function of the gain γ , sincethe barrier height increases with γ :Φ ( ¯ x ) − Φ ( − ¯ x ) = γ ( κ − κ ) v (√ − tanh − ( /√ )) . (3.36)Similarly, the MFPT is a decreasing function of the rate κ as the barrier height becomes smaller as κ approaches κ . The nonmonotonic behavior of E [ T ] for κ ≈ κ indicates a breakdown of the asymptotic analysis whenthe barrier height becomes too small. IV. DISCUSSION
In this paper we exploited the connection betweenRTPs and more general velocity jump processes in orderto calculate the MFPT for the RTP to escape from an ef-fective trapping potential in the weak noise limit. In par-ticular, following previous studies of motor-driven bidi-rectional transport, we showed how the inverse MFPTcan be identified with the principal eigenvalue λ of theCK evolution operator. We then calculated λ usingasymptotic analysis, in order to match the quasistation-ary solution in the bulk of the domain with an absorbingboundary at the escape point. We also highlighted sub-tle differences between the unimodal and cubic trappingpotentials.One issue that we did not address is to what extent onecan investigate the behavior of the RTP in the weak noiseregime using a quasi-steady-state (QSS) or adiabatic ap-proximation. It is well known that in the adiabatic limit,the CK equation of a velocity jump process or a moregeneral PDMP can be approximated by a Fokker-Planck(FP) equation for the total density p = p + p [3–7, 58–61]. The basic idea is to decompose the solution to theCK Eq. (2.6) according to p m ( x, t ) = p ( x, t ) ρ m ( x ) + ǫw n ( x, t ) , (4.1)where ∑ m = , w m ( x, t ) =
0. Using a Liapunov-Schmidtreduction one can derive the FP equation ∂p∂t = − ∂∂x ( V ( x ) p ) + ǫ ∂∂x ( D ( x ) ∂p∂x ) , (4.2)where we have dropped an O ( ǫ ) contribution to the driftterm, and D ( x ) = v α ( x ) β ( x ) α ( x ) + β ( x ) . (4.3)Under this approximation, the position of the RTPevolves according to the stochastic differential equation dX = V ( X ) dt + √ ǫD ( X ) dW ( t ) , (4.4)where W ( t ) is a Wiener process with ⟨ W ( t )⟩ = , ⟨ W ( t ) W ( t ′ )⟩ = min { t, t ′ } . (4.5)Given the specific form of the FP Eq. (4.2), the multi-plicative noise is defined according to the kinetic inter-pretation of stochastic calculus.Although the diffusion approximation is useful in cap-turing certain time-dependent aspects of the RTP, itbreaks down in the large time limit. In particular, ityields a poor estimate of the stationary density of theexact model (2.6). This point was originally highlightedwithin the context of molecular transport models [7, 8].The normalizability of the stationary density requires thecorresponding flux to be zero for all x ∈ R . In the case ofthe FP equation (4.2) this means J ( x ) = − V ( x ) p ( x ) + ∂ [ D ( x ) p ( x )] ∂x = , which yields the stationary density p ( x ) = N e − Ψ ( x )/ ǫ , (4.6)where Ψ ( x ) = ∫ x V ( y ) D ( y ) dy. (4.7)Clearly the quasipotential Ψ ( x ) differs from the exactquasipotential Φ ( x ) of Eq. (2.21), resulting in exponen-tially significant errors for small ǫ . Following [8], we canunderstand the source of this error by noting that thezero flux condition of the exact model (2.6) implies J ( x ) = ∑ n = , v n p n ( x ) = , that is, p ( x ) = p ( x ) = p ( x )/
2. The underlying as-sumption of the QSS reduction is that the solution isclose to the stationary distribution of the Markov chain,that is, p n ( x ) ∼ ρ n ( x ) . Therefore, in order to be consis-tent with the exact zero flux condition, we would require ∑ n = , v n ρ n ( x ) = V ( x ) = x . This contradicts thefact that V ( x ) only vanishes at x = x = ± ¯ x for the switching rates (2.15).The problems with the diffusion approximation for anRTP also carry over to the calculation of the MFPT.Although the diffusion approximation breaks down inthe long time limit, it can capture the behavior of a ve-locity jump process on shorter time-scales. For exam-ple, it would apply to FPT problems outside the weaknoise regime where rare events dominate. This has beenshown in a wide variety of models of motor-driven intra-cellular transport [9]. It is particularly useful when thenumber of velocity states are greater than two or trans-port occurs in more than one spatial dimension. Both ofthese latter features have been included in RTP models[18, 23, 28, 32]. A more challenging problem is extend-ing the asymptotic analysis of escape problems for RTPswith multiple internal states moving in two or more spa-tial dimensions. The first step would be to identify anappropriate mechanism for trapping. APPENDIX A: BOUNDARY LAYER ANALYSISFOR THE CUBIC POTENTIAL.
In this appendix we summarize the boundary layeranalysis of Ref. [49], which leads to the result (3.34).Again the analysis simplifies greatly by focusing on thetwo-state RTP model rather than developing the theoryfor a general PDMP, which introduces additional techni-calities. In order to deal with the blow up of the secularterm in Eq. (3.32), we introduce an additional transitionlayer between the bulk or outer solution and the bound-ary layer. The scaling of this transition layer is deter-mined by performing the change of variables x = x ∗ − ǫ θ y ,0 < θ <
1, and defining B n ( y ) = ξ ,n ( x ∗ − ǫ θ y ) . (A.1)Introduce the asymptotic expansion B n ( y ) ∼ B ( ) n ( y ) + ǫ s B ( ) n ( y ) + ǫ s B ( ) n ( y ) , s > . (A.2)Eq. (3.7) for j = ∑ m = , [ δ n,m ǫ − θ v n ddy − [ Q mn ( x ∗ ) − ǫ θ yQ ′ mn ( x ∗ ) + . . . ] ] × ( B ( ) m ( y ) + ǫ s B ( ) m ( y ) + ǫ s B ( ) m ( y )) = . (A.3)The O ( ) equation is ∑ m = , Q mn ( x ∗ ) B ( ) m ( y ) = , (A.4)which implies that B ( ) n ( y ) = a ( y ) (A.5)for some scalar function a ( y ) . The expansion (A.3) thenbecomes ǫ − θ v n a ′ ( y ) − ǫ s ∑ m = , Q mn ( x ∗ ) B ( ) m ( y ) + O ( ǫ ) + o ( ǫ s , ǫ − θ ) = . (A.6)0This suggests taking s = − θ , which yields the O ( ǫ − θ ) equation B ( ) m ( y ) = − a ′ ( y ) ̂ S m ( x ∗ ) , (A.7)where ̂ S ( x ∗ ) is the solution to Eq. (3.30). Combiningthe results so far, we have B n ( y ) ∼ a ( y ) − ǫ − θ a ′ ( y ) ̂ S n ( x ∗ ) . (A.8)The next step is to calculate a ( y ) by proceeding tohigher order. We find − ǫ ( − θ ) a ′′ ( y ) v n ̂ S n ( x ∗ ) − ǫ ( − θ ) ∑ m = , Q mn ( x ∗ ) B ( ) m ( x ∗ ) − ǫ ( − θ ) θ ya ′ ( y ) ∑ m = , Q ′ mn ( x ∗ ) ̂ S m ( x ∗ ) = . (A.9)Setting θ = / a ′′ ( y ) v n ̂ S n ( x ∗ ) + ya ′ ( y ) ∑ m = , Q ′ mn ( x ∗ ) ̂ S m ( x ∗ ) = − ∑ m = , Q mn ( x ∗ ) B ( ) m ( x ∗ ) . (A.10)Multiplying both sides by ρ n ( x ∗ ) , summing over n andapplying the Fredholm alternative theorem leads to thesolvability condition a ′′ ( y ) ∑ n = , ρ n ( x ∗ ) v n ̂ S n ( x ∗ ) + ya ′ ( y ) ∑ m,n = , ρ n ( x ∗ ) Q ′ mn ( x ∗ ) ̂ S m ( x ∗ ) = . (A.11)In addition, Q ( x ∗ ) ρ ( x ∗ ) = Q ( x ∗ ) ρ ′ ( x ∗ ) = − Q ′ ( x ∗ ) ρ ( x ∗ ) , and ∑ m,n = , ̂ S m ( x ∗ ) Q mn ( x ∗ ) ρ ′ n ( x ∗ ) = ∑ m,n = , ̂ S m ( x ∗ ) Q mn ( x ∗ ) ρ ′ n ( x ∗ ) = − ∑ n = , v n ρ ′ n ( x ∗ ) . (A.12)Therefore, Eq. (A.11) reduces to the form a ′′ ( y ) + ya ′ ( y ) ∑ n = , v n ρ ′ n ( x ∗ ) ∑ n = , ρ n ( x ∗ ) v n ̂ S n ( x ∗ ) = . (A.13)Noting that ∑ n = , v n ρ ′ n ( x ∗ ) = V ′ ( x ∗ ) and ∑ n = , ρ n ( x ∗ ) v n ̂ S n ( x ∗ ) = ρ ( x ∗ ) v [ ̂ S ( x ∗ ) − ̂ S ( x ∗ )] , it follows that the fraction on the left-hand side is equalto − Φ ′′ ( x ∗ ) so we have a ′′ ( y ) − ya ′ ( y ) Φ ′′ ( x ∗ ) = . (A.14)Exploiting the fact that Φ ′′ ( x ∗ ) <
0, the solution for a ′ ( y ) is a ′ ( y ) = ¯ c e Φ ′′ ( x ∗ ) y / and thus a ( y ) = ¯ c + ¯ c ∫ y e Φ ′′ ( x ∗ ) u / du, (A.15)where ¯ c , ¯ c are integration constants.Substituting the solution for a ( y ) into (A.8) and set-ting θ = / B n ( y ) ∼ ¯ c + ¯ c ∫ y e Φ ′′ ( x ∗ ) u / du − √ ǫ ¯ c e Φ ′′ ( x ∗ ) y / ̂ S n ( x ∗ ) , (A.16)which replaces Eq. (3.32). This solution is bounded as y → ∞ so it can be matched with the outer solution, thatis, lim y → ∞ B n ( y ) =
1. Hence,¯ c + ¯ c ∫ ∞ e Φ ′′ ( x ∗ ) u / du = ¯ c + ¯ c √ π ∣ Φ ′′ ( x ∗ )∣ = . (A.17)In addition, as y →
0, we have B n ( y ) ∼ ¯ c + ¯ c y − √ ǫ ¯ c ̂ S n ( x ∗ ) . (A.18)Matching with the boundary layer solution (3.32) thenimplies that c = ¯ c and c = − √ ǫ ¯ c . Finally, if we imposethe absorbing boundary condition B ( ) = x = x ∗ and express ¯ c in terms of ¯ c , then1 − ¯ c √ π ∣ Φ ′′ ( x ∗ )∣ − √ ǫ ¯ c ̂ S n ( x ∗ ) = . (A.19)On rearranging we recover Eq. (3.34). [1] H. C. Berg and E. M. Purcell, Physics of chemoreception.Biophys. J. E. Coli in Motion , New York, Springer(2004). [3] T. Hillen and H. Othmer, The diffusion limit of transportequations derived from velocity-jump processes. SIAM J.Appl. Math. traveling waves in linear reaction-hyperbolic equations.SIAM J. Appl. Math. ,217-246 (2005).[6] J. M. Newby and P. C. Bressloff, Quasi-steady state re-duction of molecular-based models of directed intermit-tent search. Bull. Math. Biol. , 1347-1350 (1993).[11] P. C. Bressloff and H. Kim, A search-and-capture modelof cytoneme-mediated morphogen gradient formation.Phys. Rev. E , 218103 (2008).[13] M. E. Cates and J. Tailleur, Motility-induced phase sepa-ration, Annu. Rev. Condens. Matter Phys. , 219 (2015).[14] C. Bechinger, R. Di Leonardo, H. Lowen, C. Reichhardt,G. Volpe, and G. Volpe, Active particles in complexand crowded environments, Rev. Mod. Phys. , 84 (2012).[16] G. Gradenigo and S. N. Majumdar, A first-order dynami-cal transition in the displacement distribution of a drivenrun-and-tumble particle, J. Stat. Mech. 053206 (2019).[17] P. Singh and A. Kundu, Generalised “Arcsine” laws forrun-and-tumble particle in one dimension J. Stat.Mech.083205 (2019) .[18] I. Santra, U. Basu and S. Sabhapandit, Run-and-tumbleparticles in two dimensions: Marginal position distribu-tions. Phys. Rev. E , 062120 (2020)[19] L. Angelani, R. Di Lionardo, and M. Paoluzzi, First-passage time of run-and-tumble particles, Eur. Phys. J.E , 59 (2014).[20] L. Angelani, Run-and-tumble particles, telegrapher’sequation and absorption problems with partially reflect-ing boundaries, J. Phys. A: Math. Theor. , 495003(2015).[21] K. Malakar, V. Jemseena, A. Kundu, K. Vijay Kumar, S.Sabhapandit, S. N. Majumdar, S. Redner, and A. Dhar,Steady state, relaxation and first-passage properties of arun-and-tumble particle in one-dimension, J. Stat. Mech.043215 (2018).[22] A. Scacchi and A. Sharma, Mean first passage time ofactive Brownian particle in one dimension, MolecularPhysics , 032604 (2018).[25] P. Le Doussal, S. N. Majumdar, and G. Schehr, Non- crossing run-and-tumble particles on a line, Phys. Rev.E , 012113 (2019).[26] S. N. Majumdar and M. Evans, Run and tumble particleunder resetting: a renewal approach, Journal of PhysicsA: Mathematical and Theoretical
47 (2018).[27] P. C. Bressloff, Occupation time of a run-and-tumble par-ticle with resetting. Phys. Rev. E ,032132 (2019).[30] F. J. Sevilla, A. V. Arzola, and E. P. Cital, Stationarysuperstatistics distributions of trapped run-and-tumbleparticles, Phys. Rev. E , 012145 (2019).[31] Y. Ben Dor, E. Woillez, Y. Kafri, M. Kardar, and A. P.Solon, Ramifications of disorder on active particles in onedimension, Phys. Rev. E Biophysics.
Princeton University Press,Princeton (2012).[36] P. Singh, S. Sabhapandit and S. N. Kundu, Run-and-Tumble particle in inhomogeneous media in one dimen-sion (2021)[37] P. C. Bressloff and S. D. Lawley. Temporal disorder asa mechanism for spatially heterogeneous diffusion. Phys.Rev. E R1-11 (2004)[40] M. J. I. Muller, S. Klumpp and R. Lipowsky, Tug-of-war as a cooperative mechanism for bidirectional cargotransport by molecular motors. Proc. Natl. Acad. Sci.USA ,4609-4614 (2008)[41] M. Vershinin, B. C. Carter, D. S. Razafsky, S. J. Kingand S. P. Gross, Multiple-motor based transport and itsregulation by Tau. Proc. Natl. Acad. Sci. U.S.A. α
3’ un-translated regions-directed mRNA translocation in livingneurons: Visualization by GFP linkage J. Neurosci. channel noise. Phys. Rev. Lett. Memoirs of the AMS issue 944 (2009)[56] A. Faggionato, D. Gabrielli and M. R. Crivellari, Non-equilibrium thermodynamics of piecewise deterministicMarkov Processes. J Stat Phys86