The Stretch to Stray on Time: Resonant Length of Random Walks in a Transien

First-passage times in random walks have a vast number of diverse applications in physics, chemistry, biology, and finance. In general, environmental conditions for a stochastic process are not constant on the time scale of the average first-passage time, or control might be applied to reduce noise. We investigate moments of the first-passage time distribution under a transient describing relaxation of environmental conditions. We solve the Laplace-transformed (generalized) master equation analytically using a novel method that is applicable to general state schemes. The first-passage time from one end to the other of a linear chain of states is our application for the solutions. The dependence of its average on the relaxation rate obeys a power law for slow transients. The exponent $\nu$ depends on the chain length $N$ like $\nu=-N/(N+1)$ to leading order. Slow transients substantially reduce the noise of first-passage times expressed as the coefficient of variation (CV), even if the average first-passage time is much longer than the transient. The CV has a pronounced minimum for some lengths, which we call resonant lengths. These results also suggest a simple and efficient noise control strategy, and are closely related to the timing of repetitive excitations, coherence resonance and information transmission by noisy excitable systems. A resonant number of steps from the inhibited state to the excitation threshold and slow recovery from negative feedback provide optimal timing noise reduction and information transmission.

While noise happens on the small length scales and short time scales of a system, it may trigger events on a global scale. One of the most important functions of noise for macroscopic dynamics arises from its combination with thresholds [38][39][40]. These are defined by the observation that the dynamics of a system remains close to a stable state as long as it does not cross the threshold value, and an actively amplified deviation from this state happens when the threshold is crossed. Noise drives the system across the threshold in a random manner. First passage is a natural concept to describe the timing of threshold crossings. Ignition processes are an illustrative textbook example. Although a small random spark might not be capable of igniting an inflammable material, a few of them might cause an explosion or forest fire. If the system again attains its stable state upon recovery from a deviation, such behavior is called excitable and the large deviation an excitation. The excitation is terminated by negative feedback. A forest is excitable, because it re- * Martin Falcke: Martin.Falcke@mdc-berlin.de grows after a fire. Consumption of inflammable trees acts as the negative feedback. Excitability describes not only forest fires but also the dynamics of heterogeneous catalysis [41], the firing of neurons [42], the properties of heart muscle tissue [42], and many other systems in physics, chemistry, and biology [43][44][45][46][47].
Random walks are frequently defined on a set of discrete states. The rates f i,j or waiting-time distributions Ψ i,j for transitions from state i to j set the state dwell times. The first-passage time between two widely separated states is much longer than the individual dwell times, and the conditions setting the rates and parameters of the Ψ i,j are likely to change or external control acts on the system between start and first passage. The conditions for igniting a forest fire change with the seasons or because the forest recovers from a previous fire. The occurrence of subcritical sparks during recovery has essentially no affect on the process of regrowth. More generally, noise does not affect recovery on large space and long time scales, and the random process experiences recovery as a slow deterministic change of environmental conditions. Since recovery is typically a slow relaxation process [45][46][47], it dominates event timing. Hence, first passage in an exponential transient is a natural concept through which to understand the timing of sequences of excitations. We will investigate it in this study.
We will take Markovian processes as one of the asymptotic cases of transient relaxation. There are several reasons for also considering non-Markovian formulations of continuous-time random walks. A description of a diffusing particle becomes non-Markovian whenever the particle not only passively experiences the thermal fluctuations causing its diffusive motion, but also acts back on its immediate surroundings [1,3,[5][6][7][8]48]. Non-Markovian waiting-time distributions for state transitions arise naturally in transport theory [30,[49][50][51]. Frequently in biological applications, the discrete states that arise are lumped states consisting of many "microscopic" states. Groups of open and closed states of ion channels lumped into single states are examples of this [52,53]. Transitions between lumped states are non-Markovian owing to the internal dynamics. We may also use waiting-time distributions if we lack information on all the individual steps of a process, but we do know the inter-event interval distributions. This is usually the case with in vivo measurements such as stimulation of a cell and the appearance of a response, or differentiation sequences of stem cells [37]. The state probabilities of non-Markovian processes obey generalized master equations, which we will use here [10,16,[54][55][56][57].
In Sec. II, we present a formulation of the general problem in terms of the normal and generalized master equations and give analytic solutions for both of these. These solutions apply to general state schemes. We continue with investigating first passage on linear chains of states in Sec. III. We present results on scaling of the average first-passage time with the relaxation rate of the transient γ and the chain length N in Sec. IV, and results on the phenomenon of resonant lengths in Sec. V.

A. The asymptotically Markovian master equation
In this section, we consider transition rates relaxing with rate γ to an asymptotic value λ i,j like They reach a Markov process asymptotically. The dynamics of the probability P i,j (t) to be in state j for a process that started in i at t = 0 obey the master equation In matrix notation with the vector of probabilities P i , we have with the matrices E and D defined by Eq. (2). The initial condition defines the vector r i = {δ ij }, j = 0, . . . , N . The Laplace transform of the master equation allows for a comfortable calculation of moments of the first-passage times, which we will carry out in Sec. III. The Laplace transform of Eq. (3) is the system of linear difference equations B. The generalized master equation

The waiting-time distributions
Waiting-time distributions Ψ j,k in a constant environment depend on the time t − t elapsed since the process entered state j at time t . The change in conditions causes an additional dependence on t: Ψ j,k (t, t−t ). The lumping of states, which we introduced as a major cause of dwell-time-dependent transition probabilities, often entails Ψ j, Waiting-time distributions used in transport theory exhibit similar properties [30,49,50,58]. We use a simple realization of this type of distributions by a biexponential function in t − t : The collect all factors of e −γt in Eq. (5). The functions g j,k (t − t ) describe the asymptotic part of the waitingtime distributions remaining after the transient. The Ψ j,k (t, t − t ) are normalized to the splitting probabilities C j,k = ∞ t dt Ψ j,k (t, t − t ) (the total probability for a transition from j to k given the system entered j at t ). They satisfy where N (j) out is the number of transitions out of state j.

The generalized master equation and its Laplace transform
In the non-Markovian case, the dynamics of the probabilities P i,j (t) obey a generalized master equation [59][60][61][62][63] where I i l,j (t) is the probability flux due to transitions from state l to j given that the process started at state i at t = 0, and N The second factor in the convolution is the probability of arriving in state l at time t , and the first factor is the probability of leaving toward j at time t given arrival at t . The q i l,j are the initial fluxes, with q i l,j (t) ≡ 0 for i = l.
The Laplace transform of the probability dynamics equation (8) is which contains the Laplace-transformed probability fluxesĨ i l,j . The Kronecker delta δ ij captures the initial condition.
Laplace-transforming Eq. (9) is straightforward for terms containing the asymptotic part g l,j (t − t ) of the Ψ l,j (t, t − t ) [see Eq. (6)], since they depend on t−t only and the convolution theorem applies directly. The terms containing the transient part h l,j (t − t )e −γt depend on both t − t and t and require a little more attention: This leads to the Laplace transform of Eq. (9): We write theĨ i l,j (s) andq i l,j (s) as the vectorsĨ i (s) and q i (s), respectively, to obtaiñ Solving forĨ i (s) results iñ This is again a system of linear difference equations. The entries in the matricesG(s) andH(s) are the functions g l,j (s) andh l,j (s), which are the Laplace transforms of g l,j (t − t ) and h l,j (t − t ) [Eq. (6)].

Solving the Laplace-transformed generalized and asymptotically Markovian master equations
All elements ofH(s) [Eq. (14)] vanish for s → ∞. We also expectĨ i (s) to be bounded for s → ∞. Hence, the solution for γ = ∞ is Consequently, for kγ 1, holds for the solutionĨ i (s) for all values of γ > 0 and for k a natural number. Once we knowĨ i (s + kγ), we can use Eq. (14) to findĨ i (s + (k − 1)γ): In this way, we can consecutively use Eqs. (14) and (15) to expressĨ i (s + (k − j)γ), j = 0, . . . , k, by known functions. The solution becomes exact with k → ∞. With the definitionÃ(s) as the solution of Eq. (14). Equation (17) is confirmed by verifying that it provides the correct solution in the limiting cases γ = 0 and γ = ∞. In the latter, all entries of the matrixH(s) vanish, and we obtain directly the result without transient, Eq. (15). For γ = 0, we notice that holds, which leads to the correct solution of Eq. (14): We now turn to the asymptotically Markovian master equation (2) and write its Laplace transform (4) as This equation has the same structure as Eq. (14), and since the matrix (1s − E) −1 also vanishes for s → ∞, we can calculate the Laplace transform of the P i,j completely analogously to the non-Markovian case. We de- as the solution of Eq. (19). These solutions for the Laplace transforms of the generalized and asymptotically Markovian master equations with a transient, Eqs. (14) and (20), are not restricted to state-independent waiting-time distributions or rates. They also apply to random walks with space-or statedependent waiting-time distributions and to arbitrary state networks.

III. THE PROBABILITY DENSITY OF THE FIRST-PASSAGE TIME FOR A LINEAR CHAIN OF STATES
A large class of stochastic processes, including all the examples mentioned in the introduction, are represented by state schemes like which we consider from now on. The definition of the transition probabilities for both the asymptotically Markovian system and the non-Markovian system completes its specification. We use the asymptotically Markovian rates . The process has a strong initial bias for motion toward 0, We specify the input for the generalized master equation of the non-Markovian system as The denominators in Eqs. (23) and (24) arise from the normalization to . We use identical waiting-time distributions for all transitions i → i + 1 (i > 0) and all transitions i → i − 1. The process defined by Eqs. (23) and (24) also has a strong bias toward 0 in the beginning, and is approximately symmetric in the limit t → ∞. Figure 1 illustrates the Ψ i,i±1 and their development with increasing t . The limit of γ α, β, , δ illustrates the consequences of the transient. The C i,k relax from 0 to an asymptotic value with rate γ, but the dwell times of individual states i (i = 0) are essentially constant (Fig. 1).
There is only one transition away from 0. That changes the normalization to C 0,1 (t ) = 1 and we cannot use the form of Eqs. (23) and (24) for Ψ 0,1 . We use instead The transient here causes a decrease of the dwell time in state 0 (see Fig. 1). The first-passage-time probability density F 0,N (t) provides the probability of arrival for the first time in state N in (t, t+dt) when the process started in state 0 at t=0. It is given by the probability flux out of the state range from 0 to N − 1: We denote its Laplace transform byF 0,N (s). The moments of the first-passage-time distribution are given by [15] F 0,N (t) can be determined by solving the master equations, setting state 0 as the initial condition and considering the state N as absorbing, i.e., F 0,N (t) captures not only the first-passage time 0 → N , but also, to a good approximation, transitions starting at states with indices larger than 0 (and < N ), since, owing to the initial bias, the process quickly moves into state 0 first and then slowly starts from there. With the non-Markovian waiting-time distributions, the matricesG(s) andH(s) and the vectorq i (s) of Eq. (17) specific to this problem arẽ All other entries are equal to 0. Specifically for the firstpassage problem, the matrices D and E are with all other entries being 0. The vectorr is equal to δ 1i , i = 1, . . . , N − 1. The Laplace transforms of the first-passage-time distributions arẽ with Eq. (8) and with Eq. (2). Figure 2 (a) and (b) compare analytical results for the average first-passage time T with the results of simulations. The agreement is very good, thus confirming the solutions given by Eqs. (17) and (20). This confirmation by simulations is important, since there is no method of solving difference equations that guarantees a complete solution.
IV. SCALING OF THE AVERAGE FIRST-PASSAGE TIME WITH THE RELAXATION RATE Figure 2 shows results for the average first-passage time T across four orders of magnitude of the relaxation rate γ. The results strongly suggest that T grows according to a power law γ −ν with decreasing γ, if γ is sufficiently small. The exponent exhibits a simple dependence on the number N ν of edges with identical waitingtime distributions. That number is equal to N for processes according to Eqs. (22) and to N − 1 for processes obeying Eqs. (23)- (25), since the 0 → 1 transition is different there. The exponent ν depends to leading order on N ν like N ν /(N ν + 1). This applies to both the asymptotically Markovian and non-Markovian waiting-time distributions, and to both parameter sets of waiting-time distributions simulated in the non-Markovian case.
The exponent ν is equal to 1 2 for N ν = 1, as has previously been shown analytically for f i,i±1 according to Eq. (22) (see [64], Chapter 5). The process is very unlikely to reach large N with a bias toward 0, even if this is only small. Hence, the random walk "waits" until the transient is over and symmetry of the transition rates has been reached [see Fig. 1(d)], and then goes to N . This waiting contributes a time ∝ γ −1 to T , and ν approaches 1 for large N ν .
The average first-passage time for a symmetric random walk increases with N like N (N − 1) [65], i.e., it is very long for large N . We see a contribution of the transient to T only if relaxation is slow enough for γ −1 to be comparable to this long time. Consequently, T is essentially independent of γ for large γ and large N (see N = 200 in Fig. 2(a, b)).

V. RESONANT LENGTH
The coefficient of variation CV [= standard deviation (SD)/average T ] of the first-passage time reflects the relative fluctuations. Its dependence on the chain length N is illustrated in Fig. 3. CV increases monotonically with increasing γ. Its dependence on N is not monotonic. We find a pronounced minimum of CV(N ) for small values of γ, where the minimal value is up to one order of magnitude smaller than CVs with N = 1 and with large N . Both simulations and analytical results show this behavior.
How can we get a heuristic understanding of the decrease in CV with decreasing γ and initially with increasing N ? The lingering close to state 0 shown in Fig. 1(d) means that states with index i larger than 1 are not reached before a time t ≈ γ −1 with almost certainty. This initial part of the process contributes to T but little to SD. Hence, its growth with decreasing γ and (initially) increasing N decreases CV. Additionally, the standard deviation is determined by the rates and splitting probabilities at the end of the transient, which are more favorable for reaching N than the splitting probabilities at the beginning. It is therefore less affected than the average by the γ values. This also causes a decrease in CV with decreasing γ.
To gain more insight into the N dependence of CV, we look at processes with constant transition rates and a symmetric process with transient rates but constant splitting probabilities C i,i−1 = C i,i+1 = 1 2 [see Eq. (A1)]. The permanently symmetric process exhibits only a very shallow minimum, and the minimal value of CV appears to be independent of γ [ Fig. 4(a): discrete case]. Processes with constant transition rates exhibit decreasing CV with decreasing C i,i−1 [65] [ Fig. 4(b)]. However, this decrease is minor and does not explain the data in Fig. 3 in the range C i,i−1 ≥ 1 2 covered by the transient. The comparison shows that the dynamics of the splitting probabilities provide the major part of the reduction in CV with increasing N in Fig. 3. The splitting probability C i,i−1 relaxes from 1 to its asymptotic value set by the parameters of the process. The larger the value of N , the later is the absorbing state reached and the smaller is the value of C i,i−1 when it is reached. As long as C i,i−1 (t = T ) decreases sufficiently rapidly with increasing N , so does CV. At values of N such that γT 2, the decrease in C i,i−1 (t = T ) is negligible, and CV starts to rise again with increasing N toward its large-N value.
These considerations suggests that the N with minimal CV could be fixed by a specific value of the splitting probability. This value of C i,i−1 cannot be smaller than 1 2 , since we would expect a monotonically decreasing CV in that regime [the N − 1 2 regime in Fig. 3(d)]. Since this should hold also for minima with N 1, this value of C i,i−1 also cannot be much larger than 1 2 , since T would then diverge. This is confirmed by the results shown in Fig. 3(f). The splitting probability C i,i−1 (t = T ) at the minimum changes by less than 8% over four orders of magnitude of γ, and is slightly larger than 1 2 . Hence, CV starts to rise again when the length N is so large that the average first-passage time is long enough for C i,i−1 (t = T ) to approach the symmetric limit [66].
The value of CV(N = ∞) depends on γ in the case with non-Markovian waiting-time distributions. It is 1 for large γ values, since the asymptotic splitting probability C i,i−1 (t = ∞) is larger than  The parameter values in Fig. 3 The process has a bias toward N and we see the well-known behavior CV(N ) ∝ N − 1 2 for large N . The onset of the N − 1 2 behavior moves to smaller N with increasing values of γ until finally the minimum of the CV is lost. We found minima if γ −1 (four times the average state dwell time), but we have not determined a precise critical value.
The transition from the CV of Ψ 0,1 to CV(N = ∞) for large values of γ is monotonic in Fig. 3(c) for the asymptotically non-Markovian case. The asymptotically Markovian system exhibits a minimum even for large γ values [γ ≈ λ: Fig. 3(e)]. However, it is comparably shallow, and is in line with the results of Jouini and Dallery for systems without transient [65] [ Fig. 4(b)].

A. (Generalized) master equation with exponential time dependence of rates and waiting-time distributions
In general, stochastic processes occur in changing environments or may be subjected to control. The corresponding mathematical description is provided by (generalized) master equations with time-dependent coefficients. The analytical solution of these equations is complicated even in the simplest case [see Eq. (A8)] [67]. The Laplace-transform approach is not applicable in general, but the exponential time dependences of the master equation (2) or the generalized master equation (8) with (9) do allow the transform to be performed on them, leading to difference equations in Laplace space. The solutions of these equations are given by Eqs. (17) and (20), resp. To the best of our knowledge, this is a novel (and efficient) method of solving the generalized master equation with this type of time dependence in the kernel, or the master equation with this type of time-dependent rates.
Equations (17) and (20) also apply in the case of state schemes more complicated than that in (21). It is only necessary to make appropriate changes to the definitions ofH(s) andG(s) or of D and E. Hence, we can use the method presented here to investigate the dynamics of complicated networks in a transient. This is of interest for studies of ion-channel dynamics [53,68], transport in more than one spatial dimension, gene regulatory networks, and many other applications. Many quantities of interest in these investigations are moments of first-passage type and thus can be calculated from the derivatives of the Laplace transforms. Some of these calculations might require knowledge of the probability time dependence P i,j (t), i.e., the residues ofP i,j (s). In many cases, determination of these residues will be only slightly more complicated than for the system without transient, since all terms arising from the transient involve the same factors, but with shifted argument, together with the ma-trixH(s). In general, if we know the set of residues {S 0 } of the Laplace transform of the system without transient and the matrixH(s), we can immediately write down the residues of the Laplace transform of the system with transient by shifting {S 0 } by integer multiples of γ. This allows easy generalization of many results and might be of particular interest for renewal theory [15,69].

B. The average first-passage time
The average first-passage time decreases with increasing relaxation rate according to a power law, if the transient is sufficiently slower than the time scale of the individual steps (Fig. 2). The exponent depends on the chain length to leading order like −N/(N + 1). Remarkably, this behavior has been found for substantially different parameters in both the asymptotically non-Markovian and asymptotically Markovian cases, and therefore appears to be rather universal.

C. Resonant length
Transients can substantially reduce the coefficient of variation. We have found a monotonic decrease of CV with decreasing γ when the relaxation rate is slow compared with state transition rates. CV exhibits a minimum in its dependence on the chain length N (Fig. 3) if the transient is slow compared with state dynamics. This minimum exists for both asymptotically Markovian and non-Markovian discrete systems. At a time-scale separation between the state transition rates and the transient of about 10 6 , the CV at the minimum is reduced by about one order of magnitude compared with CV(N = 1).
Exploiting our results for control purposes, we can use a transient if more precise arrival timing at the outcome of a process is desired. Applying a transient that relaxes a bias toward 0 reduces CV. In the context of charge carrier transport, such a transient could be realized by a time-dependent electric field. If we include the rising part of CV(N ) beyond the minimum in our consideration and consider all CVs smaller than CV(N = 1) as reduced, we may still see a reduction of CV even if the average first-passage time is an order of magnitude longer than γ −1 . Hence, the initial transient may have surprisingly long-term consequences and is therefore a rather robust method for CV reduction. Additionally, or if a transient is given, the state and step number can be used for optimizing the precision of arrival time. The optimum may also be at smaller lengths than that of the non-optimized process; i.e., optimization of precision may even lead to acceleration.

D. How can negative feedback robustly reduce noise in timing?
Most studies investigating the effect of negative feedback have focused on amplitude noise-we are looking at noise in timing. Timing noise substantially reduces information transmission in communication systems. Variability in timing and protein copy numbers in gene expression or cell differentiation causes cell variability [37,70]. Therefore, many studies have investigated the role of noise in gene regulatory networks and have found that immediate negative feedback is not suitable for reducing timing noise [71][72][73]. This agrees with our results for large γ values.
We have shown that slow recovery from negative feedback is a robust means of noise reduction. The recovery process corresponds to the slow transient in our present study. Negative feedback terminating the excited state is a constitutive element of many systems generating sequences of pulses or spikes, such as oscillating chemical reactions [74,75], the electrical membrane potential of neurons [76], the sinoatrial node controlling the heart beat [77], the tumor suppressor protein p53 [78,79], Ca 2+ spiking in eukaryotic cells [80][81][82], cAMP concentration spikes in Dictyostelium discoideum cells and other cellular signaling systems [83,84].
In particular, our results apply to noisy excitable spike generators in cells, where single molecular events or sequences of synaptic inputs form a discrete chain of states toward the threshold. Beside examples from membrane potential spike generation [85], Ca 2+ spiking and p53 pulse sequences have been shown to be noise-driven excitable systems [79,81,82,86]. The information content of frequency-encoding spiking increases strongly with decreasing CV of spike timing [87]. A value of the coefficient of variation between 0.2 and 1.0 has been measured for the sequence of interspike intervals of intracellular Ca 2+ signaling [81,82,[88][89][90]. Hence, CV is decreased compared with that of a Poisson process and even compared with that for first passage of a symmetric random walk. The experimental data are also compatible with the finding that the slower the recovery from negative feedback, the lower is CV [81,82,[88][89][90]. This strongly suggests that this ubiquitous cellular signaling system uses the mechanism of noise reduction described here to increase information transmission.
Optimal noise amplitudes minimize the CV of interspike intervals in noisy excitable systems [40,[91][92][93][94], which has been termed coherence resonance. Our results define the conditions of optimal CV reduction in terms of system properties-an optimal number of steps from the inhibited state to the excitation threshold during slow recovery. At the same time, our results shed light on new aspects of coherence resonance, and indeed may indicate that a more fundamental phenomenon underlies it. Since we believe excitable systems to be one of the most important applications of our results, we have chosen the term resonant length.
Coming back to the widely used graphic example, what can the drunkard learn from our results? If he chooses a pub at the right distance from home, he will arrive home sober and relatively in time for breakfast. The permanently symmetric case with a transient is defined by Its continuum limit on a domain of length L (with spatial coordinate x) satisfies The Fokker-Planck equation (A2) can be solved after applying a time transformation t → τ : The boundary and initial conditions, P (L, t) = 0, ∂P (x, t)/∂x| x=0 = 0, and P (x, 0) = δ(x) specify the firstpassage problem. We find where k n = π(2n + 1)/(2L). The ith moment of the first-passage time is given by With a n = 2λk 2 n /γ, we obtain (−1) n 2n + 1 e an a −an n Γ (a n , a n ) , where Γ(y, x) is the lower incomplete gamma function. The second moment is (−1) n (2n + 1) 5 e an × 2 F 2 ({a n , a n } ; {a n + 1, a n + 1} ; −a n ) (A8) where 2 F 2 is a hypergeometric function. Results for CV are shown in Fig. 4(a). In the continuous case, CV does not exhibit a minimum in its dependence on the length L. Interestingly, CV at small L is very close to the minimum value of the discrete case.
The a n are small for L λ/γ. Therefore, in that limit, we can approximate the incomplete Γ function and hypergeometric function by Γ (a n , a n ) ≈ a an−1 n e −an , e an 2 F 2 ({a n , a n } ; {a n + 1, a n + 1} ; −a n ) ≈ 1, and find i.e., CV approaches monotonically the value of a symmetric random walk with constant transition rates for large L.

Appendix B: Numerical methods
Evaluation of Eqs. (17) and (20): If N is large and the ratio between the fastest state transition rate and the relaxation rate γ is larger than 10 4 , high precision of the numerical calculations is required for the use of Eqs. (17) and (20). A priori known values likeĨ 0 N −1,N (s = 0) = 1 can be used to monitor the precision of the calculations. The matrix products become very large at intermediate values of jγ during the summation in Eqs. (17) and (20), and their sign alternates such that two consecutive summands nearly cancel. Intermediate summands are of order larger than 10 17 , and thus we face loss of significant figures even with the numerical floating-point number format long double. We used Arb, a C numerical library for arbitrary-precision interval arithmetic [95], to circumvent this problem. It allows for arbitrary precision in calculations with Eqs. (17) and (20). Computational speed is the only limitation with this library and has determined the parameter range for which we established analytical results. We were able to go to a time-scale separation of ≈ 10 6 using this library.
Simulation method : We use simulations for comparison with the analytic solutions and for systems with large N . The Ψ l,j (t, t − t ) are functions only of t for a given t . Since t is the time when the process moved into state l, it is always known. The simulation algorithm generates in each iteration the cumulative distribution t t dθ N (l) out i=1 Ψ l,ji (θ, θ−t ) and then draws the time t tr for the next transition. The specific transition is chosen from the relative values Ψ l,j k (t tr , t tr − t )/ N (l) out i=1 Ψ l,ji (t tr , t tr − t ) by a second draw. We verified the simulation algorithm by comparison of simulated and calculated stationary probabilities [Eq. (C4)], splitting probabilities at a variety of t values, simulations without recovery ("γ = ∞") and the comparisons in Fig. 2. We used at least 20 000 sample trajectories to calculate moments and up to 160 000 to determine the location of the minima of CV.

Appendix C: The stationary probabilities
The stationary probabilities are reached for t → ∞. This entails Ψ i,i±1 (t, t − t ) = Ψ ∞ i,i±1 (t − t ) = g i,i±1 (t − t ). We start from the idea that the stationary probability P i of being in state i is equal to the ratio of the total average time T i spent in i divided by the total average time for large t: T i is equal to the number N i of visits to state i multiplied by the average dwell time t i in i: Each visit to state i starts with a transition to i. The average number of transitions into i is equal to the number of visits to its neighboring states multiplied by the probability that the transition out of the neighboring states is toward i. With the splitting probabilities C i,i±1 (t = ∞) denoted by C i,i±1 , the N i obey Dividing by the total number of visits N j=0 N j , we get and all other entries 0. We checked for N = 2, 3, 4 that det M = 0 holds. Equation (C4) determines n only up to a common factor, which is then fixed by Eq. (C5). The average dwell times t i in the states are initially affected by the recovery from negative feedback, but are constant at large t. They have contributions t i,i±1 from both transitions to i ± 1, with weights set by C i,i±1 : (C9) C i,i±1 t i,i±1 = −(∂/∂s)g i,i±1 | s=0 holds owing to the normalization of the Ψ i,i±1 . This leads finally to To give a specific example, the stationary-state probabilities for N = 2 are P = 1 (∂/∂s) (C 1,0g0,1 +g 1,0 +g 1,2 + C 1,2g2,1 ) ×    C 1,0 (∂/∂s)g 0,1 (∂/∂s) (g 1,0 +g 1,2 ) C 1,2 (∂/∂s)g 2,1    s=0 .