Extreme-Value Statistics of Stochastic Transport Processes: Applications to Molecular Motors and Sports

TL;DR

This paper derives exact finite-time statistics for extrema in stochastic transport processes, applying the results to molecular motors and sports, revealing new symmetries and connections to random matrix theory.

Contribution

It introduces a novel analytical framework for the finite-time extrema statistics of biased random walks, connecting them to random matrix eigenvalue distributions and generalizing existing entropy production laws.

Findings

Exact formulas for extrema distributions in biased random walks.

Symmetry relations in the distribution of maxima and minima.

Application of random matrix theory to estimate extreme events.

Abstract

We derive exact expressions for the finite-time statistics of extrema (maximum and minimum) of the spatial displacement and the fluctuating entropy flow of biased random walks. Our approach captures key features of extreme events in molecular motor motion along linear filaments. Our results generalize the infimum law for entropy production and reveal a symmetry of the distribution of its maxima and minima. We also show that the relaxation spectrum of the full generating function, and hence of any moment, of the finite-time extrema distributions can be written in terms of the Mar{\v{c}}enko-Pastur distribution of random-matrix theory. Using this result, we obtain estimates for the extreme-value statistics of stochastic transport from the eigenvalue distributions of suitable Wishart and Laguerre random matrices. We confirm our results by numerical simulations of stochastic models of molecular motors and discuss as illustrative example our theory in the context of sports.

Tables

Table 1. Table 1: Number of goals scored and conceded by the teams qualified 1st at the end of the group stage in Women’s (10 teams out of 40) and Men’s (8 teams out of 32) UEFA Champions League in the 2018-2019 season, and their total number of matches. Estimation of the 1D model parameters from data of the season 2018-2019 of the soccer champions leagues (women and men). The scoring rates are given by the goals scored and goals conceded divided by the number of matches and the duration of a match ( τ match = 90 subscript 𝜏 match 90 \tau_{\rm match}=90 min), i.e. k + = GF / ( τ match ​ MP ) subscript 𝑘 GF subscript 𝜏 match MP k_{+}=\text{GF}/(\tau_{\rm match}\text{MP}) and k − = GC / ( τ match ​ MP ) subscript 𝑘 GC subscript 𝜏 match MP k_{-}=\text{GC}/(\tau_{\rm match}\text{MP}) . The rest of the parameters are obtained from our expressions corresponding to the 1D biased random walk: the scoring affinity A s subscript 𝐴 s A_{\rm s} and rate ν s subscript 𝜈 s \nu_{\rm s} [Eq. ( 2 )], and the extreme-value timescales τ 0 subscript 𝜏 0 \tau_{0} , τ 1 subscript 𝜏 1 \tau_{1} , τ ¯ ¯ 𝜏 \bar{\tau} and τ ∞ subscript 𝜏 \tau_{\infty} , see Sec.  III . The uncertainties are obtained by propagating the standard error of the mean Poissonian numbers of goals per match.

Soccer teams	Women	Men
Goals scored GF	113	108
Goals conceded GC	16	42
Number of matches MP	30	48
$k_{+}$ (goals/min)	$0.042\pm 0.006$	$0.025\pm 0.010$
$k_{-}$ (goals/min)	$0.006\pm 0.002$	$0.009\pm 0.002$
$A_s=\ln \left(k_{+}/k_{-}\right)$	$1.95\pm 0.27$	$0.94\pm 0.18$
${\nu }_s=\sqrt{k_{+}{k}_{-}}$ (goals/min)	$0.016\pm 0.002$	$0.016\pm 0.001$
${\tau }_0={\left(\sqrt{k_{+}}+\sqrt{k_{-}}\right)}^{-2}$ (min)	$13\pm 1$	$15\pm 1$
${\tau }_1={\left(k_{+}+k_{-}\right)}^{-1}$ (min)	$21\pm 2$	$29\pm 2$
$\overline{\tau }=k_{+}{\left(k_{+}-k_{-}\right)}^{-2}$ (min)	$32\pm 5$	$110\pm 30$
${\tau }_{\infty }={\left(\sqrt{k_{+}}-\sqrt{k_{-}}\right)}^{-2}$ (min)	$61\pm 13$	$280\pm 100$

Equations

$$k(x,y)=\nu (x,y)e^{A(x,y)/2},$$

$$\nu =\sqrt{k_{+}{k}_{-}},A=\ln (k_{+}/k_{-}).$$

$$P(X_{\min }=-x)$$

$$⟨X_{\min }⟩=\frac{-1}{e^A-1},⟨S_{\min }⟩=\frac{-A}{e^A-1}.$$

$$P(X_{\min }(t)=-x)$$

$$P_{\mathrm{fpt}}(\mathsf{T};x)=e^{Ax/2}\frac{|x|}{\mathsf{T}}{I}_x(2\nu \mathsf{T})e^{-2\nu \cosh (A/2)\mathsf{T}},$$

$$P_{\mathrm{fpt}}(\mathsf{T};x)/P_{\mathrm{fpt}}(\mathsf{T};-x)=e^{Ax},$$

$$P_{\mathrm{fpt}}(\mathsf{T};x)/{\tilde{P}}_{\mathrm{fpt}}(\mathsf{T};x)=e^{Ax}.$$

$$⟨S_{\min }(t)⟩=-\frac{A}{2\pi }{\int }_{-1}^1\frac{1-e^{-2\nu t(\cosh (A/2)-y)}}{{\left(\cosh (A/2)-y\right)}^2}\sqrt{1-y^2}\mathrm{d}y.$$

$$⟨S_{\max }(t)⟩-⟨S(t)⟩=-⟨S_{\min }(t)⟩.$$

$$P\left(S(0)-S_{\min }(t)=s\right)=P\left(S_{\max }(t)-S(t)=s\right),$$

$$⟨S_{\min }(t)⟩=⟨S_{\min }⟩\left(1-{\int }_{{\tau }_0}^{{\tau }_{\infty }}{e}^{-t/\tau }\rho (\tau /\bar{\tau })\frac{\mathrm{d}\tau }{\bar{\tau }}\right).$$

