Article Information

• PDF (206 kb)
• Full Text with Large Figures
• Cited by in Scopus (4)

PubMed

• Articles by Matthew R. Bennett
• Articles by Jeff Hasty

Related Articles

• A Spatial Toggle Switch Drives Boundary Formation in Develop...

  A Spatial Toggle Switch Drives Boundary Formation in Development Biophysical Journal, Volume 95, Issue 11, 1 December 2008, Pages 5111-5120 Oriol Canela-Xandri, Francesc Sagués, Ramón Reigada and Javier Buceta Abstract Herein we introduce a multicellular network motif that performs as a spatial toggle switch and explains how boundary formation can be faithfully accomplished in developmental processes. Importantly, we show that expression and activity patterns of proteins must be simultaneously characterized for a proper understanding and description of the underlying mechanism. Our in silico experiments, in agreement with in vivo results, evaluate different genetic backgrounds and shed light on the dynamics of boundary formation. In addition, we provide an estimation of relevant biological parameters and a robustness analysis. Abstract | Full Text | PDF (670 kb)

• Collective Decision Making in Bacterial Viruses

  Collective Decision Making in Bacterial Viruses Biophysical Journal, Volume 95, Issue 6, 15 September 2008, Pages 2673-2680 Joshua S. Weitz, Yuriy Mileyko, Richard I. Joh and Eberhard O. Voit Abstract For many bacterial viruses, the choice of whether to kill host cells or enter a latent state depends on the multiplicity of coinfection. Here, we present a mathematical theory of how bacterial viruses can make collective decisions concerning the fate of infected cells. We base our theory on mechanistic models of gene regulatory dynamics. Unlike most previous work, we treat the copy number of viral genes as variable. Increasing the viral copy number increases the rate of transcription of viral mRNAs. When viral regulation of cell fate includes nonlinear feedback loops, very small changes in transcriptional rates can lead to dramatic changes in steady-state gene expression. Hence, we prove that deterministic decisions can be reached, e.g., lysis or latency, depending on the cellular multiplicity of infection within a broad class of gene regulatory models of viral decision-making. Comparisons of a parameterized version of the model with molecular studies of the decision structure in the temperate bacteriophage are consistent with our conclusions. Because the model is general, it suggests that bacterial viruses can respond adaptively to changes in population dynamics, and that features of collective decision-making in viruses are evolvable life history traits. Abstract | Full Text | PDF (176 kb)

• Kinetics from Nonequilibrium Single-Molecule Pulling Experim...

  Kinetics from Nonequilibrium Single-Molecule Pulling Experiments Biophysical Journal, Volume 85, Issue 1, 1 July 2003, Pages 5-15 Gerhard Hummer and Attila Szabo Abstract Mechanical forces exerted by laser tweezers or atomic force microscopes can be used to drive rare transitions in single molecules, such as unfolding of a protein or dissociation of a ligand. The phenomenological description of pulling experiments based on Bell’s expression for the force-induced rupture rate is found to be inadequate when tested against computer simulations of a simple microscopic model of the dynamics. We introduce a new approach of comparable complexity to extract more accurate kinetic information about the molecular events from pulling experiments. Our procedure is based on the analysis of a simple stochastic model of pulling with a harmonic spring and encompasses the phenomenological approach, reducing to it in the appropriate limit. Our approach is tested against computer simulations of a multimodule titin model with anharmonic linkers and then an illustrative application is made to the forced unfolding of I27 subunits of the protein titin. Our procedure to extract kinetic information from pulling experiments is simple to implement and should prove useful in the analysis of experiments on a variety of systems. Abstract | Full Text | PDF (368 kb)

• …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 10, 3501-3512, 15 May 2007

doi:10.1529/biophysj.106.095638

Biophysical Theory and Modeling

Transient Dynamics of Genetic Regulatory Networks

Matthew R. Bennett^*, †, Dmitri Volfson^*, †, Lev Tsimring^* and Jeff Hasty^*, †, ,  

