We here provide full details of the calculations reported in the Results section. The system size-expansion which is at the heart of the analysis has to-date not been applied extensively to biological problems and thus we go into some detail in its elucidation in Sub section I, which is dedicated exclusively to Model I. For other recent applications of the general method in the context of reaction kinetics, see for example [28] and [29]. Subsections II and III (treating Model II and Model III, respectively) naturally build on the results of the first subsection and thus we only give the main steps of the calculations in these last two cases. Sub section IV has a brief discussion of the simulation methods used to verify the theoretical results.
Model I: Michaelis-Menten reaction occurring in a compartment volume of sub-micron dimensions. Substrate input into compartment is modeled as a Poisson process
The reaction scheme is
. The stochastic description of this system is encapsulated by the master equation which is a differential equation in the joint probability function π describing the system:
where π = π(n
C
, n
P
, n
S
), n
X
is the integer number of molecules of type X (where X = C, P, S), Ω is the compartment volume, and
are the step operators defined by their action on a general function g(n
X
) as:
g(n
X
) = g(n
X
± 1). Note that the relevant variables are three, not four: the integer number of molecules of free enzyme (n
E
) is not an independent variable due to the fact that the total amount of enzyme is conserved. The master equation cannot be solved exactly but it is possible to systematically approximate it by using an expansion in powers of the inverse square root of the volume of the compartments. This is generally called the system-size expansion [12].
The method proceeds as follows. The stochastic quantity, n
X
/Ω, fluctuates about the macroscopic concentrations [X]; these fluctuations are of the order of the square root of the number of particles:
Note that since n
E
+ n
C
= constant, it follows that n
E
= Ω[E] - Ω1/2ϵ
C
. The joint distribution function and the operators can now be written as functions of the new variables, ϵ
X
, giving: π = Π(ϵ
C
, ϵ
P
, ϵ
S
, t) and
; using these new variables the master equation Eq. (13) takes the form:
where
Note that in Eq. (18) terms which involve products of first and second-order derivatives, third-order derivatives or higher have been omitted - these do not affect the low-order moment equations which we will be calculating.
Analysis of Ω1/2 terms
The terms of order Ω1/2 are the dominant ones in the limit of large volumes. By equating both terms of this order on the right and left hand sides of Eq. (15) and using Eq. (16), one gets the deterministic rate equations:
These are exactly those which one would write down based on the classical approach whereby one ignores molecular discreteness and fluctuations. This is an important benchmark of the method since it shows that it gives the correct result in the limit of large volumes. On a more technical note, the cancelation of these two terms of order Ω1/2 makes Eq. (15) a proper expansion in powers of Ω-1/2. By imposing steady-state conditions we have the Michaelis-Menten (MM) equation:
where v
max
= k2 [E
T
] is the maximum reaction velocity, [E
T
] = [E] + [C] is the total enzyme concentration which is a constant at all times and K
M
= (k1 + k2)/k0 is the Michaelis-Menten constant.
Analysis of Ω0 terms
To this order, the master equation is a multivariate Fokker-Planck equation whose solution is Gaussian and thus fully determined by its first and second moments. The equations of motion for these moments can be straightforwardly obtained from the master equation to this order, leading to a set of coupled but solvable ordinary differential equations:
where,
Note that the matrices and vectors in the above equations have been reduced to a simpler form by the application of a few row operations. Note also that these equations are independent of ϵ
p
since the product-forming step is irreversible and hence the fluctuations in substrate and complex are necessarily decoupled from its fluctuations. At the steady-state, it is found that ⟨ϵ
S
,
C
⟩ → 0. From Eq. (14), it is clear that this implies that to this order the average number of substrate molecules per unit volume, ⟨n
S
/Ω⟩, is simply equal to the macroscopic concentration, [S]. The same applies for complex molecules. Hence to this order in the system-size expansion there cannot be any corrections to the macroscopic equations or to the MM equation. By writing the macroscopic concentrations in Eqs. (24) and (25) in terms of k
in
and solving, one obtains the variance and covariance of the fluctuations about the steady-state macroscopic concentrations. We here only give the result for the covariance since this will be central to our discussion later on:
where α = k
in
/v
max
is the normalized reaction velocity of the enzyme.
Analysis of Ω-1/2 terms
The system-size expansion is almost never carried out to this order because of the algebraic complexity typically involved, however it is crucial to find finite volume corrections to the deterministic rate equations and in particular to the MM equation. Using the master equation to this order, the first moment of the complex concentration is governed by the equation of motion:
Now the production of product P from complex occurs through a decay process which necessarily has to be described by a linear term of the form: k
in
= k2⟨n
C
/Ω⟩ (the steady-state condition). Since the steady-state macroscopic complex concentration is equal to [C] = k
in
/k2, then it follows that to any order in the expansion we have ⟨ϵ
C
⟩ = 0. This is always found to be the case in simulations as well. Hence it immediately follows from Eq. (27) that the average of fluctuations about the macroscopic substrate concentration are non-zero and given by:
From a physical point of view, this indicates that the steady-state concentration of substrate in the compartment is not equal to the value predicted by the MM equation (i.e. [S]) and hence the non-zero value of the average of the fluctuations about [S]. The real substrate concentration inside the compartment is obtained by substituting Eqs. (28) and (26) in Equation (14), leading to:
An alternative mesoscopic rate equation replacing the MM equation
The renormalization of the steady-state substrate concentration indicates the breakdown of the MM equation; this phenomenon occurs because of non-zero correlations between noise in the substrate and enzyme concentrations, ⟨ϵ
S
ϵ
C
⟩, which the MM equation implicity neglects. To obtain the alternative to the latter, one needs to obtain a relationship between the normalized reaction velocity, α and the real substrate concentration ⟨n
S
/Ω⟩; writing [S] in terms of α and substituting in Eq. (29), one obtains this new relation:
Note that in the limit of large volumes, the second term on the left hand side of Eq. (30) becomes vanishingly small and we are left with the MM equation. In the results section the quantity on the right hand side of Eq. (30) is referred to as α
M
since this is the normalized reaction velocity which would be predicted by the MM equation given the measured substrate concentration ⟨n
S
/Ω⟩ inside the compartment. A quick estimate of the magnitude of the error that one is bound to incur by using the conventional MM equation can be obtained by substituting α = 1/2 (i.e. enzyme is half saturated with substrate) in Eqs. (30) and (31), solving for α
M
and then using this value to compute the fractional error e = 1 - α
M
/α. This leads to the simple expression:
We finish this section by noting that Eq. (30) will be found to be valid generally and not only for the simple Michaelis-Menten scheme treated in this section; the details of the reaction network come in through the form of Eq. (31) which is reaction-specific.
Model II: Michaelis-Menten reaction occurring in a compartment volume of sub-micron dimensions. Substrate is input into compartment in groups or bursts of M molecules at a time
A natural generalization of Model I which has direct biological application is when substrate molecules are fed into the compartment M at a time with mean rate
. The total mean substrate input rate is then equal to k
in
= M
. The master equation for this process is Eq. (13) with the first term on the right hand side replaced by
. This leads to the following change in the expression for a2 (Eq. 17):
Note that since the expression for a1 (Eq. 16) is unchanged, the deterministic equations are precisely the same as those of Model I. However now the fluctuations about the macroscopic substrate concentration are enhanced by a factor M; consequently the entries in the vector B in Eq. (25) need the change k
in
→ k
in
M. The analysis proceeds in the same manner as before. The mesoscopic rate equation replacing the MM equation is now given by Eq. (30) together with:
The fractional error rate evaluated at α = 1/2 gives:
This clearly shows that generally larger deviations from the predictions of the MM equation are expected in this case compared to those computed for Model I.
Model III: Michaelis-Menten reaction with competitive inhibitor occurring in a compartment volume of sub-micron dimensions. Substrate input as in two previous models
Competitive inhibition is modeled by the set of reactions:
, where k4 =
[I] and [I] is the inhibitor concentration (similar models have been studied by Roussel and collaborators [30, 31] in the context of biochemical oscillators though these assume M = 1). In the rest of this section, we change the notation of enzyme-inhibitor complex from EI to V, just to make the math notation easier to read. The substrate input into the compartment is considered to occur as in Model II since this encapsulates that of Model I as well. The master equation for this system is:
The change of variables from n
X
to ϵ
X
is done as before, however note that now the conservation law for enzyme is different than in the two previous models. The total enzyme concentration is now equal to [E
T
] = [E] + [C] + [V] from which it follows that n
E
= Ω[E] - Ω1/2(ϵ
C
+ ϵ
V
). The description is chosen to be in terms of numbers of molecules of types C, S and V and thus E being a dependent variable does not show up explicitly in the step operators of the master equation above.
Due to the significant number of changes in the terms of the expansion from those of previous models, we will show the equivalent of Eqs. (15)-(18) in full. The master equation in the new variables ϵ
X
is given by:
where
Analysis of Ω1/2 terms
As for previous models, these terms give the macroscopic equations. Equating both terms of this order on the right and left hand sides of Eq. (37) and using Eq. (38), one obtains:
In the steady-state we have the Michaelis-Menten (MM) equation:
where β = [I]/K
i
and K
i
= k3/
is the dissociation constant of the inhibitor.
Analysis of Ω0 and Ω-1/2 terms
The equations for the first moments are easily obtained and we shall not reproduce them here; suffice to say that at steady-state, it is found that ⟨ϵ
S
,
C
,
V
⟩ → 0 which implies that to this order in the system-size expansion there cannot be any corrections to the macroscopic equations or to the MM equation. The addition of a new species, V, does however substantially increase the algebraic complexity in the equations of motion for the second moments computed using terms up to order Ω0. In particular the matrix A is now a 6 × 6 matrix, rather than the 3 × 3 matrix of the previous two models.
where,
and
In the above equations we have defined k' = k3 + k4. Note also that the system of equations has been simplified through the application of a few row operations.
Now to next order, i.e. Ω-1/2, the first moments of the concentrations of molecules of type C and V are governed by the equation of motions:
As in previous models, since the production of product P from complex occurs through a decay process, it follows that at steady-state, ⟨ϵ
C
⟩ = 0 which also implies ⟨ϵ
V
⟩ = 0 from Eq. (50). Hence it follows from Eq. (49) that ⟨ϵ
S
⟩ = [⟨ϵ
S
ϵ
C
⟩ + ⟨ϵ
S
ϵ
V
⟩]/Ω1/2 [E]. The two cross correlators can be estimated to order Ω0 by solving Eqs. (46)-(48). The non-zero value of ⟨ϵ
S
⟩ implies a renormalization of the substrate concentration inside the compartment and hence to a new rate equation replacing the MM equation. This is obtained exactly in the same manner as previously shown for Model I. The mesoscopic rate equation is found to be given by Eq. (30) together with:
where the numerator coefficients are given by:
and the denominator coefficients by:
Note that
= 0 such that at α = 0, there is no correction to the MM equation i.e. α
M
= 0 also. The case β = 0 reduces to Model II, i.e. f(α) is given by Eq. (34).
Stochastic simulation
In this section we briefly describe the simulation methods used to verify the theoretical results which are described in detail in the Results section. All simulations were carried out using Gillespie's exact stochastic simulation algorithm, conveniently implemented in the standard simulation platform, Dizzy [32].
The data points in Figure 2 were generated by iterating the following four-step procedure: (i) pick a value for α between 0 and 1. This gives the substrate input rate k
in
= α v
max
; (ii) run the simulation and measure the ensemble-averaged substrate concentration, ⟨n
S
/Ω⟩ = [S*] at long times; (iii) compute α
M
using the MM equation, α
M
= [S*]/([S*]+ K
M
); (iv) compute the absolute percentage error R
p
= 100|(1 - α
M
/α)|. The solid curves in Figure 2 were obtained by numerically solving the cubic polynomial in α given by Eqs. (7) and (8) in the Results section for given values of α
M
and then using the above expression for R
p
. Figure 3 is generated in the same manner as Figure 2, except that: in step (i) we fix M and pick a value for α between 0 and 1. Since k
in
= M
, the required simulation parameter is
= α v
max
/M; step (iv) is not needed. The solid curves were obtained by numerically solving the cubic polynomial in α given by Eqs. (7) and (9) in the Results section for given values of [S*]. The y-axis for this figure is v/v
max
= α
M
for the MM equation and v/v
max
= α for the stochastic model. Figure 4 is obtained by numerically solving the quintic polynomial in α given by Eqs. (7) and (12) in the Results section together with the coefficients given by Eqs. (52)-(62) in the present section; the inhibitor concentration, [I], is varied while the substrate concentration, [S*], is kept fixed. The substrate concentration is chosen so that at [I] = 0, v/v
max
= 0.909 in all cases. Note that for models I and II, α
M
= [S*]/([S*] + K
M
) while for Model III, α
M
= [S*]/([S*] + (1 + β)K
M
). Note that the error bars are very small on the scale of the figures and thus are not shown.