Biophysical Journal 86:1922-1945 (2004)
© 2004 The Biophysical Society
A Nonlinear Discrete Dynamical Model for Transcriptional Regulation: Construction and Properties
John Goutsias * and
Seungchan Kim 
* The Whitaker Biomedical Engineering Institute, The Johns Hopkins University, Baltimore, Maryland 21218; and Translational Genomics Research Institute, Phoenix, Arizona 85004
Correspondence: Address reprint requests to John Goutsias, E-mail: goutsias{at}jhu.edu.
 |
ABSTRACT
|
|---|
Transcriptional regulation is a fundamental mechanism of living cells, which allows them to determine their actions and properties, by selectively choosing which proteins to express and by dynamically controlling the amounts of those proteins. In this article, we revisit the problem of mathematically modeling transcriptional regulation. First, we adopt a biologically motivated continuous model for gene transcription and mRNA translation, based on first-order rate equations, coupled with a set of nonlinear equations that model cis-regulation. Then, we view the processes of transcription and translation as being discrete, which, together with the need to use computational techniques for large-scale analysis and simulation, motivates us to model transcriptional regulation by means of a nonlinear discrete dynamical system. Classical arguments from chemical kinetics allow us to specify the nonlinearities underlying cis-regulation and to include both activators and repressors as well as the notion of regulatory modules in our formulation. We show that the steady-state behavior of the proposed discrete dynamical system is identical to that of the continuous model. We discuss several aspects of our model, related to homeostatic and epigenetic regulation as well as to Boolean networks, and elaborate on their significance. Simulations of transcriptional regulation of a hypothetical metabolic pathway illustrate several properties of our model, and demonstrate that a nonlinear discrete dynamical system may be effectively used to model transcriptional regulation in a biologically relevant way.
|
INTRODUCTION
|
|---|
An emerging theme in modern biology is the development of accurate experimental techniques for monitoring cellular behavior (e.g., see Schena et al., 1996
; Brown and Botstein, 1999; Turner and Varshavsky, 2000; Zhu et al., 2001; Baldi and Hatfield, 2002). Although current techniques are mostly used to identify molecular markers for certain types of disease (e.g., cancer; see Golub et al., 1999; Bittner et al., 2000; Kobayashi et al., 2003b), it is the monitoring and modeling of cellular behavior that could mostly benefit from them.
An important cellular process under investigation is transcriptional regulation. Understanding the biological mechanisms underlying transcriptional regulation may lead to significant advances in cell biology, drug development, and medicine. It is becoming increasingly clear that, to enrich our knowledge about transcriptional regulation and understand the role it plays in cellular function, we need to construct a sufficiently predictive mathematical model for such a process, derived from basic biological principles. Moreover, experimental and computational techniques should be developed to estimate the underlying structure of the model and its parameters. Model simplicity, via reasonable biological assumptions and approximations, is important, due to limited biological knowledge of the mechanisms underlying transcriptional regulation, and difficulties of current technologies in measuring underlying parameters. If the model is sufficiently predictive, we may use it as a computational tool (even in the absence of exact parameter values) to simulate biological scenarios (e.g., steady-state analysis, mutation effects, knock-out studies, perturbation effects, homeostatic and epigenetic regulation, etc.), and generate hypotheses pertaining to the mechanisms underlying transcriptional regulation and control. This plan seems to be easier, faster, and cheaper to implement in silice (i.e., on a digital computer by simulation) than in vivo or in vitro.
There have been considerable efforts to build models for transcriptional regulation (e.g., see Thomas and D'Ari, 1990; Kauffman, 1993; Smolen et al., 2000; Gibson and Mjolsness, 2001; Hasty et al., 2001a; Savageau, 2001; de Jong, 2002; Shmulevich et al., 2002; for reviews of such models and several references). Most models can be categorized as being "qualitative" or "quantitative." The former models emphasize structural information sharing among genes and lack detailed quantitative description of transcriptional regulation. The later models focus on a quantitative description of transcriptional regulation and are often more biologically oriented than qualitative models. The Boolean network (Kauffman, 1993) is a good example of a qualitative model, whereas, transcriptional regulation models based on ordinary differential equations (ODEs) (Chen et al., 1999) are typical examples of quantitative models.
Typically, a qualitative model (like a Boolean network) is a "coarse" approximation of transcriptional regulation. It may provide some insights into the underlying mechanisms of transcriptional regulation, but it may also lead to biologically erroneous conclusions (e.g., see Hatzimanikatis and Lee, 1999). However, qualitative models may be used to predict steady-state behavior of transcriptional regulation. This is a useful property, because cells are often observed at steady state.
Cells may often transition to different states, due to environmental perturbations or genetic instability, which may result in differentiation during development, irreversible adjustments, or disease. Therefore, it is important to design transcriptional regulation models that sufficiently predict transient as well as steady-state behavior. It is believed that ODE-based models can accomplish this goal (e.g., see Hammond, 1993; Elowitz and Leibler, 2000; Gardner et al., 2000; Yildirim and Mackey, 2003). Models based on ODEs are considered to be more detailed than qualitative models, but require structural knowledge of the transcriptional machinery and of several biological parameters (e.g., identification of promoters, regulatory regions, transcription factors, mRNA decay rates, etc.). This knowledge is not currently available for most organisms, and it is thought to be the main disadvantage of ODE-based models. However, several current efforts are geared toward determining the structure of the transcriptional machinery and estimating its parameters (e.g., see Hammond, 1993; Endy et al., 1997; Arkin et al., 1998; Tavazoie et al., 1999; Akutsu et al., 2000b; Gardner et al., 2000; Turner and Varshavsky, 2000; Lee et al., 2002; Ronen et al., 2002; Wang et al., 2002). For these reasons, ODE-based models are becoming increasingly attractive as models for transcriptional regulation.
An attractive feature of a Boolean network is that it dynamically relates the state of transcriptional regulation at time t to its state at time t - 
t, for some t > 0. The state of transcriptional regulation is summarized by binary-valued variables, which are dynamically related from t - t to t by means of Boolean functions. In this formulation, the analysis and simulation of transcriptional regulation employs theoretical and computational tools from discrete dynamical systems theory (e.g., see Sandefur, 1993), specialized to the Boolean case.
On the other hand, ODE-based models represent transcriptional regulation by a (usually large) system of nonlinear ODEs. According to this formulation, the state of the system is summarized by real-valued variables, with regulatory interactions taking the form of differential and nonlinear functional relationships. Due to the size and nonlinear structure of the system, it is not in general possible to develop mathematical techniques for its analysis. In this case, analysis is done by means of numerical techniques and computer simulations. In particular, the system may be solved by a numerical technique, like a Runge-Kutta or a predictor-corrector method (e.g., see Meir et al., 2002). Although these methods lead to general analysis and simulation techniques for transcriptional regulation, they may not be efficient, and direct biological interpretation of the various terms in the resulting equations may not be possible.
As noted in Meir et al. (2002), instead of using general techniques, it may be more preferable to derive a numerical approach to transcriptional regulation by exploiting the specific nature of the problem at hand. In this article, we investigate the possibility of doing so, by replacing an ODE-based model for transcriptional regulation with a nonlinear discrete dynamical system that is "biologically transparent," in the sense that the resulting equations preserve the biological relevance and structure of the original model. This allows us to construct a biologically relevant quantitative model for transcriptional regulation that, like the Boolean network, enjoys attractive dynamical properties and is amenable to efficient simulation and analysis.
The system proposed in this article is directly obtained from a well-known model of transcriptional regulation based on ODEs. The ODE-based model is derived for a large population of cells by applying simple arguments of chemical kinetics on the processes of transcription and translation. It is required that the cell population is large, because the derivation of the ODE-based model relies on the Boltzmann distribution of statistical mechanics, which specifies how energy is distributed in a large population of identical molecules at statistical equilibrium. Because the ODE-based model is central to our work, we show in the next section how this model is derived from first principles. The purpose of our discussion is to clarify the limitations of modeling transcriptional regulation by means of ODEs, and to establish terminology and notation.
In the third section, we show how to model transcriptional regulation by means of a discrete dynamical system. We view the processes of transcription and translation as being discrete, and replace the actual transcriptional machinery with one for which the speeds of transcription and translation, as well as the delays in cis-regulation, are constant and equal to their mean values. We refer to this as an "average" transcriptional machinery. Therefore, the discrete dynamical system derived in this section models an "average" behavior of transcriptional regulation. The system is obtained by discretizing the ODE-based model discussed in the previous section. The discretization step is taken to be the time 
t that it takes the RNA polymerase II to read one nucleotide. Moreover, we assume that, for each t = t, 2t, ..., both the fraction of DNA templates committed to the transcription of a given gene and the mRNA concentration associated with that gene, remain constant within the time interval [t - t, t). The resulting dynamical system is referred to as a discrete transcriptional regulatory system. It is specified by means of parameters that characterize transcription, translation, and degradation, by functionals that characterize cis-regulation, and by time delays.
In the fourth section, we discuss the steady-state behavior of the discrete model under consideration. Our discussion is motivated by the fact that the steady-state behavior of a model for transcriptional regulation may be used to characterize the cell's phenotype, and focuses on three results. The first result shows that, at steady state, the mRNA concentration vector of the discrete model "decouples" from the steady-state protein concentration vector, in the sense that one vector can be derived as a solution of a system of (nonlinear in general) equations without knowledge of the other vector. The second result shows that there is a one-to-one correspondence between the steady-state mRNA concentration vector and the steady-state protein concentration vector. This suggests that, at steady state, mRNA expression data may be used to characterize protein activity, provided that a sufficiently good estimate of the steady-state mRNA concentration vector can be inferred from such data (this also requires that the model parameters associated with translation are known). The final result shows that the discrete model has the same steady-state behavior as the associated ODE-based model. This result, together with several computational advantages underlying the discrete model, indicates that it may be more preferable to use the proposed discrete dynamical system as a model for transcriptional regulation, than the original ODE-based model.
In the next section, and by using classical arguments from chemical kinetics, we specify the nonlinearities underlying cis-regulation and include both activators and repressors as well as the notion of regulatory modules in our formulation. The derivation is based on the assumption that regulatory proteins are free to bind at several distinct sites in a promoter's regulatory region, and on the assumption that different proteins do not interact with each other or affect each other's binding affinity. Moreover, the inclusion of repressor proteins in the formulation focuses on a specific repression mechanism by which, when a repressor protein binds on a DNA template, it either blocks the recruitment of the transcription initiation complex on the promoter or prevents the release of RNA polymerase II. Finally, we show how to model cis-regulation organized in a modular fashion. According to this organization, transcriptional activity of a given gene may be controlled by a set of distinct modules, with each module asserting its own transcriptional control, independently of other modules.
In the sixth section, we discuss several properties of the proposed discrete model, related to homeostatic and epigenetic regulation as well as to Boolean networks, and elaborate on their significance. In particular, the structure of the discrete model under consideration predicts a specific response of transcriptional regulation to changes in the cellular environment, and suggests that mRNA and protein degradation, together with the rates of mRNA and protein synthesis, may play an important role in homeostatic regulation. We show that the functional form of cis-regulation is scale-invariant. This property implies that an increase (decrease) in the rates of translation, accompanied by a proportional decrease (increase) in the affinity constants underlying the binding of proteins on a promoter's regulatory region, does not change the steady-state mRNA concentration but proportionally increases (decreases) the steady-state protein concentration. It also implies that an increase (decrease) in the rates of transcription, accompanied by a proportional decrease (increase) in the affinity constants, proportionally increases (decreases) both the steady-state mRNA and protein concentrations. These properties suggest that the rates of transcription and translation, together with the affinity constants, may play an important role in epigenetic regulation. We also discuss the problem of specifying the underlying parameters, we briefly remark on the appropriateness of the Hill function as a model for cis-regulation, and introduce a parameter that provides a trade-off between model accuracy and computational efficiency. Finally, we provide a mathematical argument that indicates a limitation of using a Boolean network as a model for transcriptional regulation.
In the seventh section, we present simulations, based on transcriptional regulation of a hypothetical metabolic pathway, that illustrate several aspects of the proposed discrete model. By varying the parameters of the model, and observing how these changes affect mRNA and protein activity, we demonstrate that the nonlinear discrete dynamical system proposed in this article may effectively be used to model transcriptional regulation in a biologically relevant way.
Finally, in the last section, we summarize our conclusions.
We believe that the main contribution of this work is to show that, by using available biological information pertaining to the processes of transcription, translation, and cis-regulation, we can derive a nonlinear discrete dynamical system that may serve as a promising and testable model for transcriptional regulation. Our theoretical discussions and simulations indicate that the proposed model is capable of sufficiently predicting basic biological function and producing biologically relevant responses. Finally, the discrete dynamical nature of the proposed model makes it very attractive for large-scale computational analysis and simulation studies of transcriptional regulation.
|
REVIEW OF A CONTINUOUS MODEL
|
|---|
To model transcriptional regulation, we consider a large population 
of genetically identical cells that express the same set of G (distinct) genes, and denote those genes by 
We take the population to be large because the derivation of the continuous model discussed here (as well as the model for cis-regulation discussed later in this article) uses the Boltzmann distribution of statistical mechanics. The Boltzmann distribution specifies how energy is distributed in a large population of identical molecules (DNA templates, mRNAs, and regulatory proteins in our case) at statistical equilibrium. We view transcriptional regulation as a complex system of interacting genes and regulatory proteins (transcription factors), whose state at time t is summarized by the G x 1 vectors r(t) and p(t), given by
where ri(t) and pi(t) are the concentrations in , at time t, of the mRNAs and regulatory proteins produced by the ith gene (measured in mol/L or molarity M; the concentrations considered in this article are with respect to the total cellular volume in ). We consider systems that are "complete," in the sense that p consists of all proteins that regulate transcription of the mRNAs in r. For ease of presentation, we focus on a transcriptional machinery that is "isolated," in the sense that it is not subject to external inputs. If necessary, our formulation can be modified to consider those cases as well (see the example depicted in Fig. 5).
Given a target gene, we need to mathematically describe how its expression level (i.e., the mRNA concentration produced by this gene) is regulated by the expression levels of other genes. Fig. 1 depicts a block diagram of a model for transcriptional regulation, in which a target gene 3 is directly regulated by two other genes, 1 and 2. By "direct regulation" we mean that changes in the expression levels of genes 1 and 2 may produce a change in the expression level of gene 3 with no mediation from other genes. According to this model, the mRNAs transcribed from genes 1 and 2, with respective concentrations r1(t) and r2(t), at time t, are translated into two regulatory proteins whose concentrations are p1(t) and p2(t). These proteins bind to the control region of gene 3 and regulate the recruitment of general transcription factors and RNA polymerase II (for eukaryotic cells) to the gene's promoter. This step is referred to as cis-regulation. After the general transcription factors and RNA polymerase II have been assembled and positioned on the promoter, the RNA polymerase II initiates transcription of gene 3, whose mRNA concentration at time t is r3(t).
View larger version (9K):
[in this window]
[in a new window]
|
FIGURE 1 Block diagram of a model for transcriptional regulation. The target gene 3 is directly regulated by two genes 1 and 2. Transcriptional regulation involves three steps: translation, cis-regulation, and transcription.
|
|
In the diagram depicted in Fig. 1, we have assumed that mRNAs and proteins do not decay, and that the tasks of translation, cis-regulation, and transcription are completed instantaneously. It is a well-known fact however that mRNAs and proteins are subject to degradation and that the time required to complete transcription and translation is not negligible. Transcription is subject to a time delay for completing RNA chain elongation, whereas, translation is subject to a time delay for completing the elongation phase of protein synthesis. Moreover, and for controlling the assembly of the transcription initiation complex (i.e., the general transcription factors and RNA polymerase II) at the promoter, appreciable time is required for the transport of proteins to the nucleus, for the binding of these proteins to the appropriate DNA regulatory sequences, and for recruiting the general transcription factors at the promoter. These effects can be accounted for, by assuming that translation, cis-regulation, and transcription are subject to time delays 
p,i, c,i, and r,i, respectively, for i 
. In general, these delays depend on the particular genes under consideration.
To obtain a model for transcriptional regulation, we need to mathematically describe the three steps of translation, transcription, and cis-regulation. To derive a mathematical model for translation, we adopt the following notation: T, absolute temperature (in degrees Kelvin, K); R, gas constant (1.9872 cal mol-1 K-1); Utr,i, activation energy of translation of the ith mRNA (in cal/mol); Udg,i, activation energy of degradation of the ith regulatory protein (in cal/mol); ri(t | U > Utr,i), concentration, at time t, of ith mRNA molecules in with energy greater than the activation energy Utr,i; pi(t | U > Udg,i), concentration, at time t, of ith regulatory protein molecules in with energy greater than the activation energy Udg,i.
The activation energy depends on the specific aspects of the underlying chemical reaction. By using standard arguments from chemical kinetics (e.g., see Moore and Pearson, 1981; Chapter 5 and Espenson, 1995; Chapter 7), we take the rate of protein synthesis (per second) during translation to be proportional to the concentration of mRNAs with energy >Utr, with proportionality constant 
tr. Similarly, we take the rate of protein degradation (per second) to be proportional to the concentration of proteins with energy >Udg, with proportionality constant dg.
By focusing on the macroscopic behavior of translation during the time interval [t, t + t], for some t > 0, we can write:
 | (1) |
for 
where tr,i and dg,i are the proportionality constants (measured in s-1) associated with the two reactions of mRNA translation and protein degradation, respectively. To obtain Eq. 1, we use a fundamental result of statistical mechanics, which states that, in a large population of identical molecules at statistical equilibrium with concentration 
, the concentration (U) of molecules with kinetic energy U is given by the Boltzmann distribution
This leads to
 | (2) |
where (U > U0) denotes the concentration of molecules in the population with kinetic energy greater than some threshold U0. From Eq. 1, we obtain
 | (3) |
By taking limits, as 
on both sides of Eq. 3, and by setting
 | (4) |
we obtain the following system of rate equations:
 | (5) |
These first-order ODEs imply that the rate of change in the concentration of the ith regulatory protein at time t is proportional to the expression level ri(t - p,i) of gene i at time t - p,i. Moreover, it implies that this protein degrades at a rate 0 < 
i < 
(in s-1), which is proportional to its concentration. In Eq. 5, 0 < 
i < (in s-1) is the rate of translation; i.e., the proteins synthesized per second from a mol of mRNA.
By following similar arguments, we can show that transcription can be modeled by the following system of rate equations:
 | (6) |
where 0 
ci(t) 1 is the fraction, at time t, of DNA templates in that are committed to the transcription of gene i, and 0 < 
i < is the transcription rate of gene i; i.e., the concentration of mRNAs synthesized per second when all DNA templates in are committed to the transcription of gene i (in M s-1). We say that a DNA template is "committed" to the transcription of a gene, if it has successfully recruited the transcription initiation complex and has anchored it at the promoter of that gene. Note that a DNA template that is committed to transcription may not necessarily lead to transcription initiation. For this to happen, the energy of the committed DNA template should be greater than the activation energy of transcription initiation. The first-order ODEs in Eq. 6 imply that the rate of change in mRNA concentration produced from gene i at time t is proportional to the fraction ci(t - r,i) of DNA templates committed to the transcription of gene i at time t - r,i. Moreover, they imply that these molecules degrade at a rate 0 < ßi < (in s-1), which is proportional to their concentration.
In general, the cis-regulation of a target gene i may be modeled by the following equations:
 | (7) |
where 
i[·] is a (nonlinear) function, which is specific to the target gene under consideration, and 
is the set of all genes in that produce proteins, which regulate the transcription of the ith gene. We refer to i[·] as the cis-regulatory function of gene i. Moreover, we refer to as the regulatory set of gene i and to the genes in as the regulating genes of gene i. We call the collection 
of all regulatory sets a transcriptional regulatory network (tRN). In Eq. 7, we assume that transcription is controlled by the protein products, at times t - c,j, 
obtained by translating the regulating genes of the target gene i.
We note here that several variations of the model governed by Eqs. 5–7 have been proposed in the literature (e.g., see Hargrove and Schmidt, 1989; Mjolsness et al., 1991; Mestl et al., 1995; Endy et al., 1997; Wolf and Eeckman, 1998; Chen et al., 1999; Hatzimanikatis and Lee, 1999; Akutsu et al., 2000b; Cherry and Adler, 2000; von Dassow et al., 2000; Elowitz and Leibler, 2000; Gardner et al. 2000; Hasty et al., 2000; Voit, 2000; Smolen et al., 2000; Gibson and Mjolsness, 2001; Mjolsness, 2001, Vohradsk
, 2001; Wahde and Hertz, 2001; de Jong, 2002; Yildirim and Mackey, 2003, and the references therein). A limitation of Eqs. 5–7 is that they only apply to a large population cells. Moreover, these equations are derived by employing a macroscopic view of the chemical reactions underlying translation, cis-regulation, and transcription. The resulting ODE-based model oversimplifies the complex structure of a cell's transcriptional activity, by ignoring several factors affecting such activity. For example, Eq. 5 ignores the effects of mRNA transport from the nucleus to the cytoplasm and mRNA localization in the cytoplasm, whereas Eq. 6 does not take into account the mechanisms of RNA processing. Eq. 7 oversimplifies transcriptional control by ignoring, for example, complex interactions, inside the cis-regulatory mechanisms, among regulatory proteins, general transcription factors, RNA polymerase II, chromatin remodeling complexes and DNA, and by ignoring the role that protein-DNA complexes play in transcriptional regulation. Finally, the model does not consider how protein folding affects transcriptional regulation and ignores several biochemical interactions among proteins and interactions between different biological and signaling pathways. Nevertheless, the ODE-based model governed by Eqs. 5–7 provides a "first-order" approximation of transcriptional activity that leads to a mathematically tractable model for transcriptional regulation.
To conclude this section, note that, by solving Eqs. 5 and 6 with respect to pi(t) and ri(t), we obtain
which in turn result in
 | (8) |
 | (9) |
for some t > 0. According to Eq. 8, the concentration ri(t) of the ith mRNA present in the cytoplasm at time t equals the concentration 
of the mRNA that survives degradation during the time interval [t - t, t), plus the concentration 
of the new mRNA that is synthesized by transcription and survives degradation during the same interval. According to Eq. 9, the concentration pi(t) of the ith protein present in the cytoplasm at time t equals the concentration 
of the protein that survives degradation during the time interval [t - t, t), plus the concentration 
of the new protein that is synthesized by translation and survives degradation during the same interval.
|
A DISCRETE MODEL
|
|---|
The previous ODE-based model provides a continuous description of transcriptional regulation. However, the processes of transcription and translation may be thought as being discrete. During transcription, the RNA polymerase II moves along the DNA in a stepwise fashion and extends the growing RNA chain by adding one nucleotide at a time (see Alberts et al., 2002, pp. 302–304). Similarly, during translation, a ribosome moves along an mRNA transcript by sequentially processing groups of three nucleotides (codons), and extends the growing polypeptide chain by adding one amino acid at a time (see Alberts et al., 2002, pp. 342–344). These observations, together with the need for solving Eqs. 5– 7 using computational techniques, motivates us to derive a discrete model for transcriptional regulation.
The motion of RNA polymerase II along a DNA molecule may not be smooth (see Alberts et al., 2002, p. 313); its speed may depend on time, the particular gene transcribed, and other factors. Moreover, ribosomes may translate with different speeds at individual codons (e.g., see Sørensen and Pedersen, 1991), whereas, for each the cis-regulation delay c,i may fluctuate. To avoid complications, we replace the actual transcriptional machinery with one for which the speeds of transcription and translation, vr and vp, and the cis-regulation delays 
are all constants, taken to be equal to their mean values. We refer to this machinery as an "average" transcriptional machinery. This implies that the discrete dynamical system to be derived in this section will model an "average" transcriptional activity. Because the transcription speed is constant, the transcription delay, r,i, will be an integer multiple of the time t that takes the RNA polymerase II to read one nucleotide. For eukaryotic cells, we may take the average transcription speed vr 
20 nucleotides/s (see Alberts et al., 2002, p. 304), in which case t 0.05 s, whereas, we may take the average translation speed vp 2 codons/s (see Alberts et al., 2002, p. 343). For this value of vp, the translation delay, p,i, is also an integer multiple of t. Finally, we assume that the "average" transcriptional machinery is also characterized by cis-regulation delays c,i, i , which are integer multiples of t as well.
For each we make the following two assumptions: 1) for each t = t, 2t, ..., the fraction ci of DNA fragments committed to the transcription of gene i remains constant in the time interval [t - t, t); and 2) for each t = t, 2t, ..., the mRNA concentration ri remains constant in the time interval [t - t, t).
In view of the small value of t, as compared to the large timescale of transcription (recall that t 0.05 s in eukaryotic cells, as compared to the duration of a typical transcription reaction, which ranges from minutes to hours), we may shift all transcription commitments within the interval (t - t, t), for t = t, 2t, ..., to time t, with negligible effects on transcription. Therefore, we may approximately assume that no new DNA templates commit to transcription within the time interval (t - t, t), which explains assumption 1. On the other hand, and due to the fact that the transcription delay r,i is an integer multiple of t, new mRNAs are synthesized only at integer multiples of t. This, together with the previous observation, implies that new mRNAs are synthesized only at times t = t, 2t, .... Moreover, experimental evidence suggests that mRNA half-lives are much larger than t (e.g., see Wang et al., 2002, and compare the value t 0.05 s for eukaryotic cells with the mRNA half-lives in yeast, which range from 
3 min to >90 min). In view of these observations, and assumption 1, we may conclude that assumption 2 is a reasonable assumption as well.
We can now employ the previous two assumptions to show that, by using t as a basic discretization step and by replacing the actual transcriptional machinery with an "average" one, Eqs. 5– 7 can be transformed into a discrete dynamical system that can effectively simulate transcriptional regulation in an iterative fashion. From assumptions 1 and 2, we have that
and
for i , which, together with Eqs. 8 and 9, with t = t, result in
 | (10) |
 | (11) |
