This Article Abstract Full Text (PDF) Supplemental Research Data Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Elf, J. Articles by Ehrenberg, M. Search for Related Content PubMed PubMed Citation Articles by Elf, J. Articles by Ehrenberg, M. Social Bookmarking                   What's this?

Methods

Fast Evaluation of Fluctuations in Biochemical Networks With the Linear Noise Approximation

Department of Cell & Molecular Biology, Uppsala University, BMC, 751 24 Uppsala, Sweden

	ABSTRACT

Top ABSTRACT RESULTS DISCUSSION METHODS REFERENCES

 
Biochemical networks in single cells can display large fluctuations in molecule numbers, making mesoscopic approaches necessary for correct system descriptions. We present a general method that allows rapid characterization of the stochastic properties of intracellular networks. The starting point is a macroscopic description that identifies the system's elementary reactions in terms of rate laws and stoichiometries. From this formulation follows directly the stationary solution of the linear noise approximation (LNA) of the Master equation for all the components in the network. The method complements bifurcation studies of the system's parameter dependence by providing estimates of sizes, correlations, and time scales of stochastic fluctuations. We describe how the LNA can give precise system descriptions also near macroscopic instabilities by suitable variable changes and elimination of fast variables.

Chemical reactions in the test tube are normally described macroscopically by (1) applying the law of mass action (Cornish-Bowden 1995; Fersht 1998) and (2) neglecting random deviations from the average number of molecules during the time evolution of the concentrations involved in the studied reactions. Such macroscopic descriptions have in the past also been used to describe intracellular chemical reactions with some notable exceptions in which stochastic reasoning has been applied (e.g., Berg 1978). More recently, however, the urgent need to take stochastic aspects of chemical reactions in single cells into account has been realized, and this insight has resulted in a burst of publications dealing experimentally as well as theoretically with the intrinsically random nature of gene expression (Becskei and Serrano 2000; Kepler and Elston 2001; Paulsson and Ehrenberg 2001; Thattai and van Oudenaarden 2001; Elowitz et al. 2002; Ozbudak et al. 2002). The focus on gene expression is explained by the fact that the copy number of individual messenger RNAs can often be very small, in particular in bacterial cells like E. coli, and this frequently leads to highly significant relative fluctuations in messenger RNA copy numbers and, sometimes, also to large fluctuations in protein levels.

To deal with the stochastic nature of these reactions, one has to abandon the macroscopic perspective of chemical reactions, where fluctuations are neglected, in favor of mesoscopic approaches (van Kampen 1997). One technique is to solve the chemical Master equation, which describes the time evolution of the probability distribution for all the different molecular components in the system (van Kampen 1997), either analytically or numerically. Alternatively, one can simulate the molecule number trajectories for such systems, using Monte Carlo methods, implemented by, for example, the Gillespie algorithm (Gillespie 1976, 1977). These well-established approaches are based on an assumption of homogeneous probability distributions in space, corresponding to much faster diffusion of all components throughout the system than the rates of all chemical reactions (van Kampen 1997). It is likely that this assumption will turn out to be untenable for many intracellular reactions, and attempts have been made to develop Monte Carlo methods to simulate space-time trajectories for stochastic reaction-diffusion-coupled systems (Stundzia and Lumsden 1996). Recently, we developed an efficient algorithm (Elf et al. 2003a) for such simulations, built on the Gillespie and the next reaction algorithms (Gibson and Bruck 2000), but in the present work we assume spatial homogeneity throughout.

Large relative fluctuations of copy numbers characterizes many mRNA levels in bacterial cells, but there are also other cases in which noise becomes highly significant, also when the average molecule copy numbers are very large. One example is systems with a continuum of macroscopic stationary states, like a limit cycle attractor (Barkai and Leibler 2000; Qian et al. 2002). Another example is systems that operate near critical points (Dogterom and Leibler 1993; Elf et al. 2003b). In all these cases, relaxation to a steady-state solution is slow or a steady state may not exist at all, meaning that very large and slowly changing fluctuations in component copy numbers can emerge.

No simple methods for the mesoscopic analysis of such systems exist today, and numerical approaches are often difficult to carry out because of the high computational demands normally associated with this type of problem. Monte Carlo methods, such as the Gillespie algorithm (Gillespie 1977), are useful only for very limited sets of parameters, because a complete simulation must be carried out for every chosen set of parameter values.

Here, we present a general method, based on the linear noise approximation (LNA; van Kampen 1997), that allows rapid characterization of the stochastic properties of intracellular networks over large parameter regions. Our analysis is useful for initial characterization of fluctuations in all kinds of intracellular reactions. We describe the method qualitatively and also give four simple, and yet realistic, examples that clarify its principles and demonstrate its applicability. The examples highlight: (1) the importance of taking stoichiometric relations into account in mesoscopic descriptions, although they enter macroscopic models in a trivial way; (2) the fundamental connection between, on one hand, the sensitivity of a system response to a change in a flow and, on the other, the relative size of random fluctuations; (3) coupled fluctuations in metabolic networks; and (4) how the LNA can be made more accurate close to near-critical points by elimination of fast variables (Gardiner 1985).

