A framework for the direct evaluation of large deviations in non-Markovian processes

The WTD can be expressed in terms of a time-dependent rate or hazard ${\beta }_{{x}_{n},{x}_{n-1}}[{t}_{n}-{t}_{n-1};w({t}_{n-1})]$, which is the probability density that there is a jump from ${x}_{n-1}$ to x_n in $[{t}_{n},{t}_{n}+{\rm{d}}t)$, conditioned on having no transitions during the interval $[{t}_{n-1},{t}_{n})$,

$$\begin{eqnarray}&&{\psi }_{{x}_{n},{x}_{n-1}}[{t}_{n}-{t}_{n-1};w({t}_{n-1})]\\ &&\quad ={\beta }_{{x}_{n},{x}_{n-1}}[{t}_{n}-{t}_{n-1};w({t}_{n-1})]\times {\phi }_{{x}_{n-1}}[{t}_{n}-{t}_{n-1};w({t}_{n-1})].\end{eqnarray} \tag{ 10 }$$

For brevity we define $\tau ={t}_{n}-{t}_{n-1}$, which is the value of the age, i.e. the time elapsed since the last jump, at which the next jump takes place (see figure 1). Roughly speaking, the hazard ${\beta }_{{x}_{n},{x}_{n-1}}[\tau ;w({t}_{n-1})]$ is the likelihood of having an almost immediate transition from a state ${x}_{n-1}$ known to be of age τ, to a state x_n, and, crucially, can also depend on the history $w({t}_{n-1})$. From equation (10), summing over ${x}_{n}\in { \mathcal S }$, and defining ${\beta }_{{x}_{n-1}}[\tau ;w({t}_{n-1})]$ $=\,{\sum }_{{x}_{n}}{\beta }_{{x}_{n},{x}_{n-1}}[\tau ;w({t}_{n-1})]$ and ${\psi }_{{x}_{n-1}}[\tau ;w({t}_{n-1})]$ $=\,{\sum }_{{x}_{n}}{\psi }_{{x}_{n},{x}_{n-1}}[\tau ;w({t}_{n-1})]$, we get

$$\begin{eqnarray}{\psi }_{{x}_{n-1}}[\tau ;w({t}_{n-1})] & = & -\displaystyle \frac{{\rm{d}}}{{\rm{d}}\tau }{\phi }_{{x}_{n-1}}[\tau ;w({t}_{n-1})]\\ & = & {\beta }_{{x}_{n-1}}[\tau ;w({t}_{n-1})]\times {\phi }_{{x}_{n-1}}[\tau ;w({t}_{n-1})].\end{eqnarray} \tag{ 11 }$$

The age-dependent sum ${\beta }_{{x}_{n-1}}[\tau ;w({t}_{n-1})]$ is also referred to as the escape rate from ${x}_{n-1}$ and, using equation (11), can be written as a logarithmic derivative

$$\begin{eqnarray}&&{\beta }_{{x}_{n-1}}[\tau ;w({t}_{n-1})]=-\displaystyle \frac{{\rm{d}}}{{\rm{d}}\tau }\mathrm{ln}{\phi }_{{x}_{n-1}}[\tau ;w({t}_{n-1})].\end{eqnarray} \tag{ 12 }$$

Integrating equation (12) with initial condition ${\phi }_{{x}_{n-1}}[0;w({t}_{n-1})]=1$ gives:

$$\begin{eqnarray}&&{\phi }_{{x}_{n-1}}[\tau ;w({t}_{n-1})]=\exp \left(-{\int }_{0}^{\tau }{\beta }_{{x}_{n-1}}[t;w({t}_{n-1})]{\rm{d}}t\right),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\psi }_{{x}_{n-1}}[\tau ;w({t}_{n-1})]={\beta }_{{x}_{n-1}}[t;w({t}_{n-1})]\exp \left(-{\int }_{0}^{\tau }{\beta }_{{x}_{n-1}}[t;w({t}_{n-1})]{\rm{d}}t\right).\end{eqnarray} \tag{ 14 }$$