^* Institute for Nonlinear Science, La Jolla, California
^† Department of Bioengineering, University of California at San Diego, La Jolla, California

Address reprint requests to J. Hasty.

As the complexity of gene regulatory networks under study increases so does the need for accurate modeling techniques 1. While exact numerical simulations are possible using Monte Carlo techniques like the Gillespie algorithm 2, such simulations can be computationally intensive. Continuous approximation schemes based on the underlying stoichiometric reactions can be used to simulate the dynamics of the average of each species in the system, but the complexity of these models can hinder both computational and theoretical analysis. Hence, many theorists have resorted to using a quasi-steady-state approximation (QSSA) 3,4 to reduce the number of dimensions in continuous models. Such reduced models do an excellent job in many cases, especially when the asymptotic state of the system is a stable equilibrium point. However, it has long been known that in some cases the QSSA does a poor job predicting the timescales over which systems equilibrate to their steady-state value 5,6,7. Moreover, interest is continually growing in gene networks that exhibit more complicated behavior, like stable limit cycles 8,9,10,11. When periodic behavior is present in a system, correct prediction of the timescales involved becomes necessary for a complete understanding of the system, and is essential for guiding experimental studies.

In this article, we present a continuous approximation scheme that reduces the number of dimensions in the system, while at the same time predicts the timescales of the full system more accurately than does the classic QSSA. By correctly applying multiple timescale analysis, the resulting reduced systems are the same as QSS approximations, but with a prefactor in front of the time derivatives of the concentrations. This method was first introduced by Kepler and Elston 6, and Bundschuh et al. 12 showed that the prefactor derived by Kepler and Elston was related to the Jacobian of a transformation relating monomer concentrations to the total concentration of a particular protein. We will examine this method in more detail, and demonstrate that the prefactors appear generally as a result of correct reduction of multiscale dynamics to a slow manifold in which fast dynamics are assumed to be in a local quasi-equilibrium.

We first introduce the approximation scheme by examining a simple example—the genetic feedback loop. We show that while both the prefactor method and the QSSA correctly predict the asymptotic behavior of the system, the transient dynamics are better modeled when the prefactor is included. Next, we look at a system with a stable limit cycle—the repressilator 8. While the QSSA correctly predicts oscillations in this system, we show that it incorrectly predicts the amplitude and frequency of those oscillations, which are better estimated with the new technique. Furthermore, the bifurcation between stable fixed points and limit cycles is more accurately estimated with the new scheme. In the last section, we will look at the lysis-lysogeny switch of bacteriophage λ. The temporal dynamics of the switch are an important aspect of its function, because the speed with which it makes transitions between the stable states will act as a limiting factor on the timescales at which the lytic cycle may be controlled. Therefore, reduced models describing the dynamics of the phage λ switch (or any switch) must accurately reproduce these timescales.

Separation of timescales is very common in the dynamics of gene regulatory networks. Many processes, like dimer- and tetramerization, occur at a much faster rate than other processes, such as transcription, translation, and degradation 4,13. Because the fast reactions are quick to equilibrate, it has been a common practice by many theorists to use a quasi-steady-state approximation to replace dynamical variables involved in these reactions with their steady-state values. This approximation reduces the number of dimensions in the resulting systems of ODEs, and hence greatly simplifies ensuing analysis.

Consider, for instance, the genetic positive feedback loop 4,14. This system involves a single gene that transcribes a single protein. Upon dimerization, the protein dimer can bind to an upstream regulatory site and stimulate transcription. Additionally, the protein monomer and the mRNA produced in transcription are subject to degradation. This situation is depicted in Fig. 1.

Display large version of this figure

Figure 1

A schematic of the genetic positive feedback loop. The protein monomers (A) bind into dimers (B) that subsequently bind to the upstream regulatory site (C), activating production of the monomer.