All examples are chosen so that analytical approximations of solutions of the Master equation can be obtained, and these approximations are then used to clarify the origin of fluctuations in each case and to trace the effects of different parameters on the stochastic properties of the system.

How the present theory can be applied to arbitrarily large biochemical networks is described in the Discussion.

	RESULTS

Top ABSTRACT RESULTS DISCUSSION METHODS REFERENCES

 
The Linear Noise Approximation
The chemical Master equation describes how the joint probability distribution of the copy numbers of different chemical species evolve in time in a spatially homogenous chemical system (van Kampen 1997). As shown in the Methods section, the Master equation can be written directly from the rate constants and stoichiometries of all the elementary reactions of a chemical system, but neither analytical nor numerical solutions are in general available. Fortunately, the Master equation can often be simplified in a linear noise approximation (LNA). This is obtained through van Kampen's -expansion of the Master equation (van Kampen 1997). In this formulation, the Master equation allows for analytical solutions that are locally valid close to macroscopic trajectories or stationary points of the system. These solutions often give very accurate descriptions of the size of molecule number fluctuations and how they are correlated. They can therefore be used for rapid, initial characterizations of the stochastic properties of chemical systems as guidelines for decisions about whether the fluctuations are so large that mesoscopic considerations are necessary or whether a macroscopic approach will suffice. The -expansion means that the Master equation is Taylor-expanded near macroscopic system trajectories or stationary solutions in powers of 1/, where is the system volume. When the Master equation is approximated near a macroscopically stable stationary solution, terms of first order in 1/ give the macroscopic rate equations, and terms of second order in 1/ give the Linear Noise Approximation (LNA). The rationale behind this approach is that for constant average concentrations, relative fluctuations will tend to decrease with the inverse of the square root of the volume. The LNA is therefore very accurate for systems for which the local linearization of the chemical rate laws is valid for component copy numbers that are frequently reached by fluctuations away from the stationary state. This is true when fluctuations are sufficiently small in relation to the component means. However, we have found that the LNA can give very good estimates of fluctuations in molecule numbers, also when they are larger than the corresponding means, and, furthermore, that the technique of separation of fast variables (see below) can be used to improve further the performance of the LNA.

Here, the LNA is derived for the stationary state of a general system of chemical reactions (Supplemental material available online at www.genome.org; Methods). The correlation matrix, C, containing the variances and covariances for the fluctuations of all the components in the system, is calculated from the linearized macroscopic rate equations near the stationary point as described by the Jacobian matrix, A, and the stoichiometries of the reactions and their rates as described by the diffusion matrix, D. The relation between the matrix C, containing information about the fluctuations, and the macroscopic matrices A and D is given by a Lyapunov matrix equation, equation 32.

(1)

How to calculate the A and D matrices from the chemical reactions and how to solve the equation for C is described in Methods. Here we focus on a few special cases in which the LNA is used to understand the origin of noise in biochemical processes.

The first example is a system that contains only one type of molecule, X. These molecules are synthesized and consumed by two different reactions with rate laws f₁() and f₂(), respectively. The reaction rates are given in molar concentration per second and is the molar concentration of X. The first reaction adds n molecules to X in a single event, and the second reaction removes m molecules from X in a single event:

(3)

When equation 3 is linearized around the steady state, , the equation for the deviation from steady state is

(5)

The variance _X² normalized to the average number (X) of molecules is often called the Fano-factor (Fano 1947). Open systems at thermodynamic equilibrium obey Poisson statistics, and for those the Fano-factor is always 1. For the creation-elimination scheme in equation 2, the Fano-factor is approximated by

(6)

To obtain this expression, we have also approximated X with the macroscopic steady-state value of the number of X molecules in the steady state. The numerical magnitude of the Fano-factor is primarily determined by the ratio between the rate f₁()/ of turnover of all the molecules in the X-pool and the rate of relaxation . The Fano-factor is also shaped by the creation (n) and elimination (m) stoichiometries. Below we use equation 6 to inspect the role of the stoichiometries (Example 1) and the role of the relaxation rate (Example 2) for the stochastic properties of the scheme in equation 2.

Example 1: The Mesoscopic Importance of Stoichiometries in Chemical Reactions
At the macroscopic description level, one can often compensate for a change in the stoichiometric order of a reaction by a change in a rate constant, so that the dynamics of the reaction remains unchanged. At the mesoscopic level, this cannot be done.

To illustrate the role of the stoichiometries, we compare two different reaction systems (A and B) with different reaction stoichometries:

(8)

From equation 5, the following approximations for the variances in the two cases are obtained:

(9)

The Fano-factors (equation 6) for systems A and B are (n + 2)/4 and 3/4, respectively. For system A, the Fano-factor increases linearly with the stoichiometry of the creation reaction, whereas the corresponding increase in the rate constant (but not the stoichiometry) for system B leaves the Fano-factor unchanged. This example is described in more detail and with general stoichiometries in the Supplemental material.

Example 2: Zero-Order Sensitivity and the Size of Mesoscopic Fluctuations
Goldbeter and Koshland Jr. (1981) pointed out that the extent of modification in certain enzymatically catalyzed modification-demodification reactions can be ultrasensitive to the amount of either the modifying or the demodifying enzyme. It was also shown (Berg et al. 2000) that such systems can display large fluctuations, which will tend to attenuate the macroscopically predicted ultrasensitivity.