These let us cast equation (10) in the more convenient form

$$\begin{eqnarray}&&{\psi }_{{x}_{n},{x}_{n-1}}[\tau ;w({t}_{n-1})]={p}_{{x}_{n},{x}_{n-1}}[\tau ;w({t}_{n-1})]\times {\psi }_{{x}_{n-1}}[\tau ;w({t}_{n-1})],\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{p}_{{x}_{n},{x}_{n-1}}[\tau ;w({t}_{n-1})]=\displaystyle \frac{{\beta }_{{x}_{n},{x}_{n-1}}[\tau ;w({t}_{n-1})]}{{\displaystyle \sum }_{{x}_{n}}{\beta }_{{x}_{n},{x}_{n-1}}[\tau ;w({t}_{n-1})]},\end{eqnarray} \tag{ 16 }$$

is the probability that the system jumps into the state x_n, given that it escapes the state ${x}_{n-1}$ at age τ. It is important to notice that the normalisation conditions ${\sum }_{{x}_{n}}{p}_{{x}_{n},{x}_{n-1}}[\tau ;w({t}_{n-1})]=1$ and ${\int }_{0}^{\infty }{\psi }_{{x}_{n-1}}[\tau ;w({t}_{n-1})]{\rm{d}}\tau =1$ are satisfied. Hence, we can sample a random waiting time τ, according to the density ${\psi }_{{x}_{n-1}}[\tau ;w({t}_{n-1})]$ and, after that, a random arrival configuration x_n, according to the probability mass ${p}_{{x}_{n},{x}_{n-1}}[\tau ;w({t}_{n-1})]$. This suggests a standard Monte Carlo algorithm for the generation of a trajectory (1):

Zoom In Zoom Out Reset image size

• (1)  
  Initialise the system to a configuration x₀ and a time t₀. Set a counter to n = 1.
• (2)  
  Draw a value τ according to the density (11) and update the time to ${t}_{n}={t}_{n-1}+\tau$.
• (3)  
  Update the system configuration to x_n, with probability given by (16).
• (4)  
  Update n to $n+1$ and repeat from (2) until t_n reaches the desired simulation time.

We now need to take into account the effect of the factor ${{\rm{e}}}^{-s{\theta }_{{x}_{n},{x}_{n-1}}}$ on the dynamics, which is to increment (if ${\theta }_{{x}_{n},{x}_{n-1}}\lt 0$) or decrement (if ${\theta }_{{x}_{n},{x}_{n-1}}\gt 0$) the 'weight' of a trajectory, within an ensemble. This can be implemented by means of the so-called 'cloning' method, which consists of assigning to each trajectory of the ensemble a population of identical clones proportional to its weight. As reviewed, e.g., in [23, 24], the idea is not new. It seems to be born in the context of quantum physics and can be traced back to Enrico Fermi [25]. Noticeably, cloning was proposed as a general scheme for the evaluation of large deviation functionals of non-equilibrium Markov processes in [10] (with further continuous-time implementation in [11]) and has also been extensively applied within equilibrium statistical physics [26–28]. Here, we refine and extend this idea for the case of non-Markovian processes. One of the devices used in [10, 11] is to define modified transition probabilities—valid only under the Markovian assumption—and a modified cloning factor, encoding the contraction or expansion of the trajectory weight. In fact, it is implicit in the original work that the redefinition of such quantities is unnecessary; in some cases it may also be inconvenient, see section 4. An arguably more natural choice, especially for non-Markovian dynamics, is to focus on the WTDs. Specifically, equations (8) and (9) suggest the following procedure:

• (1)  
  Set up an ensemble of N clones and initialise each with a given time t₀, a random configuration x₀, and a counter n = 0. Set a variable C to zero. For each clone, draw a time τ until the next jump from the density ${\psi }_{{x}_{0}}[\tau ;w({t}_{0})]$, and then choose the clone with the smallest value of t = t₀ + τ.
