Abstract
We propose a new model for describing gene regulatory networks that can capture discrete (Boolean) and continuous (differential) aspects of gene regulation. After giving some illustrations of the model, we study the problem of the reverse engineering of such networks, i.e., how to construct a network from gene expression data. We prove that for our model there exists an algorithm finding a network compatible with the given data. We demonstrate the model by simulating lambdaphage. We also describe some generalizations of the model, discuss their relevance to the realworld gene networks and formulate a number of open problems.
Keywords:
gene regulation; regulatory networks; regulatory circuits; dynamic systems; finite state automata; reverse engineeringBackground
There are many mechanisms how genes are regulated. An important role in gene regulation apparently is played by specific proteins, called transcription factors, which influence the transcription of particular genes by binding to specific parts of the DNA in the genome. In this way a product of one gene can influence the expression of another gene, and we can consider a network of gene regulation. Such regulatory networks or circuits are well studied in lambdaphage and some other viruses [1]. If the network involves only few genes, its functioning can be understood relatively directly. But what does it mean to understand a gene regulatory network of hundreds or thousands of genes? Just describing such a network may be highly nontrivial. We think that to be able to understand complex gene regulatory networks, first a formal language for describing such networks has to be developed. The language can be graph based and preferably should allow the simulation of the behaviour of the network. By simulating a network we can make predictions and compare them to experimental data. If the predictions are consistent with the data, then we can say that the model is correct (within the given accuracy limits). Such an approach is usual in physics: models (theories) are built to explain existing data, then predictions are made, which again are compared to new data. If the correspondence is good, it is claimed that the phenomenon has been understood. Preferably, the model should not be a black box, but should be interpretable, and ideally its elements should have interpretation in the real world consistent with the existing knowledge. At the same time, each model involves a simplification of the real world, which is a part of the strength of the modelling approaches.
Various models for gene regulatory networks have been proposed and studied (see for instance [2,3]). In general these models fall into two categories: boolean network based models, for instance [46], and dynamic systems described by differential or difference equations, for instance [7,8]. Each of these models have their advantages and drawbacks. The Boolean model is based on the assumption that the important aspects of gene regulation can be described by binary on/off switches, functioning in discrete time steps: the state of the network in time point n is determined by its state at timepoint n1. Even if we generalize these models to more than two discrete states they cannot describe continuous changes that happen in the cell environment. These can be described by differential equation based models, which on the other hand cannot easily describe the discrete aspects of gene regulation such as binding of a transcription factor to the DNA, which is essentially an on/off event. Also, in a differential equation model it is difficult (though not impossible) to describe nonadditive logics in gene regulation (for instance, competitive events), as well as time delays.
Models trying to combine the discrete and continuous components have been proposed, for instance in [911]. Thomas and Thieffry [12,13] describe a combined model for qualitative description of gene regulatory networks. They introduce a notion of gene state and image, the last effectively representing the substance produced by the respective gene. There is a time delay between the change of the gene state and the change of the image state. By introducing different levels of gene activity and thresholds for switching the gene states, thus they go beyond binary models. They study the qualitative behaviours of various feedback loops in their model, and show that they fall into two classes: positive loops leading to multi stable states and negative ones leading to periodicity.
The finite state linear model proposed in this paper combines the discrete and continuous aspects of gene regulation in a simple and structured way. It has a boolean network type discrete control component, and an environment of substances changing their concentrations continuously. Time is continuous, and the state of the network directly determines only the concentration change rates, while the state is affected by the concentrations themselves.
A framework (a formal language) for describing gene regulatory networks enables us to study the problem of building particular models from gene expression dataoften referred to as the reverse engineering of gene networks (e.g., [3,5]). Until recently there were little quantitative data available for building models for gene regulation. Most of the earlier gene network models, including [13] are based on observations from gene mutation data leading to phenomenological changes and not on direct observations of gene activities. This has changed with the advent of DNA microarray technology, which generates huge amounts of data characterizing gene activities under various conditions [1416] and are now being collected in various databases [17]. There can be various precise formulations for the reverse engineering problem, and there is a certain analogy between the problems of reverse engineering of gene networks and the problem of identifying finite state automata from input/output data [18].
In this paper we consider two different formulations of the reverse engineering problem. The weakest one is finding a gene network consistent with the given data. We prove that this problem is algorithmically solvable for our model. The second one involves assuming that the data have been produced by some unknown gene network, which we want to reconstruct by making experiments. This problem is still open. In the next section we describe the model, after which we study the reverse engineering problem. Then we give some informal extension of the model, and use it to describe the lambdaphage regulatory circuit. Finally we discuss some open problems.
Results and Discussion
The definition of the model
The assumptions on which our model is based are: (1) the gene activity is determined by the state of transcription factor binding sites in its promoter region; (2) each binding site can be in one of a finite number of states, characterized by having or not having bound a particular transcription factor; (3) depending on the states of the binding sites in the promoter, the gene can either be silent, or have a particular activity level; (4) if a gene is active, the concentration of the substance it produces is growing with a rate dependent on the activity level of the gene, otherwise it is decreasing (or staying 0); (5) the state of a binding site depends on the concentration of the respective transcription factor(s). To make these assumptions precise and to formalize them we have developed the model described below. We begin by describing a simpler version of the model, which we call the binary model, where each binding site and each gene have only two states: on or off. We formulate the reverse engineering problem for the binary model, before introducing the general case, though the formulation remains the same in the general case.
The binary model
Informally we assume that we have an environment of n substances 1, ..., n having concentrations c_{1}(t), ..., c_{n}(t), respectively, which may change in time t. We also assume that there are, what we call substance binding sites in the environment, each of which can attach (bind) a specific substance. In the binary case the binding site can bind only one substance. We define a binary binding site b as a triple
b = (i, a, d),
where i is the number of the substance (which can bind to b), and a and d are positive real constants 0 < d < a, called association and dissociation constants, respectively. Each binding site can be in one of two states: attached state or detached state. If binding site b = (i, a, d) is in detached state, and the concentration of substance i reaches the association constant a, i.e., c_{i}(t) ≥ a, then the b switches to attached state. If b is in attached state and the concentration c_{i}(t) falls below the dissociation constant d, i.e., c_{i}(t) ≤ d, then b switches to detached state. We denote the attached state by 1 and detached state by 0. Thus, the binding site can be described as a two state automaton in Figure 1, left. Next we define a binary gene. Each binary gene produces one substance. A binary gene can have two states on or off, depending on the state of the binding sites regulating this gene. If a gene G is on, then the respective substance is being produced and its concentration linearly increases. If G is off, the substance is being degraded by the environment, and its concentration linearly decreases (until it reaches 0, or the gene switches on). Formally a binary gene is a triple
Figure 1. Finite state automata describing a binary binding site (left) and a multilevel binding site (right).
G = (B, F, r),
where B = (b_{1}, ..., b_{k}), and b_{1}, ..., b_{k }are a subset of the binding sites, F is a boolean function called control function, and r = (i, r_{0}, r_{1}), where i is an integer denoting the number of the substance produced by the gene, r_{0 }< 0 is a real constant called degradation rate, and r_{1 }> 0 production rate. We call r a substance generator. Graphically, a gene is represented as in Figure 2, left. We can think of the binding sites and the control function, as the promoter of the gene, while the substance generator  as the coding part plus transcription machinery.
Figure 2. Left: a graphical representation of a gene. The triangles on the left represent the binding sites b_{1}, b_{2}, b_{3}. The rectangle in the middle represent the control function (in the particular example F(x_{1},x_{2},x_{3}) = x_{1 }& x_{2 }& ¬ x_{3}, meaning that the gene is on if and only if the first two binding sites are in attached state, while the third in the detached state), and the diamond on the right represents the substance generator. Right: an example gene network. In this network Γ = {G_{1}, G_{2}}, G_{1 }= ((b_{1},b_{2}), F_{1}, r_{1}), G_{2 }= ((b_{3}), F_{2}, r_{2}), b_{1 }= (1, a_{1}, d_{1}), b_{2 }= (2, a_{2}, d_{2}), b_{3 }= (1, a_{3}, d_{3}), r_{1 }= (1, r_{0,1}, r_{1,1}), and r_{2 }= (2, r_{0,2}, r_{1,2}). The solid lines can be regarded as connecting the substance produced by the gene to the respective binding sites, while the dotted lines channelling the information about the states of the binding sites and genes. Another interpretation of the lines is that the solid lines transmit real numbers, while dotted ones  boolean values.
The semantics of a gene G = (B, F, r) can be described as follows. Let q_{1}, ..., q_{k }be the states of binding sites b_{1}, ..., b_{k }where B = (b_{1}, ..., b_{k}): i.e., q_{i }= 1 if b_{i }is in attached state, and q_{i }= 0, otherwise, at some given time point t'. If F(q_{1}, ..., q_{k}) = 1, i.e., the gene is on, then the concentration c_{i}(t) of substance i (where r = (i, r_{0}, r_{1}) increases in time with rate r_{1}, i.e., c_{i}(t) = c_{i}(t') + (t  t')r_{1}. If F_{i}(q_{1}, ..., q_{k}) = 0, i.e., the gene is off, then, the concentration c_{i}(t) decreases with rate r_{0 }while it is positive, or remains equal to 0.
A binary gene network
We define a gene network as a set of genes
Γ = G_{1}, ..., G_{n}.
We can use a graphical representation of gene networks to show which gene products can attach to which binding sites. An example of such representation is given in Figure 2, right.
In general, several genes my share the same binding site (in graphical representation the dotted line coming out of a binding site can fork to several control functions). To describe the functioning of a gene network let us consider an example in Figure 3 (a more formal definition is given in Section 4.1).
Figure 3. The functioning of a simple network of two binary genes with a negative feedback loop.
Let Γ_{1 }= (G_{1}, G_{2}), G_{1 }= ((b_{1}), F_{1}, r_{1}), G_{2 }= ((b_{2}), F_{2}, r_{2}), and let us assume that the function F_{1 }is the negation (i.e., F_{1}(0) = 1 and F_{1}(1) = 0), while F_{2 }is the identity (i.e., F_{2}(0) = 0 and F_{2}(1) = 1). Gene G_{1 }produces substance 1, gene G_{2 }substance 2, and let b_{1 }= (2,a_{1},d_{1}), b_{2 }= (1,a_{2},d_{2}), r_{1 }= (1,r_{1,0},r_{1,1}), and r_{2 }= (1,r_{2,0}, r_{2,1}).
Further, we assume that at time point t_{0 }= 0 the substance 1 has some positive initial concentration c_{1}(t_{0}) > 0, while c_{2}(t_{0}) = 0, as shown in the graph in the lower part of Figure 3. We also assume that the states of both binding sites are initially equal to 0, i.e., q_{1 }= 0, q_{2 }= 0. Starting from this state at t_{0}, the network Γ_{1 }functions as follows. Since F_{1}(0) = 1, the substance 1 is produced with rate r_{1,1 }> 0, and the concentration c_{1}(t) is growing. On the other hand F_{2}(0) = 0, therefore the concentration c_{2}(t) remains 0. This linear change continues until time t = t_{1}, when c_{1}(t) = a_{2}, i.e., until the concentration of the substance 1 reaches the association constant for binding site b_{2}. At that point b_{2 }switches to attached state 1, and since F_{2}(1) = 1, gene G_{2 }switches to on state and starts producing substance 2 with rate r_{2,1}. Thus, starting from t = t_{1}, the concentration of both substances are growing. This continues until the c_{2 }reaches a_{1}, at which point b_{1 }switches to on state, switching gene G_{1 }off. The concentration c_{1}(t) starts falling, and when it reaches d_{2}, gene G_{2 }switches off and c_{2}(t) starts falling too. This continues as shown in Figure 3. The table at the bottom of Figure 3 show the states of the binding sites.
The assumption that the substance concentrations change linearly for the given state is not essential for the model. We think that linearity may be a reasonable approximation in the cases where the gene expression rates are far from saturation levels. This assumption can be relaxed by changing the linear functions to a function that behave approximately linearly while the values are relatively small, decreasing the growth rate for larger values and asymptotically approaching some given maximum. An example of such a function is the solution of the logistic differential equation dc/dt = rc(1c/k), where c is the concentration, and r and k are constants.
Another instance where the linearity may be insufficient, is if the degradation rate of a certain substance depends on the concentration of another substance (for instance, if one substance is degrading the other). Our model can be generalized to capture this situation in a straight forward manner, if there are no loops in the dependency graph describing which substances degrade which.
Although the linearity is not an essential feature of the model, in the next sections dealing with the reverse engineering, we will stick to this assumption, as we think that the properties of a simpler model should be explored first.
Reverse engineering of gene networks
Let b_{1}, ..., b_{m}, be all the binding sites in the environment, and let Q(t') = (q_{1}(t'), ..., q_{m}(t')) be their states at time point t'. We call Q(t') the binding site state vector of the network at time point t'.
Let C(t') = (c_{1}(t'), ..., c_{n}(t')) be the concentrations of all environment substances at time point t'. We call C(t') the environment concentration vector. We say that the binding site state Q(t') and concentration state C(t') are compatible, if for every binding site b_{j }= (i, a_{j}, d_{j}), if q_{j }= 0 then c_{i }< a_{j}, and if q_{j }= 1 then c_{i }> d_{j}. We define the network state vector as a pair
Σ(t') = (Q(t'), C(t'))
and we say that it is compatible if Q(t') is compatible with C(t'). We often omit t'.
Note that concentration state vector C(t') = (c_{1}(t'), ..., c_{n}(t')) at a given timepoint t' can be regarded as a concentration measurement. Let us define a measurement series as a pair of mtuples
M = ((t_{0}, t_{1}, ..., t_{m}), (C(t_{0}), C(t_{1}), ..., C(t_{m}))).
The reverse engineering problem for gene networks can be formulated as follows:
given a measurement series M = ((t_{0}, t_{1}, ..., t_{m}), (C(t_{0}), C(t_{1}), ..., C(t_{m}))), find a gene network Γ that can produce concentrations C(t_{0}), C(t_{1}), ..., C(t_{m}) at time points t_{0}, t_{1}, ..., t_{m}. In this case we say that network Γ is compatible with measurements M.
Theorem
The problem of reverse engineering is algorithmically solvable for the linear finite state gene network models, i.e., there exists an algorithm that, given a series of measurements M, outputs a gene regulatory network Γ compatible with M.
To prove the theorem, we need to introduce a few auxiliary notions. Given a network Γ and a compatible starting state Σ(t_{0}), network Γ defines the concentration change graph Δ, which is the set of all points C(t) = (c_{1}(t), ...,c_{n}(t)), for the time interval t ∈ [t_{0}, ∞]. An example of an initial part of such a graph is given in the lower part of Figure 3 and in Figure 4. Note that each concentration changes as a piecewise linear function.
Figure 4. The environment change graph
Let Γ = {G_{1}, ..., G_{n}}be a network, where G_{i }= (B_{i}, F_{i}, r_{i}). Let us consider the sets of all the binding sites in the environment and all the substance generators in the network. Each binding site and each substance generator depends on two real value constants (association and dissociation constants for binding sites, and production and degradation constants for substance generators). Let us denote the set of all binding site constants in the network by β, and the set of all substance generator constants by γ. Let α = β ∪ γ, and we call α the set of the network constants.
Let us consider an initial part Δ(t_{0},t') of a concentration change graph Δ for a network Γ in time interval [t_{0},t']. The slopes of the linear parts in the graph are determined by a subset of γ, while the transitionpoints by a subset β. We denote these subsets by γ' and β'. We call α' = β' ∪ γ' the set of reachable constants for the network Γ in [t_{0},t'] for the given starting state.
Finally, for a given network Γ, we define the network structure as the object obtained from Γ by ignoring all the network constants (formally, we can substitute all the constants in Γ, for instance, by 0). In the graphical representation the network remains the same, but the constants disappear. The control functions are a part of the structure.
Now, to prove the theorem, first, note that given an initial part of a concentration change graph Δ(t_{0},t'), we can find all reachable constants β' and γ'. We also know the number of the genes in the network, which equals n. We know the maximal number of binding sites that can switch at least once during [t_{0}, t'] from the graph. As there are only finite number of network structures for the limited number of genes and binding sites, we can enumerate them. For each structure, we can try all possible combinations of assignments of the constants from β' to the binding sites, and γ' to the substance generators and for each combination we can check the compatibility of the obtained network with the measurements. In this way, given Δ(t_{0},t'), we can construct a gene network that is compatible with it by an enumeration algorithm.
To complete the proof of the theorem, it remains to note that Δ(t_{0}, t_{m}) can be obtained from a series of measurements, for instance, by joining the points of the respective substance concentration by fragments of straight lines (i.e., c_{j}(t_{i}) is joined with c_{j}(t_{i }+ 1) for all j ∈ {1,...,n} and i ∈ {0,...,m1}). Given Δ(t_{0},t_{m}), we can construct the network by exhaustive search as described above.
Unfortunately such an enumeration algorithm needs exponential time and cannot be used in practice. We do not know if a polynomialtime reverse engineering algorithm exists for our model class. Note that even for finite state automata, the problem of finding a minimal automaton compatible with the input/output data is NPcomplete [19,20].
The theorem does not guarantee the reconstruction of the original network that has produced the concentration vectors. The method that we used in constructing the concentration change graph was very crude and can be easily improved to produce a more realistic graph (i.e., a graph that is more likely to be produced by the original network), by minimizing the number of fragments of straight lines for building the graph. Here, the notion of "more likely" is undefined. The problem of reconstructing the original network is formulated in the "open questions" section, but next, we generalize our model to nonbinary networks, and define the functioning of gene networks mathematically more precisely.
The multiple level generalization
For binary genes the control function is boolean, and consequently a gene has only two states: on or off. Also, the binding states have only two states. In the general case we assume that a binding site can bind more than one substance, and consequently has more than two states. We assume that the binding is exclusive, i.e., binding of one substance makes binding of any other substance impossible. In this way a binding site can either be in the detached state (denoted by 0), or in any of the attached states 1, 2, ..., p, characterized by the substance that is bound. For a given binding site b that can bind p substances, each substance has separate association and dissociation constants a_{h }and d_{h}, where h ∈ {1, ..., p}. In this way a generalized binding site can be described by a finite state automaton of the type given in Figure 1, right.
We also assume that a gene can have several expression levels {0, ..., k} (the 0 level usually meaning that the gene is not expressed). For this we assume that the control function F may have more than two values, i.e., instead of being a boolean, the function F maps an ntuple of finite values, to a finite value from 0 to k (i.e., F_{i}: ({0, ...,m_{1}}, ..., {0, ..., m_{n}}) → {0,..,k}). Respectively the gene can have k+1 states, and there are k+1 different concentration change rates r_{0},...,r_{k}, i.e., the substance generator has the form r = (i, r_{0},...,r_{k}). The concentration change rate of substance i is defined by the value of F(q_{1}, ..., q_{k}), where q_{1}, ..., q_{k }are the binding sites of the gene. Concretely, if F(q_{1}, ..., q_{k}) = j, then the rate equals to r_{j}.
Finally, we can also assume that genes can produce more than one substance, therefore in the general case a gene is defined as a triple G = (B, F, R), where R = {r_{1}, ..., r_{p}} and ri are the substance generators. We assume that all the substances are different (two genes cannot produce the same substance). In the graphical representation this implies that the dotted line coming out of a control function can fork to more than one substance generator (for instance, see Figure 6). In general, all the lines can fork, but they are not allowed to merge (they combine either through a control function or entering the same binding site). A dotted line leaving a binding site can enter one ore more control functions, a dotted line leaving a control function can enter one or more substance generators, and a solid line leaving a generator can enter one or more binding sites. The control functions can be regarded as defining the logics of the network, while binding sites and substance generators are mediators transforming discrete values into concentration change rates, and concentrations back into discrete values, respectively. Together with binding sites, the control function defines promoter (B, F) of gene G = (B,F,R).
Figure 6. Informal description of lambdaphage using the elements defined by our model (for further description see text)
Functioning of a gene network and simulations
The notion of binding site state vector can be generalized for multilevel networks in a straightforward way (by changing a binary vector to a vector of integers representing the states of the respective binding sites at the given moment). The notion of the compatibility of the binding site state and concentration vectors can also be easily generalized to multilevel situation. Further, we can assume that all the control functions F_{i }in the gene network have the binding site state vector Q = (q_{1}(t), ..., q_{n}(t)) as the argument (each function F_{i }can be changed to n argument function by adding dummy arguments for those binding sites which actually do not affect the gene). Let
Σ^{(i) }= (C(t_{i}), Q^{(i)})
be a compatible environment state, for i ≥ 0. We define the linear concentration change corresponding to state Σ^{(i) }as follows. For a substance j and gene G = (B,F,R), where R = {r_{1},...,r_{h},...,r_{m}} and R_{h }= (j, r_{h,1}, ..., r_{h,k}), for t ≥ t_{i }we set
c_{j}(t) = c_{j}(t_{i}) + (t  t_{i}) r_{h,j},
where j = F(Q^{(i)}). Let t = t_{i+1 }be the smallest t > t_{i}, such that (C(t), Q^{(i)}) is not a compatible state. Let b_{j1},...,b_{jp }be the binding sites the states of which are not compatible with C(t_{i+1}). Let Q^{(i+1) }be obtained from Q^{(i) }by changing the states q_{j1}, ..., q_{jp }to compatible ones. (In principle, there may be more than one way how this can be achieved  we can assume that we always change to the compatible state with the smallest number. This situation will not occur in the probabilistic generalization discussed in the next section.)
Let Σ^{(i) }= (C(t_{i+1}),Q^{(i+1)}). Then, given the initial compatible environment state Σ^{(i) }= (C(t_{0}),Q^{(0)}), the environment changes in the described manner for i = 0,1,.... The environment behaviour can be visualized as in the example in Figure 4.
We say that promoter (B,F) of gene G = (B,F,R) is active at a given time point t, if at this timepoint the concentration of the substance produced by the gene G is increasing.
Already with only a few genes the calculation of the network behaviour becomes quite laborious. Therefore we implemented a simulator ("Genenet") for these networks in JAVA. Figure 5 left shows the behaviour of a gene network consisting of only two genes, as depicted on the right of Figure 5. Both genes have a negative feedback loop to themselves. The first gene has an additional negative feedback onto the second gene, while the second gene has an additional positive feedback onto the first one (Figure 5, right). This example demonstrates that a very simple network of just two genes may show a nontrivial behaviour.
Figure 5. Left: Output of the simulation program "Genenet" (using Gnuplot for visualisation) Right: corresponding network; abbreviations: a1 stands for association constant 1, belongs to the bindingsite with the a1, d1 label, d1 is the corresponding dissociation constant; a2, d2, a3, d3, a4, d4 correspondingly
A model of lambdaphage
The model defined above was designed to describe processes involved in transcriptional regulation. Many additional cellular processes can be involved in gene regulatory networks. This makes some extensions necessary. With minor changes the model can be extended to allow the description of cellular processes like protein degradation. Some informal extensions are made to improve the readability for humans. The shaded boxes indicate how many different output states a controlfunction can have. The default value is 0,1 indicating the two possible states of the substance generator ON and OFF. But more states are possible, e.g. OFF, weak activity ON1, strong activity ON2. We demonstrate the usage of our model by describing a simplified model of lambdaphage.
lambdaphage
A lambdaphage has two modes of operating: lysis and lysogeny (for instance see [1]). During the infection of the bacterial cell by the phage a complex decision is made for either lysis or lysogeny. In the lysogenic mode the phage DNA is integrated into the bacterial genome, and the gene for lambdarepressor cI is the only expressed phage gene. External influences can trigger the switch from lysogenic to lytic behaviour. In the lytic mode the phage DNA is replicated, excised, new phage particles are produced and in the end the bacterium is broken open (lysed) to release the new phages. The lysislysogeny decision network is well studied and known to involve several cascades of events. In Figure 6 we present a simplified genetic network the lambdaphage. To make the graph more readable, we do not draw the lines between substance generators (depicted by diamonds) and the related bindingsites (depicted by triangles) but instead label them by the respective substances. We also allow more freedom to introduce connections between controlfunctions.
The mode of a lambdaphage operating is essentially determined by two proteins CI and Cro. If CI is in abundance, the phage is in lysogenic mode, if Cro is in abundance, the phage is in lytic mode. Both genes are regulated by the same DNA region (but transcribed in opposite directions), which has three binding sites: O_{R1}, O_{R2 }and O_{R3}. Each binding site can bind either Cro or CI competitively, but with different affinities. In this way each binding site can be in one of three states  unbound, Crobound, or CIbound. Depending on these states the control functions P_{R }and P_{M }have different activity levels. The circuit functions like a trigger and has two stable state: either cro is transcribed and cI is downregulated, or vice versa. The regulatory cascades of the lambdaPhage are quite complex, for reference please see [1,21]. We will now go through a simplified description (Figure 6).
On infection of the E. colicell by the lambdaPhage, only two promoters PL and PR of the lambdaGenome are active. From promoter PL the expression of N and CIII are initiated. Between both coding regions there is a leaky terminator of transcription located. Therefore CIII is produced at a lower rate than N. A second terminator is located between the coding region for CIII and Xis. This terminator is completely stopping transcription. If the concentration of N is high enough, the RNApolymerase is able to ignore the terminators and the genes are expressed at the same rate. As it will be important later, transcription from PL can be repressed by CI binding to its CI bindingsite.
The basal activity of promoter PR leads to the expression of cro and at lower level of O, P, cII, because there is also a terminator site located. Q is not expressed, because of a second terminator located upstream of it.
For the lysislysogeny decision CII is the crucial protein. It is protected by CIII from degradation by cellular enzymes. Thus, the concentration of CII depends on its rate of production, the activity of cellular proteinases and the concentration of CIII.
The promoters PE, PI, PM are active only, if enough CII is present to bind to them. Promoter PI initiates the expression of int. The Int protein is important for the integration of the phageDNA into the host genome. Promoter PE with CII leads to the production of CI, also called lambdaRepressor. Therefore the promoter is called
 P
 E
 P
 M
At this point, the lambdaDNA is integrated into the bacterial genome and cI is the only expressed lambdaPhage gene. An autoregulation circuit for controlling the concentration of CI at a high level is established. This is called the lysogenic state. Bacterial cells at this state show immunity to superinfection with lambdaphages, because they contain enough lambdaRepressor to immediately repress the expression of the newly incoming lambdaphage genes. The CI protein, however, is prone to be degraded by some bacterial enzymes, which are expressed by the bacterial cell as stress response upon e.g. UV irradiation. When the CI concentration is rapidly decreasing because of the degradation by cellular enzymes, PR is not repressed anymore. This leads to production of Cro, the counterplayer of CI in the lambdasystem. The degradation of CI triggered by stress response proteins is depicted in our model by a circular controlfunction with an input for the stress response signal, which could actually be a bindingsite for a stress response protein.
The regulatory protein Cro activates its own promoter by competing with CI for binding to O_{R1}, O_{R2}, O_{R3}. It binds to these sites with inverse preference compared to CI. Being a selfactivating system it is leading to a rapid increase of Cro protein in the cell. Cro also allows activation of PL, leading to increasing amounts of N. N is an antiterminator which binds to the terminators mentioned before. With N the expression of cIII, xis and int is increasing rapidly. Xis and Int are needed for the excision of the lambdaphageDNA from the bacterial genome. From PR not only cro is expressed, but also O, P, cII. O and P are needed for DNA replication of the lambdaPhage. With N these genes are produced at a significantly higher rate than without. N also allows the expression of Q. Q is an antiterminator for structural genes coded downstream of promoter PR'. This means, once CI is degraded to sufficiently low concentrations Cro is rapidly produced and then activating the genes necessary for excision from the host DNA, DNA replication and production of new phage particles, leading to host cell lysis and setting free new infectious phage particles ("Cro is opening Pandora"s box").
A lambdaphage simulation
In our model the promoter PL is represented by the controlfunction P_{L}, its output is 1 if the CI binding site is unbound or bound by Cro and 0 if the bindingsite is bound by CI (the controlfunction would look like "if (Crobound OR unbound) return 1 (=ON), if CIbound return 0 (=OFF)"). The first terminator is modelled by introducing a controlfunction P_{L1 }which has two inputs, one from a bindingsite for N and the other one from controlfunction P_{L}. The three different possibilities for the production rate of CIII are degradation (state 0), production at lower rate (state 1, if N is not bound to P_{L1}, 80% of full rate) and production at high rate (state 2, if the bindingsite for N at P_{L1 }is occupied, full rate). Controlfunction P_{L2 }is leads to a complete stop of transcription. The input of P_{L2 }is the used to model the second terminator site. Without N this terminator output of P_{L1 }and a bindingsite for N. The output equals the input from P_{L1 }if N is bound, or is 0 if N is not bound. The controlfunction P_{int }is used to model the transcriptional control of Int. The substance Int is generated either if P_{L2 }is active or if the CII binding site of P_{L2 }is occupied.
The implementation of the lambdaswitch in the model is achieved in a similar way. The binding sites O_{R1}, O_{R2 }and O_{R3 }can be bound by substance Cro or substance CI and are shared by the controlfunctions P_{R }and P_{M}. The association and dissociation constants for these substances to these bindingsites differ, allowing preferential binding in opposite order.
Using the simulator it is possible to run a simulation of the lambdaphage. Just using a quite arbitrary parameter set leads to the expected behaviour. In the beginning all substances are produced to a higher or lesser extend. After some time there are smaller changes of substance production, some kind of steady state is reached (we will refer to this informally as "behaviour"). Over a wide range of parameter sets we so far only found two principally different "behaviours". One possible outcome is a steady state where only CI is produced. We will refer to this as lysogeny state (Figure 7, top). The other one reaches a steady state where CI and CII are not produced but the other substances are(Figure 7, bottom). To this we will refer to as lytic state. The lytic behaviour shows downregulation of substance CI and upregulation of the other substances under control of substance Cro. Some of these are regulated by a negative feedback loop and are limited to a certain concentration. Some of the others are growing infinitely. The lysogenic behaviour is exemplified by downregulation of all substances besides CI which shows cyclic up and downregulation because of the feedback loop controlling its production/degradation. Interesting is to see, that at first the substances are upregulated and until the "decision making" has taken place. Depending on the concentration of substance Cro and substance CI either lytic or lysogenic "behaviour" is selected. By changing the starting values for the rate of production of substance CII we can trigger the model into lytic or lysogenic behaviour. This reflects some property of the "real" lambdaphage, the dependence on the number of phage particles infecting one cell. If this number is high (about 10 phage particles per cell) the preference is for lysogeny otherwise for lysis. In our model having several substance generators producing the same substance at a low rate it is equivalent to having one substance generator producing the substance at the according higher rate.
Figure 7. Simulation of lambdaphage model leading to lysogenic behaviour (top) or lytic behaviour (bottom)
The simulator allows to test for the effects of mutations easily, thus it is possible to experiment with the model and compare the simulations with the real mutants.
The potentials of the lambda model have to be examined further, for example for what range of parameter sets we get similar behaviours and how many different kind of behaviours we can find. But already using only arbitrary numbers gives promising results. What seems to be a shortcoming of the lambdaphage model is the infinite growth of some substances (e.g. Int, Q). But this might as well be some property of the lambdaphage itself, because it appears in the lytic "behaviour" and this leads finally to the lysis of the host cell. There is not strict need for a feedback control e.g. of the proteins responsible for the lysis of the cell as the major function of these proteins is to kill the cell. The next challenge would be to find parameters which are derived from experimentally measured reaction constants. But the purpose of this model and simulation is rather to illustrate how the model is working in principle than to come up with a new lambdaphage study.
It is obvious that additions to the model are necessary to get closer to the reality.
Informally we introduced in Figure 6 already a new kind of controlfunctions which are depicted by circles to stress that this is not an action which takes place on a promoter site. These controlfunctions can have the different current concentrations of substances (depicted by smaller circle labelled with the corresponding substance name) as an input and a substance generator of a different substance as an output. Thus we can model the influence of cellular components on the concentration of a substance, like for instance, a certain proteinase on the concentration of its substrate. This is depicted in our model by the circular controlfunction with input sites for CIII, CII and other cellular influences. It is important to add that this feature is not yet added to the simulator and not included in the simulation shown in Figure 7.
In the deterministic model, the state of the network is fully determined by its initial state and initial concentrations. To model the behaviour of the decisionmaking realistically [21], we need to introduce a stochastic element in the model.
Instead of setting precise thresholds for switching from detached to attached state and vice versa, we treat these switches as probabilistic events: the higher the concentration, the higher the probability of switching to attached state, and smaller to detached state, and vice versa. In this way a binary site can be defined as a triple B = (i,A,D), where as before i is the number of the substance that can bind to B, but A and D are two probability distributions, defining the probabilities of B switching from a detached to attached state and vice versa, respectively, depending on the concentration c_{i}.
Open questions
We would like to extend our model with some informal elements to allow description of the regulatory processes that may not be fully understood yet or may be too complicated for formal incorporation into the model. The extended model can be regarded as a semiformal language for depicting generegulatory networks. The goals of such a semiformal language are twofold: finding a semiformal description of a network is the first step towards building a completely formal model which can be used for simulation (i.e., to "describe" the network to a computer) and at the same time it helps to depict the regulatory network in a systematic way (to describe regulatory networks to other humans). Note that such a semiformal approach is often used in business modelling, where a formal graphbased description, which allows simulations of the given business process, are supplemented with informal comments, that can be interpreted only by humans.
As already noted, the formulation of the reverse engineering problem given in Section 3 is not entirely satisfactory, as it does not necessarily lead to the reconstruction of the "correct" network. A more satisfactory formulation involves assuming that the data have been produced by some unknown regulatory network (a black box), and the task is to find that or an equivalent network. For this, first, we need to define the equivalence of gene networks.
Let Γ_{1 }and Γ_{2 }be two gene networks and let Σ_{1}(t_{0}) and Σ_{2}(t_{0}) be their compatible starting states at time point t_{0}. Let Σ_{i}(t) = (C_{i}(t),Q_{i}(t)), for i = 1,2. We say that Γ_{1 }and Γ_{2 }are equivalent for the starting states. Σ_{1}(t_{0}) and Σ_{2}(t_{0}), if C_{1}(t_{0}) = C_{2}(t_{0}) implies C_{1}(t) = C_{2}(t) for all t > t_{0}. We say that Γ_{1 }and Γ_{2 }are equivalent, if they are equivalent for every compatible starting states Σ_{1}(t_{0}) and Σ_{2}(t_{0}), for which C_{1}(t_{0}) = C_{2}(t_{0}). We can also define an approximate equivalence, or more precisely, d  equivalence for a constant d ≥ 0. For this the requirement that C_{1}(t) = C_{2}(t) is relaxed to C_{1}(t)  C_{2}(t) ≤ d.
We define the reverse engineering problem in the strict sense in the following way. Let Γ be an unknown gene network and let Σ(t_{0}) = (C(t_{0}),Q(t_{0})) be its compatible starting state. We are allowed to measure the concentration state vector C(t) at any given timepoint t ≥ t_{0}. The task is to find time points t_{1}, t_{2}, ..., t_{n}, such that a network Γ' equivalent to Γ for the given starting state can be constructed from the measurements C(t_{1}), C(t_{2}), ..., C(t_{n}).
A generalized version of the problem is to find Γ' equivalent to Γ if we are allowed to choose arbitrary compatible starting states, and make series of concentration measurements for each of these states. Finally, a more practical problem is to find a network dequivalent to Γ, from approximate measurements.
At the moment we do not know if these problems are algorithmically solvable or not, even by an enumeration algorithm. They have a certain analogy with the problem of restoring a finite state automata from experiments[18]. This is algorithmically solvable, but is NPhard [19,20]. Despite the analogy, situation with the finite state linear networks are different form finite state automata in many respects.
Our theorem on the reverse engineering of gene networks gives us grounds for optimism that the reverse engineering problem for gene networks can be solved, still it is likely that heuristic methods will be needed for doing this in practice. To reconstruct gene networks all available background knowledge, such as knowing which binding sites belong to which gene promoters, will have to be used. Therefore systematic studies for regulatory signals in genomes, such as [22], will complement the approach followed here.
Acknowledgments
The authors benefited from discussing the gene regulation with Frank Holstege and Jaques van Helden, and discussing the model with Mathieu Louis and Jaak Vilo. Martin Vingron pointed me to logistic differential equations. A. B.s former colleagues from the Institute of Mathematics and Computer Science at the University of Latvia gave valuable insights into the problems of restoring general objects from particular examples. A conversation with Chris Sander was highly motivating. The general idea was sparked by the talk of Richard Karp in ISMB99 conference.
References

Ptashne M: A genetic switch; phage lambda and higher organisms. Oxford: Blackwell Science; 1992.

Thieffry D, Colet M, Thomas R: Formalization of regulatory networks: a logical method and its automation.

Szallasi Z: Genetic network analysis in light of massively parallel biological data acquisition.
Pac Symp Biocomput 1999, 516. PubMed Abstract

Akutsu T, Miyano S, Kuhara S: Identification of genetic networks from a small number of gene expression patterns under the Boolean network model.
Pac Symp Biocomput 1999, 1728. PubMed Abstract

Liang S, Fuhrman S, Somogyi R: Reveal, a general reverse engineering algorithm for inference of genetic network architectures.
Pac Symp Biocomput 1998, 1829. PubMed Abstract

Szallasi Z, Liang S: Modeling the normal and neoplastic cell cycle with "realistic Boolean genetic networks": their application for understanding carcinogenesis and assessing therapeutic strategies.
Pac Symp Biocomput 1998, 6676. PubMed Abstract

Chen T, He HL, Church GM: Modeling gene expression with differential equations.
Pac Symp Biocomput 1999, 2940. PubMed Abstract

D'Haeseleer P, Wen X, Fuhrman S, Somogyi R: Linear modeling of mRNA expression levels during CNS development and injury.
Pac Symp Biocomput 1999, 4152. PubMed Abstract

Smolen P, Baxter DA, Byrne JH: Modeling transcriptional control in gene networks  methods, recent results, and future directions.
Bull Math Biol 2000, 62:24792. PubMed Abstract  Publisher Full Text

Mendoza L, Thieffry D, AlvarezBuylla ER: Genetic control of flower morphogenesis in Arabidopsis thaliana: a logical analysis.
Bioinformatics 1999, 15:593606. PubMed Abstract  Publisher Full Text

Akutsu T, Miyano S, Kuhara S: Algorithms for inferring qualitative models of biological networks.
Pac Symp Biocomput 2000, 293304. PubMed Abstract

Thieffry D, Thomas R: Qualitative analysis of gene networks.
Pac Symp Biocomput 1998, 7788. PubMed Abstract

Thomas R: Regulatory Networks Seen as Asynchronous Automata: A Logical Description.

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:71728. PubMed Abstract  Publisher Full Text

DeRisi JL, Iyer VR, Brown PO: Exploring the metabolic and genetic control of gene expression on a genomic scale.
Science 1997, 278:6806. PubMed Abstract  Publisher Full Text

Eisen MB, Spellman PT, Brown PO, Botstein D: Cluster analysis and display of genomewide expression patterns.
Proc Natl Acad Sci U S A 1998, 95:148638. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Brazma A, Robinson A, Cameron G, Ashburner M: Onestop shop for microarray data.
Nature 2000, 403:699700. PubMed Abstract  Publisher Full Text

Moore EF: Gedankenexperiments on sequential machines. In In Automata Studies. Edited by Shannon CE, McCartney J. Princeton University Press; 1956:129153.

Angluin D: On the Complexity of Minimum Inference of Regular Sets.

Gold EM: Complexity of Automaton Identification from Given Data.

McAdams HH, Shapiro L: Circuit simulation of genetic networks.
Science 1995, 269:6506. PubMed Abstract

Brazma A, Jonassen I, Vilo J, Ukkonen E: Predicting gene regulatory elements in silico on a genomic scale.
Genome Res 1998, 8:120215. PubMed Abstract  Publisher Full Text