We exemplify this by treating the case with a protein that is converted from its unmodified state X₁ to its modified state X₂ by one enzyme and back to X₁ from X₂ by another. With a constant total concentration ₀ of modified and unmodified target concentrations, the state of the system is completely described by the concentration ₁ of X₁, because the concentration ₂ of X₂ is given by the conservation relation ₀ = ₁ + ₂. The reaction scheme is

(10)

The macroscopic rates of modification and demodification are f₁(₁) and f₂(₁), respectively. The modifying enzymes obey Michaelis-Menten kinetics, and when steady-state relations for their internal states are rapidly established (see discussion about elementary complex reactions in Methods), the rate equations are given by

(11)

The same K_m value is assumed for both enzymes, and the fractions of enzyme-bound target molecules are assumed to be negligible. The macroscopic dynamics is thus given by

(12)

with the steady-state condition

(13)

The macroscopic stationary value ₁ is the LNA estimate of the average number of X₁ molecules. In Figure 1A, this estimate is compared with the correct value as calculated directly from the Master equation, for different values of V_max2. The sensitivity amplification (logarithmic gain; Savageau 1976) of the system is

(14)

The sensitivity amplification is the percentage change in the concentration of ₁ in response to the percentage change in V_max2. The system can be said to be ultrasensitive if |a_1Vmax1| > 1 (Goldbeter and Koshland Jr. 1982). If ₀ >> K_m, both enzymes are close to saturation when V_max1 = V_max2 and the reaction displays ultrasensitivity (Goldbeter and Koshland Jr. 1981), because of zero-order kinetics in both the modification and demodification reactions. Zero-order kinetics means that the flow in either direction is insensitive to the concentrations of modified and unmodified substrates, so that the "force" that restores the system to its steady state is very weak.

View larger version (18K): [in this window] [in a new window]   Figure 1 Properties of the modification-demodification system (one-component system). (A) Average number of molecules in X1 state, X1. (B) The Fano-factor (C/X1). The lines correspond to the LNA equation 15, the circles to values from the numerical solution of the corresponding Master equation. The total concentrations of X1 and X2 are 10-5 M and vmax1 = 10-6 M/sec.

(15)

The equivalence between the sensitivity in equation 14 and the Fano-factor in equation 15 is seen by implicit derivation of the steady-state equation (equation 13; Elf et al. 2003b):

(16)

The maximal value of the Fano-factor, (2K_m + ₀)/(2K_m), is obtained when V_max1 = V_max2. The LNA of the Fano-factor (equation 15) is plotted in Figure 1B for increasing V_max2 and is compared with values calculated from the exact Master equation for some parameter sets. The comparison reveals that the LNA is very accurate in these cases.

The Fano-factor are much larger than 1, as seen in Figure 1B. The reason for these large fluctuations and the strong deviation from Poisson statistics (Fano-factor 1), is high turnover rates of the molecules at steady state in combination with slow relaxation back to steady state. This is exactly what is described in equation 15, where f₁/₁ is the turnover rate of the pool and f'₁(₁) - f'₂(₁) is the rate of relaxation. The high turnover rate creates a lot of chemical "noise" (large diffusion) that decays very slowly (weak restoring force), and this leads to very large fluctuations.

Example 3: Coupled Fluctuations in Multiple Component Systems
Similar principles apply to systems with more than one molecular species; that is, slow relaxations, high turnover rates, and large jumps in molecule numbers (high reaction stoichiometries) cause large fluctuations. The degrees of freedom are, however, larger in multicomponent systems, and intriguing stochastic phenomena of a novel type arise because of correlations between the fluctuations of different components. For instance, fluctuations in one species (e.g., mRNA) can induce fluctuations in another species (e.g., protein). Furthermore, fluctuations can bring the molecule number of one species to a very high value while compensating fluctuations can bring the molecule numbers of another species to a very low value in such a way that the flow properties of the system remain almost unchanged. The latter case is relevant for many biosynthetic reactions, and we will therefore discuss a simple example illustrating this type of degenerate behavior.

Although the Master equation is usually impossible to solve for multicomponent systems, the LNA can be used to estimate the covariance matrix directly from the macroscopic rate equations and the stoichiometry matrix as described in Methods.

To illustrate, we consider an anabolic system, where, respectively, molecules X and Y are synthesized from reservoirs A and B at rates k₁ and k₂. X and Y are assumed to associate irreversibly with association rate constant k₅ in the formation of a heterodimer C. The molecules X and Y also decay with first-order rate constants k₃ and k₄, respectively, as described by the scheme:

(17)

The concentrations of the A, B, and C molecules are assumed to be constant, and the concentrations of X and Y are ₁ and ₂, respectively. This general system is analyzed in the Supplemental material. Here, we focus on the most interesting parameter regime, where X and Y are synthesized at equal rates (k₁ = k₂ = ), and decay with equal and low rates (k₃ = k₄ = µ), that is,

(18)

Low decay rates imply that the degradation events are rare compared to the rate of formation of C molecules, that is, µ << . This, as we shall see, is the condition for anomalously large component number fluctuations.

The macroscopic differential equations for the concentrations and their steady-state values are

