In the past decade, advances in molecular biology such as the development of non-invasive single molecule imaging techniques have given us a window into the intricate biochemical activities that occur inside cells. In this chapter we review four distinct theoretical and simulation frameworks: (i) non-spatial and deterministic, (ii) spatial and deterministic, (iii) non-spatial and stochastic and (iv) spatial and stochastic. Each framework can be suited to modelling and interpreting intracellular reaction kinetics. By estimating the fundamental length scales, one can roughly determine which models are best suited for the particular reaction pathway under study. We discuss differences in prediction between the four modelling methodologies. In particular we show that taking into account noise and space does not simply add quantitative predictive accuracy but may also lead to qualitatively different physiological predictions, unaccounted for by classical deterministic models.
Over the past decade there has been an explosion of interest in methods for studying intracellular reactions, in particular those focusing on the determination of kinetic parameters . Typically data obtained from experimental assays must be input into a mathematical model before they can be meaningfully interpreted. The bulk of these models are based on the macroscopic and deterministic models of classical physical chemistry . A prominent example of this analysis is the characterization of enzymes using the Michaelis–Menten equation .
The use of any type of modelling methodology implicitly assumes that a number of physical, chemical or biochemical constraints are satisfied. Using a model outside of its range of application naturally implies that it is not possible to guarantee the accuracy of the model predictions. In this chapter we discuss the four main modelling strategies in the context of intracellular kinetics. In particular we show how basic experimentally obtainable parameters can be used to roughly determine which model is appropriate for studying the reaction pathway under investigation. First, we discuss the homogeneity and continuum assumptions inherent in the use of classical biochemical models. We show how the validity of these assumptions can be checked by determining the fundamental length scales which characterize a given system. Secondly, we briefly survey literature on various intracellular biochemical processes. We show that, in a significant number of cases, it is not possible to uphold both the continuum and homogeneity assumptions at the same time. These cases are amenable to study via stochastic and discrete models, which are the third topic of discussion. Fourthly, we discuss the effects of noise (stochastic compared with deterministic models) and space (homogenous compared with heterogeneous models) on physiological predictions. We conclude the chapter with a brief discussion and summary.
Classical models of chemical reaction kinetics
All reaction pathways can be decomposed into a series of elementary reaction steps, each of which is either unimolecular or bimolecular. By simulating each key step of a reaction, one can thus construct a model of the original pathway, independent of its complexity . For these reasons we limit our discussion to modelling only these two types of elementary reactions. Consider a chemical reaction occurring in a small cubic container with sides of length L. If the mixture is homogeneous throughout the container at all times, a condition which can be obtained by constant stirring, the state of the mixture at any point in time is fully described by the total number of particles of each molecular species in the container. Let nI be the total number of particles of species I in the container; then the concentration of the species, denoted by [I], is equal to . For a reversible bimolecular reaction, , the forward rate of reaction (i.e. the number of molecular events leading to a successful reaction per unit volume per unit time) is k1[A][B], whereas the backward rate of reaction is k−1[C]. The RE (rate equation) for species C is then given by the familiar classical equation (eqn 1): A crucial assumption implicit in REs is that the local concentration at each point in space inside the container must equal the global concentration at all times. In other words, all concentration gradients must be zero and, as a result, a non-spatial model is enough to describe the kinetics. This is called the homogeneity assumption. In industrial applications this is easily satisfied by constant (convective) stirring. Inside cells, homogeneity is only obtained if the effective diffusion coefficients are sufficiently large. The latter cannot be assumed a priori. The exact mathematical criterion for diffusion to suffice to maintain homogeneity inside our container is Dι≫L2, where D is the molecular diffusion coefficient and ι is the average lifetime of a reactant molecule . The distance is a measure of the distance travelled by a molecule during its lifetime. It is sometimes referred to as the Kuramoto length [5,6]. If lk≫L, then any fluctuation in concentration at a particular point inside the container quickly diffuses to all parts of the container and homogeneity is maintained at all times. On the other hand, if lk≪L, the fluctuation remains localized in the region where it occurred and distant parts of the container will have different concentrations and, as a result, homogeneity cannot be assumed.
Estimation of the Kuramoto length depends on a determination of the molecular lifetime. We illustrate how this calculation proceeds in the context of the simple reversible bimolecular reaction introduced above. Suppose that the mixture of reactant molecules is initially perfectly homogeneous throughout the container and that the reaction is at equilibrium such that the concentrations of species A, B and C are equal to [Aeq], [Beq] and [Ceq] respectively. Let δ1 present the fluctuation in the concentration of species I. The average lifetime of a molecule of type C can then be obtained by substituting [A]=[Aeq]+δA,[B]=[Beq]+δB and [C]=[Ceq]+δC in the rate equation (1) and solving this expression for δC by ignoring terms involving δC2 (i.e. assuming fluctuations are very small) and applying the equilibrium constraint δC=−δA=−δB. This gives the solution δC∝exp(−t/τ) where ι is the time for fluctuations to decay given by τ = (k1([Aeq]+[Beq])+k−1)−1. This timescale is a measure of the lifetime of a molecule of species C and can, therefore, be used to estimate the Kuramoto length of species C (another example can be found in ). Note that since ι is inversely proportional to the reaction rate, the Kuramoto length scale represents the relationship between diffusion and reaction rates.
In general, in order for homogeneity to be assumed, the Kuramoto length of each molecular species in the pathway must satisfy the restriction that it be much larger than the container size. If this is not the case, then one necessarily needs a mathematical description of the time evolution of the local concentrations, in other words, a spatial model in which diffusion is explicitly described. This can be accomplished by dividing the volume of the container into many small cubic elements of linear size δL, such that in each element homogeneity prevails, i.e. δL≪lk. The rate of change in the concentration of species C in each element is due to (i) reaction between molecules of type A and B inside the element and (ii) diffusion of the same molecules of type C from neighbouring elements. Then it follows that the dynamics in each element is described by an equation of the same form as eqn (1) with an additional diffusion term (eqn 2): This is called an RDE (reaction-diffusion equation). While REs and RDEs differ by whether the homogeneity assumption holds over the whole or parts of the volume in which the reaction occurs, they do share a common implicit assumption: both are deterministic approaches. The concentrations in these equations are mean concentrations and thus they can only be valid when the stochastic fluctuations about the mean are very small. The fluctuations in a collection of n particles are of order n−1/2 . Thus for a given concentration, the larger the volume of the cubic element in which homogeneity prevails, the larger is n and consequently the smaller are the size of the fluctuations compared with the mean concentration. This necessarily implies the conditions: li≪ δL and li ≪ L where li is the intermolecular distance. Under these conditions the discreteness of molecules is non-apparent and thus this is also called the continuum assumption.
Note that in our discussion above, we have assumed that molecules could reach any part of the container via diffusion. In this case it is clear that the size of the well-mixed region of space can only depend on the magnitude of the diffusion coefficients and the molecular lifetime (i.e. the Kuramoto length scale). However, if physical structures exist inside the container which constrain molecular movement to some specific spatial regions characterized by a length scale ls, then concentration homogeneity can only exist on a spatial scale δL if δL≪min(ls, lk).
In summary, the classical approach is only valid if the reaction volume can be subdivided into spatial elements (of linear size δL) which are small enough to ensure well-mixed conditions inside each element [the homogeneity assumption: δL≪min(ls, lk)] but large enough to ensure that molecular discreteness can be ignored (the continuum assumption: li≪δL). These are the two main assumptions implicit in the use of classical reaction kinetic equations. Diffusion-limited coalescence is a well-studied case in which the above compromise is not possible (see  for a detailed discussion). We now turn to a discussion of whether both classical assumptions hold when modelling intracellular reaction kinetics.
Can we understand intracellular reactions using classical models?
Although there is limited information on the concentration and number of molecules inside cells, diverse experimental evidence (see below) suggests that many intracellular biochemical reactions involve nanomolar (nM) or smaller concentrations of reactant molecules [8–10]. A simple calculation shows that in absolute number terms, a 1 nM concentration approx. corresponds to 2500 molecules inside an average mammalian tissue culture cell (approximated as a sphere with a diameter equal to 20 μm), and to just one molecule inside an Escherichia coli bacterial cell (approximated as a cylinder with 1 μm diameter and 2 μm in length). Absolute numbers of regulatory molecules have been reported to be approx. 1000 per cell for two endogenous yeast proteins , at most a few thousand per cell for proteins involved in the mammalian circadian clock , and to range from tens to hundreds per cell [13–15] for molecules in some bacteria. In contrast, only a few hundred are required to induce chemotaxis in amoeboid cells , to extend the growth cone in chick dorsal root ganglion cells  or to induce an intracellular calcium response in HeLa cells .
A 1 nM concentration corresponds to an average intermolecular distance, li, approximately equal to 1 μm whereas 100 nM corresponds to li approximately equal to a quarter of a μm. Note that the size of the intermolecular length scale is comparable with the typical size of many intracellular compartments inside cells (e.g. the size of a lysosome is 0.2–0.5 μm, the size of a mitochondrion is 0.5–1 μm and the size of a nucleus is 3–10 μm ).
It is important to remember that not all metabolites inside cells have concentrations in the nanomolar range. A number of glycolytic intermediates in the bacterium Lactococcus lactis are reported to be in the tens of millimolar range . Such a concentration corresponds to an average li approximately equal to 10 nm, a length scale much smaller than that of intracellular compartments.
As discussed in the previous section, the homogeneity assumption holds for length scales much smaller than both the typical size of an intracellular compartment (the latter length scale is equal to ls) and than the molecular Kuramoto lengths. Since intermolecular distances are comparable with the size of intracellular compartments, the number of reactant molecules in spatial regions of this size will typically be so small that their discrete nature cannot be neglected and the continuum approximation breaks down.
It thus appears inevitable that a reasonable model of intracellular chemical kinetics will need to relax the requirement that the continuum approximation be valid for any length scale inside the cell along with the constraint that the homogeneity approximation hold for fine enough length scales. Such models are the subject of the next section.
Discrete and stochastic models of chemical reaction kinetics
The previous section demonstrated that nanomolar concentrations of metabolites, typical in many cells, are consistent with small number of molecules in regions of the size of organelles. This implies that the continuous concentration of classical models of chemical reaction kinetics needs to be replaced with the absolute and discrete (integer) number of molecules: stochasticity cannot be neglected.
There are two types of discrete and stochastic formulations of chemical kinetics, paralleling exactly the classical RE and RDE approach discussed above. We will first discuss the case where reactants are homogeneous and then the case where they are heterogeneous in the container. In this case, the term container refers to the cell. The small elements of space into which it is subdivided are the intracellular compartments or, perhaps, even smaller spatial regions.
Stochastic kinetics in homogeneous conditions
Let us imagine a non-reversible bimolecular chemical reaction occurring in a container of size L between two different chemical species, . We assume that reactions occur in well-mixed conditions, i.e. L≪lk, but do not require that L is much larger than the typical intermolecular distance. We are, therefore, dealing with discrete (integer) numbers of molecules instead of continuous, real-valued concentrations.
Since we are assuming well-mixed conditions, the positions and velocities of individual molecules can be ignored and thus the state of the system at any time is completely specified by the number of molecules of both species (the number of molecules of type C can be ignored since they do not contribute to the reaction kinetics). Let P(nA, nB, t) be the probability that in the container there are nA molecules of species A and nB molecules of species B at time t. The state vector is (nA, nB). Applying the laws of probability yields the following equation describing the time evolution of the system (eqn 3): where T(nA′,nB′|nA,nB)Δt gives the probability of a single bimolecular reaction event occurring in a time interval Δt inside the container and causing the state vector to change from (nA, nB) to (nA′,nB′). The first term represents the state change (nA+1,nB+1)→(nA,nB) whereas the second term represents the state change (nA,nB)→(nA−1,nB−1). Eqn (3) is called the Master Equation or the CME (chemical master equation) for the bimolecular reaction. The function T is commonly referred to as the propensity function. For the reaction under consideration, T has to be proportional to both k (the familiar macroscopic rate constant) and to the number of A− B pairs in the reaction volume but inversely proportional to the container volume (since a larger volume implies rarer successful encounters). We therefore have eqn (4): This equation, together with eqn (3), completes the discrete and stochastic description of a bimolecular reaction. A similar procedure can be used to derive CMEs for any reaction process. A detailed discussion of the general derivation of propensity functions from microscopic physics can be found in .
CMEs can only be exactly solved in a few simple cases and, as a result, computer simulation is frequently essential. The most common simulation method is the SSA (stochastic simulation algorithm) of Gillespie [24,25]. A discussion of this approach and its variants in the context of reactions inside cells can be found in a review by Turner et al. . There are several readily available stochastic simulation software packages for the simulation of CMEs via the Gillespie algorithm, e.g. Dizzy , BioNetS  and Dynetica .
Many of the current models assume well-mixed conditions throughout the whole volume of space in which the intracellular processes of interest occur. These models are built by decomposing the relevant reaction pathway into a series of unimolecular and bimolecular elementary reaction steps and then writing master equations for each step. Such a procedure has now been applied to study a significant number of intracellular processes, including genetic regulatory networks , enzymatic substrate cycles (a recurring control motif in biological molecular networks) , single enzyme kinetics [32,33] and the circadian clocks of Drosophila, Neurospora  and of mammals .
We finish this section by noting that up until now we have assumed that the collision frequency is proportional to the product of the number of particle species; this is evident in the forms of eqn (1) and eqn (4). If the Maxwell velocity distribution is maintained [5,25] then the above assumption is correct. When the dominant transport process is anomalous diffusion, the molecular velocities do not obey the Maxwell distribution; this is the case when there exists significant molecular crowding effects in the region of the cell where the reactions occur (three examples of this are protein motion inside cell membranes , in the nucleoplasm  and in the cytoplasm of HeLa cells and E. coli [37,38]). The description of reaction kinetics in the presence of crowding requires more sophisticated approaches than the CME.
Stochastic kinetics in heterogeneous conditions
We now consider the case where the molecular Kuramoto lengths are much smaller than the container, thus leading to significant spatial heterogeneity in the concentrations. To derive the stochastic equations governing such conditions, we adopt the same conceptual approach used previously to derive the deterministic reaction-diffusion equations with one important distinction: the continuum approximation is not made. We divide the reaction volume into m small cubes with side length δL, such that the reaction is well-mixed in each cube, i.e. δL≪lk. The reaction inside each cube can be treated in the same way as in the previous section. We must now, however, consider diffusion between neighbouring cubes (as was also the case for RDEs). Diffusion of species A between two adjacent cubes, i and j, can be represented as where kD = D/zδL2 is the rate at which molecules diffuse between the two cubes (z is the number of nearest adjacent cubes). Hence we see that diffusion is mathematically equivalent to a set of unimolecular processes. We can therefore describe both reaction and diffusion throughout space by writing the CME for the set of reaction processes (eqn 5): This is frequently termed a multivariate master equation  or an RDME (reaction-diffusion master equation). A thorough analytical treatment of these equations can be found in the recent paper by Isaacson and Peskin .
Efficient simulations of RDMEs in biological contexts can be achieved by extensions and variants of the Gillespie algorithm, e.g. SmartCell  and MesoRD . A crucial issue to be resolved before any spatially extended system can be simulated using RDMEs is the determination of the size δL of the small cubic elements in which homogeneity is assumed to prevail. The most straightforward way to ascertain the appropriate cube size to use for a given simulation requires the comparison of the frequency of molecular transfer (due to diffusion) between adjacent cubes, fd, and the frequency of chemical reactions occurring inside the cubes, fc. When fd≫fc, homogeneity inside each cube is guaranteed; otherwise, a smaller cube size must be selected .
It is important to note that, in addition to the use of RDMEs (and their simulation via variations of the Gillespie algorithm), there are other mathematical and simulation frameworks available for modelling stochastic kinetics in heterogeneous conditions. Field-theoretic renormalization group methods (see  for a review) provide an alternative analytical approach whereas software based on Brownian dynamics (e.g. Smoldyn , Virtual Cell  and M-cell ) provides spatial stochastic simulation with single molecule detail (the latter feature is absent in software based on the Gillespie algorithm).
There have been only a few instances to date where spatial stochastic models have been used to model realistic intracellular kinetics; these include the bacterial chemotaxis signalling pathway , Min-protein oscillations in E. coli (see MesoRD website at http://mesord.sourceforge.net/), neurotransmission in the chick ciliary ganglion  and calcium waves in differentiated neuroblastoma cells .
In the above treatment we considered kD to be a constant throughout space; by allowing this parameter to be a function of space one obtains a stochastic and discrete model of intracellular reactions and the complex intracellular environment in which they occur. For example, a permeable membrane could be modelled by requiring the value of kD for each cube discretizing the membrane to be different from the value of kD for each cube discretizing the surrounding cytosol. In this manner any number of complex cytoplasmic geometries may be easily modelled. We note that intracellular processes in which anomalous diffusion is the dominant transport process can be modelled with RDMEs by introducing a spatial distribution of impenetrable obstacles (crowding agents) in the reaction volume (for a discussion of the relationship between anomalous diffusion and obstacle concentration see Saxton ).
The model-dependence of physiological predictions
As remarked in the previous section, the majority of current stochastic models of intracellular kinetics assume well-mixing in the intracellular region where the reactions occur. It is important to remember that these models are simply an approximation of spatially extended stochastic models and, as such, their predictions can be expected to be physically realistic only when the Kuramoto length of each interacting molecular species is significantly larger than the cell diameter. If the rate coefficients and the approximate number of molecules are known then the above condition can be checked by explicit calculation using the formulae discussed in the beginning of this chapter. Otherwise, one may experimentally tag key molecular species with fluorescent protein and analyse the spatial distribution over a period of time.
In general, if the reaction scheme under consideration involves one or more diffusion-limited reactions, i.e. one in which the probability of reaction upon encounter is close to one (e.g. models of the kinetics of super-efficient enzymes), then the predictions of homogeneous models may be incorrect. In these cases reactant molecules move a small distance before reacting, implying that concentrations in different regions of space will evolve with time in a different manner. For example, homogeneous stochastic simulations of the diffusion-limited bimolecular reaction with equal initial reactant concentrations via the Gillespie algorithm and analysis via the Van-Kampen expansion  show that the mean reactant concentration follows the classical decay law [A]=[B]∝1/t (in all spatial dimensions) whereas spatially extended simulations  show that the correct long-time decay law is [A] = [B]∝t−d/4, where d is the dimension of space. Note that the reaction proceeds more slowly than predicted by the corresponding non-spatial models; this is because, during the course of the reaction, alternating domains of A-only and B-only particles are formed and the reaction can thus proceed only at the domain boundaries rather than throughout the whole volume of space, as assumed by non-spatial models (Figure 1). A thorough discussion of this phenomenon can be found in .
Note that discrepancies between the concentration decay laws predicted by non-spatial and spatial models increase with decreasing dimension. This has profound implications for the modelling of intracellular biochemical networks in which some of the reactions occur via an effective lower-dimensional motion. The use of combinations of three-dimensional (in the cytosol) and one- or two-dimensional diffusion (e.g. gliding along a linear polymer or a membrane surface) is a mechanism first proposed by Adam and Delbruck . This facilitates target location  and hence it is probably quite common inside cells, e.g. the one-dimensional gliding of the tetrameric E. coli lac repressor and RNA polymerase along the DNA chain to locate specific binding sites [56,57].
Next we briefly discuss the predictions of stochastic and discrete models with those of the corresponding deterministic and continuous models for a given biological or biochemical system. Although it is now recognized that stochastic effects play important roles in cell biology, “they are often thought of as providing only moderate refinements to the behaviours otherwise predicted by the classical deterministic system description” . Indeed there are many cases in which a deterministic description suffices to capture the approximate, qualitative physiological behaviour exhibited by a real biological system. Some examples include models of the budding yeast cell cycle , of the glycolysis pathway in yeast  and L. lactis  and of feedback regulation in the lactose operon . However, there are also many instances in which it has been shown that stochastic effects lead to physiological predictions which cannot be reproduced qualitatively or quantitatively by the corresponding deterministic models. Prominent examples of such cases occur in models of biochemical oscillations (e.g. circadian rhythms, calcium oscillations and oscillations of ADP and ATP concentrations in glycolysis) where it has been shown that noise generally expands the region of parameter space in which cycles are predicted to occur by deterministic models (Figure 2) [12,31,62,63]. In some cases noise may even induce oscillations in systems in which the corresponding deterministic model predicts that no oscillations are possible over the entire parameter space. For example, noise may induce cycles in reaction schemes involving just two chemical agents whereas deterministic models predict that at least three chemical agents are necessary . In general, it has been found that fluctuations may lead to steady-states in chemical reaction-diffusion systems which are not predicted by RDEs .
In this chapter we have discussed the four main methodologies for modelling reaction kinetics in cells (Figure 3). The determination of which model is most relevant for the analysis of a particular intracellular reaction pathway is closely tied to the physical nature of the intracellular region in which it occurs and also to the typical chemical concentrations involved in the reaction.
In particular one needs to estimate the Kuramoto length scale, the size of the intracellular region in which the reaction occurs and the average intermolecular distances. These can be estimated from knowledge of the diffusion coefficients, the reaction constants and molecular concentrations. For example, consider a reaction confined to an intracellular compartment of μm size and in which the reacting species have a hydrodynamic radius equal to 4 nm (D~0.05 nm2/ns), a bimolecular rate constant equal to k = 0.5 ×106 M−1·s−1 and a typical concentration of 1 μM. This implies lk~10 μm and li~0.1 μm; since the compartment size is much smaller than lk and significantly larger than li then an RE approach may suffice. If crowding caused the molecular diffusion coefficient to be reduced by five times, the typical concentration was 100 nM and the rate constant was equal to k = 5 ×106 M−1·s−1 then lk~4.5 μm and li~0.3 μm. If this reaction forms part of a pathway which involves both the nucleus (3–10 μm) and mitochondrion (0.5–1 μm) of a mammalian cell then both homogeneity and continuum assumptions do not hold and an RDME approach will be needed.
Roughly speaking, reactions involving concentrations in the millimolar range will probably satisfy the continuum assumption whereas those involving smaller concentrations, such as nanomolar ones, will not. Reactions occurring exclusively in a small intracellular compartment will favour the homogeneity assumption; those involving multiple transport between several compartments and the cytosol or in regions where molecular crowding is significant will not satisfy this assumption.
As mentioned in the Introduction, choosing the right model for the application at hand ensures the appropriate interpretation of experimental data. An example of this is the analysis of data obtained from SMI (single-molecule imaging) techniques [64–68]. The large fluctuations inherent in this data would necessarily have to be ignored if analysis proceeds via a deterministic framework. However, analysis via the stochastic framework of master equations has shown that such fluctuations may contain important information about molecular behaviour; for example, it was recently shown that fluctuations in the rate constants of a single enzyme yield information regarding its conformational dynamics [33,69].
We conclude by noting that in general it is not possible to a priori rule out the effects of noise or space on a given intracellular reaction scheme of interest. It is important to be aware of the fact that non-spatial deterministic models (based on classical REs) constitute only a first approach to modelling intracellular dynamics; their physiological predictions should be carefully checked against those of spatial stochastic models which correspond more closely to reality.
• Intracellular reaction kinetics can be modelled using four main different approaches (RE, RDE, CME, RDME); the right model choice for a particular reaction pathway under study depends on the concentrations and Kuramoto lengths of all molecular species, the size of intracellular space in which the reaction occurs and the extent of macromolecular crowding in this region.
• In general, stochastic effects, typical of reactions involving nanomolar concentrations, do not simply add moderate quantitative refinements to the behaviours otherwise predicted by the classical deterministic models (RE, RDE). Noise may induce novel and qualitatively different physiological behaviour than that predicted by deterministic models.
We are very grateful to Michelle L. Wynn for her critical comments during the revision of this chapter. This work is based upon work supported by the NIH/NIGMS (National Institutes of Health/National Institute of General Medical Sciences) under Grant No. R01 GM076692.
- © The Authors Journal compilation © 2008 Biochemical Society