where one iteration corresponds to the time step t, nr,i = r,i/t, np,i = p,i/t, and
 | (12) |
Moreover, from Eq. 7, we have that
 | (13) |
where nc,j = c,j/t.
According to Eq. 10, the concentration ri(n) of the ith mRNA present in the cytoplasm at step n equals the concentration 
of the mRNA that survives degradation from step n - 1 to step n, plus the concentration is(ßi, t)ci(n - nr,i - 1) of the new mRNA that is synthesized by transcription and survives degradation between these two steps. According to Eq. 11, the concentration pi(n) of the ith protein present in the cytoplasm at step n equals the concentration 
of the protein that survives degradation from step n - 1 to step n plus the concentration is(i, t)ri(n - np,i - 1) of the new protein that is synthesized by translation and survives degradation between these two steps.
In the following, and to ease notation, we take the time delays p,i, c,i, and r,i to be independent of i. In this case, Eqs. 10, 11, and 13 can be written in the following compact form:
 | (14) |
 | (15) |
where 
= (r + c)/p. In Eqs. 14 and 15, 
b, c, 
, and 
are G x G diagonal matrices, given by
Moreover, 
b(n) and c(n) are G x G diagonal matrices, given by
where
and
Finally, 
[p(n - np - 1)] is an G x 1 vector-valued functional whose ith element is 
We refer to [·] as the cis-regulatory functional.
The iterations suggested by Eqs. 14 and 15 are depicted in Fig. 2, when p = t, c = 3t, and r = 2t. These iterations are initialized with an mRNA concentration vector r(0) and a protein concentration vector p(0). This implies that the fate of gene expression is determined by the initial concentrations of mRNAs and proteins. Matrices b and c model mRNA and protein degradation, respectively, whereas, matrices and model the rate of transcription and translation, respectively. For 0 n p/t = np only degradation is present. Translation of mRNAs to proteins takes place for n 
p/t + 1 = np + 1, whereas, transcription takes place for n (c + r)/t + 1 = np + 1. Note that the flow graph depicted in Fig. 2 has a modular structure; it consists of individual stages, with the nth stage being the (nonlinear in general) multi-input/multi-output system depicted in Fig. 3 (for n np + 1), where d denotes delay, such that d-mx(n) = x(n - m).
The model suggested by Eqs. 10–13 requires knowledge of the cis-regulatory functionals 
the degradation parameters 
the transcription and translation rates 
and the delays 
where i = (r,i + c,i)/p,i. We refer to the collection 
as a (discrete) transcriptional regulatory system (tRS).
|
STEADY-STATE BEHAVIOR
|
|---|
An important issue associated with a tRS is whether or not the iterations suggested by Eqs. 10–13 converge to a steady state and, if they do, to characterize that state. In most cases, and in the absence of external control, for a tRS to be biologically plausible, it is required that the mRNA and protein concentration vectors r(n) and p(n) converge, as 
to a steady-state mRNA concentration vector 
and a steady-state protein concentration vector 
In this case,
 | (16) |