When 0 < µ << , the fluctuations may still be very large, as can be seen in the Monte Carlo simulation (Fig. 2C). This is again a close-to-zero-order phenomenon, but now in two dimensions. The fluctuations in Figure 2C are anticorrelated such that the major consumption term, k₅₁₂, is close to constant, that is, both synthesis and consumption rates are near zero order.

View larger version (40K): [in this window] [in a new window]   Figure 2 Properties of a system with coupled exit flow (two-component system, equation 17). (A) The average concentrations (Y-axis) of X molecules for an increasing rate of X synthesis (k1, X-axis). Lines correspond to LNA and circles to estimates from Gillespie simulations. The synthesis rate of Y molecules is k2 = 10-6 M/sec, and the association rate constant is k5 = 106 M/sec. These parameters are also used in B-D. (B) The Fano-factor for X molecules, X2/X, as given by LNA (lines) and Gillespie simulation (circles). (C) Gillespie simulation of the variation of the number of X and Y molecules (Y-axis) over time (X-axis) for balanced synthesis (i.e., k1 = k2 = 10-6 M/sec). (D) The modified Fano-factor, 2w/X1, at the point of balanced synthesis (Y-axis) for different values of the decay rate (X-axis) as given by LNA (lines) and Gillespie simulation (circles).

In the LNA of the two-component system, given by equation 18, the fluctuations are analyzed in the variables w = ₁ - ₂ and u = ₁ + ₂. This change of variables follows automatically with our method (see Supplemental material). The advantage of using these variables is seen by subtracting the differential equations in equation 19 to get the differential equation for w = ₁ - ₂:

(20)

The difference between X and Y decays at a rate set by µ, which is small compared with the turnover rate of the X and Y pools (/). The variance in w(²_w) as calculated by the LNA, normalized to the number of molecules in one of the pools at steady state, is

(21)

This expression is exact, except for the small difference between a true concentration average and a macroscopic stationary state value. (The exact expression can be derived with moment-generating functions.) It can now be seen that it is the large separation between the turnover rate of the X and Y pools (/) and the relaxation rate µ that causes the very large fluctuations in w.

Results from the LNA of the full system (equation 17) over large parameter regions are shown in Figure 2. The LNA is compared with solutions of the Master equation at some points. A major advantage of the LNA is that the stochastic properties of a system can be rapidly scanned over a large parameter space at a very modest computational cost.

Example 4: Elimination of Fast Variables
In Example 3, large fluctuations arise because of one very slow relaxation mode in the coupled system. This behavior is typical of systems near critical points, like phase transitions (Haken 1982). In Methods, we describe a method to make the LNA more accurate in these cases. The method requires that the dynamics of the system is mainly determined by a few slow variables, which are linear combinations of the original concentrations, and that the fast variables will be at steady-state values, conditional on the present values of the slow variables. This approach gives precise estimates of the stochastic system behavior when the rate laws are approximately linear in the slow variables (but not necessarily in the fast variables), and gives a much more accurate description of the system dynamics than direct usage of the LNA.

It should be noted that this technique does not mean that one just assumes rapid equilibration among some of the molecular species. The number of relevant slow variables can, in fact, be lower than such a direct elimination would presuppose. To illustrate the principle, we exemplify with a four-component system, in which the large fluctuations are contained in just one slow variable. The example is an extended version of the system discussed in Example 3 above:

(22)

A and B molecules are synthesized by enzymes E_A and E_B, respectively. A and B can bind E_A and E_B, and thereby inhibit their own synthesis, and can also bind each other irreversibly. Lower-case letters correspond to concentrations. The stationary, marginal distribution P(B) (see Methods for a definition of marginal distribution) for the B molecules is estimated by three different methods and the results are shown in Figure 3. The distribution estimated by the Gillespie simulation (gray line) of the full system (equation 22) has similar mean and variance as the normal distribution resulting from the LNA, but the shape of the latter (black line) is very different from the more accurate distribution obtained from the Gillespie simulations. However, this failure of directly applying the LNA can be overcome by elimination of fast variables (see Supplemental material). Figure 3 shows that the elimination of all but one slow variable and application of the LNA gives a much better representation of P(B) than direct application of the LNA to the four-variable system (dashed line).

View larger version (25K): [in this window] [in a new window]   Figure 3 LNA combined with elimination of fast variables (four-component system, equation 22). The stationary marginal distribution for the number of B molecules, P(B), has been estimated by three different methods: (1) by Gillespie simulation (gray line), which converges toward the correct solution; (2) direct application of LNA (solid black line), which approximates the distribution with a normal distribution; and (3) LNA combined with elimination of fast variables (dashed black line). Parameters: = 10-15l; k1 = 200 sec-1; k2 = 1.67 x 105 M-1 sec-1; k3 = 10,000 sec-1; k4 = 100 sec-1; k5 = 1.67 x 105 M-1 sec-1; k6 = 1000 sec-1; k7 = 1.67 x 106 M-1 sec-1; eAtot = 83 nM;eBtot = 167 nM.

	DISCUSSION

Top ABSTRACT RESULTS DISCUSSION METHODS REFERENCES

