Dynamic disorder in simple enzymatic reactions induces stochastic amplification of substrate

TL;DR

This paper investigates how fluctuations in enzyme activity, known as dynamic disorder, can cause stochastic amplification of substrate levels in simple enzymatic reactions, with implications for understanding cellular processes.

Contribution

The study derives an explicit formula linking enzymatic fluctuation speed to steady-state substrate levels, revealing how dynamic disorder amplifies substrate concentration.

Findings

Fluctuation speed significantly affects mean substrate levels.

Large deviations from deterministic predictions occur due to enzyme activity fluctuations.

The connection between fluctuation speed and Markov process mixing properties is established.

Abstract

A growing amount of evidence points to the fact that many enzymes exhibit fluctuations in their catalytic activity, which are associated with conformational changes on a broad range of timescales. The experimental study of this phenomenon, termed dynamic disorder, has become possible due to advances in single-molecule enzymology measurement techniques, through which the catalytic activity of individual enzyme molecules can be tracked in time. The biological role and importance of these fluctuations in a system with a small number of enzymes such as a living cell have only recently started being explored. In this work, we examine a simple stochastic reaction system consisting of an inflowing substrate and an enzyme with a randomly fluctuating catalytic reaction rate that converts the substrate into an outflowing product. To describe analytically the effect of rate fluctuations on the average substrate abundance at steady-state, we derive an explicit formula that connects the relative speed of enzymatic fluctuations with the mean substrate level. We demonstrate that the relative speed of rate fluctuations can have a dramatic effect on the mean substrate, and lead to large positive deviations from predictions based on the assumption of deterministic enzyme activity. Our results also establish an interesting connection between the amplification effect and the mixing properties of the Markov process describing the enzymatic activity fluctuations, which can be used to easily predict the fluctuation speed above which such deviations become negligible. As the techniques of single-molecule enzymology continuously evolve, it may soon be possible to study the stochastic phenomena due to enzymatic activity fluctuations within living cells. Our work can be used to formulate experimentally testable hypotheses regarding the magnitude of these fluctuations, as well as their phenotypic consequences.

Tables

Table 1. Table 1: Estimates for ρ max subscript 𝜌 max \rho_{\textnormal{max}} , θ 𝜃 \theta and c ϵ subscript 𝑐 italic-ϵ c_{\epsilon} for various values of k s subscript 𝑘 𝑠 k_{s}

$k_s$	${\rho }_{\text{max}}$	$\theta$	$c_{ϵ}$
1	0.1703	2.1988	7.2903
5	0.1543	1.6907	8.5349
10	0.0930	0.6782	12.2383
20	0.0494	0.2550	15.4510

Equations

$$\mathbf{S}+\mathbf{E}\text{vbox}\text{ooalign}\text{raise}1.0pt\text{relbar}\text{joinrel}\rightharpoonup \text{joinrel}\text{crcr}\text{lower}1.0pt\leftharpoondown \text{joinrel}\text{relbar}\text{joinrel}S.E⟶\mathbf{P}+\mathbf{E}.$$

$$Q1=0,{\pi }^{T}Q=0^T\text{and}{\pi }^{T}1=1,$$

$$\mathbb{E}(\gamma (t))={\mathbb{E}}_{\pi }(\gamma )=\sum _{i=1}^n{\gamma }_i{\pi }_i\text{for all}t\ge 0,$$

$${\gamma }_c(t)=\gamma (ct)\text{for all}t\ge 0.$$

$$S_c(t)=S_c(0)+Y_1(k_{\text{in}}t)-Y_2\left({\int }_0^t{\gamma }_c(u)S_c(u)du\right),$$

$$\frac{d{S}_c(t)}{dt}=k_{\text{in}}-{\gamma }_c(t)S_c(t).$$

$$m_{\text{eq}}(c)=\underset{t\to \infty }{\lim }{m}_c(t).$$

$$m_c(t)=\mathbb{E}\left(S_c(0)e^{-{\int }_0^t{\gamma }_c(s)ds}\right)+k_{\text{in}}{\int }_0^t\mathbb{E}\left(e^{-{\int }_s^t{\gamma }_c(u)du}\right)ds.$$

$${\int }_0^t\mathbb{E}\left(e^{-{\int }_s^t{\gamma }_c(u)du}\right)ds={\int }_0^t\mathbb{E}\left(e^{-{\int }_0^{t-s}{\gamma }_c(u)du}\right)ds={\int }_0^t\mathbb{E}\left(e^{-{\int }_0^s{\gamma }_c(u)du}\right)ds.$$

$$m_{\text{eq}}(c)=\underset{t\to \infty }{\lim }{m}_c(t)=k_{\text{in}}{\int }_0^{\infty }\mathbb{E}\left(e^{-{\int }_0^s{\gamma }_c(u)du}\right)ds.$$

$${\int }_0^s{\gamma }_c(u)du=s\left(\frac{1}{cs}{\int }_0^{cs}\gamma (u)du\right).$$

$$\underset{c\to \infty }{\lim }{m}_{\text{eq}}(c)=k_{\text{in}}{\int }_0^{\infty }{e}^{-s{\mathbb{E}}_{\pi }(\gamma )}ds=\frac{k_{\text{in}}}{{\mathbb{E}}_{\pi }(\gamma )}:=m_{\text{eq}}^{\text{(det)}}.$$

$$\mathbb{E}\left(e^{-{\int }_0^s{\gamma }_c(u)du}\right)\ge e^{-{\int }_0^s\mathbb{E}({\gamma }_c(u))du}=e^{-s{\mathbb{E}}_{\pi }(\gamma )},$$

$$m_{\text{eq}}(c)\ge k_{\text{in}}{\int }_0^{\infty }{e}^{-s{\mathbb{E}}_{\pi }(\gamma )}ds=\frac{k_{\text{in}}}{{\mathbb{E}}_{\pi }(\gamma )}=m_{\text{eq}}^{\text{(det)}}.$$

$$\underset{c\to 0}{\lim }{m}_{\text{eq}}(c)=k_{\text{in}}\mathbb{E}\left({\int }_0^{\infty }{e}^{-s\gamma (0)}ds\right)=k_{\text{in}}\mathbb{E}\left(\frac{1}{\gamma (0)}\right)=k_{\text{in}}{\mathbb{E}}_{\pi }\left(\frac{1}{\gamma }\right):=m_{\text{eq}}^{\text{(static)}},$$

$$\inf \{x:x\in \Gamma \}\ge ϵ$$

$${\tau }_c=\inf \left\{t\ge 0:{\int }_0^t{\gamma }_c(s)ds=-\ln u\right\},$$

$$m_{\text{eq}}(c)=k_{\text{in}}\mathbb{E}({\tau }_c).$$

$$D=\text{Diag}({\gamma }_1,\dots ,{\gamma }_n)$$

$$m_{\text{eq}}(c)=k_{\text{in}}\left[{\pi }^T(D-cQ{)}^{-1}1\right].$$

$$m_{\text{eq}}^{\text{(det)}}\le m_{\text{eq}}(c)\le m_{\text{eq}}^{\text{(static)}}.$$

$$R(z)=(zI-\tilde{Q}{)}^{-1},$$

$$\tilde{Q}=\sum _{i=1}^n{\lambda }_i{u}_i{w}_i^T$$

$$R(z)=\sum _{i=1}^n\left(\frac{1}{z-{\lambda }_i}\right)u_i{w}_i^T$$

$$m_{\text{eq}}(c)=k_{\text{in}}\left[{\pi }^T(I-c\tilde{Q}{)}^{-1}{D}^{-1}1\right]=c^{-1}{k}_{\text{in}}\left[{\pi }^{T}R(c^{-1})D^{-1}1\right]$$

$${\alpha }_i=⟨\pi ,u_i⟩⟨w_i,D^{-1}1⟩=({\pi }^T{u}_i)(w_i^T{D}^{-1}1)\text{ for each }i=1,\dots ,n$$

$$m_{\text{eq}}(c)=k_{\text{in}}\sum _{i=1}^n\left(\frac{{\alpha }_i}{1-c{\lambda }_i}\right).$$

$${ϵ}_{\text{max}}=-\max \{\text{Re}({\lambda }_i):i=2,\dots ,n\},$$

$$m_{\text{eq}}(c)=k_{\text{in}}\left[{\alpha }_1+\sum _{i=2}^n\left(\frac{{\alpha }_i}{1-c{\lambda }_i}\right)\right].$$

$${\alpha }_1=\frac{m_{\text{eq}}^{\text{(det)}}}{k_{\text{in}}}\text{and}\sum _{i=2}^n{\alpha }_i=\left(\frac{m_{\text{eq}}^{\text{(static)}}-m_{\text{eq}}^{\text{(det)}}}{k_{\text{in}}}\right).$$

Full text

In this work, we examine a simple stochastic reaction system consisting of an inflowing substrate and an enzyme with a randomly fluctuating catalytic reaction rate that converts the substrate into an outflowing product. To describe analytically the effect of rate fluctuations on the average substrate abundance at steady-state, we derive an explicit formula that connects the relative speed of enzymatic fluctuations with the mean substrate level. Under fairly general modeling assumptions, we demonstrate that the relative speed of rate fluctuations can have a dramatic effect on the mean substrate, and lead to large positive deviations from predictions based on the assumption of deterministic enzyme activity. Our results also establish an interesting connection between the amplification effect and the mixing properties of the Markov process describing the enzymatic activity fluctuations, which can be used to easily predict the fluctuation speed above which such deviations become negligible. As the techniques of single-molecule enzymology continuously evolve, it may soon be possible to study the stochastic phenomena due to enzymatic activity fluctuations within living cells. Our work can be used to formulate experimentally testable hypotheses regarding the nature and magnitude of these fluctuations, as well as their phenotypic consequences.

Keywords: stochastic amplification; enzymatic fluctuations; dynamic disorder; Markov models

Mathematical Subject Classification (2010): 92C42; 92C45; 60J22; 60J28; 65C40.

1 Introduction

First made almost two decades ago, observations of enzymatic turnovers for single enzyme molecules have allowed scientists to probe enzyme behavior beyond the regime of high-copy numbers and ensemble averages [27]. Thanks to advances brought by experimental techniques such as single-molecule fluorescence spectroscopy [36, 26], the field of single-molecule enzymology developed rapidly in the subsequent years. The key observation made possible by single-molecule assays is that the catalytic rates of single enzyme molecules often display very large dynamic fluctuations over timescales much longer than the typical reaction cycle times, most likely driven by slow (spontaneous or induced) transitions in conformation [27, 25, 7, 42].