from which we have that
 | (17) |
 | (18) |
where
This shows that is a fixed-point attractor of the functional 
r[·], given by
 | (19) |
whereas, 
is a fixed-point attractor of the functional p[·], given by
 | (20) |
We refer to r[·] as the "genomic regulatory functional" and to p[·] as the "proteomic regulatory functional" because the first functional can be used to determine the steady-state mRNA concentration vector, whereas, the second functional can be used to determine the steady-state protein vector.
The fixed-point attractors and may be used to characterize the cell's phenotype. This is based on the assumption that cells may be differentiated by the concentrations of regulatory proteins synthesized at steady state (or, equivalently, by the concentrations of the corresponding mRNAs), which give each cell type its unique characteristics; e.g., see Kauffman (1993) (Chapter 12) and Alberts et al. (2002) (pp. 375–376). It is believed that the transcriptional regulatory machinery of a given organism is hardwired in its DNA. This implies that regulation of transcription is controlled by the same mechanisms, irrespective of cell type. We may however view cell differentiation as being achieved by transcriptional regulation, which guides the tRS to reach steady-state mRNA and protein concentration values and that uniquely characterize the cell type. In this case, the driving force of cell differentiation is said to be "epigenetic" regulation.
An implication of Eqs. 17 and 18 is that, at steady state, the mRNA concentration vector "decouples" from the protein concentration vector 
in the sense that can be obtained as a solution of the system of (nonlinear in general) Eq. 17, without knowledge of whereas, can be obtained as a solution of the system of (nonlinear in general) Eq. 18, without knowledge of 
It is however important to keep in mind that, despite this "decoupling," computation of the steady-state mRNA concentration vector (and the steady-state protein concentration vector ) requires knowledge of the transcription parameters 
, , and the translation parameters 
, , because Eq. 17 (and Eq. 18) depends on those parameters. Note also that there is a one-to-one correspondence between the fixed-point attractors of the genomic and proteomic regulatory functionals, because (recall Eq. 16)
 | (21) |