• (2)  
  For the chosen clone, update n to $n+1$ and then the configuration from ${x}_{n-1}$ to x_n according to the probability mass ${p}_{{x}_{n},{x}_{n-1}}[\tau ;w(t-\tau )]$.
• (3)  
  Generate a new waiting time τ for the updated clone according to ${\psi }_{{x}_{n-1}}[\tau ;w(t)]$ and increment its value of t to $t+\tau$.
• (4)  
  Cloning step. Compute $y=\lfloor {{\rm{e}}}^{-s{\theta }_{{x}_{n},{x}_{n-1}}}+u\rfloor$, where u is drawn from a uniform distribution on $[0,1)$.

  • (1)  
    If y = 0, prune the current clone. Then replace it with another one, uniformly chosen among the remaining $N-1$.
  • (2)  
    If $y\gt 0$, produce y copies of the current clone. Then, prune a number y of elements, uniformly chosen among the existing N + y.
• (5)  
  Increment C to $C+\mathrm{ln}[(N+{{\rm{e}}}^{-s{\theta }_{{x}_{n},{x}_{n-1}}}-1)/N]$. Choose the clone with the smallest t, and repeat from (2) until t − t₀ for the chosen clone reaches the desired simulation time T.

The SCGF is finally recovered as −C/T for large T. The net effect of step (4) is to maintain a constant population of samples whose mean current does not decay to $\langle j\rangle$.

The current through an ion channel in a cellular membrane can be modelled with only two states, corresponding to the gate being singly occupied (x₁) or empty (x₀); an ion can enter or leave this channel via the left (L) or right (R) boundary and non-exponential waiting times lead to a complex behaviour [43]. Specifically, we denote the WTD for a particle succeeding in entering (or leaving) through the boundary L by ${\psi }_{{x}_{1},{x}_{0},1}(\tau )$ (or ${\psi }_{{x}_{0},{x}_{1},-1}(\tau )$) with respective density ${\psi }_{{x}_{1},{x}_{0},0}(\tau )$ (or ${\psi }_{{x}_{0},{x}_{1},0}(\tau )$) for the boundary R. The rightwards current is measured by a counter that increases (decreases) by one when a particle enters (leaves) the system through the boundary L. Its exact SCGF is obtained numerically in [38] as the leading pole of the time-Laplace transform of $Z(s,t)$, for the DTI-case ${\psi }_{{x}_{i},{x}_{j},c}(\tau )={p}_{{x}_{i},{x}_{j},c}{\psi }_{{x}_{j}}(\tau )$ with ${\sum }_{c=-\mathrm{1,0}}{p}_{{x}_{0},{x}_{1},c}={\sum }_{c=\mathrm{0,1}}{p}_{{x}_{1},{x}_{0},c}=1$, and the particular choice ${\psi }_{{x}_{j}}(\tau )=g(\tau ;{k}_{j},{\lambda }_{j})$, where

$$\begin{eqnarray}&&g(\tau ;k,\lambda )={\lambda }^{k}{\tau }^{k-1}\exp (-\lambda \tau )/{\rm{\Gamma }}(k),\end{eqnarray} \tag{ 31 }$$

with ${\rm{\Gamma }}(k)$ as the Gamma function; the Markovian case is recovered for k = 1. Notably, the cloning method of section 3 can be implemented for any WTD, as only a bias of ${{\rm{e}}}^{s}$, for ions leaving the channel leftwards, and a bias of ${{\rm{e}}}^{-s}$, for ions entering from left, are needed. Figure 2(a) shows that this method reproduces, within numerical accuracy, the solution given in [38].

Zoom In Zoom Out Reset image size