Around the time when the first single enzyme molecules were observed in action, the mathematical theory of dynamic disorder was introduced by Zwanzig [43], motivated by several observations of different physico-chemical processes with a seemingly common underlying cause that boiled down to random fluctuations of key process properties. The phenomenon of dynamic disorder refers to fluctuations in enzymatic reaction rates that occur at a timescale that is either slower or comparable to the reaction timescale [4]. These fluctuations are often caused by slow transitions in the conformational state of enzymes. The most simple example of dynamic disorder first considered by Zwanzig [43] involved a so-called “rate process controlled by passage through a fluctuating bottleneck” [44]. In the language of chemical kinetics, it describes the removal of a substrate, $S$, from a system at rate $\gamma(t)S(t)$, where the time-varying rate $\gamma(t)$ is a (typically Markovian) stochastic process (see Figure 1). As the speed of $\gamma(t)$ fluctuations tends to zero, the reaction rate $\gamma(t)$ becomes a random variable which does not change with time, and we transition to the regime of static disorder [43]. On the other hand, as the speed of $\gamma(t)$ fluctuations tends to infinity, the dynamic disorder vanishes on the timescale of substrate kinetics, and we recover the classical case where the reaction rate $\gamma(t)$ becomes a deterministic constant. Our goal in this paper is to investigate the effects of $\gamma(t)$ fluctuations between these two extremities, where most realistic systems are likely to lie.

Thanks to the mathematical theory of dynamic disorder, stochastically fluctuating enzyme activities can be understood and studied within a consistent mathematical framework that can also generate testable experimental predictions [37, 20, 34, 4, 7, 27]. Already in [43] it was observed that dynamically disordered systems can give rise to macroscopic observations that differ from those expected in the absence of disorder. In subsequent years, a large body of theoretical and computational work has examined various alternative enzymatic reaction schemes, mostly focusing on the enzyme dynamics itself, e.g. on the autocorrelation of fluctuations and the distribution of waiting times between turnover events [29]. On the other hand, dynamic disorder has been observed for several biologically relevant enzymes [41], suggesting that it is ubiquitous in the cellular context. Early work [7] had already noted that enzymatic fluctuations could play an important biological role in a system containing only a small number of enzyme molecules, as often happens within a living cell, and the recent in vivo observation of fluctuating enzymatic activity confirms this claim [17].

Besides studying the intrinsic mathematical properties of dynamically disordered enzymes, it would be also highly instructive and relevant for biology to examine the consequences of dynamic disorder on substrate statistics (a first example of such a study is given in [17]). Experimental work in this area is still done in vitro using constant and large substrate pools. Here, on the contrary, we provide a mathematical treatment of how dynamic disorder alters the substrate mean abundance in the presence of substrate inflow, a condition closer to biological reality that has also been considered in [13, 40]. To this end, we analytically examine a highly simplified stochastic system with a randomly fluctuating catalytic reaction rate and describe the effect of rate fluctuations on the average substrate abundance. Under fairly general conditions, we demonstrate that the relative speed of rate fluctuations can have a dramatic effect on the mean substrate, and lead to large positive deviations from predictions based on the assumption of deterministic enzyme activity. Using a Markovian model for enzyme kinetics, we mathematically characterize this effect by deriving an explicit formula for the steady-state substrate-mean as a function of the relative speed of enzymatic fluctuations. From this formula we show that for any finite speed-value, the steady-state substrate-mean is sandwiched between the two values obtained in the static and the deterministic regimes. Furthermore, we demonstrate that the mapping between the relative speed of enzyme kinetics and the substrate-mean at steady-state can be well-approximated by a convex, monotonically decreasing function whose key shape parameter depends on the “mixing strength” of the Markov process describing enzyme kinetics. This mixing strength can be measured by computing an appropriate Dirichlet form [35] of the Markov process. Even though we consider a highly simplified situation, our analysis can serve as a guide in the case of more realistic, but analytically intractable enzymatic reaction schemes. Our results only depend on the enzymatic fluctuations, but they do not depend on the fluctuations caused by the low abundance of substate molecules (see [6]), although we account for these fluctuations by modelling the substrate kinetics as a jump Markov chain. Indeed the results we present remain unchanged even if we discard these fluctuations and describe the substrate kinetics as an ordinary differential equation (ODE) with a fluctuating rate constant $\gamma(t)$.

Enzymatic fluctuations can also arise from sources other than dynamic disorder. The abundance or the availability of enzyme molecules may also fluctuate due to gene expression noise [6, 31], and their chances of finding substrate molecules can be diffusion-limited [39]. In this work, we do not distinguish between the various sources of fluctuations and model the aggregate enzymatic activity by a Markovian stochastic process.

The biological significance of our findings is manifold since enzymatic interactions are ubiquitous is cell biology and the effects of enzymatic noise in metabolic networks have only recently started to be explored [31]. Using the relative speed of enzymatic fluctuations as parameter, our results provide a clear way to determine if the deterministic approximation is a faithful representation of reality. Our results can shed light on the timescale disparities that exist between enzyme and substrate kinetics. In particular, we see that enzyme kinetics needs to be “fast” in order to avoid any undesirable amplification of the mean substrate abundance due to inevitable variations in the enzymatic states. On the other hand, one can envisage situations where it would be beneficial for enzymes to be “slow” so that their fluctuations amplify a weak signal and enable its detection by the intracellular machinery (see Section 3.1). Such a signal-detection mechanism was the main motivation behind stochastic focusing, a sensitivity amplification phenomenon introduced in [32]. We illustrate our results on the reaction scheme of [32] in Section 3.2, where we characterize how the substrate-mean changes with the speed of the enzyme abundance dynamics. Note that in situations where enzyme kinetics is “slow”, undesirable amplification effects can be eliminated by feedback mechanisms [28]. In this context, our results can help in postulating the presence of feedback loops using experimental data. We discuss the biological importance of our results in greater detail in Section 4.

It is interesting to note that some of the expressions we derive are related to those obtained in the analysis of various physico-chemical quantum dynamical systems coupled to a randomly fluctuating environment. This theory dates back to the original work of Kubo and Anderson [1, 21] and, more recently, has been generalized to arbitrary quantum systems described by the Liouvile-von Neumann equation with Markovian and non-Markovian parametric noise [10, 12, 11]. The example system of Ref. [11] can be interpreted as a substrate decaying with a stochastically fluctuating rate. While similar to it, the system we consider here includes the inflow of substrate, which results in a non-zero steady-state and requires a different mathematical treatment.

2 Results

2.1 The model

We consider a system into which a substrate, ${\bf S}$, enters at a constant rate $k_{\textnormal{in}}$ and is degraded or (equivalently) converted into a product ${\bf P}$ that in turn leaves the system. The rate of substrate outflow depends on the activity state or abundance level of an enzyme ${\bf E}$. In turn, the catalytic activity of ${\bf E}$, denoted by $(\gamma(t))_{t\geq 0}$, is assumed to fluctuate in time $t$ according to a continuous-time Markov chain (CTMC) over a finite state-space $\Gamma=\{\gamma_{1},\dots,\gamma_{n}\}$. Here each $\gamma_{i}$ is a positive constant denoting the degradation rate constant at the $i$-th enzymatic state or abundance level. Due to fluctuations in the catalytic activity of ${\bf E}$, the degradation rate of substrate ${\bf S}$ will also fluctuate in time according to a stochastic process $(k_{d,S}(t))_{t\geq 0}$ whose value at time $t$ is given by $k_{d,S}(t)=\gamma(t)S(t)$, where $S(t)$ is the molecular count or concentration of the substrate. This model is summarised in Figure 1.

As mentioned before, the degradation reaction ${\bf S}\stackrel{{\scriptstyle\gamma(t)}}{{\longrightarrow}}{\bf\emptyset}$ can also be viewed as a conversion reaction ${\bf S}\stackrel{{\scriptstyle\gamma(t)}}{{\longrightarrow}}{\bf P}$ which is catalyzed by the enzyme. Generally this catalytic step proceeds through the reversible formation of an intermediate complex ${\bf S.E}$ which is formed when an enzyme molecule binds to a substrate molecule. In other words, the single reaction ${\bf S}\stackrel{{\scriptstyle\gamma(t)}}{{\longrightarrow}}{\bf P}$ is an abstraction for the following three reactions:

[TABLE]

If the binding/unbinding rates of ${\bf S}$ and ${\bf E}$ molecules is much higher than the rate of the conversion reaction, then we can apply the quasi-stationary assumption to conclude that the model in Figure 1 is a good approximation to the catalytic conversion dynamics (for more details see the Supplementary Material in [28]).

To describe the CTMC $(\gamma(t))_{t\geq 0}$, we need to specify its $n\times n$ transition rate matrix $Q=[q_{ij}]$ (see [30]). For any distinct $i,j\in\{1,2,\dots,n\}$, $q_{ij}\geq 0$ denotes the rate at which the process leaves state $\gamma_{i}$ and enters state $\gamma_{j}$. The diagonal entries of $Q$ are given by $q_{ii}=-\sum_{j\neq i}q_{ij}$. From now on we assume that the rate matrix $Q$ is irreducible111A matrix $Q$ is called irreducible if there does not exist a permutation matrix $P$ such that the matrix $PQP^{-1}$ is block upper-triangular. which implies that there exists a unique stationary distribution $\pi=(\pi_{1},\dots,\pi_{n})\in\mathbb{R}^{n}_{+}$ satisfying

[TABLE]

where ${\bf 0}$ and ${\bf 1}$ denote the $n\times 1$ vectors of all zeroes and ones respectively. Since the state-space is finite and the transition rate matrix $Q$ is irreducible, the CTMC $(\gamma(t))_{t\geq 0}$ is ergodic which means that the probability distribution of $\gamma(t)$ converges to the stationary distribution $\pi$ as $t\to\infty$. As we are interested in the steady-state limit, without loss of generality we can assume that the initial state $\gamma(0)$ is distributed according to $\pi$, i.e. $\mathbb{P}(\gamma(0)=\gamma_{i})$ for each $i=1,\dots,n$. This ensures that the process $(\gamma(t))_{t\geq 0}$ is a stationary stochastic process whose various statistical properties do not depend on time222For a more rigorous definition of a stationary stochastic process see the Supplementary Material.. In particular its mean $\mathbb{E}(\gamma(t))$ is equal to

[TABLE]

where $\gamma$ is a $\Gamma$-valued random variable with probability distribution $\pi$ and $\mathbb{E}_{\pi}(\cdot)$ denotes the expectation w.r.t. this distribution.

From now on, we regard $(\gamma(t))_{t\geq 0}$ as the baseline process which corresponds to enzymatic dynamics at the natural timescale. In order to study the substrate behavior, we need to view enzymatic dynamics at the timescale of substrate kinetics. For this we define a family of processes $(\gamma_{c}(t))_{t\geq 0}$ parameterised by the “relative speed” parameter $c$ as follows:

[TABLE]

Note that one time-unit of process $(\gamma_{c}(t))_{t\geq 0}$ corresponds to $c$ time-units of process $(\gamma(t))_{t\geq 0}$. In this sense, the parameter $c$ sets the speed of the fluctuation dynamics for the enzyme relative to the speed of the substrate kinetics. Like $(\gamma(t))_{t\geq 0}$, the process $(\gamma_{c}(t))_{t\geq 0}$ is also a CTMC over state-space $\Gamma=\{\gamma_{1},\dots,\gamma_{n}\}$ with transition rate matrix $Q_{c}=cQ$ and initial distribution $\pi$. Since $(\gamma(t))_{t\geq 0}$ is stationary, this process is also stationary with the same mean given by $\mathbb{E}_{\pi}(\gamma)=\mathbb{E}(\gamma_{c}(t))$ for all times $t\geq 0$. Replacing $(\gamma(t))_{t\geq 0}$ by $(\gamma_{c}(t))_{t\geq 0}$ in the model depicted in Figure 1, we will study how the steady-state mean of substrate abundance depends on the fluctuation speed $c$.