We have outlined a direct path (Methods) from macroscopic descriptions of intracellular chemical networks to the linear noise approximations (van Kampen 1997) of the Master equations associated with these systems. The starting point was to formulate the macroscopic dynamics in terms of elementary or elementary complex (Keizer 1987) reactions where the stoichiometries of the different reactions are made explicit. Special attention has been paid to the stochastic properties of chemical networks that are in macroscopically stable stationary states, where the LNA corresponds to a linear Fokker-Planck equation for the movement of a Brownian particle in a multidimensional harmonic potential (Wax 1954). Its solution (Risken 1984) evolves to a multidimensional normal distribution for deviations in molecule numbers from their macroscopic values. The stationary solution of the LNA is completely characterized by its covariance matrix and the Jacobian matrix associated with the macroscopic stationary state. We have used well-known results about the Lyapunov matrix equation (Horn and Johnson 1991) to derive an explicit expression for the covariance matrix in terms of the eigenvalues of the Jacobian matrix, from the macroscopic system at steady state, and the diffusion matrix. The analytical expressions for the cross- and autocorrelation functions of the stationary solution to the linear Fokker-Planck equation of the Master equation have also been described.

This type of analysis facilitates screening of any intracellular network with respect to its stochastic properties. The characteristics of the covariance and Jacobian matrices will decide when fluctuations are so small and fast that a system can be analyzed macroscopically. They will reveal the existence of anomalously large fluctuations that require mesoscopic approaches for precise system descriptions.

Our results show how giant fluctuations in multidimensional systems can arise owing to large terms in the diffusion matrix (large turnover rates of component pools that generate internal noise) and small eigenvalues of the Jacobian matrix (weak restoring forces). When there is a single stable, but weakly attracting, stationary state so that diffusion is almost free in its neighborhood, very large relative fluctuations can occur for systems of finite size. A striking example is the behavior of anabolic reactions in which two or more substrates are independently synthesized and then irreversibly joined by an unsaturated enzyme. When the rates of synthesis of all substrates are balanced, all combinations of substrate concentrations that are compatible with the joint exit flow are macroscopically allowed, and there is free diffusion throughout this continuum of states. In growing systems or in systems in which the rates of synthesis are attenuated by feedback or product inhibition there will be a single stationary state, but it is often so weakly constrained that highly significant fluctuations prevail (Elf et al. 2003b).

Our discussion of the LNA complements existing numerical tools (Sauro 1993; Mendes 1997) for the analysis of stability properties of stationary states and bifurcations (Strogatz 1994) in kinetic reaction networks. These tools normally make use of the Jacobian matrix and its eigenvalues. The explicit formulation of the diffusion matrix from the Liapunov equation now makes it easy to obtain also mesoscopic system properties near stationary points.

Our general analysis, exemplified by concrete examples of how anomalous fluctuations arise in intracellular networks, has also been extended beyond direct applications of the LNA to intracellular systems. When fluctuations are large compared with the range of validity of linear rate laws, then the LNA is expected to fail in many cases. We have analyzed systems in which the macroscopic stationary states are almost degenerate. In those cases, the eigenvalues of the Jacobian are very small, although the flows through the system are significant. Here, large fluctuations are associated with slow relaxation modes, and we have used this system property to systematically eliminate fast variables (Gardiner 1985) with small fluctuations. This approach improves the performance of the LNA in cases in which the slow variables can be described by linear rate laws, although the fast variables display nonlinear dynamics.

We have described several cases in which the LNA or the LNA improved by elimination of fast variables gives accurate mesoscopic system descriptions. However, in cases in which fluctuations bring a system far from the region of validity of the local LNA, a more detailed mesoscopic analysis is required. This is, for instance, the case in systems with noise-induced transitions (Horsthemke and Lefever 1984), which is exemplified (Fig. 4) by how the LNA performs for a system with noise-induced oscillations (Vilar et al. 2002). In the parameter region where an asymptotically stable macroscopic steady state exists, the LNA works, but where periodic oscillations occur because of random fluctuations, the LNA fails to reproduce the oscillations. However, the LNA correctly predicts a marked increase in the size of the fluctuations at the point where the stochastically induced oscillations begin to emerge (Fig. 4).

View larger version (13K): [in this window] [in a new window]   Figure 4 Example of failure of local linearization and LNA. This figure demonstrates the application of LNA to the system described by Vilar et al. (2002). We use the same equations and parameters as in the original article. In a certain parameter regime, the mesoscopic system displays periodic oscillation, whereas the macroscopic model describes an asymptotically stable steady state (A). As can be seen in B, the LNA well describes the Fano-factor for the repressor and repressor-inhibitor complex up to a certain value for the decay rate for the repressor (R 0.0275h-1). Above this point, the system displays large stochastic fluctuations and for even higher values (R > 0.05h-1) also periodic oscillations that cannot be described by the LNA. However, the large fluctuations predicted by LNA indicate that a more thorough stochastic analysis should be conducted. At the bifurcation, at R 0.101h-1, the macroscopic steady state loses stability and the oscillations are described also by a macroscopic model.

	METHODS