In a quest for more general classes of models to illustrate the power of our approach, we now relax the constraint of DTI and assume that each transition can be triggered independently by two mechanisms, corresponding to the two boundaries. We still assume that memory of the previous history is lost as soon as the system changes state, thus preserving the semi-Markov nature. At the instant when the gate is emptied, a particle attempts to enter the system from the left boundary after a waiting time ${T}_{0}^{L}$ with density distribution ${{\boldsymbol{\psi }}}_{0}^{L}(\tau )$, while another particle attempts to arrive from the right boundary after a time ${T}_{0}^{R}$ distributed according to ${{\boldsymbol{\psi }}}_{0}^{R}(\tau )$. The waiting times ${T}_{1}^{L}$ and ${T}_{1}^{R}$, as well as the densities ${{\boldsymbol{\psi }}}_{1}^{L}(\tau )$ and ${{\boldsymbol{\psi }}}_{1}^{R}(\tau )$ are defined similarly. In order to have a right (left) jump during the interval $[\tau ,\tau +{\rm{d}}\tau )$, we also require that the left (right) mechanism remains silent until time τ. Consequently, the WTDs are

$$\begin{eqnarray}&&{\psi }_{{x}_{1},{x}_{0},0}(\tau )={{\boldsymbol{\psi }}}_{0}^{R}(\tau ){{\boldsymbol{\varphi }}}_{0}^{L}(\tau ),\qquad {\psi }_{{x}_{1},{x}_{0},-1}(\tau )={{\boldsymbol{\psi }}}_{0}^{L}(\tau ){{\boldsymbol{\varphi }}}_{0}^{R}(\tau ),\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{\psi }_{{x}_{0},{x}_{1},0}(\tau )={{\boldsymbol{\psi }}}_{1}^{R}(\tau ){{\boldsymbol{\varphi }}}_{1}^{L}(\tau ),\qquad {\psi }_{{x}_{0},{x}_{1},1}(\tau )={{\boldsymbol{\psi }}}_{1}^{L}(\tau ){{\boldsymbol{\varphi }}}_{1}^{R}(\tau ),\end{eqnarray} \tag{ 33 }$$

where ${{\boldsymbol{\varphi }}}_{j}^{(\rho )}(\tau )={\int }_{\tau }^{\infty }{{\boldsymbol{\psi }}}_{j}^{(\rho )}(t){\rm{d}}t$ are survival probabilities, with ρ denoting the mechanism L or R. As a concrete choice, we again assign a Gamma probability distribution to the waiting time of each event

$$\begin{eqnarray}&&{{\boldsymbol{\psi }}}_{j}^{(\rho )}(\tau )=g(\tau ;{k}_{j}^{(\rho )},{\lambda }_{j}^{(\rho )}),\end{eqnarray} \tag{ 34 }$$

so that the survival probabilities are

$$\begin{eqnarray}&&{{\boldsymbol{\varphi }}}_{j}^{(\rho )}(\tau )={\rm{\Gamma }}({k}_{j}^{(\rho )},{\lambda }_{j}^{(\rho )}\tau )/{\rm{\Gamma }}({k}_{j}^{(\rho )}),\end{eqnarray} \tag{ 35 }$$

where ${\rm{\Gamma }}(k,x)$ is the upper incomplete Gamma function. The time to the next jump, given that the system just reached state x_j (i.e., its age is zero) is $\min \{{T}_{j}^{L},{T}_{j}^{R}\}$ and is associated to the total survival probability

$$\begin{eqnarray}&&{\phi }_{{x}_{j}}(\tau )={{\boldsymbol{\varphi }}}_{j}^{L}(\tau ){{\boldsymbol{\varphi }}}_{j}^{R}(\tau ).\end{eqnarray} \tag{ 36 }$$

Once the transition time is known, either the left or right trigger is chosen, according to the age-dependent rates

$$\begin{eqnarray}&&{\beta }_{j}^{(\rho )}(\tau )=g(\tau ;{k}_{j}^{(\rho )},{\lambda }_{j}^{(\rho )}){\rm{\Gamma }}({k}_{j}^{(\rho )})/\gamma ({k}_{j}^{(\rho )},{\lambda }_{j}^{(\rho )}\tau ),\end{eqnarray} \tag{ 37 }$$

where $\gamma (k,x)$ is the lower incomplete Gamma function. The SCGF of the left current is computed by biasing the WTDs ${\psi }_{{x}_{1},{x}_{0},1}(\tau )$ and ${\psi }_{{x}_{0},{x}_{1},-1}(\tau )$ with ${{\rm{e}}}^{\mp s}$, respectively.