Given a sample path of the enzyme dynamics $(\gamma_{c}(t))_{t\geq 0}$ with relative speed $c$, we regard the dynamics of substrate molecular counts as a jump Markov chain $(S_{c}(t))_{t\geq 0}$ over the set of nonnegative integers $\mathbb{N}_{0}=\{0,1,2,\dots\}$. This Markov chain can be written in the random time change representation [8] as

[TABLE]

where $Y_{1}$ and $Y_{2}$ are independent, unit rate Poisson processes. From this representation it is immediate that the substrate-production rate is constant ($k_{\textnormal{in}}$) in time, but the substrate-degradation rate is time-varying and it is equal to $\gamma_{c}(t)S_{c}(t)$ at time $t$. Here the Poisson processes $Y_{1}$ and $Y_{2}$ capture the intermittency in the firing of production and degradation reactions. This intermittency becomes unimportant if the substrate is present in high copy-numbers [22] and in this case one can regard $(S_{c}(t))_{t\geq 0}$ as the dynamics of subtrate concentration333The concentration of any species is its copy-number divided by the system volume., specified by the following ODE

[TABLE]

Note that even if the intermittency in production/degradation reactions is ignored and $(S_{c}(t))_{t\geq 0}$ is described by the ODE (2.4), the process $(S_{c}(t))_{t\geq 0}$ is still stochastic because it is driven by the stochastic process $(\gamma_{c}(t))_{t\geq 0}$ that represents enzymatic fluctuations.

Let $m_{c}(t)=\mathbb{E}(S_{c}(t))$ for each $t\geq 0$. We shall soon see that $m_{c}(t)$ does not depend on whether we use representation (2.3) or (2.4) for the substrate dynamics $(S_{c}(t))_{t\geq 0}$. Our goal in this paper is to understand the role of fluctuations in the catalytic activity of enzyme ${\bf E}$ in determining the steady-state value of the mean

[TABLE]

In particular, we study how this steady-state mean $m_{\textnormal{eq}}(c)$ depends on the relative fluctuation speed $c$ and the variability in degradation rates $\gamma_{1},\dots,\gamma_{n}$ at various enzymatic activity levels.

2.2 Expressions for $m_{\textnormal{eq}}(c)$: The general case

We can approximately find $m_{\textnormal{eq}}(c)$ by estimating $m_{c}(t)$ for a very large $t$, using simulations of the whole system. However this naive approach is highly unsatisfactory because these simulations can be computationally expensive and the approximation error incurred by replacing the steady-state mean by a finite-time mean is generally difficult to quantify. Moreover this approach does not provide us with an explicit formula for $m_{\textnormal{eq}}(c)$ that can enable us to study its dependence on the relative speed parameter $c$. In light of these difficulties, we look for alternative ways to compute $m_{\textnormal{eq}}(c)$. In this section we assume that enzymatic kinetics is given by a general stationary stochastic process with an arbitrary state-space $\Gamma\subset(0,\infty)$, and so we do not rely on the CTMC structure mentioned in Section 2.1. We specialise the results of this section to the CTMC case in Section 2.3.

Using representation (2.3) or (2.4) we can show that $m_{c}(t)=\mathbb{E}(S_{c}(t))$ is given by the following formula

[TABLE]

From the stationarity of the process $(\gamma_{c}(t))_{t\geq 0}$ we can conclude that

[TABLE]

Substituting this in (2.6) and letting $t\to\infty$, we obtain our first formula for $m_{\textnormal{eq}}(c)$, which is,

[TABLE]

From (2.2) we obtain

[TABLE]

Since the process $(\gamma(t))_{t\geq 0}$ is stationary, from Theorem 10.6 in [18] we know that as $c\to\infty$, the quantity (2.8) converges a.s. to $s\mathbb{E}_{\pi}(\gamma)$ (recall (2.1)). As a consequence $\mathbb{E}\left(e^{-\int_{0}^{s}\gamma_{c}(u)du}\right)\to e^{-s\mathbb{E}_{\pi}(\gamma)}$ and hence we get

[TABLE]

This shows that as the relative speed $c$ of enzymatic fluctuations approaches $\infty$, these fluctuations become equilibrated at the timescale of substrate kinetics, and so they do not affect the mean substrate level. In other words, from the point of view of the substrate, the enzyme kinetics is so fast that it is as if the enzyme state is constant at the equilibrium level $\mathbb{E}_{\pi}(\gamma)$. This corresponds to the classical case where there is no dynamic disorder in the enzyme activity and so this activity is well-approximated by a deterministic rate constant for the substrate degradation reaction. As the mapping $x\mapsto e^{-x}$ is convex, Jensen’s inequality tells us that

[TABLE]

where the last relation follows from the fact that $(\gamma_{c}(t))_{t\geq 0}$ is a stationary process with mean $\mathbb{E}(\gamma_{c}(t))=\mathbb{E}_{\pi}(\gamma)$ for all times $t\geq 0$. Substituting this in (2.7) we see that for any $c\geq 0$

[TABLE]

Therefore for a finite relative speed $c$, enzymatic fluctuations always amplify the mean substrate abundance, in comparison to the classical deterministic case. The natural question that now arises is - how large should speed $c$ be in order for the deterministic approximation to be acceptable within a certain tolerance level $\epsilon$? We address this question in Section 2.4.

Let us now consider the situation where the relative speed parameter $c\to 0$ and so at the timescale of substrate kinetics, the enzyme dynamics $(\gamma_{c}(t))_{t\geq 0}$ approaches a static process, i.e. $\gamma_{c}(t)=\gamma(0)$ for all $t\geq 0$. This case corresponds to the situation where the enzyme kinetics is very slow in comparison to the substrate kinetics. Hence from the point of view of the substrate, the kinetics of the enzyme is almost fixed. In this regime, we can replace $\gamma_{c}(u)$ by $\gamma(0)$ in (2.7) to obtain

[TABLE]

where we have used the fact that $\gamma(0)$ has probability distribution $\pi$ to write $\mathbb{E}(1/\gamma(0))$ as $\mathbb{E}_{\pi}(1/\gamma)$. Observe that $m^{\textnormal{(static)}}_{\textnormal{eq}}\geq m^{\textnormal{(det)}}_{\textnormal{eq}}$, which can be readily seen by letting $c\to 0$ in (2.10) or by directly using Jensen’s inequality on the convex map $f(x)=1/x$ (see Figure 2). The two extremal cases $c\to 0$ and $c\to\infty$ serve as a guide to the behavior of realistic systems with an intermediate value of $c$. In particular we can expect that for such intermediate $c$-values, the steady-state substrate mean will lie somewhere between $m^{\textnormal{(det)}}_{\textnormal{eq}}$ and $m^{\textnormal{(static)}}_{\textnormal{eq}}$ . This is precisely what happens as we shall soon see. We will also discuss how the precise value of $m_{\textnormal{eq}}(c)$ can be computed or estimated from any Markovian model of enzymatic fluctuations.

Until now, the conclusions we have drawn regarding $m_{\textnormal{eq}}(c)$ rely on the formula (2.7) that holds for any real-valued stationary stochastic process $(\gamma_{c}(t))_{t\geq 0}$ as long as its states are positive and bounded away from [math], i.e. the state-space $\Gamma$ satisfies

[TABLE]

for some $\epsilon>0$. This makes this formula very general but it is difficult to work with, because it involves an indefinite integral which is generally analytically intractable as the mapping $s\mapsto\mathbb{E}\left(e^{-\int_{0}^{s}\gamma_{c}(u)du}\right)$ does not have an explicit form. We remedy this problem in the next section by specialising this formula to the case where $(\gamma_{c}(t))_{t\geq 0}$ is a finite state-space CTMC as described in Section 2.1. Before we come to that, we provide a numerical recipe for statistically estimating $m_{\textnormal{eq}}(c)$ without the need for evaluating the indefinite integral. This scheme is based on the assumption that we can efficiently generate sample-paths of the stationary process $(\gamma_{c}(t))_{t\geq 0}$ (see [9, 2]).

Define a random variable $\tau_{c}$ by

[TABLE]

where $u$ is an independent random variable with the uniform distribution on $[0,1]$. To sample $\tau_{c}$ we can adopt the following strategy. We first sample $u$ from the uniform distribution on $[0,1]$, draw an initial condition $\gamma_{c}(0)$ from $\pi$, and then simulate the sample path $(\gamma_{c}(t))_{t\geq 0}$, keeping track of the integral $\int_{0}^{t}\gamma_{c}(s)ds$. We take $\tau_{c}$ to be first time $t$ when this integral hits the value $(-\ln u)$. From the samples of the random variable $\tau_{c}$, we can estimate its expectation $\mathbb{E}(\tau_{c})$ which gives us an estimate for $m_{\textnormal{eq}}(c)$ because it can be shown that

[TABLE]

Note that the estimator for $m_{\textnormal{eq}}(c)$ based on formula (2.14) will be unbiased but it will suffer from statistical error due to a finite sample size. However this error can be estimated and managed far more easily than the error one would incur by approximating the steady-state mean $m_{\textnormal{eq}}(c)$ by a finite-time mean $m_{c}(t)$ (recall (2.5)). This makes this formula (2.14) useful in practice (see Example 3.2).

The results from this section are collected in our next proposition which is proved in the Supplementary Material.

Proposition 2.1

Suppose $(\gamma(t))_{t\geq 0}$ is a real-valued stationary stochastic process with stationary distribution $\pi$ and state-space $\Gamma$ satisfying (2.12). Let $(\gamma_{c}(t))_{t\geq 0}$ be the speed $c$ version of this process given by (2.2) and define the substrate dynamics $(S_{c}(t))_{t\geq 0}$ either by (2.3) or by (2.4). Let $m_{c}(t)=\mathbb{E}(S_{c}(t))$ and let the steady-state limit $m_{\textnormal{eq}}(c)$ be given by (2.5). Then we have the following:

• (A)

The value $m_{\textnormal{eq}}(c)$ is well-defined (i.e. the limit in (2.5) exists) and it is given by (2.7).

• (B)

If $\tau_{c}$ is the random variable defined by (2.13) then (2.14) holds.

• (C)

The limits (2.9) and (2.11) are satisfied as $c\to\infty$ and $c\to 0$ respectively.

2.3 Expressions for $m_{\textnormal{eq}}(c)$: The finite CTMC case

In this section we specialize expression (2.7) to the case where $(\gamma(t))_{t\geq 0}$ is a stationary CTMC with a finite state-space $\Gamma=\{\gamma_{1},\dots,\gamma_{n}\}$ as described in Section 2.1. Define the CTMC $(\gamma_{c}(t))_{t\geq 0}$ by (2.2) and recall that its $n\times n$ transition-rate matrix is given by $Q_{c}=cQ$. Let $D$ be the $n\times n$ diagonal matrix

