A Explicit formula for principal eigenvalue

Interlacing Relaxation and First-Passage Phenomena in Reversible Discrete and Continuous Space Markovian Dynamics


We uncover a duality between relaxation and first passage processes in ergodic reversible Markovian dynamics in both discrete and continuous state-space. The duality exists in the form of a spectral interlacing – the respective time scales of relaxation and first passage are shown to interlace. Our canonical theory allows for the first time to determine the full first passage time distribution analytically from the simpler relaxation eigenspectrum. The duality is derived and proven rigorously for both discrete state Markov processes in arbitrary dimension and effectively one-dimensional diffusion processes, whereas we also discuss extensions to more complex scenarios. We apply our theory to a simple discrete-state protein folding model and to the Ornstein-Uhlenbeck process, for which we obtain the exact first passage time distribution analytically in terms of a Newton series of determinants of almost triangular matrices. In addition to these stand-alone results the work contains all proofs and extended technical explanations of the results reported in the accompanying letter (Hartich and Godec, submitted).

March 6, 2018

I Introduction

In his seminal work Kramers (1940) Kramers analyzed the kinetics of chemical reactions in terms of diffusive barrier crossing, assuming that the kinetic rate of a chemical reaction corresponds to the inverse of the mean first crossing time. Ever since, first passage theory is at the heart of theoretical descriptions of kinetics of chemical reactions Szabo et al. (1980); Ben-Naim et al. (1993); Oshanin et al. (1995); Mejía-Monasterio et al. (2011); Guérin et al. (2016); Li et al. (2017); see e.g. Hänggi et al. (1990); Redner (2001); Metzler et al. (2014); Bénichou and Voituriez (2014) for comprehensive reviews.

In a broader context, first passage concepts were invoked in studies of kinetics in complex media, such as reactions in fractal-like Kopelman (1988); ben Avraham and Havlin (2000) and planar domains Rupprecht et al. (2015); Grebenkov (2016), in inhomogeneous cellular environments Bressloff and Newby (2013); Godec and Metzler (2015); Vaccario et al. (2015); Godec and Metzler (2016a), in the study of neural networks Vilela and Lindner (2009); Braun et al. (2015), as well as in diverse narrow escape problems Singer et al. (2006); Schuss et al. (2007); Reingruber and Holcman (2009); Pillay et al. (2010); Isaacson et al. (2016); Grebenkov and Oshanin (2017) and so-called intermittent search strategies involving searching agents with internal dynamics Palyulin et al. (2014); Godec and Metzler (2017) (see also Bénichou et al. (2011) for a review).

First passage times play an important role in quantifying persistence properties in non-equilibrium interacting many-body systems Majumdar and Bray (2001); Majumdar and Comtet (2002); Bray et al. (2013). More recent applications of first passage concepts also include stochastic thermodynamics Jarzynski (2011); Seifert (2012); Van den Broeck and Esposito (2015), in particular, fluctuation relations for stopping time statistics and stochastic entropy production in driven molecular systems Neri et al. (2017) and in stochastic resetting processes Fuchs et al. (2016); Roldán and Gupta (2017), as well as uncertainty relations for first passage time statistics of fluctuating currents Garrahan (2017); Gingrich and Horowitz (2017) (see also Barato and Seifert (2015)).

Moreover, our current understanding of the speed and precision of transcription regulation in biological cells, and in particular of the role of the so-called proximity effect in the co-regulation of genes, Kolesov et al. (2007); Fraser and Bickmore (2007) builds on first passage time ideas. The corresponding physical principles underlying these proximity effects were explained in Bénichou et al. (2010); Meyer et al. (2011); Godec and Metzler (2016b). Notably, these works revealed the inherent insufficiency of the mean first passage time and traditional rate-based concepts for a quantitative description of biophysical dynamics in the so-called few encounter limit Godec and Metzler (2016b). As a result, a quantitative understanding of phenomena such as gene regulation Kolesov et al. (2007); Fraser and Bickmore (2007); Bénichou et al. (2010); Meyer et al. (2011); Bialek (2012); Pulkkinen and Metzler (2013); Godec and Metzler (2016b) and the misfolding-triggered pathological aggregation of proteins Dobson (2003); Chiti and Dobson (2006); Zheng et al. (2013); Yu et al. (2015); Dee and Woodside (2016); Nguyen et al. (2017), which are discussed in more detail in the accompanying letter 1, requires the consideration of the full statistics of first passage time.