While the implementation of the method of section 3.2 remains straightforward for this model, a general solution for the exact SCGF is missing. We thus specialise to the case with ${k}_{0}^{R}={k}_{1}^{R}={k}_{0}^{L}={k}_{1}^{L}=2$, for which the Laplace transform of ${\psi }_{{x}_{i},{x}_{j}}(t)={\sum }_{c}{\psi }_{{x}_{i},{x}_{j},c}(t)$ is a rational function of ν, viz.

$$\begin{eqnarray}&&{\overline{\psi }}_{{x}_{i},{x}_{j}}(\nu )=({\alpha }_{j,2}^{R}+{\alpha }_{j,2}^{L}){\left(\displaystyle \frac{{\lambda }_{j}}{\nu +{\lambda }_{j}}\right)}^{2}+({\alpha }_{j,3}^{R}+{\alpha }_{j,3}^{L}){\left(\displaystyle \frac{{\lambda }_{j}}{\nu +{\lambda }_{j}}\right)}^{3},\end{eqnarray} \tag{ 38 }$$

where we defined for convenience ${\lambda }_{j}={\lambda }_{j}^{L}+{\lambda }_{j}^{R}$ and

$$\begin{eqnarray}&&{\alpha }_{j,2}^{L}=\displaystyle \frac{{({\lambda }_{j}^{L})}^{2}}{{({\lambda }_{j})}^{2}},\quad {\alpha }_{j,3}^{L}=\displaystyle \frac{2{({\lambda }_{j}^{L})}^{2}{\lambda }_{j}^{R}}{{({\lambda }_{j})}^{3}},\quad {\alpha }_{j,2}^{R}=\displaystyle \frac{{({\lambda }_{j}^{R})}^{2}}{{({\lambda }_{j})}^{2}},\quad {\alpha }_{j,3}^{R}=\displaystyle \frac{2{({\lambda }_{j}^{R})}^{2}{\lambda }_{j}^{L}}{{({\lambda }_{j})}^{3}}.\end{eqnarray} \tag{ 39 }$$

This can be though of as resulting from a walker spending exponentially distributed times in three hidden stages, as in [44]. Hence, the two-state semi-Markov process with WTD (32) and (33) can be seen as a six-state Markov process, see online supplementary material for more details and an alternative formulation. The linearity of the Laplace transform permits the distinction of the left and right contributions in (38), hence it remains easy to bias the rates that correspond to a change in J and the SCGF can be exactly found as the leading eigenvalue of a modified Markovian stochastic generator (see again supplementary material). Figure 2(b) shows convincing agreement of our method with this exact approach.

More general non-Markovian systems are those whose WTDs depend on events which occurred during the whole observation time. Systems in this class are the 'elephant' random walk [45] and its analogues, where the transition probabilities at time t depend on the history through the time-averaged current $j(t)$. We focus here on an interacting particle system with such current-dependent rates, namely the TASEP of [46].

Non-Markovian interacting particle systems can be described by assigning a trigger for jump attempts with WTD ${{\boldsymbol{\psi }}}_{i}[\tau ;w(t)]$ and a corresponding survival function ${{\boldsymbol{\varphi }}}_{i}[\tau ;w(t)]$ to each elementary event i that controls the particle dynamics. The probability density that the next transition is of type i and occurs in the time interval $[t+\tau ,t+\tau +{\rm{d}}t)$, given that, for each j, a time τ_j has elapsed since the last event of type j, is given by⁶

$$\begin{eqnarray}&&{\psi }_{i}[\tau ;w(t)]={{\boldsymbol{\psi }}}_{i}[{\tau }_{i}+\tau ;w(t)| {\tau }_{i}]\displaystyle \prod _{j\ne i,j=1}{{\boldsymbol{\varphi }}}_{j}[{\tau }_{j}+\tau ;w(t)| {\tau }_{j}],\end{eqnarray} \tag{ 40 }$$