[TABLE]

whose entries are the degradation rates at different enzymatic states or abundance levels. One of the main results in our paper is to show that $m_{\textnormal{eq}}(c)$ can be expressed as

[TABLE]

Two alternative proofs of this result are given in the Supplementary Material. The first proof exploits some ideas from the theory of occupation measures for Markov chains [14] while the second proof is based on the Methods of Conditional Moments (MCM) approach recently developed by Hasenauer et al. [15]. Note that this formula assumes matrix $(D-cQ)$ is invertible for any $c\geq 0$ but this can be easily verified from the properties of matrix $Q$. Using formula (2.16) we can prove that for any $c\geq 0$

[TABLE]

Therefore for any relative speed $c$ of enzymatic fluctuations, the steady-state mean of the substrate is always sandwiched between the values obtained for the deterministic and the static cases. Moreover since $m_{\textnormal{eq}}(c)$ depends continuously on $c$, the limits (2.11) and (2.9) imply that for any value $m^{*}$ in the open interval $(m^{\textnormal{(det)}}_{\textnormal{eq}},m^{\textnormal{(static)}}_{\textnormal{eq}})$ there exists a relative speed value $c^{*}>0$ such that $m_{\textnormal{eq}}(c^{*})=m^{*}$. Hence the positive deviations caused by enzymatic fluctuations (in the mean substrate abundance) range from [math] to exactly $(m^{\textnormal{(static)}}_{\textnormal{eq}}-m^{\textnormal{(det)}}_{\textnormal{eq}})$.

To detemine the map $c\mapsto m_{\textnormal{eq}}(c)$ we need to evaluate $m_{\textnormal{eq}}(c)$ at several values of $c$. This can be difficult with formula (2.16) because each evaluation requires inversion of a potentially large matrix. Fortunately we can resolve this issue using simple ideas from the theory of resolvents for linear operators [19], as we now describe. Let $\widetilde{Q}=D^{-1}Q$ be the transition rate matrix of another CTMC over state-space $\Gamma=\{\gamma_{1},\dots,\gamma_{n}\}$. The difference between this new CTMC and the original CTMC is that the rates of outflow from state $\gamma_{i}$ to each state $\gamma_{j}$ (for $j\neq i$) are divided by the state value $\gamma_{i}$. Let $\mathbb{C}$ denote the field of complex numbers. The resolvent for the Markov semigroup corresponding to this CTMC is the matrix-valued function over $\mathbb{C}$ defined by

[TABLE]

where $I$ is the $n\times n$ Identity matrix. This function is well-defined for any $z$ which is not an eigenvalue of matrix $\widetilde{Q}$. Let $\lambda_{1},\dots,\lambda_{n}$ be the $n$ eigenvalues of matrix $\widetilde{Q}$, repeated according to their algebraic multiplicity. Since $\widetilde{Q}$ is the transition rate matrix of a CTMC, it has a simple444An eigenvalue is said to be simple if its algebraic multiplicity is $1$. eigenvalue (say $\lambda_{1}$) equal to [math], while its other eigenvalues have negative real parts. This implies that the resolvent function $R$ is well-defined on the positive real line $(0,\infty)$.

From now on we assume that matrix $\widetilde{Q}$ is diagonalizable555A square matrix $M$ is diagonalizable if it can be written as $M=P\Lambda P^{-1}$ for some diagonal matrix $\Lambda$ and some invertible matrix $P$. The diagonal entries of $D$ are the eigenvalues of matrix $M$. over the field $\mathbb{C}$ of complex numbers. This assumption is not very restrictive because almost every matrix is diagonalisable (see [16]) and so if $\widetilde{Q}$ is not diagonalisable, we can perturb matrix $Q$ slightly to make $\widetilde{Q}$ diagonalisable and not affect the enzyme dynamics significantly. The diagonalizability of $\widetilde{Q}$ allows us to write matrix $\widetilde{Q}$ as $\widetilde{Q}=U\Lambda U^{-1}$, where $\Lambda=\textnormal{Diag}(\lambda_{1},\dots,\lambda_{n})$ and $U$ is an invertible matrix whose columns contain the right eigenvectors for matrix $\widetilde{Q}$ corresponding to the eigenvalues $\lambda_{1},\dots,\lambda_{n}$. Similarly the rows of $U^{-1}$ contain the left eigenvectors for matrix $\widetilde{Q}$ corresponding to the eigenvalues $\lambda_{1},\dots,\lambda_{n}$. Let $u_{i}$ and $w_{i}$ be $n\times 1$ vectors denoting the $i$-th column and $i$-th row of matrices $U$ and $U^{-1}$ respectively. Therefore

[TABLE]

and we can express the resolvent function $R$ (see Chapter 5 in [19]) as

[TABLE]

Let $\langle\cdot,\cdot\rangle$ denote the standard inner product on $\mathbb{R}^{n}$. Note that formula (2.16) can be expressed as

[TABLE]

for any $c>0$. Plugging $R(c^{-1})$ from (2.20) and defining

[TABLE]

we obtain the following formula for $m_{\textnormal{eq}}(c)$:

[TABLE]

Observe that since $\alpha_{i}$-s and $\lambda_{i}$-s are independent of $c$, they only need to be computed once to construct this expression and then we can easily compute $m_{\textnormal{eq}}(c)$ for several values of $c$ without the need of evaluating the matrix inverses in (2.16). Moreover if $n$ is large, then using the values of $\alpha_{i}$ and $\lambda_{i}$ as a guide, one can derive suitable approximations of the formula (2.22) for $m_{\textnormal{eq}}(c)$. We derive one such approximation in the next section and use it as a tool to further understand the phenomenon of stochastic amplification induced by dynamic disorder in enzymatic activity.

The results from this section are collected in our next theorem which is proved in the Supplementary Material.

Theorem 2.2

Suppose $(\gamma(t))_{t\geq 0}$ is a stationary CTMC with transition rate matrix $Q$, stationary distribution $\pi$ and state-space $\Gamma=\{\gamma_{1},\dots,\gamma_{n}\}$ (see Section 2.1). Let $(\gamma_{c}(t))_{t\geq 0}$ be the speed $c$ version of this process given by (2.2) and define the substrate dynamics $(S_{c}(t))_{t\geq 0}$ either by (2.3) or by (2.4). Let the steady-state substrate mean $m_{\textnormal{eq}}(c)$ be given by (2.5) and the diagonal matrix $D$ be defined by (2.15). Then we have the following:

• (A)

The matrix $(D-cQ)$ is invertible and $m_{\textnormal{eq}}(c)$ can be expressed as (2.16).

• (B)

Suppose the matrix $\widetilde{Q}=D^{-1}Q$ is diagonalizable and let $\lambda_{1},\dots,\lambda_{n}$ be its eigenvalues. For each $i=1,\dots,n$ define $\alpha_{i}$ by (2.21). Then $m_{\textnormal{eq}}(c)$ can be expressed as (2.22).

• (C)

The relation (2.17) is satisfied for any $c\geq 0$.

2.4 Approximate formula for $m_{\textnormal{eq}}(c)$

The goal of this section is to derive an approximate formula for $m_{\textnormal{eq}}(c)$ using (2.22) and then use it to obtain some interesting insights. Recall from the previous section that $\lambda_{1},\dots,\lambda_{n}$ are the eigenvalues of matrix $\widetilde{Q}$. Among these $\lambda_{1}=0$ while the eigenvalues $\lambda_{2},\dots,\lambda_{n}$ have negative real parts. Define a positive constant $\epsilon_{\textnormal{max}}$ by

[TABLE]

where $\textnormal{Re}(z)$ denotes the real part of a complex number $z$. Setting $\lambda_{1}=0$ in (2.22) we obtain

[TABLE]

This formula is valid for any $c$ in the interval $(-\epsilon_{\textnormal{max}},\infty)$ and its form shows that the function $m_{\textnormal{eq}}(c)$ is real-analytic666A function is called real analytic at a point if it is infinitely differentiable at that point and it agrees with its Taylor series expansion around that point. at $c=0$. Therefore all the information about function $m_{\textnormal{eq}}(c)$ is contained in the value of this function and its derivatives at $c=0$.

Using limits (2.9) and (2.11) we can conclude that

[TABLE]

Let $\theta$ denote the following weighted combination of eigenvalues $\lambda_{2},\dots,\lambda_{n}$

[TABLE]

We now propose an approximate formula for $m_{\textnormal{eq}}(c)$

[TABLE]

Note this formula is much easier to use than (2.23) because it contains only one rational term. From (2.24) it is immediate that $\widehat{m}_{\textnormal{eq}}(c)$ also obeys the limits (2.9) and (2.11). Moreover it is straightforward to check that the first derivatives of $\widehat{m}_{\textnormal{eq}}(c)$ and $m_{\textnormal{eq}}(c)$ match at $c=0$. Hence the approximation error is given by a difference of second-order derivatives and we explain in Supplementary Material why this error is likely to be small. We can also view the approximation $\widehat{m}_{\textnormal{eq}}(c)$ of $m_{\textnormal{eq}}(c)$ as replacing a weighted arithmetic mean of several quantities with the corresponding harmonic mean. To see this note that from (2.23) and (2.24) we can express $m_{\textnormal{eq}}(c)$ as

[TABLE]

where

[TABLE]

is the weighted arithmetic mean of quantities $(1-c\lambda_{2})^{-1},\dots,(1-c\lambda_{n})^{-1}$ with weights $\alpha_{2},\dots,\alpha_{n}$777These weights may not be positive real numbers, as is customary in the definition of arithmetic means. However in our examples we generally find that the most significant weights indeed have a positive real part and a negligible imaginary part.. The corresponding weighted harmonic mean of these quantities is given by

[TABLE]

and observe that $\widehat{m}_{\textnormal{eq}}(c)$ can be expressed as the r.h.s. of (2.27) with arithmetic mean $\overline{x}$ replaced by the harmonic mean $\widehat{x}$.

We now illustrate the accuracy of this approximation using a couple of randomly generated $n\times n$, transition rate matrices $Q$ with $n=5$ and $n=10$ respectively. In both cases we choose the input rate to be $k_{\textnormal{in}}=1$ and the enzymatic state-values to be $\gamma_{i}=i$ for $1,2,\dots,n$. The exact function $\widehat{m}_{\textnormal{eq}}(c)$ along with its approximation $\widehat{m}_{\textnormal{eq}}(c)$ are plotted in Figure 3. The accuracy of this approximation can be easily seen. Notice that the exact function is slightly above its approximation. Assuming the significant weights ($\alpha_{i}$-s) are positive reals, this can be explained by the fact that arithmetic mean is always higher than the corresponding harmonic mean.

From (2.26) it is immediate that the shape of the function $\widehat{m}_{\textnormal{eq}}(c)$ depends crucially on the parameter $\theta$ computed according to (2.25). We now examine $\theta$ more closely and see how it is connected to an existing notion from the theory of Markov processes. Let us denote the numerator of (2.25) by

[TABLE]

Since $\alpha_{i}$-s are given by (2.21), using (2.19), $\widetilde{Q}=D^{-1}Q$ and $\lambda_{1}=0$ we can express $\Theta$ as