Existing studies of the full first passage statistics in physical systems typically focus on systems with continuous state-space dynamics, whereas much less emphasis is put on discrete-space dynamics Schnoerr et al. (2017a). Recent investigations of such discrete-state dynamics include, for example, simple models of enzyme kinetics Munsky et al. (2009); Bel et al. (2010); Grima and Leier (2017) and novel numerical approximation schemes for studying first-passage statistics based on Bayesian inference Schnoerr et al. (2017b) (see also Weber and Frey (2017) for a recent review).

Complementary to first passage processes are relaxation dynamics, which by contrast do not terminate upon reaching a given threshold for the first time. Relaxation phenomena in reversible diffusive dynamical systems are nowadays well understood in terms of the eigenmodes and eigenvalues of the underlying Fokker-Planck operators, which provide a generic and very intuitive understanding of the dynamics of complex stochastic systems Biroli and Kurchan (2001); Tănase-Nicola and Kurchan (2003, 2004). Conversely, despite for allowing an analogous spectral representation, a similar intuitive understanding of the full first passage statistics and its physical implications remains elusive. Notwithstanding, an important approximate link between the mean first passage time for escaping the deepest potential basin and the corresponding slowest relaxation mode in the potential was established in the seminal works of Matkovsky and Schuss Schuss and Matkowsky (1979); Matkowsky and Schuss (1981), which has ever since been used routinely in explaining relaxation phenomena in condensed matter systems. Nevertheless, a deeper and more generic connection between the two paradigms to date was not yet established.

Here, we present the complete duality between relaxation and first passage phenomena, which holds for all ergodic Markov processes obeying detailed balance in both, continuous and discrete state-space, in which the absorbing target is effectively one-dimensional. The duality emerges in the form of a spectral interlacing, which we prove rigorously by combining spectral-theoretic, matrix-algebraic and Greens function-theoretic concepts. On the one hand the duality allows for an intuitive generic understanding of first passage phenomena in terms of relaxation eigenmodes. On the other hand, it enables us to determine the full first passage time statistics exactly from the corresponding relaxation eigensystem. The formalism is exact and holds for all reversible Markovian systems governed by a master equation in arbitrary dimensions or by a Fokker-Planck equation, and therefore unifies the theoretical treatment of discrete and continuous space phenomena.

To illustrate the predictive power of the formalism in practice, we here predominantly focus on systems with discrete state-space dynamics, whereas continuous space dynamics are treated in more detail in the accompanying letter 2. In particular, we here apply our theory to a simple discrete-state protein folding model and to diffusion in a harmonic potential, also know as the Ornstein-Uhlenbeck process. Notably, we obtain, to the best of our knowledge, for the first time an exact analytical solution for the full first passage time distribution of the Ornstein-Uhlenbeck process.

The paper is organized as follows. In Sec. II we present a canonical formulation the first passage problem applicable to both discrete states-pace and continuous Fokker-Planck dynamics. Sections III and IV provide a step-by-step explanation of how one can exactly determine the first-passage distribution from the corresponding relaxation process, and also contain rigorous proofs of the duality in discrete and continuous state-space dynamics, respectively. We apply the duality framework in Sec. V to determine the first passage statistics for a simple protein folding model and for the Ornstein-Uhlenbeck process. A concluding perspective is provided in Sec. VI.

Ii Fundamentals

ii.1 Relaxation and first passage

We assume that the probability density to find the system in state at time upon evolving from an initial state according to microscopically reversible dynamics, , is governed by


where is a linear reversible operator, which will be specified below. We consider two classes of operators: (DS) discrete state Markov jump process, where and assume only a finite number of states, and (FP) continuous Markovian diffusion governed by a Fokker-Planck equation.

For discrete Markov state models of class (DS) the dynamics is governed by