Top ABSTRACT RESULTS DISCUSSION METHODS REFERENCES

 
Macroscopic and Mesoscopic Dynamics
Consider an intracellular biochemical system with volume and N different chemical components. The concentrations of the components are summarized in the vector x = (x₁,..., x_N)^T and the numbers of molecules in X = x = (X₁, X₂,..., X_N)^T. The state of the system is defined by X (or x), and a state change takes place by one of R elementary reactions. The probability that the elementary reaction number j will occur in a small time interval t is given by _j(x, )t, where _j(x, ) is a transition rate. By such an event the chemical component number i changes from X_i to X_i + S_ij molecules. The integers S_ij, i = 1, 2,..., N; j = 1, 2,..., R are the elements of the N x R stoichiometric matrix S of the reaction network (Heinrich and Schuster 1996).

When the system volume increases further and further, the relative fluctuations in the network will eventually become insignificant. In the limit of infinitely large , the stochastic concentration vector x will normally become deterministic, and the transition rates _j(x, ) will simplify to their macroscopic rate law counterparts . We define

(23)

Here, is the vector of macroscopic concentrations for which the dynamics is determined by N ordinary rate equations (Cornish-Bowden 1995), according to

(24)

The Master equation (van Kampen 1997) for the system determines the time evolution of the probability distribution P(X, t) for the molecule number vector X and can be written in compact form as

(25)

where E^-Sij is a step operator. When it acts on a function g(X₁,..., X_i,..., X_N), it removes S_ij molecules from X_i, that is, E^-Sijg(X₁,..., X_i,..., X_N) = g(X₁,..., X_i - S_ij,..., X_N) and E^-SijE^-Skjg(..., X_i,..., X_k,...) = g(..., X_i - S_ij,..., X_k - S_kj,...). Equation 24 is the complete macroscopic system description and equation 25 is the full Master equation for the same reaction network. Both equations 24 and 25 assume spatial homogeneity by fast diffusion of all components. When this is not the case, the Master equation has to be extended (Elf et al. 2003a). The stoichiometric matrix S appears in both equations 24 and 25, whereas the mesoscopic rate law _j in equation 25 is replaced by its macroscopic counterparts f_j in equation 24. How the macroscopic equation 24 emerges from the mesoscopic equation 25 when the system volume goes to infinity is shown in the derivation of the linear noise approximation in the Supplemental material (equation A8).

A crucial point in the mesoscopic description is correct identification of the elementary reactions in equations 24 and 25. Elementary means that the transition rates _j(x, ) should be constant as long as x does not change, that is, the rates do not depend on variables not included in the state description. A reaction can be "elementary complex" (Keizer 1987) as when it is catalyzed by an enzyme with internal states that equilibrate on a much shorter time scale than the time scale for product formation. Conversion of a substrate X₁ to a product X₂, carried out by an enzyme with Michaelis-Menten kinetics, is under such conditions elementary complex and can be characterized by the compounded transition rate (x₁)=_maxx₁/(K_m+x₁).

Linearization of Macroscopic Dynamics and the Linear Noise Approximation
In this section we first recapitulate selected parts of the standard treatment (Strogatz 1994) of linearized dynamics around a stationary state and then derive an analytical solution to the linear noise approximation of the Master equation, equation 25. If a stationary state exists for the macroscopic dynamics in equation 24, it follows from the algebraic system of equations:

(26)

Linearization of equation 24 around the stationary solution vector associated with equation 26 leads to a matrix equation for the deviations from according to

(27)

A is the Jacobian matrix, with the elements

(28)

When A can be written as A = QQ^-1, where is a diagonal matrix with the eigenvalues (₁, ₂,..., _N) of A as elements and Q is a matrix with the eigenvectors of A in the columns, the solution to equation 27 is

(29)

where is the initial condition (Strang 1988). The linear transformation , where , leads to a new set of variables (normal modes) with the property that

(30)

The Linear Noise Approximation
The mesoscopic counterpart to the macroscopic solutions in equations 29 or 30, the linear noise approximation (van Kampen 1997), follows from the Master equation, equation 25, by linearizing the rate laws as in equation 27 and by approximating its difference relations by differentials. To our knowledge, the general multidimensional LNA has not been derived in the literature and is therefore given in the Supplemental material. In the LNA, the probability distribution for the vector X, describing random deviations in molecule numbers from their macroscopic values , is obtained from a linear Fokker-Planck equation (Supplemental equation A9). This equation (Risken 1984) depends on the Jacobian A defined in equation 28 and on the noise-generating diffusion matrix BB^T, which is given by

(31)

The matrix is diagonal with the elements ; j = 1, 2,..., R. It depends on the macroscopic rate laws , rather than on the mesoscopic transition rates , because the difference between the two disappears in the LNA (see Supplemental material). The stationary distribution of the Fokker-Planck equation is a multivariate Gaussian with zero mean vector and covariance matrix C = XX^T, which is determined by the (Lyapunov) matrix equation

(32)

where A and BB^T are evaluated at the stationary state . C contains the second-order moments of the stationary Gaussian distribution for the stochastic vector X, according to

(33)

Knowledge of A and C also gives the time correlation matrix for X according to

(34)

Equation 32 can be solved for C numerically with standard methods, but analytical expressions can also be obtained. The properties of C are best explained from the analytical solution that exists for the covariance matrix for the new variables =Q^-1X, that is, fluctuations as normal modes (compare with equation 30). is associated with the linearly transformed stoichiometric matrix = Q^-1S, and it follows from equation 31 that B=Q. The covariance matrix C for X is defined from C = QCQ^* (Horn and Johnson 1985), where Q^* means Q complex conjugated and transposed, and A = QQ^-1. Insertion in equation 32 gives