[TABLE]

This relation shows that $\Theta$ (and hence $\theta$) is always real-valued even though some $\lambda_{i}$-s or $\alpha_{i}$-s may have imaginary parts. Moreover to compute $\Theta$ we do not need to compute the eigenvalues $\lambda_{1},\dots,\lambda_{n}$ of a potentially large matrix $\widetilde{Q}$. Instead we only need to evaluate the expression $\pi^{T}D^{-1}QD^{-1}{\bf 1}$ which is computationally much easier. Interestingly the definition of $\Theta$ coincides with the well-known notion of Dirichlet forms, that is extensively used in the study of mixing properties of Markov processes [35, 23]. We now discuss this connection in more detail.

Consider the CTMC $(\gamma(t))_{t\geq 0}$ with state-space $\Gamma=\{\gamma_{1},\dots,\gamma_{n}\}$ and transition rate matrix $Q=[q_{ij}]$. The generator $\mathbb{Q}$888The generator of a Markov process is an operator which specifies the rate of change of the distribution of the process. For more details see Chapter 4 in [8]. of this CTMC maps any real-valued function $f$ on $\Gamma$ to another such real-valued function $\mathbb{Q}f$ given by

[TABLE]

Define a function $f:\Gamma\to(0,\infty)$ by $f(\gamma)=1/\gamma$. Then one can see that $\Theta$ (2.36) can be expressed as

[TABLE]

In other words, if a $\Gamma$-valued random variable $\gamma$ has distribution $\pi$, then $\Theta$ is the expectation of the random variable $(-f(\gamma)\mathbb{Q}f(\gamma))$. Relation (2.38) shows that $\Theta$ is a Dirichlet form associated with the Markovian semigroup generated by $\mathbb{Q}$ (see [35]). An important consequence of this connection is that $\Theta$ is always positive (see Lemma 2.1.2 in [35]) irrespective of the entries of the rate matrix $Q$ or the state values $\gamma_{1},\dots,\gamma_{n}$. The positivity of $\Theta$ implies that $\theta$ is also positive and hence the mapping $c\mapsto\widehat{m}_{\textnormal{eq}}(c)$ is convex and monotonically decreasing from $m^{\textnormal{(static)}}_{\textnormal{eq}}$ at $c=0$ to $m^{\textnormal{(det)}}_{\textnormal{eq}}$ as $c\to\infty$. Intuitively the magnitude of Dirichlet form $\Theta$ corresponds to the mixing strength of the underlying Markov process. Therefore as $\Theta$ increases, $\theta$ also increases and the mapping $c\mapsto\widehat{m}_{\textnormal{eq}}(c)$ has a sharper “drop” to the deterministic value $m^{\textnormal{(det)}}_{\textnormal{eq}}$. Our next goal is to make this mathematically precise and quantitatively estimate the relative speed-values $c$ beyond which the deterministic assumption is acceptable.

In the rest of this section, our object of interest will be the relative stochastic amplification factor defined by

[TABLE]

which measures the difference of steady-state substrate means in the presence and absence of enzymatic fluctuations, normalized by the the steady-state substrate mean in the deterministic case. Note that $\rho(c)$ does not depend on the input rate $k_{\textnormal{in}}$ and using (2.17) we can see that $\rho(c)$ satisfies

[TABLE]

In order to study the dependence of $\rho(c)$ on $c$, we now look at its approximation $\widehat{\rho}(c)$ which is defined analogously to (2.39), with $m_{\textnormal{eq}}(c)$ replaced by $\widehat{m}_{\textnormal{eq}}(c)$. Using (2.25), (2.36) and (2.24) we see that $\theta$ is the same as the normalized Dirichlet form defined by

[TABLE]

Substituting $\widehat{\lambda}$ by $\theta$ in (2.26) and dividing by $m^{\textnormal{(det)}}_{\textnormal{eq}}$, we obtain the following formula after some simple algebraic manipulations

[TABLE]

This formula clearly indicates that as $\theta$ gets larger, the amplification factor decreases more sharply to $1$ with the relative speed parameter $c$. One can regard $\rho(c)\approx\widehat{\rho}(c)$ as the “relative error” between the actual substrate mean and the mean computed with deterministic assumption on the enzymatic kinetics. From relation (2.42) it is immediate that in order to test if this error will exceed some tolerance level $\epsilon>0$ we just need to check if the relative enzyme speed $c$ is smaller than the threshold $c_{\epsilon}$ defined by

[TABLE]

We can expect this test to be rather conservative because as we have argued before, the exact values $m_{\textnormal{eq}}(c)$ will usually lie above their approximation $\widehat{m}_{\textnormal{eq}}(c)$.

Note that $c_{\epsilon}$ is inversely proportional to $\theta$ but directly proportional to $\rho_{\textnormal{max}}$. The first parameter $\theta$ is the normalized Dirichlet form and it captures the “mixing strength” of the underlying enzymatic dynamics (see Example 3.1), while the second parameter $\rho_{\textnormal{max}}$ can be viewed as a proxy for the variance of the stationary-distribution $\pi$999To see this note that $\rho_{\textnormal{max}}$ defined by (2.40) represents the “error” in Jensen’s inequality for the convex map $x\mapsto 1/x$. It can be easily shown that this error is proportional to the variance of the distribution $\pi$ (see [3] for instance).(see Example 3.1). Generally both these parameters will increase with higher levels noise in the enzymatic dynamics. However since they affect $c_{\epsilon}$ in opposing ways, it is difficult to ascertain the overall effect of dynamical noise in setting the threshold value $c_{\epsilon}$. We explore this issue in greater detail in Section 3.2 and numerically show that increasing levels of dynamical noise in the enzymatic kinetics of that reaction network gives rise to decreasing values of $c_{\epsilon}$. This is surprising and counterintuitive because it suggests that this dynamical noise is actually beneficial in improving the accuracy of the deterministic assumption for the enzyme activity.

Finally we remark that even though most of the analysis in this paper assumes that enzymatic kinetics is described by a finite Markov chain, the formulas we derive can provide insights for a more general class of stationary stochastic processes. This is because finite Markov chains can serve as good approximations of such processes [33]. Moreover if the process is Markov, even with an arbitrary state-space, we can compute expression (2.42) for $\widehat{\rho}(c)$ by sampling its stationary distribution and using this sample to estimate $\rho_{\textnormal{max}}$ and the normalized Dirichlet form $\theta$. We illustrate this for the example network in Section 3.2 where the enzyme dynamics follows a Markov process over a countable state-space.

From (2.26) it is immediate that the shape of the function $\widehat{m}_{\textnormal{eq}}(c)$ depends crucially on the parameter $\theta$ computed according to (2.25). We now examine $\theta$ more closely and see how it is connected to an existing notion from the theory of Markov processes. Let us denote the numerator of (2.25) by

[TABLE]

Since $\alpha_{i}$-s are given by (2.21), using (2.19), $\widetilde{Q}=D^{-1}Q$ and $\lambda_{1}=0$ we can express $\Theta$ as

[TABLE]

This relation shows that $\Theta$ (and hence $\theta$) is always real-valued even though some $\lambda_{i}$-s or $\alpha_{i}$-s may have imaginary parts. Moreover to compute $\Theta$ we do not need to compute the eigenvalues $\lambda_{1},\dots,\lambda_{n}$ of a potentially large matrix $\widetilde{Q}$. Instead we only need to evaluate the expression $\pi^{T}D^{-1}QD^{-1}{\bf 1}$ which is computationally much easier. Interestingly the definition of $\Theta$ coincides with the well-known notion of Dirichlet forms, that is extensively used in the study of mixing properties of Markov processes [35, 23]. We now discuss this connection in more detail.

Consider the CTMC $(\gamma(t))_{t\geq 0}$ with state-space $\Gamma=\{\gamma_{1},\dots,\gamma_{n}\}$ and transition rate matrix $Q=[q_{ij}]$. The generator $\mathbb{Q}$101010The generator of a Markov process is an operator which specifies the rate of change of the distribution of the process. For more details see Chapter 4 in [8]. of this CTMC maps any real-valued function $f$ on $\Gamma$ to another such real-valued function $\mathbb{Q}f$ given by

[TABLE]

Define a function $f:\Gamma\to(0,\infty)$ by $f(\gamma)=1/\gamma$. Then one can see that $\Theta$ (2.36) can be expressed as

[TABLE]

In other words, if a $\Gamma$-valued random variable $\gamma$ has distribution $\pi$, then $\Theta$ is the expectation of the random variable $(-f(\gamma)\mathbb{Q}f(\gamma))$. Relation (2.38) shows that $\Theta$ is a Dirichlet form associated with the Markovian semigroup generated by $\mathbb{Q}$ (see [35]). An important consequence of this connection is that $\Theta$ is always positive (see Lemma 2.1.2 in [35]) irrespective of the entries of the rate matrix $Q$ or the state values $\gamma_{1},\dots,\gamma_{n}$. The positivity of $\Theta$ implies that $\theta$ is also positive and hence the mapping $c\mapsto\widehat{m}_{\textnormal{eq}}(c)$ is convex and monotonically decreasing from $m^{\textnormal{(static)}}_{\textnormal{eq}}$ at $c=0$ to $m^{\textnormal{(det)}}_{\textnormal{eq}}$ as $c\to\infty$. Intuitively the magnitude of Dirichlet form $\Theta$ corresponds to the mixing strength of the underlying Markov process. Therefore as $\Theta$ increases, $\theta$ also increases and the mapping $c\mapsto\widehat{m}_{\textnormal{eq}}(c)$ has a sharper “drop” to the deterministic value $m^{\textnormal{(det)}}_{\textnormal{eq}}$. Our next goal is to make this mathematically precise and quantitatively estimate the relative speed-values $c$ beyond which the deterministic assumption is acceptable.

In the rest of this section, our object of interest will be the relative stochastic amplification factor defined by

[TABLE]

which measures the difference of steady-state substrate means in the presence and absence of enzymatic fluctuations, normalized by the the steady-state substrate mean in the deterministic case. Note that $\rho(c)$ does not depend on the input rate $k_{\textnormal{in}}$ and using (2.17) we can see that $\rho(c)$ satisfies

[TABLE]

In order to study the dependence of $\rho(c)$ on $c$, we now look at its approximation $\widehat{\rho}(c)$ which is defined analogously to (2.39), with $m_{\textnormal{eq}}(c)$ replaced by $\widehat{m}_{\textnormal{eq}}(c)$. Using (2.25), (2.36) and (2.24) we see that $\theta$ is the same as the normalized Dirichlet form defined by

[TABLE]

Dividing (2.26) by $m^{\textnormal{(det)}}_{\textnormal{eq}}$, we obtain the following formula after some simple algebraic manipulations

[TABLE]

This formula clearly indicates that as $\theta$ gets larger, the amplification factor decreases more sharply to $1$ with the relative speed parameter $c$. One can regard $\rho(c)\approx\widehat{\rho}(c)$ as the “relative error” between the actual substrate mean and the mean computed with deterministic assumption on the enzymatic kinetics. From relation (2.42) it is immediate that in order to test if this error will exceed some tolerance level $\epsilon>0$ we just need to check if the relative enzyme speed $c$ is smaller than the threshold $c_{\epsilon}$ defined by