$$\rho (\lambda )\equiv \{\begin{array}{ll}\frac{1}{2\pi \delta }\frac{\sqrt{({\lambda }_{+}-\lambda )(\lambda -{\lambda }_{-})}}{\lambda }& \text{mbox}if\lambda \in [{\lambda }_{-},{\lambda }_{+}]\\ 0& \text{mbox}if\lambda \notin [{\lambda }_{-},{\lambda }_{+}],\end{array}$$

$$G_{\min }(z;t)=1+\frac{1-z}{e^A-1}{\int }_{{\tau }_0}^{{\tau }_{\infty }}\frac{1-e^{-t/\tau }}{1+f(z)\tau }\rho (\tau /\bar{\tau })\frac{\mathrm{d}\tau }{\bar{\tau }},$$

$$⟨S_{\min }(t)⟩=⟨S_{\min }⟩\left(1-\underset{m\to \infty }{\lim }\frac{1}{m}\mathrm{Tr}{e}^{-{\mathbf{M}}^{-1}t/\bar{\tau }}\right),$$

$$⟨S_{\min }(t)⟩\simeq ⟨S_{\min }⟩\left(1-\frac{1}{m}\sum _{i=1}^m{e}^{-t/(\bar{\tau }{\lambda }_i)}\right),$$

$${\nu }_{\mathrm{eff}}$$

$$A_{\mathrm{eff}}$$

$$k_{\mathrm{x}}^{\pm }$$

$$f_{\mathrm{stall}}=\frac{1}{\beta \ell }\ln \left(\frac{1+k_{\mathrm{mc}}^{-}/k_{\mathrm{m}}^{-}+k_{\mathrm{cm}}^{+}/k_{\mathrm{m}}^{-}}{1+k_{\mathrm{mc}}^{+}/k_{\mathrm{m}}^{+}+k_{\mathrm{cm}}^{-}/k_{\mathrm{m}}^{+}}\right).$$

$$⟨S_{\max }(t)⟩-⟨S(t)⟩\le 1.$$

$$\frac{\mathrm{d}{P}_x(t)}{\mathrm{d}t}=k_{+}{P}_{x-1}(t)-(k_{+}+k_{-})P_x(t)+k_{-}{P}_{x+1}(t),$$

$$P_x(t)={\left(\frac{k_{+}}{k_{-}}\right)}^{x/2}{I}_x\left(2\sqrt{k_{+}{k}_{-}}t\right)e^{-(k_{+}+k_{-})t},$$

$$P_x(t)\sim e^{-2\nu tJ(x/vt)}.$$

$$I_x(x/z)$$

$$\sim$$

$$J(u)$$

$$J(u)=\frac{\sinh (A/2{)}^2(u-1{)}^2}{2\cosh (A/2)}+\frac{\sinh (A/2{)}^4(u-1{)}^3}{6\cosh (A/2{)}^3}+O(u-1{)}^4.$$

$$⟨e^{-S(t)}|X_{[0,t^{\prime }]}⟩=e^{-S(t^{\prime })},$$

$$⟨e^{-S(t)}|X_{[0,t^{\prime }]}⟩$$

$$⟨e^{-S(\mathsf{T})}⟩=1,$$

Full text

I Introduction

Life is a non-equilibrium phenomenon characterised by fluxes of energy and matter at different scales. At the molecular level, molecular motors play a key role for the generation of movements and forces in cells. Examples are vesicle transport, muscle contraction, cell division and cell locomotion Jülicher et al. (1997); Howard et al. (2001). A molecular motor consumes a chemical fuel, adenosine triphosphate (ATP), that is hydrolysed to adenosine diphosphate (ADP) and inorganic phosphate. The chemical energy of this reaction is transduced to generate spontaneous movements and mechanical work. Single-molecule experiments have revealed that the activity of single or a few molecular motors displays strong fluctuations Noji et al. (1997); Bustamante et al. (2000); Schnitzer et al. (2000); Ishii and Yanagida (2000); Helenius et al. (2006); Ritort (2006); Vogel et al. (2009); Aathavan et al. (2009); Ananthanarayanan et al. (2013); Fakhri et al. (2014); Naganathan et al. (2016) which can be captured by the theory of stochastic processes Ndlec et al. (1997); Keller and Bustamante (2000); Qian (2005); Klumpp and Lipowsky (2005); Campas et al. (2006); Kolomeisky and Fisher (2007); Brangwynne et al. (2008); Gaspard (2016).

An important question is to understand general features and universal properties that govern the statistics of fluctuations of stochastic transport processes that include the motion of molecular motors. Universal relations for the fixed-time statistics of time-integrated currents, such as the distance traveled and the work performed, have been investigated in the framework of non-equilibrium stochastic thermodynamics Sekimoto (1997); Lau et al. (2007); Chemla et al. (2008); Toyabe et al. (2010). These results provide e.g. universal bounds for the efficiency of molecular motors given by the ratio between the mechanical power and the chemical power put in the motor Pietzonka et al. (2016). Timing statistics of enzymatic reactions, such as those powering the motion of molecular motors, have been discussed within the framework of Kramers theory Hänggi et al. (1990). Recent theory and experiments in Kinesin have revealed symmetry relations between forward and backward cycle-time distributions of enzymatic reactions Qian and Xie (2006); Nishiyama et al. (2002); Carter and Cross (2005). Related results have been derived in the context of waiting times of active molecular processes Neri et al. (2017) and transition-path times in folding transitions of DNA hairpins Gladrow et al. (2019).

When discussing dynamic processes, it is also sufficient to study averages and small fluctuations. However, rare events and large fluctuations play an important role when resilience and reliability of a system are investigated. In this context, the statistics of extreme values and of extreme excursions from the average play an important role, as has been discussed in fields ranging from statistical mechanics to climate Cheng et al. (2014); Diffenbaugh et al. (2005) and finance Rocco (2014); Brodin and Klüppelberg (2014); Novak (2011). Extreme events are also important in biophysics and are key to understand the robustness of biological processes. Illustrating examples are microtubule catastrophes or a sperm winning a race against a billion competitors. Here we consider the extreme value statistics of transport processes. In spite of the significant progress in research on steady-state currents and time fluctuations, little is known yet about extreme-value statistics of stochastic processes such as active biomolecular processes. Extreme-value theory has provided useful insights for e.g. long-range correlations of DNA sequences Peng et al. (1992); Arneodo et al. (1995) and DNA replication statistics in frogs’ embryonic cells Bechhoefer and Marshall (2007). However, extreme-value statistics of molecular motors have not been discussed so far. For example, what is the maximal excursion of a stochastic motor against or in the average direction of its motion within a given time? How long does it take a motor to reach its maximum excursion against the chemical bias? What is the entropy production associated with an extreme fluctuation of a molecular motor?

In this article, we provide novel insights on the aforementioned questions by deriving exact results for extreme-value statistics of simple models of stochastic transport given in Eqs. (3-4) and (9-15). We discuss the statistics of the maximum and minimum excursion (with respect to its initial location) and the associated extremal entropy changes. Moreover, we investigate the timescales associated with those extrema, combining concepts from stochastic thermodynamics, random walks and random-matrix theory. As we show below, our results provide insights on extreme-value statistics beyond recently derived inequalities for the finite-time infimum statistics of entropy production Chetrite and Gupta (2011); Neri et al. (2017), and relate to record statistics of correlated stochastic processes Krug (2007); Majumdar and Ziff (2008); Wergen et al. (2011); Godreche et al. (2014); Harris (2015); Bénichou et al. (2016); Hartich and Godec (2019); Sabhapandit (2019).

The article is organized as follows: Sec. II describes the stochastic model of stochastic transport used in this paper and provides exact extreme-value statistics for one dimensional (1D) biased random walks. Sec. III discusses the connection between extrema of 1D biased random walks and random matrix theory. In Sec. IV we apply our theory to two-dimensional stochastic models of molecular motors and discuss its implications using sports as illustrative example. Sec. VI concludes the paper. Details on the mathematical derivations and numerical simulations are provided in the Appendices.

II Extreme-value statistics of 1D biased random walks

Many nonequilibrium phenomena at mesoscales can be described at a coarse-grained level as a continuous-time Markov jump process between discrete states $x,y,z$ etc, with exponential waiting times. The transition rate from $x$ to $y$ can be written as Van Kampen (1992); Bressloff (2014); Baiesi and Maes (2018)

[TABLE]

with $\nu(x,y)=\nu(y,x)$ symmetric and $A(x,y)=-A(y,x)$ antisymmetric with respect to the exchange $x\to y$. If local detailed balance holds $A(x,y)=\beta[W(x)-W(y)]$ with $W(x)$ the potential energy of state $x$, $\beta=(k_{\rm B}T)^{-1}$ with $k_{\rm B}$ Boltzmann’s constant and $T$ the temperature of a thermostat. First, we consider the simple case of a 1D biased random walk on a line with discrete states denoted by the integer $x\in\mathbb{Z}$ (see Fig. 1a for an illustration). The forward and backward jump rates are given by $k_{\pm}\equiv k(x,x\pm 1)=\nu e^{\pm A/2}$. Here $\nu>0$ is a rate and $A>0$ the affinity, which satisfy

[TABLE]

This biased random walk describes e.g. the motion of a molecular motor along a periodic track fuelled by ATP Kolomeisky and Fisher (2007). In the simplest case, $A=\beta\Delta\mu$, with $\Delta\mu=\mu_{\rm ATP}-(\mu_{\rm ADP}+\mu_{P})$ the chemical potential difference of ATP hydrolysis, often of the order of $20k_{\rm B}T$, and the rate $\nu$ depends on ATP concentration and internal timescales that determine the dwell-time statistics of the motor.

An individual trajectory of a motor starting from a reference state $X(0)=0$ at time $t=0$ is denoted by $X_{[0,t]}=\{X(s)\}_{s=0}^{t}$. It contains jumps $j=1,2,\dots$ from state $x_{j}^{-}$ to state $x_{j}^{+}$ that occur at random times $t_{j}$. The entropy production in units of $k_{\rm B}$ associated with this trajectory is $S(t)=\ln[\mathcal{P}(X_{[0,t]})/\mathcal{P}(\tilde{X}_{[0,t]})]=AX(t)$ Neri et al. (2017). Here $\mathcal{P}$ is the path probability and $\tilde{X}_{[0,t]}=\{X(t-s)\}_{s=0}^{t}$ is the time reversed path. Thus, the entropy production $S(t)$ is a stochastic variable that undergoes a biased random walk of step size $A$ with trajectories $S_{[0,t]}=AX_{[0,t]}$. For $A$ positive, both the average velocity $v=\langle X(t)\rangle/t=(k_{+}-k_{-})=2\nu\sinh(A/2)$ and the average rate of entropy production $\sigma=\langle S(t)\rangle/t=vA$ are positive. Here and in the following we denote by $\langle\,\cdot\,\rangle$ averages over many realizations of the process $X(t)$. However, due to fluctuations, the stochastic variables $X(t)$ and $S(t)$ can in principle take any value with finite probability and even become negative.

We now derive exact expressions for the statistics of the minimum $X_{\rm min}(t)=\text{min}_{\tau\in[0,t]}X(\tau)$ and the maximum $X_{\rm max}(t)=\text{max}_{\tau\in[0,t]}X(\tau)$ of the position of the motor with respect to its initial position, see Fig. 1b for illustrations. We also discuss the global minimum and maximum of the stochastic entropy production $S(t)=AX(t)$ denoted by $S_{\rm min}(t)$ and $S_{\rm max}(t)$, respectively. We first discuss the statistics of the global extrema of the position $X_{\rm min}\equiv\lim_{t\to\infty}X_{\rm min}(t)$ and $X_{\rm max}\equiv\lim_{t\to\infty}X_{\rm max}(t)$, and of the entropy production, $S_{\rm min}$ and $S_{\rm max}$. The probability that the global minimum of the discrete position is $-x$, with $x\geq 0$, is $P(X_{\rm min}=-x)=P_{\rm abs}(-x)-P_{\rm abs}(-x-1)$, where $P_{\rm abs}(-x)=e^{-Ax}$ Neri et al. (2017) is the probability that $X(t)$ reaches an absorbing site in $-x$ at a finite time. Thus, the global minimum follows a geometric distribution

[TABLE]

for $x\geq 0$ and $P(X_{\rm min}=x)=P(S_{\rm min}=Ax)=0$ for $x>0$. From Eq. (3) we obtain the mean global minimum of a 1D biased random walk and of its associated entropy production:

[TABLE]

Therefore, the global minimum of the position diverges in the limit of small bias $A$ whereas the entropy production minimum is bounded for all $A\geq 0$ and obeys the infimum law $\langle S_{\rm min}\rangle\geq-1$ Neri et al. (2017). This bound is saturated in the limit of small affinity, which corresponds to the diffusion limit Pigolotti et al. (2017). Because $S(t)$ and $X(t)$ have positive drift, the average global maxima of entropy production and displacement are not defined. However the difference $\lim_{t\to\infty}[\langle S_{\rm max}(t)\rangle-\langle S(t)\rangle]=A/(e^{A}-1)$ is finite and obeys symmetry properties that we discuss below.

Finite-time extrema statistics of the 1D biased random walk may be obtained from the finite-time absorption probabilities

[TABLE]

where $P_{\rm abs}(-x;t)$ is the probability that $X$ reaches an absorbing site at $-x$ at any time smaller or equal than $t$. The absorption probability $P_{\rm abs}(x;t)=\delta_{x,0}+\int_{0}^{t}P_{\rm fpt}(\mathsf{T};x){\rm d}\mathsf{T}$, with $\delta_{i,j}$ Kronecker’s delta and

[TABLE]

is the first-passage time probability for the motor to first reach an absorbing site in $x$, with $|x|\geq 1$, at time $\mathsf{T}\geq 0$ Redner (2001); Roldán et al. (2016), see Appendix A. Here $I_{x}$ denotes the $x-$th order modified Bessel function of the first kind. Note that $\int_{0}^{\infty}\text{d}\mathsf{T}P_{\rm fpt}(\mathsf{T};x)=P_{\rm abs}(x)\leq 1$. We identify in Eq. (6) two timescales. The smaller timescale $\tau_{1}=(k_{+}+k_{-})^{-1}=(2\nu\cosh(A/2))^{-1}$ is the average waiting time between two jumps, and $\tau_{2}=(2\sqrt{k_{+}k_{-}})^{-1}=(2\nu)^{-1}$ is inversely proportional to the geometric mean of the transition rates; their ratio $\tau_{2}/\tau_{1}=\cosh(A/2)\geq 1$ increases with the bias strength. Normalizing (6) by $P_{\rm abs}(x)$, we obtain the mean $\langle\mathsf{T}\rangle=|x|A/\sigma$ and variance $\text{Var}[\mathsf{T}]=(\coth(A/2)/|x|)\langle\mathsf{T}\rangle^{2}$ of the first-passage time, in agreement with the first-passage time uncertainty relation $\text{Var}[\mathsf{T}]/\langle\mathsf{T}\rangle\geq 2/\sigma$ Gingrich and Horowitz (2017). Furthermore, the first-passage time probability density (6) obeys the following symmetry properties. First, the ratio

[TABLE]

is independent on $\mathsf{T}$, as follows from the stopping-time fluctuation theorem Neri et al. (2017); Qian and Xie (2006); Singh et al. (2019). Second, the “conjugate” first-passage time probability $\tilde{P}_{\rm fpt}(\mathsf{T};x)$, obtained exchanging $k_{+}$ by $k_{-}$ (i.e. $A$ by $-A$), obeys

[TABLE]

These two properties imply $\tilde{P}_{\rm fpt}(\mathsf{T};x)=P_{\rm fpt}(\mathsf{T};-x)$, which has interesting consequences for random walks Krapivsky and Redner (2018) and for the extrema statistics of $S(t)$, see below.

In order to derive exact finite-time extrema statistics, it is often convenient to use generating functions and Laplace transforms. The generating functions of the distributions of finite-time entropy production extrema are defined as $G_{\rm min/max}(z;t)=\sum_{x=-\infty}^{\infty}z^{x}P\left(S_{\rm min/max}(t)=Ax\right)\Theta(\mp x)$, with $\Theta$ the Heaviside function. Their Laplace transforms are given by $\hat{G}_{\rm min/max}(z;s)=\int_{0}^{\infty}\text{d}te^{-st}G_{\rm min/max}(z;t)=s^{-1}[1-\hat{P}_{\rm fpt}(s;-1)]/[1-z^{\mp 1}\hat{P}_{\rm fpt}(s;-1)]$. Here, $\hat{P}_{\rm fpt}(s;n)=e^{Ax/2}e^{-|x|\,\mathrm{cosh}^{-1}\left[s/2\nu+\cosh(A/2)\right]}-\delta_{x,0}$ is the Laplace transform of the first-passage-time probability density at site $x$. These expressions enable the computation of the Laplace transform of all the moments of the extrema from successive derivatives of the generating functions with respect to $\ln z$. In particular, the Laplace transform of the average minimum of entropy production reads $s\langle\hat{S}_{\min}(s)\rangle=-A/[\hat{P}_{\rm fpt}(s;-1)^{-1}-1]$. In the time domain, we may write this equality as $\langle S_{\min}(t)\rangle=-A\int_{0}^{t}\text{d}\mathsf{T}\sum_{x=1}^{\infty}P_{\rm fpt}(\mathsf{T};-x)$ Hartich and Godec (2019) which can be written as (see Appendix B):

[TABLE]

Numerical simulations of the 1D biased random walk are in excellent agreement with Eq. (9) (Fig. 2, blue symbols). Note that Eq. (9) can also be expressed in terms of the Kampé de Fériet function F as $\langle S_{\min}(t)\rangle/A=-k_{-}t+{}^{2+0}F_{1+1}\left[\left.\genfrac{}{}{0.0pt}{}{[2\;1]\;,\;\varnothing\;,\;\varnothing}{3\;,\;2\;,\;2}\right\rvert-k_{-}t\;,\;-k_{+}t\right](k_{-}k_{+}t^{2})/2$ (see Appendix E). Interestingly, our simulations reveal (Fig. 2, red symbols) that the average maximum of entropy production minus the average entropy production at time $t$ equals to minus the right-hand side of Eq. (9):

[TABLE]

Equation (10) follows from the symmetry relation of the 1D biased random walk $\tilde{P}_{\rm fpt}(\mathsf{T};x)=P_{\rm fpt}(\mathsf{T};-x)$ and the relation $\langle S_{\max}(t)\rangle-\langle S(t)\rangle=A\int_{0}^{t}\text{d}\mathsf{T}\sum_{x=1}^{\infty}\tilde{P}_{\rm fpt}(\mathsf{T};x)$. Moreover, this symmetry extends to the distribution of minima and maxima of entropy production

[TABLE]

where, in this case, $S(0)=0$. Figure 3 shows empirical distributions of entropy-production minima and maxima obtained from numerical simulations, which fulfil the symmetry relation (11).

III Random-matrix approach

We now explore a connection between entropy-production extrema and random-matrix theory. More precisely, we relate the previously derived expressions for the average and distribution of extrema with eigenvalue distribution of specific random matrices. Equation (9) can also be written as (see Appendix C)

[TABLE]

Here $\tau_{0}\equiv\left(\sqrt{k_{+}}+\sqrt{k_{-}}\right)^{-2}$ is the minimal relaxation time of the extreme value statistics and $\tau_{\infty}\equiv\left(\sqrt{k_{+}}-\sqrt{k_{-}}\right)^{-2}$ is the maximal extrema relaxation time. Here $\rho$ is the Marčenko-Pastur distribution, where times are normalized by $\bar{\tau}\equiv k_{+}/(k_{+}-k_{-})^{2}$. The Marčenko-Pastur distribution is given by Marčenko and Pastur (1967):

[TABLE]

with $\int_{0}^{\infty}\rho(\lambda){\rm d}\lambda=1$, where $\delta=k_{-}/k_{+}=e^{-A}$ and $\lambda_{\pm}=\left(1\pm\sqrt{\delta}\right)^{2}=(\sqrt{k_{+}}\pm\sqrt{k_{-}})^{2}/k_{+}$. Interestingly, $\rho(\lambda)$ is the distribution of eigenvalues in the large size limit of Hermitian matrices drawn from the ensemble of the Wishart-Laguerre random matrices Nica and Speicher (2006); Fridman et al. (2012); Livan et al. (2018), whose structure is explained below. The average (12) follows from the generating function of the distribution of the minimum (see Appendix C)

[TABLE]

where $f(z)=k_{+}(z-1)+k_{-}(z^{-1}-1)$.

Equations (12-14) imply that the time at which the distribution of the extrema relax to their long time limit is given by the largest timescale of the relaxation spectrum $\tau_{\infty}=\bar{\tau}\lambda_{+}$. We demonstrate this result in an example shown in the inset of Fig. 2 which shows this for the case of $\langle S_{\rm min}(t)\rangle$. Furthermore, it implies that the following trace formula holds:

[TABLE]

where $\mathbf{M}$ is a $m\times m$ random matrix. Eq. (15) holds for $\mathbf{M}$ drawn from one of two random matrix ensembles: (i) Wishart matrices $\mathbf{M}=n^{-1}\mathbf{R}\mathbf{R}^{T}$, where $\mathbf{R}$ is a $m\times n$ rectangular random matrix, $\mathbf{R}^{T}$ its transpose where $n=\lceil e^{A}m\rceil$. Its entries $R_{ij}$ are independent and identically distributed Gaussian random numbers with zero mean and unit variance; (ii) Laguerre matrices $\mathbf{M}=n^{-1}\mathbf{R}\mathbf{R}^{T}$, where $\mathbf{R}$ is a $m\times m$ square matrix with $R_{ij}$ independent but not identically distributed random variables, drawn from specific $\chi$ distributions Dumitriu and Edelman (2002), see Appendix D.

Eq. (15) can be approximated numerically using random matrices with finite but sufficiently large $m$. In practice, we use the following estimate:

[TABLE]

where $\lambda_{i}$ is the $i$-th eigenvalue of $\mathbf{M}$ drawn from either the Wishart or Laguerre ensembles.

In Fig. 4 we illustrate the validity of our random-matrix approach. More precisely, we compare the theoretical value of the average minimum of entropy production with the estimate given by the right-hand side of Eq. (15). For this purpose, we compute numerically the eigenvalues of a single random matrix of the Wishart ensemble and of a matrix drawn from the $\beta-$Laguerre ensemble with parameter $\beta=2$. Notably, using a single $64\times 64$ random matrix from the Wishart or Laguerre ensembles, we obtain an estimate of the average entropy production minimum that differs with respect to the exact value only by about $2\%$ at all times.

IV Extrema statistics of molecular motors

We now investigate whether similar results also hold for more complex stochastic models of molecular motors. We consider a biochemical process where a molecular motor’s fluctuating motion is described by a continuous-time Markov jump process on a potential energy surface in two dimensions $x$ and $y$ (Fig. 5a). Here $x$ denotes the spatial displacement of the motor along a discrete track of period $\ell$, and $y$ is a chemical reaction coordinate denoting the net number of fuel molecules spent by the motor.

The motion of the motor is biased along the track by a mechanical force $f_{\rm ext}$ applied to the motor. In addition, the motor hydrolyzes ATP with chemical potential difference $\Delta\mu$. We consider that both $f_{\rm ext}$ and $\Delta\mu$ to be independent of the state of the motor, which corresponds to the limits where the external force and the concentration of fuel molecules are stationary. States $(x,y)$ of the motor are in local equilibrium at temperature $T=\beta^{-1}$. The dynamics of the motor is as follows. From a given state the motor can perform, at a random time, a jump to eight adjacent states corresponding to the following four transitions and their reversals: (i) sliding along the track by a distance $\ell$ ($-\ell$) without consuming fuel but generating work in (against) the direction of the force, at a rate $k_{\rm m}^{+}$ ($k_{\rm m}^{-}$), with $k_{\rm m}^{-}=k_{\rm m}^{+}e^{-\beta f_{\rm ext}\ell}$; (ii) consumption of one ATP (ADP) molecule without generating work, at a rate $k_{\rm c}^{+}$ ($k_{\rm c}^{-}$), with $k_{\rm c}^{-}=k_{\rm c}^{+}e^{-\beta\Delta\mu}$; (iii) work generation in (against) the external force using ATP (ADP) at a rate $k_{\rm mc}^{+}$ ($k_{\rm mc}^{-}$), with $k_{\rm mc}^{-}=k_{\rm mc}^{+}e^{-\beta(f_{\rm ext}\ell+\Delta\mu)}$ and (iv) work generation against (in) the external force using ATP (ADP) at a rate $k_{\rm cm}^{+}$ ($k_{\rm cm}^{-}$), with $k_{\rm cm}^{-}=k_{\rm cm}^{+}e^{\beta(f_{\rm ext}\ell-\Delta\mu)}$. We use transition rates of the form $k_{\alpha}^{\pm}=\nu_{\alpha}e^{\pm A_{\alpha}/2}$, where $\nu_{\alpha}$ give different weights to each transition type.

A single trajectory of the motor is a 2D random walk containing snapshots $(X(t),Y(t))$ of the state of the motor at time $t$. Here $X(t)$ is the spatial coordinate of the motor (with respect to its initial position) and $Y(t)$ is the reaction coordinate representing the net number of ATP molecules consumed up to time $t$. Note that when $Y(t)$ is negative, the motor has consumed more ADP than ATP molecules. The entropy production associated with a single trajectory of the molecular motor is $S(t)=A_{\rm m}X(t)+A_{\rm c}Y(t)$, where $A_{\rm m}=\beta f_{\rm ext}\ell$, $A_{\rm c}=\beta\Delta\mu$ are the mechanical and chemical affinities. Thus $S(t)$ is a random walk with four different step lengths $A_{\rm m}$, $A_{\rm c}$, $A_{\rm mc}\equiv A_{\rm m}+A_{\rm c}$ and $A_{\rm cm}\equiv-A_{\rm m}+A_{\rm c}$ corresponding to the jumps along the $X$, $Y$ and the diagonal directions, respectively.

We perform numerical simulations of this 2D stochastic model of the molecular motor using Gillespie’s algorithm, and evaluate the entropy flow associated with different trajectories of the motor. Obtaining exact extreme-value statistics in this model is challenging. However the following simple approximation provides good estimates. The finite-time average (Fig. 5b) and the distribution (Fig. 5c) of the entropy production extrema obtained from simulations can be approximated by Eq. (9) and Eq. (3) replacing $A$ and $\nu$ by the effective parameters

[TABLE]

where the index $\alpha$ runs over the four types of transitions $\alpha=\text{m}$, c, mc and cm. This approximation based on effective parameters follows from considering effective 1D models with jumping rate $k^{+}_{\rm eff}+k^{-}_{\rm eff}=\sum_{\alpha}k^{+}_{\alpha}+k^{-}_{\alpha}$. Moreover, the symmetry between the time-dependent distribution of the extrema (11) is yet satisfied with high accuracy in our numerical simulations, even though these distributions can have very irregular shapes depending on the ratios of the affinities (Fig. 5c).

The average extrema of the mechanical and chemical currents are shown in Fig. 6 as a function of the external force. The behavior of these currents can be understood from a mapping to 1D biased random walks $X(t)$ and $Y(t)$ with effective forward and backward hopping rates given, respectively, by

[TABLE]

This yielding effective affinities $A_{\rm x/y}$ and rates $\nu_{\rm x/y}$ defined as in Eq. (2). This allows us to calculate the extreme value statistics of the numbers of step $X(t)$ and $Y(t)$, as illustrated in Fig. 6 (a) as solid lines. Note that $A_{\rm x}=\beta(f_{\rm ext}-f_{\rm stall})\ell$ is related to the mechanical affinity $A_{\rm m}=\beta f_{\rm ext}\ell$ and the stall force

[TABLE]

In the examples shown in Figs. 5,6, $f_{\rm stall}\simeq-1.2$pN in the simulations. $A_{\rm y}$ and the chemical affinity $A_{\rm c}$ obey a similar relationship. We observe in Fig. 6 that the largest average extreme values of entropy production occur when the external force applied to the motor is near the stall force. Interestingly, the minima and maxima of the displacements, $X_{\min}$ and $X_{\max}$, show diverging averages when the stall force is approached from above or below respectively. Such behavior could be observable in experiments on the statistics of the stepping of molecular motors near stall forces. Note that in Fig. 6 the long time limit of the average extrema is already reached at $0.1$s, see Fig. 5.

V Illustration of extreme value statistics using data from sports

Our analysis of extreme values in Sec. II applies to any stochastic process described by a biased random walk, e.g. scores during a sports match. Let us consider a game where two teams (Team 1 and Team 2) play against each other. We denote the number of points scored by each team up to time $t$ by $N_{1}(t)$ and $N_{2}(t)$, respectively. We can define the average rate of scores for Team 1 and Team 2 as $k_{+}$ and $k_{-}$, where we consider without loss of generality $k_{+}>k_{-}$. Here $k_{\pm}$ are the expected number of points scored per unit of time (e.g. goals per minute). Their difference $X(t)=N_{1}(t)-N_{2}(t)$ is the net score in favor of Team 1, i.e. it is positive when Team 1 is winning, zero when the match is tied, and negative when Team 2 is winning. When applying our minimal model we assume for simplicity that the score of each team follows Poissonian statistics, i.e. we do not take into account correlations and non-stationary effects Karlis and Ntzoufras (2003).

We now discuss our results in the context of this simplified model of sports. At finite time $t$, the stronger team is expected to lead by $\langle X(t)\rangle=(k_{+}-k_{-})t$, with variance $\text{Var}[X(t)]=(k_{+}+k_{-})t$. The following questions then appear naturally: “By how many points do we expect the weaker team to lead over the stronger one during the game?”, “How long should a match last in order to see an exciting game with comebacks of the weaker team?”. We provide insights on these questions in terms of $X_{\min}(t)$ i.e. the maximum net score against the stronger team up to time $t$. Knowledge on the characteristic time of the score’s extrema statistics is useful for the design of an exciting game that takes our breath away as long as possible: the longer the game, the more certain the stronger team is to win but the more predictable is the game, whereas the shorter the game, the more random it is. Our theory reveals a continuous spectrum of the characteristic relaxation times of the score extrema, following the Marčenko-Pastur distribution (12-13). An exciting game with comebacks should last at least $\tau_{0}=\left(\sqrt{k_{+}}+\sqrt{k_{-}}\right)^{-2}$. This is because $\tau_{0}$ is the timescale for the first action in the game. It should not be longer than $\tau_{\infty}=\left(\sqrt{k_{+}}-\sqrt{k_{-}}\right)^{-2}$ because after that time the weaker team is not expected to lead anymore. Indeed, $\tau_{\infty}$ corresponds to the time after which the weaker time falls behind for the rest of the match. Interestingly the smallest timescale $\tau_{0}$ is shorter than the expected time for the first point $\tau_{1}=(k_{+}+k_{-})^{-1}$, which defines the maximum of the Marčenko-Pastur density function.

We illustrate these results using data from soccer, corresponding to the 2018-2019 season of Men’s and Women’s UEFA Champions League, see Table 1. We collect data of the total number of goals scored (GF) and goals conceded (GC) by all teams that qualified 1st at the end of the group stage. We estimate the parameters $k_{+}$ and $k_{-}$ for a typical match between the group leader against an average team, see Appendix.

We find that the effective scoring affinity $A_{\rm s}=\ln(k_{+}/k_{-})$, which measures the strength difference between group leaders and average soccer teams, is twice larger for Women’s than for Men’s soccer. This striking difference is not likely to be a fundamental difference in the way soccer is played for each gender, as suggested by their equal scoring rate $\nu_{\rm s}=0.016$goals/min. It rather indicates more heterogeneity between strong and average teams within the Women’s Champions League as compared to Men’s. This difference however changes drastically the relaxation timescales of the extreme statistics, as revealed in Fig. 8. The distribution of extreme-value timescales peaks close to the first half of the match both for Women’s ($21$min) and Men’s ($29$min) soccer. On the other hand, their largest timescale for comebacks differ strongly, $\tau_{\infty}\simeq 61$min for Women’s whereas $\tau_{\infty}\simeq 280$min for Men’s. For Women’s soccer matches one that therefore does not expect a comeback of the weaker after about 60 minutes, whereas for Men’s soccer such a comeback is possible until the end of the game. Indeed, only 50% of the relaxation timescales in Men’s soccer can take place in a match of 90min (Fig. 8 inset). This suggests that shorter matches of less than 70min for Women’s soccer would remain exciting until the final whistle blows, because the weaker team would yet have chances to win the match.

It is interesting to compare soccer to other sports such as handball and basketball, see Fig. 8. We estimate the parameters $k_{+}$ and $k_{-}$ in a similar manner. For both handball and basketball, men and women, the scoring affinities are close to $A\simeq 0.14$. The scoring rates $\nu$ are given in the caption of Fig. 8, in which we can appreciate the difference compared to soccer. For handball and basketball the larger tails towards shorter timescales $\tau_{0}$ in the Marčenko-Pastur distributions indicate that the first action or comebacks can happen much more rapidly than for soccer. Handball, basketball and mens soccer share a game duration which falls within the relaxation time spectrum. The long timescale $\tau_{\infty}$ exceeds significantly the game duration, for all cases except women’s soccer.

VI Discussion

We have derived analytical expressions for the distribution and moments of the finite-time minimum and maximum values of continuous-time biased random walks. Such stochastic processes provide minimal models to describe the fluctuating motion of molecular motors and cyclic enzymatic reactions that take place in a thermal reservoir and under non-equilibrium conditions induced by e.g. external forces and/or chemical reactions.

Our key results are: (i) exact statistics of the extrema of the position and the entropy production of a biased random walk (ii) a novel connection between extreme-value statistics of biased random walks and the Marčenko-Pastur distribution of random matrix theory; (iii) symmetry relations between distributions of extrema of stochastic entropy production; (iv) estimates of extreme value statistics from spectral properties of random matrices.

For biased random walks, our results provide insights beyond the infimum law for nonequilibrium steady states, $\langle S_{\rm min}(t)\rangle\geq-1$, which states that the entropy production of a mesoscopic system plus its environment cannot be reduced on average by more than the Boltzmann constant. For continuous systems this bound is approached at large times, $\lim_{t\rightarrow\infty}\langle S_{\rm min}(t)\rangle=-1$. Here we have shown that the effects of discreteness are very important. At large times we find that a model dependent bound above $-1$ is reached, see Eq. (4) and Fig. 2. Equations (10) and (11) reveal that, for this class of models, the maximum difference between the entropy production and its initial value has the same statistics as the maximum difference between its maximum and its final value, for any given time interval $[0,t]$. Moreover, this result reveals that a ”supremum law” for entropy production bounds the average of the difference between entropy production’s maximum and its value at a fixed time $t\geq 0$, i.e.

[TABLE]

The inequality (21) follows applying the results of Sec. V in Ref. Neri et al. (2017) to the process $R(t)=\mathcal{P}(X_{[0,t]})/\mathcal{P}(\tilde{X}_{[0,t]})$, which is a martingale with respect to the time-reversed measure. We have shown that the inequalities for the average extrema of entropy production and displacement of a biased random walk saturate in the limit of small affinity $A\ll 1$. This limit corresponds to systems that exchange a small amount of heat with their environment — below the thermal energy $k_{\rm B}T$ — in each forward or backward step of the walker. For larger values of $A$, our analytical expressions reveal that details on the discreteness in the walker’s motion have strong influence in the extrema statistics. At large times, the time-asymmetric parameter $A$ fully determines the distribution of extrema, whereas extrema at finite times are determined by both $A$ and the time-symmetric rate constant $\nu$, see Eq. (2). Moreover, our results for 1D biased random walks also apply to homogeneous 2D biased random walks describing the motion of molecular motors, upon suitable definitions of effective affinity and rate parameters given by Eqs. (17) and (18).

We have also revealed a connection between the statistics of finite-time extrema of biased random walks with random-matrix theory, Eqs. (12) and (14). We show that the average of finite-time extrema (minimum and maximum) of a 1D biased random walk has a relaxation time spectrum given by Marčenko-Pastur distribution of eigenvalues of a random matrix ensemble, see Eq. (13). Therefore, efficient estimates for the extreme-value statistics of biased random walks can be obtained from the eigenvalue distributions of suitable Wishart and Laguerre random matrices Nica and Speicher (2006); Fridman et al. (2012); Livan et al. (2018). Our numerical simulations show that random-matrix estimates can overperform the accuracy and convergence of Monte Carlo simulations to determine extrema statistics. We remark that the connection between extrema and random-matrix theory developed here is of distinct nature to the relation between statistical properties of stochastic processes (e.g. interface growth) and extreme eigenvalues of random matrices described by a Tracy-Widom distribution Tracy and Widom (1994, 1996); Prähofer and Spohn (2000); Dhar et al. (2017). We expect that this surprising link of stochastic thermodynamics with random matrix theory may turn out to be more general, through the universal convergence of Wishart-Laguerre random matrices in the Marchenko-Pastur theorem, and could hold for larger classes of systems in a similar way.

Our work has important consequences for the theory of nonequilibrium fluctuations of active molecular processes and biomolecules. For example, the statistics of the maximum excursion of a motor against its net motion along a track provides insights on the physical limits of pernicious effects of fluctuations at finite times, which can be relevant in e.g. the finite-time efficiency of enzymatic reactions responsible of polymerization processes, muscle contraction by molecular motors, etc. We have shown that the displacement of motors with small cycle affinity exhibits large extreme values on average. However, the associated extreme entropy flows are on average always bounded in absolute value by the Boltzmann constant. Insights of our theory could be also discussed in the context of more complex biomolecular stochastic processes (e.g. microtubule growth Helenius et al. (2006); Redemann et al. (2017) and transport in actin networks Richard et al. (2019)). It will be interesting to extend our theory to Markovian and non-Markovian processes with time-dependent driving Harris and Schütz (2007); Strasberg et al. (2016); Harris and Shreshtha (2019); Proesmans et al. (2019); Loos and Klapp (2019); Oberreiter et al. (2019), stochastic processes with hidden degrees of freedom Polettini and Esposito (2017); Martínez et al. (2019), and also to explore whether extrema statistics from single-molecule experimental data reveal relaxation spectra described by the Marčenko-Pastur distribution.

We have illustrated our results in simple stochastic models of sports (soccer, handball and basketball) yielding the following conclusions. The Marčenko-Pastur distributions of the extreme-value relaxation times gives an interesting representation of the thrill of the game across timescales, which is both dependent on the rate $\nu$, quite characteristic of the sport, and the strength difference $A$ between the teams. The duration of the game acts as a cursor within this distribution that balances the predictability versus the randomness of the game, favoring the possibility for the weaker team to win from a great coup. Women’s soccer appears as an outlier in the comparison of several champions leagues for different sports and genders, that makes the outcome of the game more deterministic. In order to guarantee exciting games, women’s soccer matches duration may thus be optimized to shorter durations (60 to 75min, based on the women to men ratio of $\tau_{0}$ or $\tau_{1}$) without affecting the comeback possibilities.

We acknowledge fruitful discussions with Izaak Neri, Simone Pigolotti, Ken Sekimoto, Carlos Mejía-Monasterio, and enlightening discussions with Pierpaolo Vivo. AG acknowledges MPIPKS and ICTP for their hospitality, and Bertrand Fourcade for his guidance towards MPIPKS, and EUR Light S&T for funding.

Appendix

Appendix A First-passage-times, large deviations, and absorptions probabilities of biased random walks

In this section we review some knowledge of random walk theory (e.g. first-passage statistics Redner (2001)) to derive Eq. (6) in the Main Text i.e. the exact formula for the first-passage time distribution of a 1D continuous-time biased random walk. We also discuss large-deviation properties of this model, absorption probabilities using martingale theory, as well as parameters estimation from observed trajectories.

A.1 Model and solution of the Master equation

We consider a continuous-time biased random walk in a discrete one-dimensional lattice, where $X(t)=\{0,\pm 1,\pm 2,\dots\}$ denotes the position of the walker at time $t\geq 0$. We assume that $X(0)=0$ and that the walker can jump from state $X(t)=x$ to $x+1\,(x-1)$ at a rate $k_{+}\,(k_{-})$.

The waiting time at any site is exponentially distributed with the rate parameter $k_{+}+k_{-}$. The probability $P_{x}(t)=P\left(X(t)=x\right)$ to find the walker at the lattice site $x$ at time $t$ obeys the master equation

[TABLE]

with initial condition $P_{x}(0)=\delta_{x,0}$. From this evolution equation, a velocity $v=k_{+}-k_{-}$ and a diffusion coefficient $D=(k_{+}+k_{-})/2$ can be defined. Its solution is the Skellam distribution Skellam (1946) $P_{\rm Sk}(x;\mu_{1},\mu_{2})=(\mu_{1}/\mu_{2})^{x/2}I_{x}(2\sqrt{\mu_{1}\mu_{2}})e^{-(\mu_{1}+\mu_{2})}$ with parameters $\mu_{1}=k_{+}t$ and $\mu_{2}=k_{-}t$

[TABLE]

where $I_{x}(y)$ denotes the $x-$th order modiffied Bessel function of the first kind. Equation (23) follows from (22) and the exact expression for the generating function of the modified Bessel function of the first kind $\sum_{x=-\infty}^{\infty}z^{x}I_{x}(y)=e^{y(z+z^{-1})/2}$.

A.2 Large deviation and diffusion limit

We now discuss and review large deviation properties of the 1D biased random walk Touchette (2009); Speck et al. (2012) and relate them to the statistics of 1D drift diffusion process. For this purpose we consider the scaling limit of $P_{x}(t)$ given by Eq. (23) for large $x\sim vt$, with $v=k_{+}-k_{-}$ the net velocity of the walker. We assume a large deviation principle for $P_{x}(t)$ of the form

[TABLE]

In order to derive an analytical expression for the rate function $J$, we first approximate the modified Bessel function in Eq. (23) using a saddle point approximation

[TABLE]

where the saddle point is given by $i\sin\theta_{0}=z$ i.e. $i\theta_{0}=\sinh^{-1}z$ and we have used $\cos i\theta_{0}=\cosh\theta_{0}$ and $z^{-1}\cosh\sinh^{-1}z=\sqrt{1+z^{-2}}$. The rate function $J(u)$ with $u=x/vt$ can be evaluated from the leading term of Eq. (23) which is found using Eq. (26):

[TABLE]

where $z=u\sinh(A/2)=x/2\nu t$ and the change of variables $(k_{+},k_{-})\to(A,\nu)$ have been used for convenience.

In the vicinity of the minimum where $u\sim 1$ and thus $z\sim\sinh(A/2)$, the large deviation function behaves as

[TABLE]

Interestingly, the ratio between the second and the leading term of (28), given by $\tanh(A/2)^{2}(u-1)/3$, vanishes for small deviations $u\sim 1$ but also for large deviations in the limit of a small bias $A\ll 1$. The continuum limit of the biased random walk for $A$ small simplifies to the drifted Brownian motion $P_{x}(t)=e^{-(x-vt)^{2}/4Dt}/\sqrt{4\pi Dt}$, where the polynomial prefactor is recovered by normalization, and we have used the expressions for the velocity $v=2\nu\sinh(A/2)$ and the diffusion coefficient $D=\nu\cosh(A/2)$. In this regard, the 1D biased random walk can be seen as a generalization of the drifted Brownian motion for any bias. A finite bias modifies occurrences of extreme large deviations with respect to those occurring in the drift diffusion process. Consequently, the bias $A$ is expected to affect the extreme value statistics of the process, as shown below.

A.3 Martingales and absorption probabilities

In this subsection we employ martingale theory to derive an analytical expression of the absorption probability $P_{\rm abs}(-x)$ for a 1D biased random walk starting at $x=0$ to ever reach an absorbing boundary located at $-x<0$.

We first show explicitly that $e^{-S(t)}$ is a martingale process with respect to $X(t)$, i.e.

[TABLE]

for $t\geq t^{\prime}$. In words, the average of $e^{-S(t)}$ over all trajectories with common history $X_{[0,t^{\prime}]}$ up to time $t^{\prime}\leq t$ equals to its value at the last time of the conditioning $e^{-S(t^{\prime})}$. The proof is as follows:

[TABLE]

In the first and second lines we have used the additive property and the Markov property of entropy production, respectively. In the third line we have used the definitions $\Delta t\equiv t-t^{\prime}$ and $S(t)=X(t)\ln(k_{+}/k_{-})$. In the fourth line we have used the identity $\sum_{x=-\infty}^{\infty}z^{x}I_{x}(y)=e^{y(z+z^{-1})/2}$.

We remark that the proof sketched above can be simplified using the integral fluctuation relation $\langle e^{-S(t-t^{\prime})}\rangle=1$ in the second line, which holds for any $t\geq t^{\prime}$ Jarzynski (1997); Seifert (2005). It has been shown Neri et al. (2019, 2017) that the martingality of $e^{-S(t)}$ implies a set of integral fluctuation relations at stopping times

[TABLE]

where $\mathsf{T}$ is any bounded stopping time, i.e. a stochastic time at which the process $X(t)$ satisfies for the first time a certain criterion. In particular, Eq. (31) holds for the first-passage time $\mathsf{T}_{2}$ of $X(t)$ to reach any of two absorbing barriers located at $-x_{-}$ and $x_{+}$, with $x_{+}$ and $x_{-}$ two arbitrary positive integer numbers. When applying Eq. (31) to this particular stopping time, we can unfold the average in the left-hand side using the absorption probabilities

[TABLE]

where in the second line we have used the fact that $e^{-S(\mathsf{T}_{2})}=e^{-Ax_{+}}$ with probability $P_{\rm abs}(x_{+})$ and $e^{-S(\mathsf{T}_{2})}=e^{Ax_{-}}$ with probability $P_{\rm abs}(x_{-})$. In the third line we have used $P_{\rm abs}(x_{+})+P_{\rm abs}(x_{-})=1$, and in the fourth line Eq. (31). Solving the third line Eq. (32) for the absorption probability we obtain

[TABLE]

Taking the limit $x_{+}\to\infty$ in Eq. (33) we obtain the well-known analytical expression for the absorption probability

[TABLE]

which we used to derive the analytical expressions $-$ Eqs. (3) and (4) $-$ of the distribution and mean of the global minimum of entropy production in the biased random walk.

A.4 First-passage-time distribution

The first-passage-time density $P_{\textrm{fpt}}(t^{\prime};x)$ can be derived from the solution of the Master equation (22) with an absorbing boundary at site $x\neq 0$, with $x$ an integer number Redner (2001). It can also be derived from the distribution of the walker using Laplace transforms through the renewal equation:

[TABLE]

where $P_{0}(t)$ is the probability to be at a state at time $t$ when the system was at the same state at $t=0$. This convolution integral becomes a product in the Laplace domain, for any $x\neq 0$:

[TABLE]

where

[TABLE]

and

[TABLE]

We thus obtain, using Eqs. (37) and (38) in (36)

[TABLE]

In the above equations and in the following we will use the variables $\nu=(k_{+}k_{-})^{1/2}$ and $A=\ln(k_{+}/k_{-})$, see Eq. (2) in the Main Text. The inverse Laplace transform of Eq. (39) implies $P_{\textrm{fpt}}(t;x)=(|x|/t)P_{x}(t)$:

[TABLE]

which is Eq. (6) in the Main Text.

Appendix B Exact extrema statistics

for the 1D biased random walk

In this section, we use generating functions to derive the statistics of the finite-time extrema of entropy production. In particular, we focus on the generating function for the probability $G_{\min}$ of the minimum and the generating function for the absorption probability $G_{\rm abs}$, which are defined respectively as

[TABLE]

These two generating functions are related by

[TABLE]

Moments and probabilities follow taking derivatives of the generating functions (41) with respect to $z$

[TABLE]

For instance, inserting (43) into (44) to compute the first moment reduces to the simple expression

[TABLE]

B.1 Finite time statistics in the Laplace domain

We rewrite the probability-generating function $G_{\min}$ using the first-passage-time density derived above (39). First we use the fact the the first-passage-time density is the derivative of the absorption probability, $\hat{P}_{\textrm{fpt}}(s;x)=s\hat{P}_{\textrm{abs}}(x;s)-\delta_{x,0}$, which holds for any $x$. Equation (40) implies that the Laplace transform $\hat{P}_{\textrm{fpt}}(s;-x)$ of the first-passage-time probability to reach an absorbing site located at $-x$, with $x\geq 1$, can be expressed as the $x-$th power of the Laplace transform $\hat{P}_{\textrm{fpt}}(s;-1)$ of the first-passage time probability to reach $x=-1$:

[TABLE]

Consequently, the statistics of the minimum (maximum) can be expressed in terms of just $\hat{P}_{\textrm{fpt}}(s;-1)$ ($\hat{P}_{\textrm{fpt}}(s;1)$). Using Eqs. (41-43) and (47), we derive the generating functions in terms of $\hat{P}_{\textrm{fpt}}(s;-1)$

[TABLE]

Because $\hat{P}_{\textrm{fpt}}(s;-1)=e^{-A/2}(h(s)+\sqrt{h(s)^{2}-1})$ is algebraic, all the moments and probabilities, obtained from Eqs. (44-45), are algebraic expressions. For instance, the Laplace transform of the mean minimum is obtained directly from (46) and (48):

[TABLE]

B.2 Integral representations of extreme value statistics

We start from the first-passage-time density formula (40) and we exploit two properties of the modified Bessel function of the first kind. This allows us to rewrite the absorption probability as a definite integral of trigonometric and hyperbolic functions and the parameters $A$ and $\nu$:

[TABLE]

In (B.2) we have used Eq. (40). In (B.2) we have used the definition $I_{x}(y)=(1/\pi)\int_{0}^{\pi}e^{y\cos\theta}\cos(x\theta){\rm d}\theta$ and the property $I_{x}(y)=(y/2x)[I_{x-1}(y)-I_{x+1}(y)]$. Finally in (53) we have performed the integration over $t$. Using Eq. (53), we express the generating function of the absorption probability as an integral:

[TABLE]

where we have used in the above equation the identity

[TABLE]

which follows from the generating function of Chebyshev polynomials of the first kind $T_{x}(\cos\theta)\equiv\cos(x\theta)$. Performing the change of variable $y=\cos\theta$ in Eq. (54), and setting $z=1$ (we recall $\langle S_{\min}(t)\rangle=-AG_{\rm abs}(1,t)$, see Eq. (46)) we obtain Eq. (9) in the Main Text:

[TABLE]

Appendix C Extrema statistics and Marčenko-Pastur distribution

In this Section, we derive analytical expressions for the extrema statistics of 1D biased random walks in terms of the Marčenko-Pastur distribution (13) of random matrix theory, copied here for convenience:

[TABLE]

where $\lambda$ is a positive random variable, $\delta\leq 1$ a parameter, and $\lambda_{\pm}=\left(1\pm\sqrt{\delta}\right)^{2}$. Note that this distribution is normalized with mean $\int_{\lambda_{-}}^{\lambda_{+}}\lambda\rho(\lambda)\text{d}\lambda=1$.

C.1 Average value of the finite-time minimum of entropy production

Performing the changes of variable $k=2\nu(\cosh(A/2)-y)$, as well as $\tau=1/k$, in Eq. (LABEL:intrpz) we obtain

[TABLE]

where we have introduced the variables

[TABLE]

Thus, the average entropy production minimum can be expressed as exponential relaxation process with a spectrum of relaxation times distributed according to Marčenko-Pastur distributions (57)

[TABLE]

where

[TABLE]

and

[TABLE]

Note that Eq. (63) provides Eq. (12) of the Main Text.

C.2 Generating functions of the absorption probability and the minimum

The generating function for the absorption probability can also be expressed using Marčenko-Pastur distributions. Using the same method as described above for Eq. (54), we find

[TABLE]

where

[TABLE]

Using Eqs. (68) and (43) we obtain the generating function of the distribution of minima given by Eq. (14) in the Main Text

[TABLE]

C.3 Laplace transforms

Taking the Laplace transform of Eq. (67), we obtain

[TABLE]

Equations (71) and (46) imply that the Laplace transform of the average minimum can be written as a Stieltjes-like transform of the Marčenko-Pastur distribution for the variable $k$

[TABLE]

and similar relations hold for moments of any order. Notably, Eqs. (71) and (72) have a similar mathematical structure as the Laplace transform of the first-passage-time density of Markovian stochastic processes found in Hartich and Godec (2019), where instead $\hat{P}_{\rm fpt}(s;x)$ is expressed as a weighted discrete sum of relaxation modes.

Appendix D Random-matrix estimates

of extreme-value statistics

In this section we discuss the connection between the relaxation spectrum of first-passage and extrema statistics in the 1D biased random walk with random-matrix theory. We now describe how one can estimate finite-time statistics of the minimum entropy production from the spectrum of suitable random matrices. For this purpose, we use a celebrated result by Marčenko and Pastur Marčenko and Pastur (1967). Consider a real $m\times m$ Wishart matrix defined as

[TABLE]

where $\mathbf{R}$ (its transpose $\mathbf{R}^{T}$) is a $m\times n$ random matrix, with $n\geq m$. The random matrix $\mathbf{R}$ is filled with independent identically distributed (i.i.d.) random variables drawn from a normal distribution of zero mean and unit variance, i.e. $\mathbf{R}_{ij}\sim\mathcal{N}(0,1)$, for all $i,j\leq m$. The resulting positive definite random matrix follows the Wishart distribution of degree of freedom $n$ and density $c_{n,m}\text{det}\mathbf{w}^{(n-m-1)/2}e^{\operatorname{Tr}\mathbf{w}\,n/2}$ (where $c_{n,m}$ is a normalization factor). Following Marčenko and Pastur, the eigenvalues $\lambda$ of the Wishart random matrix $\mathbf{W}$ are asymptotically distributed according to the distribution (57) in the limit $n,m\to\infty$ with finite rectangularity $m/n=\delta<1$. It has been shown that this asymptotic result also holds when all $\mathbf{R}_{ij}$ are i.i.d. random variables drawn from any distribution of zero mean and unit variance Livan et al. (2018).

We now put in practice Marčenko and Pastur’s result, namely we find random matrices whose spectral density matches with that of the relaxation spectrum of the average minimum of entropy production. This can be achieved e.g. by using a Wishart random matrix of rectangularity $\delta=m/n$ identified as $k_{-}/k_{+}=e^{-A}$ in terms of the bias $A$ of the walker, i.e. we draw a real $m\times n$ Wishart random matrix $\mathbf{W}$ with $m,n\gg 1$ and $m/n\simeq e^{-A}$ (for instance $n=\lceil e^{A}m\rceil$). Then we evaluate the $m$ eigenvalues $\lambda_{i}$ of the matrix $\mathbf{W}$ and we give them a dimension using Eqs. (64-65) and performing the changes of variables $k=\lambda\bar{k}$ in Eq. (62) and $\tau=\lambda\bar{\tau}$ in Eq. (63) respectively. We thus obtain the following two estimators, $\langle\tilde{S}_{\min}(t)\rangle_{k}$ and $\langle\tilde{S}_{\min}(t)\rangle_{\tau}$, for the average entropy production minimum:

[TABLE]

where $\bar{k}=k_{+}$ and $\bar{\tau}=k_{+}/(k_{+}-k_{-})^{2}$ as identified previously. These estimators converge respectively to the exact result in the limit of a large matrix size. Using Eq. (67-68), the same procedure can be applied to estimate the generating function and any order moment of the distribution of entropy production extrema.

To test the convergence of these estimators, we define their relative error $\epsilon_{k}(t)$ and $\epsilon_{\tau}(t)$ as the relative difference

[TABLE]

which is a random real quantity for both $k$ and $\tau$ estimates. Their limiting values are related and can be calculated analytically:

[TABLE]

which vanishes in the limit of a large random matrix because $\langle 1/\lambda\rangle_{\rho}=1/(1-e^{-A})$, with $\langle...\rangle_{\rho}$ denoting an average over the Marčenko-Pastur distribution (57). From the limits (77), we conclude that the estimator (74) is advantageous to study the short-time behavior whereas (75) is most suited for large-time asymptotics. Our numerical results show that $|\epsilon_{k}(t)|\leq|\epsilon_{\max}|$ and also $|\epsilon_{\tau}(t)|\leq|\epsilon_{\max}|$ for all tested parameter values and for all times $t$. Therefore we will use $\epsilon_{\max}$ given by Eq. (78) as a conservative bound for the relative error of the random-matrix estimates at any time $t$.

The estimates introduced above rely on the fact that one can achieve a rectangularity $\delta=m/n\simeq e^{-A}$ with large enough random matrices. Because $e^{-A}$ is in general a real number, it is desirable to develop random-matrix estimators that achieve the Marčenko-Pastur distribution accurately for any value of $A$. Following Nica and Speicher (2006), the $\beta$-Laguerre matrices are an alternative ensemble whose spectral density tends asymptotically to the distribution $\text{MP}(e^{-A})$ in the large size limit. A $\beta$-Laguerre $m\times m$ random matrix $\mathbf{L}$ is defined as

[TABLE]

where $n=m\beta/\delta$. Here, $\mathbf{R}$ is a $m\times m$ random matrix with all entries equal to zero except the $m\times(m-1)$ diagonal and sub-diagonal elements. The non-zero entries $R_{ij}$ are drawn following $\chi(d_{ij})$ distributions of $d_{ij}$ degrees of freedom:

[TABLE]

The random variable $Z_{ij}=\sqrt{\sum_{k=1}^{d_{ij}}X_{k}^{2}}$, where the $X_{k}\sim\mathcal{N}(0,1)$ are independent and identically distributed (i.i.d.), follows the $\chi(d_{ij})$ distribution. Equivalently, one can obtain a random variable that follows the $\chi(d_{ij})$ distribution by taking the square root of a random variable drawn from a chi-square distribution $\chi^{2}(d_{ij})$. The degrees of freedom $d_{ij}$ of the $\chi$ distributions in the $\beta$-Laguerre random matrix are:

[TABLE]

We recall that in this case $n=m\beta/\delta$ is not the dimension of the matrix $\mathbf{R}$, which is a $m\times m$ square matrix, but a positive real number. Therefore, the rectangularity parameter does not need to be approximated in this method.

Figure 9 shows numerical results of the random-matrix estimators of the average entropy-production minimum for the 1D biased random walk with bias $A=1$. We draw $m\times m$ random matrices from the $\beta=1,2,1000$ Laguerre ensembles. Note that $\beta=1,2,4$ Laguerre ensembles are equivalent to Wishart random matrices with $\mathbf{R}_{ij}$ given respectively by real, complex and quaternionic normal random variables, and use the appropriate conjugate transpose of $\mathbf{R}$ Dumitriu and Edelman (2002). We plot the maximum relative error $\epsilon_{\rm max}$ (78) as a function of the size of the random matrix $m$ (Figure 9a). The Wishart ($1$-Laguerre) ensemble provides a biased overestimate of the real value with maximum relative error $2.3\%$ for small random matrices of sizes larger than $64\times 64$. We observe that the $2$-Laguerre ensemble provides an estimator that is practically unbiased, even using small matrices (except in the limit $A\ll 1$). Furthermore, $\beta$-Laguerre matrices with large values of $\beta$ (e.g. $\beta=1000$) yield small dispersion in the relative difference but a bias (underestimation) for small matrix sizes. The mean and the standard deviation of $\epsilon_{\rm max}$ obtained from a large population of computer-generated random matrices are observed to converge to zero with the matrix size $m$ as $\sim 1/m$. This fast convergence (compared to the usual $1/\sqrt{m}$) is a consequence of the correlation between the $m$ eigenvalues of the $\beta$-Laguerre random matrices. The convergence of the estimators is revealed in the difference between the spectral density of the random matrices and the Marčenko-Pastur distribution for $m=2^{6}$ (Fig. 9b) and $m=2^{12}$ (Fig. 9c). Remarkably, even though for $m=2^{6}$ the eigenvalue distribution is a rough approximation to the Marčenko-Pastur distribution, the relative error of the estimator is smaller than $\pm$2.3% for a single $1$-Laguerre and $\pm$1.5% for a single $2$-Laguerre random matrix.

Eventually, we notice that all the expressions that involve a sum over the eigenvalues can be recast into random matrix traces. For instance, the two minimum entropy estimators and their associated maximum relative error read:

[TABLE]

where $\mathbf{M}$ is either a Wishart random matrix or the Laguerre random matrix of size $m\times m$, scaled by its degree of freedom $n=e^{A}m$ (see Eqs. (73,79)).

Appendix E Explicit expression for the average entropy production minimum

We now employ exact expressions for the moments of the Marčenko-Pastur distribution $\rho$ to derive an analytical expression for the average minimum of the 1D biased random walk. The $n-$th moment ($n\geq 1$) of the Marčenko-Pastur distribution $\rho(\lambda)$, can be expressed as a sum involving binomial coefficients:

[TABLE]

where $\langle...\rangle_{\rho}$ denotes an average over the Marčenko-Pastur distribution (57). Then, it is convenient to expand the exponential in Eq. (62), and express the integrals in terms of the moments (85). By manipulating the indices, we obtain two infinite series:

[TABLE]

This form can be identified as the Kampé de Fériet function $F$ Exton (1978), a two-variable generalization of hypergeometric functions. It is defined for integer vectors $\bf a,b,b^{\prime},c,d,d^{\prime}$, of lengths $p,q,q,r,s$ and $s$ respectively, as follows:

[TABLE]

where $(m)_{n}=\prod_{k=0}^{n-1}(m-k)=\frac{(m+n-1)!}{(m-1)!}$ denotes the rising factorial. The translations of the binomial coefficient into rising factorials yields the following expression of the mean minimum entropy:

[TABLE]

Appendix F Sports data

We provide additional details of the data used in Sec. VI about sports, selected from European competitions during the season 2018-2019. For soccer, we use data from the UEFA Champions League in the group stage. We estimate the goal marking rate $k_{+}$ (goal conceding rate $k_{-}$) of the group leaders by summing the total number of goals scored (conceded) by all the group stage leaders divided by the cumulated duration of all the matches played in the group stage. There are 8 groups in Men’s and 10 groups in Women’s Champions League with 4 teams each. Each Men’s team plays 6 matches, each Women’s team plays 3 matches, each match lasts 90 minutes. Therefore, we select 48 matches of 8 teams out of 32 for Men’s soccer and 30 matches of 10 teams out of 40 for Women’s soccer.

For handball, we select the 2 best teams of each group (40 matches of 4 teams out of 12) in the main round of the Women’s EHF Champions League, 2 best teams of each group (56 matches of 8 teams out of 28) in the group stage of the Men’s EHF Champions League, each match lasts 60 minutes. For basketball, we select the 2 best teams of each group (56 matches of 4 teams out of 16) of the regular season of the Women’s FIBA EuroLeague (the equivalent of the Champions League which does not exists for Women’s basketball), and the best team or equally-ranked of each group (70 matches of 5 teams out of 32) of the regular season of the Men’s FIBA Champions League, each match lasts 40 minutes.