(35)

The diagonal elements of , that is, the variances , are the summed squares of the transformed stoichiometric coefficients _ij (which can be complex valued) multiplied with their corresponding rate laws (flows), divided by the real part of the rate of relaxation for the normal mode in equation 30:

(36)

It is seen that internal noise is boosted by high rates f_j and large stoichiometry coefficients|_ij| and that fluctuations are attenuated by a large (negative) real part of the eigenvalue of the normal mode in equation 30. When A is symmetric, Q can be made unitary and then = Q^TS, which simplifies analysis and interpretation (see Example 3 in the Supplemental material).

Single-Component Systems
When there is only one chemical species X₁, then equation 32 is scalar and the LNA of the stationary variance becomes

(37)

The Fano-Factor and the Correlation Coefficient
The Fano-factor (Fano 1947) is a measure of the relative size of the internal fluctuations. It is defined as the variance divided by the average value, that is, , which in LNA is approximated by . F_X = 1 when X has Poisson statistics. As both the average value and variance of internal fluctuations are proportional to the volume of the system, the Fano-factor is independent of volume. Another volume-independent stochastic measure is the correlation coefficient , which in the LNA becomes . It is zero for uncorrelated fluctuations and +1 or -1 if fluctuations are perfectly linearly correlated or anticorrelated, respectively.

Elimination of Fast Variables Can Make the LNA More Accurate Near Critical Points
Starting from the full system description in terms of S and , we first identify the correlation matrix from the LNA in equation 35, using equations 28 and 31. We identify the variables (), which correspond to the slow relaxation modes where Re_1...j << Re_{j + 1...N}. The macroscopic dynamics of the system (equation 24) are then formulated in terms of transformed variables similar to equation 30:

(38)

Next, the fast variables are locked to their macroscopic steady-state values, conditional on , by setting the differential equations g_{j + 1}...g_N equal to zero and solving the following system of equations:

(39)

Subsequently, the LNA is applied only to the slow variables to get their stationary probability density function P₁(₁,..., _j). This is a normal distribution centered at the macroscopic steady state with a covariance matrix _{1...j, 1...j}. For each value of (₁,..., _j), the fast variables are locked to their conditional stationary values , and the overall stationary distribution function is therefore:

(40)

where (_k - _hk(₁,..., _j)) is the Dirac function and _i = ^-1_i. () is not a normal distribution any more because the h = (h_{j + 1},..., h_N) functions generally are nonlinear. Finally, we get the probability distribution function in the original variables as

(41)

where is a normalization constant. For a multidimensional P(X), it is convenient to consider only the marginal distributions, defined as

(42)

	Acknowledgements

 
This work was supported by the National Graduate School of Scientific Computing (NGSSC). We thank Martin Lovmar for comments on the manuscript and Lars Elden for pointing out the Horn and Johnson reference.

The publication costs of this article were defrayed in part by payment of page charges. This article must therefore be hereby marked "advertisement" in accordance with 18 USC section 1734 solely to indicate this fact.

	Footnotes

 
¹ Corresponding authors.
E-MAIL Johan.elf{at}icm.uu.se; FAX 46184714262. E-MAIL ehrenberg{at}xray.bmc.uu.se; FAX 46184714262.

Article and publication are at http://www.genome.org/cgi/doi/10.1101/gr.1196503.

[Supplemental material is available online at www.genome.org.]

	REFERENCES

Top ABSTRACT RESULTS DISCUSSION METHODS REFERENCES

 

Barkai, N. and Leibler, S. 2000. Circadian clocks limited by noise. Nature 403: 267-268.[Medline]

Becskei, A. and Serrano, L. 2000. Engineering stability in gene networks by autoregulation. Nature 405: 590-593.[CrossRef][Medline]

Berg, O.G. 1978. A model for the statistical fluctuations of protein numbers in a microbial population. J. Theor. Biol. 71: 587-603.[CrossRef][Medline]

Berg, O.G., Paulsson, J., and Ehrenberg, M. 2000. Fluctuations and quality of control in biological cells: Zero-order ultrasensitivity reinvestigated. Biophys. J. 79: 1228-1236.[Abstract/Free Full Text]

Cornish-Bowden, A. 1995. Fundamentals of enzyme kinetics. Portland Press, London.

Dogterom, M. and Leibler, S. 1993. Physical aspects of the growth and regulation of microtubule structures. Phys. Rev. Lett. 70: 1347-1350.[CrossRef][Medline]

Donachie, W.D. and Robinson, A. 1987. Cell division: Parameter values and the process. In Escherichia coli and Salmonella typhimurium: Cellular and molecular biology (eds. F. Neidhardt et al.), pp. 1578-1593. ASM Press, Washington, DC.

Elf, J., Doncic, A., and Ehrenberg, M. 2003a. Mesoscopic reaction-diffusion in intracellular signaling. In Fluctuations and noise in biological, biophysical and biomedical systems (eds. S. Bezrukov et al.), pp. 114-124. SPIE, Bellingham, WA.