The second equation above implies that, at steady state, may be determined from 
provided that the underlying translation parameters , are known. This observation suggests that mRNA expression data, obtained by means of microarray gene expression profiling, may be used to characterize protein activity at steady state, provided that the translation parameters , are known, and a sufficiently good estimate of the steady state mRNA concentration vector can be inferred from such data.
Because Eqs. 10–13 have been obtained by discretizing Eqs. 5–7, it is of interest to investigate how the steady-state behavior of the discrete tRS is related to the steady-state behavior of the continuous tRS. From Eqs. 10–12, we can show that
 | (22) |
On the other hand, if 
and 
are the steady states of Eqs. 5 and 6, then (by setting the derivatives in Eqs. 5 and 6 equal to zero), we obtain
 | (23) |
Eqs. 22 and 23 verify that the discrete tRS has the same steady states as the continuous tRS, and show that the steady-state behavior of the discrete tRS is identical to that of the continuous tRS.
Besides fixed-point attractors, the tRS may be subject to limit-cycle attractors, which lead to oscillatory behavior. A tRS with limit-cycle attractors may be useful for modeling periodic cellular behavior, such as cell cycle control or circadian rhythms; see Kauffman (1993) (Chapter 12), and Elowitz and Leibler (2000); Smolen et al. (2000); Goldbeter et al. (2001); Hasty et al. (2001a,b); Tyson et al. (2001). In this article, we do not consider limit-cycle attractors (however, see Fig. 11 d).
|
A MODEL FOR cis-REGULATION
|
|---|
The cis-regulatory functions i[·] in Eq. 7 are at the core of a tRS, because these functions specify how proteins regulate transcription. In this section, we derive a form for these functions by using simple arguments from chemical kinetics (see also Hill, 1985; Wang et al., 1999). Keep in mind however that the resulting model oversimplifies cis-regulation, because cis-regulation is controlled by rather complicated biochemical interactions (e.g., see Holstege et al., 1998).
To model cis-regulation, we consider again a large population of cells, and assume at the moment that the promoter of a given target gene is controlled by two regulatory proteins P1 and P2, with concentrations p1 and p2, respectively. Moreover, we assume that protein P1 is free to bind at anyone of S1 distinct sites of the promoter's regulatory region, whereas, protein P2 is free to bind at anyone of S2 distinct sites, with the binding sites of P1 being different than that of P2. Let D[s1, s2] be a DNA template with s1 out of the S1 sites being occupied by P1 and s2 out of the S2 sites being occupied by P2. The binding of proteins P1 and P2 at the promoter's regulatory region can be described by means of the following reversible reactions:
 | (24) |
