# Reverse engineering of gene regulatory networks: a finite state linear model

- Alvis Brazma
^{1}and - Thomas Schlitt
^{1}Email author

**4**:P5

**DOI: **10.1186/gb-2003-4-6-p5

© BioMed Central Ltd 2003

**Received: **14 April 2003

**Published: **29 April 2003

## 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 lambda-phage. We also describe some generalizations of the model, discuss their relevance to the real-world gene networks and formulate a number of open problems.

### Keywords

gene regulation regulatory networks regulatory circuits dynamic systems finite state automata reverse engineering## Background

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 lambda-phage 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 [4–6], 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 time-point n-1. 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 non-additive 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 [9–11]. 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 feed-back 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 data-often 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 [14–16] 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 lambda-phage 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),

*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

G = (B, F, r),

_{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.

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.

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(1-c/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 time-point t' can be regarded as a concentration measurement. Let us define a *measurement series* as a pair of m-tuples

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.

_{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.

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 transition-points 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,...,m-1}). 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 polynomial-time 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 NP-complete [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 non-binary 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 n-tuple 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}.

_{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).

### Functioning of a gene network and simulations

The notion of binding site state vector can be generalized for multilevel networks in a straight-forward 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 time-point the concentration of the substance produced by the gene G is increasing.

### A model of lambda-phage

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 control-function 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 lambda-phage.

### lambda-phage

A lambda-phage 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 lambda-repressor *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 lysis-lysogeny decision network is well studied and known to involve several cascades of events. In Figure 6 we present a simplified genetic network the lambda-phage. 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 control-functions.

The mode of a lambda-phage 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, *Cro*-bound, or *CI*-bound. 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 down-regulated, or vice versa. The regulatory cascades of the lambda-Phage are quite complex, for reference please see [1, 21]. We will now go through a simplified description (Figure 6).

On infection of the *E. coli*-cell by the lambda-Phage, only two promoters PL and PR of the lambda-Genome 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 RNA-polymerase 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 lysis-lysogeny 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 phage-DNA into the host genome. Promoter PE with *CII* leads to the production of *CI*, also called lambda-Repressor. Therefore the promoter is called Promoter repressor Early (PE). *CI* binds to the operator sites O_{R1} and O_{R2} in promoter PR and to PL, thus blocking transcription from PR, PL, PE. But it activates its own synthesis via promoter PM (Promoter for repressor Maintenance). Thus the single gene for *cI* can be either transcribed from PM or PE. Actually these promoters are serially organized on DNA level. The promoters PM and PR are sharing the operator sites O_{R1}, O_{R2}, O_{R3}. These sites are bound by increasing concentrations of *CI*. Binding to O_{R1} and O_{R2} leads to inactivation of PR and activation of PM. However, binding to O_{R3} at even higher *CI* concentrations leads to inactivation of PR and PM, thus down-regulating its own expression.

At this point, the lambda-DNA is integrated into the bacterial genome and *cI* is the only expressed lambda-Phage gene. An auto-regulation 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 super-infection with lambda-phages, because they contain enough lambda-Repressor to immediately repress the expression of the newly incoming lambda-phage 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 counter-player of *CI* in the lambda-system. The degradation of *CI* triggered by stress response proteins is depicted in our model by a circular control-function 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 self-activating 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 anti-terminator 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 lambda-phage-DNA 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 lambda-Phage. With *N* these genes are produced at a significantly higher rate than without. *N* also allows the expression of *Q*. *Q* is an anti-terminator 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 lambda-phage simulation

In our model the promoter PL is represented by the control-function 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 control-function would look like "if (Cro-bound OR unbound) return 1 (=ON), if CI-bound return 0 (=OFF)"). The first terminator is modelled by introducing a control-function P_{L1} which has two inputs, one from a bindingsite for N and the other one from control-function 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). Control-function 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 control-function 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 lambda-switch 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 control-functions P_{R} and P_{M}. The association and dissociation constants for these substances to these bindingsites differ, allowing preferential binding in opposite order.

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 lambda-phage model is the infinite growth of some substances (e.g. Int, Q). But this might as well be some property of the lambda-phage 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 lambda-phage 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 control-functions which are depicted by circles to stress that this is not an action which takes place on a promoter site. These control-functions 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 control-function 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 decision-making 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 semi-formal language for depicting gene-regulatory networks. The goals of such a semi-formal language are twofold: finding a semi-formal 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 semi-formal approach is often used in business modelling, where a formal graph-based 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 time-point 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 d-equivalent 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 NP-hard [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.

## Declarations

### 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.

## Authors’ Affiliations

## References

- Ptashne M: A genetic switch; phage lambda and higher organisms. 1992, Oxford: Blackwell ScienceGoogle Scholar
- Thieffry D, Colet M, Thomas R: Formalization of regulatory networks: a logical method and its automation. Math Model Sci Comput. 1993, 55: 144-151.Google Scholar
- Szallasi Z: Genetic network analysis in light of massively parallel biological data acquisition. Pac Symp Biocomput. 1999, 5-16.Google Scholar
- 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, 17-28.Google Scholar
- Liang S, Fuhrman S, Somogyi R: Reveal, a general reverse engineering algorithm for inference of genetic network architectures. Pac Symp Biocomput. 1998, 18-29.Google Scholar
- 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, 66-76.Google Scholar
- Chen T, He HL, Church GM: Modeling gene expression with differential equations. Pac Symp Biocomput. 1999, 29-40.Google Scholar
- D'Haeseleer P, Wen X, Fuhrman S, Somogyi R: Linear modeling of mRNA expression levels during CNS development and injury. Pac Symp Biocomput. 1999, 41-52.Google Scholar
- Smolen P, Baxter DA, Byrne JH: Modeling transcriptional control in gene networks - methods, recent results, and future directions. Bull Math Biol. 2000, 62: 247-92. 10.1006/bulm.1999.0155.PubMedView ArticleGoogle Scholar
- Mendoza L, Thieffry D, Alvarez-Buylla ER: Genetic control of flower morphogenesis in Arabidopsis thaliana: a logical analysis. Bioinformatics. 1999, 15: 593-606. 10.1093/bioinformatics/15.7.593.PubMedView ArticleGoogle Scholar
- Akutsu T, Miyano S, Kuhara S: Algorithms for inferring qualitative models of biological networks. Pac Symp Biocomput. 2000, 293-304.Google Scholar
- Thieffry D, Thomas R: Qualitative analysis of gene networks. Pac Symp Biocomput. 1998, 77-88.Google Scholar
- Thomas R: Regulatory Networks Seen as Asynchronous Automata: A Logical Description. J theor Biol. 1991, 153: 1-23.View ArticleGoogle Scholar
- 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: 717-28.PubMedView ArticleGoogle Scholar
- DeRisi JL, Iyer VR, Brown PO: Exploring the metabolic and genetic control of gene expression on a genomic scale. Science. 1997, 278: 680-6. 10.1126/science.278.5338.680.PubMedView ArticleGoogle Scholar
- Eisen MB, Spellman PT, Brown PO, Botstein D: Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci U S A. 1998, 95: 14863-8. 10.1073/pnas.95.25.14863.PubMedPubMed CentralView ArticleGoogle Scholar
- Brazma A, Robinson A, Cameron G, Ashburner M: One-stop shop for microarray data. Nature. 2000, 403: 699-700. 10.1038/35001676.PubMedView ArticleGoogle Scholar
- Moore EF: Gedanken-experiments on sequential machines. In Automata Studies. Edited by: Shannon CE, McCartney J. 1956, Princeton University Press, 129-153.Google Scholar
- Angluin D: On the Complexity of Minimum Inference of Regular Sets. Inform Control. 1978, 39: 337-350.View ArticleGoogle Scholar
- Gold EM: Complexity of Automaton Identification from Given Data. Inform Control. 1978, 37: 302-320.View ArticleGoogle Scholar
- McAdams HH, Shapiro L: Circuit simulation of genetic networks. Science. 1995, 269: 650-6.PubMedView ArticleGoogle Scholar
- Brazma A, Jonassen I, Vilo J, Ukkonen E: Predicting gene regulatory elements in silico on a genomic scale. Genome Res. 1998, 8: 1202-15.PubMedPubMed CentralGoogle Scholar