where denote the discrete states, is the rate of jumping from state to state () and is the total rate of leaving state guaranteeing conservation of probability (). In order to have reversible dynamics we need to additionally impose detailed balance, i.e. the constraint (see., e.g., van Kampen (2007)), which assures that the system will relax to a Boltzmann distribution in a potential on ergodic timescales , where is the inverse thermal energy. We call such a reversible ergodic process that conserves probability a relaxation process. If we add an absorbing point at we call the resulting process a first passage process or in short absorption, which we introduce in the following way. First, we modify the generator () such that all transitions corresponding to jumps out of the absorbing state are removed, i.e., the elements of the first passage generator read


Using a bra-ket matrix notation Kadanoff and Swift (1968) we rewrite this equation as


where is a vector with all entries except the th one; consequently, we identify . The first passage time density to reach state at time starting from is then formally defined by


which is nothing but the normalized probability flux into state with . Note that with Eq. (4) we use the convention that is the unique stationary solution with .

For a continuous space Markovian diffusion the transition probability density function (the ’propagator’) instead obeys the Fokker-Planck equation (2)


where is the probability current, is the diffusion constant, is a force field generated by the potential at position , and is the inverse temperature, which we set to to express energies in units of from now on. The scenario with reflecting barriers at with Redner (2001), we term a relaxation process, where correspond to so-called natural boundary conditions 3.

Conversely, an absorbing boundary at enters the Fokker-Planck equation via the Dirichlet boundary condition , without altering the partial differential equation (6), i.e., the first passage operator still reads . However, here the first passage time density becomes the probability flux into state . For convenience we use the operator as shorthand for Eq. (6) under the boundary condition .

We note that without an absorbing point both, dynamics governed by the master equation (2) and the Fokker-Planck equation (6) relax to the Boltzmann distribution , whereas with the absorbing boundary condition the particle will eventually reach the target with probability 1.

ii.2 Eigendecomposition

Since is assumed to generate a reversible Markov process, we can expand the generator in a bi-orthogonal eigenbasis Gardiner (2004). Denoting the eigenvalues of the relaxation process by and the corresponding left (right) eigenvectors by (), respectively, the generators from Eqs. (2) and (6) become in the respective eigenbases


where , and . We assume the eigenvalues to be ordered such that , and the generator to be irreducible , which means that there is a unique equilibrium state van Kampen (2007). Note that for a Fokker-Planck equation with reflecting barriers at (relaxation) the eigenfunction must satisfy the zero flux condition , with .

The generator with the absorbing point at state , can similarly be expanded in a bi-orthogonal set of eigenfunctions


where is the -th eigenvalue and () denote the corresponding left (right) eigenfunctions of the first passage process. Without loss generality we use an ordered labeling such that , where .

The left and right eigenvectors of the absorption (at position ) as well as of the relaxation process are related via and , respectively. In the case of a discrete number of states, the lowest eigenvalue of the generator (4) will be with the right eigenfunction , whereas for Fokker-Planck dynamics one imposes the boundary condition .

In our previous work an explicit Newton series expression for in terms of a series of almost triangular matrices was derived Godec and Metzler (2016b), which corresponds to a large deviation limit . One of our main goals here is to obtain the full first passage statistics explicitly in terms of relaxation eigenmodes. Our theory builds on the renewal theorem, which we briefly review in the following subsection.

ii.3 Renewal theorem

The classical renewal theorem provides a well known implicit connection between first passage and relaxation processes. It relates the probability density of the freely propagating system to be in state at time upon starting from a state , to the first passage distribution from to :


where both and admit a spectral representation




respectively. In other words, a system starting from state must pass through state before reaching the final state , which for an effectively 1-dimensional Fokker-Planck necessarily means or . In Eq. (11) we introduced in the first passage weights


for discrete state (DS) and Fokker-Planck (FP) dynamics, respectively, which must satisfy with the first nonzero weight being strictly positive , and where we introduced . Note that the first line of Eq. (12), i.e. the DS case, is equivalent to with from Eq. (4). In the case of FP dynamics the second line of Eq. (12) is equivalent to , which follows from a partial integration using both Eq. (6) and Eq. (8).

In the case of the renewal theorem (9) has the simple interpretation: a system being in state at time must have arrived at that point at some earlier time for the first time (), and then returned to the same position again at time , where corresponds to the time of first arrival.

Laplace transforming the renewal theorem (9), where a generic function is transformed according to , yields Siegert (1951)


Based on this well known renewal theorem we construct in the following section a method that allows to determine explicitly the first passage time statistics exactly in terms of the relaxation process.