for 
If we assume that P1 and P2 do not interact with each other or affect each other's binding activity, then Eq. 24 can be sequentially written as
 | (25) |
 | (26) |
In the following, d[s1, s2] denotes the concentration of DNA templates in with s1 out of the S1 sites being occupied by P1 and s2 out of the S2 sites being occupied by P2. Moreover, Ubd,i denotes the activation energy (in cal/mol) for a regulatory protein Pi to bind on a DNA template, whereas, Uds,i denotes the activation energy (in cal/mol) for Pi to dissociate itself from the template.
At equilibrium, the concentration of free regulatory proteins P1 that bind (per second) on DNA templates D[s1, s2] to produce DNA templates D[s1 + 1, s2] by means of the forward reaction in Eq. 25 must equal the concentration of regulatory proteins P1 freed (per second) by the backward reaction. By using molecular collision theory (e.g., see Moore and Pearson, 1981, Chapter 4), it can be shown that the first concentration is proportional to the concentration of those proteins P1 with kinetic energy >Ubd,1 times the concentration of sites available for P1 to bind to, with proportionality constant bd,1 (measured in M-1s-1). Because each DNA template D[s1, s2] has S1 - s1 sites available for P1 to bind to, the concentration of available binding sites for P1 is (S1 - s1)d[s1, s2]. In this case, the concentration of free regulatory proteins P1 that bind (per second) on DNA templates D[s1, s2] to produce DNA templates D[s1 + 1, s2] by means of the forward reaction in Eq. 25, is given by
where we have used Eq. 2.
On the other hand, the concentration of regulatory proteins P1 freed (per second) by the backward reaction in Eq. 25 is proportional to the concentration of bound P1 molecules on the DNA template D[s1 + 1, s2] with kinetic energy >Uds,1, with proportionality constant ds,1 (measured in s-1). Because each DNA template D[s1 + 1, s2] contains s1 + 1 bound P1 molecules, this concentration is given by
where we have used again Eq. 2. Therefore, at equilibrium, we have that
 | (27) |
 | (28) |