[TABLE]

We can expect this test to be rather conservative because as we have argued before, the exact values $m_{\textnormal{eq}}(c)$ will usually lie above their approximation $\widehat{m}_{\textnormal{eq}}(c)$.

Note that $c_{\epsilon}$ is inversely proportional to $\theta$ but directly proportional to $\rho_{\textnormal{max}}$. The first parameter $\theta$ is the normalized Dirichlet form and it captures the “mixing strength” of the underlying enzymatic dynamics (see Example 3.1), while the second parameter $\rho_{\textnormal{max}}$ can be viewed as a proxy for the variance of the stationary-distribution $\pi$111111To see this note that $\rho_{\textnormal{max}}$ defined by (2.40) represents the “error” in Jensen’s inequality for the convex map $x\mapsto 1/x$. It can be easily shown that this error is proportional to the variance of the distribution $\pi$ (see [3] for instance).(see Example 3.1). Generally both these parameters will increase with higher levels noise in the enzymatic dynamics. However since they affect $c_{\epsilon}$ in opposing ways, it is difficult to ascertain the overall effect of dynamical noise in setting the threshold value $c_{\epsilon}$. We explore this issue in greater detail in Section 3.2 and numerically show that increasing levels of dynamical noise in the enzymatic kinetics of that reaction network gives rise to decreasing values of $c_{\epsilon}$. This is surprising and counterintuitive because it suggests that this dynamical noise is actually beneficial in improving the accuracy of the deterministic assumption for the enzyme activity.

Finally we remark that even though most of the analysis in this paper assumes that enzymatic kinetics is described by a finite Markov chain, the formulas we derive can provide insights for a more general class of stationary stochastic processes. This is because finite Markov chains can serve as good approximations of such processes [33]. Moreover if the process is Markov, even with an arbitrary state-space, we can compute expression (2.42) for $\widehat{\rho}(c)$ by sampling its stationary distribution and using this sample to estimate $\rho_{\textnormal{max}}$ and the normalized Dirichlet form $\theta$. We illustrate this for the example network in Section 3.2 where the enzyme dynamics follows a Markov process over a countable state-space.

3 Examples

In this section we present a couple of examples to illustrate our results. Our first example of a two-state switching enzyme is such that all calculations can be easily done analytically allowing us to clearly understand the stochastic amplification effect. We also see how enzymes can utilize their fluctuations to serve as high-gain amplifiers. Our second example is the reaction network of Paulsson et al. [32] which displays stochastic focusing. We apply our results to this network and demonstrate that in some cases dynamical fluctuations can actually be beneficial in reducing the unwanted stochastic amplification effects.

3.1 Stochastic amplification induced by a two-state switching enzyme

Consider a simple instance of the system in Figure 1, in which a single enzyme molecule is present and can fluctuate between two states of activity: with enzyme ${\bf E}$ in the low-activity (“0”) state, the degradation rate of substrate ${\bf S}$ is assumed to be $\gamma_{0}$, while it is equal to $\gamma_{1}$ when ${\bf E}$ is highly active (“1”). The 0-to-1 and 1-to-0 rates are given by $k_{\textnormal{on}}$ and $k_{\textnormal{off}}$ respectively. A schematic representation of this model is presented in Figure 4.

The time-varying degradation rate $(\gamma(t))_{t\geq 0}$ induced by this fluctuating enzyme ${\bf E}$ is a CTMC with state-space $\Gamma=\{\gamma_{0},\gamma_{1}\}$ and transition-rate matrix

[TABLE]

One can check that the unique stationary distribution $\pi=(\pi_{0},\pi_{1})$ for this CTMC is simply given by

[TABLE]

For each $i=0,1$, we can regard $\pi_{i}$ as the steady-state probability of the enzyme being in state $i$. Due to the Ergodic Theorem (see Theorem 10.6 in [18]) we can also view $\pi_{i}$ as the proportion of time that the enzyme spends in state $i$ in the long-run. Let $\gamma$ be a $\Gamma$-valued random variable with probability distribution $\pi$. Then its mean and variance can be computed as

[TABLE]

Suppose that the speed of enzymatic kinetics relative to the substrate is $c$ and so the degradation rate is given by the process $(\gamma_{c}(t))_{t\geq 0}$ defined by (2.2). Let $m_{\textnormal{eq}}(c)$ (2.5) be the steady-state substrate mean in this case. Using part (A) of Theorem 2.2 we obtain

[TABLE]

This formula involves the inverse of a $2\times 2$ matrix, which can be easily computed explicitly. Substituting this inverse along with the expressions for $\pi_{0}$ and $\pi_{1}$ (see (3.44)) we get

[TABLE]

It is interesting to point our the formal similarity of (3.51) with the formula for the mean transfer time of a relaxation process whose rate is modeled by a two-state CTMC [11]. From (3.51) it can be readily seen that the steady-state substrate means in the static and deterministic cases are given by

[TABLE]

Hence we can compute the maximum relative amplification factor $\rho_{\textnormal{max}}$ (2.40) as

[TABLE]

Recall the formula for $\textnormal{Var}_{\pi}(\gamma)$ from (3.45) and oberve that $\rho_{\textnormal{max}}$ can be expressed as

[TABLE]

which reinforces the point we made in Section 2.4 that $\rho_{\textnormal{max}}$ serves as a proxy for the variance of the stationary distribution.

The Dirichlet form $\Theta$ (2.38) for this CTMC is given by

[TABLE]

This yields the following formula for the normalized Dirichlet form $\theta$ (2.41)

[TABLE]

which determines the shape of our approximate formula $\widehat{m}_{\textnormal{eq}}(c)$ (2.26) for the steady-state substrate mean. It is straightforward to check that for this example, this approximate formula is exact because the expression (3.51) for $m_{\textnormal{eq}}(c)$ can be written as

[TABLE]

Note that $\theta$ is a measure of the mixing strength of the enzymatic kinetics and so it is not surprising that it increases linearly with the transition rates $k_{\textnormal{on}}$ and $k_{\textnormal{off}}$.

Define the relative amplification factor by (2.39). It can be exactly expressed as

[TABLE]

We now consider the situation when the degradation rate induced by the enzyme E in the low-activity state (“0”) is negligible. In this case $\gamma_{0}\approx 0$ and $\rho(c)$ simplifies to

[TABLE]

which shows that the relative amplification factor is proportional to $1/c$ and the proportionality constant is simply the product of the proportion of time ($\pi_{0}$) the enzyme spends in the low-activity state, the degradation rate ($\gamma_{1}$) at the high-activity state and the reciprocal of the sum of transition rates $k_{\textnormal{on}}$ and $k_{\textnormal{off}}$. In particular as $c$ approaches [math], the relative amplification factor $\rho(c)$ can be enormous, thereby indicating that such a switching enzyme ${\bf E}$ can exploit its fluctuations to function as a biological amplifier with a very high gain.

3.2 Stochastic focusing Network

In this section we apply our results to the famous stochastic focusing network given in [32]. This network involves three species: substrate S, product P and enzyme E121212In [32], S was called I and E was called S. We have changed the notation to ensure consistency with the notation in this paper.. The molecules of substrate S are produced constitutively at rate $k_{\textnormal{in}}$ and converted into product P through a first-order reaction with rate constant $k_{p}$. Both subtrate and product molecules degrade spontaneously at rates $k_{a}e$ and $\delta_{p}$ respectively, where $e$ denotes the current state or abundance level of enzyme E. The schematic representation of these reactions is as follows

[TABLE]

The enzymatic dynamics in this example is given by the Markovian birth-death process with birth-rate $k_{s}$ and death-rate $k_{d}$:

[TABLE]

This process evolves on state-space $\mathbb{N}_{0}$, which is the set of all nonnegative integers and its unique stationary distribution is Poisson with mean $k_{s}/k_{d}$. We assume that the initial enzymatic state is a random variable with this stationary distribution.

Multiplying the rate constants $k_{s}$ and $k_{d}$ by $c$, we obtain enzymatic kinetics whose speed relative to the substrate is $c$. Let $m_{\textnormal{eq}}^{({\bf S})}(c)$ and $m_{\textnormal{eq}}^{({\bf P})}(c)$ denote the steady-state means of substrate and product respectively when the relative enzyme speed is $c$. From the first-order moment equations for the network (3.58) one can easily show (see Supplementary Material) that for any $c\geq 0$

[TABLE]

To study the amplification of steady-state means due to enzymatic fluctuations we use the relative stochastic amplification factor defined by (2.39). Due to the linear relationship (3.60) between $m_{\textnormal{eq}}^{({\bf S})}(c)$ and $m_{\textnormal{eq}}^{({\bf P})}(c)$, these amplification factors are same for both product and subtrate. Therefore we can understand the amplification phenomenon by replacing network (3.58) with our simplified scheme (see Figure 1), where the degradation rate at time $t$ is given by

[TABLE]

and $E_{c}(t)$ denotes the state at time $t$ of enzymatic kinetics with relative speed $c$. Let $\gamma_{i}=k_{p}+k_{a}i$ for each $i=0,1,\dots$. Note that $(\gamma_{c}(t))_{t\geq 0}$ is a CTMC with state-space

[TABLE]

and stationary distribution131313This stationary distribution is obtained by applying the linear change of variables $\gamma=k_{p}+k_{a}e$ on the Poisson distribution with mean $k_{s}/k_{d}$.

[TABLE]

This CTMC transitions from state $\gamma_{i}$ to state $\gamma_{i+1}$ at rate $ck_{s}$ and from state $\gamma_{i}$ to state $\gamma_{i-1}$ at rate $cik_{d}$. In other words, the generator for this CTMC is given by

[TABLE]

for any bounded function $f:\Gamma\to\mathbb{R}$.

In the rest of this section, we denote the steady-state substrate mean by $m_{\textnormal{eq}}(c)$ instead of $m_{\textnormal{eq}}^{({\bf S})}(c)$. Since the state-space $\Gamma$ is not finite, we cannot use the results from Section 2.3 to compute $m_{\textnormal{eq}}(c)$. However we can easily simulate the paths of process $(\gamma_{c}(t))_{t\geq 0}$ with Gillespie’s Algorithm [9], and obtain samples of the random variable $\tau_{c}$ defined by (2.13). The corresponding sample mean then serves as an estimator for $m_{\textnormal{eq}}(c)$ (see part (C) of Proposition 2.1). Note that the steady-state substrate mean in the absence of enzymatic fluctations is simply given by

[TABLE]

Dividing $m_{\textnormal{eq}}(c)$ by $m_{\textnormal{eq}}^{\textnormal{(det)}}$ and subtracting $1$, we obtain an estimate for the relative stochastic amplification factor $\rho(c)$ (see (2.39)). This factor only depends on four rate constants $k_{p},k_{a},k_{s}$ and $k_{d}$ which we now set as

[TABLE]