In this system there are nine reactions occurring among five chemical species. A list of these reactions is given in Table 1. Here x and y are the concentrations of protein monomers and dimers, respectively; d₀ is the concentration of promoter sites that are free of the dimer; d_r is the concentration of promoter sites that are bound to the protein dimer; and m is the concentration of mRNA strands.

Table 1 Reactions in the genetic positive feedback network

	Number		Reaction
	1	x+x		y
	2	y		x+x
	3	y+d0		dr
	4	dr		y+d0
	5	d0		d0+m
	6	dr		dr+m
	7	m		m+x
	8	x		Ø
	9	m		Ø

The first four reactions in Table 1 represent the dimerization of the protein, the binding of the dimer to the upstream regulatory site, and their reverse processes. Reactions 5–7 represent transcription and translation, while the last two reactions are the degradation of the protein monomers and the mRNA. Reactions 1–4 typically occur at a much faster timescale than reactions 5–9.

From the reactions given in Table 1, we can write down a system of differential equations that represent the time evolution of the average concentration of each species. These equations are

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

The problem with this method of reduction is in treating x as a slow variable. It is true that depends on slow reactions (namely translation and degradation), but it also depends on two fast reactions (dimerization and dissociation). Therefore, x is not a slow variable, but a mixture of both slow and fast. While rigorous multiple timescale analysis is possible for systems such as these, it is not always necessary. In many cases, the variables can be transformed into “slow” and “fast” variables, so that the differential equations for the transformed variables contain either slow or fast reaction terms, but not both. Note that this does not imply that there are only two timescales in the problem, since there are typically many timescales involved in such networks. Instead, by “slow” and “fast” we mean a partitioning of all timescales into two sets that are separated by at least an order of magnitude. When this occurs, the classic QSSA can be applied to the transformed system. Afterwards, the system must be transformed back into the original variables. For systems in which suitable transformations are not forthcoming, we provide in the Appendix a rigorous reduction scheme for a class of gene networks.

For the positive feedback loop in question, we can make the proper transformation by noting that both dimerization and dissociation keep the total number of protein molecules constant, while translation and degradation do not. Therefore, we can track the truly slow variable n_x=x+2y+2d_r, representing the total concentration of protein molecules (in any form), and write

(11)

(12)

(13)

(14)

Combining Eq. (11) with Eq. (14), we arrive at a new system of equations for the time evolution of x and m, namely,

(15)

(16)

Fig. 2 shows the dynamics of the protein monomer number in both the QSSA (dashed curve) and corrected version with the prefactor (dash-dot curve). Also shown is the result of integrating the full system, without any dimensional reduction (solid curve). The reduced system with the prefactor does a much better job in predicting the correct timescale over which the system relaxes to the fixed point. Furthermore, the prefactor method estimates the relaxation time of the positive feedback loop throughout a wider range of parameter values more accurately than does the QSSA. Let t_1/2 be the time it takes the system to go from a zero concentration of proteins and mRNAs to the time at which the monomer concentration comes to one-half its steady-state value. Fig. 3 shows a comparison of the values of t_1/2 for the QSSA (dashed curve), prefactor method (solid curve), and the nonreduced system (circles) for a range of monomer degradation rates, γ_p. The prefactor method predicts t_1/2 more accurately than does the QSSA. Note that t_1/2 is not a monotonic function of the degradation rate. This is due to our definition of t_1/2, and the existence of nonlinear positive feedback in the system. As the degradation rate is changed, so too is the position of the closest stable fixed point relative to our choice of initial conditions. Furthermore, because the relaxation is nonlinear, a true “half-life” for the relaxation can only be obtained very near the fixed points, making alternative measures of the relaxation necessary.

Display large version of this figure

Figure 2

A comparison between the QSSA (dashed curve) and the ODE system with the prefactor correction (dash-dot curve) for a positive feedback gene network. Also shown (solid curve) is the nonreduced system of ODEs, Eqs. (1). Here κ₊=k₊=50, κ_–=k_–=1000, α=1, β=10, σ=3, γ_p=1, γ_m=6, and n_c=20. The initial conditions are x(0)=10, x₂(0)=0, m(0)=0, and d₀(0)=d.