with 1 = bd,1/ds,1 and U1 = Ubd,1 - Uds,1 being the binding free energy. The parameter 
1 (measured in M-1) is characteristic to the binding sites and is referred to as "affinity constant." At equilibrium, and when Ubd,1 = Uds,1, the values of p1 (S1 - s1) d[s1, s2] and (s1 + 1) d[s1 + 1, s2] must be equal; therefore, 1 = 1 M-1. In addition, because Ubd,1 Uds,1, we have that 1 1 .
A similar argument applies to Eq. 26 and leads to
 | (29) |
with U2 = Ubd,2 - Uds,2 and 2 = 1 M-1.
From Eqs. 27 and 29, it can be shown that
 | (30) |
are the so-called Binomial coefficients. If we ignore additional processes underlying cis-regulation (see Wang et al., 1999 for such processes) and assume that, for transcription to be initiated, it is necessary (but not sufficient) that a DNA template is bound by at least one P1 protein or one P2 protein, then the fraction c of DNA templates in committed to the transcription of the target gene will be given by
 | (31) |
The previous assumption agrees with the fact that transcription in eukaryotic cells can be initiated only in the presence of activator proteins (e.g., see Alberts et al., 2002, pp. 312–313). Eq. 31, together with Eq. 30, leads to
where
 | (32) |