where ${{\boldsymbol{\psi }}}_{i}[{\tau }_{i}+\tau ;w(t)| {\tau }_{i}]$ $=\,{{\boldsymbol{\psi }}}_{i}[{\tau }_{i}+\tau ;w(t)]/{{\boldsymbol{\varphi }}}_{i}[{\tau }_{i};w(t)]$ and ${{\boldsymbol{\varphi }}}_{i}[{\tau }_{i}+\tau ;w(t)| {\tau }_{i}]$ $=\,{{\boldsymbol{\varphi }}}_{i}[{\tau }_{i}+\tau ;w(t)]/{{\boldsymbol{\varphi }}}_{i}[{\tau }_{i};w(t)]$. With exact expressions for these WTDs, we can implement the algorithms of section 3.

The TASEP consists of a one-dimensional lattice of length L, where each lattice site l, $1\leqslant l\leqslant L$, can be either empty (${\eta }_{l}=0$) or occupied by a particle (${\eta }_{l}=1$). Particles on a site $l\lt L$ are driven rightwards: they attempt a bulk jump to site $l+1$ with WTD ${{\boldsymbol{\psi }}}_{b}[\tau ;w(t)]$, the attempt being successful if ${\eta }_{l+1}=0$, as in [47–49]. With open boundaries, a particle that reaches the rightmost site L leaves the system with WTD ${{\boldsymbol{\psi }}}_{L}[\tau ;w(t)]$. Also, as soon as ${\eta }_{1}=0$, a further boundary mechanism turns on and particles arrive on the leftmost site with WTD ${{\boldsymbol{\psi }}}_{0}[\tau ;w(t)]$. The special choice ${{\boldsymbol{\psi }}}_{0}[\tau ;w(t)]=\alpha {{\rm{e}}}^{-\alpha \tau }$, ${{\boldsymbol{\psi }}}_{b}[\tau ;w(t)]=p{{\rm{e}}}^{-p\tau }$, and ${{\boldsymbol{\psi }}}_{L}[\tau ;w(t)]=\beta {{\rm{e}}}^{-\beta \tau }$ corresponds to the standard Markovian TASEP with constant left, bulk and right rates α, p, and β.

We now assume that only the left boundary has a non-exponential WTD, while the particle triggers have exponential WTDs with rate 1 for free particles in the bulk, and rate β for the particle on the rightmost site. Consequently the inter-event time density distribution, conditioned on a time τ₀ having elapsed since the last arrival, is

$$\begin{eqnarray}\psi [\tau ;w(t)| {\tau }_{0}] & = & \left(\displaystyle \frac{{{\boldsymbol{\psi }}}_{0}[\tau +{\tau }_{0};w(t)]}{{{\boldsymbol{\varphi }}}_{0}[{\tau }_{0};w(t)]}(1-{\eta }_{0})+{\mathsf{n}}+\beta {\eta }_{L}\right)\\ & & \times \exp \left\{\mathrm{ln}\left[\displaystyle \frac{{\varphi }_{0}[\tau +{\tau }_{0};w(t)]}{{\varphi }_{0}[{\tau }_{0};w(t)]}\right](1-{\eta }_{0})-({\mathsf{n}}+\beta {\eta }_{L})\tau \right\}.\end{eqnarray} \tag{ 41 }$$

The probability mass distribution, conditioned on an age τ and elapsed time τ₀ is