Display large version of this figure

Figure 3

The relaxation times (t_1/2) as a function of γ_p for the QSSA (dashed curve), prefactor method (solid curve), and the nonreduced system (circles). Note that the values of t_1/2 predicted by the prefactor method are nearly identical with the true values obtained from the nonreduced system. Here κ₊=k₊=50, κ_–=k_–=1000, α=1, β=10, σ=3, γ_m=6, and d=20. The initial conditions are x(0)=0, x₂(0)=0, m(0)=0, and d₀(0)=d.

The correct estimation of the transient dynamics of genetic feedback loops may be of little consequence, since in most cases we are only concerned with the final state of the system. However, for such systems as oscillators, the final state will not be a stable fixed point. If accurate estimation of the timescales of the oscillations (i.e., periods) is desired, we cannot use the QSSA, since it incorrectly predicts these. To examine this issue, let us do another example, that of the repressilator system 8,16, as shown in Fig. 4. The repressilator is a three-element gene network based on a ring architecture in which each of the elements represses the next one down the line. In other words, gene G₁ (after transcription and translation of the resulting mRNA) produces protein x₁, which upon dimerization inhibits transcription of gene G₂. Similarly, the protein dimer y₂ represses the gene G₃, whose protein product, y₃, represses transcription of G₁.

Display large version of this figure

Figure 4

A schematic of the repressilator system. Gene G₁ produces protein x₁, which upon dimerization inhibits transcription of gene G₂. Similarly, the protein dimer y₂ represses the gene G₃, whose protein product, y₃, represses transcription of G₁.

The reactions of the repressilator are shown in Table 2. The first two equations represent the dimerization and dissociation of the protein monomers (x_i) and dimers (y_i). Reactions 3 and 4 are the binding and dissociation of the protein dimers to/from the free (d_0,i) and repressed (d_r,i) promoter sites. If the promoter of gene G_i is free (unrepressed) it can transcribe its associated mRNA (m_i), which in turn can translate its associated protein (reactions 5 and 6). Furthermore, the protein monomers and mRNAs will degrade with time (reactions 7 and 8).

Table 2 Reactions in the repressilator; here i ∈ {1, 2, 3} while k ∈ {3, 1, 2}

	Number		Reaction
	1	xi+xi		yi
	2	yi		xi+xi
	3	d0,i+yk		dr,i
	4	dr,i		d0,i+yk
	5	d0,i		d0,i+mi
	6	mi		mi+xi
	7	xi		Ø
	8	mi		Ø

As before, we can use the reactions given in Table 2 to write a system of ODEs that represent the average concentrations of each species in the thermodynamic limit. These equations are

(17)

(18)

(19)

(20)

(21)

As with the feedback system, we assume that dimerization and dissociation of the proteins (to themselves and the promoters) are fast compared to the other processes. If we take these reactions to be in equilibrium, we find that