Our discussion so far has been based on the assumption that regulatory proteins activate transcription; i.e., binding of regulatory proteins on a DNA template recruits the transcription initiation complex and initiates transcription. However, cis-regulation may also be controlled by regulatory proteins that repress transcription. Although eukaryotic genes employ several mechanisms for repressing transcription, we only consider here the mechanism by which a repressor protein binds on the DNA template, and either blocks the recruitment of the transcription initiation complex to the promoter, or prevents the release of the RNA polymerase II (see Alberts et al., 2002, pp. 405–406). This implies that, once a repressor protein binds on a DNA template, transcription cannot be initiated by that template, and leads to a simple model for the repression of transcription. Keep in mind however that, if necessary, it may be possible to derive models for other repression mechanisms as well.
If we assume that the previous protein P2 is a repressor, then the repression mechanism under consideration implies that transcription may be initiated only if a DNA template is free of protein P2 and there is at least one activator protein P1 bound to it. Then, the fraction c of DNA templates committed to the transcription of the target gene will be given by
which leads to
as opposed to Eq. 31. In general, if the promoter of a given target gene is controlled by activators P1, P2,..., Pk and repressors Pk+1, Pk+2,..., PJ, it can be shown that
 | (33) |
It is believed that the organization of cis-regulation is modular; see Davidson (2001) (Chapter 1), Alberts et al. (2002) (pp. 408–413), Arnone and Davidson (1997), and Bolouri and Davidson (2002). This means that the regulatory region of a target gene may be partitioned into several entities (modules), with each entity being associated with different sets of regulatory proteins, which may assert a different type of control on transcription. This modular structure allows a gene to express itself under different conditions and different contexts.
We now derive a model for modular cis-regulation. For the purpose of
TOP
ABSTRACT
INTRODUCTION
denotes an activator, whereas, 
denotes a repressor.