We estimate $\rho(c)$ for several values of $c$ in the interval $(0,20)$ and plot these estimates in Figure 5. For each value of $c$, $\rho(c)$ was estimated using $10^{5}$ samples of $\tau_{c}$ and the resulting standard error141414The standard error is simply the standard deviation of the distribution of the sample mean. is also displayed in Figure 5. In Section 2.4 we develop an approximate expression $\widehat{\rho}(c)$ (2.42) for the relative amplification factor which is likely to hold even though the state-space $\Gamma$ is not finite. Using the stationary distribution $\pi$ (3.61) and the generator $\mathbb{Q}_{c}$ (with $c=1$) we estimate the maximum amplification factor $\rho_{\textnormal{max}}$ (2.40) and the normalized Dirichlet form $\theta$ (2.41) as $\rho_{\textnormal{max}}=0.1703$ and $\theta=2.1988$ respectively. With these values we evaluate the map $c\mapsto\widehat{\rho}(c)$ and plot it in the interval $(0,20)$ in Figure 5. The close agreement between the estimated and the approximate values of the relative amplification factor can be easily seen. In Figure 5 we also indicate the threshold speed $c_{\epsilon}$ (2.43) for the $1\%$ threshold level (i.e. $\epsilon=0.01$). This threshold speed is $c_{\epsilon}=7.2903$ which indicates that if $c<7.2903$ then enzymatic fluctuations will amplify the steady-state substrate mean by more than $1\%$ in comparison to the deterministic case. In other words, the relative error in assuming that the enzyme activity is deterministic exceeds $1\%$ if $c<7.2903$.

We now explore the effects of changing the levels of noise in the enzymatic activity, on the relative amplification factor for the steady-state substrate mean. This noise can be measured using the coefficient of variation (CV)151515The coefficient of variation of a probability distribution is its standard deviation divided by its mean. It measures the dispersion of a distribution relative to its mean. of the stationary distribution for the enzyme abundance. Since this distribution is Poisson with mean $k_{s}/k_{d}$, the CV is $(\sqrt{k_{d}}/\sqrt{k_{s}})$ which shows that for a fixed $k_{d}$, we can decrease the relative noise level by simply increasing $k_{s}$. With this in mind we repeat the above computations (see Figure 5) for three additional values of $k_{s}$: 5, 10 and 20, and the results are provided in Figure 6. For each value of $k_{s}$, the corresponding estimates for the maximum relative amplification factor $\rho_{\textnormal{max}}$, the normalized Dirichlet form $\theta$ and the threshold speed $c_{\epsilon}$ (for $\epsilon=0.01$) are given in Table 1.

Recall the discussion at the end of Section 2.4 on the effects of noise in the enzymatic dynamics. From Table 1 it is clear that as expected, decreasing noise (or increasing $k_{s}$) results in the decline of both $\rho_{\textnormal{max}}$ and $\theta$. These parameters influence the threshold speed $c_{\epsilon}$ (see (2.43)) in opposite ways, but their overall effect is to increase $c_{\epsilon}$, indicating that as the noise levels go down, the relative enzyme speed needs to be higher and higher for the assumption of deterministic enzymatic activity to be acceptable. In other words, even though noise in enzyme activity causes the stochastic amplification effect it also helps in eliminating it.

4 Discussion

We examined the mathematical properties of a system consisting of a substrate that is degraded through an enzyme with stochastically fluctuating activity levels. Our analysis focused on the effect of enzymatic fluctuations on the mean substrate abundance and its deviations from the deterministic model predictions. It should be pointed out that even if the substrate inflow rate is assumed to be an independent stationary stochastic process with mean $k_{\textnormal{in}}$, our results will not be affected.

Whereas a stochastically varying production rate would leave the mean substrate level unaffected and equal to that of the deterministic model, fluctuations in the removal rate of the substrate result in a system that behaves very differently in the stochastic and deterministic regimes due to the product term in the degradation rate of $S$. Our formulas help quantify this discrepancy and study its behavior as the speed of enzymatic fluctuations varies from zero to infinity. They also provide an interesting connection between the amplification effect and the mixing properties of the Markov process describing the enzymatic activity fluctuations, which allow us to determine the speed above which this amplification becomes negligible for a given system parametrization. Note that the study of such systems through the use of approximate stochastic models such as the Linear Noise Approximation [5] is particularly challenging, since these methods typically fail to capture the very strong negative correlations between enzyme activity and substrate that can arise at slow enzyme fluctuations (see also discussion in [28]). On the contrary, the results presented here are valid under much milder simplifying assumptions, and can thus accurately reveal the magnitude of the discrepancy between stochastic and deterministic descriptions of the system.

Given the prevalence of enzymatic interactions in cell biology, stochastic fluctuations in enzyme activity and/or abundance are expected to play a large role in shaping the mean intracellular abundances of substrates [17], which could potentially also deviate significantly from the deterministically predicted amounts. Since many enzymes are allosterically regulated [24] by their products, substrates or other small signaling molecules, it would be very interesting to also study the effects of this regulation on the statistics of substrates and products, and examine potential noise reduction [28] or signal amplification strategies. As the sensitivity of single-molecule enzymology experimental techniques increases, it may soon be possible to study the phenomena described theoretically in this work within living cells.

Supplementary Material

Appendix S1 The model

In this paper we consider a system where the substrate ${\bf S}$ enters at a constant rate $k_{\textnormal{in}}$ and is degraded at a rate that depends on the activity state or abundance level of an enzyme ${\bf E}$. This activity state is assumed to fluctuate in time $t$ according to a continuous-time Markov chain (CTMC) $(\gamma(t))_{t\geq 0}$ over a finite state-space $\Gamma=\{\gamma_{1},\dots,\gamma_{n}\}$. The system can be written as

[TABLE]

The CTMC $(\gamma(t))_{t\geq 0}$ is described by its $n\times n$ transition rate matrix $Q=[q_{ij}]$ (see [30]). For any distinct $i,j\in\{1,2,\dots,n\}$, $q_{ij}\geq 0$ denotes the rate at which the process leaves state $\gamma_{i}$ and enters state $\gamma_{j}$. The diagonal entries of $Q$ are given by $q_{ii}=-\sum_{j\neq i}q_{ij}$. From now on we assume that the rate matrix $Q$ is irreducible161616A matrix $Q$ is called irreducible if there does not exist a permutation matrix $P$ such that the matrix $PQP^{-1}$ is block upper-triangular. which implies that there exists a unique stationary distribution $\pi=(\pi_{1},\dots,\pi_{n})\in\mathbb{R}^{n}_{+}$ satisfying

[TABLE]

where ${\bf 0}$ and ${\bf 1}$ denote the $n\times 1$ vectors of all zeroes and ones respectively. Since the state-space is finite and the transition rate matrix $Q$ is irreducible, the CTMC $(\gamma(t))_{t\geq 0}$ is ergodic which means that the probability distribution of $\gamma(t)$ converges to the stationary distribution $\pi$ as $t\to\infty$. As we are interested in the steady-state limit, without loss of generality we can assume that the initial state $\gamma(0)$ is distributed according to $\pi$, i.e. $\mathbb{P}(\gamma(0)=\gamma_{i})$ for each $i=1,\dots,n$. This ensures that the process $(\gamma(t))_{t\geq 0}$ is a stationary stochastic process whose finite-dimensional distributions are invariant under time-shifts. This means that for any finite collection of time-points $t_{1},t_{2},\dots,t_{n}$ the joint distribution of the random vector $(\gamma(t_{1}+s),\dots,\gamma(t_{n}+s))$ remains the same for all $s\geq 0$. This also implies that various statistical properties of this process do not depend on time. In particular its mean $\mathbb{E}(\gamma(t))$ is equal to

[TABLE]

where $\gamma$ is a $\Gamma$-valued random variable with probability distribution $\pi$ and $\mathbb{E}_{\pi}(\cdot)$ denotes the expectation w.r.t. this distribution.

In what follows, we need to view enzymatic dynamics at the timescale of substrate kinetics. For this we define a family of processes $(\gamma_{c}(t))_{t\geq 0}$ parameterised by the “relative speed” parameter $c$ as

[TABLE]

Like $(\gamma(t))_{t\geq 0}$, the process $(\gamma_{c}(t))_{t\geq 0}$ is also a CTMC over state-space $\Gamma=\{\gamma_{1},\dots,\gamma_{n}\}$ with transition rate matrix $Q_{c}=cQ$ and initial distribution $\pi$. Since $(\gamma(t))_{t\geq 0}$ is stationary, this process is also stationary with the same mean given by $\mathbb{E}_{\pi}(\gamma)=\mathbb{E}(\gamma_{c}(t))$ for all times $t\geq 0$. Replacing $(\gamma(t))_{t\geq 0}$ by $(\gamma_{c}(t))_{t\geq 0}$ in (S1), we will study how the steady-state mean of substrate abundance depends on the fluctuation speed $c$.

Given a sample path of the enzyme dynamics $(\gamma_{c}(t))_{t\geq 0}$ with relative speed $c$, we regard the dynamics of substrate molecular counts as a jump Markov chain $(S_{c}(t))_{t\geq 0}$ over the set of nonnegative integers $\mathbb{N}_{0}=\{0,1,2,\dots\}$. This Markov chain can be written in the random time change representation [8] as

[TABLE]

where $Y_{1}$ and $Y_{2}$ are independent, unit rate Poisson processes. Here the Poisson processes $Y_{1}$ and $Y_{2}$ capture the intermittency in the firing of production and degradation reactions. This intermittency becomes unimportant if the substrate is present in high copy-numbers [22] and in this case one can regard $(S_{c}(t))_{t\geq 0}$ as the dynamics of substrate concentration171717The concentration of any species is its copy-number divided by the system volume., specified by the following ODE

[TABLE]

Let $m_{c}(t)=\mathbb{E}(S_{c}(t))$ for each $t\geq 0$, where $(S_{c}(t))_{t\geq 0}$ evolves according to either (S4) or (S5). Our goal in this paper is to understand the role of fluctuations in the catalytic activity of enzyme ${\bf E}$ in determining the steady-state value of the mean

[TABLE]

Appendix S2 Expressions for $m_{\textnormal{eq}}(c)$: The general case

In this section we prove Proposition 2.1 in the main paper. For convenience, we restate this proposition below.

Proposition S2.1

Suppose $(\gamma(t))_{t\geq 0}$ is a real-valued stationary stochastic process with stationary distribution $\pi$ and state-space $\Gamma$ satisfying

[TABLE]

for some $\epsilon>0$. Let $(\gamma_{c}(t))_{t\geq 0}$ be the speed $c$ version of this process given by (S3) and define the substrate dynamics $(S_{c}(t))_{t\geq 0}$ either by (S4) or by (S5). Let $m_{c}(t)=\mathbb{E}(S_{c}(t))$ and let the steady-state limit $m_{\textnormal{eq}}(c)$ be given by (S6). Then we have the following:

• (A)

The value $m_{\textnormal{eq}}(c)$ is well-defined (i.e. the limit in (S6) exists) and it is given by

[TABLE]

• (B)

Let $\tau_{c}$ is the random variable defined by

[TABLE]

where $u$ is an independent random variable with the uniform distribution on $[0,1]$. Then we have

[TABLE]

• (C)