Iii Principal result for discrete state systems

Starting from the renewal theorem (13), we now derive an expression for the first passage time density for discrete state Markov processes in terms of relaxation modes in the following three steps. The first step involves a crucial relation between the eigenvalues of the relaxation process and absorption process , which are here shown to interlace


for . For effectively one dimensional finite lattice models with the target at an outer edge these inequalities become strict


which will also apply identically to Fokker-Planck dynamics discussed in Sec. IV. In the second step we exactly express the first passage eigenvalues in the form of a Newton series of determinants of almost triangular matrices, which generalizes the result for the slowest mode from Godec and Metzler (2016b) to all first passage modes. The third and final step corresponds to a straightforward application of the residue theorem, which is used to determine the first passage weights .

iii.1 Interlacing of eigenmodes (step 1)

For a discrete system with states the eigenvalues and correspond to the roots of the respective characteristic polynomials


i.e., and . Inserting Eq. (4), which is , into the second characteristic polynomial (16) and using the matrix determinant lemma establishes a link between the two characteristic polynomials


where is called the adjugate of a matrix A satisfying Cramer’s rule . We note that the same mathematical concepts have been used recently to determine the stalling distribution of irreversibly driven systems (cf. deletion-contraction formula in Polettini and Esposito (2017); Bisker et al. (2017)).

The adjugate of a diagonal matrix D with elements ( if ) is diagonal as well, with elements . Consequently, the bi-orthogonal expansion (7) implies


which inserted into Eq. (17) gives


where we used the eigenvalue equation . Eq. (19) constitutes an essential step in our calculations, which allows us to express the diagonal of the relaxation propagator solely in terms of eigenvalues and (see the following subsection for more details).

Moreover, the characteristic polynomials of the first passage process and relaxation process change sign one after the other, since for all , which proves that the eigenvalues of the first passage process and eigenvalues of the relaxation process interlace according to Eq. (14). We note that this result is directly related to the interlacing of eigenvalues generated from a ”lumping” of states which is proven in Grone et al. (2008). In the following paragraph we briefly discuss the scenario, in which the interlacing of eigenvalues becomes strict (15), which will be the case for systems with Fokker-Planck dynamics discussed in Sec. IV.2.



Figure 1: Characteristic polynomials (solid blue line) and (dash-dotted red line) for simple three state models (see insets) along real axis in . (a) Linear chain of states with the absorbing state at the border . crosses the axis at which correspond to , respectively. (b) The same model as in (a) but with the absorbing state at . The first eigenmode vanishes at the target. (c) Fully connected three state model, in which . All rates in (a)-(c) are set to 1 (relaxation process). The transition rates away from absorbing point (dashed arrows) are set to zero (absorption). We note this special choice of rates deliberately generates a multiplicity of the first relaxation mode in (b) and the first passage mode in (c).

The stronger condition (15) holds if all eigenfunctions are nonzero at the target and all relaxation eigenvalues are non-degenerate, that is, for all . One can show that this condition is always trivially satisfied for 1-dimensional models ( if ), in which the target is placed at the border (e.g., or ); see inset of Fig. 1a for such an exemplary 3-state system.

Inserting the relaxation eigenvalues into the characteristic polynomial of the first passage process (19) yields


where we used the relations for all and for all , as well as . Since for each eigenvalue is positive (), the characteristic polynomial of the first passage process is proportional to . Consequently, changes sign exactly once between any two consecutive relaxation modes . The fact that and are polynomials of the same degree forbids more than a single root, and hence implies the strict interlacing of eigenvalues from Eq. (15), which completes the proof. The aforementioned reasoning is illustrated in Fig. 1a for a simple three state model in which the vertical arrows represent Eq. (20).

For fine-tuned systems in which the target is not located at the very outer position (see e.g., Fig. 1b) or systems that are not effectively one dimensional (see e.g., Fig. 1c) the strict interlacing theorem (15) can be violated, whereas the ”slightly weaker” interlacing condition (14) still holds.

iii.2 Diagonal of the relaxation propagator in terms of bare eigenvalues

Using the results from the previous subsection we are now in the position to represent (i.e. Eq. (10) with ), using only the eigenvalues of both the first passage process and the relaxation process, and , respectively. Laplace transforming the eigenmode expansion in Eq. (10) assuming yields


