Investigating the robustness of the classical enzyme kinetic equations in small intracellular compartments
aa r X i v : . [ q - b i o . S C ] O c t Investigating the robustness of the classical enzyme kineticequations in small intracellular compartments
Ramon Grima
School of Biological SciencesCentre for Systems Biology at Edinburgh, University of Edinburgh, UK • Corresponding author: Ramon Grima, email: [email protected] bstract
Background : Classical descriptions of enzyme kinetics ignore the physical nature of the in-tracellular environment. Main implicit assumptions behind such approaches are that reactionsoccur in compartment volumes which are large enough so that molecular discreteness can be ig-nored and that molecular transport occurs via diffusion. Though these conditions are frequentlymet in laboratory conditions, they are not characteristic of the intracellular environment, whichis compartmentalized at the micron and submicron scales and in which active means of transportplay a significant role.
Results : Starting from a master equation description of enzyme reaction kinetics and assumingmetabolic steady-state conditions, we derive novel mesoscopic rate equations which take intoaccount (i) the intrinsic molecular noise due to the low copy number of molecules in intracel-lular compartments (ii) the physical nature of the substrate transport process, i.e. diffusionor vesicle-mediated transport. These equations replace the conventional macroscopic and de-terministic equations in the context of intracellular kinetics. The latter are recovered in thelimit of infinite compartment volumes. We find that deviations from the predictions of classicalkinetics are pronounced (hundreds of percent in the estimate for the reaction velocity) for en-zyme reactions occurring in compartments which are smaller than approximately 200nm, for thecase of substrate transport to the compartment being mediated principally by vesicle or granuletransport and in the presence of competitive enzyme inhibitors.
Conclusions : The derived mesoscopic rate equations describe subcellular enzyme reaction ki-netics, taking into account, for the first time, the simultaneous influence of both intrinsic noiseand the mode of transport. They clearly show the range of applicability of the conventional de-terministic equation models, namely intracellular conditions compatible with diffusive transportand simple enzyme mechanisms in several hundred nanometre-sized compartments. An activetransport mechanism coupled with large intrinsic noise in enzyme concentrations is shown to leadto huge deviations from the predictions of deterministic models. This has implications for thecommon approach of modeling large intracellular reaction networks using ordinary differential2quations and also for the calculation of the effective dosage of competitive inhibitor drugs.3 ackground
The inside of a cell is a highly complex environment. In the past two decades, detailed mea-surements of the chemical and biophysical properties of the cytoplasm have established thatthe conditions in which intracellular reactions occur are, by and large, very different than thosetypically maintained in laboratory conditions. One of the outstanding differences between invivo and in vitro conditions, is that in the former, biochemical reactions typically occur in mi-nuscule reaction volumes [1]. For example, in eukaryotic cells, many biochemical pathways aresequestered within membrane-bound compartments, ranging from ∼ −
400 mg/ml which implies that between 5% and 40% of the total intracellular volumeis physically occupied by these molecules [5]. The concentration of these crowding molecules ishighly heterogeneous (see for example [6]), meaning that typically one will find small pocketsof intracellular space, characterized by low macromolecular crowding, surrounded by a “sea”of high crowding; such pockets of space may serve as effective compartments where reactionsmay occur more easily than in the rest of the cytosol. Analysis of experimental data for thedependence of diffusion coefficients with molecular size suggests the length scale of such effectivecompartments is in the range 35-50nm [7], a size comparable to that of the smallest vesicles. Thesignificant crowding also suggests that frequently an active means of transport such as vesicle-mediated transport, may be more desirable than simple diffusion as a means of intracellulartransport.The volume of a spherical cavity of space of diameter 50nm is merely ∼ . × − liters,an extremely small number compared to the typical macroscopic reaction volumes of in vitro experiments (experimental attolitre biochemistry is still in its infancy - see for example [8]).These very small reaction volumes imply that at physiologically relevant concentrations (nanoto millimolar), the copy number of a significant number of intracellular molecules is very small[1] and consequently that intrinsic noise cannot be ignored; for example 255 µM corresponds toan average of just 10 molecules in a 50nm vesicle and fluctuations about this mean of the orderof 3 molecules [9]. 4he traditional mathematical framework of physical chemistry ignores the basic physicalproperties of the intracellular environment. Kinetics are described by a set of coupled ordinarydifferential equations which implicitly assume (i) that the reaction compartment is so largethat molecular discreteness can be ignored and that hence integer numbers of molecules perunit volume can be replaced by a continuous variable, the molar concentration. Since thenumber of molecules is assumed to be very large, stochastic fluctuations are deemed negligibleand the equations are hence deterministic; (ii) the reaction compartment is well-stirred so thathomogeneous conditions prevail throughout [9]. Both assumptions can be justified for reactionsoccurring in a constantly stirred reactor of macroscopic dimensions. However if diffusion is thedominant transport process inside the compartment then the homogeneity assumption holdsonly if the volume is small enough so that in the time between successive reactions, a moleculewill diffuse a distance much larger than the size of the compartment. This comes at the expenseof the first assumption. It hence appears natural that for intracellular applications, the firstassumption, namely that of deterministic kinetics cannot be justified a priori . The secondassumption can be justified if reactions are localized in sufficiently small parts of the cell and inparticular for reaction-limited processes i.e. those for which the typical time for two moleculesto meet each other via diffusion is much less than the typical time for them to react if they arein close proximity. For such conditions, a molecule will come within reaction range several timesbefore participating in a successful reaction, in the process sampling the compartment manytimes which naturally leads to well-mixed conditions [9, 10, 11].In this article we seek to understand the magnitude of deviations from the classical ki-netic equations in small intracellular compartment volumes. We specifically focus on the caseof reaction-limited enzyme reactions which allows us to relax the first assumption of physicalchemistry while keeping the second one; this makes the mathematics tractable. We quantifydeviations from classical kinetics in the context of the Michaelis-Menten (MM) equation; this isthe cornerstone of present day enzyme kinetics and is a derivative of the traditional determinis-tic mathematical framework based on ordinary differential equations. In steady-state metabolicconditions, it is predicted to be exact. Thus this equation is ideal as a means to accurately testthe effects of small-scale compartmentation on chemical kinetics. We consider three successive5iological models of intracellular enzyme kinetics, each one building on the biological detail andrealism from the previous one (Figure 1). The models incorporate the intrinsic noisiness ofkinetics in small compartments, the details of the substrate transport process to the compart-ment (diffusion or active transport) and the presence of intra-compartmental molecules otherthan substrate molecules which may modulate the enzyme-catalyzed process e.g. inhibitors. Onthe macroscopic level, i.e. for large volumes, the steady-state kinetics of all models conformto the MM equation. We test whether this equation holds on the on the length scale of smallintracellular compartments by deriving the dependence of the ensemble averaged rate of productformation on the ensemble-averaged substrate concentration from the corresponding stochasticmodels in the steady-state. It is shown via both calculation and stochastic simulation that atthese small length scales the MM equation breaks down, being replaced by a new more generalequation. Practical consequences of this breakdown are illustrated in the context of the calcu-lation of the effective dosage of enzyme inhibitor drug needed to suppress intra-compartmentalenzyme activity by a given amount. To make our approach accessible to readers not familiarwith stochastic equations and their analysis, in the Results/Discussion sections we mainly focuson the biological/biophysical context and implications of the models together with the mainmathematical results which are verified by simulation. Detailed mathematical derivations andthe methods of simulation are relegated to the Methods section. Results
Model I: Michaelis-Menten reaction occurring in a compartment volume ofsub-micron dimensions. Substrate input into compartment occurs via a Pois-son process
This is the simplest, biologically-relevant case (Figure 1(A)). The reaction scheme is k in −−→ S + E k − ⇀↽ − k C k −→ E + P . Substrate molecules ( S ) are continuously supplied inside the compartmentat some rate k in , they reversibly bind to enzyme molecules ( E ) with rate constants k (forwardreaction) and k (backward reaction) to form transitory enzyme-substrate complex molecules( C ) which then decay with rate k into enzyme and product molecules ( P ). The substrate inputis assumed to be governed by a Poisson process with mean k in ; this is consistent with substratetransport to the compartment being dominated by normal diffusion. The enzyme acts as a cat-6lyst, effectively speeding up the reaction by orders of magnitude. It is assumed that diffusioninside the compartment is normal and not rate-limiting on the catalytic process i.e. well-mixedconditions or rate-limited kinetics inside the compartment. Given these conditions we ask our-selves what is the relationship between the reaction velocity and the substrate concentrationinside the compartment. The simplest approach consists of writing down the rate equations oftraditional physical chemistry: [ E T ] = [ E ] + [ C ] = constant, (1) d [ S ] /dt = k in − k [ E ][ S ] + k [ C ] , (2) d [ C ] /dt = k [ E ][ S ] − ( k + k )[ C ] , (3) d [ P ] dt = k [ C ] . (4)By imposing steady-state conditions we get the sought-after relationship which is simply thewell-known MM equation: d [ P ] dt = v = v max [ S ] K M + [ S ] , (5)where K M = ( k + k ) /k is the MM constant, v max = k [ E T ] is the maximum reaction velocityand square brackets indicate the macroscopic concentrations. We note that steady-state condi-tions for substrate necessarily require that k in ≤ v max otherwise the substrate will continuouslyaccumulate with time. Though this approach is simple and straightforward, as mentioned in theintroduction, the assumptions behind the formulation of the rate equations are not consistentwith the known physical properties of the cytoplasm. In particular it is clear that if the volumeof our compartment is very small (as is the case), the numbers of particles is also quite small,meaning that the concept of a continuous variable such as the average macroscopic concentrationhas little meaning. Rather we require a mathematical description in terms of discrete, integernumbers of particles and which is stochastic. The natural description of such cases is a masterequation which is a differential equation in the joint probability function π describing the system[12, 13, 14, 15]: dπdt = k in Ω(Θ − S − π + k Ω (Θ S Θ − C − n S n E π + k (Θ C Θ − S − n C π + k (Θ C Θ − P − n C π, (6)7here π = π ( n S , n C , n P ), n Y is the integer number of molecules of type Y , Ω is the compartmentvolume, and Θ ± X are step operators defined in the Methods section. This equation cannot besolved exactly. However it can be solved perturbatively using the system-size expansion due tovan Kampen [12]. This expansion is one in powers of the inverse square root of the compartmentvolume. In the Methods section, we calculate the first three terms of the expansion, namelythose proportional to Ω / , Ω and Ω − / . The first term, being the dominant one for largevolumes, gives back as expected, the rate equations Eqs. (1)-(4). The second term gives themagnitude of stochastic fluctuations about the macroscopic concentrations. Corrections to therate equations and to the MM equation (due to small compartment volumes or equivalentlydue to intrinsic noise) are found by considering the third term. In the rest of the article,instead of using the reaction velocity v , we use the normalized reaction velocity, α , which issimply the velocity of the reaction, v , divided by the maximum reaction velocity, v max . Givensome measured intracompartmental substrate concentration, [ S ∗ ] = h n S / Ω i (angled bracketsimply average), the relationship between the normalized reaction velocity predicted by the MMequation ( α M = [ S ∗ ] / ( K M + [ S ∗ ])) and the actual normalized reaction velocity ( α ), as predictedby our theory, is given by: α + (1 − α M ) f ( α )Ω − = α M , (7)where, f ( α ) = α K M + [ E T ](1 − α ) . (8)Hence the prediction of the MM equation is only correct, i.e. α = α M , in the limit of infinitelylarge compartment volumes, in which case the second term on the left hand side of Eq. (7)will become vanishingly small and can be neglected. For finite compartment volumes, the MMequation is not exact (except in the two limiting cases of α M → α M →
1) but is at bestan approximation, even though steady-state conditions are imposed; this is at odds with theprediction of the conventional deterministic theory. An inspection of Eqs. (7) and (8) shows thatthe magnitude of the deviations from the MM equation depends on the two non-dimensionalquantities: (i) K M Ω, a measure of the rate at which enzyme-substrate combination events occurrelative to the rate of decay of complex molecules; (ii) [ E T ]Ω, the total integer number of enzymemolecules in the compartment. 8s shown in the Methods section, the MM equation is found to implicity assume that thenoise about the macroscopic substrate and enzyme concentrations is uncorrelated (this assump-tion has generally been found to be at the heart of many macroscopic models - for examplesee [16]); properly taking into account these non-zero correlations leads to the corrections en-capsulated by Eqs. (7) and (8). These correlations are expected to be small in two particularcases: (i) if K M is large; in this case when substrate molecules combine with an enzyme toform a complex, the latter dissociates very quickly back into free enzyme and thus successiveenzyme-substrate events to the same enzyme molecule are bound to be almost independent ofeach other. The opposite situation of small K M would imply that the bottleneck in the catalyticprocess is the decay of complex rather than enzyme-substrate combination; if a successful com-bination occurs, the next substrate to arrive to the same enzyme molecule would have to waituntil the complex decays, naturally leading to correlations between successive enzyme-substratecombination events. (ii) if the total number of enzyme molecules is large; in such a case, at anyone time, the noise about the macroscopic concentrations will be the sum total from a largenumber of enzymes, each at a different stage in the catalytic process and each independent fromall others, which naturally dilutes any temporal correlations.To estimate the magnitude of the deviations from the MM equation inside cells, we usethe above two equations, Eqs. (7) and (8), to compute the absolute percentage error R p =100 | − α M /α | . These estimates are also computed by stochastic simulation of the master Eq.(6), using the exact stochastic simulation algorithm of Gillespie [10] (see Methods for detailsregarding the method of simulation); this provides a direct test of the theory. Figure 2 showsthe typical dependence of R p on α M , as predicted by both theory (solid lines) and simulation(data points). Generally the agreement between the two is found to be very good; discrepanciesincrease as K M and compartment volume decrease but are small for parameter values realisticfor intracellular conditions. The maxima of such plots gives the maximum absolute percentageerror which is a measure of the maximum expected deviations from the MM equation. Table 1summarizes these estimates (theory and simulation) over wide ranges of the parameters typicalof in vivo conditions: K M = 10 µM − µM [17], enzyme copy numbers of ten and one hundredper compartment which correspond to enzyme concentrations ranging from 4 µM to 2.5 m M and9ompartment diameters ranging from 50nm to 200nm. Note that the maximum deviations fromthe MM equation are estimated to be less than approximately 20% and typically just a fewpercent over large ranges of parameter values – this robustness of the MM equation with respectto intrinsic molecular noise is indeed surprising, since strictly speaking it is only valid for infinitecompartment volumes.The theory is always found to underestimate the actual deviations predicted by simulations;hence the theoretical expressions provide a quick, convenient way by which one can generallyestimate a lower bound on the deviations to be expected from the MM equation without theneed to perform extensive stochastic simulation.
Model II: Michaelis-Menten reaction occurring in a compartment volume ofsub-micron dimensions. Substrate is input into compartment in groups orbursts of M molecules at a time
Model I captures the basics of a general enzyme-catalyzed process occurring in a small intracel-lular compartment. In this section we build upon this model to incorporate further biologicalrealism. In particular, in the previous model we assumed that substrate input can be well de-scribed by a Poisson process, where one molecule at a time is fed into the compartment withsome average rate k in . This is the simplest possible assumption and approximates well the sit-uation in which molecules are brought to the compartment via normal diffusion. However thereare many situations where this may not be the case; we now describe two such cases.The intracellular condition of macromolecular crowding limits the Brownian motion of moleculesin the cytoplasm, this being reflected in the relatively small diffusion coefficients measured invivo compared to those known in vitro for moderately to relatively large molecules. Experimentswith inert tracer particles in the cytoplasm of Swiss 3T3 cells show that the in vivo diffusioncoefficient is an order of magnitude less than that in vitro for molecules with hydrodynamicradius 14nm and diffusion becomes negligibly small for molecules larger than approximately25nm [7]; similar results have been obtained in Xenopus neurons [18] and skeletal muscle my-otubes [19]. If diffusion is considerably hindered, one expects active transport to become a moredesirable mode of transport. Indeed there exists ample evidence for the active transport ofmacromolecules: they are typically packaged in a vesicle or a granule which is then transported10long microtubules or by some other means. It is also found that each vesicle or granule typ-ically contains several of these molecules (examples are: mRNA molecules - several estimatedper granule [20, 21]; cholesterol molecules which are transported in low-density lipoproteins [2]- approximately 1500 per lipoprotein).Generally an active means of transport is not exclusively linked with the transport of largesubstrate molecules. The cell being a highly compartmentalized and dynamic entity requires forits survival the precise transport of certain molecules from one compartment to another and aregulation of this transport depending on its current physiological state. Brownian motion leadsto an isotropic movement of molecules down the concentration gradient and to a consequentdamping of the substrate concentration with distance. In contrast active transport provides adirected (anisotropic) means of transport with little or no loss of substrate with distance, isindependent of the concentration gradient and it is also easily amenable to modulation.Hence it follows that a more general process modeling molecular entry into an intracellularcompartment is one in which M molecules are fed into the compartment at a rate k in ; the latterrate constant is the rate at which vesicles or granules arrive to the site of the compartment(Figure 1(B)). The total mean substrate input rate is then k in = M k in . The special case M = 1 corresponds to Model I. We construct the relevant master equation and employ thesystem-size expansion as for the previous model (see Methods for details); it is found that thedeterministic rate equations are exactly Eqs. (1)-(4) i.e. at the macroscopic level, given tworeactions occurring in two different compartments, both with the same total mean substrateinput rate k in but one occurring via diffusion (e.g. M = 1 , k in = 1) and the other via activetransport (e.g. M = 10 , k in = 0 .
1) , cannot be distinguished. However if the compartmentvolumes become small, then once again we find corrections to the MM equation and interestinglythese corrections are sensitive to the mode of transport. The relationship between the normalizedreaction velocity predicted by the MM equation ( α M ) and the actual normalized reaction velocity( α ), as predicted by our theory, is given by Eq. (7) together with: f ( α ) = α [ α + ( M − K m + [ E T ](1 − α ) . (9)This suggests that generally deviations from the predictions of the MM equation increase withthe carrying capacity, M , of the vesicle or granule. To compare the effects of active transport11nd diffusion on the kinetics, we set M = 50 and adjusted k in so that in all cases, the total meansubstrate input rate for model II is equal to k in , the input rate of Model I (i.e. the two modelswould be indistinguishable from a macroscopic point of view). Using the same procedure as forModel I, we computed the maximum percentage error using Eqs. (7) and (9) and also fromsimulations. The results are summarized in Table 2. Notice that now the deviations from theMM equation are much larger than before, running into hundreds of percent rather than the tensas for Model I. Because of the increase in substrate fluctuations, the quantitative accuracy ofthe theory is now less than before; it generally fares very well for compartments with diameterslarger than ∼ K M larger than ∼ µM . Nevertheless in all cases theory doescorrectly predict a large increase in discrepancy between the reaction velocities given by thedeterministic MM equation and those from stochastic simulation compared to the case of ModelI. The intuitive reason behind these increases in discrepancy is that substrate which is input inbursts enhances correlations between successive enzyme-substrate events.The explicit dependence of the reaction velocity on substrate concentration is complex andgenerally requires the solution of the cubic polynomial encapsulated by Eqs. (7) and (9). How-ever for small substrate concentrations, the equations simplify to a simple linear equation: α = [ S ∗ ] (cid:18) K M (cid:20) M − K M + [ E T ]) (cid:21)(cid:19) − (10)Note that if the MM equation was correct, one would expect α = [ S ∗ ] /K M . Indeed Eq. (10)reduces to the latter prediction in the limit of large volumes. Note also that this renormalizationof the proportionality constant occurs only if the substrate input occurs in bursts, i.e. M >
Model III: Michaelis-Menten reaction with competitive inhibitor occurringin a compartment volume of sub-micron dimensions. Substrate input as inprevious models
In this last section, we further build on the previous two models by adding competitive inhibitorsto the intracellular compartment in which enzymes are localized. A competitive inhibitor, I ,is one which binds reversibly to the active site of the enzyme (forming a complex EI ), thuspreventing substrate molecules from binding to the enzyme and slowing down catalysis (Figure1(C)). In standard textbooks and in the literature, this is typically modeled by the set of reactions12see for example [22]): k in −−→ S + E k − ⇀↽ − k C k −→ E + P, E k − ⇀↽ − k EI , where k = k [ I ] and [ I ] is theinhibitor concentration. Note that it is implicitly assumed that inhibitor is in such abundancethat the second-order bimolecular reaction between inhibitor and enzyme can be replaced by apseudo first-order reaction with constant inhibitor concentration. We shall assume the same inour model. Substrate input into the compartment is considered to occur as in Model II sincethis encapsulates that of Model I as well. The deterministic model of this set of reactions leadsto a MM equation of the form: d [ P ] dt = v max [ S ] K m (1 + β ) + [ S ] , (11)where β = [ I ] /K i and K i = k /k is the dissociation constant of the inhibitor. The perturbativesolution of the master equation describing the system is now significantly more involved than inprevious models; the underlying reason for this is that the computation of the noise correlatorsto order Ω requires the inversion of a 6 × × f ( α ) = 1 + βK m [ E T ] P i =0 c i (1 − α ) i P i =0 d i (1 − α ) i , (12)where c i and d i are coefficients with a complex dependence on the various enzyme parameters(these are given in full in the Methods Section). Table 3 shows the maximum percentage errorcomputed using Eqs. (7) and (12) and also from simulations for the cases in which substrateinput occur a molecule at a time and in bursts of 50 at a time. The parameter values chosenin the simulations and calculations (see caption of Table 3) are typical for many enzymaticprocesses: the bimolecular rate coefficients span the range 10 − s − M − [22], the backwarddecay processes are in the middle of the range 10 − s − [22], the inhibitor concentration is tentimes larger than the total enzyme concentration (satisfying the implicit assumption that theinhibitor is in significantly larger concentration than free enzyme), and the intracompartmentalenzyme concentrations are in the range 4 − µM . The deviations from the MM equation inthis case are more severe than in the previous two models, this being due to non-zero correlationsbetween substrate and the complex EI in addition to the already present correlations betweensubstrate and complex C . Note that the agreement between theory and simulations is overall13etter than in previous models, even when the burst size is large, M = 50. As mentioned in thesection for Model I, discrepancies between theory and simulation are generally found to decreasewith increasing K M ; for the case of competitive inhibition, the effective K M of the reactionis significantly larger than that of the enzyme (see Eq. (11)), which explains the increasedagreement between theory and simulations for Model III compared to the previous two models.A significant number of drugs suppress a chain of biochemical reactions by reducing theactivity of key enzymes in the pathway via competitive inhibition [17]. The conventional methodto estimate the required concentrations of these inhibitors involves plotting the variation ofthe enzyme activity with inhibitor concentration, [ I ], using the MM equation; the substrateconcentration is kept fixed and is chosen so that at [ I ] = 0, the reaction velocity is close tothe maximum, v max . Since there are significant corrections to the MM equation when reactionsoccur in intracellular compartments, it is not clear how accurate are estimates of [ I ] basedupon it. Figure 4 compares the enzyme activity curve based on the MM equation with thetheoretical predictions for the corrected enzyme activity curves based on the mesoscopic rateequation embodied by Eqs. (7) and (12), for compartments of diameter 50nm and 100nm (inset)and for substrate input burst sizes of M = 20 and 50. The substrate concentration is chosen sothat at [ I ] = 0, v/v max = 0 .
909 in all cases. We find that generally as the burst size increases,the actual inhibitor concentration needed to suppress enzyme activity by a given amount islarger than that estimated by the MM equation; this discrepancy decreases with increasingcompartment volume. For the example in Figure 4, for the case in which substrate is input intothe compartment in bursts of M = 50, the actual inhibitor concentration needed to decrease theenzyme activity from 0 .
909 to 0 . iscussion and Conclusion In this last section we discuss some fine points regarding: (i)the assumptions behind the use ofmaster equations which throws light on the range of use of the derived mesoscopic equations,(ii) the use of the system-size expansion to perturbatively solve the master equation and (iii)the assumption of steady-state metabolic conditions. We conclude by placing our work in thecontext of previous recent studies of stochastic enzyme kinetics and discuss possible experimentsto verify some of the conclusions we have reached.We have implicitly assumed throughout the article that a single (global) master equationmodel suffices to capture the deviations from classical kinetics due to fluctuations in chemicalconcentrations inside a single subcellular compartment. As noted by Baras and Mansour [23],“the global master equation selects the very limited class of exceptionally large fluctuationsthat appear at the level of the entire system, disregarding important nonequilibrium featuresoriginated by local fluctuations.” Hence the results presented here necessarily underestimatethe possible deviations from classical kinetics, in particular the local fluctuations due to dif-fusion of molecules inside the compartment. These local fluctuations are typically small forreaction-limited processes (as in this article) but significant for diffusion-limited ones. To cap-ture them effectively, one would be required to spatially discretize the compartment into manysmall elements and describe the reaction-diffusion processes between these elements by means ofa multivariate master equation [12, 23]. The latter is known as a reaction-diffusion master equa-tion; typically it does not allow detailed analytical investigation as for a global master equationand one is limited to stochastic simulation. Use of the global master equation is also restrictedfor compartments which are not too small: in particular the linear dimensions of the compart-ment should be larger than the average distance traveled by a molecule before undergoing asuccessful reaction with another molecule i.e. the length scale is much larger than that inherentin molecular dynamics simulation [23].We have applied the systematic expansion due to van Kampen to perturbatively solve themaster equation. It is sometimes a priori assumed that because this expansion is about themacroscopic concentrations, it cannot give information regarding the stochastic kinetics of fewparticle / small volume systems. This is true if one restricts oneself to the expansion to order15 i.e. the linear-noise approximation; this is commonly the case found in the literature sincethe algebra becomes tedious if one considers more terms. However we note that as argued andshown by van Kampen himself [12], terms beyond the linear-noise approximation in the system-size expansion add terms to the fluctuations that are of order of a single particle relative tothe macroscopic quantities and are essential to understanding how fluctuations are affected bythe presence of non-linear terms in the macroscopic equation (substrate-enzyme binding in ourcase). In our theory we went beyond the linear-noise approximation. We find that the predictedtheoretical results are in reasonable agreement, in many cases (comparison of bold and italicvalues in Tables 1, 2 and 3), with stochastic simulations of just a few tens of enzyme moleculesin sub-micron compartments, which justifies our methodology.We have also imposed metabolic steady-state conditions inside the subcellular compartment.Technically this is convenient since in such a case one does not deal with complex transients.Also since under such conditions the MM equation is exact from a deterministic point of view,it provides a very useful reference point versus which to accurately compute deviations due tointrinsic noise. In reality one may not always have steady-state conditions inside cells, thisdepending strongly on the rate of substrate input relative to the maximum rate at which theenzyme can process substrate. Another possibility is that one is dealing with a batch reactioni.e. one in which a number of substrate molecules are transported at one go and just once tothe subcellular compartment (e.g. via vesicle-mediated transport) and the reaction proceedsthereafter without any further substrate replenishment. This latter scenario is compatible withthe presentation of the MM equation typical in standard physical chemistry textbooks. TheMM equation is then an approximation (not exact as in steady-state case) to the deterministickinetics, when substrate is present in much larger concentration than enzyme. This case iscurrently under investigation using the same perturbative framework used in this article.We note that this is not the first attempt to study stochastic enzyme kinetics. The bulk ofrecent studies [24, 25, 26, 27] have focused on understanding the kinetics of a Michaelis-Mententype reaction catalyzed by a single enzyme molecule. Deviations from classical kinetics werefound to be most pronounced when one takes into account substrate fluctuations [26]. Thesepioneering studies were restricted to a single-enzyme assisted reaction which reduces complexity16hereby making it ideal from a theoretical perspective; since the reaction is dependent on justa single enzyme molecule one also finds maximum deviations from deterministic kinetics. Inreality, it is unlikely to find just one enzyme molecule inside a subcellular compartment - asmentioned in the introduction a physiological concentration of just a few hundred micromolarwould correspond to few tens inside the typically smallest subcellular compartment. It is also thecase that diffusion may not always be the main means of substrate transport to the compartmentand that the reaction maybe more complex than the simple Michaelis-Menten type reaction ofthese previous studies. The present study fills in these gaps by using a systematic methodto derive approximate and relatively simple analytic expressions for mesoscopic rate equationsdescribing the kinetics of the general case of N enzyme molecules in a subcellular compartmentwith or without active transport of substrate and in the presence of enzyme inhibitors. Mostimportantly our approach shows the effects of intrinsic noise on the kinetics can be capturedvia effective ordinary differential equations. This enables quick estimation of the magnitude ofstochastic effects on reaction kinetics and thus gives insight into whether a model or parts ofa model should be designed to be stochastic or deterministic without the need for extensivestochastic simulation. In the present study, this approach enabled us to readily compute, forthe first time, the deviations from deterministic kinetics for a broad range of realistic in vivo parameter constants (Tables 1, 2 and 3), a task which would be considerably lengthy if one hadto rely solely on data obtained from ensemble-averaged stochastic simulations.We conclude by briefly discussing possible experiments which can verify the predictionsmade in this article. It is arguably not an easy task to perform the required experiments in real-time in a living cell. A viable alternative would consist of monitoring reaction kinetics insidesingle artificially-made vesicles. Pick et al [8] have shown that the addition of cytochalasin tomammalian cells induces them to extrude from their plasma membrane minuscule vesicles ofattolitre volume with fully functional cell surface receptors and also retaining cytosolic proteinsin their interior. The change in the intra-vesicular calcium ion concentration in response tosurface ligand binding was measured using fluorescence confocal microscopy (FCM). Since thevesicle sizes are of typical small sub-cellular compartment dimensions (1 attolitre correspondsto a spherical vesicle of approximate diameter 120nm) and FCM allows the measurement of the17oncentration of a fluorescent probe (via a calibration procedure), this experimental techniqueappears ideal to verify the predictions of Model I and of Model III for the case of diffusivesubstrate transport. Model II and Model III with vesicle-transport of substrate are probablymuch more challenging to verify since one then needs to construct the in vitro equivalent ofmicrotubules. This is within the scope of synthetic biology and may be a possibility in the nextfew years. Methods
We here provide full details of the calculations reported in the Results section. The systemsize-expansion which is at the heart of the analysis has to-date not been applied extensively tobiological problems and thus we go into some detail in its elucidation in Subsection I, whichis dedicated exclusively to Model I. For other recent applications of the general method in thecontext of reaction kinetics, see for example [28] and [29]. Subsections II and III (treating ModelII and Model III, respectively) naturally build on the results of the first subsection and thus weonly give the main steps of the calculations in these last two cases. Subsection IV has a briefdiscussion of the simulation methods used to verify the theoretical results.
Model I: Michaelis-Menten reaction occurring in a compartment volume ofsub-micron dimensions. Substrate input into compartment is modeled as aPoisson process
The reaction scheme is k in −−→ S + E k − ⇀↽ − k C k −→ E + P . The stochastic description of this systemis encapsulated by the master equation which is a differential equation in the joint probabilityfunction π describing the system: dπdt = k in Ω(Θ − S − π + k Ω (Θ S Θ − C − n S n E π (13)+ k (Θ C Θ − S − n C π + k (Θ C Θ − P − n C π, 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 Θ ± X are the step operators defined by their action on ageneral function g ( n X ) as: Θ ± X g ( n X ) = g ( n X ± n E ) is not an independent variable18ue to the fact that the total amount of enzyme is conserved. The master equation cannotbe solved exactly but it is possible to systematically approximate it by using an expansion inpowers of the inverse square root of the volume of the compartments. This is generally calledthe system-size expansion [12].The method proceeds as follows. The stochastic quantity, n X / Ω, fluctuates about the macro-scopic concentrations [X]; these fluctuations are of the order of the square root of the numberof particles: n X = Ω[ X ] + Ω / ǫ X . (14)Note that since n E + n C = constant , it follows that n E = Ω[ E ] − Ω / ǫ C . The joint distributionfunction and the operators can now be written as functions of the new variables, ǫ X , giving: π = Π( ǫ C , ǫ P , ǫ S , t ) and Θ ± X = 1 ± Ω − / ∂/∂ǫ X + Ω − ∂ /∂ǫ X + O (Ω − / ); using these newvariables the master equation Eq. (13) takes the form: ∂ Π ∂t − Ω / (cid:18) d [ C ] dt ∂ Π ∂ǫ C + d [ P ] dt ∂ Π ∂ǫ P + d [ S ] dt ∂ Π ∂ǫ S (cid:19) = Ω / a Π + Ω a Π + Ω − / a Π + O (Ω − )(15)where a = − ( k in + k [ C ] − k [ E ][ S ]) ∂∂ǫ S + (( k + k )[ C ] − k [ E ][ S ]) ∂∂ǫ C − k [ C ] ∂∂ǫ P , (16) a = 12 k in ∂ ∂ǫ S + 12 (cid:18) ∂∂ǫ S − ∂∂ǫ C (cid:19) ( k [ S ][ E ] + k [ C ]) + k (cid:20) ∂∂ǫ C − ∂∂ǫ P (cid:21) ǫ C + (cid:20) ∂∂ǫ S − ∂∂ǫ C (cid:21) [ k ( ǫ S [ E ] − ǫ C [ S ]) − k ǫ C ] + 12 k (cid:18) ∂∂ǫ P − ∂∂ǫ C (cid:19) [ C ] , (17) a = 12 (cid:18) ∂∂ǫ S − ∂∂ǫ C (cid:19) ( k ǫ S [ E ] − k ǫ C [ S ] + k ǫ C ) − k (cid:20) ∂∂ǫ S − ∂∂ǫ C (cid:21) ǫ S ǫ C + 12 k (cid:18) ∂∂ǫ P − ∂∂ǫ C (cid:19) ǫ C . (18)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 momentequations which we will be calculating. 19 nalysis of Ω / terms The terms of order Ω / are the dominant ones in the limit of large volumes. By equating bothterms of this order on the right and left hand sides of Eq. (15) and using Eq. (16), one gets thedeterministic rate equations: d [ S ] /dt = k in − k [ E ][ S ] + k [ C ] , (19) d [ C ] /dt = k [ E ][ S ] − ( k + k )[ C ] , (20) d [ P ] /dt = k [ C ] . (21)These are exactly those which one would write down based on the classical approach whereby oneignores molecular discreteness and fluctuations. This is an important benchmark of the methodsince it shows that it gives the correct result in the limit of large volumes. On a more technicalnote, the cancelation of these two terms of order Ω / makes Eq. (15) a proper expansion inpowers of Ω − / . By imposing steady-state conditions we have the Michaelis-Menten (MM)equation: d [ P ] dt = v max [ S ] K M + [ S ] , (22)where v max = k [ E T ] is the maximum reaction velocity, [ E T ] = [ E ] + [ C ] is the total enzymeconcentration which is a constant at all times and K M = ( k + k ) /k is the Michaelis-Mentenconstant. Analysis of Ω terms To this order, the master equation is a multivariate Fokker-Planck equation whose solution isGaussian and thus fully determined by its first and second moments. The equations of motionfor these moments can be straightforwardly obtained from the master equation to this order,leading to a set of coupled but solvable ordinary differential equations: ∂ t (cid:20) h ǫ S ih ǫ C i (cid:21) = (cid:18) − k [ E ] k + k [ S ] k [ E ] − k ( K M + [ S ]) (cid:19) (cid:20) h ǫ S ih ǫ C i (cid:21) (23) ∂ t h ǫ S ih ǫ C ih ǫ S ǫ C i = A · h ǫ S ih ǫ C ih ǫ S ǫ C i + B, (24)20here, A = − k [ E ] 0 2( k + k [ S ])0 − k ( K M + [ S ]) 2 k [ E ]0 − k − k , B = k in + k [ C ] + k [ S ][ E ] k ([ S ][ E ] + K M [ C ]) k in + k [ C ] . (25)Note that the matrices and vectors in the above equations have been reduced to a simpler formby the application of a few row operations. Note also that these equations are independentof ǫ p since the product-forming step is irreversible and hence the fluctuations in substrate andcomplex are necessarily decoupled from its fluctuations. At the steady-state, it is found that h ǫ S,C i →
0. From Eq. (14), it is clear that this implies that to this order the average number ofsubstrate molecules per unit volume, h n S / Ω i , is simply equal to the macroscopic concentration,[ S ]. The same applies for complex molecules. Hence to this order in the system-size expansionthere cannot be any corrections to the macroscopic equations or to the MM equation. By writingthe macroscopic concentrations in Eqs. (24) and (25) in terms of k in and solving, one obtains thevariance 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 lateron: h ǫ C ǫ S i = K M [ E T ] α K M + [ E T ](1 − α ) , (26)where α = k in /v max is the normalized reaction velocity of the enzyme. Analysis of Ω − / terms The system-size expansion is almost never carried out to this order because of the algebraiccomplexity typically involved, however it is crucial to find finite volume corrections to the de-terministic rate equations and in particular to the MM equation. Using the master equation tothis order, the first moment of the complex concentration is governed by the equation of motion: d h ǫ C i /dt = − k ([ S ] + K M ) h ǫ C i + k [ E ] h ǫ S i − k Ω − / h ǫ S ǫ C i . (27)Now the production of product P from complex occurs through a decay process which necessarilyhas to be described by a linear term of the form: k in = k h n C / Ω i (the steady-state condition).Since the steady-state macroscopic complex concentration is equal to [ C ] = k in /k , then itfollows that to any order in the expansion we have h ǫ C i = 0. This is always found to be the21ase in simulations as well. Hence it immediately follows from Eq. (27) that the average offluctuations about the macroscopic substrate concentration are non-zero and given by: h ǫ S i = h ǫ S ǫ C i Ω / [ E ] . (28)From a physical point of view, this indicates that the steady-state concentration of substrate inthe compartment is not equal to the value predicted by the MM equation (i.e. [S]) and hencethe non-zero value of the average of the fluctuations about [S]. The real substrate concentrationinside the compartment is obtained by substituting Eqs. (28) and (26) in Equation (14), leadingto: D n S Ω E = [ S ] + K M α (1 − α )[ K M + [ E T ](1 − α ) ]Ω (29) An alternative mesoscopic rate equation replacing the MM equation
The renormalization of the steady-state substrate concentration indicates the breakdown of theMM equation; this phenomenon occurs because of non-zero correlations between noise in thesubstrate and enzyme concentrations, h ǫ S ǫ C i , which the MM equation implicity neglects. Toobtain the alternative to the latter, one needs to obtain a relationship between the normalizedreaction velocity, α and the real substrate concentration h n S / Ω i ; writing [ S ] in terms of α andsubstituting in Eq. (29), one obtains this new relation: α + (cid:18) − h n S / Ω i K M + h n S / Ω i (cid:19) f ( α )Ω − = h n S / Ω i K M + h n S / Ω i , (30) f ( α ) = α K M + [ E T ](1 − α ) (31)Note that in the limit of large volumes, the second term on the left hand side of Eq. (30) becomesvanishingly small and we are left with the MM equation. In the results section the quantity onthe right hand side of Eq. (30) is referred to as α M since this is the normalized reaction velocitywhich would be predicted by the MM equation given the measured substrate concentration h n S / Ω i inside the compartment. A quick estimate of the magnitude of the error that oneis bound to incur by using the conventional MM equation can be obtained by substituting α = 1 / α M andthen using this value to compute the fractional error e = 1 − α M /α . This leads to the simple22xpression: e = [1 + Ω([ E T ] + 4 K M )] − (32)We finish this section by noting that Eq. (30) will be found to be valid generally and not only forthe simple Michaelis-Menten scheme treated in this section; the details of the reaction networkcome in through the form of Eq. (31) which is reaction-specific. Model II: Michaelis-Menten reaction occurring in a compartment volume ofsub-micron dimensions. Substrate is input into compartment in groups orbursts of M molecules at a time. A natural generalization of Model I which has direct biological application is when substratemolecules are fed into the compartment M at a time with mean rate k in . The total meansubstrate input rate is then equal to k in = M k in . The master equation for this process is Eq.(13) with the first term on the right hand side replaced by Ω(Θ − MS − k in . This leads to thefollowing change in the expression for a (Eq. 17):12 k in ∂ ∂ǫ S → k in M ∂ ∂ǫ S . (33)Note that since the expression for a (Eq. 16) is unchanged, the deterministic equations areprecisely the same as those of Model I. However now the fluctuations about the macroscopicsubstrate concentration are enhanced by a factor M ; consequently the entries in the vector Bin 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) togetherwith: f ( α ) = α [ α + ( M − K M + [ E T ](1 − α ) . (34)The fractional error rate evaluated at α = 1 / e = MM + Ω([ E T ] + 4 K M ) (35)This clearly shows that generally larger deviations from the predictions of the MM equation areexpected in this case compared to those computed for Model I.23 odel III: Michaelis-Menten reaction with competitive inhibitor occurring ina compartment volume of sub-micron dimensions. Substrate input as in twoprevious models. Competitive inhibition is modeled by the set of reactions: k in −−→ S + E k − ⇀↽ − k C k −→ E + P, E k − ⇀↽ − k EI ,where k = k [ I ] and [ I ] is the inhibitor concentration (similar models have been studied byRoussel 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 compartmentis considered to occur as in Model II since this encapsulates that of Model I as well. The masterequation for this system is: dπdt = k in Ω(Θ − MS − π + k Ω (Θ S Θ − C − n S n E π + k (Θ C Θ − S − n C π + k (Θ C Θ − P − n C π + k (Θ V − n V + k (Θ − V − n E . (36)The change of variables from n X to ǫ X is done as before, however note that now the conservationlaw for enzyme is different than in the two previous models. The total enzyme concentrationis now equal to [ E T ] = [ E ] + [ C ] + [ V ] from which it follows that n E = Ω[ E ] − Ω / ( ǫ 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 masterequation above.Due to the significant number of changes in the terms of the expansion from those of previousmodels, we will show the equivalent of Eqs. (15)-(18) in full. The master equation in the newvariables ǫ X is given by: ∂ Π ∂t − Ω / (cid:18) d [ C ] dt ∂ Π ∂ǫ C + d [ P ] dt ∂ Π ∂ǫ P + d [ S ] dt ∂ Π ∂ǫ S + d [ V ] dt ∂ Π ∂ǫ V (cid:19) =Ω / a Π + Ω a Π + Ω − / a Π + O (Ω − ) (37)where a = − ( k in + k [ C ] − k [ E ][ S ]) ∂∂ǫ S + (( k + k )[ C ] − k [ E ][ S ]) ∂∂ǫ C (38)+ ( k [ V ] − k [ E ]) ∂∂ǫ V − k [ C ] ∂∂ǫ P , = 12 k in M ∂ ∂ǫ S + 12 (cid:18) ∂∂ǫ S − ∂∂ǫ C (cid:19) ( k [ S ][ E ] + k [ C ]) + k (cid:20) ∂∂ǫ C − ∂∂ǫ P (cid:21) ǫ C + (cid:20) ∂∂ǫ S − ∂∂ǫ C (cid:21) [ k ( ǫ S [ E ] − ( ǫ C + ǫ V )[ S ]) − k ǫ C ] + k (cid:18) ∂∂ǫ V ǫ V + 12 ∂ ∂ǫ V [ V ] (cid:19) + k (cid:18)
12 [ E ] ∂ ∂ǫ V + ∂∂ǫ V ( ǫ C + ǫ V ) (cid:19) + 12 k (cid:18) ∂∂ǫ P − ∂∂ǫ C (cid:19) [ C ] , (39) a = 12 (cid:18) ∂∂ǫ S − ∂∂ǫ C (cid:19) ( k ǫ S [ E ] − k ( ǫ C + ǫ V )[ S ] + k ǫ C ) − k (cid:20) ∂∂ǫ S − ∂∂ǫ C (cid:21) ǫ S ( ǫ C + ǫ V )+ 12 k (cid:18) ∂∂ǫ P − ∂∂ǫ C (cid:19) ǫ C + 12 k ∂ ∂ǫ V ǫ V − k ∂ ∂ǫ V ( ǫ C + ǫ V ) . (40) Analysis of Ω / terms As for previous models, these terms give the macroscopic equations. Equating both terms ofthis order on the right and left hand sides of Eq. (37) and using Eq. (38), one obtains: d [ S ] /dt = k in − k [ E ][ S ] + k [ C ] , (41) d [ C ] /dt = k [ E ][ S ] − ( k + k )[ C ] , (42) d [ P ] /dt = k [ C ] , (43) d [ V ] /dt = k [ E ] − k [ V ] . (44)In the steady-state we have the Michaelis-Menten (MM) equation: d [ P ] dt = v max [ S ] K M (1 + β ) + [ S ] , (45)where β = [ I ] /K i and K i = k /k is the dissociation constant of the inhibitor. Analysis of Ω and Ω − / 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 h ǫ S,C,V i → V , does however substantially increase thealgebraic complexity in the equations of motion for the second moments computed using terms25p to order Ω . In particular the matrix A is now a 6 × × ∂ t h ǫ S ih ǫ C ih ǫ V ih ǫ S ǫ V ih ǫ C ǫ V ih ǫ S ǫ C i = A · h ǫ S ih ǫ C ih ǫ V ih ǫ S ǫ V ih ǫ C ǫ V ih ǫ S ǫ C i + B, (46)where, A = − k [ E ] 0 0 2 k [ S ] 0 2( k + k [ S ])0 − k ( K M + [ S ]) 0 0 − k [ S ] 2 k [ E ]0 0 − k ′ − k − k − k ′ − ( k + k ′ ) − k − k − k [ S ] k [ E ] − k ( K M + [ S ]) − k ′ − k − k , (47)and B = k in M + k [ C ] + k [ S ][ E ] k ([ S ][ E ] + K M [ C ]) k [ E ] + k [ V ]00 ( k in M + k [ C ]) . (48)In the above equations we have defined k ′ = k + k . Note also that the system of equations hasbeen simplified through the application of a few row operations.Now to next order, i.e. Ω − / , the first moments of the concentrations of molecules of type C and V are governed by the equation of motions: d h ǫ C i /dt = − k ([ S ] + K M ) h ǫ C i − k [ S ] h ǫ V i + k [ E ] h ǫ S i − k Ω − / ( h ǫ S ǫ C i + h ǫ S ǫ V i ) (49) d h ǫ V i /dt = − k h ǫ V i − k ( h ǫ C i + h ǫ V i ) (50)As in previous models, since the production of product P from complex occurs through a decayprocess, it follows that at steady-state, h ǫ C i = 0 which also implies h ǫ V i = 0 from Eq. (50).Hence it follows from Eq. (49) that h ǫ S i = [ h ǫ S ǫ C i + h ǫ S ǫ V i ] / Ω / [ E ]. The two cross correlatorscan be estimated to order Ω by solving Eqs. (46)-(48). The non-zero value of h ǫ S i implies arenormalization of the substrate concentration inside the compartment and hence to a new rate26quation replacing the MM equation. This is obtained exactly in the same manner as previouslyshown for Model I. The mesoscopic rate equation is found to be given by Eq. (30) together with: f ( α ) = (1 + β ) K M [ E T ] P i =0 c i (1 − α ) i P i =0 d i (1 − α ) i , (51)where the numerator coefficients are given by: c = + k ( β + 1) K M k [ E T ] , (52) c = + K M ( β + 1) [( β + 1)[ E T ] k − (3 β + 2)[ E T ] k K M k + k βv max K M ] , (53) c = − K M ( β + 1)[2(2 β + β + 1)[ E T ] k − ((3 β + 4 β + 1)[ E T ] k K M + (54) − β ( β + 1) v max − k [ E T ] ) k + β (1 + 2 β ) k v max K M − ( β + 1)[ E T ] k v max ] ,c = + [(1 + 3 β + β + 3 β )[ E T ] K M k − ( β ( β + 1) [ E T ] k K M + (55)( − β ( β + 1) v max − β )[ E T ] k ) K M + β ( β + 1)[ E T ] v max ) k + β (1 + β ) k v max K M − (2 + 3 β + 2 β )[ E T ] k v max K M ] ,c = − [( − ( β + 1)[ E T ] k K M + β ( β + 1)[ E T ] v max ) k + (56) − [ E T ]( β + β + 1) k v max K M ] , and the denominator coefficients by: d = + K M k k (1 + β ) , (57) d = + K M k ( β + 1) [ β ( k − k K M ) + k ] , (58) d = + K M k ( β + 1) [ k [ E T ]( β + 2) + v max ] , (59) d = + ( β + 1)[ k β [ E T ] − k β k K M [ E T ] + 2 k β [ E T ] (60) − k k βK M [ E T ] − k βv max K M + k [ E T ]] , (61) d = + [ E T ] k [ k [ E T ] + k β [ E T ] + v max ] . (62)Note that P i =0 c i = 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 resultswhich are described in detail in the Results section. All simulations were carried out using27illespie’s exact stochastic simulation algorithm, conveniently implemented in the standardsimulation 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) runthe simulation and measure the ensemble-averaged substrate concentration, h n s / Ω i = [ S ∗ ] atlong times; (iii) compute α M using the MM equation, α M = [ S ∗ ] / ([ S ∗ ] + K M ); (iv) compute theabsolute percentage error R p = 100 | (1 − α M /α ) | . The solid curves in Figure 2 were obtained bynumerically solving the cubic polynomial in α given by Eqs. (7) and (8) in the Results sectionfor given values of α M and then using the above expression for R p . Figure 3 is generated in thesame manner as Figure 2, except that: in step (i) we fix M and pick a value for α between 0and 1. Since k in = M k in , the required simulation parameter is k in = αv max /M ; step (iv) is notneeded. The solid curves were obtained by numerically solving the cubic polynomial in α givenby Eqs. (7) and (9) in the Results section for given values of [ S ∗ ]. The y-axis for this figureis v/v max = α M for the MM equation and v/v max = α for the stochastic model. Figure 4 isobtained by numerically solving the quintic polynomial in α given by Eqs. (7) and (12) in theResults section together with the coefficients given by Eqs. (52)-(62) in the present section; theinhibitor concentration, [ I ], is varied while the substrate concentration, [ S ∗ ], is kept fixed. Thesubstrate concentration is chosen so that at [ I ] = 0, v/v max = 0 .
909 in all cases. Note that formodels 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. Acknowledgments
It is a pleasure to thank Arthur Straube and Philipp Thomas for interesting discussions. Theauthor gratefully acknowledges support from SULSA (Scottish Universities Life Sciences Al-liance).
References [1] Luby-Phelps K:
Cytoarchitecture and Physical properties of Cytoplasm: Volume,Viscosity, Diffusion, Intracellular Surface Area.
Int Rev Cytology :189–221282] Alberts B et al: Molecular Biology of the Cell. Garland Publishing; 1994[3] Minton AP:
How can biochemical reactions within cells differ from those in testtubes?
J Cell Sci :2863-2869[4] Trepat X et al:
Universal physical responses to stretch in the living cell.
Nature :592-596[5] Schnell S, Turner TE:
Reaction kinetics in intracellular environments with macro-molecular crowding: simulations and rate laws.
Prog Biophys Mol Biol :235–260[6] Medalia O et al: Macromolecular Architecture in Eukaryotic Cells Visualized byCryoelectron tomography.
Science :1209–1213[7] Luby-Phelps K, Castle PE, Taylor DL, Lanni F:
Hindered diffusion of inert tracerparticles in the cytoplam of mouse 3T3 cells.
Proc Natl Acad Sci USA :4910–4913[8] Pick H et al: Investigating Cellular Signaling Reactions in Single Attoliter vesi-cles.
J Am Chem Soc :2908–2912[9] Grima R, Schnell S:
Modelling reaction kinetics inside cells.
Essays in Biochemistry :41–56[10] Gillespie DT: Exact stochastic simulation of coupled chemical reactions , J PhysChem :2340–2361[11] Ben-Avraham D, Havlin S: Diffusion and Reactions in Fractals and Disordered Systems .Cambridge University Press; 2000[12] van Kampen NG:
Stochastic processes in physics and chemistry . Amsterdam: Elsevier; 2007[13] Bartholomay AF:
A Stochastic Approach to Statistical Kinetics with Applicationto Enzyme Kinetics , Biochemistry :223–2302914] Bartholomay AF: Enzymatic Reaction-Rate Theory - A Stochastic Approach , An-nals of the New York Academy of Sciences :897[15] Jachimowski CJ, McQuarrie DA, Russell ME: A Stochastic Approach to Enzyme-Substrate Reactions , Biochemistry :1732–1736[16] Grima R: Multiscale modeling of biological pattern formation
Curr Top Dev Biol :435–460[17] Berg JM, Tymoczko JL, Stryer L Biochemistry th edition. New York: Freeman; 2002[18] Popov S, Poo MM: Diffusional transport of macromolecules in developing nerveprocesses.
J Neurosci :77–85[19] Arrio-Dupont M, Cribier S, Foucault G, Devaux PF, d’Albis A: Diffusion of fluorescentlylabelled macromolecules in cultured muscle cells.
Biophys J :2327–2332[20] Ainger K et al: Transport and localization of exogenous myelin basic proteinmRNA microinjected into Oligodendrocytes.
J Cell Biol :431-441[21] Bassell G, Singer RH: mRNA and cytoskeletal filaments.
Curr Opin Cell Biol :109–115[22] Fersht A Structure and mechanism in protein science . New York: Freeman; 1998[23] Baras F, Mansour MM:
Reaction-diffusion master equation: A comparison withmicroscopic simulations.
Phys Rev E :6139–6148[24] English BP, Min W, van Oijen AM, Lee KT, Luo G et al: Ever-fluctuating single enzymemolecules: Michaelis-Menten equation revisited.
Nat Chem Biol :87–94[25] Kou SC, Cherayil BJ, Min W, English BP, Xie SX: Single-molecule Michaelis-Mentenequations.
J. Phys. Chem B :19068–19081[26] Stefanini MO, McKane AJ, Newman TJ:
Single enzyme pathways and substratefluctuations.
Nonlinearity :1575–15953027] Qian H, Elson EL: Single-molecule enzymology: stochastic Michaelis-Menten ki-netics.
Biophys Chem :565–576[28] Grima R:
Noise-induced breakdown of the Michaelis-Menten equation in steady-state conditions.
Phys Rev Letts :218103[29] McKane AJ, Nagy J, Newman TJ, Stefanini M:
Amplified biochemical oscillations incellular systems
J Stat Phys :165–191[30] Ngo LG, Roussel MR:
A new class of biochemical oscillators based on competitivebinding
Eur J Biochem :182–190[31] Davis KL, Roussel MR:
Optimal observability of sustained stochastic competitiveinhibition oscillations at organellar volumes
FEBS journal :84–95[32] Ramsey S, Orrell D, Bolouri H:
Dizzy: Stochastic simulation of large-scale geneticregulatory networks
J Bioinform Comput Biol : 415–43631able 1: Maximum Percentage error in reaction velocity from prediction of the MM equation forModel I. The copy number indicates the total number of enzyme molecules per compartment.Values in bold and in square brackets are those estimated by simulation; the italic values areobtained from the derived theoretical expressions, Eqs. (7) and (8).D/nm K M = 10 µM µM µM Copy No.50 [17.00] [4.33] [5.00] [0.74] [5.33] [2.02] [2.23] [0.52] K M = 10 µM µM µM Copy No.50 [291.56] [331.66] [58.39] [6.99] [8.50] [61.66] [66.03] [24.52] [6.06] [6.91]
CS P+S SS vesicle / granule MT
E CS P+S SS
E CS P+S SS I EI (A)(B)(C)
Figure 1:
Schematic illustrating the three models considered in this article. (A) ModelI: Michaelis-Menten reaction occurring in a compartment volume of sub-micron dimensions(shown by dashed rectangle). Substrate input into compartment occurs via a Poisson processi.e. diffusion-mediated substrate transport. (B) Model II: As for Model I but now substrateis input into compartment in groups or bursts of M molecules at a time i.e. vesicle-mediatedsubstrate transport along microtubules (MT). (C) Model III: Michaelis-Menten reaction withcompetitive inhibitor ( I ) occurring in a small subcellular compartment. Substrate transport asin previous two models. 33 p α Μ (A) α Μ (B) R p Figure 2:
Deviations from the predictions of the MM equation for diffusion-mediatedsubstrate transport. (Model I) Plot of the Percentage Error in reaction velocity, R p =100 | − α M /α | , versus the normalized reaction velocity of the MM equation, α M for 10 enzymes(green) and 100 enzymes (red) with K M = 10 µM in compartments with diameter 100nm (A)and 50nm (B) . The solid lines show the theoretical predictions, as encapsulated by Eqs. (7)and (8); the data points are obtained by stochastic simulation (see Methods for details).34 (A)(B)(C) [S*] / µ M v v max Figure 3:
Deviations from the predictions of the MM equation for vesicle-mediatedsubstrate transport. (Model II) Testing the validity of the MM relationship at small substrateconcentrations for the case in which substrate input into compartments occurs in bursts. Thedata is for 10 enzymes with K M = 100 µM in compartments of diameter (A) 200nm (circles),(B) 100nm (diamonds) and (C) 50nm (crosses); substrate is input M = 50 molecules at a time.The deterministic prediction for all three cases is the same MM equation shown by the greencurve. In contrast, the stochastic models, [Eqs. (7) and(9)], predict different rate equations foreach case (red solid lines). Data points are obtained by stochastic simulation (see Methods fordetails). Note that v/v max = α M and α for solid green and red lines respectively.35 v max [I]/[E ] T v v max [I]/[E ] T Figure 4:
Effects of intrinsic noise on the inhibition of enzyme activity in smallcompartments. (Model III) Plots of normalized enzyme activity versus normalized inhibitorconcentration (measured in units of the total enzyme concentration [ E T ]) for 10 enzymes with K M = 100 µM in compartments of diameter 50nm and 100nm (inset). The colors correspondto: (red) MM equation; (green) stochastic model, M = 20; (blue) stochastic model, M = 50.The latter two curves are those predicted by theory [Eqs. (7) and(12)]. Parameters same asmentioned in caption of Table 3 (except for [I], which is a variable in the present case). Substrateconcentrations chosen so that at [ I ] = 0, v/v max = 0 .
909 in all cases. Black dashed lines contrastthe inhibitor concentration required to decrease enzyme activity from 0.909 to 0.1 as predictedby the MM equation and the stochastic models. Note that v/v max = α M and α for solid redand blue/green lines respectively. 36able 3: Maximum Percentage error in reaction velocity from prediction of the MM equation forModel III. The total number of enzyme molecules per compartment is ten in all cases. Valuesin bold and in square brackets are those estimated by simulation; the italic values are obtainedfrom the theoretical expressions, Eqs. (7) and (12). The parameters are: k = 10 s − M − , k = k = 1000 s − , k = 10 s − M − , and [ I ] = 10[ E T ].D/nm K M = 10 µM µM µM M (burst size)50 [76.8] [76.5] [26.4] [26.1] [169.4] [345.5] [75.2] [31.5] [11.5]3.6