Elf, J., Paulsson, J., Berg, O.G., and Ehrenberg, M. 2003b. Near-critical phenomena in intracellular metabolite pools. Biophys. J. 84: 154-170.[Abstract/Free Full Text]

Elowitz, M.B., Levine, A.J., Siggia, E.D., and Swain, P.S. 2002. Stochastic gene expression in a single cell. Science 297: 1183-1186.[Abstract/Free Full Text]

Fano, U. 1947. Ionization yield of ratios. II. The fluctuations of the number of ions. Phys. Rev. 72: 26-29.[CrossRef]

Fersht, A.R. 1998. Structure and mechanism in protein science. Freeman, New York.

Gardiner, C. 1985. Handbook of stochastic methods. Springer Verlag, Berlin.

Gibson, M. and Bruck, J. 2000. Efficient exact stochastic simulation of chemical systems with many species and channels. J. Phys. Chem. A 104: 1876-1889.[CrossRef]

Gillespie, D. 1976. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22: 403-434.[CrossRef]

____. 1977. Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81: 2340-2361.[CrossRef]

Glad, T. and Ljung, L. 2000. Control theory: Multivariable and nonlinear methods. Taylor & Francis, London.

Goldbeter, A. and Koshland Jr., D.E. 1981. An amplified sensitivity arising from covalent modification in biological systems. Proc. Natl. Acad. Sci. 78: 6840-6844.[Abstract/Free Full Text]

____. 1982. Sensitivity amplification in biochemical systems. Q. Rev. Biophys. 15: 555-591.[Medline]

Guptasarama, P. 1995. Does replication-induced transcription regulate the myriad of low copy number proteins of Escherichia coli? Bioessays 17: 987-997.[CrossRef][Medline]

Haken, H. 1982. Synergetics. Nonequilibrium phase transitions and self-organization in biological systems. In Thermodynamics and kinetics in biological processes (eds. I. Lamprecht and A.I. Zotin). Walter de Gruyter, Berlin.

Heinrich, R. and Schuster, S. 1996. The regulation of cellular systems. Chapman and Hall, New York.

Horn, R.A. and Johnson, C.R. 1985. Matrix analysis. Cambridge University Press, Cambridge, UK.

____. 1991. Topics in matrix analysis. Cambridge University Press, Cambridge, UK.

Horsthemke, W. and Lefever, R. 1984. Noise-induced transitions. Theory and applications in physics, chemistry, and biology. Springer Verlag, Berlin.

Keizer, J. 1987. Statistical thermodynamics of nonequilibrium processes. Springer-Verlag, Berlin.

Kepler, T.B. and Elston, T.C. 2001. Stochasticity in transcriptional regulation: Origins, consequences, and mathematical representations. Biophys. J. 81: 3116-3136.[Abstract/Free Full Text]

Mendes, P. 1997. Biochemistry by numbers: Simulation of biochemical pathways with Gepasi 3. Trends Biochem. Sci. 22: 361-363.[CrossRef][Medline]

Ozbudak, E.M., Thattai, M., Kurtser, I., Grossman, A.D., and van Oudenaarden, A. 2002. Regulation of noise in the expression of a single gene. Nat. Genet. 31: 69-73.[CrossRef][Medline]

Paulsson, J. and Ehrenberg, M. 2001. Noise in a minimal regulatory network: Plasmid copy number control. Q. Rev. Biophys. 34: 1-59.[CrossRef][Medline]

Qian, H., Saffarian, S., and Elson, E.L. 2002. Concentration fluctuations in a mesoscopic oscillating chemical reaction system. Proc. Natl. Acad. Sci. 99: 10376-10381.[Abstract/Free Full Text]

Risken, H. 1984. The Fokker-Planck equation. Spinger Verlag, Berlin.

Sauro, H.M. 1993. SCAMP: A general-purpose simulator and metabolic control analysis program. CABIOS 9: 441-450.

Savageau, M.A. 1976. Biochemical systems analysis: A study of function and design in molecular biology. Addison-Wesley, Reading, MA.

Schrödinger, E. 1944. What is life? Cambridge University Press, Cambridge, UK.

Strang, G. 1988. Linear algebra and its applications. Harcourt Brace Jovanovich Collage Publishers, Fort Worth, TX.

Strogatz, S.H. 1994. Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering. Addison-Wesley, Cambridge, MA.

Stundzia, A.B. and Lumsden, C.J. 1996. Stochastic simulation of coupled reaction-diffusion processes. J. Comput. Phys. 127: 196-207.[CrossRef]

Thattai, M. and van Oudenaarden, A. 2001. Intrinsic noise in gene regulatory networks. Proc. Natl. Acad. Sci. 98: 8614-8619.[Abstract/Free Full Text]

van Kampen, N.G. 1997. Stochastic processes in physics and chemistry, 2nd ed. Elsevier, Amsterdam.

Vilar, J.M., Kueh, H.Y., Barkai, N., and Leibler, S. 2002. Mechanisms of noise-resistance in genetic oscillators. Proc. Natl. Acad. Sci. 99: 5988-5992.[Abstract/Free Full Text]

Wax, N. 1954. Selected papers on noise and stochastic processes. Dover, New York.

Received January 20, 2003; accepted in revised format July 22, 2003.
CiteULike    Connotea