(22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

(32)

When the prefactor, p(x_i), is set to unity, Eqs. (30) are exactly the equations produced by the QSSA for the repressilator system 8. Both sets of equations still have regions of stable limit cycles, but the prefactor version more accurately reproduces the full dynamics of the system. This can be seen in Fig. 5. Figure 5a shows the results of a simulation of the full, unreduced equations, Eqs. (34), while Figure 5bc, are from the QSSA and the prefactor method, respectively. Notice that the prefactor method does a much better job in estimating both the period and the amplitude of the oscillations. This can also be seen in Figure 6ab, which show the period and amplitude of each of the three systems as a function of γ_p.

Display large version of this figure

Figure 5

An example of the oscillations seen in the repressilator system for (a) the unreduced equations, Eqs. (34), (b) the QSSA—i.e., Eqs. (30) with p(x_i)=1, and (c) the prefactor method approximation, Eqs. (30). Here γ_p=6, γ_m=1, κ₊=k₊=5, κ_–=k_–=100, α=10, σ=20, and d=20. The initial conditions were x_i=10, y_i=0, d_i=d, and m_i=0.

Display large version of this figure

Figure 6

A comparison of the period (a) and amplitude (b) of the oscillations as a function of the protein degradation rate in the repressilator for the unreduced equations (black curves), the prefactor method (circles), and the QSSA (dashed curves). Here, we calculate the amplitude by max(x₁)-min(x₁) on the limit cycle. The parameters for these figures are the same as those in Fig. 5.

To guide experiments, it is necessary to estimate the regions in parameter space where one should expect stable limit cycles. The prefactor, Eq. (32), contains the parameter r which does not appear in the QSSA and as the value of r is changed (which amounts to changing the ratio of c_d to c_p while keeping their product fixed), the bifurcation curve separating stable fixed points from stable limit cycles (in κ–β space) changes. This is shown in Figure 7a. The solid line represents the bifurcation between stable fixed points (to the left of the curve) and stable limit cycles (to the right of the curve) according to the QSSA. As long as the product c_dc_p remains fixed, the bifurcation curve will not change, no matter what the ratio is. This is in contrast to the predictions of the prefactor method. While keeping all the parameters the same as in the QSSA, we change the parameter r to see how the bifurcation curve changes. The colored curves in Figure 7a represent the bifurcation curve for r=10 (dash-dot), r=1 (dotted), and r=0.1 (dashed). As r→0, the prefactor approximation approaches the QSSA, and as r increases, there is a significant divergence between the regions in parameter space in which stable limit cycles are present for the two methods. Because the timescales of the full system are more accurately predicted with the prefactor method, it does a better job in estimating the region of oscillations than does the QSSA.

Display large version of this figure

Figure 7

Comparison of the parameter space for each model. (a) Bifurcation diagram for the repressilator. The solid line represents the bifurcation curve of the QSSA between stable fixed points (to the left of the curve) and stable limit cycles (to the right of the curve). The remaining curves are the same bifurcation curves for the prefactor method with r=10 (dash-dot), r=1 (dotted), and r=0.1 (dashed). For each of the prefactor curves the parameters are such that, while r differs, the corresponding parameters for the QSSA are the same. The curves for the full system virtually coincide with the curves for prefactor model (μ=ν=10⁴) and are not shown for clarity. For all curves, d′=1. (b) Comparison of the Hopf curves for the prefactor model with the results for the full system for different values of μ=ν. The difference becomes significant as μ and ν are decreased and the assumption of the separation of timescales breaks down. All curves are for r=10 and d′=1. For smaller values of r and/or d′, the instability threshold predicted with the prefactor model converges to the threshold for the full system at smaller values of μ (not shown). (c) Comparison of the Hopf curves for the prefactor model computed numerically for r=0.1 and r=10 (for d′=1, solid lines) with the asymptotic approximations (dashed lines), Eqs. (39). (d) The frequency of the unstable mode along the threshold as a function of β. Circles and squares correspond to numerical results for the full system (μ=ν=10⁴) and for the prefactor method, respectively. The solid line corresponds to the function ω_H specified by Eq. (43). The dashed line predicted by QSSA systematically overestimates the frequency of the oscillatory instability. All curves are for r=1 and d′=1.

To compare these results with the full model, we rescale the additional variables that do not appear in the reduced systems. If we use the same rescalings as before, the dimensionless form the full system for repressilator becomes

(34)

(35)

(36)

(37)

When μ=ν is large enough, it becomes possible to predict the onset of the instability leading to the limit cycle. The system of equations, Eqs. (30), always has a symmetric equilibrium parameterized by which is the unique real solution of the equation

(38)

(39)

(40)

(41)

(42)

The Hopf frequency may be expressed in a compact form, for both the prefactor method and the QSSA, with result

(43)

One important class of gene networks is the toggle switch 17. These networks are designed to respond to an external signal in such a way that they are either on (fully induced) or off (no transcription). A quantitative understanding of the dynamics of gene switches is a crucial first step toward a modular description of gene regulation. In addition, the ability to rapidly switch between multiple stable states is important to the development of sophisticated cellular control schemes. Nonlinearities giving rise to two stable states suggest the possibility of using these states as digital signals to be processed in cellular-level computations 18,19. One may eventually be able to produce systems in which sequences of such switching events are combined to control gene expression in complex ways. In any such application, the speed with which systems make transitions between their stable states will act as a limiting factor on the timescales at which cellular events may be controlled. With this in mind, it becomes important that mathematical models of gene switches correctly predict the timescales over which they travel from one state to the other.

The genetic network of λ-phage switches its host bacterium from the dormant lysogenous state to the lytic growth state in ∼20min 13,20. As shown in Fig. 8, it consists of two promoters, P_R and P_RM, which share three operator sites, O_R1, O_R2, and O_R3. The product of the left promoter (P_RM) is the protein cI, which (upon dimerization) can bind to the promoter sites. When cI is bound, its purpose is twofold. First, it activates transcription of the left promoter, causing a positive feedback loop. Second, it represses the right promoter (P_R), blocking its transcription. The system will remain in this state until some external signal (such as UV radiation) causes the rapid degradation of cI. This is done by an activated form of the protein RecA, which cleaves cI monomers, rendering them incapable of dimerizing. At this point, the concentration of cI becomes low enough to free up the promoter sites, releasing the P_R so that it may produce its protein, Cro. Cro can then bind to the promoter sites, repressing the production of cI.

Display large version of this figure

Figure 8

Genetic network diagram of the phage λ switch. The left promoter produces cI, which activates itself and represses the right promoter. The right promoter produces Cro, which represses the left promoter. Both of the proteins, cI and Cro, can bind to the three promoter sites once they have dimerized.

Transcription of repressor (cI) takes place when there is no protein (of either type) bound to O_R3. When repressor is bound to O_R2, the rate of repressor transcription is enhanced, and Cro is transcribed only when O_R3 is either vacant, or has a Cro dimer bound to it. If either repressor or Cro is bound to either O_R1 or O_R2, the production of Cro is halted. For brevity, we omit the full derivation of the equations of motion, and ignore the intermediate step of mRNA translation. This step can easily be put into the model, but it will not affect the timescale analysis, or the prefactor. For a listing of the relevant chemical reactions in the phage λ-switch and a derivation of the equations without the prefactor, we refer the interested reader to Hasty et al. 4. Letting x and y represent the concentrations of cI and Cro, the competition for operator sites leads to deterministic rate equations of the form

(44)

(45)

(46)

(47)

(48)

(49)

Fig. 9 shows simulations of Eqs. (44), with and without the prefactor. Before irradiation with UV light, cI is in its high state, while the concentration of Cro is near zero. At time t=0 the irradiation occurs, and the degradation rate of cI is increased. Subsequently, the concentration of cI becomes very low, while the concentration of Cro becomes very high, beginning the lytic process in phage λ. Notice that the timescales over which the switch changes states is much longer when the prefactor is included.

Display large version of this figure

Figure 9

Dynamics of the phage λ switch. Here, the parameters are m=20, α=11, ρ_y=62.92, σ₁=2, σ₂=0.8, β₁=β₂=β₃=0.08, β₄=β₅=1, K₁=K₃=5×10⁷ M⁻¹, K₂=K₄=3.3×10⁸ M⁻¹, and γ_y=0.03min⁻¹. For t<0, γ_x=0.03min⁻¹. At time t=0 the degradation rate of cI is increased to γ_x=30min⁻¹, representing cleavage of cI by the activated form of RecA due to an increase in UV radiation.

We have introduced a method for approximating and reducing the full systems of ODEs for genetic regulatory networks. Like the QSSA, the prefactor method simplifies the analysis of such systems by reducing the number of dimensions. Unlike the QSSA, however, the timescales of the system are preserved. Timescale problems have long been known to exist when using the QSSA to derive Michaelis-Menten type reaction equations for enzyme kinetics 21,22,23,24. It is not surprising, then, that timescale problems arise when using the QSSA to derive reduced equations for genetic regulatory networks.

The timescale problems can be overcome with the correct use of timescale separation. While we have shown that a rigorous approach can be derived, one based on multiple timescale analysis, we have also shown that one can derive prefactors by plugging the equilibrium values of the “fast” reactants into differential equations for the slow variables. Generally, the slow variables are not the protein concentrations, but instead are total concentrations of all versions of the specific protein (free or bound into dimers or larger complexes), since the total number of proteins in any form is affected only by the slow reactions of translation and degradation.

In general, the inclusion of the prefactors into the analysis slows down the resulting dynamics of the reduced system. Whereas standard QSSA generally has faster dynamics than the nonreduced system, the prefactor method estimates the timescales very accurately. This accuracy is of paramount importance when modeling gene networks designed to exhibit temporal dynamics 25. For instance, the time that it takes for a genetic switch 26,27 to move from one stable fixed point to another can have important consequences. Furthermore, when designing genetic oscillators, it becomes necessary to correctly predict the regions in parameter space in which oscillations are expected. As we saw with the repressilator, the prefactor method predicts these regions more accurately than does the QSSA. In many cases, the separation of timescales is the key to nontrivial behavior, and correctly modeling such systems is therefore necessary.

It should be noted that in this article we have ignored the stochastic fluctuations due to the intrinsic randomness of the underlying biochemical reactions 2,6,28. These fluctuations can have profound effects on networks, and the timescales of their correlations may determine the nature of their influence 29. A procedure similar to the one given above that does include stochasticity involves the reduction of the full multidimensional master equation onto a manifold of slowly evolving variables 6,12. One can then derive corresponding Langevin equations for protein concentrations that will have similar prefactor terms to the ones derived here. However, these equations ignore the stochasticity inherent in the fast reactions. It is still an open question, then, how to project a full system of Langevin equations (containing both slow and fast reactions) onto a slow manifold while preserving the stochastic effects of the fast reactions. We are planning to address this issue in our future work.

Appendix

General method of timescale analysis for gene networks

While the procedures given above for reducing systems of ODEs for gene regulatory networks are illustrative, a more rigorous approach may be desired. To this end, we present in this Appendix a method based on multiple timescale analysis 30 for projecting the full dynamics onto a slow manifold for a general regulatory network.

Let us consider a general framework for a genetic regulatory network. The notations for the components of the network are listed in Table A1. We assume that the network consists of N genes each with a concentration d_i (i=1, …, N). These genes transcribe N corresponding mRNAs, denoted by m_i. The mRNAs, in turn, are translated at a rate σ_i to produce n_i protein monomers, x_i. Each protein can form an associated dimer, y_i, with rates κ_i and κ_–i for dimerization and dissociation, respectively. Accompanying each gene is a promoter that may contain one or more binding sites to which the protein dimers may bind, and thereby enhance or diminish the transcription rate of the gene.

Table A1 Notation for the generalized genetic network model

	Symbol	Description
	di	Gene concentrations.
	xi	Monomer concentrations.
	yi	Dimer concentrations.
	mi	mRNA concentrations.
	σi	Translation rates.
	γp,i	Degradation rates of the monomers.
	γm,i	Degradation rates of the mRNA.
	pij	Rate at which protein i is degraded by protein j.
		Unoccupied binding site .
	Di	Matrix of occupied binding site concentrations for the ith gene.
		The jth, nth element of Di (binding site occupied by protein n).
	ni	Number of protein molecules translated per one mRNA.
		Label for the jth binding site of the ith gene.
	N	Number of genes.
	Mi	Number of binding sites of the ith gene.
	κ±i	Dimerization/dissociation rates of the proteins to/from themselves.
	αi	Transcription rate of the ith gene.
	K±i	Rate matrix for the binding/dissociation of the dimers to the promoter sites.
		The jth, nth element of K±i (binding/dissociation rate of the nth dimer to the jth binding site of the ith gene).
	Hi	The ratio κi/κ–i.
		The ratio .

We assume that the promoter for the i^th gene contains M_i binding sites (j=1, …, M_i) and any protein can bind to any of the M_i binding sites of the i^th promoter. We denote by the concentration of unoccupied binding sites and by the concentration of the binding sites occupied by n^th protein. Note that the sum

(A1)

is the constant concentration of the i^th gene. In what follows, if we refer to we will be referring to bound promoter sites (n>0) unless we explicitly refer to or unless it is used in a sum such as .

The transcription rate of the i^th gene is α_i, which is a function of the current state of the promoter (i.e., which regulatory proteins are bound to which promoter sites). In general, we can include this dependency by making the transcription rates functions of the bound promoter site concentrations, i.e., α_i=α_i(D_i), where D_i is an M_i×N matrix with elements As written, α_i(D_i) can be thought of as the transcriptional rate of the i^th gene, averaged over each copy of that gene. In other words, we can write

(A2)

where is an M_i-tuple of integers representing a possible configuration of a gene, consisting of protein ν_k occupying the k^th promoter site, and 〈…〉_si represents averaging over all possible states of the promoter. If ν_k=0, then that particular binding site is free. The binding/dissociation rates of the proteins to the i^th promoter are described by two M_i×N matrices, K_i and K_–i, with components and respectively. We will assume that there is no cooperativity between the binding of the proteins to the promoter sites, so that K_±i is not a function of D_i. This assumption is made to simplify the analysis, but is not necessary. We also include into the consideration protein-protein interactions: protein j may play a role of protease for the protein i, with p_ij denoting the corresponding degradation rate. Finally, both the protein monomers and the mRNA will degrade at rates γ_p,i and γ_m,i, respectively. The chemical reactions occurring in the network are listed in Table A2.

Table A2 Reactions in the generalized genetic network model

	NumberReaction			Rate
	1	xi+xi		yi	κi
	2	yi		yi+yi	κ−iyi
	3	xi		Ø	γp, ixi
	4	xi+xj		xj	pijxixj
	5	di		di+mi	αi(Di)di
	6	+yn
	7
	8	mi		mi+nixi	σimi
	9	mi		Ø	γm,imi

From the reactions given in Table A2 we can write down a system of ODEs corresponding to the time evolution of the concentrations of each of the chemical species. As a result, we obtain

(A3)

(A4)

(A5)

(A6)

(A7)

Because of conservation of the total concentration of the genes, Eq. (A6) can be eliminated, and replaced by giving us

(A8)

(A9)

(A10)

(A11)

Among the reactions in Eqs. (A8), dimerization and protein binding/dissociation reactions are usually fast, while transcription, translation, and degradation are slow. This is not to say that there are only two timescales present in the system. On the contrary, gene networks can possess a diverse spectrum of timescales. However, these timescales can usually be partitioned into two classes—those that occur on fast timescales, and those that are slow in comparison. When such a partitioning occurs, it becomes possible to characterize the ratio of the characteristic time constants of fast and slow reactions by the small parameter ϵ, and introduce scaled kinetic constants for the slow reactions, and In the spirit of multiple timescale analysis we introduce the fast and slow times t and T=ϵt, respectively. We assume that all concentrations are functions of these two independent variables and expand them in a power series in the small parameter, ϵ,

(A12)

where z=z(t, T) stands for variables x_i, y_i, m_i, and We replace the time derivatives in Eqs. (A8) by ∂/∂t+ϵ∂/∂T. Then, collecting terms of equal power of ϵ we obtain:

(A13)

(A14)

(A15)