The limits below are satisfied as $c\to\infty$ and $c\to 0$ respectively:

[TABLE]

Proof. Let $\{\mathcal{F}_{t}\}$ be the filtration generated by the process $(\gamma_{c}(t))_{t\geq 0}$ and let $\mathcal{F}_{\infty}=\lim_{t\to\infty}\mathcal{F}_{t}$ be its limiting value. Given the information $\mathcal{F}_{\infty}$, the random path of the process $(\gamma_{c}(t))_{t\geq 0}$ is completely known. Hence we can formulate the ODE for the conditional first-moment $\mathbb{E}(S_{c}(t)|\mathcal{F}_{\infty})$ as follows

[TABLE]

This equation remains unchanged whether we use representation (S4) or (S5) for the substrate dynamics $(S_{c}(t))_{t\geq 0}$. Using $\exp(\int_{0}^{t}\gamma_{c}(s)ds)$ as the integrating factor we can write (S13) as

[TABLE]

Finally, integrating both sides w.r.t. time $t$ we obtain

[TABLE]

Taking expectations we get

[TABLE]

where the last equation follows from the Fubini’s theorem. Since the process $(\gamma_{c}(t))_{t\geq 0}$ is stationary, the distribution of the random variable $\int_{s}^{t}\gamma_{c}(u)du$ is same as the distribution of $\int_{0}^{t-s}\gamma_{c}(u)du$. Hence using a simple change of variables we can write

[TABLE]

which gives us the following formula for $m_{c}(t)$

[TABLE]

Using (S7) we see that

[TABLE]

and using the dominated convergence theorem we can conclude that

[TABLE]

Taking the limit $t\to\infty$ in (S15) we obtain the formula (S8) for $m_{\textnormal{eq}}(c)$. This finishes the proof of part (A) of the proposition.

Let $\tau_{c}$ be the random variable defined by (S9). Since it is a continuous random variable with range $[0,\infty)$ we can write

[TABLE]

Let $\{\mathcal{F}_{t}\}$ be the filtration as defined before. As $u$ in (S9) is a uniform $[0,1]$ random variable independent of the process $(\gamma_{c}(t))_{t\geq 0}$ we have

[TABLE]

Taking expectations both sides we get

[TABLE]

Substituting this in (S16) we obtain

[TABLE]

Replacing this integral in (S8) by $\mathbb{E}(\tau_{c})$ proves formula (S10). This finishes the proof of part (B) of the proposition. The proof of part (C) is already outlined in the main paper. $\Box$

Appendix S3 Expressions for $m_{\textnormal{eq}}(c)$: The finite CTMC case

In this section we prove Theorem 2.2 in the main paper. For convenience, we restate this theorem below.

Theorem S3.1

Suppose $(\gamma(t))_{t\geq 0}$ is a stationary CTMC with transition rate matrix $Q$, stationary distribution $\pi$ and state-space $\Gamma=\{\gamma_{1},\dots,\gamma_{n}\}$ (see Section S1). Let $(\gamma_{c}(t))_{t\geq 0}$ be the speed $c$ version of this process given by (S3) and define the substrate dynamics $(S_{c}(t))_{t\geq 0}$ either by (S4) or by (S5). Let the steady-state substrate mean $m_{\textnormal{eq}}(c)$ be given by (S6) and the diagonal matrix $D$ be defined by

[TABLE]

Then we have the following:

• (A)

The matrix $(D-cQ)$ is invertible and $m_{\textnormal{eq}}(c)$ can be expressed as

[TABLE]

• (B)

Suppose the matrix $\widetilde{Q}=D^{-1}Q$ is diagonalizable and let $\lambda_{1},\dots,\lambda_{n}$ be its eigenvalues. For each $i=1,\dots,n$ define $\alpha_{i}$ by

[TABLE]

Then $m_{\textnormal{eq}}(c)$ can be expressed as

[TABLE]

• (C)

The following relation is satisfied for any $c\geq 0$

[TABLE]

Proof. For each $i=1,\dots,n$ and $t\geq 0$ define

[TABLE]

Let $\beta(t)$ denote the vector $\beta(t)=(\beta_{1}(t),\dots,\beta_{n}(t))$. Note that

[TABLE]

Note that $\beta(0)=\pi$ because for each $i$ we have $\beta_{i}(0)=\mathbb{E}({\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}_{\{\gamma_{c}(0)=s_{i}\}})=\mathbb{P}(\gamma_{c}(0)=s_{i})=\pi_{i}$. From Proposition 4.1 in [14] we see that $\beta(t)$ satisfies the following ODE:

[TABLE]

Since $\beta(0)=\pi$, the solution to this ODE is

[TABLE]

Therefore

[TABLE]

On integrating from $t=0$ to $t=\infty$ we get

[TABLE]

This relation along with (S8) proves (S18). Since the matrix $Q$ is irreducible and Metzler with non-positive eigenvalues, the matrix $(D-cQ)$ is invertible and its inverse $(D-cQ)^{-1}$ is a matrix with only nonnegative real entries (see Theorem 2.6 in [38]).

We now provide another proof of (S18) using the Methods of Conditional Moments (MCM) approach of Hasenauer et al. [15]. We thus define $m_{i}(t):=\mathbb{E}[S(t)|\gamma(t)=\gamma_{i}]$ and consider the algebraic equations that describe the system of conditional moments for $m_{1},\dots,m_{n}$:

[TABLE]

If we define $M_{i}=\pi_{i}m_{i}$ and consider the limit as $t\to\infty$ in (S22), we get the following system of linear equations:

[TABLE]

Letting $M:=\begin{bmatrix}M_{1}&\dots&M_{n}\end{bmatrix}$, we get

[TABLE]

Therefore,

[TABLE]

The proof of part (B) is outlined in the main paper. For part (C) note that one of the inequalities $m^{\textnormal{(det)}}_{\textnormal{eq}}\leq m_{\textnormal{eq}}(c)$ was already shown in the main paper. Hence it suffices to prove that $m_{\textnormal{eq}}\leq m^{\textnormal{(static)}}_{\textnormal{eq}}$ which is equivalent to

[TABLE]

Let $\Pi$ be the $m\times m$ diagonal matrix with entries $\pi_{1},\dots,\pi_{m}$. Define another $n\times n$ matrix by

[TABLE]

Since $\Pi D^{-1}$ is a positive diagonal matrix and $\Pi(D-Q)^{-1}$ is a componentwise nonnegative matrix, we can conclude that $M$ is also a Metzler matrix. In order to prove (S24) we just need to show that

[TABLE]

which is equivalent to proving that

[TABLE]

Note that $(M+M^{T})$ is a symmetric Metzler matrix. Relation (S26) will hold if this matrix is nonpositive definite, which is same as saying that all its eigenvalues have nonpositive real parts. Theorem 2.6 in [38] shows that this matrix is nonpositive definite if we can find a componentwise positive vector $v$ such that

[TABLE]

Let $v=(\gamma_{1},\dots,\gamma_{n})$. Note that since $Q{\bf 1}={\bf 0}$ we have $(D-cQ){\bf 1}=D{\bf 1}=v$ and so $(D-cQ)^{-1}v={\bf 1}$. Therefore

[TABLE]

Similarly since $Q^{T}\pi={\bf 0}$ we have $(D-cQ^{T})\pi=D\pi$ and so $\pi=(D-cQ^{T})^{-1}D\pi=(D-cQ^{T})^{-1}\Pi v$. Hence

[TABLE]

Combining (S28) and (S29) proves (S27) and shows that $(M+M^{T})$ is a symmetric nonpositive definite matrix. Therefore (S26) holds and this finishes the proof of part (C) of this theorem. $\Box$

Appendix S4 Approximate formula for $m_{\textnormal{eq}}(c)$

In this section we examine the approximate formula for $m_{\textnormal{eq}}(c)$ and discuss why the approximation error is likely to be small. Recall from Theorem S3.1 that $\lambda_{1},\dots,\lambda_{n}$ are the eigenvalues of matrix $\widetilde{Q}$. Among these $\lambda_{1}=0$ while the eigenvalues $\lambda_{2},\dots,\lambda_{n}$ have negative real parts. Using part (B) of Theorem S3.1 gives us this exact formula for $m_{\textnormal{eq}}(c)$:

[TABLE]

From limits (S11) and (S12) we can conclude that

[TABLE]

Let $\theta$ denote the following weighted combination of eigenvalues $\lambda_{2},\dots,\lambda_{n}$

[TABLE]

Recall from the main paper that $\theta$ is the normalized Dirichlet form which is always positive. We proposed the following approximate formula for $m_{\textnormal{eq}}(c)$

[TABLE]

We define the relative error between the exact value $m_{\textnormal{eq}}(c)$ and its approximation $\widehat{m}_{\textnormal{eq}}(c)$ by

[TABLE]

Due to limits (S11) and (S12), we know that this error function is contained between [math] and $1$. Our goal is to argue that this relative error function has a relatively small magnitude. For this purpose we define a change of variables as

[TABLE]

Note that as $c$ goes from [math] to $\infty$, $x$ goes from [math] to $1$. Define a function $f$ by

[TABLE]

where

[TABLE]

for each $i=2,\dots,n$. Replacing $c$ by $x/((1-x)\theta)$ in (S30) and using (S31) we get

[TABLE]

Similarly $\widehat{m}_{\textnormal{eq}}(c)$ can be written as

[TABLE]

which allows us to express the relative error as

[TABLE]

Note that

[TABLE]

We can compute the first and second order derivatives of function $f(x)$ as

[TABLE]

Therefore from (S35) we obtain

[TABLE]

We can see that $\kappa=f^{\prime\prime}(0)/2$ measures the weighted “spread” of the eigenvalues $\lambda_{2},\dots,\lambda_{n}$ around $-\theta$, where the weights are given by $\beta_{2},\dots,\beta_{n}$. Numerical experiments indicate that generally only a few of these weights are significant while the others are negligibly small. This is precisely the situation where this spread is small since $\theta=-\sum_{i=2}^{n}\lambda_{i}\beta_{i}$. Hence we can safely assume that $\kappa$ is small. The function $f(x)$ is real-analytic at $x=0$ and its Taylor series expansion is given by

[TABLE]

Using (S34) we can conclude that the relative error behaves like $\kappa x^{2}$ which is small since $\kappa$ is small.

Appendix S5 Stochastic Focusing example

Consider the stochastic focusing network given in [32]. It involves three species: substrate S, product P and enzyme E. The molecules of substrate S are produced constitutively at rate $k_{\textnormal{in}}$ and converted into product P with rate constant $k_{p}$. Both subtrate and product molecules degrade spontaneously at rates $k_{a}e$ and $\delta_{p}$ respectively, where $e$ denotes the current state or abundance level of enzyme E. These reactions can be expressed as

[TABLE]

Let $m^{({\bf S})}(t)$ and $m^{({\bf P})}(t)$ denote the expected abundance level at time $t$ of substrate and product molecules respectively. Furthermore assume that the limit

[TABLE]

exists. From the reaction network (S36) it is immediate that

[TABLE]

Solving this ODE we get

[TABLE]

Therefore using limit (S37) we can conclude that

[TABLE]