$$\begin{eqnarray}{p}_{0}[\tau ;w(t)| {\tau }_{0}] & = & \displaystyle \frac{{\psi }_{0}[\tau +{\tau }_{0};w(t)]}{{\varphi }_{0}[{\tau }_{0};w(t)]}(1-{\eta }_{0})\\ & & \times {\left(\displaystyle \frac{{\psi }_{0}[\tau +{\tau }_{0};w(t)]}{{\varphi }_{0}[{\tau }_{0};w(t)]}(1-{\eta }_{0})+{\mathsf{n}}+\beta {\eta }_{L}\right)}^{-1},\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{p}_{i}[\tau ;w(t)| {\tau }_{0}]={\left(\displaystyle \frac{{\psi }_{0}[\tau +{\tau }_{0};w(t)]}{{\varphi }_{0}[{\tau }_{0};w(t)]}(1-{\eta }_{0})+{\mathsf{n}}+\beta {\eta }_{L}\right)}^{-1},\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{p}_{L}[\tau ;w(t)| {t}_{0}]=\beta {\delta }_{{\eta }_{L}}{\left(\displaystyle \frac{{\psi }_{0}[\tau +{t}_{0};w(t)]}{{\varphi }_{0}[{t}_{0};w(t)]}(1-{\eta }_{0})+{\mathsf{n}}+\beta {\eta }_{L}\right)}^{-1},\end{eqnarray} \tag{ 44 }$$

where η₀ and η_L encode the exclusion rules at boundaries, $i=1,2,\ldots ,{\mathsf{n}}$, and n is the number of free particles in the bulk, which depends on the lattice configuration before the jump.

Let us impose now that the arrival rate α depends linearly on the input current $j(t)$, i.e., $\alpha (j)={\alpha }_{0}+{aj}$, which defines a time-dependent rate ${\beta }_{0}(t):= \alpha [j(t)]$. Similar functional dependence (but on the instantaneous output current) has been used to model ribosome recycling in protein translation [50, 51]. Generically, such rates describe a simple form of positive feedback (for $a\gt 0$), whose effect on the stationary state of the TASEP is to shrink the low-density phase [46, 51]. The current fluctuations are also altered; the rate function $\hat{e}(j)$ in this phase has already been computed, for our model, by means of the temporal additivity principle [8, 46], hence this provides a testing ground for the cloning method of section 3. The particle arrival mechanism starts when the leftmost site is emptied, when we set an age of $\tau =0$. Denoting by q the current immediately after the last arrival, which occurred at $t-{\tau }_{0}$, the value $j(t+\tau )$ at age τ can be expressed as $q(t-{\tau }_{0})/(t+\tau )$, hence the trigger hazard is

$$\begin{eqnarray}&&{\beta }_{0}(t+\tau )={\alpha }_{0}+{aq}(t-{\tau }_{0})/(t+\tau ),\end{eqnarray} \tag{ 45 }$$

where τ is the trigger age. Initial values of τ₀ and $q$ are chosen to be 1 and 0, respectively. This allows us to derive the trigger survival probability and trigger WTD (see also [46]) which are, respectively

$$\begin{eqnarray}&&{{\boldsymbol{\varphi }}}_{0}[\tau ;w(t)]={\left(\displaystyle \frac{t}{t+\tau }\right)}^{{aq}(t-{t}_{0})}{{\rm{e}}}^{-{\alpha }_{0}\tau },\end{eqnarray} \tag{ 46 }$$

and ${{\boldsymbol{\psi }}}_{0}[\tau ;w(t)]={\beta }_{0}(t+\tau ){{\boldsymbol{\varphi }}}_{0}[\tau ;w(t)]$. Using these in equations (41) and (44) allows us to generate the trajectories, and the large deviation function can be evaluated by applying a bias ${{\rm{e}}}^{-s}$ to the arrival WTD and following the algorithm of section 3. The results are plotted in figure 3 and validated by the exact numerical calculation of [46] (which assumes the temporal additivity principle).

Zoom In Zoom Out Reset image size

It is worth noting that, for large negative values of s, the SCGF displays a linear branch with slope ${j}^{* }$. If $e(s)$ remains linear with the same slope for $s\lt -4$, its Legendre–Fenchel transform will be defined only for $j\leqslant \,{j}^{* }$. This appears to be related to the dynamical phase transition seen in the Markovian TASEP, where large current fluctuations require correlations on the scale of the system size and the rate function, in the corresponding regime, diverges with L [52]. Indeed, space–time diagrams (not shown) of the density profile from the cloning simulations seem to suggest that the correlation length increases as s becomes more negative. It would be interesting to further reduce the finite-ensemble errors (as discussed in section 6) in order to probe larger negative values of s.