where we identified the equilibrium probability density in the first term. Comparing now in (21) with from Eq. (19) and from Eq. (16) yields after some algebra


The second equality in Eq. (22) follows from Eq. (16). Hence, encodes the eigenvalues of both, the relaxation and the first passage processes. Due to Eq. (21) contains only simple poles and decreases monotonically in between any two consecutive poles since . If (e.g., 1d models with the target at at the border) each root of represents a first passage eigenvalue (), which is located in between two relaxation modes , thus providing an alternative proof of relation (15) Keilson (1964). In the following section we determine the roots of the autocorrelation function explicitly, which due to Eq. (22) correspond to first passage eigenvalues .

Let us briefly reformulate in a way that can also be applied to continuous systems with an infinite number of states. Isolating the equilibrium probability, which is the first term in Eq. (21), from the product formula (22) yields


Since increase monotonically with we will later be able to adopt these results to systems governed by Fokker-Planck dynamics, which formally corresponds to the limit for which the product in Eq. (23) still converges. Eq. (23) was used to derive Eq. (10) in the accompanying letter.

iii.3 From the relaxation spectrum to the first passage time spectrum (steps 2 and 3)

Based on the interlacing theorem presented in Eq. (14), which is also given in Eq. (4) of the accompanying letter Note1 (), we can determine the full first passage time spectrum from the corresponding relaxation spectrum, . For simplicity we first consider the eigenvalues to be both ordered and non-degenerate, and also assume that holds for all values of . The extension to situations with , which also includes degenerate eigenvalues, for some is straightforward and will be dealt with at the end of this subsection.

Before determining the weights , we first determine the first passage eigenvalues , which were shown to be encoded in the roots of in Eq. (22). We introduce the th ”modified autocorrelation function”


which still encodes all of the first passage eigenvalues according to Eq. (22), i.e., it has exactly the same roots as . However, in contrast to the modified function is strictly concave within the interval , which can easily be confirmed by taking the second derivative and realizing that holds within the region of interest .

For and the modified functions and both are strictly concave within the interval and, consequently, also locally concave around the th first passage eigenvalue , i.e., and . Moreover, both functions and allow a Taylor expansion around the midpoint that converges within the whole interval including the root at which .

The method we present in the following is an analytical technique based on the principles of Newton iteration, which is a simple root finding algorithm that is guaranteed to work for functions that are both negative and concave between the starting point and the first root. Hence, to determine the th eigenvalue we accordingly choose the modified function


such that


which guarantees both negativity and concavity between and .

According to the interlacing theorem (15) is the only zero within the interval . With the midpoint starting condition the th first passage eigenvalue can be represented exactly in a series of determinants of almost triangular matrices


where is the th derivative of as defined in (25) with respect to at , and stands for an almost triangular matrix with elements Godec and Metzler (2016b)


with denoting the Heaviside step function ( if ) and . Moreover, we adopt the convention . We note that this method generalizes the method recently derived to determine the slowest first passage mode Godec and Metzler (2016b) to all first passage eigenmodes .

Let us briefly repeat the two crucial steps towards Eq. (27). First, the interlacing theorem (15) guarantees that the Taylor series around the midpoint coverges in the entire spectral interval , which also contains the first passage eigenvalue . Second, due to in Eqs. (24)-(26) the function is strictly concave and negative between and , which in turn guarantees the convergence of the explicit Newton series (27).

Eqs. (24)-(26) provide a universal method for determining explicitly first passage eigenvalues from the corresponding relaxation spectrum and constitute the central result of this work. We show in the Appendix A a simpler derivation of as well as a compact approximation of the principal first passage eigenvalue , which is particularly useful in the case of time scale separation (or ).

In the following we briefly comment on the practical implementation of the exact result for to render Eqs. (24)-(27) fully explicit. The weights will be determined afterwards in this subsection. The th derivative of with respect to at , , can be written explicitly as


where or is chosen according to Eq. (26). Note that condition (26) is equivalent to the condition , implying the first line of Eq. (29) to be either negative for or for , i.e. one has to evaluate the first line of Eq. (29) for : if one must to change to and reevaluate . Once one has determined and one can proceed with the second line of Eq. (29) to determine and insert the result in the almost triangular matrix (28). The determinant of almost triangular matrices can be calculated elegantly using the simple recursion relation from Cahill et al. (2002), see also Jia (2018) for an efficient numerical implementation.

