German Cavelier | Freelancer Systems Biology

Systems Biology

PROTEINS: Structure, Function, and Bioinformatics 55:339 –350 (2004) Data-Based Model and Parameter Evaluation in Dynamic Transcriptional Regulatory Networks German Cavelier1 and Dimitris Anastassiou1,2* 1 Genomic Information Systems Laboratory, Department of Electrical Engineering, Columbia University, New York, New York 10027 2 Columbia Genome Center and Center for Computational Biology and Bioinformatics, Columbia University, New York, New York 10027 ABSTRACT Finding the causality and strength of connectivity in transcriptional regulatory networks from time-series data will provide a powerful tool for the analysis of cellular states. Presented here is the design of tools for the evaluation of the network’s model structure and parameters. The most effective tools are found to be based on evolution strategies. We evaluate models of increasing complexity, from lumped, algebraic phenomenological models to Hill functions and thermodynamically derived functions. These last functions provide the free energies of binding of transcription factors to their operators, as well as cooperativity energies. Optimization results based on published experimental data from a synthetic network in Escherichia coli are presented. The free energies of binding and cooperativity found by our tools are in the same physiological ranges as those experimentally derived in the bacteriophage lambda system. We also use time-series data from high-density oligonucleotide microarrays of yeast meiotic expression patterns. The algorithm appropriately ﬁnds the parameters of pairs of regulated regulatory yeast genes, showing that for related genes an overall reasonable computation effort is sufﬁcient to ﬁnd the strength and causality of the connectivity of large numbers of them. Proteins 2004;55:339 –350. © 2004 Wiley-Liss, Inc. Key words: binding free energy; evolutionary algorithms; high-density oligonucleotide microarrays; regulatory networks; transcription; transcription factor INTRODUCTION Transcriptional activation of genes in living cells must follow an orderly pattern of steps in order to appropriately regulate cellular functions. Together with experimental data, a systems perspective would be useful in describing such a complex pattern.1– 6 The increasing availability of complete genome data from a growing number of organisms necessitates computational algorithms and tools1,6 – 8 with which to decipher the cell’s schema for activating and deactivating genes under various conditions.1,6,9 A modelindependent instrument that is able to obtain such characteristics would be of great value in the assessment of different competing models that are compatible with a © 2004 WILEY-LISS, INC. given data set.6 This article describes the design, construction, and testing of such an instrument. Interactions between genes in the lambda bacteriophage lysis/lysogeny decision circuit have been studied with the theory of neural networks.10 However, a biophysics- and thermodynamics-based model11,12 yielded more insight into the transcription factor/promoter interactions of this system than an abstract nonlinear algebraic function. Among eukaryotes, the yeast Saccharomyces cerevisiae is the most extensively studied in terms of its transcriptional regulation. Several efforts have been directed toward the elucidation of connections and causality in the network of transcriptional regulators in this organism.9,13–15 Genomewide analyses have been performed on it,16,17 yielding not only the network structure but also new functional classiﬁcations and a description of the motifs present in the structure.18 The interactions between transcriptional regulatory genes and their causal relationships can be obtained by other methods from available data.9,19 In this initial effort, we limit ourselves to deterministic, continuoustime, ordinary differential equation models of genetic networks.3,6,10,20 –23 The modeling issues themselves are addressed only tangentially, because they have been well described elsewhere.6,11,12 MATERIALS AND METHODS Modeling The mathematical model chosen here to represent the dynamical system of a network of interconnected genes and gene products is based on a state-variables description of the system with nonlinear ordinary differential equations3,20: drᠬ ⫽f共rᠬ 兲⫺Vrᠬ , dt (1) where n is the total number of genes in the network, P is the total number of gene products (proteins) in the netAbbreviations: ES, evolution strategy; GFP, green ﬂuorescent protein; RNAP, RNA polymerase *Correspondence to: D. Anastassiou. E-mail: anastas@ee. columbia.edu. Received 15 July 2003; Revised 28 September 2003; Accepted 20 October 2003 Published online 5 March 2004 in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/prot.20056 340 G. CAVELIER AND D. ANASTASSIOU work, rᠬ is the mRNA and protein concentrations [(n ⫹ P)-dimensional vector, function of time t], f(rᠬ ) is the transcription/translation function, and V is the degradation rate of mRNA/protein [nondegenerate, constant, (n ⫹ P) ⫻ (n ⫹ P), diagonal matrix]. The function f(rᠬ ) may be a lumped6,12 representation of the mRNA/protein variables in the system or may be simply a phenomenological6 description, such as a saturable function of the sum of positively and negatively weighted molecular concentrations. We treat the function f(rᠬ ) for transcription as a nonlinear function, whereas translation is described with a linear law. We used three increasingly complex representations for the nonlinear function. The ﬁrst one tested was a continuous, sigmoidtype function of the sum of weighted components of rᠬ , similar to the activation function used in training artiﬁcial neural networks.22 This corresponds to a phenomenological6 function. The nonlinear function in this case is therefore 1 2 再冋冉 tanh ⫺S 䡠冊册冎 P ⫺1 THR ⫹1 , (2) where P is the protein (unbound repressor dimer) concentration; THR is a threshold for P, where half-maximal repression occurs; and S is a slope factor (similar to the Hill coefﬁcient below). A more biophysically based approach would be the Hill function. The Hill function, as used in deterministically modeling the synthetic oscillatory network of transcriptional regulators,3 was therefore the second representation used here: dm i mi ⫹␣0 ⫹ ␣ 䡠 f共Pj兲, ⫽⫺ dt ␶m used and analyzed,24 and here we restate it using our notation. The ﬁrst step is to write the biochemical reaction for the complex formation between the transcription factor (e.g., an activator with concentration [A]) and its DNA binding site. This binding reaction can be written as ka A⫹DNA | 0 A 䡠 DNA, d关A 䡠 DNA兴 ⫽k a 䡠关Afree兴䡠关DNAfree兴 ⫺ kd关A 䡠 DNA兴, dt (7) where [A 䡠 DNA] is the (average) concentration of the complex between A and DNA, ka is the association rate of A to DNA, kd is its dissociation rate, [Afree] is the concentration of unbound activator in its molecular form for binding (sometimes these transcription factors dimerize or oligomerize before binding, but we do not include such detail), and [DNAfree] is the (average) concentration of unbound DNA. The equilibrium dissociation constant for the reaction in Eqs. (6) and (7) is deﬁned at a stationary steady state. The use of the stationary steady state is justiﬁed by noting that these biochemical reactions are much faster than any of the other mechanisms involved in Eqs. (3)–(5) and they can be considered comparatively instantaneous. Therefore, at stationary steady state (ss) we have d关A 䡠 DNA兴 ⬅ 0, dt (3) 关K diss兴eq ⬅ Kdiss ⬅ dP j Pj ⫽ ⫺ ⫹ kt 䡠 mi, dt ␶P (4) tetR, cI 冉 i⫽lacI, j⫽CI, LACI, TETR 冊, (5) (6) kd 册 kd ka ⫽ ss 关Afree兴䡠关DNAfree兴 . 关A 䡠 DNA兴 (8) (9) The fractional saturation or fractional (average) occupancy of DNA by A is deﬁned by where mi is the mRNA concentration, Pj is the corresponding (translated) protein concentration, ␣0 is the “leakiness” of the promoter (mRNA concentration per unit time produced from a given promoter type during continuous growth in the presence of saturating amounts of repressor), ␣ ⫹ ␣0 is the mRNA concentration per unit time produced from a given promoter type during continuous growth in the absence of a repressor, kt is the protein translation coefﬁcient, ␶P is the protein degradation time constant, and ␶m is the mRNA degradation time constant. Note that in Eqs. (3)–(14) the times are in minutes and the concentrations are in nanomoles. In the following paragraphs we formulate a biochemical framework for deriving the nonlinear function for transcription f(P). A biochemical approach to modeling the system is comprehensive because from its results one can derive, by including more or less detail as desired, several different models as discussed below. Each of them is appropriate for a different case according to the degree of complexity. This biochemical approach to modeling has been previously fa ⬅ 关A 䡠 DNA兴关A 䡠 DNA兴 ⬅ . 关DNAtot兴关A 䡠 DNA兴 ⫹ 关DNAfree兴 (10) Obtaining [DNAfree] from Eq. (9) and replacing it in Eq. (10) we get the fractional saturation or fractional (average) occupancy of DNA by the activator A: fa ⫽ 关Afree兴 . Kdiss ⫹ 关Afree兴 (11) For large amounts of A, the concentration of unbound A is practically equal to the total concentration of A: [A] ⬵ [Afree]. Therefore, fa ⫽ 关A兴 ␳a , ⫽ Kdiss ⫹ 关A兴 1 ⫹ ␳a (12) where ␳a ⫽ [A]/Kdiss is a normalized activator concentration. The form in Eq. (12) is the nonlinear form to which we assume that the transcription rate is proportional. A similar form would be valid for the binding of the repressor [R], but the transcription rate would then be proportional 341 GENOMIC TRANSCRIPTION BIOPHYSICS TABLE I. States and Energies for System With Two Operators (Op1, Op2), One Promoter (Pr), One Transcription Factor (x), and [RNAP] ⴝ 30 nM ⴝ x2 Site a (Op1) — x — — x x — x Site b (Op2) Site c (Pr) State s energy ⌬Gs Partition function term — — x — x — x x — — — x2 — x2 x2 x2 ⌬G1 ⫽ 0 ⌬G2 ⌬G3 ⌬G4 ⌬G5 ⫽ (⌬G2 ⫹ ⌬G3 ⫹ ⌬G23) ⌬G6 ⫽ (⌬G3 ⫹ ⌬G4 ⫹ ⌬G34) ⌬G7 ⫽ (⌬G2 ⫹ ⌬G4 ⫹ ⌬G24) ⌬G8 ⫽ (⌬G2 ⫹ ⌬G3 ⫹ ⌬G4 ⫹ ⌬G234) 1 (k2) 䡠 x (k3) 䡠 x (k4) 䡠 x2 (k2k3k23) 䡠 x2 (k3k4k34) 䡠 x2x (k2k4k24) 䡠 x2x (k2k3k4k234) 䡠 x2x2 to the fraction of unoccupied DNA or (1 ⫺ fr), because a bound repressor represses transcription: 关R兴 Kdiss 1 , 共1 ⫺ f r兲 ⫽ 1⫺ ⫽ ⫽ Kdiss ⫹ 关R兴 Kdiss ⫹ 关R兴 1 ⫹ ␳r (13) where ␳r ⫽ [R]/Kdiss is the normalized repressor concentration. The same expression can be found by applying statistical thermodynamical analysis.25 It is known that when measurements of these binding reactions are performed, the exact sigmoid form represented by Eq. (12) or Eq. (13) is sometimes not closely followed by the experimental data. For example, the slope may be steeper.25 This has been explained by calling into play additional details in the mechanism, for example, dimerization or oligomerization of the transcription factors before binding, cooperativity upon binding, or multiple binding sites.24,25 When these additional considerations are taken into the modeling of the nonlinear function11,12,24 –27 by writing and solving a much more detailed and complete set of biochemical reactions, it is usually found that the end effect is an increase in the switching slope of the nonlinear sigmoid in Eqs. (12) and (13). This has been modeled phenomenologically by introducing the Hill, or similar, exponential coefﬁcients,24 as shown below. Because the results of the parameter identiﬁcation that we performed with experimental data generally produced a Hill coefﬁcient much larger than one, we did not use the forms in Eqs. (12) and (13) but instead the Hill function shown in Eq. (14). The normalized Hill nonlinear (repression) function is therefore introduced in Eq. (3) as follows: f共P兲 ⫽ 1⫹ 1 P THR 冉冊 nHill , (14) where nHILL is the Hill coefﬁcient, THR is the repressor threshold concentration necessary to half-maximally repress a promoter [equivalent to Kdiss in Eq. (9)], and P is the concentration of unbound repressor dimers. The form of f(P) in Eq. (14) corresponds to the case where P is a repressor. In the case that P is an activator, replace f(P) by [1 ⫺ f(P)]. The effects of multiple operator binding sites and of cooperativity may alternatively be included using the thermodynamically derived model explained below. The case of multiple binding sites is relevant to the gene regulatory organization that is usually present in (mostly eukaryotic) genes, and therefore we will discuss it in some detail. The Hill function is a simpliﬁcation of the more realistic thermodynamic and biophysically based description of the process of transcription, as described elsewhere.6,11,12,28 –32 Therefore, a possible improvement would be a thermodynamically derived model for such a process, and this served as the third and most complex representation of the nonlinear function used in this work. As a ﬁrst approximation, we use a simpliﬁed form of the thermodynamically derived function, the so-called “Omodel.”33 In this model, o operators are assumed for repressor binding. In an experimental implementation of the system o ⫽ 2. Under simplifying assumptions,33 a nonlinear repression curve that will replace f(P) in Eq. (3) is obtained, and it constitutes the O-model nonlinear (repression) function: 冉 1 冊 Pd 1⫹ THR o , (15) where Pd is the concentration of unbound repressor dimers. A second approximation would be to use the thermodynamically derived function without simpliﬁcations.11,12,32 The equilibrium statistical thermodynamical probability is used to determine the occupancy of operator sites. The probability of an operator in conﬁguration s is12 具c s 典 ⫽ exp共⫺⌬Gs/RT兲 ⫻ 关R2兴i关RNAP兴j , ¥s,i,j exp共⫺⌬Gs/RT兲 ⫻ 关R2兴i关RNAP兴j (16) where ⌬Gs is the Gibbs free energy of the s conﬁguration, R is the gas constant, T is the absolute temperature, [R2] is the concentration of unbound repressor dimers, and [RNAP] is the concentration of RNA polymerase molecules. Indices i, j indicate the number of each protein bound to an operator in the s conﬁguration. The summation is taken over the total s states whereas i, j have the values 0, 1, or 2 corresponding to the state s. This particular description corresponds to the case of a twooperator site, one-promoter, one-transcription factor system modeling that we chose for the purposes of illustration in the synthetic3 system. Other possible conﬁgurations have been analyzed.32 The ⌬Gs is a sum of all of the free energies of binding to the operator sites, plus cooperative 342 G. CAVELIER AND D. ANASTASSIOU interaction terms (see Table I). The fractional saturation of the promoter by RNAP is the sum of the fractional saturations 具cs典 of those states with RNAP bound to the promoter. Because of the presence of several binding sites, the binding reactions leading at equilibrium to Eq. (16) may be multiple and simultaneous, which leads to a combinatorial explosion in the number of states. An enumeration of the different physiologically relevant states is therefore the ﬁrst step in this modeling (Table I). The stationary steady state at equilibrium for these binding reactions is described in a more expeditious form by the equilibrium statistical thermodynamical treatment described above25 in Eq. (16) than by solving a series of simultaneous biochemical and algebraic binding reactions. As in the case of Eqs. (8)–(13), the underlying assumption is that the concentration of the transcription factors is high enough that the concentration of unbound transcription factor is practically equal to the total concentration of that transcription factor. As mentioned above, we analyzed the twooperator, one-promoter, one-transcription factor model described in detail by Wolf and Eeckman32 but without considering transcription factor dimerization in order to keep the analysis simpler. This system is similar to the one studied by Ackers et al.11 and is a subsystem of that studied by Shea and Ackers.12 It considers only one transcription factor (repressor x ⫽ [R2]) and RNAP binding to the two-operator, one-promoter set of sites in DNA. The eight distinct microscopic conﬁgurations or states possible for this system are shown in Table I, where 冉 k s ⫽ exp ⫺ 冊 ⌬Gs RT 共s ⫽ 1, 2, . . . , Nstates ⫽ 8兲, (17) Parameter Identiﬁcation This problem has been solved in engineering and in biological systems by means of gradient descent, simulated annealing,22 or genetic/evolutionary programming algorithms. The results we obtained using the gradient descent techniques were not satisfactory unless the initial point in the search was very near the optimum. Therefore, the more powerful and global genetic programming approach was taken instead. Because the parameters of the system are real numbers with continuous values, binary coding is not appropriate and the best approach is a type of evolutionary algorithm called an evolution strategy (ES). We developed, calibrated, and tested a self-adaptive34 parameter-identiﬁcation algorithm based on the theory of ES as described by Back and colleagues.35 There were previous genetic regulatory network simulations in this area, although the modeling and evolutionary algorithm used36 were different and no experimental data were used for testing the method. The basic evolutionary algorithm is simple and closely follows the description given by Back et al.35 An optimization problem requires that a set yᠬ 僆 M of free parameters (or object variables) be found for the system under consideration. (For simplicity, we explain the components of the algorithm with a vector yᠬ , although in the systems evaluated below the parameters form a matrix.) A quality or ﬁtness criterion f : M 3 ℜ, which is also called an objective function, must be maximized (or minimized) at the same time: f共yᠬ 兲 3 max. The solution to the global optimization problem is given by a vector yᠬ * such that @yᠬ 僆M : f共yᠬ 兲 ⱕ f共yᠬ *兲. and where the interaction (cooperativity) energies between bound transcription factors are represented by 冉冊 ⌬Gkl k kl ⫽ exp ⫺ . RT (18) To determine the rate of transcription, two additional assumptions were made, namely, that transcription initiation is the rate-limiting step12 and it is proportional to the fractional saturation of the corresponding promoter by RNAP. We also assumed12,32 that the rate of transcriptional initiation is proportional to a single value (the same for each transcription factor) of the rate constant that governs the isomerization of RNAP from a closed to an open complex. The nonlinear function is then in this case k 4 䡠关RNAP兴 ⫹ 共k6 ⫹ k7兲䡠关RNAP兴䡠 x ⫹ k8 䡠关RNAP兴䡠 x2 , 共k1 ⫹ k4 䡠关RNAP兴兲 ⫹ 共k2 ⫹ k3 ⫹ 2 共k6 ⫹ k7兲䡠关RNAP兴兲䡠 x ⫹ 共k5 ⫹ k8 䡠关RNAP兴兲䡠 x (19) where ks ⫽ exp(⫺⌬Gs/RT) from Table I, the concentration of unbound repressor dimers is x ⫽ [R2], and [RNAP] ⫽ 30 nM ⫽ x2 is the concentration of RNAP. (20) (21) ES algorithm In this method, evolution is deﬁned as the result of interactions between the creation of new genetic information and its evaluation and selection. The general structure of an evolutionary algorithm35 is as follows: Evolutionary Algorithm 1: t ⫽ 0 initialize P(t); evaluate P(t); while not terminate do P⬘(t) ⫽ variation [P(t)] ; evaluate [P⬘(t)] ; P(t ⫹ 1) ⫽ select from P⬘(t) ; t ⫽ t ⫹ 1 ; end where P(t) is a population of ␮ individuals at generation t. An offspring population P⬘(t) of size ␭ is generated from the population P(t) using operators such as recombination and mutation. These offsprings are then evaluated by calculating the objective function values f(yᠬ k) for each of the solutions yᠬ k (k ⫽ 1, 2, . . . , ␭) and a selection is performed based on some ﬁtness criterion to drive the process toward better solutions. 343 GENOMIC TRANSCRIPTION BIOPHYSICS The main constraints of the problem at hand are that the parameters cannot lead to negative values of the problem variables (molecular concentrations) and that some parameters cannot be negative. The data include the nonnegativity constraint on the problem variables, and negative parameters usually give exponentially growing variables (which lead to a high error). Therefore, the two constraints can be handled only by performing the optimization procedure, as long as this procedure approaches the data with a small error. Tuning some parameters of the algorithm (such as ␶⬘ and ␶ and others, see below) is also necessary, but this is usually accomplished empirically. Few rules for the design and parameterization of evolutionary algorithms exist, because they are part of a relatively new area of research.35 We used a particular implementation of the ES called (␮, ␭)-ES, where ␮ parents create ␭ ⬎ ␮ offsprings by recombination and mutation; subsequently the best ␮ offsprings are deterministically selected to replace the parents, even if the previous minimum has not been lowered.35,37 This permits us to escape local minima. It has been found by theoretical studies, as well as in practice, that a ratio of ␭/␮ ⫽ 7 is the best for this type of ES, and in most cases the population should be ␭ ⬇ 100. Recombination This operator must be applied as the ﬁrst one in the variation part of evolutionary algorithm 1. Recombination generates a new intermediate population of ␭ individuals [each individual yᠬ k (k ⫽ 1, 2, . . . , ␭) is a complete set of parameters of the model] by ␭-fold application to the parent population, creating one individual per application from ␮ individuals. The recombination types for object variables and for the so-called strategy parameters in self-adaptation (see below) differ; typical examples are discrete recombination (random choices of single variables from parents) and intermediary recombination (often arithmetic averaging, or other variants35). After testing several schemes, we found the best-performing option for the object variables in our system was discrete recombination (i.e., essentially tossing a coin and choosing a randomly selected parameter among the ␮ parent individuals). For the strategy parameters in our system, the best choice after several trials with different formulas was found to be intermediary recombination, as depicted by Back.37 These choices in our particular implementation can be described in pseudocode as follows: discrete recombination for object variables (problem parameters) from 1 to ␭ and for each parameter: toss ⫽ random number between 0 and 1 if toss ⬎ 0.5 recombined ⫽ randomly chosen from ␮ individuals; else recombined ⫽ randomly chosen from ␮ individuals (a different random choice than the previous) ; intermediary recombination for strategy parameters from 1 to ␭ and for each parameter: p1 ⫽ randomly chosen parameter from ␮ individuals; p2 ⫽ randomly chosen parameter from ␮ individuals; (a different random choice than the previous) ; recombined parameter ⫽ p1 ⫹ (p2 ⫺ p1)/2. Mutation The targets of the optimization routine are the object variables (or components of each individual yᠬ k) yi 僆 ℜ (1 ⱕ i ⱕ no). In our realization, mutation adds a uniformly distributed random value between zero and ␴ to each yi. The notation ␴Ui( 䡠 , 䡠 ) indicates such a uniformly distributed random value and that the random variable is sampled anew for each i. y⬘i ⫽ y i ⫹ ␴U i 共0, 1兲. (22) The step size ␴ should ideally be adaptable to the particular conﬁguration of the search surface for each iteration, leading to the so-called self-adaptive strategies.34 Self-adaptation This is a mechanism that controls the step size ␴ by incorporating it as a new vector parameter to be optimized.34 The ES is then applied to an enlarged set of parameters including the object variables yi, in addition to the step sizes ␴i, which are also called strategy parameters. This results in the self-adaptive mutation shown below. The notation N(0, 1) indicates a normally distributed random value from a distribution with mean 0 and standard deviation 1. The notation Ni(0, 1) indicates that the normally distributed random variable is sampled anew for each i: ␴⬘i ⫽ ␴ i exp关␶⬘ 䡠 N共0, 1兲 ⫹ ␶䡠Ni共0, 1兲兴 (23) y⬘i ⫽ y i ⫹ ␴⬘i U i 共0, 1兲 (24) where ␶⬘ ⬀ 共冑 2n o 兲 ⫺1 and ␶ ⬀ 共冑2 冑no兲⫺1 are the recommended formulas for these quantities. However, there is ﬂexibility in choosing the proportionality constant and future research may yield more appropriate values.35 Practice appears to indicate that values for the self-adaptive mutation should be chosen empirically for each particular case. In our case, after testing several mutation schemes, we found the best-performing algorithm for our systems to be the following, written here for reference in pseudocode: self-adaptive mutation for strategy parameters and object variables For each individual yᠬ k(k ⫽ 1, 2, . . . , ␭) calculate: exp1 ⫽ exp(␶⬘ 䡠 N(0, 1)) from 1 to ␭ and for each parameter: calculate: exp2 ⫽ exp(␶ 䡠 Ni(0, 1)) mutated strategy parameter ⫽ (recombined strategy parameter) 䡠 exp1 䡠 exp2 mutated object variable ⫽ recombined object variable ⫹ 344 G. CAVELIER AND D. ANASTASSIOU (mutated strategy parameter)ⴱ (normally distributed random number, 0 ⱕ r ⱕ 1). Deleting nonappropriate individuals In addition to the steps of an ES algorithm as described by Back,37 we added a further discriminating function. This step is usually not included in an evolutionary algorithm. Because of the nature of genetic regulatory networks, an on/off or high/low sequence in time is characteristically found in the time evolution of molecular concentrations. This requires deleting, or penalizing, those responses that are “ﬂat” in time. Therefore, following the recombination/mutation operations of the ES, we subtract the mean value from the simulated response and test the number of zero crossings in the resulting time series. If the number of zero crossings is lower than say, 3 zero crossings, we have a nonappropriate individual (nonoscillating) and its parameters are deleted from the pool of recombined/ mutated parameters. New recombination/mutation operations are then performed until an acceptable individual is found for this criterion. The result of this operation is an appreciable speeding of the optimization process because spurious results are continuously eliminated before they can enter into the pool of individuals to be tested Selection The selection operator is based on a ﬁtness criterion applied to the individuals. The best ␮ offspring individuals should be deterministically selected to replace the parents, even if the previous minimum in the ﬁtness function has not been lowered. We tried several different functional forms of the “sum of squared errors” as a ﬁtness criterion. It was found that, for simple cases with a small number of parameters, the plain sum of squared errors works quite well. The experimental data are presented to the algorithm as the sampled measured concentrations for each variable. For each variable and for each sample, the difference (or error) between such a measured value and the model output is calculated and squared. The squared errors are then summed over all the samples for each variable and over all variables, and a mean value may be calculated by dividing by the number of variables and the number of samples. The resulting mean square error is then used as the ﬁtness criterion. Brieﬂy, the rules for developing such an algorithm must be followed closely and for each case a judicious choosing of parameters and initial conditions of the algorithm is essential. The best way to perform a correct choice is by trial and error, but some general guides help: 1. Generate a random initial population of individuals. In the integration of the differential equations, the initial conditions have to be chosen while taking care that they are in the physiologically relevant range. The standard deviations for the mutation of the parameters must be chosen near a value of 3, although for some cases a lower value (e.g., 0.1) is better. The initial values of the parameters must also be chosen with the physiological ranges in mind. 2. The self-adaptive parameters ␶ and ␶⬘ explained in the references also have to be chosen carefully. Particularly useful sets have been found to be ␶ ⫽ 1 and ␶⬘ ⫽ 1, ␶ ⫽ 0.1 and ␶⬘ ⫽ 0.1; and, even in extreme cases, where there are many variables or when otherwise the convergence is very slow, ␶ ⫽ 3 and ␶⬘ ⫽ ⫺3. 3. After ranking the best individuals in a population, if the best of them has a lower error than the previous lowest error, choose that new population as the new minimum. In any case, allow the new population to go over the next iteration, so as to allow escaping from local minimums. The new initial condition, if it is being optimized, should be chosen randomly among the best individuals. 4. If possible, use an ordinary differential equation integrator for stiff systems. 5. Delete nonappropriate individuals generated by recombination and mutation before evaluating them, and replace them by appropriate individuals (“appropriate” meaning here that the time course of the variables follows the data in its general ﬂuctuations, e.g., oscillations). RESULTS Simulation and Parameter Identiﬁcation As detailed above, we used the gradient descent method in the initial trials of simulation and parameter identiﬁcation of a genetic network. In the model of Eq. (1) the tested nonlinear functions were sigmoids of the sum of weighted activator/repressor concentrations. If the initial search point was close enough to the optimum, the sum of squared errors diminished, but it usually only approached a local minimum (data not shown). To improve the optimization, global search methods were investigated. An initial trial of an ES adopted the model that yielded the basis for the experimental design of a synthetic oscillatory network of transcriptional regulators.3 This model is based on the Hill function approximation as the nonlinear (repression) transcription function [Eqs. (3)–(5) and (14)]. Simulated data from the model using the parameters that lead to oscillation3 were used in an ES intended to identify these parameters from the data. The results were quite satisfactory, giving a low sum of squared errors and the identiﬁed parameters within the order of magnitude of the real parameters (Fig. 1). The values identiﬁed for the parameters in the scaled version of the model3 are in the range termed “oscillatory” by Elowitz and Leibler3 and are ␣ ⫽ 11, ␣0 ⫽ 0.5, and ␤ ⫽ 3. In order to keep the number of parameters small, the Hill coefﬁcient was not optimized in these preliminary trials of the method with simulated data. The value for the Hill coefﬁcient (nHILL) was instead ﬁxed in these preliminary trials at a value of 2, based on the value employed in published simulations of this system.3 The correctness of this value for the Hill coefﬁcient was corroborated when the method was applied to the experimental data: this time we also optimized the Hill coefﬁcient with the other parameters, and we found that the optimized value for nHILL is 1.9665 (see below). These encouraging results GENOMIC TRANSCRIPTION BIOPHYSICS Fig. 1. (a) Simulated normalized output (six state variables) sampled uniformly 15 times (linear interpolation between samples) as a data set from this network [Eqs. (3)–(5) and (14)] for parameter ﬁtting. (b) Optimized normalized output of the same network. The time in the graph is in units of time intervals between samples. using a global search method led us to seek improved versions for use with larger networks, as well as for use with real data. Trials with simulated data using different types of models proved that the ES algorithm is effective for systems with 10 –20 variables and 30 – 440 parameters, obtaining an average error per sample at the minimum as low as 2.39% (not shown). Parameter Identiﬁcation From Experimental Data The ES algorithm was then improved by introducing self-adaptive features.34 The designed synthetic oscillatory network of transcriptional regulators3 provided the ﬁrst test of the improved ES algorithm, using the data published in ﬁgure 2 of Elowitz and Leibler3 with the ﬂuorescence data proportional to the reporter green ﬂuorescent protein (GFP) concentration. In the experimental setup for obtaining the measurements in the published data,3 the repressor protein lifetimes where brought closer to that of mRNA (about 2 min, on average, in Escherichia coli) by inserting a carboxy-terminal tag at the 3⬘ end of each repressor gene. Proteases in E. coli recognize this tag and target the attached protein for destruction. These tags therefore created the lite variety of repressor proteins, which have a much smaller lifetime (around 4 min) than the wild-type variety. Once their concentration has been increased by an input, it decreases rapidly as the input disappears. In our case it should approach the zero baseline, as found in the simulations. The intermediate stability variant of GFP used in the published experimental setup3 (GFP-AAV) has a half-life of 90 min and the oscillations found in the experiments 345 have a period of about 180 min. Therefore, beginning at a minimum and completing one cycle of oscillation, the concentration of GFP-AAV cannot return to the previous level. This causes an increase of the baseline level of measured ﬂuorescence intensity (superimposed to the actual oscillations) in the chosen cell. This increase is not due to an increase in the baseline concentration of the gene product (repressor protein) whose concentration we want to measure. It is actually an accumulation of GFP-AAV because this protein is degraded more slowly than the repressor proteins. The baseline increase is approximately linear in the time range of the measurements. Therefore, we subtracted this baseline by subtracting a ramp with the corresponding baseline slope in order to better represent the actual concentration of the gene product (repressor protein) in the synthetic oscillatory network of transcriptional regulators.3 The resulting ramp baseline-subtracted oscillation data have now a near zero baseline, which agrees more with real, rapidly decaying concentrations of lite repressor proteins in the synthetic system. It is important for future uses of the method to have measurements of the mRNA concentrations (and eventually protein concentrations) in molar (nanomolar) concentration units because in this way the appropriate threshold parameter in nanomolar units can be ﬁtted in a unique way. This corresponds to the Kdiss of Eq. (9) in a simple binding reaction, Kdiss thus being a measure of the necessary concentration for achieving a switch from low to high occupancy and from low to high activation or repression of transcription. In the case of the thermodynamically derived function, this permits the ﬁnding of the free energies of binding (see below and Fig. 2). The free energies of binding are important because they can be directly linked to the structural details of the repressor protein binding to DNA. Posttranslational structural modiﬁcations of the transcription factor proteins or covalent modiﬁcations of DNA are often used by the cell to control transcription, for example, protein phosphorylations, DNA methylations, and other modiﬁcations of the transcription machinery. Changes in the free energies of binding may also be correlated with the effect of mutations and singlenucleotide polymorphisms. It could also be possible that some other structural modiﬁcations like redox modiﬁcations caused by redox active molecules such as nitric oxide and reactive oxygen species (e.g., thiol residue modiﬁcations, nitrosylation, or free radicals in some amino acids) might be identiﬁed by this procedure as effectors of disease (manifest in malfunction of the intracellular network). Because in the synthetic system data3 and in the usual microarray data there is always an arbitrary scale factor for the magnitude of the measurements, we analyzed the problem by incorporating this scaling factor as a parameter in the identiﬁcation procedure and getting a measurement of the sensitivity of the parameter identiﬁcation to it. It was found that generally for the ranges studied, the results of the identiﬁcation show low sensitivity to the scaling factor. However, in the future it would be better to have some way of exactly calibrating the measurements to nanomolar concentrations, for example, performing (at 346 G. CAVELIER AND D. ANASTASSIOU Fig. 2. (---) Experimental baseline-subtracted data of the synthetic3 network (six state variables, only one measured) sampled 60 times at equal intervals as a data set for parameter ﬁtting. The network is modeled as in Eqs. (3)–(5) with the following nonlinear functions: (a) hyperbolic tangent [Eq. (2)]; (b) Hill function [Eq. (14)] also ﬁtting the Hill coefﬁcient; (c) O-model function [Eq. (15)]; (d) thermodynamically derived function [Eq. (19)]: (—) Optimized output of the same networks. The time is the real time in minutes. least for a small subset of the mRNAs) complementary quantitative real-time polymerase chain reaction measurements at the same time that the microarray measurements are done and with the same samples. The usual mRNA concentrations are expected to be in the range of 100 nM, and we had in the available published data3 ﬂuorescence units in the numerical range of 100s. We therefore assumed that the units were in nanomoles. These data were obtained experimentally by other groups, so we do not have ways of performing the calibration and can only assume that they are in the range of nanomolar units and test if there is sensitivity to their actual magnitude. We empirically tested several scaling factors and found that, in general, the sensitivity of the parameter identiﬁcation to such scaling factors is not very high, so we decided that the units used for ﬂuorescence roughly corresponded to the value for the nanomolar concentration. Proper calibration will remove this particular uncertainty from the parameter-identiﬁcation problem in future work. We then carried out several parameter identiﬁcations in succession: an algebraic sigmoid, the Hill model, the so-called O-model,33 and the thermodynamically derived model, using the same published data in all cases.3 The results are shown in Figure 2(a– d). In Figure 2(a) the hyperbolic tangent [Eq. (2)] is used as a normalized, nonlinear repression function of the protein concentrations, rather than the function f(P) of Eq. (14), but the remaining parameters are the same as in Eqs. (3)–(5). In Figure 2(b) the optimization with the Hill [Eq. (14)] nonlinear function is also shown, this time optimizing among the parameters, the Hill coefﬁcient as well. The algorithm found an optimized nHILL value of 1.9665. The same optimization results are shown in Figure 2(c) with the O-model function [Eq. (15)] and in Figure 2(d) with the thermodynamically derived nonlinear function [Eq. (19) and Table I]. The values obtained for the free energies were remarkably well positioned in the range of free energies expected for the binding of transcription factors (repressor dimers) to DNA promoters, as well as for cooperativity interactions between adjacent DNA-bound repressors.11,12 The full set of parameters (the initial conditions were found as parameters by the identiﬁcation routine, because they are not available in the published data,3 but they are not shown) extracted by the optimization for the curves in Figure 2 is shown in Table II. It should be noted that we did not include the equations for transcription and translation of GFP, and the model here uses only one (ﬁxed, optimized previously by trial and error) value for the following four parameters: the three protein’s translation constants, the three mRNA’s degradation constants, the three proteins’ degradation constants, and the three thresholds for half-maximal repression in the nonlinear functions. In addition, we did not take the issue of dimerization of repressor monomers into account.33 Some of these restrictions were eliminated in the processing of the yeast data discussed below. Figure 3 compares the normalized nonlinear function identiﬁed for the Hill model [Fig. 2(b)] with nHILL ⫽ 1.9665 and the thermodynamically derived nonlinear function with two operators after parameter identiﬁcation [Fig. 2(d)]. We ﬁrst ﬁtted the Hill model and found by trial and error an initial estimate for the free energies that brought the thermodynamically derived function very close to the ﬁtted Hill function with nHILL ⫽ 1.9665. These free energy estimates provide an ideal starting point for the parameter identiﬁcation with the thermodynamically derived nonlinear function of Eq. (19). They were slightly corrected for 347 GENOMIC TRANSCRIPTION BIOPHYSICS TABLE II. Parameters Found in Different Models for Synthetic Network 再冋冉冊册冎 1 P ⫺1 ⫹1 tanh ⫺S䡠 2 THR E min ⫽ 18.28% THR ⫽ 10 nM ␣ ⫽ 5.32 ␣0 ⫽ 0.03 ␶m ⫽ 10.00 ␶P ⫽ 17.86 S ⫽ 2.58 1 P 1⫹ THR 冉冊 nHILL E min ⫽ 16.25% THR ⫽ 10 nM ␣ ⫽ 7.76 ␣0 ⫽ 0.02 ␶m ⫽ 14.00 ␶P ⫽ 15.22 nHILL ⫽ 1.97 1 Pd 1⫹ THR 冉冊 o E min ⫽ 19.29% THR ⫽ 10 nM ␣ ⫽ 9.95 ␣0 ⫽ 0.02 ␶m ⫽ 14.00 ␶P ⫽ 13.46 o ⫽ 2.74 Fig. 3. (—) Thermodynamically derived nonlinear function [Eq. (19), Fig. 2(d)] calculated with the free energies of binding and cooperativity determined by the ﬁtting (kcal/mol, Table II). (---) Nonlinear function used in the optimization shown in Figure 2(b) (Hill function with Hill coefﬁcient nHill ⫽ 1.9665). k4䡠关RNAP兴⫹共k6⫹k7兲䡠关RNAP兴䡠x⫹k8䡠关RNAP兴䡠x2 共k1⫹k4䡠关RNAP兴兲⫹共k2⫹k3⫹共k6⫹k7兲䡠关RNAP兴兲䡠x⫹共k5⫹k8䡠关RNAP兴兲䡠x2 E min ⫽ 21.87% ⌬G2 ⫽ ⫺8.59 kcal/mol THR ⫽ 10 nM ⌬G3 ⫽ ⫺9.24 kcal/mol ␣ ⫽ 12.46 ⌬G4 ⫽ ⫺5.71 kcal/mol ␣0 ⫽ 0.10 ⌬G23 ⫽ ⫺4.86 kcal/mol ␶m ⫽ 17.50 ⌬G24 ⫽ ⫺0.29 kcal/mol ␶P ⫽ 10.12 ⌬G34 ⫽ ⫺0.35 kcal/mol ⌬G234 ⫽ ⫺1.27 kcal/mol to four pairs of yeast genes gives the results shown in Table III and the graphs of Figure 4. The four respective pairs of genes and the causality of the connectivity found are YBR083W activator of YHR084W, YBL008W repressor of YML010W, YAL021C activator of YHR119W, and YBL005W activator of YGL013C. Repeated application of the algorithm for one pair of genes gives on average the same values for the parameters. We analyzed more than 10 pairs of yeast genes in this way and the results are consistently good (not shown). It was found that some of the paths in the connectivity are in the opposite direction to that shown in the published data.9 We believe that the optimization done here gives a more accurate representation of causality in these cases, because the time course of some processes represented in detail by the transcription and translation events in the equations cannot be reversed. Attempts to do so produced a very bad ﬁt to the data (not shown). DISCUSSION the interaction free energies by the optimization routine to get the ﬁnal values shown in Table II. As a further proof that the algorithm works with experimental data, we tested it using results from high-density oligonucleotide microarrays of yeast meiotic expression time-series patterns. The data used here were from the yeast strain SK1 MATa/alpha (http://re-esposito.bsd. uchicago.edu/), which was published by Primig et al.13 The connectivities used to begin the parameter identiﬁcation were downloaded from the http://www.nature.com/ng/ journal/v31/n1/suppinfo/ng873_S1.html and the method for ﬁnding them can be found in Ref. 9. These connectivities can be found equivalently by other methods.16,17,19,38 The model we used was similar to the Hill model of Eqs. (3)–(5) and (14). The difference in this case is that one of the experimental mRNA concentrations in time, namely m1(t), is taken as independent input for translation of its gene product (P1) using a form similar to Eq. (4). Using f(P1) in a form similar to Eq. (3), the output of this model (m2, solid red curves, Fig. 4) is then found. It is correspondingly ﬁtted to the yeast strain SK1 MATa/alpha sampled data mRNA time series (m2, dashed red curves, Fig. 4) with the ES algorithm, ﬁnding the whole set of parameters and P1 (t ⫽ 0). The time constants ␶m and ␶P are in minutes, and THR and P1 (t ⫽ 0) are in nanomoles. The other units are the same as previously. This ﬁtting applied Meaningful biophysical parameters were derived for a model of transcriptional regulatory networks from experimental data found with the tool described here. It is signiﬁcant that the free energies of binding of transcription factors to their promoters can be determined, as can any cooperativity effects between distant, operator-bound transcription factors. It will be particularly interesting to probe these free energies of binding and cooperativity in cases where different cellular functions, mutations, covalent modiﬁcations, or misfolding of proteins are involved. Models with differing architecture and complexity may be evaluated and compared via the same procedure. The Hill coefﬁcient in the Hill model or the o coefﬁcient in the O-model also clearly indicate how to construct a more detailed, thermodynamically based model that is best suited for each particular case. An analysis of the parameters found (Table II) for the synthetic network experimental data shows that the various test models generally yield values of the same order of magnitude. The relatively small differences are most likely due to the different nonlinear functions that were used. The Hill coefﬁcient and its analogous quantities S and o are consistently close to 2, as expected for a system with some cooperativity.12,33 We did not include the possibility of transcription factor dimerization in the modeling of the Hill function, the O-model, and the thermodynamically derived function. 348 G. CAVELIER AND D. ANASTASSIOU Fig. 4. Yeast strain SK1 MATa/alpha (---) time-series data and (—) its ﬁtted-model curves for four pairs of genes. Time-series data for (---) m1 ⫽ mRNA1 and (---) m2 ⫽ mRNA2. The ﬁrst is independent input for translation of the corresponding gene product, protein P1. Protein P1 is in turn transcriptional activator or repressor of (—) m2 ⫽ mRNA2 (optimized output of the model). The Hill model is used [Eqs. (3)–(5) and (14)]. Pairs of yeast genes: (a) BR083W activator of YHR084W, (b) YBL008W repressor of ML010W, (c) YAL021C activator of YHR119W, (d) YBL005W activator of GL013C. TABLE III. Hill Model Parameters for Yeast Genes ␣ ␣0 ␶m ␶P kt nHILL THR P1(0) Emin YBR083W YHR084W YBL008W YML010W YAL021C YHR119W YBL005W YGL013C -% -% -% -% For more accurate results, this should be included and may be accomplished easily.33 Dimerization may increase the cooperativity effect of the overall function; thus, in detailed models such as the thermodynamically derived function its effect should be evaluated, and we plan to do so in future work. The value of the threshold at which there is half-maximal repression/activation by the transcription factor in the cases shown in Figures 1–3 was optimized, but this was done empirically. In the last runs of the method with yeast data (Table III and Fig. 4) the threshold was indeed optimized by the ES tool, as can be seen in the values that were found. The average error per sample shown in Tables II and III is close to 20%, which is a low value considering the small number of samples available and other difﬁculties such as noise, lack of protein concentration measurements, lack of initial conditions of many variables, lack of precise calibra- tion to nanomolar concentrations, and inherent imprecision because of simpliﬁcations of the model. The free energies of binding and cooperativity found by our tools are in the same physiological ranges as the experimentally derived energies of binding and cooperativity in the bacteriophage lambda system.12 In terms of the cooperativity free energies obtained from the thermodynamically derived model in the studied case, only one (⌬G23) appears to have a meaningful value and the others are too low for any real effect. Although much less detail is present on the model description of degradation, the time constants ␶m and ␶P found in the different cases are in the physiological ranges as well. It has been found that intrinsic noise cannot be separated from normal functioning of genetic and biochemical networks.39,40 Our models can be described by a set of biochemical reactions from the dimerization of transcription factors to the DNA binding reaction of these transcription factors and RNAP. Transcription and translation kinetics are in turn represented by ordinary differential equations in which the rate of transcription is proportional to the fractional binding occupancy of the promoter by RNAP. In the steady state of these binding reactions, the fractional occupancy will be a multidimensional function of the different activator and repressor concentrations, with parameters according to the model chosen among the ones discussed here or more complex ones. Therefore, the whole description of transcription/translation in genetic regulatory networks amounts to a biochemical set of equations plus the kinetic equations. The kinetic equations can also be interpreted as a special case of biochemi- 349 GENOMIC TRANSCRIPTION BIOPHYSICS cal equations for the corresponding reactions of (enzymatic) mRNA and protein production and degradation. As discussed above, many of these are usually taken as linear, with the exception of the nonlinear functions for mRNA transcription. Algorithms for simulation of biochemical equations with intrinsic noise are well known.40,41 Among them, the StochSim approach41 is interesting because it relies only on the biochemical equations for solving the system behavior stochastically. Because our parameter identiﬁcation is independent of the model chosen for the system, we can use the set of biochemical and kinetic equations discussed above in a StochSim stochastic model with intrinsic noise as a better representation of Eqs. (3)–(5) in order to solve for the time-series simulations in the model. The parameter-identiﬁcation routine can then be applied to ﬁnd the parameters. An analysis of this approach together with the deterministic continuous approach based on the same equations (as performed here) will permit a comparison of both cases and determine if the stochastic model with intrinsic noise yields similar parameters, has the same or different sensitivities to changes, and is better in cases involving few molecules and small volumes, as should be expected.40,41 CONCLUSION We can ﬁnd the model parameters from available yeast (or other genomes) time-series data by choosing appropriate subsystems of the whole genome and by using connectivity data from independent algorithms as a primer. The simplest of these subsystems correspond to regulated regulatory genes, giving a pair of mRNAs and the protein product of the ﬁrst gene, which regulates the second gene. By repeating this procedure with other subnetworks across the entire genome of an organism, one could construct a database of parameters and a corresponding model of increasingly complex subnetworks from the accumulated results with a reasonable computation effort. Comparing the healthy and diseased states of cells with the parameters identiﬁed here will be useful in drug design and gene therapy applications. Analyses of the changes in degradation rates, free energies of binding, or cooperativity effects are important for disease mechanism elucidation, as well as for diagnosis. In addition, because the network equations are predictive of gene concentrations in time under different conditions, the model analysis displays the capacity for disease prognosis. The instrument developed here also provides a detailed description of the causality and time dependence of the network’s protein and mRNA concentrations, as opposed to the static pictures obtained with other tools. ACKNOWLEDGMENTS The authors thank C. Zukowski, I. Tagkopoulos, A. Rzhetsky, H. Bussemaker, and C. Wiggins for their assistance and useful discussions and the anonymous reviewers for very helpful comments that improved the manuscript. REFERENCES 1. Hasty J, McMillen D, Isaacs F, Collins JJ. Computational studies of gene regulatory networks: in numero molecular biology. Nat Rev Genet 2001;2:268 –279. 2. Hasty J, McMillen D, Collins JJ. Engineered gene circuits. Nature 2002;-):224 –230. 3. Elowitz MB, Leibler S. A synthetic oscillatory network of transcriptional regulators. Nature 2000;-):-. Gardner TS, Cantor CR, Collins JJ. Construction of a genetic toggle switch in Escherichia coli. Nature 2000;-):339 –342. 5. McAdams HH, Shapiro L. Circuit simulation of genetic networks. Science 1995;-):650 – 656. 6. Gilman A, Arkin AP. Genetic “code”: representations and dynamical models of genetic components and networks. Annu Rev Genom Hum Genet 2002;3:-. Quackenbush J. Computational analysis of microarray data. Nat Rev Genet 2001;2:418 – 427. 8. Brazma A, Vilo J. Gene expression data analysis. Microbes Infect 2001;3:-. Guelzim N, Bottani S, Bourgine P, Kepes F. Topological and causal structure of the yeast transcriptional regulatory network. Nat Genet 2002;31:60 – 63. 10. Vohradsky J. Neural model of the genetic network. J Biol Chem 2001;276:36168 -. Ackers GK, Johnson AD, Shea MA. Quantitative model for gene regulation by lambda phage repressor. Proc Natl Acad Sci USA 1982;79:1129 –1133. 12. Shea MA, Ackers GK. The OR control system of bacteriophage lambda. A physical– chemical model for gene regulation. J Mol Biol 1985;181:-. Primig M, Williams RM, Winzeler EA, Tevzadze GG, Conway AR, Hwang SY, Davis RW, Esposito RE. The core meiotic transcriptome in budding yeasts. Nat Genet 2000;26:-. Chu S, DeRisi J, Eisen M, Mulholland J, Botstein D, Brown PO, Herskowitz I. The transcriptional program of sporulation in budding yeast. Science 1998;-):699 –705. 15. Holstege FC, Jennings EG, Wyrick JJ, Lee TI, Hengartner CJ, Green MR, Golub TR, Lander ES, Young RA. Dissecting the regulatory circuitry of a eukaryotic genome. Cell 1998;95:-. Bussemaker HJ, Li H, Siggia ED. Regulatory element detection using correlation with expression. Nat Genet 2001;27:-. Lee TI, Rinaldi NJ, Robert F, Odom DT, Bar-Joseph Z, Gerber GK, Hannett NM, Harbison CT, Thompson CM, Simon I, Zeitlinger J, Jennings EG, Murray HL, Gordon DB, Ren B, Wyrick JJ, Tagne JB, Volkert TL, Fraenkel E, Gifford DK, Young RA. Transcriptional regulatory networks in Saccharomyces cerevisiae. Science 2002;-):799 – 804. 18. Wyrick JJ, Young RA. Deciphering gene expression regulatory networks. Curr Opin Genet Dev 2002;12:130 –136. 19. Kholodenko BN, Kiyatkin A, Bruggeman FJ, Sontag E, Westerhoff HV, Hoek JB. Untangling the wires: a strategy to trace functional interactions in signaling and gene networks. Proc Natl Acad Sci USA 2002;99:-. Chen T, He HL, Church GM. Modeling gene expression with differential equations. Proc of Pac Symp Biocomput 1999:29 – 40. 21. Vohradsky J. Neural network model of gene expression. Faseb J 2001;15:846 – 854. 22. Reinitz J, Sharp DH. Mechanism of eve stripe formation. Mech Dev 1995;49:-. D’Haeseleer P, Liang S, Somogyi R. Genetic network inference: from co-expression clustering to reverse engineering. Bioinformatics 2000;16:-. Bolouri H, Davidson EH. Modeling transcriptional regulatory networks. Bioessays 2002;24:1118 –1129. 25. van Holde KE, Johnson WC, Shing Ho P. Principles of physical biochemistry. Upper Saddle River, NJ: Prentice Hall; 1998. 26. Ferrell JE, Machleder EM. The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes. Science 1998;-):-. Huang CYF, Ferrell JE. Ultrasensitivity in the mitogen-activated protein kinase cascade. Proc Natl Acad Sci USA 1996;93:10078 -. Pray TR, Burz DS, Ackers GK. Cooperative non-speciﬁc DNA binding by octamerizing lambda cI repressors: a site-speciﬁc thermodynamic analysis. J Mol Biol 1998;282:947–958. 350 G. CAVELIER AND D. ANASTASSIOU 29. Perlmutterhayman B. Cooperative binding to macromolecules—a formal approach. Acc Chem Res 1986;19:90 –96. 30. Darling PJ, Holt JM, Ackers GK. Coupled energetics of lambda cro repressor self-assembly and site-speciﬁc DNA operator binding II: cooperative interactions of cro dimers. J Mol Biol 2000;302:-. Rusinova E, Ross JB, Laue TM, Sowers LC, Senear DF. Linkage between operator binding and dimer to octamer self-assembly of bacteriophage lambda cI repressor. Biochemistry 1997;36:12994 -. Wolf DM, Eeckman FH. On the relationship between genomic regulatory element organization and gene regulatory dynamics. J Theor Biol 1998;195:-. Elowitz ME. Transport, assembly, and dynamics in systems of interacting proteins. Princeton, NJ: Princeton University Press; 1999. 34. Beyer HG, Deb K. On self-adaptive features in real-parameter evolutionary algorithms. IEEE Trans Evol Comput 2001;5:250 – 270. 35. Back T, Hammel U, Schwefel H-P. Evolutionary computation: comments on the history and current state. IEEE Trans Evol Comput 1997;1:3–17. 36. Iba H, Mimura A. Inference of a gene regulatory network by means of interactive evolutionary computing. Inform Sci 2002;145:-. Back T. Evolutionary algorithms in theory and practice: evolution strategies, evolutionary programming, genetic algorithms. New York: Oxford University Press; 1996. p xii. 38. Tegner J, Yeung MK, Hasty J, Collins JJ. Reverse engineering gene networks: integrating genetic perturbations with dynamical modeling. Proc Natl Acad Sci USA 2003;100. 39. McAdams HH, Arkin A. It’s a noisy business! Genetic regulation at the nanomolar scale. Trends Genet 1999;15:65-69. 40. Arkin A, Ross J, McAdams HH. Stochastic kinetic analysis of developmental pathway bifurcation in phage lambda-infected Escherichia coli cells. Genetics 1998;149:-. Morton-Firth CJ, Bray D. Predicting temporal ﬂuctuations in an intracellular signalling pathway. J Theor Biol 1998;192:117–128.