Having obtained the first passage eigenvalues, the weights of the first passage time distribution can be calculated using the standard residue theorem. The Laplace transform of the spectral expansion of the first passage time density (11) reads


Using the residue theorem to invert the Laplace transformed renewal theorem (13) yields


where is taken at . The explicit Newton series (27) along with the first passage weights (31) fully characterize the first passage time distribution in terms of relaxation eigenmodes . This completes our third and final step, which allows, for the first time, to analytically deduce first passage time statistics directly from relaxation eigenmodes. We call this relation the explicit forward duality between first passage and relaxation. This completes the central result of this paper.

The spectral representation is very useful for determining the moments of the first passage time, . Moreover, as explained in more detail in the accompanying letter Note1 (), the full spectral expansion is required for a correct explanation of kinetics in the so-called few encounter limit, where molecules starting from position are searching for the target at . The probability density that the first molecule out of arrives at time at for the first time for this case becomes , which can be understood as follows. The probability that the first molecule have not yet reached the target will be given by , while the th particle arrives at with a rate ; hence the probability density that any particle out of molecules arrives at the target for the first time according to . Further details of the -particle problem and in particular the physical implications of the few-encounter limit are discussed in the accompanying letter Note1 ().

Let us now briefly generalize the method to systems with degenerate eigenvalues or vanishing relaxation modes. An eigenfunction that vanishes at the target will have a vanishing spectral weight as a result of Eq. (31). Hence, ’manually’ removing such modes will not affect the first passage time distribution . Moreover, if a relaxation eigenvalue is degenerate we define


and replace as well as and take the sums in Eq. (13) over all different values of . After renumbering all distinct contributing eigenvalues we obtain a strict interlacing (15). Therefore, we can apply our standard forward duality also to degenerate eigensystems. In the next subsection we will briefly derive a formal backward duality after which we reformulate the results from this subsection to continuous Fokker-Planck dynamics.

iii.4 Backward duality

In contrast to the explicit forward duality, which was presented in the previous subsection, an explicit reverse relation in the time-domain could not be established. In Laplace space, however, the forward duality can be inverted to give a backward duality as follows. Inserting the first passage generator from (4) into the Laplace transform of the propagator and using the Sherman-Morrison-Woodbury formula yields


Let us now insert the expression for the first passage time distribution from Eq. (5), which can be written as , into Eq. (33) to obtain


where is the generator of the relaxation process. Notably, this is expression corresponds to the backward duality and is the formal inverse of the renewal theorem, where is the Laplace transform of the cumulative first passage time distribution .

Iv Principal result for Fokker-Planck dynamics

iv.1 Greens function with natural boundaries

We restrict our discussion to effectively 1-dimensional dynamics, which include diffusion in dimensions in an isotropic potential as discussed in Godec and Metzler (2016b), where may also be fractal. Introducing an absorbing target at position splits the first passage problem into two cases (I) and (II) . Case (I) corresponds to an absorption from the left, and case (II) to an absorption from the right. In the following paragraph we demonstrate that all first passage modes of both distinct cases (I) and (II) are entirely encoded in , which allows to formulate the results from Sec. III.3 also for systems with Fokker-Planck dynamics.

Laplace transforming the Fokker-Planck equation (6) yields


where and is the initial position of the relaxation process. Eq. (35) is a inhomogeneous linear differential equation which can be solved using the standard Green’s function approach. First, we find the two independent solutions of the homogeneous problem , where we use the label ”” and ”” for the solution satisfying the left and right boundary condition, respectively. That is, a diffusion process within an interval imposes the probability current to vanish at the boundaries, i.e. . The special case of so-called natural boundary conditions correspond to the limit or analogously , that is, . The full solution of (35) is a continuous function in with a discontinuity of its first derivative at . Using the scaled Wronskian 4


the propagator, which satisfies the proper jump condition of the first derivative (current function) at , becomes


We note that the Wronskian (36) is proportional to the Boltzmann factor (see, e.g., Ref. Keilson (1964)), i.e., . Hence using the renewal theorem (13) and as well as from Eq. (37) yields the Laplace transform of the first passage time distribution